蘇博美利奴羊毛囊發(fā)育相關(guān)lncRNA與mRNA共表達(dá)網(wǎng)絡(luò)的構(gòu)建

陳華峰，田可川，黃錫霞，阿布來提·蘇來曼，何軍敏，田月珍，徐新明，付雪峰，趙冰茹，朱樺，哈尼克孜·吐拉甫

蘇博美利奴羊毛囊發(fā)育相關(guān)lncRNA與mRNA共表達(dá)網(wǎng)絡(luò)的構(gòu)建

陳華峰1,2，田可川2，黃錫霞1，阿布來提·蘇來曼1，何軍敏3，田月珍2，徐新明2，付雪峰2，趙冰茹4，朱樺1，哈尼克孜·吐拉甫2

（1新疆農(nóng)業(yè)大學(xué)動(dòng)物科學(xué)學(xué)院，烏魯木齊 830052；2新疆畜牧科學(xué)院畜牧研究所，烏魯木齊 830011；3甘肅農(nóng)業(yè)大學(xué)動(dòng)物科學(xué)技術(shù)學(xué)院，蘭州 730070；4中國農(nóng)業(yè)大學(xué)動(dòng)物科技學(xué)院，北京 100083）

【】隨著高通量測序技術(shù)的發(fā)展和完善，海量的轉(zhuǎn)錄組測序數(shù)據(jù)隨之涌現(xiàn)，基因網(wǎng)絡(luò)的方法也越來越多被用于研究中；細(xì)毛羊毛囊發(fā)育受多個(gè)基因及l(fā)ncRNA共同調(diào)控，單獨(dú)研究某一分子并不足以發(fā)現(xiàn)其調(diào)控機(jī)制，本研究旨在構(gòu)建毛囊發(fā)育相關(guān)lncRNA和mRNA的共表達(dá)網(wǎng)絡(luò)，挖掘細(xì)毛羊毛囊發(fā)育的潛在候選基因。以蘇博美利奴羊胎齡65 d（D65）、85 d（D85）、105 d（D105）、135 d（D135）的胎兒以及出生后7 d（A7）、30 d（A30）的羔羊肩胛部皮膚組織，每個(gè)時(shí)期有3個(gè)生物學(xué)重復(fù)，共18個(gè)樣本，進(jìn)行轉(zhuǎn)錄組測序以獲得6個(gè)不同發(fā)育時(shí)期的lncRNA與mRNA表達(dá)譜數(shù)據(jù)，篩選出相鄰時(shí)期的差異表達(dá)lncRNA與mRNA，運(yùn)用加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析（weighted gene co-expression network analysis，WGCNA）方法構(gòu)建共表達(dá)模塊，使用DAVID在線工具進(jìn)行GO（gene ontology）和 KEGG pathway富集分析以找到與毛囊發(fā)育相關(guān)的模塊，最后從目標(biāo)模塊篩選高連通度的lncRNA與mRNA使用Cytoscape軟件進(jìn)行網(wǎng)絡(luò)可視化。從表達(dá)譜數(shù)據(jù)中篩選得到9 070個(gè)差異表達(dá)的lncRNA與mRNA，運(yùn)用WGCNA方法得到11個(gè)模塊，使用DAVID富集分析發(fā)現(xiàn)honeydew1、paleturquoise、skyblue2模塊中的基因富集在皮膚發(fā)育、毛囊發(fā)育、毛囊形態(tài)發(fā)生、Wnt信號通路負(fù)調(diào)控、細(xì)胞粘附等生物過程，以及Wnt信號通路、TGF-β信號通路、Hedgehog信號通路、緊密連接、MAPK信號通路、ECM-受體相互作用等毛囊發(fā)育相關(guān)的信號通路，并篩選這3個(gè)模塊中連通度高的lncRNA和mRNA構(gòu)建了子網(wǎng)絡(luò)，得到包括、、、、、、、、在內(nèi)的毛囊發(fā)育相關(guān)基因，預(yù)測出ENSOART00000029117、TCONS_00489976、TCONS_00376759等24個(gè)lncRNA的可能靶基因。首次運(yùn)用WGCNA方法構(gòu)建了蘇博美利奴羊毛囊發(fā)育相關(guān)lncRNA與mRNA的共表達(dá)網(wǎng)絡(luò)，鑒別出3個(gè)毛囊發(fā)育相關(guān)的共表達(dá)模塊，發(fā)現(xiàn)多個(gè)與毛囊發(fā)育相關(guān)的潛在候選基因，并預(yù)測出24個(gè)lncRNA可能的靶基因。

蘇博美利奴羊；WGCNA；mRNA；lncRNA；毛囊發(fā)育

0 引言

【研究意義】隨著高通量測序技術(shù)的高速發(fā)展，尤其是基因芯片技術(shù)和第二代測序技術(shù)的發(fā)展，基因表達(dá)數(shù)據(jù)迅速增長，利用系統(tǒng)生物學(xué)網(wǎng)絡(luò)的方法，對高通量測序技術(shù)下的基因表達(dá)數(shù)據(jù)進(jìn)行分析和挖掘，成為了生物信息學(xué)中重要的研究方向[1]。毛囊發(fā)育及羊毛品質(zhì)受多個(gè)分子共同調(diào)控，單一的個(gè)體分子研究不足以發(fā)現(xiàn)其機(jī)制，因此，利用網(wǎng)絡(luò)的方法構(gòu)建與毛囊發(fā)育相關(guān)的基因網(wǎng)絡(luò)勢在必行?！厩叭搜芯窟M(jìn)展】蘇博美利奴羊是在2014年培育出來的一個(gè)超細(xì)型細(xì)毛羊的新品種，該品種以澳洲美利奴羊公羊?yàn)楦副?，以中國美利奴羊、新吉?xì)毛羊、敖漢細(xì)毛羊?yàn)槟副具M(jìn)行雜交所培育的，具有著羊毛品質(zhì)好、羊毛產(chǎn)量高、羊毛平均纖維直徑為17.0—19.0 μm、羊毛平均長度為8.5—11.0 cm、抗逆性強(qiáng)等特點(diǎn)[2-3]。張立春等[4]利用PCR-SSCP方法研究了蘇博美利奴羊和東北細(xì)毛羊FGF5基因多態(tài)性對毛用性狀的影響；付雪峰等[5]研究發(fā)現(xiàn)蘇博美利奴羊MYL6基因在超細(xì)型細(xì)毛羊皮膚組織表達(dá)量顯著高于細(xì)型細(xì)毛羊；朱樺等[6]發(fā)現(xiàn)隨著蘇博美利奴羊毛囊的逐步成熟，DSG4基因的表達(dá)量呈上調(diào)趨勢。加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析（weighted gene co-expression network analysis，WGCNA）是ZHANG等[7]于2005年提出的，LANGFELDER等[8]在2008年編寫發(fā)布了WGCNA的R包。它是一種用于研究分子作用機(jī)制以及網(wǎng)絡(luò)關(guān)系的準(zhǔn)確性高的系統(tǒng)生物學(xué)方法，目前已被廣泛應(yīng)用于構(gòu)建基因共表達(dá)網(wǎng)絡(luò)。【本研究切入點(diǎn)】羊毛是紡織工業(yè)重要的原材料，研究與羊毛生長息息相關(guān)的皮膚毛囊發(fā)育相關(guān)長鏈非編碼RNA（long noncoding RNA，lncRNA）與mRNA具有重要價(jià)值，近年來，關(guān)于毛囊的研究也逐漸成為皮膚生物學(xué)領(lǐng)域中的熱點(diǎn)。毛囊是控制羊毛生長的皮膚附屬器官，其形成是由初級毛囊發(fā)育至次級毛囊的一個(gè)復(fù)雜過程，受多個(gè)基因聯(lián)合調(diào)控，目前關(guān)于毛囊發(fā)育的研究多關(guān)注個(gè)體分子，并不足以發(fā)現(xiàn)毛囊發(fā)育的分子調(diào)控機(jī)制，同時(shí)關(guān)于細(xì)毛羊毛囊發(fā)育的網(wǎng)絡(luò)研究尚未見報(bào)道?！緮M解決的關(guān)鍵問題】本研究采用基于R語言的WGCNA軟件包，挖掘蘇博美利奴羊18個(gè)皮膚樣本的轉(zhuǎn)錄組表達(dá)數(shù)據(jù)，構(gòu)建與毛囊發(fā)育相關(guān)lncRNA與mRNA的共表達(dá)模塊，對所有模塊進(jìn)行富集分析，找到與毛囊發(fā)育相關(guān)的模塊并構(gòu)建了共表達(dá)網(wǎng)絡(luò)，并確定這些模塊中的核心基因。

1 材料與方法

1.1 數(shù)據(jù)來源

本研究在課題組前期的轉(zhuǎn)錄組測序基礎(chǔ)上進(jìn)行，根據(jù)阿布來提·蘇來曼等[9]于2015年在新疆科創(chuàng)繁育中心進(jìn)行的試驗(yàn)，所采集的蘇博美利奴羊胎齡65 d（D65）、85 d（D85）、105 d（D105）、135 d（D135）的胎兒以及出生后7 d（A7）、30 d（A30）的羔羊肩胛部皮膚組織，每個(gè)時(shí)期有3個(gè)生物學(xué)重復(fù)，共18個(gè)樣本。將采集好的皮膚樣品送至上海伯豪生物技術(shù)有限公司進(jìn)行轉(zhuǎn)錄組測序（RNA-Seq），包括lncRNA測序和mRNA測序，以獲得lncRNA和mRNA的表達(dá)譜數(shù)據(jù)。后續(xù)分析于2018年在新疆畜牧科學(xué)院畜牧研究所進(jìn)行。

1.2 篩選差異基因

為了使不同基因和不同樣本之間基因的表達(dá)水平能夠具有可比性，將每個(gè)基因的表達(dá)水平標(biāo)準(zhǔn)化為FPKM（Fragments Per Kilobase of exon model per Million mapped reads）[10]。將6個(gè)發(fā)育時(shí)期的lncRNA和mRNA的表達(dá)譜分別以后一個(gè)時(shí)期與前一個(gè)時(shí)期對比，根據(jù)各基因的FPKM值，應(yīng)用edgeR包[11]進(jìn)行樣本間差異基因分析，當(dāng)q-value≤0.05且Fold-change≥2時(shí)，認(rèn)定該基因?yàn)椴町惐磉_(dá)基因。

1.3 共表達(dá)網(wǎng)絡(luò)的構(gòu)建

從表達(dá)矩陣中提取出差異表達(dá)基因的表達(dá)值，并使用WGCNA軟件包進(jìn)行共表達(dá)網(wǎng)絡(luò)的分析。使用WGCNA建立檢測到的差異表達(dá)基因和差異表達(dá)的lncRNA的共表達(dá)網(wǎng)絡(luò)，步驟：

計(jì)算兩兩基因間的Pearson相關(guān)系數(shù)，構(gòu)建表達(dá)相關(guān)矩陣，使用pick-SoftThreshold函數(shù)計(jì)算無尺度拓?fù)鋽M合指數(shù)選取合適的β值作為構(gòu)建網(wǎng)絡(luò)的軟閾值，將相關(guān)矩陣轉(zhuǎn)換為鄰接矩陣，根據(jù)鄰接矩陣建立無尺度的網(wǎng)絡(luò)關(guān)系。

在所得鄰接矩陣的基礎(chǔ)上，使用拓?fù)渲丿B方法評估兩個(gè)基因間表達(dá)模式的相關(guān)性，并構(gòu)建拓?fù)渲丿B矩陣（topological overlap matrix，TOM），根據(jù)TOM矩陣計(jì)算出基因間相異度矩陣dissTOM=1-TOM。

在相異度矩陣上進(jìn)行層次聚類以對具有高度相似的共表達(dá)關(guān)系的基因進(jìn)行分組，設(shè)置模塊至少包含30個(gè)基因，使用混合動(dòng)態(tài)樹剪切法進(jìn)行分配，得到基因共表達(dá)模塊，隨后計(jì)算每個(gè)模塊的模塊特征值（module eigengene，ME），ME表示模塊的整體表達(dá)水平，計(jì)算模塊ME間的Pearson相關(guān)系數(shù)和模塊ME間的平均距離，基于該平均距離，使用average-linkage層次聚類對各個(gè)模塊進(jìn)行聚類，合并相似度高的模塊，確定最終的共表達(dá)模塊。

1.4 模塊的富集分析

使用注釋、可視化和綜合發(fā)現(xiàn)數(shù)據(jù)庫（database for annotation，visualization and integrated discovery，DAVID）[12]分別對各個(gè)模塊中的基因進(jìn)行基因本體論（gene ontology，GO）和京都基因和基因組百科全書（kyoto encyclopedia of genes and genomes，KEGG）富集分析，設(shè)置＜0.05作為顯著性閾值，篩選出與毛囊發(fā)育相關(guān)的共表達(dá)模塊，并利用R語言的ggplot2軟件包[13]對富集結(jié)果進(jìn)行可視化。

1.5 網(wǎng)絡(luò)可視化

由WGCNA生成毛囊發(fā)育相關(guān)共表達(dá)模塊網(wǎng)絡(luò)的節(jié)點(diǎn)和邊屬性文件，輸出為Cytoscape文件。拓?fù)渲丿B度量（topological overlap measurement，Tom）是計(jì)算兩個(gè)基因之間連通性的參數(shù)，如果兩個(gè)基因間具有高的Tom值，則這兩個(gè)基因與相同的基因集有高度的連接性[14]；基因的連通度是與其相連接的基因的邊屬性之和，與網(wǎng)絡(luò)中度的概念類似。篩選Tom＞0.15且由Cytoscape 3.6.1軟件[15]保留連通度排序前60內(nèi)的mRNA和lncRNA，進(jìn)行可視化和網(wǎng)絡(luò)模塊分析。

2 結(jié)果

2.1 差異基因的篩選

通過RNA-Seq測定了lncRNA和mRNA的表達(dá)譜，利用edgeR包篩選出了差異表達(dá)（differential expression，DE）lncRNA和mRNA，去除各個(gè)對比時(shí)期重復(fù)出現(xiàn)的差異表達(dá)lncRNA與mRNA使其保留唯一，得到140個(gè)差異表達(dá)lncRNA，4 372個(gè)差異表達(dá)已知mRNA以及4 558個(gè)差異表達(dá)novel mRNA（表1），用于構(gòu)建lncRNA與mRNA共表達(dá)網(wǎng)絡(luò)。

2.2 網(wǎng)絡(luò)模塊的構(gòu)建

2.2.1 軟閾值的選擇 為了構(gòu)建符合無尺度分布的網(wǎng)絡(luò)模塊，同時(shí)為了降低計(jì)算的復(fù)雜性并減少計(jì)算時(shí)間，選擇2＞0.85時(shí)，β=18作為構(gòu)建共表達(dá)模塊的軟閾值較為合理（圖1），同時(shí)對應(yīng)的平均連通度為70.5。

2.2.2 網(wǎng)絡(luò)模塊的鑒定 在選擇加權(quán)系數(shù)β之后，獲得基因之間的相異度矩陣dissTOM，對該相異度矩陣dissTOM進(jìn)行層次聚類得到基因聚類樹（圖2中Cluster Dendrogram），使用混合動(dòng)態(tài)樹剪切法從該基因聚類樹中識別出65個(gè)模塊（圖2中Dynamic Tree Cut），對這65個(gè)模塊以0.75的相關(guān)性合并相似度高的模塊后得到11個(gè)模塊（圖2中Merged dynamic）。最后，將9 070個(gè)lncRNA和mRNA分成11個(gè)共表達(dá)網(wǎng)絡(luò)模塊，各個(gè)網(wǎng)絡(luò)模塊所包含的基因數(shù)及l(fā)ncRNA數(shù)見表2，其中g(shù)rey模塊是無法分配到其他模塊的基因集合。

表1 差異表達(dá)lncRNA與mRNA數(shù)的統(tǒng)計(jì)

圖1 軟閾值的篩選

圖2 基因聚類樹與模塊構(gòu)建

表2 共表達(dá)網(wǎng)絡(luò)模塊中基因數(shù)與lncRNA數(shù)

2.3 模塊的富集分析

為了了解各個(gè)共表達(dá)模塊中基因的生物學(xué)功能，對各個(gè)模塊中的基因進(jìn)行富集分析，找到與毛囊發(fā)育相關(guān)的共表達(dá)模塊。

利用DAVID 6.8在線分析軟件分別對11個(gè)模塊中的基因進(jìn)行富集分析，以＜0.05作為顯著性閾值，在GO富集分析中有發(fā)現(xiàn)honeydew1模塊中的基因顯著富集于皮膚發(fā)育過程（GO：0043588）和毛囊發(fā)育過程（GO：0001942）等相關(guān)的生物學(xué)過程（圖3 honeydew1），paleturquoise模塊中的基因顯著富集于毛囊形態(tài)發(fā)生（GO：0031069）、皮膚發(fā)育（GO：0043588）以及毛囊發(fā)育（GO：0001942）等相關(guān)的生物學(xué)過程（圖3 paleturquoise），而skyblue2模塊中的基因顯著富集于Wnt信號通路負(fù)調(diào)控（GO：0090090）、細(xì)胞粘附（GO：0007155）等生物過程（圖3 skyblue2）。3個(gè)目標(biāo)模塊部分GO富集結(jié)果及相關(guān)基因見表3。

表3 目標(biāo)模塊部分GO富集

在KEGG富集分析中，發(fā)現(xiàn)honeydew1模塊中的基因富集在Wnt信號通路（oas04310）、TGF-β信號通路（oas04350）、Hedgehog信號通路（oas04340）、粘著力（oas04510）、ECM-受體相互作用（oas04512）等毛囊發(fā)育相關(guān)的通路（圖4 honeydew1），paleturquoise模塊中的基因富集在緊密連接（oas04530）、MAPK信號通路（oas04010）等與毛囊發(fā)育相關(guān)的通路（圖4 paleturquoise），skyblue2模塊中的基因富集于Wnt信號通路（oas04310）、ECM-受體相互作用（oas04512）、粘著力（oas04510）等毛囊發(fā)育相關(guān)的通路（圖4 skyblue2）。3個(gè)目標(biāo)模塊部分GO富集結(jié)果及相關(guān)基因見表4。

表4 目標(biāo)模塊部分KEGG富集結(jié)果

2.4 網(wǎng)絡(luò)可視化

富集分析顯示在honeydew1、paleturquoise、skyblue2模塊中的基因與毛囊發(fā)育相關(guān)，根據(jù)功能關(guān)聯(lián)原則可知這3個(gè)模塊中的基因與毛囊發(fā)育相關(guān)，通過WGCNA導(dǎo)出3個(gè)毛囊發(fā)育相關(guān)共表達(dá)模塊網(wǎng)絡(luò)的node和邊edge文件，通過篩選Tom＞0.15且由Cytoscape軟件保留連通度排序前60內(nèi)的mRNA和lncRNA，進(jìn)行網(wǎng)絡(luò)可視化。各模塊中連通度排序前5的mRNA與lncRNA見表5。

表5 連通度前5的mRNA與lncRNA

在honeydew1模塊中新基因XLOC_170356連通度最高，為該模塊中的核心基因，與876個(gè)mRNA和lncRNA相連，已知基因中ANGPT1基因連通度最高，而連通度最高的lncRNA是ENSOART00000029117（表5）。通過連通度的高低，篩選出該模塊中連通度高的38個(gè)mRNA和4個(gè)lncRNA構(gòu)建共表達(dá)網(wǎng)絡(luò)（圖5），在該網(wǎng)絡(luò)中發(fā)現(xiàn)一個(gè)有趣的基因轉(zhuǎn)化生長因子-β2（Transforming Growth Factor Beta 2，），該基因是honeydew1模塊子網(wǎng)絡(luò)中唯一一個(gè)與4個(gè)lncRNA均存在著共表達(dá)關(guān)系的基因，同時(shí)該基因顯著富集在毛囊發(fā)育這一生物過程中（表3）。

點(diǎn)的大小代表連通度高低The size of the dot represents the level of connectivity

在paleturquoise模塊中TEC基因共與1 462個(gè)mRNA和lncRNA相連接，連通度最高，為該模塊中的核心基因，而連通度最高的lncRNA是TCONS_00489976（表5）。篩選出paleturquoise模塊中連通度高的19個(gè)lncRNA和41個(gè)mRNA，進(jìn)行該模塊子網(wǎng)絡(luò)的可視化（圖6），該子網(wǎng)絡(luò)中的19個(gè)lncRNA與41個(gè)mRNA均存在著共表達(dá)關(guān)系。

在skyblue2模塊中新基因XLOC_001041共與 1 462個(gè)mRNA和lncRNA相連接，連通度最高，為該模塊中的核心基因，連通度最高的已知基因?yàn)镃DH11基因，而連通度最高的lncRNA是TCONS_00041508（表5）。篩選skyblue2模塊中連通度高的18個(gè)lncRNA和21個(gè)mRNA，進(jìn)行該模塊子網(wǎng)絡(luò)的可視化（圖7），發(fā)現(xiàn)該子網(wǎng)絡(luò)中的有3個(gè)lncRNA（TCONS_ 00376759、TCONS_00571105、TCONS_00479936）與21個(gè)mRNA均存在著共表達(dá)關(guān)系。

LncRNA Trans作用靶基因預(yù)測認(rèn)為lncRNA的功能與其共表達(dá)的基因相關(guān)，因此在本研究中l(wèi)ncRNA與mRNA的WGCNA分析同時(shí)預(yù)測出lncRNA的靶基因。在honeydew1模塊子網(wǎng)絡(luò)中發(fā)現(xiàn)ENSOART00000029117和TCONS_00236664的38個(gè)靶基因，在paleturquoise模塊子網(wǎng)絡(luò)中鑒定出TCONS_00489976等19個(gè)lncRNA的41個(gè)靶基因，在skyblue2模塊子網(wǎng)絡(luò)中發(fā)現(xiàn)3個(gè)lncRNA（TCONS_00376759、TCONS_00571105、TCONS_00479936）的21個(gè)靶基因。

3 討論

基因共表達(dá)網(wǎng)絡(luò)分析是基于系統(tǒng)生物學(xué)方法，將表達(dá)模式相似的基因進(jìn)行聚類，以網(wǎng)絡(luò)模塊的形式表現(xiàn)基因的共調(diào)控特性，使用模塊作為基因分類和功能識別的框架，使用功能關(guān)聯(lián)原則以及模塊內(nèi)研究較為透徹的基因，借助其與其他基因之間的網(wǎng)絡(luò)關(guān)系，推測其他基因的功能，有利于快速鑒定新的目標(biāo)基因[16]。

點(diǎn)的大小代表連通度高低The size of the dot represents the level of connectivity

點(diǎn)的大小代表連通度高低 The size of the dot represents the level of connectivity

WGCNA是基因共表達(dá)網(wǎng)絡(luò)的常用算法，它是一種高效、準(zhǔn)確的生物信息學(xué)、生物數(shù)據(jù)挖掘方法。與基于相關(guān)系數(shù)的共表達(dá)網(wǎng)絡(luò)算法相比，WGCNA特點(diǎn)是在閾值選擇上采用的是軟閾值的方法代替之前的硬閾值。之前的基因共表達(dá)網(wǎng)絡(luò)分析方法存在一個(gè)很明顯的局限，如將閾值定為0.85時(shí)，甚至是相關(guān)系數(shù)為0.84的基因?qū)σ矊⒈环诸愔痢盁o關(guān)”的組中，而WGCNA使用基于軟閾值的決策方法來避免此問題。

WGCNA興起于人類醫(yī)學(xué)的相關(guān)研究中，而近年來在畜禽研究中也已有所應(yīng)用，例如：MIAO等[17]為了鑒定與綿羊繁殖力有關(guān)的新lncRNA并更好地了解它們的分子機(jī)制，在鑒定差異表達(dá)mRNA和lncRNA后，利用WGCNA對重要模塊中的基因和lncRNA進(jìn)行功能富集分析，確定了與綿羊繁殖力相關(guān)的幾個(gè)基因模塊，以及由hub基因和lncRNA組成的網(wǎng)絡(luò)，這些網(wǎng)絡(luò)可能通過調(diào)節(jié)與TGF-β和催產(chǎn)素信號傳導(dǎo)途徑有關(guān)的目標(biāo)mRNA來促成綿羊多產(chǎn)性；XU等[18]為了研究在兩個(gè)綿羊品種的不同胎兒階段背最長肌有關(guān)的基因表達(dá)變化的動(dòng)態(tài)，使用時(shí)間序列表達(dá)分析在各個(gè)胎兒階段鑒定了1 472個(gè)差異表達(dá)的基因，并利用WGCNA進(jìn)行分析，識別與胎兒體重顯著相關(guān)的基因模塊，并使用具有高度重要性的基因構(gòu)建網(wǎng)絡(luò)，鑒定了，，等與胎兒體重相關(guān)的基因；PACHOLEWSKA等[19]使用RNA-seq分析了LPS誘導(dǎo)的馬外周血單核細(xì)胞中基因表達(dá)的差異，并利用加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析，識別了10個(gè)共同表達(dá)的基因模塊受LPS體外刺激調(diào)節(jié)。但是應(yīng)用WGCNA構(gòu)建細(xì)毛羊毛囊相關(guān)的研究尚未見報(bào)道。

本研究使用蘇博美利奴羊不同發(fā)育時(shí)期的轉(zhuǎn)錄組數(shù)據(jù)，利用WGCNA方法得到了11個(gè)模塊，使用DAVID進(jìn)行GO、KEGG富集分析確定了3個(gè)與毛囊發(fā)育相關(guān)的共表達(dá)模塊，其中honeydew1模塊包含1 853個(gè)mRNA和14個(gè)lncRNA，paleturquoise模塊中包含1 748個(gè)mRNA和38個(gè)lncRNA，skyblue2模塊中含有1 870個(gè)mRNA和23個(gè)lncRNA，分別從這3個(gè)模塊中篩選連通度高的基因構(gòu)建了子網(wǎng)絡(luò)。

子網(wǎng)絡(luò)分析發(fā)現(xiàn)在3個(gè)模塊中均存在著與毛囊發(fā)育相關(guān)的基因。在honeydew1模塊中，顯著富集在毛囊發(fā)育生物學(xué)過程上，并且已有研究發(fā)現(xiàn)在小鼠中能夠誘導(dǎo)靜止毛囊干細(xì)胞增殖并進(jìn)入組織再生模式[20]。在paleturquoise模塊中，發(fā)現(xiàn)41個(gè)mRNA與19個(gè)lncRNA均存在著共表達(dá)關(guān)系，其中組織蛋白酶B（cathepsin B，）據(jù)報(bào)道是小鼠毛囊中潛在的半胱氨酸蛋白酶抑制劑M/E的靶標(biāo)[21]；14-3-3σ（stratifin，）基因是一種細(xì)胞周期調(diào)節(jié)因子，在小鼠中發(fā)現(xiàn)對于毛囊發(fā)育是至關(guān)重要的，特別是對于毛干的形成[22]，而卜然等的研究表明基因可能參與了細(xì)毛羊羊毛生長的調(diào)控[23]；絲氨酸肽酶抑制劑，Kunitz 1型（serine peptidase inhibitor, Kunitz type 1，）與肝細(xì)胞生長因子激活劑抑制劑1型（hepatocyte growth factor activator inhibitor type 1，）形成一種膜結(jié)合的絲氨酸蛋白酶抑制劑，在小鼠中發(fā)現(xiàn)/在表皮角質(zhì)化和毛發(fā)發(fā)育中具有重要作用[24]；序列相似性家族83成員G（family with sequence similarity 83 member G，）在狗和小鼠的研究中發(fā)現(xiàn)參與了毛發(fā)形態(tài)發(fā)生，毛囊分化和毛囊周期[25]；任亮在藏系綿羊中的研究結(jié)果表明Gasdermin A（）基因可能通過Wnt信號通路影響毛囊的發(fā)育以及周期變化[26]，楊雪的研究表明牦牛GSDMA基因在毛生長周期中參與了毛囊細(xì)胞凋亡與增值的調(diào)控[27]，而LEI等的研究發(fā)現(xiàn)家族成員GSDMA3基因突變可以阻止小鼠毛囊周期中退行期的發(fā)生[28]；LEIVA等的研究表明髓磷脂蛋白零樣3（myelin protein zero like 3，）在控制皮膚發(fā)育，毛發(fā)生長和脂肪細(xì)胞功能中的復(fù)雜作用[29]。在skyblue2模塊中，發(fā)現(xiàn)18個(gè)lncRNA中有3個(gè)lncRNA與21個(gè)mRNA均存在著共表達(dá)關(guān)系，其中波形蛋白（vimentin，）在內(nèi)蒙古絨山羊的研究中發(fā)現(xiàn)可能通過影響外根鞘參與毛囊生長周期的調(diào)節(jié)[30]；Collins等的研究發(fā)現(xiàn)細(xì)胞維甲酸結(jié)合蛋白1（cellular retinoic acid binding protein 1，）在發(fā)育形態(tài)發(fā)生，成人毛囊周期和皮膚腫瘤形成期間被動(dòng)態(tài)調(diào)節(jié)[31]，表明CRABP1基因與毛囊相關(guān)。

子網(wǎng)絡(luò)分析中找到了部分與毛囊相關(guān)的基因，根據(jù)功能關(guān)聯(lián)原則，與這些基因存在著共表達(dá)關(guān)系的基因以及l(fā)ncRNA也與毛囊發(fā)育相關(guān)。子網(wǎng)絡(luò)的可視化結(jié)果根據(jù)連通度的高低僅顯示了排名靠前的部分基因，所篩選的3個(gè)模塊中還包含著許多已報(bào)道的與毛囊發(fā)育相關(guān)的基因，根據(jù)這些已知功能的基因，可以構(gòu)建相應(yīng)的子網(wǎng)絡(luò)，快速找到與之共表達(dá)的基因，并推測未知基因的功能。

在lncRNA Trans作用靶基因預(yù)測的基本原理是lncRNA的功能與其共表達(dá)的蛋白編碼基因相關(guān)[32]，當(dāng)lncRNA與一些距離較遠(yuǎn)的基因在表達(dá)量上存在正相關(guān)或負(fù)相關(guān)的情況時(shí)，能夠通過WGCNA共表達(dá)分析來預(yù)測lncRNA的靶基因。因此本研究的網(wǎng)絡(luò)分析中預(yù)測出ENSOART00000029117、TCONS_00489976、TCONS_00376759等24個(gè)lncRNA的靶基因，為后續(xù)進(jìn)行功能驗(yàn)證提供了參考依據(jù)。

本研究構(gòu)建了蘇博美利奴羊毛囊發(fā)育相關(guān)lncRNA與mRNA的共表達(dá)網(wǎng)絡(luò)，找到了lncRNA與mRNA間的共表達(dá)關(guān)系，但其具體調(diào)控方式還需要進(jìn)一步深入研究。

4 結(jié)論

在本研究中，利用蘇博美利奴羊皮膚組織的RNA-seq數(shù)據(jù)，通過使用加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析方法構(gòu)建共表達(dá)網(wǎng)絡(luò)模塊，將篩選出的差異基因劃分為11個(gè)模塊，利用注釋、可視化和綜合發(fā)現(xiàn)數(shù)據(jù)庫進(jìn)行富集分析發(fā)現(xiàn)honeydew1、paleturquoise、skyblue2這3個(gè)與毛囊發(fā)育相關(guān)的模塊，鑒定出許多與毛囊發(fā)育相關(guān)的潛在候選基因及l(fā)ncRNA，并預(yù)測出24個(gè)lncRNA可能的靶基因。

[1] 汪濤, 蔣慶華, 彭佳杰, 王亞東. 基因共表達(dá)網(wǎng)絡(luò)的構(gòu)建及分析方法研究綜述. 智能計(jì)算機(jī)與應(yīng)用, 2014, 4(6): 47-50.

WANG T, JIANG Q H, PENG J J, WANG Y D. A review of the construction and analysis of gene co-expression network., 2014, 4(6): 47-50.(in Chinese)

[2] ZHANG L C, SUN F L, JIN H G, DALRYMPLE B P, CAO Y, WEI T, VUOCOLO T, ZHANG M X, PIAO Q L, INGHAM A B. A comparison of transcriptomic patterns measured in the skin of Chinese fine and coarse wool sheep breeds.2017, 7(1): 14301-14312.

[3] 柳楠, 田可川, 石國慶, 賀建寧, 劉積鳳, 狄江, 楊永林. 蘇博美利奴羊核心群級進(jìn)雜交不同世代對產(chǎn)毛性狀影響的研究. 中國畜牧雜志, 2015, 51(15): 6-10.

LIU N, TIAN K C, SHI G Q, HE J N, LIU J F, DI J, YANG Y L. Effect of different generations on wool traits during upgrading cross stages in subo merino nuleus herd.2015, 51(15): 6-10.(in Chinese)

[4] 張立春, 尹峰, 柳儉強(qiáng), 李夢姝, 王春昕, 曹陽, 樸慶林, 金海國, 張明新. 綿羊FGF5基因Exon 3多態(tài)性對毛用性狀的影響. 中國草食動(dòng)物科學(xué), 2018, 38(06): 1-5.

ZHANG L C, YIN F, LIU J Q, LI M S, WANG C X, CAO Y, PIAO Q L, JIN H G, ZHANG M X. Relationship between polymorphisms of FGF5 gene Exon 3 and wool traits., 2018, 38(06): 1-5.(in Chinese)

[5] 付雪峰, 黃錫霞, 劉春潔, 徐新明, 張艷花, 吳偉偉, 于麗娟, 石剛, 田可川. 不同羊毛纖維直徑細(xì)毛羊MYL6基因克隆及差異表達(dá)研究. 中國畜牧雜志, 2016, 52(13): 5-11.

FU X F, HUANG X X, LIU C J, XU X M, ZHANG Y H, WU W W, YU L J, SHI G, TIAN K C. Study on clone and differential expression of MYL6 gene in fine-wool sheep with different fiber diameter.2016, 52(13): 5-11.(in Chinese)

[6] 朱樺, 潘家華, 徐新明, 付雪峰, 田月珍, 阿布來提·蘇來曼, 黃錫霞, 田可川. Dsg4基因在蘇博美利奴羊毛囊發(fā)育不同時(shí)期及不同組織中的表達(dá)差異. 中國畜牧雜志, 2018, 54(10): 47-51.

ZHU H, PAN J H, XU X M, FU X F, TIAN Y Z, ABLAT S, HUANG X X, CHUANTIAN K. Differential expression of Dsg4 gene in different stages of hair follicles development and different tissues in subo merino sheep., 2018, 54(10): 47-51.(in Chinese)

[7] ZHANG B, HORVATH S. A general framework for weighted gene co-expression network analysis., 2005, 4(1): 1-37.

[8] LANGFELDER P, HORVATH S. WGCNA: an R package for weighted correlation network analysis., 2008, 9(1): 559-571.

[9] 阿布來提·蘇來曼, 田月珍, 許琴, 何軍敏, 趙冰茹, 徐新明, 付雪峰, 熱娜古麗·艾尼瓦爾江, 黃錫霞, 田可川. 蘇博美利奴羊胎兒皮膚毛囊結(jié)構(gòu)及形態(tài)發(fā)育. 中國農(nóng)業(yè)科學(xué), 2017, 50(16): 3226-3235.

ABLAT S, TIAN Y Z, XU Q, HE J M, ZHAO B R, XU X M, FU X F, RANAGUL A, HUANG X X, TIAN K C. Study on fetal skin hair follicle structure and morphological development of Subo Merino sheep., 2017, 50(16): 3226-3235.(in Chinese)

[10] MORTAZAVI A, WILLIAMS B A, MCCUE K, SCHAEFFER L, WOLD B. Mapping and quantifying mammalian transcriptomes by RNA-Seq.2008, 5(7): 621-628.

[11] ROBINSON M D, MCCARTHY D J, SMYTH G K. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data., 2010, 26(1): 139-140.

[12] HUANG D W, SHERMAN B T, LEMPICKI R A. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources.2009, 4(1): 44-57.

[13] 哈德利?威克姆. ggplot2: 數(shù)據(jù)分析與圖形藝術(shù). 西安: 西安交通大學(xué)出版社, 2013.

HADLEY W. ggplot2.. Xi'an: Xi'an Jiaotong University Press, 2013.(in Chinese)

[14] CHEN Y C, GUO Y F, HE H, LIN X, WANG X F, ZHOU R, LI W T, PAN D Y, SHEN J, DENG H W. Integrative analysis of genomics and transcriptome data to identify potential functional genes of BMDs in females., 2016, 31(5): 1041-1049.

[15] SHANNON P, MARKIEL A, OZIER O, BALIGA N S, WANG J T, RAMAGE D, AMIN N, SCHWIKOWSKI B, IDEKER T. Cytoscape: a software environment for integrated models of biomolecular interaction networks., 2003, 13(11): 2498-2504.

[16] 魏凱, 張婷婷, 馬磊. 豬基因共表達(dá)網(wǎng)絡(luò)模塊的構(gòu)建及功能分析. 畜牧獸醫(yī)學(xué)報(bào), 2017, 48(11): 2205-2215.

WEI K, ZHANG T T, MA L. Construction and functional analysis of gene co-expression network modules.2017, 48(11): 2205-2215.(in Chinese)

[17] MIAO X Y, LUO Q M, ZHAO H J, QIN X Y. Co-expression analysis and identification of fecundity-related long non-coding RNAs in sheep ovaries.2016, 6(1): 39398-39407.

[18] XU L Y, ZHAO F P, REN H X, LI L, LU J, LIU J S, ZHANG S F, LIU G E, SONG J Z, ZHANG L, WEI C H, DU L X. Co-expression analysis of fetal weight-related genes in ovine skeletal muscle during mid and late fetal development stages., 2014, 10(9): 1039-1050.

[19] PACHOLEWSKA A, MARTI E, LEEB T, JAGANNATHAN V, GERBER V. LPS-induced modules of co-expressed genes in equine peripheral blood mononuclear cells., 2017, 18(1): 34-45.

[20] OSHIMORI N, FUCHS E. Paracrine TGF-beta signaling counterbalances BMP-mediated repression in hair follicle stem cell activation., 2012, 10(1): 63-75.

[21] OORTVELD M A W, VAN VLIJMEN WILLEMS I M J J, KERSTEN F F J, CHENG T, VERDOES M, VAN ERP P E J, VERBEEK S, REINHECKEL T, HENDRIKS W J A J, SCHALKWIJK J, ZEEUWEN P L J M. Cathepsin B as a potential cystatin M/E target in the mouse hair follicle.2017, 31(10): 4286-4294.

[22] HAMMOND N L, HEADON D J, DIXON M J. The cell cycle regulator protein 14-3-3sigma is essential for hair follicle integrity and epidermal homeostasis.2012, 132(6): 1543-1553.

[23] 卜然, 李晶, 劉開東, 劉積鳳, 柳楠. 細(xì)毛羊肩頸部與腹部皮膚組織SFN基因差異表達(dá)的研究. 黑龍江畜牧獸醫(yī), 2014(19): 80-82.

BU R, LI J, LIU K D, LIU J F, LIU N. Study on the differential expression of the SFN gene in the skin tissues of the shoulder /neck and abdomen in fine-wool sheep.2014(19): 80-82.(in Chinese)

[24] NAGAIKE K, KAWAGUCHI M, TAKEDA N, FUKUSHIMA T, SAWAGUCHI A, KOHAMA K, SETOYAMA M, KATAOKA H. Defect of hepatocyte growth factor activator inhibitor type 1/serine protease inhibitor, Kunitz type 1 (Hai-1/Spint1) leads to ichthyosis- like condition and abnormal hair development in mice., 2008, 173(5): 1464-1475.

[25] BALMER P, FELLAY A K, SAYAR B S, HARITON W V J, WIENER D J, GALICHET A, MULLER E J, ROOSJE P J. FAM83G/ Fam83g genetic variants affect canine and murine hair formation., 2018, 28(4): 13729-13738.

[26] 任亮. 藏系綿羊GSDMA基因皮膚特異表達(dá)載體的構(gòu)建 [D]. 成都: 西南民族大學(xué), 2015.

REN L. Construction of skin-specific expression vector for GSDMA gene of Tibetan sheep [D]. Chengdu: Southwest Minzu University, 2015.(in Chinese)

[27] 楊雪. 牦牛皮膚毛囊周期性結(jié)構(gòu)變化規(guī)律及其相關(guān)調(diào)控因子的研究 [D]. 蘭州: 甘肅農(nóng)業(yè)大學(xué), 2017.

YANG X. Study on histological characteristics of hair follicle and expression pattern of related factors in the skin during hair cycle in yak [D]. Lanzhou: Gansu Agricultural University, 2017.(in Chinese)

[28] LEI M X, XIANG G, LI Y, TIAN Y, LIAN X H. Gsdma3 gene is needed for the induction of apoptosis-driven catagen during mouse hair follicle cycle.2011, 136(3): 335-343.

[29] LEIVA A G, CHEN A L, DEVARAJAN P, CHEN Z, DAMANPOUR S, HALL J A, BIANCO A C, LI J, BADIAVAS E V, ZAIAS J, MITEVA M, ROMANELLI P, NOURI K, CAO WIKRAMANAYAKE T. Loss of Mpzl3 function causes various skin abnormalities and greatly reduced adipose depots.2014, 134(7): 1817-1827.

[30] RILE N, LIU Z H, GAO L X, QI J K, ZHAO M, XIE Y C, SU R, ZHANG Y J, WANG R J, LI J, XIAO H M, LI J Q. Expression of vimentin in hair follicle growth cycle of inner Mongolian cashmere goats.2018, 19(1): 38-45.

[31] COLLINS C A, WATT F M. Dynamic regulation of retinoic acid-binding proteins in developing, adult and neoplastic skin reveals roles for beta-catenin and Notch signalling.2008, 324(1): 55-67.

[32] YAN P X, LUO S, LU J Y Y, SHEN X H. Cis- and trans-acting lncRNAs in pluripotency and reprogramming.2017, 46: 170-178.

Construction of Co-expression Network of lncRNA and mRNA Related to Hair Follicle Development of Subo Merino Sheep

CHEN HuaFeng1,2, TIAN KeChuan2, HUANG XiXia1, Ablat Sulayman1, HE JunMin3, TIAN YueZhen2, XU XinMing2, FU XueFeng2, ZHAO BingRu4, ZHU Hua1, Hanikezi Tulafu2

(1College of Animal Science, Xinjiang Agricultural University, Urumqi 830052;2Institute of Animal Science, Xinjiang Academy of Animal Science, Urumqi 830011;3College of Animal Science and Technology, Gansu Agricultural University, Lanzhou 730070;4College of Animal Science and Technology, China Agricultural University, Beijing 100083)

【】With the development and improvement of high-throughput sequencing technology, massive transcriptome sequencing data has emerged, and more and more methods of gene network have been used in research; the hair follicles development of fine wool sheep is regulated by multiple genes and lncRNAs. It is not enough to study a certain molecule alone to discover its regulatory mechanism. The aim of this study was to construct a co-expression network of lncRNAs and mRNAs related to hair follicle development, and to explore potential candidate genes for hair follicle development. 【】The study was performed on fetal skin tissues of the 65th, 85th, 105th, and 135th days and lambs skin tissues of the 7th and 30th days from Subo Merino sheep, with 3 biological replicates in each period, a total of 18 samples were sequenced by transcriptome to obtain lncRNAs and mRNAs expression profiles of 6 different developmental stages, and the differential expressed lncRNAs and mRNAs in adjacent periods was screened to construct a co-expression module by Weighted Gene Co-expression Network Analysis (WGCNA) method, and GO (Gene Ontology) and KEGG pathway enrichment analysis of DAVID online tool was used to find hair follicle related modules. Finally, the high-interconnect lncRNAs and mRNAs were screened from the target module for network visualization using Cytoscape software.【】From the expression profile data, 9070 differentially expressed lncRNAs and mRNAs were screened, and 11 modules were obtained by WGCNA method. The DAVID enrichment analysis revealed that the genes in the honeydew1 module, paleturquoise module, and skyblue2 module were involved in biological processes such as skin development, hair follicle development, hair follicle morphogenesis, negative regulation of Wnt signaling pathway, and cell adhesion, et al; and the signaling pathways involved in hair follicle development such as Wnt signaling pathway, TGF-β signaling pathway, Hedgehog signaling pathway, tight junction, MAPK signaling pathway, ECM-receptor interaction, et al. The sub-networks of lncRNAs and mRNAs with high connectivity were screened in these 3 modules, and the hair follicle development related genes including,,,,,,,andwere obtained, and possible target genes of 24 lncRNAs including ENSOART00000029117, TCONS_00489976, TCONS_00376759 were predicted.【】In this study, we first constructed a co-expression network of lncRNAs and mRNAs related to the hair follicles development of Subo merino sheep using WGCNA method, identified 3 co-expression modules related to hair follicle development, and found several potential candidate genes related to hair follicle development,and predicted possible target genes of 24 lncRNAs.

Subo Merino sheep; WGCNA; mRNA; lncRNA; hair follicle development

10.3864/j.issn.0578-1752.2019.19.016

2018-12-19；

2019-07-03

國家自然科學(xué)基金（31460593，31360543）、國家現(xiàn)代農(nóng)業(yè)產(chǎn)業(yè)技術(shù)體系建設(shè)項(xiàng)目（CARS-40）、新疆絨毛用羊遺傳育種與繁殖自治區(qū)重點(diǎn)實(shí)驗(yàn)室項(xiàng)目（XJYS1105）、新疆農(nóng)業(yè)大學(xué)研究生科研創(chuàng)新項(xiàng)目（XJAUGRI-2017013）

陳華峰，E-mail：1504159731@qq.com。

黃錫霞，E-mail：au-huangxixia@163.com。通信作者田可川，E-mail：tiankechuan@163.com

（責(zé)任編輯 林鑒非）