利用WGCNA鉴定非生物胁迫相关基因共表达网络

李旭凯 李任建 张宝俊,*

1 山西农业大学生命科学学院,山西太谷030801；2 山西农业大学农学院,山西太谷 030801

干旱、冷害和高盐等非生物逆境会对植物的生长和发育产生不利影响,造成植物生长发育缓慢,严重时可导致植物死亡[1]。干旱、冷害和高盐导致的植物应激可破坏其细胞结构和影响关键的生理功能[1],使植物细胞产生渗透胁迫,水分外流,蛋白失活或变性,产生大量的活性氧(reactive oxygen species,ROS)进而损伤细胞,最终导致植物光合作用被抑制,代谢功能紊乱[2]。在农业方面,这些逆境反应造成了作物产量和品质的降低。植物的耐受性会随着环境的改变而改变[3],同时受许多基因控制[4],基因工程为发展多耐性作物提供了方法,可将多个基因组装导入作物,实现作物更高的耐受性[5]。高通量测序技术的飞速发展,也为基因的挖掘提供了快捷的方法。研究人员利用转录组学和代谢组学分析了烟草对低温胁迫的反应[6]。王涛研究团队通过建立黄花苜蓿非生物胁迫转录组数据库,揭示了植物膜锚定转录因子在干旱胁迫下重定位进入细胞核调控机体氧化还原平衡的分子机制[7]。

加权共表达网络分析(Weighted Gene Co-expression Network Analysis,WGCNA)利用植物体内各个生命活动互相关联的特点,通过高通量测序技术获得基因的表达情况,对基因间的连接系统进行幂次处理,使网络内部的基因符合无尺度网络拓扑结构,即网络中存在的少数基因与大多数基因存在联系的特点[8],该中心基因即为起关键调控的基因(hub基因)。该方法将具有相似表达模式的基因聚类到同一个模块,多被用于研究共表达网络与植物性状之间的生物学关系,挖掘与性状之间高度关联的关键基因。本研究利用水稻的根、叶、花蕊、种子等材料,以及在干旱、冷、盐胁迫下的不同处理时间共47 份转录组数据,构建了共表达网络,得到了15 个模块。从国家水稻数据中心(http：//www.ricedata.cn/)获取已克隆的干旱、冷、盐相关的406 个基因(附表1),通过查看已知基因所在的模块,并通过对冷、干旱、盐胁迫处理的不同时间段的转录组数据的差异表达分析,根据已报道相关胁迫基因在模块中的数量和各胁迫不同处理中差异表达基因存在的模块数量,确定并对数目较多的模块进行生物学功能分析。通过对干旱、冷、盐胁迫共同存在的模块进行了功能分析,预测水稻耐受性关联的关键基因,为提高水稻耐受性方面的研究提供了有效参考,打开了该研究方向的新思路。

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 来构建无尺度网络,之后按照混合动态剪切的标准划分模块,随之依次计算每个模块的特征向量值(eigengenes),对距离较近的模块进行合并,最终获得15个共表达模块(图1),每个模块用不同的颜色来表示聚类到一起的基因。其中pink 模块聚类到的基因数目最多,为7386 个,bisque4 模块包含的基因最少,有286 个。

将水稻数据库中已报道与干旱、冷、盐胁迫相关的基因与模块中的基因比较。已报道的干旱基因在共表达网络中有156 个,其中lightgreen 模块和green 模块存在较多,分别为32 个和29 个；已报道的冷胁迫相关基因共有164 个,其中lightgreen 和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%左右。在对盐胁迫的差异表达分析中发现,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 模块。对这些已报道基因有关联的基因进行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)

(续表1)

植物对逆境的响应是多个基因调控的,在目前已报道与逆境胁迫相关的基因中,仍有部分潜在的基因其在逆境胁迫中发挥的具体作用机制尚不明确。本研究对干旱胁迫、冷胁迫、盐胁迫中已报道逆境相关基因存在较多的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 干旱、冷、盐胁迫相关基因功能富集分析

对不同胁迫在共表达网络分析中得到的基因进行韦恩图分析,发现干旱、冷、盐胁迫的相关报道基因大量被聚类到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

(续表2)

(续表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),和一些与免疫相关的功能,包括免疫系统过程(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

(续表3)

(续表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)也得到了富集,也与一些细胞程序性死亡负调节(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

(续表5)

3 讨论

本研究通过对干旱、冷、盐胁迫不同处理下47个不同材料建立共表达网络,根据聚类情况划分为15 个模块,对已克隆报道的水稻数据库中与干旱、冷、盐胁迫相关基因进行统计,并与划分的模块进行比较,发现各个模块都能得到与各胁迫相关的基因。在对每个胁迫下存在的基因数目最多的2 个模块进行了GO 功能富集分析,各个模块都得到了具有相关生物学意义的功能富集结果,并且发现不同的胁迫下寄主的响应也有所不同。如与干旱胁迫相关的 lightgreen 模块中的渗透胁迫响应(GO：0009737)、对水的响应(GO：0006970),与冷胁迫相关的茉莉酸介导的信号转导通路的调控(GO：2000022),darkmagenta 模块中的对温度刺激响应(GO：0009266)通路等；盐胁迫相关的green 模块中的盐胁迫响应(GO：0009651),这些富集到的功能与样品的处理方式生物功能基本吻合,证明了WGCNA 方法在对大数据分析时具有较高的可靠性。在对模块green 中3种胁迫下各自特有的基因进行富集也表现出功能上的明显不同,其中对冷胁迫的响应较敏感。对共表达网络中的核心基因分析发现,在darkmagenta 中,OsPDS(LOC_Os03g08570)编码的八氢番茄红素脱氢酶,在水稻中报道为对干旱敏感性增加,而耐冷性和氧化胁迫性增强[31]。OsCOIN(LOC_ Os01g01420)编码的锌指蛋白,过量表达可增强水稻对冷、盐和干旱的抗性,对冷胁迫较为敏感[32]。在lightgreen 中,处于表达网络中心的OsMSR15(LOC_ Os03g41390)是一类具有转录激活活性锌指蛋白,在冷、干旱、盐胁迫中均能强烈响应,是植物响应干旱的重要调节因子[33]。SRWD1(LOC_Os08g 38880)编码盐应答WD40 蛋白,受盐胁迫调节,具有多种调节模式[34]。对预测的连通性较高的基因进行注释发现,LOC_Os01g01420在水稻中功能尚不明确,拟南芥中注释为编码过氧化氢酶,参与多重胁迫反应[35]。对green 模块分析和共有的基因富集分析中发现,酰胺生物合成、细胞酰胺代谢被显著富集,而在拟南芥中,IAA 酰胺合成酶GH3-6 负调控干旱和盐胁迫的反应[36]。苯丙烷代谢中的木质素生物合成在非生物胁迫下被诱导,通过积累保护寄主细胞[37]。在构建的共表达网络中,处于枢纽位置的基因被注释到了DUF26 激酶,该激酶在拟南芥中被报道参与防御和程序化细胞的调节[38]。

本研究对干旱、冷、盐这3 种胁迫下共同存在的调控机制进行了分析,但鉴定到的结果都是直接或间接地参与胁迫调控,需进一步挖掘和验证这些基因的具体生物功能。总之,利用WGCNA 共表达网络对干旱、冷、盐胁迫相关功能基因进行挖掘,对水稻抗逆的研究具有重要的意义。

4 结论

通过对冷、干旱、盐胁迫数据的共表达网络分析,得到了15 个模块。对这3 种胁迫下各自的生物学意义进行了描述,对参与2 种或3 种胁迫的基因进行了鉴定,构建了具有生物学意义的共表达网络。本研究结果可为后续的水稻多抗逆性和综合抗逆性相关的调控机制研究提供重要的参考依据。

附表请见网络版：1)本刊网站http：//zwxb.chinacrops.org/；2)中国知网http：//www.cnki.net/；3)万方数 据 http：//c.wanfangdata.com.cn/Periodical-zuowxb.aspx。