李旭凱 李任建 張寶俊,*

1 山西農業(yè)大學生命科學學院,山西太谷030801；2 山西農業(yè)大學農學院,山西太谷 030801

干旱、冷害和高鹽等非生物逆境會對植物的生長和發(fā)育產生不利影響,造成植物生長發(fā)育緩慢,嚴重時可導致植物死亡[1]。干旱、冷害和高鹽導致的植物應激可破壞其細胞結構和影響關鍵的生理功能[1],使植物細胞產生滲透脅迫,水分外流,蛋白失活或變性,產生大量的活性氧(reactive oxygen species,ROS)進而損傷細胞,最終導致植物光合作用被抑制,代謝功能紊亂[2]。在農業(yè)方面,這些逆境反應造成了作物產量和品質的降低。植物的耐受性會隨著環(huán)境的改變而改變[3],同時受許多基因控制[4],基因工程為發(fā)展多耐性作物提供了方法,可將多個基因組裝導入作物,實現(xiàn)作物更高的耐受性[5]。高通量測序技術的飛速發(fā)展,也為基因的挖掘提供了快捷的方法。研究人員利用轉錄組學和代謝組學分析了煙草對低溫脅迫的反應[6]。王濤研究團隊通過建立黃花苜蓿非生物脅迫轉錄組數據庫,揭示了植物膜錨定轉錄因子在干旱脅迫下重定位進入細胞核調控機體氧化還原平衡的分子機制[7]。

加權共表達網絡分析(Weighted Gene Co-expression Network Analysis,WGCNA)利用植物體內各個生命活動互相關聯(lián)的特點,通過高通量測序技術獲得基因的表達情況,對基因間的連接系統(tǒng)進行冪次處理,使網絡內部的基因符合無尺度網絡拓撲結構,即網絡中存在的少數基因與大多數基因存在聯(lián)系的特點[8],該中心基因即為起關鍵調控的基因(hub基因)。該方法將具有相似表達模式的基因聚類到同一個模塊,多被用于研究共表達網絡與植物性狀之間的生物學關系,挖掘與性狀之間高度關聯(lián)的關鍵基因。本研究利用水稻的根、葉、花蕊、種子等材料,以及在干旱、冷、鹽脅迫下的不同處理時間共47 份轉錄組數據,構建了共表達網絡,得到了15 個模塊。從國家水稻數據中心(http：//www.ricedata.cn/)獲取已克隆的干旱、冷、鹽相關的406 個基因(附表1),通過查看已知基因所在的模塊,并通過對冷、干旱、鹽脅迫處理的不同時間段的轉錄組數據的差異表達分析,根據已報道相關脅迫基因在模塊中的數量和各脅迫不同處理中差異表達基因存在的模塊數量,確定并對數目較多的模塊進行生物學功能分析。通過對干旱、冷、鹽脅迫共同存在的模塊進行了功能分析,預測水稻耐受性關聯(lián)的關鍵基因,為提高水稻耐受性方面的研究提供了有效參考,打開了該研究方向的新思路。

1 材料與方法

1.1 水稻轉錄組數據的獲取與分析

水稻干旱脅迫、冷脅迫、鹽脅迫處理的100 個轉錄組數據來源于NCBI的SRA(Sequence Read Archive)數據庫(附表2)。首先利用SRAtoolkit[9]軟件的fast- dump 將SRA 格式文件轉換為序列文件fastq,利用FastQC[10]軟件對原始測序數據進行質量評估,緊接著用Trimmonatic[11]軟件去除接頭、過濾低質量的reads,得到clean data。利用Hisat2[12]軟件將clean data 比對到水稻基因組。水稻基因組序列和注釋信息下載自Rice Genome Annotation Project(http：//rice.plantbiology.msu.edu)數據庫。用featureCounts[13]獲得數據的reads計數,之后根據基因表達的TPM 值(transcripts per million)公式編寫計算TPM 值的R語言代碼,構建基因表達矩陣,R語言代碼如下：

1.2 共表達網絡構建和已知基因模塊的選取

基因的表達量矩陣來自不同脅迫處理下的基因表達量,帶重復的樣本進行均值預處理。使用R 軟件(R version 3.4.4)和WGCNA(R version 1.6.6)包,構建加權基因共表達網絡與劃分相關模塊[14]。使用WGCNA 包中的pickSoftThreashold 計算權重值,本實驗選取power 值為12,使用blockwiseModules 構建無尺度網絡,參數按照默認設置[14]。利用edgeR[15]分別對數據庫中下載到的冷脅迫處理3 h、6 h、24 h的水稻幼苗和干旱脅迫處理1 h、6 h、24 h 的水稻幼苗以及鹽脅迫處理1 h、5 h、24 h 的水稻幼苗進行差異表達分析,分析差異表達所在的模塊,并結合水稻已克隆報道的脅迫基因,確定已報道基因所在模塊及數量,確定各脅迫下的特異性模塊。

1.3 模塊的GO 富集分析

首先提取模塊的中基因,利用TBtools(A Toolkit for Biologists integrating various biological data handling tools with a user-friendly interface)進行GO富集分析[16]。使用水稻基因組作為參考數據庫,P值小于0.005。利用Cytoscape(version 3.6.1)對模塊中的網絡進行可視化[17]。

1.4 轉錄因子注釋分析

通過PlantTFDB 數據庫(http：//planttfdb.cbi.pku.edu.cn/index.php?sp=Osj),獲得水稻中的轉錄因子共2409 個,分屬于52 個家族[18]。與模塊中預測的基因進行比較,得到共表達模塊中預測的轉錄因子信息(附表3)。

2 結果與分析

2.1 水稻加權基因共表達網絡構建

對水稻表達矩陣中表達量較低的基因進行過濾,獲得30,339 個高表達基因。選擇權重值β=12 來構建無尺度網絡,之后按照混合動態(tài)剪切的標準劃分模塊,隨之依次計算每個模塊的特征向量值(eigengenes),對距離較近的模塊進行合并,最終獲得15個共表達模塊(圖1),每個模塊用不同的顏色來表示聚類到一起的基因。其中pink 模塊聚類到的基因數目最多,為7386 個,bisque4 模塊包含的基因最少,有286 個。

將水稻數據庫中已報道與干旱、冷、鹽脅迫相關的基因與模塊中的基因比較。已報道的干旱基因在共表達網絡中有156 個,其中l(wèi)ightgreen 模塊和green 模塊存在較多,分別為32 個和29 個；已報道的冷脅迫相關基因共有164 個,其中l(wèi)ightgreen 和green 模塊存在基因數目最多,分別為31 個和32 個；鹽脅迫的已報道基因在共表達模塊中共有70 個,在green 模塊和darkmagenta 模塊中分別存在18 個和17 個(圖2)。干旱、冷、鹽脅迫相關基因被聚在各個模塊中,其中在green 模塊中大量的聚集。

為了驗證和確定模塊選取的準確性,同時對數據中冷脅迫處理3 h、6 h、12 h 處理的水稻幼苗,干旱處理1 h、6 h、24 h 的水稻幼苗和鹽脅迫處理1 h、5 h、24 h 的水稻幼苗3 份轉錄組數據進行了差異表達基因分析。在對冷脅迫的差異表達分析中,模塊darkmagenta 中不同處理下的差異基因數目占所有基因總數的 17%左右,green 模塊中占18%左右。在干旱脅迫差異表達分析中,模塊green 中的差異基因占總差異基因數目的13%左右,lightgreen 模塊差異基因約占總數目的16%左右。在對鹽脅迫的差異表達分析中發(fā)現(xiàn),green 模塊中的差異基因約占總差異基因數目的 13%左右,lightgreen 模塊中5 h 處理下的差異基因數目約占40%左右(附表4)。

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

圖2 報道基因分布情況 Fig.2 Distribution of opented genes

2.2 模塊的GO 富集分析

選取已報道基因所在最多的模塊,結合不同脅迫下差異表達基因在模塊中的分布結果(附表4),干旱選取lightgreen 模塊和green 模塊,冷脅迫選取green 模塊和 darkmagenta 模塊,鹽脅迫選取lightgreen 模塊和green 模塊。對這些已報道基因有關聯(lián)的基因進行GO 功能富集分析,富集結果涉及到了生物學過程(biological process,BP)、分子功能(molecular function,MF)和細胞組分(cellular component,CC)。

與干旱相關的lightgreen 和green 模塊都富集到了與干旱相關的調控通路,主要是與脂肪酸分解過程(GO：0009062)、活性氧的響應(GO：0000302)、脫落酸響應(GO：0000302)、滲透脅迫響應(GO：00009737)等相關。與冷相關的green 模塊和darkmagenta 模塊富集到了與鈣依賴性蛋白激酶活性(GO：0010857)、細胞酰胺代謝過程(GO：0044106)等,以及茉莉酸介導的信號轉導通路的調控(GO：2000022)。與鹽脅迫相關的green 模塊和lightgreen 模塊也富集到了鹽脅迫 響 應(GO：0009651)、內源性激素的響應(GO：0009719)和激素應答(GO：0009725)的一些響應脅迫相關的功能(表1和附表5)。

表1 模塊GO 富集情況(部分)Table1 GO enrichments of network modules(part)

(續(xù)表1)

植物對逆境的響應是多個基因調控的,在目前已報道與逆境脅迫相關的基因中,仍有部分潛在的基因其在逆境脅迫中發(fā)揮的具體作用機制尚不明確。本研究對干旱脅迫、冷脅迫、鹽脅迫中已報道逆境相關基因存在較多的6 個模塊,通過權重共表達網絡分析,并使用Cytoscape 軟件對網絡進行可視化(圖3)。其中,在darkmagenta 中的WLP1[19](LOC_ Os01g54540)、OsPDS[20](LOC_Os03g08570)、OsCYL4[21](LOC_Os09g02270)、hbd2[22](LOC_Os02g40860)；lightgreen 模塊中的OsDREB1G[23](LOC_Os08g 43210)、OCPI2[24](LOC_Os01g42860)、OsMSR15[25](LOC_Os03g41390)、SRWD1[26](LOC_ Os08g38880)、OsNHX1[27](LOC_Os07g47100)都處于網絡的核心位置(圖3)。對處于網絡中心權重值較高的基因在水稻中進行了注釋和在擬南芥中進行了同源基因的注釋(表2)。

2.3 干旱、冷、鹽脅迫相關基因功能富集分析

對不同脅迫在共表達網絡分析中得到的基因進行韋恩圖分析,發(fā)現(xiàn)干旱、冷、鹽脅迫的相關報道基因大量被聚類到green 模塊(圖3)。冷脅迫共有3006 個基因被富集到了green 模塊,其中與干旱共有的基因數目為2680 個,與冷脅迫相關的基因數目為2716 個,特有的基因有209 個；干旱脅迫在綠色模塊富集到的基因數目多達4226 個,與鹽脅迫的共有基因比與冷脅迫共有的基因數目要多,有3663 個,特有的基因數目為462 個；鹽脅迫總共有4057 個基因被富集到綠色模塊；三者共有的基因數目為2599個(圖4)。

(圖3)

圖3 冷、干旱、鹽脅迫相關的共表達網絡 Fig.3 Coexpression network related to cold,drought,and salt stress

表2 模塊核心基因功能注釋 Table2 Functional annotation of modular hub genes

(續(xù)表2)

(續(xù)表2)

對冷脅迫富集到了焦磷酸酶活性(GO：0016462)、核苷三磷酸酶活性(nucleoside-triphosphatase activity)等分子功能(molecular function,MF),與樣本具有較高生物學關系的響應冷脅迫(GO：0009409)等生物學過程(biological process,BP),以及有機生物合成過程(GO：1901576)、細胞糖類代謝過程(GO：0044262)等生物學過程(圖5)。

對干旱脅迫特有的基因,富集到一些與蛋白質相關的功能,包括蛋白質結合(GO：0005515)、作用于蛋白質的催化活性(GO：0140096)等；在細胞組分中,顯著富集到4 個與細胞膜相關的功能；在生物學過程中,富集到蛋白質修飾(GO：0006464)、蛋白質磷酸化(GO：0006468)等功能,也富集到與信號轉導相關的細胞表明受體信號通路(GO：0007166),脫落酸刺激的細胞反應(GO：0071215),和一些與免疫相關的功能,包括免疫系統(tǒng)過程(GO：0002376)、免疫應答(GO：0006955)、防御反應(GO：0006952)、固有免疫反應(GO：0045087)等(圖5)。

對鹽脅迫特有的基因,最顯著富集到的是信號受體活性(GO：0038023)、分子傳感性活性(GO：0060089)；細胞組分中富集到一些與膜組分相關的功能,包括膜的整體組份(GO：0016021)、膜部分(GO：0044425)；生物學功能富集到較多與鹽脅迫相關的功能,包括滲透脅迫反應(GO：0006970)、鹽脅迫 反應(GO：0009651)、對水的反應(GO：0009415)、對水楊酸的反應(GO：0009751)、對乙烯的反應(GO：0009723)(圖5)。

圖4 綠色模塊中與3 種脅迫相關的共表達基因 Fig.4 Modular network gene in co-expression network

在3 種脅迫共有的基因中有20 個注釋到了已報道與冷脅迫相關的基因,35 個注釋到了已報道與干旱脅迫相關的基因,35 個注釋到了已報道與鹽脅迫 相關的基因,其中已報道的同時參與3 種脅迫的基因有3 個,分別是編碼AT-hook 蛋白的OsAHL1[28]、編碼絲裂原活化蛋白激酶的OsMAPK5[29]和NAC 轉錄因子OsNAC6[30]；同時參與調控干旱和鹽脅迫的基因有17 個,同時參與調控干旱和冷脅迫的基因有6 個,同時參與調控冷和鹽脅迫的基因有2 個(表3)。

圖5 綠色模塊中3 種脅迫特有共表達基因的GO 功能富集 Fig.5 GO functional enrichment of three kinds of stress-specific co-expressed genes in the green module

表3 多種脅迫調控基因功能注釋 Table3 Annotation of functions of various stress regulating genes

(續(xù)表3)

(續(xù)表3)

通過GO 功能富集分析,干旱、冷、鹽脅迫在分子功能中富集到了與氧化反應相關的氧化還原酶活性(GO：0016491)、過氧化物酶活性(GO：0004601)、氧化還原酶活性對過氧化物作為受體的作用(GO：0016684)、抗氧化活性(GO：0016209)等。在受到干旱、冷、鹽脅迫后,酰胺的生物合成過程(GO：0043604)、細胞酰胺代謝(GO：0043603),過氧化氫分解代謝過程(GO：0042744)、氧化還原過程(GO：0055114)、過氧化氫代謝過程(GO：0072593)、毒素代謝過程(GO：0009404)；苯丙烷代謝過程(GO：0009698)和木質素代謝過程(GO：00009808)也得到了富集,也與一些細胞程序性死亡負調節(jié)(GO：0043069)和細胞分解代謝過程(GO：0044248)顯著富集(表4和附表6)。

表4 模塊GO 富集情況(部分)Table4 GO enrichments of network modules(part)

以已報道的對 3 種脅迫都具有調控作用的OsAHL1、DSM1和ONAC095為核心基因,構建共表達網絡,利用Cytoscape 對網絡進行可視化。其中LOC_Os02g30900 編碼的蛋白激酶結構域含蛋白質是NAC 轉錄因子(LOC_Os01g66120)和絲裂原活化蛋白激酶(LOC_Os03g17700)的樞紐基因(圖6)。對這些預測到的樞紐基因在玉米數據庫中進行檢索,得到同源基因,并進行了功能注釋(表5)。

圖6 Green 模塊中3 種脅迫下共有基因的共表達網絡 Fig.6 Co-expression network of shared genes under three stresses in green module

表5 候選基因功能注釋 Table5 Annotation of candidate gene

(續(xù)表5)

3 討論

本研究通過對干旱、冷、鹽脅迫不同處理下47個不同材料建立共表達網絡,根據聚類情況劃分為15 個模塊,對已克隆報道的水稻數據庫中與干旱、冷、鹽脅迫相關基因進行統(tǒng)計,并與劃分的模塊進行比較,發(fā)現(xiàn)各個模塊都能得到與各脅迫相關的基因。在對每個脅迫下存在的基因數目最多的2 個模塊進行了GO 功能富集分析,各個模塊都得到了具有相關生物學意義的功能富集結果,并且發(fā)現(xiàn)不同的脅迫下寄主的響應也有所不同。如與干旱脅迫相關的 lightgreen 模塊中的滲透脅迫響應(GO：0009737)、對水的響應(GO：0006970),與冷脅迫相關的茉莉酸介導的信號轉導通路的調控(GO：2000022),darkmagenta 模塊中的對溫度刺激響應(GO：0009266)通路等；鹽脅迫相關的green 模塊中的鹽脅迫響應(GO：0009651),這些富集到的功能與樣品的處理方式生物功能基本吻合,證明了WGCNA 方法在對大數據分析時具有較高的可靠性。在對模塊green 中3種脅迫下各自特有的基因進行富集也表現(xiàn)出功能上的明顯不同,其中對冷脅迫的響應較敏感。對共表達網絡中的核心基因分析發(fā)現(xiàn),在darkmagenta 中,OsPDS(LOC_Os03g08570)編碼的八氫番茄紅素脫氫酶,在水稻中報道為對干旱敏感性增加,而耐冷性和氧化脅迫性增強[31]。OsCOIN(LOC_ Os01g01420)編碼的鋅指蛋白,過量表達可增強水稻對冷、鹽和干旱的抗性,對冷脅迫較為敏感[32]。在lightgreen 中,處于表達網絡中心的OsMSR15(LOC_ Os03g41390)是一類具有轉錄激活活性鋅指蛋白,在冷、干旱、鹽脅迫中均能強烈響應,是植物響應干旱的重要調節(jié)因子[33]。SRWD1(LOC_Os08g 38880)編碼鹽應答WD40 蛋白,受鹽脅迫調節(jié),具有多種調節(jié)模式[34]。對預測的連通性較高的基因進行注釋發(fā)現(xiàn),LOC_Os01g01420在水稻中功能尚不明確,擬南芥中注釋為編碼過氧化氫酶,參與多重脅迫反應[35]。對green 模塊分析和共有的基因富集分析中發(fā)現(xiàn),酰胺生物合成、細胞酰胺代謝被顯著富集,而在擬南芥中,IAA 酰胺合成酶GH3-6 負調控干旱和鹽脅迫的反應[36]。苯丙烷代謝中的木質素生物合成在非生物脅迫下被誘導,通過積累保護寄主細胞[37]。在構建的共表達網絡中,處于樞紐位置的基因被注釋到了DUF26 激酶,該激酶在擬南芥中被報道參與防御和程序化細胞的調節(jié)[38]。

本研究對干旱、冷、鹽這3 種脅迫下共同存在的調控機制進行了分析,但鑒定到的結果都是直接或間接地參與脅迫調控,需進一步挖掘和驗證這些基因的具體生物功能?？傊?利用WGCNA 共表達網絡對干旱、冷、鹽脅迫相關功能基因進行挖掘,對水稻抗逆的研究具有重要的意義。

4 結論

通過對冷、干旱、鹽脅迫數據的共表達網絡分析,得到了15 個模塊。對這3 種脅迫下各自的生物學意義進行了描述,對參與2 種或3 種脅迫的基因進行了鑒定,構建了具有生物學意義的共表達網絡。本研究結果可為后續(xù)的水稻多抗逆性和綜合抗逆性相關的調控機制研究提供重要的參考依據。

附表請見網絡版：1)本刊網站http：//zwxb.chinacrops.org/；2)中國知網http：//www.cnki.net/；3)萬方數 據 http：//c.wanfangdata.com.cn/Periodical-zuowxb.aspx。