基于WGCNA的棉花子葉抗冷相關共表達模塊鑒定

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

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

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

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

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

1 材料與方法

1.1 試驗材料及處理

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

1.2 權重基因共表達網絡構建

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

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

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

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

1.4 核心基因的qRT-PCR驗證

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

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

2 結果與分析

2.1 構建權重基因共表達網絡

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

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

2.1.2 基因共表達網絡中基因的聚類 根據基因間表達量的相關性構建聚類樹，采用動態切割法將產生的聚類樹切割，把表達模式相似的基因合并在同一分支上，每個分支代表1個共表達模塊，根據模塊相似度（0.8）對表達模式相似的模塊合并后進行模塊劃分，最終獲得9個共表達模塊（圖2），每一種顏色代表一個模塊，灰色（Grey）模塊代表無法歸入任何一個模塊的基因。藍色（Blue）模塊包含的基因數量最多，為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 共表達模塊中基因數量分布Fig.3 Distribution of number of genes in co-expression modules

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

圖4 基因共表達網絡模塊與不同樣本的關聯熱圖Fig.4 Association analysis of gene co-expression network modules with different samples

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

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

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

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

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

圖7 Brown模塊的基因共表達網絡及其核心基因Fig.7 Gene co-expression network and hub genes in Brown module

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

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

2.3 核心基因的RT-qPCR驗證

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

圖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個冷脅迫處理下與棉花品種和處理時間相關聯的共表達模塊，利用各個模塊的特征向量基因對模塊進行相關性分析，發現Brown和Blue模塊間存在極顯著的相關性且關聯度最高，可作為子葉抗冷性狀機理研究的目標模塊。對這2個模塊內的基因進行GO功能富集分析，均富集到了抗逆相關的具有生物學意義的調控途徑。例如，Blue和Brown模塊富集到轉錄因子活性（GO：0003700）、對脅迫的反應（GO：0006950）等通路，該結果與謝勇軍［11］的研究結果一致。

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

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