家豬長(zhǎng)鏈非編碼RNA的組織特異性表達(dá)及其對(duì)肉質(zhì)的影響

王浩杰， 張珍華， 李琰， 張霞， 彭銳

（四川大學(xué)生命科學(xué)學(xué)院，生物資源與生態(tài)環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室，成都610065）

哺乳動(dòng)物的基因組中約2/3可轉(zhuǎn)錄，而僅約1.9%的部分可合成蛋白（Mattick，2001；Djebaliet al.，2012）。多數(shù)非編碼基因的轉(zhuǎn)錄產(chǎn)物稱為非編碼RNA（non-coding RNA，ncRNA），其中長(zhǎng)度超過(guò)200 nt的RNA被定義為長(zhǎng)鏈非編碼RNA（long non-coding RNA，lncRNA）。與 mRNA相同的是，lncRNA也是由RNA聚合酶Ⅱ轉(zhuǎn)錄，具有5’-和3’-poly A尾部的結(jié)構(gòu)（Mercer & Mattick，2013；Marcheseet al.，2017）。除了無(wú)法編碼蛋白外，與mRNA相比，lncRNA有表達(dá)量更低、更高的組織特異性表達(dá)模式等的特點(diǎn)（Mongelliet al.，2019）。LncRNA曾被視為轉(zhuǎn)錄組的“噪音”，但它可通過(guò)與蛋白或蛋白編碼基因（protein-coding gene，PCG）等生物大分子作用，從而參與一系列重要的生物過(guò)程，如轉(zhuǎn)錄調(diào)控、細(xì)胞分化和細(xì)胞自噬等（Gil & Ulitsky，2020；Sunet al.，2020）。

在包括人類和家畜在內(nèi)的多個(gè)物種中，lncRNA都可呈現(xiàn)出組織特異性的表達(dá)特點(diǎn)（Tsoiet al.，2015；Kernet al.，2018；Cavaet al.，2019）。一些組織特異性lncRNA（tissue-specific lncRNA，TS lncRNA）已被證實(shí)對(duì)其所在的組織發(fā)育有著重要的影響，Grote等（2013）在小鼠Mus musculus的胚胎細(xì)胞中發(fā)現(xiàn)了特異性表達(dá)于側(cè)中胚層的TS lncRNA Fendrr，并證實(shí)對(duì)Fendrr的沉默會(huì)導(dǎo)致小鼠的心臟和體壁發(fā)育異常。Morán等（2012）篩選了人胰島內(nèi)的TS lncRNA，并發(fā)現(xiàn)它們是胰島β細(xì)胞分化過(guò)程中的重要組成，還發(fā)現(xiàn)了lncRNA HILNC25對(duì)糖尿病相關(guān)基因GLIS3的正調(diào)控作用。

豬肉不僅在肉類消費(fèi)市場(chǎng)上占有重要地位，也是人類攝入蛋白的重要來(lái)源。因此，提升豬肉的品質(zhì)對(duì)生豬養(yǎng)殖業(yè)十分重要。肌肉內(nèi)脂肪（intramuscular fat，IMF）是影響豬肉品質(zhì)的重要因素，較高的IMF含量可以提升豬肉的多汁度、柔軟度和風(fēng)味（Daszkiewiczet al.，2005）。影響IMF含量的因素包括品種、年齡、飼養(yǎng)環(huán)境、飼料中蛋白質(zhì)含量等（Parket al.，2018；Malgwiet al.，2022），然 而lncRNA-mRNA的相互作用對(duì)IMF的影響還不清楚。

本研究利用公共數(shù)據(jù)庫(kù)中的RNA-seq數(shù)據(jù)，重新組裝了家豬Sus scrofa domestica的轉(zhuǎn)錄組并鑒定了lncRNA，篩選了在肌肉、脂肪、脾臟、腎臟、肝臟和卵巢中的TS lncRNA并預(yù)測(cè)了其功能。此外，還鑒定了影響其IMF含量的lncRNA并建立了相關(guān)的調(diào)控網(wǎng)絡(luò)。這些數(shù)據(jù)為生豬養(yǎng)殖業(yè)和未來(lái)lncRNA的研究提供了參考。

1 材料與方法

1.1 數(shù)據(jù)來(lái)源

參考基因組（ftp：//ftp.ensembl.org/pub/release-98/fasta/sus_scrofa/）和注釋文件（ftp：//ftp.ensembl.org/pub/release-98/gtf/sus_scrofa/）均 來(lái) 自 Ensembl數(shù)據(jù)庫(kù)。RNA-seq數(shù)據(jù)均下載于基因表達(dá)（Gene Expression Omnibus，GEO）數(shù)據(jù)庫(kù)，包括6個(gè)組織或器官（肌肉、脂肪、脾臟、肝臟、腎臟和卵巢）共36組用于 lncRNA 鑒定（Fushanet al.，2015；Draget al.，2017；Liet al.，2017；Limet al.，2017；TangQZet al.，2017；TangZLet al.，2017；Cheet al.，2018；Kimet al.，2018；Oczkowiczet al.，2018；Liuet al.，2019；Wanget al.，2019）。3組高IMF含量和3組低IMF含量的轉(zhuǎn)錄組測(cè)序數(shù)據(jù)（Limet al.，2017）用于差異表達(dá)分析。

表1 數(shù)據(jù)集信息匯總Table 1 The summary of datasets

1.2 轉(zhuǎn)錄組組裝

RNA-seq 由 Trim Galore 0.6.4（Martin，2011）去除接頭和低質(zhì)量的測(cè)序片段后，使用Hisat2 2.1.0（Kimet al.，2019）將數(shù)據(jù)比對(duì)到參考基因組上并得到sam格式的比對(duì)文件。最后，使用StringTie 2.0（Perteaet al.，2015）將比對(duì)的結(jié)果進(jìn)行組裝并定量，合并后得到gtf格式的轉(zhuǎn)錄組注釋文件。

1.3 LncRNA的鑒定

首先將注釋文件中被標(biāo)記為“l(fā)ncRNA”的轉(zhuǎn)錄本視作已知的lncRNA。然后使用Cuffcompare 2.2.1工具將所有新鑒定的轉(zhuǎn)錄本進(jìn)行分類，獲取其中被標(biāo)記為“i”“j”“u”“x”和“o”5 種類型之一的轉(zhuǎn)錄本，并過(guò)濾掉長(zhǎng)度小于200 nt以及外顯子數(shù)量＜2的轉(zhuǎn)錄本。隨后，分別使用CNCI（Sunet al.，2013）、PLEK 1.2（Liet al.，2014）和 Pfam（Batemanet al.，2004）對(duì)上一步所得轉(zhuǎn)錄本進(jìn)行編碼能力預(yù)測(cè)，獲取在3種方法中都顯示無(wú)法編碼蛋白的轉(zhuǎn)錄本。

1.4 TS lncRNA的鑒定

使用 Tspex 0.6.1（https：//tspex.lge.ibi.unicamp.br/）計(jì)算每個(gè)lncRNA的組織特異性指數(shù)（tissue specific index，TSI），TSIi=，式中，N為組織的總數(shù)，xi為轉(zhuǎn)錄本x在第i個(gè)組織中經(jīng)歸一化的表達(dá)量（Julienet al.，2012）。TSI的值在0～1之間，越接近1表明組織特異性程度越高。將TSI＞0.9且在所有樣本中的平均表達(dá)量＞0.1 FPKM的lncRNA視作具有組織特異性。

1.5 加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析

使 用 R 包 WGCNA（Langfelder & Horvath，2008）構(gòu)建加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析（weighted gene co-expression network analysis，WGCNA）。首先，根據(jù)表達(dá)量計(jì)算每對(duì)轉(zhuǎn)錄本之間的Pearson系數(shù)，選擇5為軟閾值并構(gòu)建鄰接矩陣。其次，根據(jù)鄰接矩陣計(jì)算拓?fù)渲丿B度并得到拓?fù)渲丿B矩陣，再將其轉(zhuǎn)換為相異矩陣并據(jù)其構(gòu)建層次聚類樹，通過(guò)動(dòng)態(tài)樹切割將表達(dá)模式相似的轉(zhuǎn)錄本聚集到同一模塊。隨后，計(jì)算出每個(gè)組織與模塊之間的相關(guān)性，值越接近1或-1表示該組織與模塊之間的相關(guān)性越高。最后，將與每個(gè)TS lncRNA相關(guān)性最高的10個(gè)mRNA視作其潛在的反式作用位點(diǎn)。

1.6 功能富集分析

首先使用R包的BiomaRt 2.42.1（Smedleyet al.，2009）將每個(gè)轉(zhuǎn)錄本轉(zhuǎn)換為與之對(duì)應(yīng)的Entrez id，隨后使用 R 包的 ClusterProfiler 3.14.0（Yuet al.，2012）進(jìn)行基因本體GO和KEGG富集分析，結(jié)果中P-value≤0.05的注釋條目視為有效結(jié)果。

1.7 影響IMF的lncRNA-mRNA篩選

首先利用 R 包的 DESeq2 1.26.0（Loveet al.，2014）對(duì)不同IMF水平的RNA-seq進(jìn)行差異分析。結(jié)果中l(wèi)og2（Fold Change）≥1且FDR≤0.05的mRNA被視作顯著差異，而log2（Fold Change）≥1且P-value≤0.05的lncRNA被視作顯著差異。其次，計(jì)算所有差異表達(dá)的lncRNA和mRNA之間的Pearson相關(guān)系數(shù)，選擇Pearson系數(shù)絕對(duì)值＞0.9的轉(zhuǎn)錄本作為共表達(dá)的lncRNA-mRNA作用對(duì)。隨后，利用Lnc-Tar（Liet al.，2015）估算這些lncRNA-mRNA之間的歸一化結(jié)合自由能（ndG），以-0.13作為閾值篩選出可能相互作用的lncRNA-mRNA。最后進(jìn)行功能富集分析，將與脂質(zhì)代謝相關(guān)的lncRNA-mRNA視作潛在影響IMF的lncRNA-mRNA作用對(duì)。

1.8 LncRNA的二級(jí)結(jié)構(gòu)預(yù)測(cè)

利用 ViennaRNA 2.4.18（Lorenzet al.，2011）基于成環(huán)結(jié)構(gòu)以及動(dòng)態(tài)規(guī)劃算法預(yù)測(cè)目標(biāo)lncRNA的二級(jí)結(jié)構(gòu)：最小自由能結(jié)構(gòu)和質(zhì)心二級(jí)結(jié)構(gòu)。

1.9 RNA干擾與過(guò)表達(dá)實(shí)驗(yàn)

首先將豬PK-15細(xì)胞培養(yǎng)于DEME（高糖）培養(yǎng)基，置于37 ℃，5%CO2環(huán)境。在PLKO.1質(zhì)粒中插入shRNA構(gòu)建目標(biāo)lncRNA的干擾載體，對(duì)照組質(zhì)粒中插入scramble shRNA。過(guò)表達(dá)實(shí)驗(yàn)中，將目標(biāo)lncRNA的全長(zhǎng)cDNA插入pcDNA-3.1質(zhì)粒構(gòu)建過(guò)表達(dá)載體，對(duì)照組使用空質(zhì)粒。均使用High-Gene轉(zhuǎn)染試劑將載體轉(zhuǎn)染入細(xì)胞。

1.10 實(shí)時(shí)熒光定量PCR（RT-qPCR）

總RNA由TRIzol試劑提取并反轉(zhuǎn)錄為cDNA，隨后使用 2×Taq SYBR?Green qPCR Premix對(duì)MSTRG.7256.5、PPARγ、SCD、AQP7和SORBS1進(jìn)行RT-qPCR實(shí)驗(yàn)，選用管家基因GADPH為內(nèi)參，每個(gè)樣本3組重復(fù)，采用2-ΔΔCt法計(jì)算每個(gè)基因的相對(duì)表達(dá)水平（表2）。

表2 實(shí)驗(yàn)所用引物Table 2 Primers used in this study

1.11 油紅O染色實(shí)驗(yàn)

首先將收集的細(xì)胞用PBS緩沖液清洗3次，再用10%甲醛溶液固定5 min，使用60%異丙醇溶液清洗5 min后，將細(xì)胞用油紅O染液染色10 min，最后將細(xì)胞用蒸餾水洗滌并置于顯微鏡下觀察。

2 結(jié)果

2.1 LncRNA的鑒定與分析

共獲得23 678條lncRNA，其中，新鑒定的lncRNA有14 309條，包括3 061條反義lncRNA，3 051條基因間lncRNA，2 212條內(nèi)含子lncRNA和5 985條正義lncRNA。與mRNA相比，lncRNA的長(zhǎng)度更短，且這一點(diǎn)在新鑒定的lncRNA中更顯著（圖1：A）。在每條染色體上，mRNA的數(shù)量都明顯多于lncRNA，僅在第10條、第11條、第16條和Y染色體上，二者的數(shù)量差距不大（圖1：B）。此外，lncRNA還有著更少的外顯子數(shù)量和更低的表達(dá)量（圖1：C，D）。

圖1 LncRNA的特征Fig. 1 Characterization of lncRNAs

2.2 組織特異性lncRNA的鑒定與分析

通過(guò)計(jì)算每條lncRNA的TSI指數(shù)，共有1 218條被鑒定為具有組織特異性，其中，卵巢的TS lncRNA數(shù)量最多（429條），而肌肉的最少（85條）（圖2：A）。將TS lncRNA的位置映射到染色體上，第1條上的數(shù)量最多，Y染色體上為0，TS lncRNA的數(shù)量基本與其所在染色體的長(zhǎng)度呈正相關(guān)關(guān)系（r=4.25×10-7，P=4.62×10-6）（圖2：B，C）。TS lncRNA的表達(dá)量經(jīng)歸一化后的中位數(shù)均＜1（圖2：D）。盡管卵巢中TS lncRNA的數(shù)量最多，但表達(dá)水平最低（圖2：A，D）。

圖2 TS lncRNA的特征Fig. 2 Characteristics of TS lncRNAs

2.3 TS lncRNA的功能預(yù)測(cè)

首先，查找TS lncRNA上、下游50 kb內(nèi)的mRNA，將其視作順式作用的潛在位點(diǎn)，并進(jìn)行GO和KEGG富集分析：僅有腎臟、肝臟和脾臟中TS lncRNA的功能與對(duì)應(yīng)的組織相關(guān)。在肝臟中，顯著的GO條目有“甘油三酸酯穩(wěn)態(tài)”“化學(xué)內(nèi)穩(wěn)態(tài)”“有機(jī)酸代謝過(guò)程”等（圖3：A）。在脂肪中，顯著的GO條目與mRNA剪接、mRNA處理等相關(guān)（圖3：B）。

圖3 通過(guò)順式作用的TS lncRNA功能預(yù)測(cè)Fig. 3 Functional prediction of TS lncRNAs by cis-acting

為獲取TS lncRNA潛在的反式作用位點(diǎn)，選取所有l(wèi)ncRNA和mRNA構(gòu)建加權(quán)基因共表達(dá)分析。最終表達(dá)模式相近的轉(zhuǎn)錄本被分到相同的模塊，每個(gè)模塊被標(biāo)記為單獨(dú)的顏色，共獲得42個(gè)模塊（圖4：A）。所有TS lncRNA都位于與其所在組織高度相關(guān)的模塊（P＜0.05）。選取與每個(gè)TS lncRNA相關(guān)程度最高的10個(gè)mRNA，并對(duì)其進(jìn)行功能富集分析。結(jié)果顯示，每個(gè)組織中TS lncRNA的功能都與該組織相關(guān)，肌肉中TS lncRNA功能預(yù)測(cè)中顯著的GO條目有“肌肉組織發(fā)育”“橫紋肌細(xì)胞分化”等（圖4：B），而脾臟中對(duì)應(yīng)的GO條目有“淋巴細(xì)胞激活”“免疫系統(tǒng)過(guò)程的調(diào)節(jié)”等（圖4：C）。

圖4 通過(guò)反式作用的TS lncRNA功能預(yù)測(cè)Fig. 4 Functional prediction of TS lncRNAs by trans-action

2.4 調(diào)控IMF的lncRNA-mRNA篩選

從GEO數(shù)據(jù)庫(kù)中下載有關(guān)IMF差異的家豬RNA-seq數(shù)據(jù)，通過(guò)差異表達(dá)分析，共獲得239條顯著差異表達(dá)的mRNA和245條lncRNA（圖 5：A，B）。共篩選出 120個(gè) lncRNA-mRNA作用對(duì)，功能富集分析發(fā)現(xiàn)23對(duì)lncRNA-mRNA與脂質(zhì)代謝相關(guān)的通路有關(guān)，包括“PPAR信號(hào)通路”“脂肪酸代謝”“甘油酯代謝”等（圖5：C）。其中，新鑒定的lncRNA——MSTRG.7256.5與脂質(zhì)相關(guān)基因的連接最多，其中部分參與PPAR通路（圖5：D）。因此，推測(cè)MSTRG.7256.5可能是調(diào)控IMF含量的因素。

MSTRG.7256.5位于第12條染色體，具有2個(gè)外顯子，全長(zhǎng)300 nt，其對(duì)應(yīng)的最小自由能結(jié)構(gòu)和質(zhì)心二級(jí)結(jié)構(gòu)如圖所示（圖5：E，F(xiàn)）。

圖5 與IMF相關(guān)的lncRNA篩選Fig. 5 Identification of lncRNAs responsible for IMF

2.5 MSTRG.7256.5功能的驗(yàn)證

RT-qPCR的結(jié)果顯示，構(gòu)建的shRNA載體有效降低MSTRG.7256.5的水平后（圖6：A），AQP7、SORBS和PPARγ的表達(dá)量增加，其中，PPARγ的增加顯著，而SCD的表達(dá)量顯著下降（圖6：B）。隨后構(gòu)建并轉(zhuǎn)染了MSTRG.7256.5的過(guò)表達(dá)載體，MSTRG.7256.5的表達(dá)水平升高后（圖6：C），AQP7、SORBS和PPARγ的表達(dá)量顯著下降，而SCD的表達(dá)量顯著上升，整體上與干擾實(shí)驗(yàn)相反（圖6：D）。

圖6 MSTRG.7256.5作用的RT-qPCR驗(yàn)證Fig. 6 Correlation validation of MSTRG.7256.5 by RT-qPCR

油紅O染色實(shí)驗(yàn)結(jié)果顯示，MSTRG.7256.5干擾后的細(xì)胞內(nèi)脂滴數(shù)量高于對(duì)照組（圖7：A）。過(guò)表達(dá)了MSTRG.7256.5的細(xì)胞脂滴數(shù)量則更低（圖7：B）。這些結(jié)果表明，MSTRG.7256.5可能通過(guò)PPAR途徑調(diào)控家豬的IMF含量。

圖7 MSTRG.7256.5干擾（A）與過(guò)表達(dá)后（B）的細(xì)胞內(nèi)脂滴變化（油紅O染色實(shí)驗(yàn)）Fig. 7 Change of lipid droplets after RNA interference（A） and overexpression （B）（oil red O staining experiment）

3 討論

由于其營(yíng)養(yǎng)價(jià)值和經(jīng)濟(jì)價(jià)值，家豬已經(jīng)成為許多研究中的重要生物模型。lncRNA在許多生物過(guò)程中發(fā)揮著關(guān)鍵作用（Zhuet al.，2019），然而lncRNA-mRNA的相互作用，及其對(duì)組織特異性表達(dá)和IMF沉積影響的研究還有待深入。在本研究中，利用公共數(shù)據(jù)庫(kù)中的RNA-seq數(shù)據(jù)，綜合分析了lncRNAs的組織特異性表達(dá)模式以及對(duì)豬IMF含量的影響，為未來(lái)相關(guān)養(yǎng)殖業(yè)的研究提供了數(shù)據(jù)和參考。

在組裝了豬的轉(zhuǎn)錄組后，共鑒定出23 678個(gè)lncRNA，其中14 309個(gè)為新鑒定的。經(jīng)過(guò)特征分析，發(fā)現(xiàn)lncRNA與編碼蛋白的mRNA相比，具有長(zhǎng)度更短、表達(dá)更低、外顯子更少等特征，這與之前的研究一致（Zouet al.，2018）。新鑒定的與已注釋的lncRNA也存在明顯的差異，在其他研究中也有類似的現(xiàn)象（Liet al.，2020），這可能是由目前l(fā)ncRNA鑒定流程不夠成熟造成的。有研究報(bào)道，lncRNA與蛋白編碼基因相比具有更高的組織特異性（Cabiliet al.，2011；Kornienkoet al.，2016）。因此，本文通過(guò)計(jì)算TSI指數(shù)篩選了所有的TS lncRNA，其中，卵巢組織中的TS lncRNA數(shù)量明顯高于其他組織。曾有研究報(bào)道lncRNA在哺乳動(dòng)物卵巢和睪丸中特異性表達(dá)的數(shù)量高于其他組織（Iwakiriet al.，2017；Cavaet al.，2019），因此，推測(cè)lncRNA可能在生殖系統(tǒng)中有更高的組織特異性。

lncRNA可通過(guò)順式或反式作用調(diào)控其靶基因（Gil & Ulitsky，2020），因此通常使用2種方式預(yù)測(cè)lncRNA的功能：一是根據(jù)其基因座附近的編碼基因預(yù)測(cè)；二是根據(jù)與其共表達(dá)的基因預(yù)測(cè)。因此，本文通過(guò)尋找TS lncRNA的鄰近基因（＜50 kb）和與其共表達(dá)的mRNA來(lái)預(yù)測(cè)lncRNA的功能。富集分析的結(jié)果顯示，與順式作用相比，lncRNA的反式作用靶點(diǎn)與其所在組織功能的相關(guān)性更高。卵巢TS ln-cRNA的鄰近基因在功能富集分析中未獲得顯著結(jié)果，但與其共表達(dá)的mRNA在生殖系統(tǒng)相關(guān)的GO條目中顯著富集。推測(cè)TS lncRNA傾向于通過(guò)反式作用發(fā)揮功能。

最后通過(guò)一系列流程篩選可能調(diào)控IMF含量的lncRNA-mRNA作用對(duì)，并構(gòu)建了一個(gè)lncRNA-mRNA相互作用網(wǎng)絡(luò)。發(fā)現(xiàn)1條新鑒定的lncRNA——MSTRG.7256.5與脂質(zhì)代謝相關(guān)的PCG連接數(shù)最多，其中一些PCG富集在PPAR通路中，已有報(bào)道稱，PPAR通路與脂質(zhì)調(diào)節(jié)相關(guān)（Walczak & Tontonoz，2002；Li & Glass，2004）。 因此 ，本 文 對(duì)MSTRG.7256.5分別進(jìn)行了RNA干擾和過(guò)表達(dá)實(shí)驗(yàn)，以驗(yàn)證其與潛在靶基因的關(guān)系。結(jié)果表明，MSTRG.7256.5表達(dá)趨勢(shì)與AQP7、SORBS和PPARγ的表達(dá)呈負(fù)相關(guān)，與SCD的呈正相關(guān)，但SORBS的表達(dá)量在干擾實(shí)驗(yàn)中的變化不顯著，AQP7的表達(dá)量在干擾實(shí)驗(yàn)與過(guò)表達(dá)實(shí)驗(yàn)中的變化均不顯著。此前已有研究表明，AQP7和SORBS與脂質(zhì)代謝相關(guān)（Yanget al.，2003；Skowronskiet al.，2007），而SCD可以調(diào)節(jié)豬的脂肪酸沉積和IMF含量（Yokotaet al.，2012；Ros-Freixedeset al.，2016）。此外，油紅O染色實(shí)驗(yàn)的結(jié)果顯示，MSTRG.7256.5的表達(dá)水平與細(xì)胞內(nèi)的脂滴數(shù)量負(fù)相關(guān)，但過(guò)表達(dá)實(shí)驗(yàn)的結(jié)果不如RNA干擾實(shí)驗(yàn)的結(jié)果明顯。

綜上所述，認(rèn)為MSTRG.7256.5可能通過(guò)PPAR途徑負(fù)向調(diào)節(jié)家豬的脂肪沉積。然而，由于此步驟目的在于驗(yàn)證MSTRG.7256.5能否對(duì)其靶基因和脂肪沉積產(chǎn)生影響，因此僅選用了PK-15細(xì)胞用于實(shí)驗(yàn)，而PK-15細(xì)胞較低的脂肪含量可能導(dǎo)致實(shí)驗(yàn)結(jié)果不顯著?？傊?，MSTRG.7256.5對(duì)其靶基因具體的作用機(jī)制和對(duì)肌肉間IMF含量的影響還須通過(guò)實(shí)驗(yàn)進(jìn)一步探究。