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

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

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

1 材料與方法

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

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

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

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

2 結果與分析

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

表1 不同時期miRNA-lncRNA 靶向預測

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

圖1 不同時期lncRNA-miRNA-mRNA 調控網絡

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

圖2 不同時期靶基因GO 富集分析

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

圖3 不同時期靶基因KEGG 通路分析

3 討論

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

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

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

4 結論

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