利用WGCNA進(jìn)行棉花果枝節(jié)間伸長(zhǎng)相關(guān)基因共表達(dá)模塊鑒定

巨飛燕，張思平，劉紹東，馬慧娟，陳靜，葛常偉，沈倩，張小萌，劉瑞華，趙新華，張永江，龐朝友*

（1. 棉花生物學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室河北基地/河北農(nóng)業(yè)大學(xué)農(nóng)學(xué)院，河北保定071001；2. 棉花生物學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室/中國(guó)農(nóng)業(yè)科學(xué)院棉花研究所，河南安陽(yáng)455000）

棉花（Gossypium hirsutumL.）是全世界主要纖維作物之一，也是我國(guó)重要的經(jīng)濟(jì)作物，在國(guó)民經(jīng)濟(jì)中有著非常重要的地位[1-2]。合理密植是棉花生產(chǎn)中提高綜合產(chǎn)量的一項(xiàng)關(guān)鍵栽培技術(shù)[3]，適當(dāng)提高棉花的種植密度是提高產(chǎn)量的有效途徑之一，然而寬大松散的株型結(jié)構(gòu)已成為限制棉花種植密度的瓶頸。因此合理的植物株型和群體結(jié)構(gòu)對(duì)棉花的高產(chǎn)種植具有重要的意義。棉花的株型結(jié)構(gòu)主要由主莖高度、根的數(shù)量及長(zhǎng)度、主莖和果枝的節(jié)間長(zhǎng)度、棉鈴生長(zhǎng)結(jié)構(gòu)及分布等組成[4]，其中果枝類(lèi)型是棉花株型表現(xiàn)的決定性狀[5]。開(kāi)展棉花果枝伸長(zhǎng)調(diào)控的分子機(jī)理研究就顯得意義重大。

在高通量技術(shù)時(shí)代，由芯片技術(shù)到測(cè)序技術(shù)產(chǎn)生的大量數(shù)據(jù)驅(qū)使人們需要客觀地提取有意義的信息。而著眼于單個(gè)基因的差異表達(dá)分析技術(shù)或聚類(lèi)分析等手段在面對(duì)由大量DNA、RNA及蛋白質(zhì)所構(gòu)成的復(fù)雜網(wǎng)絡(luò)時(shí)，往往錯(cuò)過(guò)很多有效信息，因此基于全局基因表達(dá)模式的生物網(wǎng)絡(luò)分析法應(yīng)運(yùn)而生。當(dāng)前應(yīng)用較廣泛的是權(quán)重基因共表達(dá)網(wǎng)絡(luò)分析（Weighted gene co-expression network analysis,WGCNA）方法。WGCNA 是一種基于大樣本轉(zhuǎn)錄組數(shù)據(jù)的研究方法，其首先假定基因網(wǎng)絡(luò)服從無(wú)尺度分布，根據(jù)功能相似的基因往往具有類(lèi)似的表達(dá)變化構(gòu)建基因模塊，然后對(duì)模塊進(jìn)行深入挖掘，比如篩選模塊的樞紐基因（Hub gene）、模塊與性狀的關(guān)聯(lián)分析、模塊富集分析和建立基因表達(dá)網(wǎng)絡(luò)等[6-7]。自該技術(shù)問(wèn)世以來(lái)，因其強(qiáng)大的分析效能，在多個(gè)研究領(lǐng)域中得到廣泛的運(yùn)用。Farber[8]通過(guò)整合GWAS（Genome-wide association study）和WGCNA 方法，成功發(fā)掘出在常規(guī)GWAS 分析中沒(méi)有被發(fā)現(xiàn)的與骨密度相關(guān)的TNF基因，同時(shí)驗(yàn)證了TNF基因在骨代謝中的核心地位。Greenham等[9]為研究油菜（Brassica napusL.）基因組的干旱脅迫下轉(zhuǎn)錄水平表達(dá)模式的變化，利用干旱和對(duì)照處理下的48個(gè)轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行了共表達(dá)網(wǎng)絡(luò)構(gòu)建，鑒定到與干旱相關(guān)的生物學(xué)模塊，并挖掘了目標(biāo)模塊內(nèi)的核心基因。Hollender等[10]利用K-mean 聚類(lèi)方法和WGCNA 網(wǎng)絡(luò)分析方法，將30個(gè)樣本成千上萬(wàn)的基因劃分為12個(gè)組織和發(fā)育時(shí)期特異的表達(dá)模塊，鑒定到27個(gè)與花托發(fā)育相關(guān)的核心基因。鞠正等[11]結(jié)合RNA-Seq（RNA-Sequencing）數(shù)據(jù)和WGCNA 方法構(gòu)建番茄加權(quán)基因共表達(dá)網(wǎng)絡(luò)，并通過(guò)對(duì)果實(shí)成熟階段組織的特異性模塊分析，挖掘出58個(gè)潛在的與番茄果實(shí)成熟相關(guān)的轉(zhuǎn)錄因子。鮮小華等[12]結(jié)合甘藍(lán)型油菜種子發(fā)育過(guò)程中8個(gè)時(shí)期的轉(zhuǎn)錄組數(shù)據(jù)，采用WGCNA 方法構(gòu)建了8個(gè)共表達(dá)模塊，挖掘出與黃籽表型相關(guān)的微效作用基因。楊宇昕等[13]利用玉米14個(gè)不同組織的RNA-Seq 數(shù)據(jù)，通過(guò)WGCNA 方法鑒定到14個(gè)組織特異性模塊，鑒定了開(kāi)花相關(guān)的調(diào)控基因，為解析玉米發(fā)育的遺傳機(jī)理提供了重要的參考依據(jù)。

本研究利用棉花18個(gè)不同長(zhǎng)度果枝節(jié)間樣品的RNA-Seq 數(shù)據(jù)，通過(guò)WGCNA 方法構(gòu)建了加權(quán)基因共表達(dá)網(wǎng)絡(luò)，得到了13個(gè)基因共表達(dá)模塊，通過(guò)對(duì)模塊內(nèi)基因的富集分析鑒定到具有生物學(xué)意義的特異性模塊，預(yù)測(cè)到13個(gè)與棉花果枝節(jié)間伸長(zhǎng)調(diào)控相關(guān)的樞紐基因，并挖掘到80個(gè)潛在的候選基因，為棉花果枝節(jié)間的伸長(zhǎng)調(diào)控機(jī)理研究和棉花株型育種提供了有價(jià)值的數(shù)據(jù)支持和理論依據(jù)。

1 材料與方法

1.1 材料與數(shù)據(jù)來(lái)源

本研究所用的轉(zhuǎn)錄組數(shù)據(jù)來(lái)源于實(shí)驗(yàn)室前期的試驗(yàn)，該試驗(yàn)于2017年在河南省安陽(yáng)市白壁鎮(zhèn)中國(guó)農(nóng)業(yè)科學(xué)院棉花研究所試驗(yàn)基地進(jìn)行，試驗(yàn)材料為棉花品種魯棉研28（L28，長(zhǎng)果枝品種）和新陸中77（XLZ77，短果枝品種）。4月20日播種，株距28cm，行距70cm，行長(zhǎng)8 m，大田常規(guī)管理。于7月6日選取10 株長(zhǎng)勢(shì)一致的棉花對(duì)其倒一、倒二、倒三果枝的第一節(jié)間進(jìn)行取樣，3個(gè)生物學(xué)重復(fù)，共18組樣品，對(duì)其進(jìn)行RNA-Seq，得到各樣本的基因表達(dá)量數(shù)據(jù)，使用FPKM(Fragments per kilobase of exon per million fragments mapped)值來(lái)衡量基因的表達(dá)水平。

1.2 模塊的劃分

在進(jìn)行WGCNA 分析之前，需對(duì)本試驗(yàn)選用的基因集進(jìn)行篩選過(guò)濾，把在超過(guò)半數(shù)樣本中沒(méi)有表達(dá)的低質(zhì)量的基因去掉，從而提高網(wǎng)絡(luò)構(gòu)建的精度。利用R 軟件中的WGCNA（v1.47）包完成權(quán)重基因共表達(dá)網(wǎng)絡(luò)的構(gòu)建[6]。首先根據(jù)各樣本中基因差異表達(dá)規(guī)律進(jìn)行聚類(lèi)，以分析樣本關(guān)系。根據(jù)無(wú)尺度網(wǎng)絡(luò)原則確定軟閾值（Soft thresholding power），取相關(guān)系數(shù)達(dá)到平臺(tái)期時(shí)最小的Power 值作為后續(xù)分析的參數(shù)，同時(shí)統(tǒng)計(jì)在不同Power 值下，基因平均連通性的變化。根據(jù)基因間表達(dá)量的相關(guān)性構(gòu)建基因聚類(lèi)樹(shù)，利用動(dòng)態(tài)剪枝法識(shí)別共表達(dá)模塊，設(shè)定模塊最少基因數(shù)為50個(gè)，然后根據(jù)模塊特征值的相似度（0.75）對(duì)表達(dá)模式相近的模塊進(jìn)行合并。

將模塊基因在各個(gè)樣本中的表達(dá)模式用模塊特征值MEmagenta 來(lái)展示，橫坐標(biāo)為樣本，縱坐標(biāo)為模塊，用模塊特征值來(lái)繪制樣本表達(dá)模式熱圖。通過(guò)該圖可以找出與特定樣本顯著相關(guān)的模塊，從而后續(xù)可選擇相應(yīng)的模塊進(jìn)行進(jìn)一步的研究。

1.3 模塊基因的功能富集分析

對(duì)各模塊基因進(jìn)行GO 功能分析和KEGG Pathway 分析。通過(guò)WEGO（http://wego.genomics.org.cn/cgi-bin/wego/index.pl）進(jìn)行GO 分類(lèi)，結(jié)果顯示各模塊基因可以歸為分子功能（Molecular function）、生物過(guò)程（Biological process）和細(xì)胞組分（Cellular component）三個(gè)本體，通過(guò)GO 分析能確定差異表達(dá)基因行使的主要生物學(xué)功能。Pathway 顯著性富集分析以KEGG Pathway 為單位，應(yīng)用超幾何檢驗(yàn)，找出與整個(gè)基因組背景相比在差異表達(dá)基因中顯著性富集的Pathway，通過(guò)Pathway 顯著性富集能確定差異表達(dá)基因參與的最主要生化代謝途徑和信號(hào)轉(zhuǎn)導(dǎo)途徑。經(jīng)過(guò)多重檢驗(yàn)校正后，以FDR（False discovery rate）≤0.05 為閾值，滿(mǎn)足此條件的GO term 和Pathway定義為在差異表達(dá)基因中顯著富集的GO term和Pathway。

1.4 棉花果枝節(jié)間伸長(zhǎng)調(diào)控相關(guān)模塊的鑒定

根據(jù)每個(gè)模塊的GO 和KEGG Pathway 富集分析結(jié)果，并結(jié)合樣本表達(dá)模式熱圖篩選出與果枝節(jié)間長(zhǎng)度性狀顯著相關(guān)的目的模塊，根據(jù)對(duì)目的模塊的富集分析及模塊內(nèi)基因連通性（Gene connectivity）的篩選得到參與調(diào)控棉花果枝節(jié)間伸長(zhǎng)過(guò)程的樞紐基因。

1.5 基因互作網(wǎng)絡(luò)的構(gòu)建

篩選出模塊間相關(guān)系數(shù)排名靠前的目標(biāo)基因后，提取與其表達(dá)相關(guān)性排名前10的weight值對(duì)應(yīng)的基因，利用Cytoscape_3.3.0 進(jìn)行加權(quán)基因共表達(dá)網(wǎng)絡(luò)圖的繪制，網(wǎng)絡(luò)中每一個(gè)節(jié)點(diǎn)代表一個(gè)基因，邊代表基因相互間的調(diào)控關(guān)系?；蛘{(diào)控關(guān)系網(wǎng)絡(luò)圖能幫助我們準(zhǔn)確篩選潛在與目標(biāo)基因存在調(diào)控關(guān)系的候選基因，并且可以利用已知基因的功能預(yù)測(cè)未知基因功能。

2 結(jié)果與分析

2.1 兩棉花品種的株型差異

于2017年7月7日對(duì)兩個(gè)棉花品種的果枝第一節(jié)間長(zhǎng)度進(jìn)行測(cè)量，如圖1 所示，XLZ77的果枝節(jié)間長(zhǎng)度較短，自頂部向下各果枝節(jié)間長(zhǎng)度變化幅度較小，而L28 品種的果枝節(jié)間長(zhǎng)度顯著高于XLZ77，在倒6 果枝處長(zhǎng)度達(dá)到14cm，由此反映出兩品種株型結(jié)構(gòu)之間差異顯著。為進(jìn)一步解釋不同果枝節(jié)間長(zhǎng)度差異形成的原因，前期的試驗(yàn)中對(duì)兩品種棉花倒一、倒二和倒三果枝的第一節(jié)間進(jìn)行了轉(zhuǎn)錄組測(cè)序，每個(gè)節(jié)位取3 次生物學(xué)重復(fù)，共得到18份果枝節(jié)間樣品的RNA-Seq 數(shù)據(jù)。

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

對(duì)18個(gè)果枝節(jié)間樣品的RNA-Seq 數(shù)據(jù)進(jìn)行篩選，過(guò)濾掉在超過(guò)半數(shù)樣本中沒(méi)有表達(dá)的低質(zhì)量的基因后，獲得34 559個(gè)有效基因用于WGCNA 分析。利用WGCNA 包計(jì)算權(quán)重值，以使網(wǎng)絡(luò)符合無(wú)尺度網(wǎng)絡(luò)分布。根據(jù)圖2的結(jié)果，選擇軟閾值β＝8 構(gòu)建共表達(dá)網(wǎng)絡(luò)。根據(jù)基因間表達(dá)量的相關(guān)性構(gòu)建基因聚類(lèi)樹(shù)，如圖3 所示，樹(shù)的一個(gè)分支對(duì)應(yīng)一簇表達(dá)量高度相關(guān)的基因集。

圖1 兩棉花品種果枝第一節(jié)間長(zhǎng)度的變化Fig. 1 Changes in the first internode length of fruiting branches of two cotton varieties

圖2 軟閾值的確定Fig. 2 Determination of soft threshold power

利用動(dòng)態(tài)剪切樹(shù)法對(duì)分支進(jìn)行剪切區(qū)分，產(chǎn)生不同的模塊，根據(jù)模塊相似度（0.75）對(duì)表達(dá)模式相似的模塊進(jìn)行合并后的模塊劃分，最終獲得13個(gè)共表達(dá)模塊（圖4），每一種顏色代表一個(gè)模塊，灰色的模塊代表無(wú)法歸入任何一個(gè)模塊的基因。其中Green 模塊包含的基因數(shù)量最多，為11 004個(gè)基因；Firebrick4 模塊所含基因最少，為98個(gè)，各模塊基因數(shù)量見(jiàn)表1。

2.3 果枝節(jié)間伸長(zhǎng)調(diào)控相關(guān)模塊的篩選

根據(jù)對(duì)每個(gè)模塊GO 和KEGG Pathway的富集分析結(jié)果對(duì)13個(gè)模塊進(jìn)行篩選，過(guò)濾掉沒(méi)有顯著富集或僅富集到少數(shù)幾個(gè)GO Term 和Pathway的模塊，如Yellow 模塊僅顯著富集到GO:0071289細(xì)胞對(duì)鎳離子的反應(yīng)及GO:0010045 鎳陽(yáng)離子響應(yīng)兩條GO Term，并且未顯著富集到任何KEGG Pathway。而Antiquewhite 模塊未顯著富集到任何GO Term，KEGG 富集分析也僅顯著富集到Spliceosome 剪接體這一條通路。Darkgrey模塊顯著富集到GO: 0008643 碳水化合物運(yùn)輸和GO:0009959 負(fù)向重性?xún)蓷lGO Term，且僅顯著富集到Citrate cycle（TCA cycle）三羧酸循環(huán)這一條KEGG Pathway。此外，F(xiàn)irebrick4 這一模塊未富集任何KEGG 通路，而GO 富集分析共獲得23個(gè)GO Term，大多是蛋白質(zhì)修復(fù)、蛋白質(zhì)代謝等，與本研究方向無(wú)參考價(jià)值，故將此類(lèi)模塊排除掉。通過(guò)篩選最終得到6個(gè)模塊，分別是Tan、Plum、Cyan、Sienna3、Lightcoral、Green 模塊。

圖3 基因聚類(lèi)樹(shù)和模塊劃分Fig. 3 Gene cluster dendrograms and module division

表1 共表達(dá)模塊中基因數(shù)量分布Table 1 Distribution of number of genes in co-expression modules

2.4 共表達(dá)模塊的功能富集分析及目標(biāo)模塊的鑒定

為進(jìn)一步了解各模塊的生物學(xué)功能，本研究對(duì)各模塊的KEGG 通路和GO 富集分析中的生物過(guò)程進(jìn)行深度解析。其中Tan 模塊顯著富集到氮素代謝、半胱氨酸和蛋氨酸代謝、糖和葡萄糖醛酸轉(zhuǎn)換等6條KEGG 通路，并顯著富集到GO:0071554細(xì)胞壁組織或生物發(fā)生、GO: 0009809木質(zhì)素合成及代謝、GO:0071669 植物型細(xì)胞壁組織或生物發(fā)生等152條GO Term。Plum 模塊僅富集到氨基糖和核苷酸代謝一條KEGG Pathway，以及GO:0009723 乙烯反應(yīng)、GO:0009725激素應(yīng)答、GO: 0009753 茉莉酸反應(yīng)、GO:0009873 乙烯介導(dǎo)的信號(hào)通路等59條GO Term。Sienna3 模塊顯著富集到苯丙素的生物合成、苯丙氨酸代謝、脂肪酸伸長(zhǎng)等8條KEGG 通路，并顯著富集到GO:0042546細(xì)胞壁生物發(fā)生、GO:0009698 苯丙烷代謝、GO: 0010410 半纖維素代謝、GO:0010413 葡糖醛酸木聚糖代謝等166條GO Term。Lightcoral 模塊中顯著富集到核苷酸切除修復(fù)、嘧啶代謝、嘌呤代謝、堿基修復(fù)、花生四烯酸代謝等KEGG 通路，且共有GO: 0006259 DNA 代謝、GO: 0051276染色質(zhì)組裝、GO:0006281 DNA 修復(fù)、GO: 0016568染色質(zhì)修復(fù)、GO: 0007049細(xì)胞周期、GO: 0016570組蛋白修飾等649條GO Term 得到了顯著富集。而Green模塊則有核糖體、光合作用、氧化磷酸化等4條KEGG 通路得到顯著富集，并有GO:0006334 核小體裝配、GO: 0031497染色質(zhì)組裝、GO:0009735細(xì)胞分裂素反應(yīng)、GO: 0044843細(xì)胞周期G1/S 相變等228條GO Term 得到顯著富集。其中，Cyan 模塊不僅顯著富集到亞麻酸代謝、植物激素信號(hào)轉(zhuǎn)導(dǎo)、類(lèi)固醇生物合成等7條KEGG通路，還顯著富集到GO:0009725 激素反應(yīng)、GO:0009723 乙烯反應(yīng)、GO: 0009753 茉莉酸反應(yīng)、GO: 0009864 茉莉酸介導(dǎo)的信號(hào)通路等414條GO Term。

鑒于Cyan 模塊的KEGG Pathway 通路富集分析及GO Category 富集分析均顯著富集得到植物激素信號(hào)調(diào)控這一通路，并結(jié)合樣本表達(dá)模式熱圖（圖4）可知Cyan 模塊基因在兩品種3個(gè)果枝節(jié)間中的表達(dá)量呈逐步升高趨勢(shì)，且XLZ77品種基因的表達(dá)量顯著高于L28 品種，其表達(dá)趨勢(shì)相較于其他模塊更有規(guī)律可循，因此將Cyan模塊作為目標(biāo)模塊以進(jìn)行下一步對(duì)樞紐基因的挖掘與分析，如圖5 所示為該模塊的KO 富集氣泡圖。

圖4 樣本表達(dá)模式熱圖Fig. 4 Sample expression heatmap

圖5 Cyan 模塊的KEGG 富集分析氣泡圖Fig. 5 Analysis of KEGG enrichment in Cyan module

2.5 樞紐基因的挖掘及互作網(wǎng)絡(luò)的構(gòu)建

模塊內(nèi)基因的連通性（Connectivity）代表該基因與其他基因之間的調(diào)控關(guān)系，連通性越高，說(shuō)明該基因在模塊中的調(diào)控作用越大，越有可能成為潛在的樞紐基因。根據(jù)這個(gè)原理，首先對(duì)Cyan 模塊KEGG Pathway 富集得到的植物激素信號(hào)轉(zhuǎn)導(dǎo)通路中的121個(gè)基因以及GO Category富集分析得到的GO: 0009725 激素反應(yīng)通路中的649個(gè)基因按照連通性進(jìn)行降序排序，發(fā)現(xiàn)有13個(gè)編碼茉莉酸ZIM 結(jié)構(gòu)域蛋白的JAZ（Jasmonate-zim-domain protein）基因的連通性顯著高于所在通路中的其他基因。對(duì)于Cyan 模塊所有基因而言，該類(lèi)JAZ基因的連通性排名也非?？壳?，其中有8個(gè)JAZ基因的連通性在該模塊內(nèi)排名前10%，因此綜合KEGG Pathway 富集分析、GO Category 富集分析和模塊內(nèi)基因連通性排序三個(gè)方面的數(shù)據(jù)結(jié)果，認(rèn)定這13個(gè)JAZ基因作為本研究的樞紐基因（表2），即在棉花果枝節(jié)間伸長(zhǎng)過(guò)程中發(fā)揮著重要調(diào)控作用的核心基因。

為了進(jìn)一步挖掘參與棉花果枝節(jié)間伸長(zhǎng)調(diào)控的候選基因，本研究以選定的13個(gè)JAZ類(lèi)基因作為樞紐基因，篩選與其表達(dá)相關(guān)性排名前10的權(quán)重值對(duì)應(yīng)的基因作為調(diào)控棉花果枝節(jié)間伸長(zhǎng)過(guò)程的候選基因，利用13個(gè)核心基因及其對(duì)應(yīng)的80個(gè)候選基因通過(guò)Cytoscape_3.3.0 進(jìn)行加權(quán)基因共表達(dá)網(wǎng)絡(luò)調(diào)控圖的繪制（圖6），網(wǎng)絡(luò)中每一個(gè)節(jié)點(diǎn)代表一個(gè)基因，連接線代表基因相互間的調(diào)控關(guān)系，處于連接線兩端的基因被認(rèn)為具有相同的生物學(xué)功能。通過(guò)基因功能注釋發(fā)現(xiàn)這80個(gè)候選基因分別被注釋為WRKY家族轉(zhuǎn)錄因子、WRKYDNA 結(jié)合蛋白、蛋白激酶超家族蛋白、木葡聚糖水解酶等，而本研究通過(guò)基因共表達(dá)網(wǎng)絡(luò)調(diào)控圖發(fā)現(xiàn)這80個(gè)基因與樞紐基因存在緊密的聯(lián)系，故推論這80個(gè)候選基因可能與JAZ類(lèi)蛋白存在功能上的相似性，繼而通過(guò)協(xié)同調(diào)控作用來(lái)調(diào)控棉花果枝節(jié)間的伸長(zhǎng)生長(zhǎng)過(guò)程。

3 討論

3.1 WGCNA 分析的重要意義

在大樣本條件下，傳統(tǒng)的趨勢(shì)分析無(wú)法對(duì)基因進(jìn)行有效的分類(lèi)，不能充分有效地挖掘海量數(shù)據(jù)中蘊(yùn)含的生物學(xué)意義，依賴(lài)數(shù)據(jù)的功能分析亦無(wú)法推測(cè)新的調(diào)控關(guān)系。而新興的共表達(dá)網(wǎng)絡(luò)分析方法的出現(xiàn)則能打破傳統(tǒng)趨勢(shì)分析的局限，它可以將復(fù)雜的數(shù)據(jù)進(jìn)行歸納和整理，高效研究基因整體表達(dá)規(guī)律，同時(shí)能夠系統(tǒng)地反饋樣本中的基因間相互作用模式。WGCNA 分析的意義在于通過(guò)基因的模塊分類(lèi)，呈現(xiàn)樣本中基因的全局表達(dá)規(guī)律，尋找與樣本性狀相關(guān)的基因模塊和核心基因，并利用已知基因推測(cè)未知基因的功能，為下游生物學(xué)實(shí)驗(yàn)提供線索。因此，WGCNA 在分析多樣本的RNA-Seq 數(shù)據(jù)中發(fā)揮著重要的作用。

表2 Cyan 模塊中核心基因的功能注釋Table 2 Annotation of hub genes in Cyan module

圖6 Cyan 模塊果枝節(jié)間伸長(zhǎng)調(diào)控相關(guān)基因共表達(dá)網(wǎng)絡(luò)圖Fig. 6 Co-expression network of genes related to internode elongation regulation of cotton fruiting branches in Cyan module

3.2 棉花果枝節(jié)間伸長(zhǎng)相關(guān)共表達(dá)模塊的鑒定

本研究針對(duì)XLZ77 和L28 兩品種果枝節(jié)間長(zhǎng)度差異顯著這一性狀展開(kāi)分析，通過(guò)對(duì)18個(gè)不同節(jié)位果枝產(chǎn)生的34 559個(gè)有效差異基因進(jìn)行模塊劃分，共獲得13個(gè)共表達(dá)模塊，模塊內(nèi)的基因具有相同或相似的表達(dá)趨勢(shì)，而模塊間的表達(dá)趨勢(shì)各不相同。通過(guò)對(duì)各模塊KEGG Pathway和GO Category 富集分析結(jié)果的篩選，排除掉Yellow、Antiquewhite、Darkgrey、Grey等7個(gè)未顯著富集或極少富集到通路的無(wú)參考意義的模塊，最終獲得Tan、Plum、Cyan、Sienna3、Lightcoral、Green 6個(gè)模塊。鑒于Cyan 模塊無(wú)論在KEGG Pathway 還是在GO Category 富集分析中，均有植物激素信號(hào)轉(zhuǎn)導(dǎo)通路得到了顯著富集，且該模塊內(nèi)基因的表達(dá)趨勢(shì)在兩品種內(nèi)具有明顯且有規(guī)律的差異，因此本研究將Cyan 模塊作為目標(biāo)模塊，認(rèn)定該模塊基因的差異表達(dá)可能是導(dǎo)致兩品種果枝節(jié)間長(zhǎng)度差異的主要原因。

3.3 樞紐基因JAZ 與候選基因間互作調(diào)控棉花果枝節(jié)間的伸長(zhǎng)

根據(jù)對(duì)Cyan 模塊所有基因及KEGG 激素信號(hào)轉(zhuǎn)導(dǎo)通路和GO 富集激素反應(yīng)通路中基因連通性的排序，本研究發(fā)現(xiàn)共有13個(gè)編碼JAZ的基因在KEGG 和GO 分析中均得到了顯著富集，且有8個(gè)JAZ基因在該模塊內(nèi)連通性排名前10 %，故認(rèn)定這13個(gè)JAZ類(lèi)基因作為Cyan 模塊的樞紐基因。已知JAZ蛋白作為茉莉酸響應(yīng)的抑制因子，在高水平JA-Ile 存在時(shí)可促進(jìn)JAZ蛋白和冠菌素不敏感因子復(fù)合體SCFCOI1結(jié)合使JAZ蛋白泛素化，泛素化的JAZ蛋白由26s 蛋白酶體途徑降解，從而釋放轉(zhuǎn)錄因子MYC2，啟動(dòng)JA 應(yīng)答基因的轉(zhuǎn)錄[14]。而在前人的研究中曾指出茉莉酸抑制植物的生長(zhǎng)[15]，能夠抑制細(xì)胞周期相關(guān)蛋白的表達(dá)從而抑制細(xì)胞分裂[16]。也有人指出茉莉酸通過(guò)影響細(xì)胞分裂和細(xì)胞伸長(zhǎng)而抑制擬南芥主根生長(zhǎng)[17]，GhJAZ通過(guò)抑制GhMYB25-like的轉(zhuǎn)錄激活功能而抑制下游基因的表達(dá)，從而抑制了棉花纖維的突起[18]。擬南芥JAZ7 過(guò)量表達(dá)突變體植株地上部與野生型相比表現(xiàn)出生長(zhǎng)矮小的性狀，而番茄JAZ7 過(guò)表達(dá)轉(zhuǎn)基因植株與野生型相比，也表現(xiàn)出株高矮小、幼葉莖稈卷縮、果實(shí)發(fā)育進(jìn)程緩慢的特點(diǎn)[19]。此外，在本實(shí)驗(yàn)室前期工作中，我們測(cè)定了兩個(gè)品種果枝節(jié)間內(nèi)的激素含量，結(jié)果顯示XLZ77 內(nèi)茉莉酸含量顯著高于L28。而利用13個(gè)JAZ基因作為Cyan 模塊的核心基因，選取了與其相關(guān)性排名前10的共80個(gè)基因作為候選基因，進(jìn)行基因互作調(diào)控圖的繪制，發(fā)現(xiàn)這80個(gè)基因被注釋為WRKY家族轉(zhuǎn)錄因子、WRKYDNA 結(jié)合蛋白、蛋白激酶超家族蛋白、木葡聚糖水解酶等，根據(jù)共表達(dá)網(wǎng)絡(luò)基因的生物學(xué)功能相似性原理可以認(rèn)定候選基因與樞紐基因具有相同或相似的功能，因此推測(cè)本研究中XLZ77 品種內(nèi)較高的茉莉酸能夠激活JAZ基因的表達(dá)，而高表達(dá)的JAZ蛋白可能通過(guò)與候選基因的協(xié)同調(diào)控作用介導(dǎo)了茉莉酸反應(yīng)，從而減弱了棉花果枝節(jié)間細(xì)胞的分裂與伸長(zhǎng)，減緩了棉花的營(yíng)養(yǎng)生長(zhǎng)，出現(xiàn)果枝節(jié)間較短的現(xiàn)象。

以上結(jié)果表明，WGCNA 分析能夠?qū)?fù)雜的數(shù)據(jù)劃分為具有共表達(dá)關(guān)系的模塊，更能反映真實(shí)生物網(wǎng)絡(luò)中的有效信息，能夠更好地挖掘與棉花果枝節(jié)間伸長(zhǎng)相關(guān)的核心基因，并預(yù)測(cè)潛在候選基因的功能，為研究不同性狀的基因調(diào)控網(wǎng)絡(luò)提供新的研究思路。

4 結(jié)論

利用權(quán)重基因共表達(dá)網(wǎng)絡(luò)分析（WGCNA）的方法將34 559個(gè)差異基因劃分成13個(gè)共表達(dá)模塊，完成了模塊內(nèi)基因的功能富集分析，揭示了模塊所蘊(yùn)含的生物學(xué)意義，鑒定到Cyan 模塊是與棉花果枝節(jié)間伸長(zhǎng)相關(guān)的特異性模塊，并篩選得到植物激素信號(hào)轉(zhuǎn)導(dǎo)通路中的13個(gè)JAZ基因?yàn)樵撃K內(nèi)的樞紐基因。利用樞紐基因挖掘到80個(gè)潛在的候選基因，并用其構(gòu)建了基因互作調(diào)控網(wǎng)絡(luò)。本研究結(jié)果可為棉花果枝節(jié)間伸長(zhǎng)的調(diào)控機(jī)理研究提供數(shù)據(jù)支持，為進(jìn)一步調(diào)控棉花株型結(jié)構(gòu)奠定理論基礎(chǔ)。