任一帆，孟 婕，楊 晴，張 貝，李可強(qiáng)，張傳生*，耿立英，趙雪艷，李祥龍

(1 河北科技師范學(xué)院 河北省特色動(dòng)物種質(zhì)資源挖掘與創(chuàng)新重點(diǎn)實(shí)驗(yàn)室，河北 秦皇島，066600；2 山東省農(nóng)業(yè)科學(xué)院畜牧獸醫(yī)研究所)

卵巢作為雞繁殖活動(dòng)的基礎(chǔ)，近年來受到研究人員的廣泛關(guān)注，特別是對(duì)雞胚卵巢轉(zhuǎn)錄組的研究。例如：曹頂國(guó)[1，2]對(duì)高、低產(chǎn)汶上蘆花雞的卵巢組織轉(zhuǎn)錄組分析，分別獲得有功能注釋的基因4 675個(gè)和4 472個(gè)，并以低產(chǎn)汶上蘆花雞為對(duì)照，在高產(chǎn)個(gè)體的卵巢組織中獲得了158個(gè)差異表達(dá)基因。前人在雞胚卵巢后期生長(zhǎng)發(fā)育的轉(zhuǎn)錄組研究取得了階段性研究成果，但對(duì)雞胚胎期的卵巢發(fā)育研究還很匱乏，尤其缺乏利用系統(tǒng)生物學(xué)網(wǎng)絡(luò)分析的方法對(duì)雞胚胎期卵巢組織轉(zhuǎn)錄組基因的研究。此外，雞卵巢發(fā)育受多個(gè)分子共同調(diào)控，單一的個(gè)體分子研究不足以發(fā)現(xiàn)其機(jī)制，而加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析作為一種系統(tǒng)生物學(xué)方法，目前在多個(gè)領(lǐng)域被廣泛應(yīng)用[3～6]。通過WGCNA可系統(tǒng)的分析基因之間的相互關(guān)系，發(fā)現(xiàn)共表達(dá)基因的調(diào)控中心，并將基因模塊與性狀相聯(lián)系，從而找出目標(biāo)性狀的關(guān)鍵基因[7～10]。為此，筆者以不同發(fā)育階段的轉(zhuǎn)錄組測(cè)定的卵巢組織雞的RNA-Seq數(shù)據(jù)為基礎(chǔ)，利用WGCNA分析雞卵巢發(fā)育過程中的基因表達(dá)譜，挖掘出與雞卵巢發(fā)育高度相關(guān)的基因，并試圖篩選出關(guān)鍵的調(diào)控基因，為雞卵巢發(fā)育具體機(jī)制的研究提供幫助。

1 材料與方法

1.1 數(shù)據(jù)獲取

雞卵巢在不同發(fā)育階段的轉(zhuǎn)錄組測(cè)序數(shù)據(jù)來自文獻(xiàn)[11]，樣品的RNA-Seq數(shù)據(jù)在NCBI的登記號(hào)碼為：SRP081829。該數(shù)據(jù)以白來航雞為材料，分別選取孵化第6天(E6)、E12和出殼第1天(D1)的左右側(cè)性腺，進(jìn)行高通量測(cè)序，每組包括2個(gè)生物學(xué)重復(fù)。所有基因的表達(dá)量通過FPKM標(biāo)準(zhǔn)化，并去除標(biāo)準(zhǔn)差小于5的基因。

1.2 加權(quán)基因共表達(dá)網(wǎng)絡(luò)構(gòu)建

通過R軟件(R version 3.6.2)的WGCNA包(版本1.6.9)提供的函數(shù)構(gòu)建加權(quán)基因共表達(dá)網(wǎng)絡(luò)[12，13]。為了盡可能的滿足無尺度網(wǎng)絡(luò)分布要求，調(diào)用WGCNA中的函數(shù)pick Soft Threshold，設(shè)置合適的軟閾值β，并采用WGCNA自帶的動(dòng)態(tài)切割法劃分基因模塊，設(shè)定模塊的最小基因總數(shù)為100，切割高度為0.85，其他參數(shù)參考函數(shù)的默認(rèn)值。

1.3 組織特異性模塊鑒定

計(jì)算出模塊基因的模塊特征向量(module eigengene, ME)，模塊的特征向量與不同組織性狀之間的相關(guān)性系數(shù)r>0時(shí)表示模塊表達(dá)量與性狀呈正相關(guān)[14]。繪出模塊與性狀的相關(guān)性熱圖，本次研究定義相關(guān)系數(shù)r>0.65且顯著性P<0.01的模塊為組織特異性模塊。

1.4 模塊的富集分析

通過R軟件提供的clusterProfiler，enrichplot等軟件包，對(duì)所有8個(gè)模塊進(jìn)行GO(gene ontology，GO)和KEGG (Kyoto Encyclopedia of Genes and Genomes)分析。GO分析主要是將每個(gè)模塊的基因分為生物過程(Biological Process，BP)、分子功能(Molecular Function，MF)以及細(xì)胞組分(Cellular Component，CC)等3個(gè)部分，經(jīng)過GO分析可以明確每個(gè)模塊基因具有的生物學(xué)作用。KEGG分析通過信號(hào)通路的顯著性富集，確定差異表達(dá)基因參與的主要生化代謝途徑和信號(hào)轉(zhuǎn)導(dǎo)途徑。設(shè)置閾值p<0.1作為顯著性閾值，篩選出與雞卵巢發(fā)育相關(guān)的基因，并利用R語言的ggplot2軟件包[15]對(duì)富集結(jié)果進(jìn)行可視化。

1.5 hub gene的篩選

根據(jù)模塊中各基因的模塊連通性和臨床性狀關(guān)系選擇候選hub基因。模塊連通性定義為基因間皮爾遜相關(guān)性的絕對(duì)值(模塊成員，即MM)，臨床性狀關(guān)系定義為各基因與性狀的皮爾遜相關(guān)絕對(duì)值(基因顯著性,即GS)。設(shè)置候選hub基因的MM>0.8，GS>0.6。同時(shí)，將hub模塊中的候選基因，使用搜索工具For the Retrieve of Interactive gene(String；https://string-db.org/)數(shù)據(jù)庫(kù)構(gòu)建PPI網(wǎng)絡(luò)，對(duì)基因的功能進(jìn)行注釋。將候選基因?qū)隒ytoscape篩選連接值前30的基因，并呈現(xiàn)網(wǎng)絡(luò)進(jìn)行下一步分析。

2 結(jié)果與分析

2.1 篩選表達(dá)譜基因和樣本

將所有的基因通過FPKM標(biāo)準(zhǔn)化并去除低表達(dá)量的數(shù)值之后，得到了5 557個(gè)表達(dá)基因的FPKM值。通過對(duì)樣本進(jìn)行聚類分析，結(jié)果表明，12個(gè)樣本分為2組，一組包含8個(gè)樣本，另一組為4個(gè)樣本，本次研究數(shù)據(jù)的樣本聚集度相對(duì)較好，沒有進(jìn)行離群值的刪除(圖1)。

圖1 雞胚卵巢基因樣本聚類分析

2.2 軟閾值的確定

在確定無尺度網(wǎng)絡(luò)擬合指數(shù)R2>0.85時(shí)，軟閾值β確定為16比較合適，在β=16時(shí)所對(duì)應(yīng)的平均連通度為257，最大的連接系數(shù)為882，中位數(shù)為171(圖2)。

注：圖中的橫軸均代表軟閾值(β)。 A: 縱坐標(biāo)對(duì)應(yīng)的是無尺度網(wǎng)絡(luò)模型指數(shù); B: 每一個(gè)軟閾值對(duì)應(yīng)的網(wǎng)絡(luò)平均連接度。 圖2 雞胚卵巢加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析軟閾值(β)的確定

2.3 基因共表達(dá)網(wǎng)絡(luò)構(gòu)建

在軟閾值為16的條件下，采用動(dòng)態(tài)剪切法構(gòu)建加權(quán)基因共表達(dá)網(wǎng)絡(luò)，共得到8個(gè)模塊(圖3)，每個(gè)模塊都代表一組相關(guān)度比較高的基因，灰色模塊代表不能被納入其他模塊的基因。各個(gè)模塊都用不同的顏色表示，各模塊包含的基因數(shù)量相差較大，最少的是黑色模塊，包含109個(gè)基因；最多的是青色模塊，包含2 118個(gè)基因，其余6個(gè)模塊包含的基因數(shù)分別為：534(黃色),602(綠色),854(黑色)和143(紅色)。

注：基因聚類樹的每一個(gè)顏色對(duì)應(yīng)一個(gè)模塊分支1～7代表不同的模塊，1：黑色模塊，2：藍(lán)色模塊，3：紅色模塊，4：綠色模塊，5：黃色模塊，6：棕色模塊，7：青色模塊，灰色代表灰色模塊代表不能被納入其他模塊的基因圖3 基于拓?fù)渲丿B構(gòu)建的雞胚卵巢基因聚類樹

注：橫軸表示不同的性狀, 縱軸表示每一個(gè)模塊的特征向量。紅色的格子代表性狀與模塊具有正相關(guān)性, 綠色的格子代表性狀與模塊具有負(fù)相關(guān)性。矩形框里的數(shù)字代表模塊與性狀之間的相關(guān)系數(shù)及相應(yīng)p值。(E6L：孵化6天雞胚左側(cè)卵巢；E6R：孵化6天雞胚右側(cè)卵巢；E12L：孵化12天雞胚左側(cè)卵巢；E12R孵化12天雞胚右側(cè)卵巢；D1L：出殼后第1天雛雞左側(cè)卵巢；D1R：出殼后第1天雛雞右側(cè)卵巢)圖 4 雞胚卵巢基因共表達(dá)網(wǎng)絡(luò)模塊與性狀相關(guān)性

2.4 組織特異性模塊鑒定

在8個(gè)模塊中，除去灰色模塊共有6個(gè)模塊與組織存在高度特異性(圖4)，分別為：黑色(r=0.76，P=4×10-3)、藍(lán)色(r=0.95，P=2×10-6)、棕色(r=0.97，P=3×10-7)、青色(r=0.88，P=2×10-4)、紅色(r=0.81，P=10-3)與黃色(r=0.96，P=10-6)。大多數(shù)模塊都有與其高度相關(guān)聯(lián)的組織,例如棕色模塊與D1L卵巢組織存在高度相關(guān)性,并且連接度也最高。而綠色模塊沒有與任何一個(gè)組織存在高度相關(guān)性，通過篩選最終得到6個(gè)模塊進(jìn)行下一步的富集分析以及目標(biāo)模塊的鑒定。

2.5 雞卵巢發(fā)育相關(guān)模塊的篩選

為進(jìn)一步了解各模塊的生物學(xué)功能，本次研究對(duì)篩選的6個(gè)模塊的KEGG通路和GO富集分析進(jìn)行深度解析。

在孵化期的第6天，雞卵巢特征顯現(xiàn)，只有右側(cè)卵巢組織有與之顯著相關(guān)的模塊：黑色，該模塊富集到(GO:0002443)白細(xì)胞介導(dǎo)的免疫反應(yīng)以及跟免疫相關(guān)的94條GO 條目，并富集到糖異生，丙酮酸代謝等9條KEGG通路。在孵化進(jìn)行到第12天時(shí)，與左側(cè)卵巢組織顯著相關(guān)的青色模塊可以富集的GO條目有：(GO:0008585)雌性性腺發(fā)育，(GO:0061458)生殖系統(tǒng)發(fā)育，(GO:0030111) Wnt信號(hào)通路的調(diào)控等共328條，而KEGG富集的通路有：細(xì)胞周期，卵母細(xì)胞減數(shù)分裂，p53信號(hào)通路，TGF-β信號(hào)通路等共25條；與右側(cè)卵巢高度相關(guān)的模塊有2個(gè),分別是紅色模塊和黃色模塊，紅色模塊顯著富集到：(GO:0006342)染色質(zhì)沉默，(GO:0016458)基因沉默，(GO:0004857)酶抑制劑的活動(dòng)等122條GO條目，富集到的KEGG通路有：氨基酸的生物合成等共17條。黃色模塊共富集到 151 個(gè) GO 條目, 其中多個(gè)GO條目與體液免疫反應(yīng)以及防御反應(yīng)相關(guān)，KEGG共富集到39條通路，大多是酪氨酸代謝和精氨酸生物合成等與氨基酸代謝相關(guān)的通路。

出殼第1天，左側(cè)卵巢組織與棕色模塊存在最高的相關(guān)性，棕色模塊富集到(GO:0001838)胚胎上皮管形成等166條GO條目，并且顯著富集到孕酮介導(dǎo)的卵母細(xì)胞成熟和GnRH信號(hào)通路等9條KEGG通路。而右側(cè)卵巢組織與藍(lán)色模塊相關(guān)性最高，藍(lán)色模塊主要富集于(GO:0055085)跨膜轉(zhuǎn)運(yùn)等248條GO條目，以及氨基酸合成與代謝等21條KEGG信號(hào)通路。

雖然6個(gè)模塊各自都有相關(guān)的性狀，但根據(jù)富集分析結(jié)果發(fā)現(xiàn)，與雞卵巢發(fā)育相關(guān)的通路主要集中在青色模塊，并且該模塊富集到本試驗(yàn)領(lǐng)域已經(jīng)報(bào)道的Wnt通路。右側(cè)卵巢的發(fā)育停滯以及萎縮主要與紅色模塊的基因相關(guān)，所以選擇青色模塊與紅色模塊作為特征性模塊進(jìn)行下一步的分析。青色模塊與紅色模塊的GO富集分析氣泡見圖5。

圖5 雞胚卵巢基因共表達(dá)網(wǎng)絡(luò)模塊中青色模塊與紅色模塊的GO富集分析

2.6 樞紐基因的篩選

模塊內(nèi)基因的連通性(Connectivity)代表該基因與其他基因之間的調(diào)控關(guān)系，連通性越高，說明該基因在模塊中的調(diào)控作用越大，而樞紐基因的連通性通常情況下是在該模塊中連通性較大的基因。參考WGCNA 的標(biāo)準(zhǔn)，對(duì)特征性進(jìn)行樞紐基因篩選，確定了2個(gè)模塊的候選樞紐基因數(shù)量，分別為：1 863(青色)，442(紅色)。每個(gè)模塊MM值排名前10的候選樞紐基因見附表1。將候選樞紐基因利用STRING計(jì)算出其互作關(guān)系，輸出為Cytoscape文件，并借助STRING網(wǎng)站注釋了基因的功能(附表2)。

附表1 雞胚卵巢基因共表達(dá)網(wǎng)絡(luò)模塊中的6個(gè)模塊前10個(gè)樞紐基因

附表2 雞胚卵巢基因共表達(dá)網(wǎng)絡(luò)模塊候選樞紐基因的功能注釋

將特征性模塊的候選基因利用Cytoscape軟件的cytoHubba插件根據(jù)程度值篩選出連接度前30的基因作為中心網(wǎng)絡(luò)，候選基因的連線分別為：青色模塊205條，紅色模塊328條。

使用Cytoscape將結(jié)果可視化，青色模塊和紅色模塊的中心網(wǎng)絡(luò)見圖6。

注：在網(wǎng)絡(luò)中的每一個(gè)基因代表一個(gè)節(jié)點(diǎn)，顏色由深到淺代表連接度從大到小。A:青綠色模塊的中心網(wǎng)絡(luò)，B：紅色模塊的中心網(wǎng)絡(luò)。圖6 雞胚卵巢基因共表達(dá)網(wǎng)絡(luò)模塊中青色模塊和紅色模塊的中心網(wǎng)絡(luò)結(jié)構(gòu)

青色模塊的核心基因CDC5L與紅色模塊的核心基因PLG在卵巢的表達(dá)量變化見圖7?？梢钥闯鯟DC5L基因與PLG基因都在胚胎左側(cè)卵巢發(fā)育的第12天表達(dá)量激增，不同的是CDC5L在左側(cè)卵巢表達(dá)量增加，PLG基因在右側(cè)卵巢表達(dá)量增加。此外，PLG基因在胚胎期第12天右側(cè)卵巢的表達(dá)量遠(yuǎn)大于其他時(shí)期，但CDC5L基因在同時(shí)期卵巢左側(cè)的表達(dá)量與出殼第1天的數(shù)值相差不大。

A:青色模塊的核心基因表達(dá)量，B：紅色模塊的核心基因表達(dá)量圖7 雞胚卵巢基因共表達(dá)網(wǎng)絡(luò)模塊中青色模塊和紅色模塊的核心基因表達(dá)量

3 討 論

3.1 WGCNA 分析的重要意義

傳統(tǒng)的生物學(xué)研究側(cè)重于從分子水平闡釋單個(gè)功能元件(如DNA, mRNA和蛋白質(zhì))對(duì)生命活動(dòng)的影響,只能局部地解釋某一生命活動(dòng)發(fā)生的原因。而新興的共表達(dá)網(wǎng)絡(luò)分析方法可以將復(fù)雜的數(shù)據(jù)進(jìn)行歸納和整理，高效研究基因整體表達(dá)規(guī)律。相比于其他調(diào)控網(wǎng)絡(luò),WGCNA可以特異地篩選出與性狀相關(guān)的基因，進(jìn)行模塊化分類,得到具有高度生物學(xué)意義的共表達(dá)模塊,已經(jīng)被證明是一種高效的數(shù)據(jù)挖掘手段[16]。

3.2 雞胚卵巢發(fā)育相關(guān)共表達(dá)模塊的鑒定

本次研究根據(jù)雞卵巢不同發(fā)育時(shí)期的轉(zhuǎn)錄組數(shù)據(jù)，利用 WGCNA 方法進(jìn)行分析，最終獲得了8個(gè)共表達(dá)模塊，通過對(duì)6個(gè)組織特異性模塊進(jìn)行富集分析,發(fā)現(xiàn)這些模塊均可以得到具有生物學(xué)意義的調(diào)控通路。例如：青色模塊富集到(GO:0007049)細(xì)胞周期等細(xì)胞過程，藍(lán)色模塊可富集到(GO:0015267)通道活性，(GO:0055085)跨膜轉(zhuǎn)運(yùn)等通路。富集結(jié)果提示了參與雞卵巢發(fā)育生理活動(dòng)的差異基因除了集中于細(xì)胞內(nèi)，還廣泛分布于細(xì)胞膜上，該結(jié)論與蘭道亮等[17]對(duì)牦牛卵巢轉(zhuǎn)錄組研究獲得的結(jié)果相一致，這說明細(xì)胞膜組分在卵巢生理生殖活動(dòng)中可能也扮演了重要角色。

3.3 雞胚卵巢發(fā)育節(jié)點(diǎn)特征

雞胚的原始生殖腺位于中腎嵴腹內(nèi)側(cè)，由中胚層發(fā)育而來，在孵化期第3.5天左右出現(xiàn)[18]；性腺分化在孵化期第5.5天出現(xiàn)，到孵化期第6天時(shí)通過形態(tài)學(xué)方法能夠分辨出性別[19]，性腺分化在孵化期第7天時(shí)差異更顯著；原始生殖細(xì)胞在孵化期第11天左右完全分化[20]。卵巢皮質(zhì)部在孵化期第11～12天時(shí)開始出現(xiàn)原始卵泡，并且數(shù)量不斷增加。卵原細(xì)胞(oogonia)減數(shù)分裂生成初級(jí)卵母細(xì)胞出現(xiàn)在孵化期第15.5天[21，22]。到了孵化期第16～18天，雌性左側(cè)卵巢繼續(xù)發(fā)育，外觀似腔卵泡樣結(jié)構(gòu)；而右側(cè)性腺退化，外觀呈睪丸樣結(jié)構(gòu)。雌雞性腺一般情況下為左側(cè)正常發(fā)育，右側(cè)退化并逐漸萎縮[23，24]。但有學(xué)者研究發(fā)現(xiàn)，在雌雞性腺發(fā)育過程中會(huì)出現(xiàn)雙側(cè)功能卵巢閉鎖，并證實(shí)該現(xiàn)象并非巧合，原因之一是等位基因的組合配對(duì)導(dǎo)致出現(xiàn)這樣的表型[25]。

3.4 雞胚卵巢發(fā)育基因的篩選

本次研究根據(jù)WGCNA樞紐基因的篩選標(biāo)準(zhǔn)：GS>0.6，MM>0.8，并結(jié)合富集分析結(jié)果，篩選出2個(gè)特定模塊的候選樞紐基因；并將候選基因利用在線網(wǎng)絡(luò)工具：STRING，根據(jù)已有的研究數(shù)據(jù)建立基因之間的連通性，最后將基因間的連通性信息導(dǎo)入到 Cytoscape 軟件中，篩選每個(gè)模塊連接度前30的基因，并根據(jù)STRING對(duì)每個(gè)基因的注釋，進(jìn)一步識(shí)別出每個(gè)模塊的樞紐基因。

在孵化期第12天，卵巢顯著相關(guān)的模塊所參與的生物調(diào)控通路差別較大。左側(cè)卵巢與青色模塊高度相關(guān)，參照富集分析結(jié)果，該模塊主要調(diào)控卵母細(xì)胞分裂并調(diào)控卵巢的發(fā)育；模塊的核心基因CDC5L基因大量存在于細(xì)胞核，在mRNA剪接前的第2個(gè)催化步驟起重要作用；CDC5L在癌癥以及腫瘤方面的研究成果較多，研究表明：CDC5L對(duì)于豬卵母細(xì)胞減數(shù)分裂和早期胚胎發(fā)育至關(guān)重要，但是該基因在雞胚發(fā)育方面研究較少，本次試驗(yàn)結(jié)果可能提示CDC5L在雞胚胎卵巢發(fā)育早期起到重要作用。右側(cè)卵巢與紅色模塊顯著相關(guān)，富集分析顯示紅色模塊主要調(diào)控基因以及染色質(zhì)的沉默，并參與酶抑制劑的合成活動(dòng)，右側(cè)性腺在該時(shí)期處于停滯發(fā)育并出現(xiàn)萎縮的時(shí)期，富集分析結(jié)果與該時(shí)期節(jié)點(diǎn)特征相一致。紅色模塊的核心基因PLG在胚胎右側(cè)卵巢發(fā)育的第12天表達(dá)量激增，該基因編碼蛋白為纖維溶酶原，在包括胚胎發(fā)育，組織重塑在內(nèi)的多種過程中充當(dāng)?shù)鞍姿庖蜃?，但目前在雞胚卵巢發(fā)育方面對(duì)PLG的研究資料很少，本次試驗(yàn)結(jié)果可能提示該蛋白在胚胎右側(cè)卵巢發(fā)育停滯時(shí)期起到關(guān)鍵作用。此外，PLG大量存在于細(xì)胞外，這個(gè)結(jié)果提示核心基因不僅在細(xì)胞內(nèi)，在細(xì)胞外也同樣存在。

4 結(jié) 論

利用權(quán)重基因共表達(dá)網(wǎng)絡(luò)分析(WGCNA)的方法將5 558個(gè)差異基因劃分成8個(gè)共表達(dá)模塊，將6個(gè)有意義的模塊完成了模塊內(nèi)基因的功能富集分析，揭示了模塊所蘊(yùn)含的生物學(xué)意義，得到不同發(fā)育階段卵巢組織差異表達(dá)基因的功能和代謝通路，并篩選出2個(gè)特征模塊的樞紐基因：CDC5L和PLG。