家豬長鏈非編碼RNA的組織特異性表達及其對肉質的影響

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

（四川大學生命科學學院，生物資源與生態環境教育部重點實驗室，成都610065）

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

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

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

本研究利用公共數據庫中的RNA-seq數據，重新組裝了家豬Sus scrofa domestica的轉錄組并鑒定了lncRNA，篩選了在肌肉、脂肪、脾臟、腎臟、肝臟和卵巢中的TS lncRNA并預測了其功能。此外，還鑒定了影響其IMF含量的lncRNA并建立了相關的調控網絡。這些數據為生豬養殖業和未來lncRNA的研究提供了參考。

1 材料與方法

1.1 數據來源

參考基因組（ftp：//ftp.ensembl.org/pub/release-98/fasta/sus_scrofa/）和注釋文件（ftp：//ftp.ensembl.org/pub/release-98/gtf/sus_scrofa/）均 來 自 Ensembl數據庫。RNA-seq數據均下載于基因表達（Gene Expression Omnibus，GEO）數據庫，包括6個組織或器官（肌肉、脂肪、脾臟、肝臟、腎臟和卵巢）共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含量的轉錄組測序數據（Limet al.，2017）用于差異表達分析。

表1 數據集信息匯總Table 1 The summary of datasets

1.2 轉錄組組裝

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

1.3 LncRNA的鑒定

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

1.4 TS lncRNA的鑒定

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

1.5 加權基因共表達網絡分析

使 用 R 包 WGCNA（Langfelder & Horvath，2008）構建加權基因共表達網絡分析（weighted gene co-expression network analysis，WGCNA）。首先，根據表達量計算每對轉錄本之間的Pearson系數，選擇5為軟閾值并構建鄰接矩陣。其次，根據鄰接矩陣計算拓撲重疊度并得到拓撲重疊矩陣，再將其轉換為相異矩陣并據其構建層次聚類樹，通過動態樹切割將表達模式相似的轉錄本聚集到同一模塊。隨后，計算出每個組織與模塊之間的相關性，值越接近1或-1表示該組織與模塊之間的相關性越高。最后，將與每個TS lncRNA相關性最高的10個mRNA視作其潛在的反式作用位點。

1.6 功能富集分析

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

1.7 影響IMF的lncRNA-mRNA篩選

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

1.8 LncRNA的二級結構預測

利用 ViennaRNA 2.4.18（Lorenzet al.，2011）基于成環結構以及動態規劃算法預測目標lncRNA的二級結構：最小自由能結構和質心二級結構。

1.9 RNA干擾與過表達實驗

首先將豬PK-15細胞培養于DEME（高糖）培養基，置于37 ℃，5%CO2環境。在PLKO.1質粒中插入shRNA構建目標lncRNA的干擾載體，對照組質粒中插入scramble shRNA。過表達實驗中，將目標lncRNA的全長cDNA插入pcDNA-3.1質粒構建過表達載體，對照組使用空質粒。均使用High-Gene轉染試劑將載體轉染入細胞。

1.10 實時熒光定量PCR（RT-qPCR）

總RNA由TRIzol試劑提取并反轉錄為cDNA，隨后使用 2×Taq SYBR?Green qPCR Premix對MSTRG.7256.5、PPARγ、SCD、AQP7和SORBS1進行RT-qPCR實驗，選用管家基因GADPH為內參，每個樣本3組重復，采用2-ΔΔCt法計算每個基因的相對表達水平（表2）。

表2 實驗所用引物Table 2 Primers used in this study

1.11 油紅O染色實驗

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

2 結果

2.1 LncRNA的鑒定與分析

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

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

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

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

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

2.3 TS lncRNA的功能預測

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

圖3 通過順式作用的TS lncRNA功能預測Fig. 3 Functional prediction of TS lncRNAs by cis-acting

為獲取TS lncRNA潛在的反式作用位點，選取所有lncRNA和mRNA構建加權基因共表達分析。最終表達模式相近的轉錄本被分到相同的模塊，每個模塊被標記為單獨的顏色，共獲得42個模塊（圖4：A）。所有TS lncRNA都位于與其所在組織高度相關的模塊（P＜0.05）。選取與每個TS lncRNA相關程度最高的10個mRNA，并對其進行功能富集分析。結果顯示，每個組織中TS lncRNA的功能都與該組織相關，肌肉中TS lncRNA功能預測中顯著的GO條目有“肌肉組織發育”“橫紋肌細胞分化”等（圖4：B），而脾臟中對應的GO條目有“淋巴細胞激活”“免疫系統過程的調節”等（圖4：C）。

圖4 通過反式作用的TS lncRNA功能預測Fig. 4 Functional prediction of TS lncRNAs by trans-action

2.4 調控IMF的lncRNA-mRNA篩選

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

MSTRG.7256.5位于第12條染色體，具有2個外顯子，全長300 nt，其對應的最小自由能結構和質心二級結構如圖所示（圖5：E，F）。

圖5 與IMF相關的lncRNA篩選Fig. 5 Identification of lncRNAs responsible for IMF

2.5 MSTRG.7256.5功能的驗證

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

圖6 MSTRG.7256.5作用的RT-qPCR驗證Fig. 6 Correlation validation of MSTRG.7256.5 by RT-qPCR

油紅O染色實驗結果顯示，MSTRG.7256.5干擾后的細胞內脂滴數量高于對照組（圖7：A）。過表達了MSTRG.7256.5的細胞脂滴數量則更低（圖7：B）。這些結果表明，MSTRG.7256.5可能通過PPAR途徑調控家豬的IMF含量。

圖7 MSTRG.7256.5干擾（A）與過表達后（B）的細胞內脂滴變化（油紅O染色實驗）Fig. 7 Change of lipid droplets after RNA interference（A） and overexpression （B）（oil red O staining experiment）

3 討論

由于其營養價值和經濟價值，家豬已經成為許多研究中的重要生物模型。lncRNA在許多生物過程中發揮著關鍵作用（Zhuet al.，2019），然而lncRNA-mRNA的相互作用，及其對組織特異性表達和IMF沉積影響的研究還有待深入。在本研究中，利用公共數據庫中的RNA-seq數據，綜合分析了lncRNAs的組織特異性表達模式以及對豬IMF含量的影響，為未來相關養殖業的研究提供了數據和參考。

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

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

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

綜上所述，認為MSTRG.7256.5可能通過PPAR途徑負向調節家豬的脂肪沉積。然而，由于此步驟目的在于驗證MSTRG.7256.5能否對其靶基因和脂肪沉積產生影響，因此僅選用了PK-15細胞用于實驗，而PK-15細胞較低的脂肪含量可能導致實驗結果不顯著?？傊?，MSTRG.7256.5對其靶基因具體的作用機制和對肌肉間IMF含量的影響還須通過實驗進一步探究。