基于WGCNA的谷子苗期冷脅迫應答基因網絡構建與核心因子發掘

張會， 王越越， 趙波， 張麗玲， 郄倩茹， 韓淵懷， 李旭凱

（1.山西農業大學生命科學學院，山西 太谷 030801； 2.山西農業大學山西省后稷實驗室，太原 030031；3.山西農業大學農學院，山西 太谷 030801）

谷子（Setaria italica）古稱稷、粟、粱，是起源于我國黃河流域的古老作物，一萬多年前由青狗尾草（Setaria viridis）馴化而來[1]。谷子是二倍體（2n=2X=18），基因組較?。s440 Mb）[2-4]，具有耐旱、耐貧瘠、適應性強等特點，被譽為“五谷之首”，與C3 模式作物水稻（Oryza sativa）“優勢互補”“水旱對照”“南北呼應”，成為旱生C4 禾谷類作物分子育種研究的模式植物[5]。本實驗室通過對‘晉谷21’進行EMS（ethylmethane sulfonate）誘變獲得矮稈超早熟突變體材料xiaomi，xiaomi生育期短，約60 d，其全基因組測序的完成為本研究奠定了基礎[4]。

植物生存缺乏自主移動性，易遭受低溫、干旱、鹽堿等非生物脅迫，嚴重影響植物的生長發育[6]。低溫脅迫對作物整個生育過程都會造成不同程度的影響，如種子萌發、植株生長發育、產量和品質形成等[7]。近年來，國內外多個研究團隊在水稻耐冷相關的數量性狀位點（quantitative trait locus，QTL）克隆與調控機制解析等方面取得了一定進展[8]。水稻低溫感受器基因COLD1編碼膜定位蛋白，可與G 蛋白α 亞基RGA1 互作以感知低溫，激活Ca2+通道，并增強G 蛋白GTP 酶活性；COLD1能夠正向調控下游轉錄因子OsDREB1B和OsDREB1C的表達，并偶聯葉綠體中的維生素E-維生素K1 代謝調控網絡，進而調控水稻的耐冷性[9-10]。Ca2+作為第二信使，在冷信號的激活和轉導中發揮重要作用。研究發現，過表達膜定位的鈣通道蛋白OsCNGC9 可以增強水稻幼苗的耐冷性：在轉錄水平，OsCNGC9受轉錄因子OsDREB1A的直接轉錄激活調控；在蛋白水平，激酶OsSAPK8 可與OsCNGC9 互作并使其磷酸化，進而激活Ca2+的內流[11]；AP2 轉錄因子家族的CBF/DREBs 轉錄因子在單、雙子葉植物中高度保守，可以促進多個冷相關基因（COR）的表達，增強植株的耐冷性[12]。此外，OsWKRY71/94/45/76、OsMYB3R-2/30、OsTCP-5/6/8/21、OsbZIP38/52/71/73、OsSPL3/14/17、OsMADS57和OsSLR1等轉錄因子基因均直接參與水稻的低溫脅迫調控[13]。玉米中也有較多耐冷相關QTL 定位和候選基因預測的報道，但克隆的基因較少[14]。谷子中耐冷基因定位及克隆的研究較為匱乏。因此，挖掘谷子中的耐冷基因有助于進一步探索谷子的抗逆機理，闡明其分子遺傳機制，有助于推進谷子耐冷遺傳改良進程。

加權基因共表達網絡分析（weighted gene coexpression network analysis, WGCNA）是在鑒定出表達模式相似的基因集合（module）的基礎上，解析不同基因集合與不同分組樣品表型之間的關聯性，繪制基因集合之間的調控網絡并鑒定關鍵調控基因的方法[15]。該方法多被用于研究共表達網絡與植物性狀之間的生物學關系，能夠輔助挖掘與目標性狀高度關聯的核心基因。

本研究利用谷子的根、葉、穗、莖等不同組織、不同冷處理條件下的33 份轉錄組數據，利用WGCNA 方法構建谷子共表達網絡模塊，并進行差異表達分析，將冷脅迫處理組與對照組進行比較，設置錯誤發現率（false discovery rate，FDR）閥值為0.01 篩選差異，與WGCNA 模塊相結合，在模塊中找出谷子抗冷基因，預測與谷子耐冷相關的候選基因，為進一步研究谷子的抗逆機理奠定基礎。

1 材料與方法

1.1 谷子RNA-Seq數據的獲取與分析

谷子品種‘豫谷1 號’（Yugu 1）冷脅迫處理的部分轉錄組數據中來源于NCBI（national center for biotechnology information）的SRA（sequence read archive）數據庫（https://www.ncbi.nlm.nih.gov/sra/），從‘豫谷1 號’12 d 的幼苗期開始，分別在6 ℃冷脅迫處理0.0、0.5、1.0、3.0、6.0、16.0、24.0 h，在每個時間點分別從對照組和冷脅迫處理組同時取樣，進行RNA 測序[16]。xiaomi突變體的11 個不同組織、時期的轉錄組數據來自于本實驗室[4]，數據已上傳NGDC（national genomics data center）的GSA（genome sequence archive）數 據 庫（https://bigd.big.ac.cn/）；其他部分為本實驗室RNA-seq 測序所得。

將‘晉谷21’（JG21）、‘牛毛白’（NMB）和‘平定大谷’（PD）種子播種于山西省太谷縣（37°25′13″ N，112°35′26″ E）大田。取‘晉谷21’出苗2周整株幼苗的混合葉、種植30 d的倒2葉及灌漿期S2、S4 時期未成熟的種子；對于‘牛毛白’及‘平定大谷’均只取灌漿期S2、S4 時期未成熟的種子。將所有樣品取樣后立即置于液氮中，隨后采用RNAprep 純植物試劑盒（天根生化科技有限公司）分離總RNA，并送樣測序，每個樣本3 個生物學重復。

對于NCBI 下載的數據首先將SRA 格式文件利用SRAtoolkit 軟件的fast-dump 命令轉換為fastq序列文件[17]。所有的測序原始數據利用FastQC對其進行質控[18]，然后利用Trimmomatic 軟件（軟件中參數threads 設為10，phred 設為33）去除接頭，過濾掉低質量的reads，得到clean data[19]。利用Hisat2（參數phred 設為12）將clean data 比對到本團隊構建的谷子超早熟突變體xiaomi高質量基因組[20]。利用featurecounts[21]對數據的reads計數，之后根據基因表達的TPM（transcripts per million）值，編寫計算TPM 值的R 語言代碼，構建基因表達矩陣[22]

式中，xi和xj分別表示比對到基因i上的reads數和比對到任一基因j上的reads 數；Li表示基因i的外顯子長度的總和；Lj表示任一基因j的外顯子長度總和。

1.2 共表達網絡構建

基因的表達量矩陣來自于谷子不同組織不同處理下的基因表達量。參照Langfelder[23]的方法，使用R 軟件（R version 3.4.4）中的WGCNA 包（R version 1.6.6）構建加權基因共表達網絡。利用WGCNA 包中的pickSoftThreshold 函數計算權重值，根據谷子轉錄組數據不同軟閾值下的網絡拓撲結構分析結果，選擇軟閾值power=12，使用默認參數利用函數blockwiseModules 構建無尺度網絡。

1.3 谷子耐冷基因檢索及差異表達基因分析

通過國家水稻數據中心（http://www.ricedata.cn/ontology/）檢索已報道的水稻耐冷相關基因，通過blastp 將其蛋白序列與谷子xiaomi的蛋白庫進行比對[24]，獲取與水稻耐冷相關基因的蛋白序列相似性最高的谷子基因，統計這些基因分布最多的3個模塊。

將count 數據文件在百邁克云平臺進行基因差異表達分析，采用FDR 來篩選差異，設置FDR閥值為0.01 且|log2FC|>2 進行差異表達基因分析，用在線軟件（http://jvenn.toulouse.inra.fr/app/example.html）制作韋恩圖[25]，并統計差異表達基因在上述3個模塊中的分布數目。

1.4 模塊的GO富集分析

提取出含冷脅迫基因最多的3 個模塊，利用TBtools 軟件對3 個模塊中的所有基因進行GO 富集分析[26]；之后選出模塊內與冷GO注釋相關的基因，利用Cytoscape（version 3.6）對模塊中的這些基因進行可視化[27]。

1.5 冷脅迫候選基因功能注釋

將核心基因在水稻（https://rapdb.dna.affrc.go.jp）及玉米數據庫（https://www.maizegdb.org）中進行同源比對，設置閾值E-value<10-5且Identity>85%進行篩選，參照水稻和玉米中的同源基因功能進行注釋。

1.6 冷脅迫應答候選基因的實時熒光定量PCR（RT-qPCR）分析

為了驗證冷脅迫應答候選基因的表達模式，對其中6個候選基因進行了RT-qPCR分析。將‘預谷1 號’播種12 d 后對其進行6 ℃冷脅迫處理，分別在處理0、3、6、24 h時，從對照組和冷脅迫處理組取樣（地上部），每組樣本取3個生物學重復，使用博邁德Green qPCR Master Mix 熒光定量試劑盒進行RT-qPCR 驗證，檢測候選基因對冷脅迫的響應情況。用NCBI 的Primer-BLAST（https://www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi?LINK_LOC=BlastHome）進行引物設計，如表1 所示。內參基因為谷子Actin基因，并采用2-ΔΔCt方法[28]計算基因的相對表達量。

表1 RT-qPCR引物序列Table 1 RT-qPCR primer sequence

2 結果與分析

2.1 谷子加權基因共表達網絡構建

利用WGCNA 包中的goodSamplesGenes 函數檢測并過濾掉表達量偏低的基因，最終獲得32 017 個高表達基因。利用WGCNA 包內的pickSoftThreshold 函數計算并選擇權重值β=12（R2=0.8），以實現所有基因構成的關系網絡符合無尺度網絡分布。按權重值β=12，計算并獲得每2 個基因間的拓撲重疊矩陣（topological overlap matrix，TOM）；之后采用動態剪切算法進行共表達網絡模塊劃分，最終得到44個共表達模塊，如圖1所示，圖中每種顏色代表聚成一簇的基因所形成的模塊。其中darkturquoise 模塊包含的基因最多，有7 415 個；navajowhite 模塊包含的基因最少，僅37個。

圖1 基因聚類樹和樣品分割Fig. 1 Gene cluster dendrogram and module detecting

2.2 冷脅迫相關模塊提取

利用水稻已報道的與冷脅迫直接相關的基因在谷子數據庫中進行檢索比對，共鑒定到68 個冷脅迫候選基因。通過映射查找，這些基因分布于18 個共表達網絡模塊，其中，darkturquoise、turquoise、blue 模塊中所含的耐冷基因數量較多，分別為13、10、8個。

為了評估冷脅迫相關基因映射選取模塊的準確性，對谷子冷脅迫處理0.5、1.0、3.0、6.0、16.0和24.0 h 的轉錄組數據進行差異表達基因分析。結果顯示，turquoise 模塊中不同處理下的差異基因數量最多，占所有差異基因數的12.69%；darkturquoise 模塊中的差異基因數次之，占12.06%；blue模塊中差異基因數占7.68%。

對6 種冷脅迫處理的差異表達基因進行分析（圖2），發現有9 個基因在6 種冷脅迫時間下均差異表達。通過同源比對，這9 個基因在水稻和擬南芥（Arabidopsis thaliana）中的同源基因主要涉及TPR、AP2、WRKY、bHLH、TIFY、PAO2 和MYB家族以及赤霉素、冷調控等相關功能（表2），這些家族均與植物非生物脅迫應答相關[29-31]。

圖2 谷子6個冷脅迫時間差異表達基因分析Fig. 2 Analysis of differentially expressed genes in foxtail millet under cold stress

表2 谷子6個冷脅迫處理條件下同時差異表達基因Table 2 Differentially expressed genes under 6 cold stress conditions in foxtail millet

進一步比較這9 個基因在6 種冷脅迫處理條件下的表達模式，結果（圖3）顯示，與各處理的對照組相比，除Si6g03220在各處理下均呈上調表達外，其他8個基因在0.5和1.0 h冷處理下呈下調表達，而在3.0、6.0、16.0 和24.0 h 冷處理下呈持續上調表達。其中Si3g01950與Si3g02940以及Si2g38970與Si9g50340隨著處理時間延長，上調表達的強度呈遞增趨勢（圖3A）。同時，對這9 個基因在xiaomi全生育期的表達模式進行分析，發現，Si6g03220、Si3g01950、Si3g02940和Si9g50340在全生育時期的表達量均較高（圖3B）。

圖3 差異表達基因的表達模式Fig. 3 Expression patterns of differentially expressed genes

2.3 模塊的GO富集分析

選取冷響應差異表達基因數量最多的模塊，并結合谷子中與水稻已報道耐冷基因同源最多的模塊，進行GO 富集分析。對冷脅迫基因數最多的3個模塊中的所有基因進行GO 功能富集分析，富集結果涉及到生物學過程（biological process，BP）、分子功能（molecular function，MF）和細胞組 分（cellular component，CC）。由 圖4 可 知，darkturquoise 模塊富集到了對寒冷的響應（GO:0009409）、低溫馴化（GO:0009631）等相關的調控通路；turquoise 模塊富集到了半胱氨酸類型的肽酶活性（GO:0008234）、作用于蛋白質的催化活性（GO:0140096）及脫落酸激活的信號通路（GO:0009738）等相關的調控通路；blue模塊富集到了多種激酶活性的調控通路，如氧化還原酶活性（GO:0016491）、碳水化合物激酶活性（GO:0019200）。

圖4 模塊特有的基因GO富集結果Fig. 4 GO enrichment results of module specific gene

從上述GO 富集結果中挑選與冷脅迫相關的基因，進行權重共表達網絡分析，并利用Cytoscape 軟件對網絡進行可視化，如圖5 所示。同時對模塊darkturquoise、turquoise 及blue 內權重較高的核心基因分別在水稻和擬南芥數據庫中進行同源基因注釋（部分見表3）。結果表明，這些基因涉及ERF、bHLH 和WRKY 等家族，而這些家族基因均與植物非生物脅迫應答相關[29-32]，可作為谷子耐冷脅迫研究的候選基因。

圖5 3個模塊中與冷脅迫相關基因的共表達網絡Fig. 5 Co-expression network of genes relative to cold stress in 3 modules

表3 3個模塊中與冷脅迫相關的核心基因的功能注釋Table 3 Functional annotation of hub genes relative to cold stress in 3 modules

2.4 目標模塊中冷脅迫候選基因的篩選及功能分析

篩選與核心基因高連通性的基因作為候選基因，為了進一步了解這些候選基因的功能，通過功能注釋及文獻報道發現候選基因都與冷脅迫相關。Si8g19810和Si7g12250在水稻和玉米中的同源基因均屬于bHLH 家族；Si4g26860和Si7g26280在水稻和玉米中的同源基因屬于bZIP家族；Si1g19560和Si9g44570屬于ERF 家族；這幾個家族的轉錄因子都與植物非生物脅迫應答相關[33-39]。因此，通過WGCNA 構建冷脅迫基因共表達網絡，并進一步篩選具有冷脅迫相關生物學意義的候選基因是可行的。從目標模塊中的篩選獲得的候選基因可作為谷子冷脅迫相關基因挖掘的可用資源，可作為后續研究中重點關注的對象。

2.5 候選基因冷脅迫應答的RT-qPCR分析

為進一步驗證候選基因對冷脅迫的響應模式，對其中6 個核心脅迫應答基因進行RT-qPCR驗證，結果（圖6）顯示，與對照相比， 6個基因在不同冷脅迫處理時間后的表達量均呈現出上調趨勢，其表達模式與轉錄組分析結果（圖3）基本一致，證明利用轉錄組數據分析得到的冷脅迫響應差異表達基因是可靠的。

圖6 6個核心基因在冷脅迫不同時期的表達Fig. 6 Expression of 6 hub genes in different times under cold stress

3 討論

3.1 共表達網絡分析富集到多個冷脅迫相關的基因

本研究通過對谷子冷脅迫處理組及對照組多個組織、不同時間的33 份RNA-seq 數據進行分析，構建了共表達網絡，根據聚類情況將其劃分為44 個模塊。利用水稻數據庫中已報道與冷脅迫相關的基因，在谷子蛋白數據庫中比對尋找谷子中的同源基因，并與劃分的模塊進行比較，發現在18 個模塊中都能找出與冷脅迫相關的基因。對冷脅迫相關基因數量最多的3 個模塊進行了GO功能富集分析，各個模塊都得到了與冷脅迫應答相關的功能富集結果，并且發現不同模塊內的響應存在差異。其中，darkturquoise 模塊富集到了對寒冷的響應（GO:0009409）、低溫馴化（GO:0009631）等相關的調控通路；turquoise模塊富集到了半胱氨酸類型的肽酶活性（GO:0008234）、作用于蛋白質的催化活性（GO:0140096）及脫落酸激活的信號通路（GO:0009738）等相關的調控通路；blue模塊富集到了多種激酶活性的調控通路，如氧化還原酶活性（GO:0016491）、碳水化合物激酶活性（GO:0019200）。

對turquoise模塊共表達網絡中的核心基因進行分析發現，核心基因Si4g02480注釋到了對寒冷的響應（GO:0009409），且該基因在水稻蛋白數據庫中的同源基因為脫水應答元件結合蛋白，屬于DREB轉錄因子的ERF亞族基因。此亞族基因通過與啟動子區域的順式作用元件特異性結合來調控與干旱、高鹽和低溫等非生物脅迫相關的應答基因表達[33]，因此在植物抵抗非生物脅迫過程中起重要作用[34]；此外，ERF 基因家族可以促進其下游抗逆相關基因的表達，提高植物的抗逆境脅迫能力[35]。本研究對各模塊的候選基因進行注釋發現，Si9g05520在水稻和玉米中的同源基因LOC_Os03g58890和GRMZM2G054224都 屬 于 氧化還原酶，其酶活性受溫度影響；另有2 個候選基因屬于bHLH 轉錄因子；有4 個基因屬于ERF 轉錄因子。研究表明，bHLH 轉錄因子家族參與植物的生長發育，其表達受鹽、冷等非生物脅迫誘導[36]，還可以響應低溫，進而調控植物生長發育[37]；AP2/ERF 家族受低溫誘導，同時激活或抑制下游基因的表達[38]。本研究除檢測到ERF轉錄因子外，還有MYB、bZIP、WRKY、bHLH、TCP 轉錄因子。據報道，MYB 轉錄因子對高溫具有耐受性[39]；bZIP在植物應對鹽害、干旱、冷害、機械損傷以及滲透脅迫的非生物脅迫中發揮重要作用[40]；WRKY 轉錄因子可結合目標基因啟動子上的作用元件W-BOX（TGACCA/T），從而調控下游目的基因的表達，調節植物適應低溫脅迫，其編碼蛋白質通常表現為抑制作用，在低溫脅迫下，下游靶基因的表達量下調[41]。由此表明，多個轉錄因子家族與冷脅迫相關，因此，本研究中發現的候選基因為后續谷子耐冷育種提供了基因資源。

3.2 Si6g03220 可能是谷子冷脅迫響應的關鍵候選基因

利用‘豫谷1 號’在6 個不同冷處理時間下的差異表達基因集，并結合全生育期的表達模式進行分析，發現Si6g03220在全生育期都有較高表達，且在不同冷脅迫處理下均上調表達，該基因被注釋為乙烯反應相關轉錄因子11，屬AP2 基因家族。研究表明，當植物受到冷害和凍害刺激時，體內的乙烯含量顯著增加[42-43]；而乙烯在植物冷脅迫應答中具有一定的調控作用[44]。因此，該基因可作為研究谷子耐冷機制的重要候選基因。

本研究通過對谷子多時期、多組織的冷脅迫數據進行共表達網絡分析，得到了44 個模塊，其中有3 個模塊的冷脅迫響應基因數較多，進一步對他們進行功能富集分析，篩選出了與冷脅迫相關的候選基因；通過對‘豫谷1號’在6個不同冷脅迫處理中的差異表達基因進行分析，獲得了9 個高置信度的候選基因；轉錄組數據和RT-qPCR 驗證表明，這些基因與冷脅迫應答密切相關。以上研究結果為谷子冷脅迫響應機理研究奠定了基礎，得到的候選基因為谷子抗逆育種研究提供重要資源。