張小偉，何軍敏，徐新明，田月珍，關(guān)鳴軒，張 漫，付雪峰，哈尼克孜?吐拉甫，黃錫霞，田可川（1.新疆農(nóng)業(yè)大學(xué)動(dòng)物科學(xué)學(xué)院，新疆烏魯木齊 83005；.新疆畜牧科學(xué)院畜牧研究所，新疆烏魯木齊 830011；3.山東省農(nóng)業(yè)科學(xué)院畜牧獸醫(yī)研究所，山東濟(jì)南 50100）

羊毛品質(zhì)由密度、細(xì)度、強(qiáng)度、長度等多個(gè)數(shù)量性狀共同影響，這些數(shù)量性狀受多個(gè)基因或者多個(gè)基因座共同調(diào)控，且性狀間存在著十分復(fù)雜的相互作用。羊毛的品質(zhì)和產(chǎn)量是影響羊毛生產(chǎn)經(jīng)濟(jì)效益的重要因素，要提高動(dòng)物毛發(fā)的產(chǎn)量和品質(zhì)，首先應(yīng)了解毛囊的發(fā)生發(fā)育，闡明其發(fā)育的分子調(diào)控機(jī)制，解析與毛囊發(fā)育相關(guān)基因的功能。過去人們多對(duì)單個(gè)基因的功能作用進(jìn)行研究，近年來，毛囊發(fā)育和毛發(fā)生長的分子遺傳學(xué)研究進(jìn)展迅速，許多調(diào)控毛囊形態(tài)發(fā)生和毛發(fā)周期的信號(hào)通路以及基因相繼被發(fā)現(xiàn)和鑒定。生物體是一個(gè)復(fù)雜的系統(tǒng)，各基因間存在相互調(diào)控的關(guān)系，進(jìn)而形成一個(gè)復(fù)雜的調(diào)控網(wǎng)絡(luò)。

近年來的研究表明，ncRNA 參與了多種生物過程，lncRNA 和miRNA 在動(dòng)物的生長發(fā)育中起著重要作用，在各種組織和器官中發(fā)揮著多種生物學(xué)功能。已報(bào)道了幾種與毛囊生長相關(guān)的lncRNA 和miRNA，如lnc-000133 參與山羊毛囊生長，Qian 等發(fā)現(xiàn)lncRNA參與絨山羊的毛囊生長。Gao 等通過研究探索候選miRNA 與毛囊發(fā)育特征之間的相關(guān)性，發(fā)現(xiàn)miR-143會(huì)調(diào)節(jié)湖羊毛囊生長。Zhai 等通過鑒定毛囊發(fā)育中miRNA-21 的靶基因，發(fā)現(xiàn)其可為調(diào)節(jié)細(xì)胞發(fā)育提供新的研究途徑。何玉龍等通過比較miRNA 在小尾寒羊和灘羊皮膚毛囊中的表達(dá)差異，對(duì)差異表達(dá)miRNA 靶基因進(jìn)行富集分析發(fā)現(xiàn)其主要參與代謝過程、催化活性、細(xì)胞進(jìn)程和細(xì)胞組分等過程。毛囊的生長發(fā)育和周期性變化是受多種因子調(diào)控的復(fù)雜分子調(diào)控過程，miRNA、lncRNA 和mRNA 之間相互作用構(gòu)成了一個(gè)復(fù)雜的調(diào)控過程，越來越多的研究報(bào)道了miRNA和lncRNA在毛囊發(fā)育及周期性生長中的作用，但是對(duì)毛囊生長發(fā)育調(diào)控網(wǎng)絡(luò)的研究還較少。因此，本研究以蘇博美利奴羊?yàn)閷?shí)驗(yàn)動(dòng)物，通過生物信息學(xué)方法挖掘蘇博美利奴羊皮膚組織樣本轉(zhuǎn)錄組數(shù)據(jù)，分析lncRNA、miRNA 與mRNA 之間存在的相互調(diào)控關(guān)系，構(gòu)建lncRNA-miRNA-mRNA 分子調(diào)控網(wǎng)絡(luò)，揭示了與毛囊發(fā)育相關(guān)的ncRNA 網(wǎng)絡(luò)調(diào)控機(jī)制，為后續(xù)毛囊發(fā)育的研究提供了一定的參考。

1 材料與方法

1.1 實(shí)驗(yàn)材料 課題組前期于新疆科創(chuàng)畜牧繁育中心采集蘇博美利奴羊胎齡65 d（G1）、85 d（G2）、105 d（G3）、135 d（G4）及出生后7 d（G5）、30 d（G6）6 個(gè)時(shí)期羔羊肩胛部皮膚組織，每個(gè)時(shí)期各采集3 個(gè)樣本，6 個(gè)時(shí)期共有18 個(gè)樣本。提取RNA 構(gòu)建cDNA 文庫，并利用Illumina Hiseq 2000 進(jìn) 行RNA-seq 測序。采用Edge R 對(duì)蘇博美利奴羊不同時(shí)期轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行差異分析（fold change ≥ 2，q value ≤ 0.05），得到不同時(shí)期的差異表達(dá)lncRNA（DELs）、差異表達(dá)miRNA（DEMs）及差異表達(dá)基因（DEGs）。

1.2 靶向關(guān)系預(yù)測 結(jié)合課題組前期lncRNA、miRNA和mRNA 測序和篩選毛囊發(fā)育過程中差異表達(dá)的mRNA及ncRNA結(jié)果，通過生物信息學(xué)方法分析蘇博美利奴羊毛囊發(fā)育中l(wèi)ncRNA、miRNA 與mRNA間存在的調(diào)控作用關(guān)系。使用DIANA Tools 的LncBase Predicted v.2 網(wǎng)站預(yù)測miRNA 靶向調(diào)控的lncRNA，運(yùn)用Venny 2.1 軟件與差異lncRNA 取交集，獲得miRNAlncRNA 靶向關(guān)系。運(yùn)用TargetScan和PITA 網(wǎng)站在線預(yù)測miRNA 靶向調(diào)控的mRNA，運(yùn)用Venny 2.1軟件與差異mRNA 取交集，獲得miRNA-mRNA 靶向關(guān)系。

1.3 lncRNA-miRNA-mRNA 調(diào)控網(wǎng)絡(luò)的構(gòu)建 通過以上生物信息學(xué)方法預(yù)測獲得miRNA-lncRNA 及miRNAmRNA 靶向調(diào)控關(guān)系，考慮到mRNA 大多通過miRNA的負(fù)調(diào)控而發(fā)揮重要生物學(xué)作用，所以本研究特別關(guān)注了在表達(dá)水平上與miRNA 表達(dá)量呈負(fù)相關(guān)的靶基因，即上調(diào)的miRNA 對(duì)應(yīng)下調(diào)的mRNA，而下調(diào)的miRNA對(duì)應(yīng)上調(diào)的mRNA。運(yùn)用Cytoscape 軟件（版本3.8.0）構(gòu)建lncRNA-miRNA-mRNA 調(diào)控關(guān)系網(wǎng)絡(luò)，網(wǎng)絡(luò)由節(jié)點(diǎn)（node）和邊（edge）組成，節(jié)點(diǎn)代表DELs、DEMs 和DEGs，連線代表節(jié)點(diǎn)之間存在調(diào)控關(guān)系。

1.4 功能分析 運(yùn)用DAVID v6.8 對(duì)調(diào)控網(wǎng)絡(luò)中的靶基因進(jìn)行GO 富集分析（Gene Ontology 基因本體論，GO），主要從生物學(xué)過程（Biological Process，BP）、細(xì)胞組分（Cellular Component，CC）以及分子功能（Molecular Function，MF）3 個(gè)方面探索靶基因在蘇博美利奴羊毛囊發(fā)育過程中所起的作用。進(jìn)一步對(duì)靶基因進(jìn)行KEGG（Kyoto Encyclopedia of Genes and Genomes 京都基因與基因組百科全書，KEGG）通路分析，以研究其生物學(xué)功能。并利用R 語言（版本4.0.0）ggplot2 軟件包對(duì)富集結(jié)果進(jìn)行繪圖。

2 結(jié)果與分析

2.1 靶向關(guān)系預(yù)測 通過生物信息學(xué)方法預(yù)測獲得miRNAlncRNA 及miRNA-mRNA 靶向調(diào)控關(guān)系，分析發(fā)現(xiàn)在G3vs.G2 組TCONS_00008348 與miR-150、miR-154b-3p、miR-485-3p、miR-539-3p 及miR-655-5p 有靶向關(guān)系，在G6vs.G1 組TCONS_00008348、TCONS_00001791 和TCONS_00022645 與miR-10b、miR-26b 等有靶向關(guān)系，如表1 所示。

表1 不同時(shí)期miRNA-lncRNA 靶向預(yù)測

2.2 lncRNA-miRNA-mRNA 調(diào)控網(wǎng)絡(luò) 根據(jù)lncRNAmiRNA 及miRNA-mRNA 之間的靶向作用關(guān)系，運(yùn)用Cytoscape 軟件構(gòu)建lncRNA-miRNA-mRNA 調(diào)控關(guān)系網(wǎng)絡(luò)，考慮到miRNA 大多通過負(fù)調(diào)控mRNA 而發(fā)揮重要生物學(xué)作用，所以本研究特別關(guān)注了在表達(dá)水平上與miRNA 表達(dá)量呈負(fù)相關(guān)的靶基因，即上調(diào)的miRNA對(duì)應(yīng)下調(diào)的mRNA，而下調(diào)的結(jié)果如圖1 所示，不同時(shí)期mRNA 與lncRNA 及miRNA 具有復(fù)雜的網(wǎng)絡(luò)連接。在G6vs.G1 組中，由于時(shí)期間隔較長，其調(diào)控網(wǎng)絡(luò)也較復(fù)雜。

圖1 不同時(shí)期lncRNA-miRNA-mRNA 調(diào)控網(wǎng)絡(luò)

2.3 功能分析 通過GO 富集分析發(fā)現(xiàn)，G2vs.G1 組靶基因共富集到48 個(gè)GO Term 上，G3vs.G2 組靶基因共富集到59 個(gè)GO Term 上，G4vs.G3 組靶基因共富集到41 個(gè)GO Term 上，G5vs.G4 組靶基因共富集到20個(gè)GO Term 上，G6vs.G5 組靶基因共富集到17 個(gè)GO Term 上，G6vs.G1 組靶基因共富集到97 個(gè)GO Term 上，差異的靶基因主要參與細(xì)胞增殖負(fù)調(diào)控、RNA 聚合酶II 啟動(dòng)子轉(zhuǎn)錄的正調(diào)控、基因表達(dá)的正調(diào)控、凋亡過程正調(diào)控等生物學(xué)過程，參與細(xì)胞外空間、細(xì)胞外泌體、細(xì)胞質(zhì)、轉(zhuǎn)錄因子復(fù)合體、質(zhì)膜的整體成分等細(xì)胞組分；參與轉(zhuǎn)錄因子活性、結(jié)構(gòu)分子活性、細(xì)胞外基質(zhì)結(jié)構(gòu)組成、染色質(zhì)結(jié)合及鋅鈣離子結(jié)合等分子功能，并且分析發(fā)現(xiàn)G2vs.G1 組GO:0031069 為毛囊形態(tài)發(fā)生Term，富集到該Term 的基因?yàn)楹停桓鹘MGO富集分析Top30 如圖2 所示。

圖2 不同時(shí)期靶基因GO 富集分析

通過KEGG 富集分析發(fā)現(xiàn)，G2vs.G1 組的靶基因共富集到7 個(gè)通路上，G3vs.G2 組的靶基因共富集到11 個(gè)通路上，G4vs.G3 組的靶基因共富集到11 個(gè)通路上，G5vs.G4 組的靶基因共富集到7 個(gè)通路上，G6vs.G5 組的靶基因共富集到6 個(gè)通路上，G6vs.G1 組靶基因共富集到12 個(gè)通路上，如圖3 所示，差異的靶基因主要參與TGF-beta 信號(hào)通路、癌癥途徑、MAPK 信號(hào)通路、PI3K-Akt 信號(hào)通路、細(xì)胞因子與細(xì)胞因子受體的相互作用等多種途徑。

圖3 不同時(shí)期靶基因KEGG 通路分析

3 討論

本研究采用生物信息學(xué)方法，從蘇博美利奴羊皮膚組織樣本轉(zhuǎn)錄組數(shù)據(jù)篩選出差異表達(dá)的基因，構(gòu)建了差異基因的lncRNA-miRNA-mRNA 調(diào)控網(wǎng)絡(luò)，解析毛囊發(fā)生發(fā)育的分子調(diào)控機(jī)制，從轉(zhuǎn)錄組學(xué)層面對(duì)毛囊生長發(fā)育進(jìn)行探討。

本實(shí)驗(yàn)通過靶基因預(yù)測發(fā)現(xiàn)差異lncRNA 中TCONS_00008348、TCONS_00022645 及TCONS_00001791與miR-150、miR-154b-3p 等miRNA 有著靶向調(diào)控的關(guān)系，通過GO 富集分析發(fā)現(xiàn)差異的靶基因主要參與細(xì)胞增殖負(fù)調(diào)控、基因表達(dá)的正調(diào)控、凋亡過程正調(diào)控、細(xì)胞外空間、細(xì)胞外泌體、轉(zhuǎn)錄因子復(fù)合體等生物過程，通過KEGG 富集分析發(fā)現(xiàn)差異的靶基因主要參與TGF-beta 信號(hào)通路、MAPK 信號(hào)通路、PI3K-Akt 信號(hào)通路等多種途徑。且在G2vs.G1 組中靶基因和參與毛囊形態(tài)發(fā)生（GO:0031069）的生物學(xué)過程，G2vs.G1 組為胚胎期65 d 與85 d 比較組，此前課題組阿布來提? 蘇來曼等發(fā)現(xiàn)在蘇博美利奴羊胎齡65 d 時(shí)，胎兒表皮結(jié)構(gòu)完整，各部位均出現(xiàn)初級(jí)毛囊，在胎齡85 d 時(shí)，初級(jí)毛囊發(fā)生形成毛芽，初級(jí)毛囊基底部可見次級(jí)毛囊的囊泡結(jié)構(gòu)，所以胚胎期65～85 d 為蘇博美利奴羊毛囊發(fā)育的關(guān)鍵時(shí)期。吳萍發(fā)現(xiàn)Notch信號(hào)通路的、等基因涉及控制毛囊發(fā)育和生長的分子信號(hào)途徑來改良毛用動(dòng)物的生產(chǎn)性能，因此基因在毛囊發(fā)育過程中起著關(guān)鍵作用。TGF-及其受體在毛囊中的表達(dá)具有部位和毛發(fā)生長周期的特異性，且有多種信號(hào)途徑參與TGF-對(duì)毛囊生長發(fā)育的調(diào)控。

lncRNA-miRNA-mRNA 調(diào)控機(jī)制在各種疾病中得到了廣泛的研究，而在綿羊毛囊發(fā)育調(diào)控網(wǎng)絡(luò)中的研究還相對(duì)較少。本實(shí)驗(yàn)分析發(fā)現(xiàn)，在G2vs.G1 組調(diào)控網(wǎng)絡(luò)中參與毛囊形態(tài)發(fā)生生物學(xué)過程的基因和為miR-133、miR-433-5p 的靶基因，因此構(gòu)成miR-133-、miR-133-及miR-133-、miR-433-5p-的調(diào)控關(guān)系。這與Zhao 等分析胚胎90 d、胚胎120 d 及出生日的敖漢細(xì)毛羊ncRNA 和mRNA 表達(dá)轉(zhuǎn)錄數(shù)據(jù)所構(gòu)建的與毛囊生長調(diào)控網(wǎng)絡(luò)關(guān)系相類似，其通過調(diào)控網(wǎng)絡(luò)發(fā)現(xiàn)miR-21 與lncRNA 作為競爭的內(nèi)源性RNA 參與毛囊的發(fā)育。Wang等對(duì)毛囊循環(huán)周期中生長期和靜止期的ncRNAs進(jìn)行分析，發(fā)現(xiàn)調(diào)控網(wǎng)絡(luò)中miR-221-5p-lnc_000679-WNT3、miR-34a-lnc_000181-GATA3 和miR-214-3plnc_000344-SMAD3 在絨山羊毛囊生物學(xué)中發(fā)揮著重要作用。Han 等通過對(duì)敖漢細(xì)毛羊的關(guān)鍵靶基因進(jìn)行分析，獲得了7 個(gè)可能調(diào)控的候選lncRNA，并構(gòu)建了一個(gè)lncRNA-mRNA 共表達(dá)網(wǎng)絡(luò)，推測這些候選lncRNAs 可能通過調(diào)節(jié)靶基因的表達(dá)發(fā)揮作用。因此，ncRNA 調(diào)控網(wǎng)絡(luò)在調(diào)節(jié)毛囊生長中發(fā)揮作用，但其在毛囊發(fā)育中的具體功能仍需進(jìn)一步研究。

4 結(jié)論

本研究通過構(gòu)建蘇博美利奴羊毛囊發(fā)育相關(guān)lncRNA-miRNA-mRNA 調(diào)控網(wǎng)絡(luò)，揭示了與毛囊發(fā)育相關(guān)的ncRNA 網(wǎng)絡(luò)調(diào)控機(jī)制，通過富集分析發(fā)現(xiàn)靶基因參與細(xì)胞增殖的負(fù)調(diào)控、毛囊形態(tài)發(fā)生、TGF-beta信號(hào)通路等多種途徑，并在網(wǎng)絡(luò)中篩選出與毛囊發(fā)育的關(guān)鍵基因和參與毛囊形態(tài)發(fā)生的生物學(xué)過程，為后續(xù)毛囊發(fā)育的研究提供了一定的參考。