周雨青， 楊永飛， 葛常偉， 沈倩， 張思平， 劉紹東， 馬慧娟，陳靜， 劉瑞華， 李士叢， 趙新華， 李存東， 龐朝友*

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

棉花是一種喜溫作物，起源于熱帶和亞熱帶地區(qū)，近年來，我國主要棉區(qū)逐步“西進(jìn)、東移、北上”［1］，其栽培范圍已逐步擴(kuò)展至較冷區(qū)域。新疆植棉面積較大，該地區(qū)屬于典型的大陸性氣候，春季經(jīng)常發(fā)生“倒春寒”，因而在播種出苗期及苗期常常遭受低溫的危害，其中苗期冷害發(fā)生較頻繁，大幅度降溫、下霜并伴隨降水的強(qiáng)烈災(zāi)害性天氣對(duì)棉苗生長危害極大［2］，每年因爛種、爛芽、死苗或晚發(fā)而導(dǎo)致部分棉田重播［3］，嚴(yán)重影響了我國棉花的產(chǎn)量及纖維品質(zhì)。因此，揭示棉花苗期抗冷機(jī)制和培育抗冷棉花品種具有重要意義。

加權(quán)基因表達(dá)網(wǎng)絡(luò)分析（weighted gene coexpression network analysis，WGCNA）是基于大樣本轉(zhuǎn)錄組數(shù)據(jù)的生物信息學(xué)分析方法，其首先假定基因網(wǎng)絡(luò)服從無尺度分布，將表達(dá)模式相似的基因聚類形成不同的模塊，分析模塊與特定性狀或表型之間的關(guān)聯(lián)性，通過基因聚類和關(guān)聯(lián)分析結(jié)果構(gòu)建共表達(dá)調(diào)控網(wǎng)絡(luò)。位于調(diào)控網(wǎng)絡(luò)中心的基因被稱為核心基因（hub gene），這類基因通常是關(guān)鍵的調(diào)控基因，值得深入挖掘和分析［4］。在棉花研究中，巨飛燕等［5］利用2個(gè)株型結(jié)構(gòu)差異顯著的棉花品種的18份不同長度果枝節(jié)間轉(zhuǎn)錄組數(shù)據(jù)構(gòu)建共表達(dá)網(wǎng)絡(luò)，鑒定到與棉花果枝節(jié)間伸長相關(guān)的特異性模塊，并發(fā)現(xiàn)植物激素信號(hào)轉(zhuǎn)導(dǎo)通路中的JAZ基因?yàn)樵撃K的樞紐基因。傅明川等［6］利用21份黃萎病菌侵染不同時(shí)間點(diǎn)的海島棉幼苗根尖轉(zhuǎn)錄組數(shù)據(jù)構(gòu)建共表達(dá)網(wǎng)絡(luò)，鑒定到5個(gè)與抗黃萎病相關(guān)特異性模塊，并挖掘出了網(wǎng)絡(luò)中的核心基因。在抗冷研究中，李旭凱等［7］利用47份正常水稻組織轉(zhuǎn)錄組數(shù)據(jù)，通過冷脅迫、干旱脅迫、鹽脅迫不同的處理方式，使用WGCNA方法挖掘到2 599個(gè)與3種脅迫都相關(guān)的基因，并預(yù)測出25個(gè)抗逆相關(guān)的關(guān)鍵基因，為水稻的綜合抗逆能力等研究提供了新思路。秦夢(mèng)凡等［8］利用2個(gè)抗凍響應(yīng)有差別的甘藍(lán)型油菜材料經(jīng)不同低溫處理的36個(gè)轉(zhuǎn)錄組數(shù)據(jù)，鑒定到了共同響應(yīng)凍害和耐寒的特異性模塊，并對(duì)其調(diào)控機(jī)制進(jìn)行了分析，為油菜的耐寒調(diào)控機(jī)制研究提供重要的參考依據(jù)。而關(guān)于棉花抗冷研究的WGCNA分析未見報(bào)道。

本研究以4℃低溫處理不同時(shí)間的棉花子葉轉(zhuǎn)錄組數(shù)據(jù)為材料，對(duì)其進(jìn)行差異表達(dá)分析；通過構(gòu)建加權(quán)基因共表達(dá)網(wǎng)絡(luò)劃分基因模塊，篩選出抗冷相關(guān)的特異性模塊；經(jīng)GO富集分析探究模塊功能，根據(jù)基因在相應(yīng)網(wǎng)絡(luò)中的連通性鑒定出模塊內(nèi)的核心基因。本研究可為進(jìn)一步理解棉花低溫冷害的分子機(jī)制提供理論基礎(chǔ)，并為棉花抗冷育種提供新的基因資源。

1 材料與方法

1.1 試驗(yàn)材料及處理

試驗(yàn)材料新陸中16為抗冷品種（cold tolerance，CT），新陸中 32為冷敏感品種（cold sensitivity，CS），均由中國農(nóng)業(yè)科學(xué)院棉花研究所種質(zhì)資源中期庫提供。棉花材料先于28℃光照培養(yǎng)室進(jìn)行育苗（光照16 h，黑暗8 h），采用營養(yǎng)土與蛭石體積3∶1的混合基質(zhì)培養(yǎng)，并保證水分充足，待子葉平展時(shí)將材料放入4℃冷室處理0、1、3、6、9、12 h，選取長勢一致的棉苗對(duì)其子葉進(jìn)行取樣，每個(gè)時(shí)期3次重復(fù)，每次重復(fù)子葉數(shù)量為8～10片，共36個(gè)樣本，對(duì)其進(jìn)行轉(zhuǎn)錄組測序，得到各樣本的基因表達(dá)量數(shù)據(jù)，使用FPKM（fragments per kilobase of exon per million fragments mapped）值來衡量基因的表達(dá)水平［9］。新陸中16和新陸中32的測序樣品組設(shè)為CT和CS，根據(jù)冷處理時(shí)間0、1、3、6、9、12 h，新陸中16處理組命名為 CT0、CT1、CT3、CT6、CT9、CT12，新陸中 32 處理組命名為CS0、CS1、CS3、CS6、CS9、CS12。

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

對(duì)選用的基因集進(jìn)行篩選過濾，需同時(shí)滿足以下條件：基因在樣本中的表達(dá)量均值要大于1，超過一半的樣本表達(dá)量大于0，變異系數(shù)大于0.2。去掉不符合條件的低質(zhì)量基因，提高網(wǎng)絡(luò)構(gòu)建的精度。利用R軟件中的WGCNA（v1.47）包計(jì)算權(quán)重值，完成權(quán)重基因共表達(dá)網(wǎng)絡(luò)的構(gòu)建。首先根據(jù)所有基因的表達(dá)量對(duì)所有樣本進(jìn)行聚類，以分析樣本關(guān)系。根據(jù)無尺度網(wǎng)絡(luò)原則確定軟閾值（soft thresholding power），取相關(guān)系數(shù)達(dá)到平臺(tái)期（或大于0.8）時(shí)最小power值作為后續(xù)分析參數(shù)，同時(shí)統(tǒng)計(jì)在不同power值下基因平均連通性的變化。采用動(dòng)態(tài)切割法（dynamic tree cut）對(duì)基因進(jìn)行聚類及模塊劃分。根據(jù)基因間表達(dá)量的相關(guān)性構(gòu)建基因聚類樹，并根據(jù)基因間的聚類關(guān)系劃分基因模塊，然后根據(jù)模塊特征值的相似度合并表達(dá)模式相近的模塊。設(shè)定模塊最少基因數(shù)為50，相似模塊合并閾值為0.8，將模塊基因在各個(gè)樣本中的表達(dá)模式用模塊特征值來展示，并繪制樣本表達(dá)模式熱圖。通過模塊特征向量基因（module eigengene，ME）分析確定與抗冷顯著相關(guān)的特異性模塊，后續(xù)選擇相應(yīng)的模塊進(jìn)行深入研究。

1.3 核心基因的篩選和功能分析

為進(jìn)一步挖掘特異性模塊的功能，對(duì)特異性模塊基因進(jìn)行GO功能分析。通過GO數(shù)據(jù)庫（http：//www.geneontology.org/）進(jìn)行GO分類，結(jié)果顯示目標(biāo)模塊基因可以歸為分子功能（molecular function）、生物過程（biological process）和細(xì)胞組分（cellular component）3個(gè)大類。經(jīng)過多重檢驗(yàn)校正后，以P值＜0.05為閾值，滿足此條件的定義為在基因中顯著富集的GO term。

為獲得特異性模塊中的核心基因，利用Cytoscape3.7.1軟件對(duì)基因互作網(wǎng)絡(luò)進(jìn)行可視化處理，根據(jù)模塊中基因的連通性（connectivity）及轉(zhuǎn)錄組數(shù)據(jù)中基因的FPKM值，篩選出連通性高且表達(dá)量高的基因，將這些基因作為該模塊中的核心基因。在這些網(wǎng)絡(luò)中，每個(gè)節(jié)點(diǎn)代表1個(gè)基因，處于連線兩端的基因通常被認(rèn)為具有相同的生物功能?；蛘{(diào)控關(guān)系網(wǎng)絡(luò)圖能準(zhǔn)確篩選與核心基因存在調(diào)控關(guān)系的候選基因，并利用已知基因的功能預(yù)測未知基因功能。

1.4 核心基因的qRT-PCR驗(yàn)證

將用于轉(zhuǎn)錄組測序的備份RNA反轉(zhuǎn)錄，利用SYBR PremixEx Taq（DRR041A）熒光定量試劑盒對(duì)篩選得到的核心基因進(jìn)行qRT-PCR（quantitative real-time polymerase chain reaction）驗(yàn)證，檢測核心基因在CT和CS中對(duì)冷脅迫的響應(yīng)情況。通過網(wǎng)站qPrimerDB（https：//biodb.swu.edu.cn/qprimerdb/）設(shè)計(jì)qRT-PCR引物（表1）。內(nèi)參基因?yàn)镚hUBQ7（Gh_A11G0969），每個(gè)反應(yīng)設(shè)置3個(gè)重復(fù)。用2-ΔΔCT方法［10］計(jì)算基因的相對(duì)表達(dá)量。

表1 qRT-PCR所用基因引物序列Table 1 Specific primers for the selected genes

2 結(jié)果與分析

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

2.1.1 軟閾值確定 從圖1可以看出，當(dāng)軟閾值β=8時(shí)，無尺度網(wǎng)絡(luò)擬合指數(shù)R2＞0.8，平均連通性趨近于0，表明用此值進(jìn)行冪處理可以得到符合要求的無尺度網(wǎng)絡(luò)，因此選擇β=8構(gòu)建無尺度網(wǎng)絡(luò)。

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

2.1.2 基因共表達(dá)網(wǎng)絡(luò)中基因的聚類 根據(jù)基因間表達(dá)量的相關(guān)性構(gòu)建聚類樹，采用動(dòng)態(tài)切割法將產(chǎn)生的聚類樹切割，把表達(dá)模式相似的基因合并在同一分支上，每個(gè)分支代表1個(gè)共表達(dá)模塊，根據(jù)模塊相似度（0.8）對(duì)表達(dá)模式相似的模塊合并后進(jìn)行模塊劃分，最終獲得9個(gè)共表達(dá)模塊（圖2），每一種顏色代表一個(gè)模塊，灰色（Grey）模塊代表無法歸入任何一個(gè)模塊的基因。藍(lán)色（Blue）模塊包含的基因數(shù)量最多，為16 025個(gè)，其次是棕色（Brown）模塊，包含1 893個(gè)基因，黃色（Yellow）模塊包含1 818個(gè)基因，灰色（Grey）模塊包含655個(gè)基因，黑色（Black）模塊包含616個(gè)基因，綠色（Green）模塊包含489個(gè)基因，洋紅色（Magenta）模塊包含337個(gè)基因，粉色（Pink）模塊包含184個(gè)基因，黃綠色（Greenyellow）模塊所含基因最少，為66個(gè)（圖3）。

圖2 基因聚類樹和模塊劃分Fig.2 Gene cluster dendrograms and module division

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

2.1.3 基因共表達(dá)網(wǎng)絡(luò)中模塊與樣本的關(guān)聯(lián)分析 將所獲模塊與樣本進(jìn)行關(guān)聯(lián)分析，得到9個(gè)與不同品種和處理時(shí)間相關(guān)聯(lián)的模塊，部分模塊與品種和處理時(shí)間高度關(guān)聯(lián)。例如，Brown模塊的基因表達(dá)模式隨著低溫處理時(shí)間的增加在2個(gè)品種中均存在由負(fù)到正的趨勢，與低溫處理9和12 h時(shí)間點(diǎn)正相關(guān)，與低溫處理0、1、3和6 h時(shí)間點(diǎn)負(fù)相關(guān)；Blue模塊的基因表達(dá)模式隨著低溫處理時(shí)間的增加在兩個(gè)品種中存在由正到負(fù)的趨勢，與低溫處理0、1和3 h時(shí)間點(diǎn)正相關(guān)，與低溫處理6、9和12 h時(shí)間點(diǎn)負(fù)相關(guān)（圖4）。

圖4 基因共表達(dá)網(wǎng)絡(luò)模塊與不同樣本的關(guān)聯(lián)熱圖Fig.4 Association analysis of gene co-expression network modules with different samples

模塊中具有代表性的基因即為特征向量基因（ME），它們之間的相關(guān)性越高，其所在模塊的相關(guān)性也就越高，通過ME相關(guān)性分析發(fā)現(xiàn)，Blue模塊與Brown模塊中ME之間的相關(guān)性達(dá)到-0.86（圖5）。因此將這兩個(gè)模塊作為抗冷相關(guān)特異性模塊進(jìn)行深入研究，挖掘模塊中的核心基因。

圖5 不同模塊兩兩之間ME的相關(guān)性Fig.5 ME correlation between different modules

2.2 核心基因的篩選和功能分析

對(duì)Blue和Brown模塊進(jìn)行GO功能富集分析發(fā)現(xiàn)，這2個(gè)模塊都可以顯著富集到生物學(xué)過程（P）、分子功能（F）以及細(xì)胞組分（C）的若干GO通路（表2）。在Blue模塊中篩選出了5個(gè)核心基因（圖6），Brown模塊中篩選出了5個(gè)核心基因（圖7），利用棉花基因數(shù)據(jù)庫（https：//cottonfgd.org/）和NCBI數(shù)據(jù)庫（https：//www.ncbi.nlm.nih.gov/）獲取這些核心基因的功能信息，并借助TAIR數(shù)據(jù)庫（https：//www.arabidopsis.org/）注釋這些核心基因在擬南芥中的同源基因的功能（表3），結(jié)合已有研究進(jìn)展，本研究發(fā)現(xiàn)這些篩選出的核心基因基本都與冷脅迫相關(guān)，比如Blue模塊中的C2H2型鋅指蛋白Gh_A13G2112富集到刺激反應(yīng)（GO：0050896）、幾丁質(zhì)反應(yīng)（GO：0010200）、冷反應(yīng)（GO：0009409）、溫度刺激反應(yīng)（GO：0009266）、脅迫反應(yīng)（GO：0006950）等多個(gè)通路中；Brown模塊中參與冷調(diào)節(jié)的基因Gh_D09G1773和Gh_A05G1 554均富集到刺激反應(yīng)（GO：0050896）、幾丁質(zhì)反應(yīng)（GO：0010200）、溫度刺激反應(yīng)（GO：0009266）、冷反應(yīng)（GO：0009409）、外部刺激反應(yīng)（GO：0009605）、內(nèi)源性刺激反應(yīng)（GO：0009719）等多個(gè)通路。

圖6 Blue模塊的基因共表達(dá)網(wǎng)絡(luò)及其核心基因Fig.6 Gene co-expression network and hub genes in Blue module

圖7 Brown模塊的基因共表達(dá)網(wǎng)絡(luò)及其核心基因Fig.7 Gene co-expression network and hub genes in Brown module

表2 特異性模塊的部分GO富集分析結(jié)果Table 2 Partial GO enrichment analysis of target module

表2 特異性模塊的部分GO富集分析結(jié)果Table 2 Partial GO enrichment analysis of target module 續(xù)表Continued

2.3 核心基因的RT-qPCR驗(yàn)證

從相關(guān)性最高的2個(gè)基因模塊中篩選到10個(gè)可能與冷脅迫相關(guān)的核心基因（表3），對(duì)這10個(gè)基因進(jìn)行qRT-PCR驗(yàn)證，最終基因表達(dá)量的變化趨勢與轉(zhuǎn)錄組結(jié)果基本一致（圖8和9），進(jìn)一步證明了轉(zhuǎn)錄組測序結(jié)果的可靠性。這10個(gè)核心基因的表達(dá)量在2個(gè)抗冷性差異品種中隨著低溫處理時(shí)間的增加而增加，均呈現(xiàn)上調(diào)表達(dá)的趨勢，驗(yàn)證了這些核心基因響應(yīng)冷脅迫。

圖8 核心基因的表達(dá)模式Fig.8 Expression pattern of hub genes

表3 目標(biāo)模塊中核心基因的功能注釋Table 3 Annotation of hub genes in target module

圖9 核心基因的RT-qPCR驗(yàn)證Fig.9 qRT-PCR validation of hub genes

3 討論

本研究通過WGCNA富集到9個(gè)冷脅迫處理下與棉花品種和處理時(shí)間相關(guān)聯(lián)的共表達(dá)模塊，利用各個(gè)模塊的特征向量基因?qū)δK進(jìn)行相關(guān)性分析，發(fā)現(xiàn)Brown和Blue模塊間存在極顯著的相關(guān)性且關(guān)聯(lián)度最高，可作為子葉抗冷性狀機(jī)理研究的目標(biāo)模塊。對(duì)這2個(gè)模塊內(nèi)的基因進(jìn)行GO功能富集分析，均富集到了抗逆相關(guān)的具有生物學(xué)意義的調(diào)控途徑。例如，Blue和Brown模塊富集到轉(zhuǎn)錄因子活性（GO：0003700）、對(duì)脅迫的反應(yīng)（GO：0006950）等通路，該結(jié)果與謝勇軍［11］的研究結(jié)果一致。

對(duì)這2個(gè)模塊的核心基因進(jìn)行分析發(fā)現(xiàn)，核心基因里也有部分基因與冷脅迫相關(guān)，如Gh_A05G1931、 Gh_A13G2112、 Gh_A12G2357、Gh_D09G1773、Gh_A05G1554。篩選出的核心基因在模式植物擬南芥和其他作物（如玉米、水稻等）中有較多的研究，但在棉花中未見報(bào)道，可進(jìn)行深入探究。Gh_A05G1931（AT4G08950）參與了對(duì)油菜素類固醇（brassinosteroids，BR）的反應(yīng)。遏藍(lán)菜（Thlaspi arvense）比其近緣植物擬南芥（Arabidopsis thaliana）或甘藍(lán)型油菜（Brassicanapus）具有更高的抗凍性，原因可能是一些基因在表達(dá)水平和表達(dá)時(shí)間上存在差異，試驗(yàn)證明EXO（AT4G08950）基因在甘藍(lán)型油菜和遏藍(lán)菜中均顯著上調(diào)表達(dá)，但在擬南芥中沒有觀察到變化，且在甘藍(lán)型油菜中表達(dá)量的增加最為明顯［12］，說明EXO基因能增強(qiáng)遏藍(lán)菜的抗凍性。玉米中，用油菜素類固醇進(jìn)行浸種處理，可增加幼苗可溶性糖含量，減輕細(xì)胞膜傷害，提高脯氨酸含量，從而增強(qiáng) 植 物 對(duì) 低 溫 的 抗 性［13］。 Gh_A13G2112（AT1G27730）屬于C2H2型鋅指蛋白家族。C2H2鋅指蛋白含有ERF相關(guān)的兩親性抑制（EAR）結(jié)構(gòu)域［14-16］，在調(diào)節(jié)植物對(duì)非生物脅迫的防御反應(yīng)中起著關(guān)鍵作用。ZAT10/STZ（AT1G27730）最初被鑒定為鹽和冷反應(yīng)蛋白［17］。ZAT10轉(zhuǎn)錄受脫落酸、干旱、氯化鈉、低溫和強(qiáng)光等多種非生物脅迫的高度誘導(dǎo)。過表達(dá)ZAT10的轉(zhuǎn)基因植物表現(xiàn)出生長遲緩，同時(shí)它們表現(xiàn)出對(duì)多種非生物脅迫的耐受性增強(qiáng)，并且敲除該基因也可以增強(qiáng)對(duì)鹽、滲透脅迫的耐受性［15，17-20］。此外，ZAT10為擬南芥絲裂原活化蛋白激酶（MAPKs或MPKs）的直接靶向底物，其活性可能通過MPK直接磷酸化來調(diào)節(jié)［21］。植物MPKs參與了對(duì)各種生物和非生物脅迫的脅迫耐受反應(yīng)，如病原體感染、低溫、低濕度、鹽、傷害、紫外線、臭氧和重金屬［22-27］。Gh_A12G2357（AT5G51990）是C-重復(fù)結(jié)合因子/脫水反應(yīng)結(jié)合元件（CBF/DREB1）家族的成員。CBF基因家族主要包括4個(gè)成員：CBF1、CBF2、CBF3和CBF4［28］。冷處理下CBF1和CBF3基因的表達(dá)既不上調(diào)也不下調(diào)，但激活了CBF2和CBF4（AT5G51990）基因，尤其是低溫脅迫下CBF4基因的表達(dá)高度上調(diào)［29］。CBF4過表達(dá)使轉(zhuǎn)基因擬南芥更耐寒和耐旱［30］。Gh_D09G1773（AT3G49530）是植物特異轉(zhuǎn)錄因子NAC家族的成員，在發(fā)育過程和脅迫反應(yīng)中均起作用［31］。研究發(fā)現(xiàn)，NAC062（AT3G49530）受冷誘導(dǎo)活化，進(jìn)入細(xì)胞核調(diào)控抗病相關(guān)基因，提高植物抗病能力，冷害脅迫期間NAC062的活化可能由膜組成和流動(dòng)性改變導(dǎo)致［32］。Gh_A05G1554（AT1G20440）屬于脫水蛋白家族，與脫落酸途徑有關(guān)，當(dāng)其過表達(dá)時(shí)，植物具有耐寒性，是一種冷調(diào)節(jié)基因。脫水蛋白基因通常在植物脫水時(shí)表達(dá)，脫水可能是由干旱、滲透脅迫或低溫引起的［33］。有研究表明，COR47（AT1G20440）是在低溫下積累的主要脫氫酶（DHNs），有助于冷應(yīng)激反應(yīng)，在體外已實(shí)現(xiàn)COR47的離子和水結(jié)合、低溫保護(hù)活性、類囊體膜結(jié)合和金屬結(jié)合等功能，且COR47的過表達(dá)與低溫條件下擬南芥抗冷性的提高有關(guān)［34］。

以上結(jié)果表明，利用WGCNA分析可挖掘到與目標(biāo)性狀的生物學(xué)意義高度關(guān)聯(lián)的基因模塊和核心基因，為挖掘目標(biāo)基因提供新的研究思路，為解析復(fù)雜的農(nóng)藝性狀提供參考。本研究重點(diǎn)關(guān)注了Blue和Brown這2個(gè)相關(guān)性高的調(diào)控模塊，其余的基因模塊雖然沒有被詳細(xì)論述，但也可能包含與子葉抗冷相關(guān)的通路，可進(jìn)一步挖掘其蘊(yùn)含的生物學(xué)意義。