利用WGCNA鉴定甘薯耐盐相关共表达网络及核心基因

张 毅，吴万亿，刘霞宇，张 洁，唐锐敏，贾小云

(1.山西农业大学 生命科学学院，山西 太谷 030801； 2.山西农业大学 农学院，山西 太谷 030801)

随着全球人口的快速增加，人类对粮食的需求量日益提高，但全球土地资源有限，并且约有20%的土地正遭受盐碱的危害。土地盐碱化是一个全球性问题，造成大片耕地不能正常使用，据估计，到2050年，将有50%的耕地盐碱化[1]。在盐碱化土地上种植的植物会遭受离子毒害、渗透胁迫和氧化胁迫等，进而抑制植物生长甚至产量[2]。甘薯[Ipomoeabatatas(L.)Lam]是重要的粮食作物，其块根富含淀粉、膳食纤维、维生素和花青素等，具有抗癌和保健作用[3-4]。甘薯对环境有较好的适应快，但盐碱胁迫仍然是影响甘薯正常生长的主要因素。因此，挖掘甘薯盐胁迫响应基因将为解析甘薯适应盐环境的机制和培育甘薯耐盐新品种提供重要的理论依据。

随着高通量测序技术的发展和测序价格的降低，海量的转录组数据产生[5]。传统的2个样本间转录组数据的差异分析已经满足不了需求，而且传统的转录组数据差异分析只能对有功能注释的基因进行分析。加权基因共表达网络分析(Weighted gene co-expression network analysis，WGCNA)可以利用多个转录组数据构建基因共表达网络，通过软阈值来鉴定表达模式相似的基因集合模块(Module)，并特异性筛选出与性状相关的共表达模块，进而挖掘出共表达模块中的关键基因。表达模式相似的基因通常被认为可能参与相似的生物学过程，WGCNA是研究表达模式相似基因集的重要工具[6]。WGCNA通过构建无尺度分布的网络，可以对参与生物学过程的没有功能注释的新基因进行预测，并结合已有的基因功能注释发现参与特定代谢途径的共表达模块、预测新基因的功能。WGCNA在植物学研究中广泛应用。前人利用WGCNA分析水稻不同时间点的镉胁迫转录组数据，鉴定到了271个差异表达基因参与镉的调控[7]；利用WGCNA分析了红皮甘薯及其突变体，鉴定出2个调节甘薯皮色的基因模块[8]；利用WGCNA分析了红皮苹果及其突变体黄皮苹果的转录组数据，发现1个与花青素含量高度相关的共表达模块，同时发现MdMYB10基因的启动子甲基化是导致果皮颜色差异的原因[9]。但是目前尚未见关于利用WGCNA鉴定甘薯耐盐相关共表达网络及核心基因的报道。为此，利用甘薯盐胁迫转录组数据结合WGCNA方法，鉴定甘薯耐盐相关共表达网络及核心基因，为进一步研究甘薯耐盐分子机制及耐盐甘薯新品种的培育奠定基础。

1 材料和方法

1.1 数据获取与分析

甘薯盐胁迫转录组数据(SRP092215)是从EBI(European Bioinformatics Institute)的ENA(European Nucleotide Archive)数据库获得的[10]。该转录组数据是以栗子香和ND98两个甘薯品种为材料，对生长30 d的甘薯幼苗利用 200 mmol/L NaCl溶液进行盐胁迫处理，分别在处理后0、12、48 h取根系进行转录组测序，每个处理2个重复，共12个样本[11]。首先用Aspera软件下载甘薯盐胁迫转录组(SRP092215)数据的fastq文件，然后利用FastQC[12]软件对原始数据进行质量评估处理，利用fastp软件[13]去除接头和低质量碱基，利用Hisat2软件[14]将过滤后的clean data与甘薯参考基因组[15](https://ipomoea-genome.org)比对，然后利用featureCount软件[16]计算每个基因reads数，最后使用TPM(Transcripts per kilobase of exon model per million mapped reads)体现基因表达量。

1.2 加权基因共表达网络的构建

利用R 3.61软件中WGCNA 1.69软件包[17]构建加权基因共表达网络，以标准化后的12个盐胁迫转录组数据中的基因表达矩阵作为输入数据，使用R软件中genefiliter 软件包筛选在各个样本间表达量变异最大的前50%基因构建加权共表达网络，利用WGCNA 1.69软件中的pickSoftThreshold 函数计算共表达网络的软阈值，并利用powerEstimate函数估计最优power值，最终选择power=16对原始有尺度关系矩阵进行幂处理，得到无尺度化邻接矩阵。为更好评估基因间表达模式的相关性, 进一步将邻接矩阵转化为拓扑重叠矩阵(Topological overlap matrix,TOM)，并利用拓扑相异矩阵(dissTOM=1-TOM)，采用动态剪切算法进行基因聚类及模块划分。使用WGCNA自动网络构建函数 blockwiseModules 构建网络,其中，最小模块基因数为30，相似模块的合并阈值为0.25(mergeCutHeight=0.25)，其他参数按照默认设置。

1.3 盐胁迫特异性模块的鉴定

对每个共表达模块中的所有基因进行主成分分析(Principle component analysis,PCA)，将主成分1(PC1)称为此模块的特征向量(Module eigengene,ME)，为了筛选盐胁迫相关特异性模块，计算各个模块的ME值与不同性状(此处为不同盐胁迫时间)的相关系数r和P值。r>0表示正相关，r<0表示负相关。在本研究中，若模块ME值和性状相关系数|r|>0.70且P<0.05，则认为该模块是特异性模块。

1.4 GO功能富集分析

利用TBtools (A tool kit for biologists integrating various biological data handling tools with a user-friendly interface)软件[18]对每个共表达模块中的基因进行GO(Gene ontology)富集分析，阈值为P<0.05。

1.5 耐盐特异性模块核心基因的鉴定及基因互作网络构建

将共表达模块中每个基因的表达量与模块ME值进行关联，得到模块关系值(Module membership，MM)，其可以体现基因在模块中的重要性，MM值越高，基因在模块中越重要。选取耐盐正相关特异性模块中MM值排名前5的基因为模块中的核心基因，并利用Cytoscape软件展示核心基因互作网络。

1.6 转录因子分析

将各个模块中的蛋白质序列提交到plantTFDB数据库[19](http://planttfdb.gao-lab.org/)进行转录因子分析预测，获得各个模块中预测的转录因子。

2 结果与分析

2.1 甘薯加权共表达网络的构建

本研究过滤掉变异较小的基因，最终选择32 147个基因用于构建加权基因共表达网络。使用WGCNA软件包中powerEstimate函数估计最适power值为16，当power值为16时，无尺度网络拟合指数R2>0.85，平均连通性趋于0(图1)。该power值可以获得符合要求的无尺度网络。采用动态剪切法划分模块，合并相似度大于75%的模块，最终构建了20个共表达模块，不同颜色表示不同模块(图2)。各模块含基因数目如图3，其中，Turquoise模块所含基因数最多，为10 383个；Royalblue模块中基因数最少，为124;其他模块的基因数介于242～6 745；Grey模块中的基因是一组没有分配到其他模块中的基因集。

a.不同软阈值对应的无尺度网络拟合指数R2，直线代表R2=0.85; b.不同软阈值对应的平均连通性

上部为基因聚类树，下部是按树的分枝切割得到的模块，相同的模块用同一种颜色表示

图3 模块中基因数目分布

2.2 甘薯盐胁迫相关特异性模块的鉴定

由图4可知，本研究共鉴定到6个与盐胁迫特异相关的基因共表达模块(|r|>0.70，P<0.05)。其中，Green模块(r=0.71，P=0.01)和Red模块(r=0.85，P=5e-04)与盐胁迫48 h正相关；Black模块(r=0.83，P=9e-04)和Yellow模块(r=0.79，P=0.002)与盐胁迫12 h正相关；Midnightblue模块(r=-0.74,P=0.006)和Magenta模块(r=-0.73,P=0.007)与盐胁迫48 h负相关。为了验证模块划分的可靠性，进一步选择与盐胁迫正相关的4个特异性模块(Green模块、Black模块、Yellow模块和Red模块)作为目标模块做深入分析。

括号外数据、括号内数据分别表示模块和性状之间的相关系数、P值

2.3 甘薯盐胁迫正相关特异性模块GO富集分析

由图5可知，选取4个与甘薯盐胁迫正相关的特异性模块(Green模块、Black模块、Yellow模块和Red模块)进行GO富集分析，发现其均可以富集到生物学过程、分子功能和细胞组分。Green模块富集到139个生物学过程、26个细胞组分、41个分子功能，主要富集到干旱响应(GO:0009414)、叶绿体(GO:0009507)、非生物刺激的反应(GO:0009628)、高渗响应(GO:0006972)等；Black模块富集到84个生物学过程、51个细胞组分、24个分子功能，主要富集到跨膜运输相关通路，如活性离子跨膜转运活性(GO:0022853)、质子跨膜运输(GO:1902600)、跨膜运输(GO:0055085 )等；Yellow模块富集到265个生物学过程、24个细胞组分、53个分子功能，主要富集到信号传递(GO:0023052)、质膜(GO:0005886)、脱落酸激活的信号通路(GO:0009738)、离子结合(GO:0043167)等；Red模块富集到63个生物学过程、25个细胞组分、22个分子功能，主要富集到抗氧化活性(GO:0010035)、光合作用膜(GO:0034357)、水杨酸刺激的细胞反应(GO:0071446)等。

图5 甘薯耐盐相关特异性模块GO富集分析

2.4 甘薯盐胁迫显著正相关共表达模块中核心基因的鉴定及基因互作网络构建

本研究选取与甘薯盐胁迫显著正相关的共表达模块中MM 值排名前5的核心基因，并将其与拟南芥蛋白质序列进行比对，预测核心基因的功能(表1)。核心基因主要有HAK5(High affinity K+transporter 5)、TLP5(Tubby like protein 5)、NAGS2(N-acetyl-L-glutamate synthase 2)、bHLH115(Basic helix-loop-helix 115)、SK11(Shaggy-related kinase 11)、ABCG28(ABC transporter G family member 28)、DUF699(GNAT acetyltransferase)等，其中部分核心基因功能已经明确，如Green模块中核心基因bHLH115是一个参与植物缺铁响应的正调控转录因子基因，与参与响应植物逆境胁迫的基因共表达;Yellow模块中的核心基因SK11可以提高植物的耐盐性；Red模块中的核心基因ABCG28有物质转运功能。进一步探寻与核心基因互作关系强的基因，利用Cytoscape软件展示互作网络(图6)，发现与核心基因互作的基因也有相应的功能并且已经被报道，如Green模块中的核心基因TLP5(g42340)与NAC[NAM(no apical meristem)，ATAF1(Arabidopsistranscription activation factor 1)，ATAF2，CUC2(cup-shaped cotyledon 2)]基因家族NAC2(g28394)基因互作，核心基因g48306与HsfB2a(Heat shock transcription factor B2a)(g20568)互作；Black模块中的核心基因g14417与FAR1(Far-red-impaired response 1)(g54849)基因存在互作；Yellow模块中核心基因g23125与DREB2C(Dehydration responsive element binding protein 2C)(g1860)基因互作;Red模块中核心基因g11018与ERF1(Ethylene response factor 1)(g46237)基因互作。

表1 甘薯盐胁迫显著正相关共表达模块中的核心基因及其功能

灰色圆形代表核心基因，三角形代表转录因子，灰色三角形代表既是核心基因又是转录因子，圆形代表其他基因

3 结论与讨论

植物在长期适应盐胁迫的过程中形成了应对盐胁迫的分子机制。本研究利用甘薯盐胁迫的转录组数据，构建了加权基因共表达网络，挖掘出与盐胁迫显著相关的6个特异性模块，并对其中的4个与盐胁迫显著正相关的基因模块进行GO富集分析，发现其可以显著富集到与耐盐相关的具有生物学意义的条目，主要包括跨膜运输、水杨酸刺激的细胞反应、非生物刺激的反应及高渗响应等条目。其中，Green模块可以显著富集到非生物刺激的反应(GO:0009628)和干旱响应(GO:0009414)等GO条目。研究表明，植物在盐胁迫过程中由于渗透胁迫引起细胞水势过低，植物无法吸收足够的水分满足其正常生命活动，进而引起植物生理干旱[20]。Black模块可以显著富集到活性离子跨膜转运活性(GO:0022853)和质子跨膜运输(GO:1902600)等相关通路。研究表明，植物在响应盐胁迫的过程中进行渗透调节，植物生物膜上的质子泵通过水解ATP(5′-adenylate triphosphate)产生能量，定向运输氢离子形成跨膜电势梯度，从而可以使得盐离子进行跨膜运输[21]。植物的质子泵与离子通道等决定了其对离子的选择性吸收。植物在盐胁迫下，质子泵提供能量，把钠离子排出细胞外，使得植物体内钠离子含量减少，进而使细胞处于稳态和代谢正常[22]。Red模块可以显著富集到水杨酸刺激的细胞反应(GO:0071446)、抗氧化活性(GO:0010035)等条目。Yellow模块可以显著富集到脱落酸激活的信号通路(GO:0009738)。研究表明，脱落酸、茉莉酸、水杨酸等植物激素在植物响应盐胁迫中起着重要的作用，其中脱落酸最关键，其可以参与调控耐盐相关基因的表达、气孔关闭和离子平衡等[23]。

本研究选择模块中MM值排名前5的基因作为核心基因进行研究，模块中核心基因为HAK5、TLP5、NAGS2、bHLH115、TBL19、SK11、DUF699等。由于核心基因在网络中具有较高的连通性，所以其可能参与重要的生物学过程[24]。通过对目标模块中的核心基因进行功能分析，结果表明，其主要功能为跨膜运输、蛋白质磷酸化、转录调控等，然而发现大部分核心基因的功能在甘薯中鲜有报道，但是它们在拟南芥、番茄等植物中的同源基因已经被报道过可响应盐胁迫。如Green模块中的核心基因g34447在拟南芥中的同源基因bHLH115在维持体内铁稳态中起着关键作用[25],g1032在番茄中的同源基因NAGS1可提高拟南芥的耐盐性和抗旱性[26]；Black模块中的核心基因g57672在拟南芥中的同源基因HAK5，属于高亲和K+转运体HAK，可以调节K+离子浓度，进而提高植物的耐盐性[27]；Yellow模块中的核心基因g23125在拟南芥中的同源基因SK11在盐胁迫下被激活，随后又激活G6PD基因(Glucose-6-phosphate dehydrogenase)，提高了耐盐性[28]；Red模块中的核心基因g33775在拟南芥中的同源基因SRO5与P5CDH(Pyrroline-5-carboxylate dehydrogenase)共同控制活性氧(ROS)产生和响应胁迫反应[29]。植物耐盐性不是由单个基因决定的，而是多个抗逆基因共同作用的结果。与核心基因互作的基因也可能参与重要的生物学过程。本研究通过对与核心基因互作的基因进行分析，发现部分基因功能已经被报道。如与Green模块中核心基因g42340互作的g28394基因在拟南芥中的同源基因NAC2响应多种非生物胁迫和灰霉病菌侵染[30]。与Green模块中核心基因g48306互作的g20568基因在辣椒中的同源基因HsfB2a参与辣椒对青枯病菌的免疫和对高温高湿的耐受性[31]。与Black模块中核心基因g14417互作的g54849基因在拟南芥中的同源基因FAR1通过整合叶绿素生物合成和水杨酸信号通路，在调节植物免疫中发挥重要作用[32]。与Yellow模块中核心基因g23125互作的g1860基因在拟南芥中的同源基因DREB2C在氧化应激耐受性方面发挥重要作用[33]。与Yellow模块中核心基因g56260互作的g58671基因在番茄中的同源基因ZFP1(Zinc-finger protein 1)可以提高番茄幼苗对非生物胁迫的耐受能力[34]。与Red模块中核心基因g11018互作的g46237基因在拟南芥中的同源基因ERF1通过整合茉莉酸、乙烯和脱落酸信号调控胁迫相关基因，在高盐、干旱和高温胁迫中发挥积极作用[35]。

本研究侧重分析了与盐胁迫显著正相关的模块，其他模块中也可能有一些重要基因参与甘薯的耐盐胁迫。本研究所分析的核心基因在甘薯中的功能尚未明确，后续需要通过分子生物学试验进一步验证。