綿羊組織特異性長鏈非編碼RNA的鑒定與功能分析

趙 琛，劉曉懿，王喜宏

( 西北農(nóng)林科技大學(xué) 動物科技學(xué)院，陜西 楊凌 712100)

長鏈非編碼RNA(Long non-coding RNA, lncRNA)是指長度超過200個核苷酸，且不具備蛋白編碼能力的RNA[1]。lncRNA在染色質(zhì)調(diào)控、轉(zhuǎn)錄調(diào)控和轉(zhuǎn)錄后調(diào)控中具有重要作用[2]，主要通過招募染色質(zhì)重塑復(fù)合物，介導(dǎo)染色體間的相互作用，作為轉(zhuǎn)錄因子“向?qū)А被颉罢T餌”，充當(dāng)miRNA的“海綿”等方式發(fā)揮抑制或激活作用[3]。近年來高通量測序方法的迅速發(fā)展使得對家畜lncRNA的研究增加，Bakhtiarizadeh等[4]分析綿羊8種組織的RNA-seq數(shù)據(jù)并報告了第一個基因間lncRNA(lincRNA)數(shù)據(jù)集(n=325)。Bush等[5]根據(jù)綿羊3.1版本參考基因組(Oar_v3.1)在11個綿羊組織中鑒定出11646個新的lncRNA。但是與其他哺乳動物(例如人，小鼠，牛和豬)相比，綿羊基因組中注釋的lncRNA很少。NONCODE數(shù)據(jù)庫是一個比較全面的ncRNA相關(guān)注釋的數(shù)據(jù)庫，包含人類，小鼠，大鼠，牛，雞和豬等17個物種的ncRNA信息[6]，但是該數(shù)據(jù)庫中缺少綿羊的lncRNA數(shù)據(jù)集。本研究以綿羊為研究對象，以高通量測序為手段，篩選并鑒定新型的lncRNA及組織性lncRNA，旨在發(fā)現(xiàn)對家畜重要經(jīng)濟性狀有調(diào)控作用的lncRNA，豐富家畜lncRNA數(shù)據(jù)庫，為lncRNA在育種工作中的應(yīng)用提供借鑒和參考。

1 材料與方法

1.1 試驗數(shù)據(jù)

綿羊(Ovisaries)RNA-seq測序數(shù)據(jù)包括特克賽爾羊♂×哈薩克羊♀雜交綿羊F1代個體、哈薩克羊、特克塞爾羊♂×蘇格蘭黑臉羊♀的雜交F1代個體的192個組織樣本。由肌肉(半腱肌、半膜肌、背最長肌)、皮膚(腹側(cè)皮膚、體側(cè)皮膚)、脂肪(尾脂、皮下脂肪)、心臟、肺、肝臟、腎臟、小腸、脾臟構(gòu)成，包含新生、出生后1周、8周、3個月、5個月、6個月與成年的數(shù)據(jù)。原始數(shù)據(jù)下載自NCBI SRA數(shù)據(jù)庫，BioProject號為PRJNA485657[7]。

1.2 序列比對與轉(zhuǎn)錄本組裝

使用Trimmomatic(版本0.36)[8]對每個樣本的轉(zhuǎn)錄組測序原始數(shù)據(jù)進行質(zhì)控，具體參數(shù)如下：LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:40。去除接頭序列和過濾低質(zhì)量的讀段后，通過STAR軟件(版本2.5.1)[9]將樣本讀段比對到綿羊參考基因組(Oar_v4.0)[10]。利用StringTie軟件(版本1.3.4)[11]將比對結(jié)果進行樣本的組裝，并通過String Tie --merge將多個樣本組裝合并，得到綿羊的轉(zhuǎn)錄本集合。

1.3 lncRNA的鑒定

lncRNA的鑒定步驟：(1)去除長度小于200 bp，外顯子數(shù)目小于2的轉(zhuǎn)錄本[12]；(2)篩去蛋白編碼基因；(3)利用CPC2軟件[13]，CNCI軟件[14]預(yù)測候選lncRNA的編碼能力，去除評分不合格的轉(zhuǎn)錄本；(4)使用blastx與UniProt數(shù)據(jù)庫(版本2016_04)[15]和Pfam數(shù)據(jù)庫(版本29)[16]中的已知蛋白質(zhì)進行比對，認(rèn)為E值≤10-3，比對長度≥10個氨基酸和一致性≥95%的轉(zhuǎn)錄本為“潛在的編碼基因”，將此類轉(zhuǎn)錄本去除[17]。保留篩選下來的轉(zhuǎn)錄本作為綿羊候選lncRNA轉(zhuǎn)錄本。統(tǒng)計lncRNA與編碼基因的外顯子數(shù)目、GC含量、表達(dá)量等特征。

1.4 表達(dá)量評估

采用StringTie軟件(版本1.3.4)[11]計算lncRNA在每一個組織中的表達(dá)量，使用每千堿基外顯子百萬片段數(shù)(fragments per kilobase of exons per million fragments mapped，F(xiàn)PKM)來評估表達(dá)水平，并保留至少在任意一個組織中有20%以上的樣本FPKM表達(dá)量>1的轉(zhuǎn)錄本[7]，作為綿羊lncRNA的組織表達(dá)譜。

1.5 組織特異性分析

選取每種組織的表達(dá)量的平均值代表該組織，對候選lncRNA轉(zhuǎn)錄本計算τ組織特異性指數(shù)。τ指數(shù)是基因特異或廣泛表達(dá)的指標(biāo)，取值范圍為0-1，其中1表示僅在一個組織中特異表達(dá)，而0表示所有組織的廣泛表達(dá)[18]。給定基因的τ指數(shù)可以用以下公式計算:

式中n是被研究組織的數(shù)量，xi是給定組織的表達(dá)譜。基于τ計算組織特異性分?jǐn)?shù)，認(rèn)為τ≥0.8且在所有組織中表達(dá)量排名前兩位的為具有組織特異性的lncRNA，分別計算肌肉、脂肪、皮膚特異性的lncRNA。

1.6 靶基因功能預(yù)測

lncRNA的順式作用主要通過影響鄰近蛋白編碼基因的表達(dá)發(fā)揮調(diào)控功能[19]。lncRNA的反式作用主要通過識別與候選lncRNA共表達(dá)的mRNA，調(diào)控遠(yuǎn)距離mRNA的轉(zhuǎn)錄激活與表達(dá)[20]。(1)順式作用預(yù)測：皮爾森相關(guān)系數(shù)是用來反映兩個變量線性相關(guān)程度的統(tǒng)計量。提取位于lncRNA上游10 kb，下游100 kb內(nèi)的編碼基因作為靶基因[21]。(2)反式作用預(yù)測：采用皮爾森相關(guān)系數(shù)法分析樣本間lncRNA與蛋白編碼基因的相關(guān)性，取相關(guān)性較高的蛋白編碼基因作為靶基因(r≥0.90或r≤-0.90, Pvalue<0.01)[22]。

1.7 功能富集分析

使用R軟件中的ClusterProfiler軟件包[23]對基因進行KEGG通路富集分析，設(shè)置adjust.p<0.05作為顯著性閾值。使用Metascape[24]對候選靶基因進行基因本體論(gene ontology, GO)富集分析，設(shè)置Pvalue<0.05作為顯著性閾值。

1.8 lncRNA-mRNA互作網(wǎng)絡(luò)的構(gòu)建

利用STRING數(shù)據(jù)庫(版本11.0)[22]分別對肌肉特異性lncRNA、皮膚特異性lncRNA、脂肪特異性lncRNA反式作用的靶基因數(shù)據(jù)集進行蛋白質(zhì)相互作用(Protein-Protein Interaction, PPI)網(wǎng)絡(luò)預(yù)測。整合PPI網(wǎng)絡(luò)與共表達(dá)的lncRNA-mRNA來構(gòu)建最終的互作關(guān)系。利用Cytoscape軟件(版本3.7)[25]中的ClusterONE插件(版本1.0)[26]檢測網(wǎng)絡(luò)互作中的顯著模塊。設(shè)置Pvalue≤0.05，表達(dá)模塊中最少的基因數(shù)量為3。使用Cytoscape軟件[25]對lncRNA-mRNA互作網(wǎng)絡(luò)進行可視化。

1.9 數(shù)量性狀位點(QTL)分析

綿羊QTL數(shù)據(jù)(Oar_v4.0)來自AnimalQTLdb[27]。通過bedtools intersect[28]比較QTL和lncRNA在全基因組中的位置，預(yù)測lncRNA的功能。起始位置和結(jié)束位置都在QTL區(qū)域內(nèi)的lncRNA被認(rèn)為成功注釋在QTL中[29]。

2 結(jié)果與分析

2.1 lncRNA的預(yù)測結(jié)果與組織特異性分析

利用軟件CNCI、CPC2與蛋白質(zhì)數(shù)據(jù)庫Pfam和UniProt分別預(yù)測出94 721條、75 580條、58 136條和94 987條lncRNA，四者的交集為36 340條lncRNA，認(rèn)為其為高可信度的綿羊lncRNA(圖1A)。對lncRNA和蛋白編碼基因的基本特征進行比較分析，lncRNA的GC含量為47.17%，低于編碼基因的GC含量48.79%(圖1B)；lncRNA外顯子平均數(shù)目為3.45個，低于蛋白編碼基因的12.8個；但外顯子個數(shù)在2～10區(qū)間中l(wèi)ncRNA為35 423個，多于編碼基因的23 466個(圖1C)；lncRNA和蛋白編碼基因的平均表達(dá)值分別為0.38和4.49，lncRNA的平均表達(dá)量低于蛋白編碼基因(圖1D)。

圖1 綿羊lncRNA與蛋白編碼基因的基本特征比較

基于RNA-Seq數(shù)據(jù)的主成分分析結(jié)果表明(圖2)，相似類型的組織樣本更傾向于聚集在一起，如半腱肌、半膜肌、背最長肌三種肌肉組織聚集，體側(cè)皮膚、腹側(cè)皮膚兩種皮膚組織聚集，尾脂和皮下脂肪兩種脂肪組織聚集，呈現(xiàn)出明顯的組織特異性。試驗共獲得肌肉特異性lncRNA 213個，皮膚特異性lncRNA 467個，脂肪特異性lncRNA 295個。

圖2 組織樣本的主成分分析

2.2 肌肉特異性lncRNA

基于綿羊的213個肌肉特異性lncRNA預(yù)測到了255個上下游基因。為了評估上下游基因的調(diào)控能力，從AnimalTFDB3.0數(shù)據(jù)庫[30]中提取了綿羊轉(zhuǎn)錄因子集合，將其與順式作用靶基因進行比較，發(fā)現(xiàn)有14個靶基因為綿羊轉(zhuǎn)錄因子。對肌肉特異性lncRNA上下游靶基因進行GO分析，發(fā)現(xiàn)顯著富集在肌肉系統(tǒng)過程(19，GO:0003012)(P<0.05)和肌肉結(jié)構(gòu)發(fā)育(21，GO:0007517)(P<0.05)。Six1為綿羊轉(zhuǎn)錄因子，是MSTRG.71513.3(7:69708802-69711134)的鄰近基因，對骨骼肌發(fā)育和肌纖維類型轉(zhuǎn)化具有調(diào)控作用[31]。肌肉特異性lncRNA共預(yù)測出63個反式作用的靶基因，GO富集注釋結(jié)果顯示，可注釋到7個功能條目，顯著富集是肌肉系統(tǒng)發(fā)育(16，GO:0003012)(P<0.05)，肌肉器官發(fā)育(11，GO:0007517)(P<0.05)，鈣介導(dǎo)信號(7，GO:0019722)(P<0.05)，肌動蛋白絲運動的調(diào)節(jié)(3，GO:1903115)(P<0.05)。括號內(nèi)包含注釋到該條目的基因數(shù)目與GO編號(圖3A)。

圖3 肌肉與皮膚特異性lncRNA反式作用靶基因的GO富集分析

2.3 皮膚特異性lncRNA

基于綿羊的467個皮膚特異性lncRNA預(yù)測到了1006個上下游基因，有44個為綿羊轉(zhuǎn)錄因子。對靶基因進行GO分析，共富集到11個功能條目，其中最顯著富集的為皮膚發(fā)育(28，GO:0043588)(P<0.05)，包括表皮生長因子家族基因AREG，骨形態(tài)發(fā)生蛋白受體1A型(BMPR1A)，角蛋白基因(KRT)，細(xì)胞間粘附分子(ICM)等。AREG位于MSTRG.68599.1(6:88881341-88903389)鄰近，KRT(KRT9，KRT15，KRT17)位于MSTRG.12764.8(11:41165824-41274496)鄰近，ICM(ICAM1，ICAM4，ICAM5)位于MSTRG.62954.1(5:12366881-12385388)鄰近。皮膚特異性lncRNA共預(yù)測出334個反式作用的靶基因。GO富集注釋結(jié)果顯示(圖3B)，lncRNA反式作用的靶基因可注釋到32個功能條目，顯著富集的前7位分別是角質(zhì)化(28，GO:0070268)(P<0.05)，蛻皮周期(14，GO:0042303)(P<0.05)，牙齒發(fā)生(11，GO:0042476)(P<0.05)，表皮發(fā)育調(diào)節(jié)(10，GO:0045682)(P<0.05)，長鏈脂肪?；o酶A代謝(7，GO:0035336)(P<0.05)，蛋白交聯(lián)(7，GO:0018149)(P<0.05)，皮膚屏障的建立(6，GO:0061436)(P<0.05)。進一步選取更顯著的lncRNA-mRNA對(r≥0.95或r≤-0.95)進行共表達(dá)分析，發(fā)現(xiàn)可分為4個表達(dá)模塊，黃色圓圈表示ClusterONE[26]識別的顯著表達(dá)模塊，紅色圓圈表示組織特異性lncRNA，藍(lán)色圓圈表示lncRNA反式作用靶基因。若僅關(guān)注包含lncRNA的表達(dá)模塊，最顯著的為MSTRG.39053.2(2:228846018-228866160)所在模塊(P<0.05)(圖4)，共10個基因，包含TGM1、IL36RN、PLA2G4E、THEM5、CERS3、PNPLA1、ALOX12B、KRT28、DSG3編碼基因，對該模塊富集分析發(fā)現(xiàn)這些基因顯著富集在角質(zhì)化(4，GO:0070268)(P<0.05)。

圖4 皮膚lncRNA-mRNA共表達(dá)網(wǎng)絡(luò)圖

2.4 脂肪特異性lncRNA

基于綿羊的295個脂肪特異性lncRNA預(yù)測到了787個上下游基因，有43個為綿羊轉(zhuǎn)錄因子。對上下游基因進行GO分析，顯著富集到脂質(zhì)定位(22，GO:0010876)(P<0.05)。轉(zhuǎn)錄因子HOXC基因HOXC11、HOXC12、HOXC13位于LNC12062(LOC105611004，3：132242058-132246141)鄰近。脂肪特異性lncRNA共預(yù)測出79個反式作用的靶基因。KEGG數(shù)據(jù)庫注釋結(jié)果顯示，注釋到8條通路(P<0.05)，注釋基因數(shù)最多的前5位分別為甘油脂代謝(5)(P<0.05)、脂肪酸伸長率(3)(P<0.05)、不飽和脂肪酸的生物合成(3)(P<0.05)、纈氨酸，亮氨酸和異亮氨酸的降解(3)(P<0.05)、脂肪的消化吸收(3)(P<0.05)。GO富集注釋結(jié)果顯示，注釋到11個功能條目，顯著富集的是脂質(zhì)代謝(14，R-HSA-556833)(P<0.05)，DAG和TAG的酰基鏈重塑(14，R-HSA-556833)(P<0.05)，跨化學(xué)突觸的傳遞(7，R-HSA-112315)(P<0.05)，膽汁鹽和有機酸，金屬離子和胺化合物的轉(zhuǎn)運(4，R-HSA-425366)(P<0.05)。lncRNA-mRNA共表達(dá)分析發(fā)現(xiàn)了1個顯著表達(dá)模塊(圖5)，包含8個基因，其中包含2個預(yù)測的lncRNA——MSTRG.57079.1(3:166886406-167025566)、MSTRG.68891.1(6:100780317-100820136)。對該模塊進行GO富集分析發(fā)現(xiàn)主要富集在?；o酶A代謝過程(4，GO:0006637)(P<0.05)，KEGG通路最顯著富集于甘油三酯的生物合成(P<0.05)，這兩個過程都與脂肪合成相關(guān)。括號內(nèi)的數(shù)字表示注釋到該條目的基因數(shù)。

圖5 脂肪lncRNA-mRNA共表達(dá)網(wǎng)絡(luò)圖

2.5 QTL分析

綿羊QTL數(shù)據(jù)庫[27]包含肌肉相關(guān)QTL 681個，脂肪相關(guān)QTL 231個。對組織特異性lncRNA與綿羊QTL數(shù)據(jù)庫的重合關(guān)系進行檢索，共發(fā)現(xiàn)14個lncRNA可能作為綿羊肌肉與脂肪相關(guān)性狀的候選lncRNA。存在6個肌肉特異性lncRNA與肌肉相關(guān)QTL區(qū)域重合，其中有3個位于胴體肌肉重量QTL區(qū)域，1個位于肌肉密度QTL區(qū)域，2個位于體重QTL區(qū)域(表1)。11個脂肪特異性lncRNA與脂肪相關(guān)QTL區(qū)域重合，其中有5個位于尾脂沉積區(qū)域，6個位于全部脂肪區(qū)域(表1)。MSTRG.77930.1(9:93629571-93634360)與尾脂沉積QTL重合，其順式作用基因之一為與脂質(zhì)合成相關(guān)的DEPTOR，可通過抑制胰島素信號mTORC1介導(dǎo)的反饋抑制作用來促脂肪合成[32]。

表1 位于QTL區(qū)域的組織特異性lncRNA

3 討 論

本研究共鑒定36 340個綿羊lncRNA，發(fā)現(xiàn)肌肉特異性lncRNA 213個，皮膚特異性lncRNA 467個，脂肪特異性lncRNA 295個，擴展了綿羊基因組lncRNA集合，為進一步分析提供了良好的基礎(chǔ)。與編碼基因相比，lncRNA的GC含量、表達(dá)豐度較低，但平均長度更長，并且組織特異性高。這些基本特征與前人報道的哺乳動物lncRNA相似[33]。

肌肉質(zhì)量是決定畜禽經(jīng)濟價值的重要因素之一，肌肉發(fā)育和生長過程受遺傳、環(huán)境和營養(yǎng)的影響。lncRNA已被證實通過順式或反式作用調(diào)節(jié)肌肉生長發(fā)育[34]，很可能是重要的潛在的肌肉生長調(diào)節(jié)因子。Cai等[35]發(fā)現(xiàn)位于雞Six1基因上游432 bp的lncRNA-Six1通過編碼一種微肽激活Six1基因，發(fā)揮順式調(diào)節(jié)作用，促進細(xì)胞增殖并參與肌肉生長。Mueller等[36]發(fā)現(xiàn)小鼠成肌轉(zhuǎn)錄因子MyoD基因的轉(zhuǎn)錄起始位點上游5 kb處編碼的lncRNA MUNC，可以促進MyoD在骨骼肌生成。MSTRG.46395.2的靶基因是肌球蛋白1(MYOM1)，位于MYOM1轉(zhuǎn)錄下游區(qū)域。MYOM1在橫紋肌的細(xì)胞分化中編碼一種高度特異性的蛋白質(zhì)標(biāo)記，并可能在肌原纖維的組裝或維持中起作用[37]。本研究發(fā)現(xiàn)MSTRG.46395.2在肌肉組織中特異表達(dá)，可能表明MSTRG.46395.2與MYOM之間可能存在調(diào)控關(guān)系，但是需要進一步研究以闡明確切的作用機理。Six1是參與骨骼肌發(fā)育的重要轉(zhuǎn)錄因子[38]，肌肉特異性lncRNA MSTRG.71513.3位于其下游7 kb。有研究表明Six1不僅可以調(diào)節(jié)灘羊的肌肉纖維類型，而且還可以調(diào)節(jié)肉的嫩度，可以作為選擇肉嫩度的候選基因[39]。MSTRG.71513.3在肌肉組織中特異表達(dá)，推測MSTRG.71513.3與Six1之間可能存在調(diào)控關(guān)系，調(diào)控肌肉發(fā)育。但是確切的作用機理需要進一步研究。

毛囊發(fā)育和周期性生長過程中受多種信號通路調(diào)控，主要有WNT通路、轉(zhuǎn)化生長因子(TGF)家族、絲裂原活化蛋白激酶(MAPK)家族、Notch信號通路、骨形成蛋白(BMP)信號通路等[40]。有研究發(fā)現(xiàn)綿羊胎兒皮膚初級毛囊中下調(diào)的lncRNA與角蛋白(KRT14和KRT15)、WNT信號通路(WNT16和SFRP1)等通路基因相互作用，進而影響表皮和羊毛發(fā)育[41]。本研究發(fā)現(xiàn)，皮膚特異性lncRNA的上下游基因與反式作用共表達(dá)基因主要與角質(zhì)化、毛囊發(fā)育相關(guān)。Notch信號通路和毛囊的形態(tài)發(fā)生具有直接的關(guān)系，可能參與了毛囊發(fā)育的早期環(huán)節(jié)，并在毛干分化過程中起重要作用[40]。DLL4(Notch配體)位于MSTRG.70607.8(7:33364804-33478104)鄰近，猜測該皮膚特異性lncRNA可能對毛囊發(fā)育具有調(diào)控作用。角蛋白是一種中間纖維形成蛋白，是復(fù)雜上皮細(xì)胞的主要成分[42]。有研究表明具有無效的KRT17基因的小鼠在出生后的第一周內(nèi)無法形成完整的皮毛[43]，MSTRG.12764.8鄰近基因包括多個KRT基因KRT9、KRT15、KRT17，猜測其對表皮角質(zhì)化及毛囊發(fā)育有重要作用，但具體作用機制待進一步分析。

lncRNA在脂質(zhì)代謝中的調(diào)節(jié)功能已經(jīng)在不同生物體中得到了證實，包括人類[44]，小鼠[45]，豬[46]，牛[47]等。有研究表明，皮下脂肪組織中HOXC簇的表達(dá)上調(diào)[48]。還有研究發(fā)現(xiàn)，HOXC11，HOXC12，HOXC13等發(fā)育基因與綿羊尾部脂肪沉積有關(guān)[49]。本研究發(fā)現(xiàn)LOC105611004分別位于HOXC11，HOXC12，HOXC13附近，猜測該lncRNA可能對綿羊尾脂沉積有調(diào)控作用。還發(fā)現(xiàn)MSTRG.77930.1順式調(diào)控的靶基因是DEPTOR。DEPTOR通過激活脂肪細(xì)胞中的過氧化物酶體增殖物激活受體γ自主促進脂肪生成，是脂肪形成的新調(diào)節(jié)劑[50]。lncRNA-mRNA共表達(dá)分析發(fā)現(xiàn)了2個lncRNA(MSTRG.57079.1和MSTRG.68891.1)出現(xiàn)在與脂肪相關(guān)的表達(dá)模塊中，與?；o酶A代謝和甘油三酯的生物合成功能相關(guān)，有研究表明骨骼肌脂肪儲存的增加通常伴隨脂質(zhì)代謝產(chǎn)物(例如甘油二酯和長鏈酰基輔酶A)的逐步積累[51]，推測這兩個lncRNA可能對脂肪代謝具有調(diào)控作用。

4 結(jié) 論

本研究通過分析綿羊不同轉(zhuǎn)錄組的RNA-seq數(shù)據(jù)，獲得了新的綿羊lncRNA數(shù)據(jù)集，并篩選出了具有肌肉、皮膚、脂肪組織特異性的lncRNA?；陧樖阶饔煤头词阶饔脙蓚€方面，對lncRNA靶基因的功能進行分析，發(fā)現(xiàn)其可能在組織代謝、發(fā)育等多個功能中發(fā)揮調(diào)控作用，推測肌肉特異性lncRNA可能與肌肉生長功能相關(guān)，脂肪特異性lncRNA可能與脂肪合成相關(guān)，皮膚特異性lncRNA可能與表皮形成和毛囊發(fā)育相關(guān)。同時通過QTL篩選，獲得了與綿羊肌肉、皮膚與脂肪性狀相關(guān)的候選長鏈非編碼RNA。本研究為探究lncRNA在綿羊重要經(jīng)濟性狀的功能提供了新見解。