基于WGCNA的谷子苗期冷脅迫應(yīng)答基因網(wǎng)絡(luò)構(gòu)建與核心因子發(fā)掘

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

（1.山西農(nóng)業(yè)大學(xué)生命科學(xué)學(xué)院，山西 太谷 030801； 2.山西農(nóng)業(yè)大學(xué)山西省后稷實(shí)驗(yàn)室，太原 030031；3.山西農(nóng)業(yè)大學(xué)農(nóng)學(xué)院，山西 太谷 030801）

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

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

加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析（weighted gene coexpression network analysis, WGCNA）是在鑒定出表達(dá)模式相似的基因集合（module）的基礎(chǔ)上，解析不同基因集合與不同分組樣品表型之間的關(guān)聯(lián)性，繪制基因集合之間的調(diào)控網(wǎng)絡(luò)并鑒定關(guān)鍵調(diào)控基因的方法[15]。該方法多被用于研究共表達(dá)網(wǎng)絡(luò)與植物性狀之間的生物學(xué)關(guān)系，能夠輔助挖掘與目標(biāo)性狀高度關(guān)聯(lián)的核心基因。

本研究利用谷子的根、葉、穗、莖等不同組織、不同冷處理?xiàng)l件下的33 份轉(zhuǎn)錄組數(shù)據(jù)，利用WGCNA 方法構(gòu)建谷子共表達(dá)網(wǎng)絡(luò)模塊，并進(jìn)行差異表達(dá)分析，將冷脅迫處理組與對(duì)照組進(jìn)行比較，設(shè)置錯(cuò)誤發(fā)現(xiàn)率（false discovery rate，F(xiàn)DR）閥值為0.01 篩選差異，與WGCNA 模塊相結(jié)合，在模塊中找出谷子抗冷基因，預(yù)測(cè)與谷子耐冷相關(guān)的候選基因，為進(jìn)一步研究谷子的抗逆機(jī)理奠定基礎(chǔ)。

1 材料與方法

1.1 谷子RNA-Seq數(shù)據(jù)的獲取與分析

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

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

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

式中，xi和xj分別表示比對(duì)到基因i上的reads數(shù)和比對(duì)到任一基因j上的reads 數(shù)；Li表示基因i的外顯子長(zhǎng)度的總和；Lj表示任一基因j的外顯子長(zhǎng)度總和。

1.2 共表達(dá)網(wǎng)絡(luò)構(gòu)建

基因的表達(dá)量矩陣來自于谷子不同組織不同處理下的基因表達(dá)量。參照Langfelder[23]的方法，使用R 軟件（R version 3.4.4）中的WGCNA 包（R version 1.6.6）構(gòu)建加權(quán)基因共表達(dá)網(wǎng)絡(luò)。利用WGCNA 包中的pickSoftThreshold 函數(shù)計(jì)算權(quán)重值，根據(jù)谷子轉(zhuǎn)錄組數(shù)據(jù)不同軟閾值下的網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)分析結(jié)果，選擇軟閾值power=12，使用默認(rèn)參數(shù)利用函數(shù)blockwiseModules 構(gòu)建無尺度網(wǎng)絡(luò)。

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

通過國家水稻數(shù)據(jù)中心（http://www.ricedata.cn/ontology/）檢索已報(bào)道的水稻耐冷相關(guān)基因，通過blastp 將其蛋白序列與谷子xiaomi的蛋白庫進(jìn)行比對(duì)[24]，獲取與水稻耐冷相關(guān)基因的蛋白序列相似性最高的谷子基因，統(tǒng)計(jì)這些基因分布最多的3個(gè)模塊。

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

1.4 模塊的GO富集分析

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

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

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

1.6 冷脅迫應(yīng)答候選基因的實(shí)時(shí)熒光定量PCR（RT-qPCR）分析

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

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

2 結(jié)果與分析

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

利用WGCNA 包中的goodSamplesGenes 函數(shù)檢測(cè)并過濾掉表達(dá)量偏低的基因，最終獲得32 017 個(gè)高表達(dá)基因。利用WGCNA 包內(nèi)的pickSoftThreshold 函數(shù)計(jì)算并選擇權(quán)重值β=12（R2=0.8），以實(shí)現(xiàn)所有基因構(gòu)成的關(guān)系網(wǎng)絡(luò)符合無尺度網(wǎng)絡(luò)分布。按權(quán)重值β=12，計(jì)算并獲得每2 個(gè)基因間的拓?fù)渲丿B矩陣（topological overlap matrix，TOM）；之后采用動(dòng)態(tài)剪切算法進(jìn)行共表達(dá)網(wǎng)絡(luò)模塊劃分，最終得到44個(gè)共表達(dá)模塊，如圖1所示，圖中每種顏色代表聚成一簇的基因所形成的模塊。其中darkturquoise 模塊包含的基因最多，有7 415 個(gè)；navajowhite 模塊包含的基因最少，僅37個(gè)。

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

2.2 冷脅迫相關(guān)模塊提取

利用水稻已報(bào)道的與冷脅迫直接相關(guān)的基因在谷子數(shù)據(jù)庫中進(jìn)行檢索比對(duì)，共鑒定到68 個(gè)冷脅迫候選基因。通過映射查找，這些基因分布于18 個(gè)共表達(dá)網(wǎng)絡(luò)模塊，其中，darkturquoise、turquoise、blue 模塊中所含的耐冷基因數(shù)量較多，分別為13、10、8個(gè)。

為了評(píng)估冷脅迫相關(guān)基因映射選取模塊的準(zhǔn)確性，對(duì)谷子冷脅迫處理0.5、1.0、3.0、6.0、16.0和24.0 h 的轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行差異表達(dá)基因分析。結(jié)果顯示，turquoise 模塊中不同處理下的差異基因數(shù)量最多，占所有差異基因數(shù)的12.69%；darkturquoise 模塊中的差異基因數(shù)次之，占12.06%；blue模塊中差異基因數(shù)占7.68%。

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

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

表2 谷子6個(gè)冷脅迫處理?xiàng)l件下同時(shí)差異表達(dá)基因Table 2 Differentially expressed genes under 6 cold stress conditions in foxtail millet

進(jìn)一步比較這9 個(gè)基因在6 種冷脅迫處理?xiàng)l件下的表達(dá)模式，結(jié)果（圖3）顯示，與各處理的對(duì)照組相比，除Si6g03220在各處理下均呈上調(diào)表達(dá)外，其他8個(gè)基因在0.5和1.0 h冷處理下呈下調(diào)表達(dá)，而在3.0、6.0、16.0 和24.0 h 冷處理下呈持續(xù)上調(diào)表達(dá)。其中Si3g01950與Si3g02940以及Si2g38970與Si9g50340隨著處理時(shí)間延長(zhǎng)，上調(diào)表達(dá)的強(qiáng)度呈遞增趨勢(shì)（圖3A）。同時(shí)，對(duì)這9 個(gè)基因在xiaomi全生育期的表達(dá)模式進(jìn)行分析，發(fā)現(xiàn)，Si6g03220、Si3g01950、Si3g02940和Si9g50340在全生育時(shí)期的表達(dá)量均較高（圖3B）。

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

2.3 模塊的GO富集分析

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

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

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

圖5 3個(gè)模塊中與冷脅迫相關(guān)基因的共表達(dá)網(wǎng)絡(luò)Fig. 5 Co-expression network of genes relative to cold stress in 3 modules

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

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

篩選與核心基因高連通性的基因作為候選基因，為了進(jìn)一步了解這些候選基因的功能，通過功能注釋及文獻(xiàn)報(bào)道發(fā)現(xiàn)候選基因都與冷脅迫相關(guān)。Si8g19810和Si7g12250在水稻和玉米中的同源基因均屬于bHLH 家族；Si4g26860和Si7g26280在水稻和玉米中的同源基因?qū)儆赽ZIP家族；Si1g19560和Si9g44570屬于ERF 家族；這幾個(gè)家族的轉(zhuǎn)錄因子都與植物非生物脅迫應(yīng)答相關(guān)[33-39]。因此，通過WGCNA 構(gòu)建冷脅迫基因共表達(dá)網(wǎng)絡(luò)，并進(jìn)一步篩選具有冷脅迫相關(guān)生物學(xué)意義的候選基因是可行的。從目標(biāo)模塊中的篩選獲得的候選基因可作為谷子冷脅迫相關(guān)基因挖掘的可用資源，可作為后續(xù)研究中重點(diǎn)關(guān)注的對(duì)象。

2.5 候選基因冷脅迫應(yīng)答的RT-qPCR分析

為進(jìn)一步驗(yàn)證候選基因?qū)涿{迫的響應(yīng)模式，對(duì)其中6 個(gè)核心脅迫應(yīng)答基因進(jìn)行RT-qPCR驗(yàn)證，結(jié)果（圖6）顯示，與對(duì)照相比， 6個(gè)基因在不同冷脅迫處理時(shí)間后的表達(dá)量均呈現(xiàn)出上調(diào)趨勢(shì)，其表達(dá)模式與轉(zhuǎn)錄組分析結(jié)果（圖3）基本一致，證明利用轉(zhuǎn)錄組數(shù)據(jù)分析得到的冷脅迫響應(yīng)差異表達(dá)基因是可靠的。

圖6 6個(gè)核心基因在冷脅迫不同時(shí)期的表達(dá)Fig. 6 Expression of 6 hub genes in different times under cold stress

3 討論

3.1 共表達(dá)網(wǎng)絡(luò)分析富集到多個(gè)冷脅迫相關(guān)的基因

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

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

3.2 Si6g03220 可能是谷子冷脅迫響應(yīng)的關(guān)鍵候選基因

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

本研究通過對(duì)谷子多時(shí)期、多組織的冷脅迫數(shù)據(jù)進(jìn)行共表達(dá)網(wǎng)絡(luò)分析，得到了44 個(gè)模塊，其中有3 個(gè)模塊的冷脅迫響應(yīng)基因數(shù)較多，進(jìn)一步對(duì)他們進(jìn)行功能富集分析，篩選出了與冷脅迫相關(guān)的候選基因；通過對(duì)‘豫谷1號(hào)’在6個(gè)不同冷脅迫處理中的差異表達(dá)基因進(jìn)行分析，獲得了9 個(gè)高置信度的候選基因；轉(zhuǎn)錄組數(shù)據(jù)和RT-qPCR 驗(yàn)證表明，這些基因與冷脅迫應(yīng)答密切相關(guān)。以上研究結(jié)果為谷子冷脅迫響應(yīng)機(jī)理研究奠定了基礎(chǔ)，得到的候選基因?yàn)楣茸涌鼓嬗N研究提供重要資源。