基于WGCNA的棉花子葉抗冷相關(guān)共表達模塊鑒定

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

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

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

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

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

1 材料與方法

1.1 試驗材料及處理

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

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

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

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

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

1.4 核心基因的qRT-PCR驗證

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

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

2 結(jié)果與分析

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

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

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

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

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

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

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

圖4 基因共表達網(wǎng)絡模塊與不同樣本的關(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)性達到-0.86（圖5）。因此將這兩個模塊作為抗冷相關(guān)特異性模塊進行深入研究，挖掘模塊中的核心基因。

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

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

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

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

圖7 Brown模塊的基因共表達網(wǎng)絡及其核心基因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驗證

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

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

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

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

3 討論

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

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

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