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

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

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

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

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

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

1 材料與方法

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

谷子品種‘豫谷1 號’（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 號’12 d 的幼苗期開始，分別在6 ℃冷脅迫處理0.0、0.5、1.0、3.0、6.0、16.0、24.0 h，在每個時間點分別從對照組和冷脅迫處理組同時取樣，進行RNA 測序[16]。xiaomi突變體的11 個不同組織、時期的轉(zhuǎn)錄組數(shù)據(jù)來自于本實驗室[4]，數(shù)據(jù)已上傳NGDC（national genomics data center）的GSA（genome sequence archive）數(shù) 據(jù) 庫（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 下載的數(shù)據(jù)首先將SRA 格式文件利用SRAtoolkit 軟件的fast-dump 命令轉(zhuǎn)換為fastq序列文件[17]。所有的測序原始數(shù)據(jù)利用FastQC對其進行質(zhì)控[18]，然后利用Trimmomatic 軟件（軟件中參數(shù)threads 設為10，phred 設為33）去除接頭，過濾掉低質(zhì)量的reads，得到clean data[19]。利用Hisat2（參數(shù)phred 設為12）將clean data 比對到本團隊構建的谷子超早熟突變體xiaomi高質(zhì)量基因組[20]。利用featurecounts[21]對數(shù)據(jù)的reads計數(shù)，之后根據(jù)基因表達的TPM（transcripts per million）值，編寫計算TPM 值的R 語言代碼，構建基因表達矩陣[22]

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

1.2 共表達網(wǎng)絡構建

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

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

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

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

1.4 模塊的GO富集分析

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

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

將核心基因在水稻（https://rapdb.dna.affrc.go.jp）及玉米數(shù)據(jù)庫（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 驗證，檢測候選基因?qū)涿{迫的響應情況。用NCBI 的Primer-BLAST（https://www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi?LINK_LOC=BlastHome）進行引物設計，如表1 所示。內(nèi)參基因為谷子Actin基因，并采用2-ΔΔCt方法[28]計算基因的相對表達量。

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

2 結(jié)果與分析

2.1 谷子加權基因共表達網(wǎng)絡構建

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

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

2.2 冷脅迫相關模塊提取

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

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

對6 種冷脅迫處理的差異表達基因進行分析（圖2），發(fā)現(xiàn)有9 個基因在6 種冷脅迫時間下均差異表達。通過同源比對，這9 個基因在水稻和擬南芥（Arabidopsis thaliana）中的同源基因主要涉及TPR、AP2、WRKY、bHLH、TIFY、PAO2 和MYB家族以及赤霉素、冷調(diào)控等相關功能（表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 種冷脅迫處理條件下的表達模式，結(jié)果（圖3）顯示，與各處理的對照組相比，除Si6g03220在各處理下均呈上調(diào)表達外，其他8個基因在0.5和1.0 h冷處理下呈下調(diào)表達，而在3.0、6.0、16.0 和24.0 h 冷處理下呈持續(xù)上調(diào)表達。其中Si3g01950與Si3g02940以及Si2g38970與Si9g50340隨著處理時間延長，上調(diào)表達的強度呈遞增趨勢（圖3A）。同時，對這9 個基因在xiaomi全生育期的表達模式進行分析，發(fā)現(xiàn)，Si6g03220、Si3g01950、Si3g02940和Si9g50340在全生育時期的表達量均較高（圖3B）。

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

2.3 模塊的GO富集分析

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

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

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

圖5 3個模塊中與冷脅迫相關基因的共表達網(wǎng)絡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 目標模塊中冷脅迫候選基因的篩選及功能分析

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

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

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

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

3 討論

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

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

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

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

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

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