利用WGCNA鉴定谷子内源脱落酸响应禾生指梗霉胁迫的共表达基因

常国蓉，李任建，张琦，张育铭，韩渊怀，张宝俊

（山西农业大学植物保护学院，山西太谷 030801）

0 引言

【研究意义】脱落酸（ABA）作为一类“逆境激素”，在响应干热、高盐、低温、重金属及辐射等非生物胁迫过程中起着重要作用。同时，脱落酸可直接或通过与水杨酸（SA）、茉莉酸（JA）和乙烯（ETH）等激素的互作或调控等参与寄主的生物胁迫[1]。脱落酸信号转导依赖其受体蛋白 PYR/PYL/PCAR及负调控因子2C类蛋白磷酸酶（PP2C）和正调控因子SNF1相关的蛋白激酶（SnRK2）组成的双重负调控系统的调控，激活其下游基因的表达[2-3]。据报道，在拟南芥（Arabidopsis thaliana）、水稻（Oryza sativa）等作物中，约 10%的基因受脱落酸的调控[3]。谷子（Setaria italica）是我国北方地区一种重要的杂粮作物，具有较强的抗旱、抗逆能力，且生育期短，经济效益高[4]。但谷子生产上常受到多种生物胁迫及非生物胁迫因素的影响，特别是近年来，随着谷子种植规模扩大，由卵菌纲禾生指梗霉（Sclerospora graminicola）引起的谷子白发病在各地频发，对谷子产量和品质的提高造成严重威胁，分析脱落酸在谷子响应禾生指梗霉侵染中的调控作用，有利于解析谷子白发病的发生特点及规律，并为其防治提供新思路。【前人研究进展】稻瘟病菌（Magnaporthe oryzae）在水稻感染部位能够激活脱落酸信号，抑制水稻中最初入侵的免疫信号，从而负调控水稻对稻瘟病菌的抗病性[5-6]。SONG等[7]研究发现，对番茄叶面外施脱落酸，可有效降低番茄早疫病菌（Alternaria solani）对番茄植株的病害程度。ARF10可诱导甘蓝型油菜对脱落酸的敏感性，提高对黑斑病菌（Alternaria brassicicola）的抗性[8]。PYL和SnRK2是脱落酸信号转导途径中两个关键的正反馈调节的家族基因。当植物体内脱落酸大量积累时，脱落酸与其受体PYR/PYLs/RCARs结合后，其以单体形式与磷酸酶PP2Cs结合，通过抑制PP2Cs，从而释放SnRK2s蛋白激酶，SnRK2s磷酸化激活ABFs及RAV1等转录因子、激活下游脱落酸响应基因[9-14]。反之当脱落酸含量较低时，PYR/PYLs/RCARs以二聚体的形式存在，SnRK2的活性受到抑制，阻止脱落酸信号传导[15]。ZHANG等[16]研究发现，SnRK2家族的OsSAPK9可与OsSGT1相互作用，提高水稻的抗白叶枯病能力；YAO等[17]通过GWAS群体和转录关联的方法，鉴定到PYL和SnRK2在玉米早期抗穗腐病中显著上调。【本研究切入点】现有研究表明，脱落酸在寄主响应病菌侵染、扩展、致病过程中发挥重要作用[9]，且PYL和SnRK2是脱落酸信号转导途径中两类关键因子，在脱落酸代谢过程中发挥着重要作用。然而目前尚未有脱落酸代谢在谷子抗病中的相关研究。【拟解决的关键问题】通过基因家族、加权基因共表达网络（WGCNA）、转录组分析对PYL和SnRK2家族基因进行鉴定和分析，探索谷子响应禾生指梗霉侵染过程中脱落酸发挥的作用及关键的调控因子，为后续谷子抗白发病或抗其他生物胁迫的调控机制研究提供参考。

1 材料与方法

1.1 试验材料

感病品种晋谷 21号种子由山西农业大学生物工程研究所种质资源库提供，病菌孢子于2016年在山西农业大学试验田的病株上分离获得。2017年5月初将种子拌菌处理后种植于山西农业大学谷子杂粮育种基地，以消毒不接种谷子为对照组，采取随机区组设计方法划分试验小区，每个小区面积为3 m×2 m。在谷子叶片刚表现出“灰背”症状后，每隔7 d取材一次，样品经液氮速冻后于-80℃保存，用于脱落酸含量测定和RNA提取，共取材5次，以健康叶片为对照，3次重复。

1.2 谷子PYL和SnRK2家族全基因筛选

PYL家族筛选：在 Phytozome数据库（https://phytozome.jgi.doe.gov/pz/portal.html）中获取谷子（v2.2）、水稻（v7_JGI）和拟南芥（TAIR10）的基因组数据。在Pfam数据库（http://pfam.xfam.org/）中下载PYL（Polyketide_cyc2）HMMER模型（PF10604），通过HMMER软件[18]中的hmmsearch命令搜索谷子全基因蛋白质数据，选取比对值 E-value＜0.01。SnRK2家族筛选：利用蛋白激酶模型（PF00069）通过hmmsearch对谷子基因组进行搜索，随后利用拟南芥数据库和水稻数据库下载的SnRK2家族基因对上一步具有蛋白激酶结构的基因进行同源比对，选取比对值E-value＜0.01，利用MEGA计算筛选的基因与水稻、拟南芥SnRK2的遗传距离，筛选谷子SnRK2家族基因。

使用 Pfam中的 Batch search对检索到PYL和SnRK2蛋白序列结构域进行验证。使用TBtools软件[19]对谷子的基因组注释gff3文件的基因定位结果可视化分析。

1.3 谷子PYL和SnRK2蛋白理化性质和基因定位分析

采用在线工具 ExPASY-ProtParam（https://web.expasy.org/protparam/）对谷子PYL和SnRK2家族基因的氨基酸数目、分子量、等电点和总平均疏水指数进行分析[20]。

1.4 谷子PYL和SnRK2的系统进化分析

采用MEGA 7.0[21]中的ClustalW将拟南芥（9个）、水稻（8个）的PYL序列与谷子PYL序列进行多重比对，将拟南芥（5个）、水稻（9个）SnRK2序列与谷子SnRK2序列进行多重比对，构建系统进化树，矫正后，通过Test Neighbor-Joining Tree（邻近法）构建进化树。利用进化树在线美化网站ITOL进行进化树美化。

1.5 谷子PYL和SnRK2的顺式作用元件分析

采用TBtools软件[19]提取谷子PYL和SnRK2家族基因的启动子序列，将提取的启动子序列到在线数据库 PlantCARE（http://bioinformatics.psb.gent.e/eboolslantcare/html）[22]分析谷子PYL和SnRK2家族基因顺式作用元件。

1.6 脱落酸含量的测定

取0.5 g叶片材料在液氮中研磨成粉末，将粉末与3 mL 80%甲醇在4℃下混合4 h，4℃，12 000×g离心20 min。收集上清液，用作激素提取液[23]。采用ELISA试剂盒（上海酶联）测定内源激素脱落酸含量[24]。用酶标仪（SP-Max2300A2）在450 nm处测定吸光度，用标准曲线法计算待测激素的含量。以同时期健康叶片为对照，3次重复。

1.7 共表达网络的构建和模块鉴定

利用同时期转录组表达矩阵，过滤 FPKM值＜1的基因。使用R软件（R version 3.4.4）和WGCNA（R version 1.6.6）包，构建加权基因共表达网络与划分相关模块[25]。使用 WGCNA包中的 pickSoftThreashold计算软阈值，选取最佳网络构建的 power值。使用blockwiseModules构建无尺度网络，参数按照默认值设置[25]。

1.8 脱落酸及其下游信号转导关键模块鉴定和功能分析

利用moduleTraitCor＜-cor（mergedMEs，datTraits，use=“p”）对脱落酸含量和共表达模块进行相关性分析，鉴定关键模块。分析PYL和SnRK2基因家族在共表达网络中的分布情况，鉴定其信号转导调控的关键模块。提取模块中的基因，利用TBtools进行GO富集[19]，以谷子基因组为参考数据库，P-value＜0.05，利用Cytoscape对模块进行可视化[26]。

1.9 谷子PYL和SnRK2在禾生指梗霉胁迫下的表达模式分析

对转录组数据进行分析后，利用pheatmap对谷子PYL和SnRK2家族成员的表达水平进行可视化。

1.10 实时荧光定量PCR（qRT-PCR）

选取4个基因的转录本进行qRT-PCR验证。使用primer 3设计定量PCR引物，利用M-MuLV第一链cDNA合成试剂盒（TaKaRa）对提取的RNA进行纯化并反转录成cDNA。采用2×SG快速qPCR混合试剂盒（TaKaRa）进行荧光定量分析。qRT-PCR扩增的反应体积为15 μL，包括7.5 μL 2×SG快速qPCR主混合物、0.6 μL每个引物、40 ng cDNA模板、3.3 μL ddH2O。PCR流程：95℃变性30 s，95℃变性5 s，60℃退火/延伸30 s，40个循环，以actin为内对照基因，用2-ΔΔCT法计算相关基因表达。3次重复。

2 结果

2.1 谷子PYL和SnRK2的筛选

基因家族分析共鉴定到PYL和SnRK2家族基因各11个（表1）。谷子PYL和SnRK2家族基因除6号和8号染色体外，其余染色体均有分布，其中4号染色体上只有PYL基因Seita.4G239500，2号和7号染色体上只有SnRK2，分别为 Seita.2G394500和 Seita.7G100500，3号和9号染色体上数目最多，分别有6个。PYL家族基因氨基酸长度介于175—220 aa，等电点介于 4.39—8.88，分子量较小，分子量最大为23.67 kD。SnRK2家族基因氨基酸数目在333—454 aa，等电点介于 4.73—8.49，分子量较大，最高可达51.77 kD。

表1 谷子PYL和SnRK2基因家族成员基本信息Table 1 The basic information of PYL and SnRK2 gene family members in S.italica

2.2 谷子PYL和SnRK2的系统进化分析

谷子PYL家族在系统进化中分为4类，橙色分支中的 Seita.3G072600、Seita.5G302400和 Seita.5G140800较为保守，蓝色分支中的谷子PYL与水稻PYL聚集在一起，表明这部分谷子PYL可能与水稻PYL亲缘性较高。谷子SnRK2家族在系统进化中可分为3类，在灰色、紫色和黄色3个分支中，水稻SnRK2均与谷子SnRK2位于同一小分支上，而拟南芥SnRK2在分支上较独立，表明SnRK2在禾本科间较为保守（图1）。

2.3 顺式作用元件预测

谷子PYL和SnRK2家族基因顺式作用元件分析表明（图2），PYL和SnRK2家族基因中在127个区域预测到光响应元件。其中，22个家族基因都预测到脱落酸响应元件；谷子PYL家族基因Seita.1G013900、Seita.9G311900和SnRK2家族基因Seita.3G230400和Seita.3G369900被预测到防御和胁迫响应元件；PYL和SnRK2家族基因同时预测到其他激素类响应元件，其中与茉莉酸甲酯响应有关的基因最多，PYL家族9个，SnRK2家族10个。

2.4 脱落酸含量的测定

受侵染不同时期叶片中脱落酸含量测定表明，脱落酸在病菌侵染前期显著积累，TG1和TG2时期脱落酸含量显著高于对照组，分别为 22.50和 18.08 ng·mL-1，随后含量显著减少。在TG3时期降低到与对照组同等水平13.79 ng·mL-1，在TG4时期达到最低值13.18 ng·mL-1，且含量显著低于对照组 21.82 ng·mL-1（图3）。

2.5 谷子脱落酸共表达网络鉴定

利用 gsg=goodSamplesGenes（datExpr，verbose=3），gsg$allOK参数测试数据的合理性，在所测RNA-seq数据中筛选软阈值β=12构建共表达网络（图4-A），利用18 535个基因共得到34个基因共表达模块（图4-B）。其中MEturquoise模块的基因数目最多有3 924个，MEdarkolivegreen模块的基因数目最少，仅有43个。

图1 谷子PYL和SnRK2家族系统进化树Fig.1 Phylogenetic tree of PYL and SnRK2 family genes in S.italica

图2 谷子PYL和SnRK2家族基因的启动子顺式作用元件分析Fig.2 Functions of promoter cis-acting elements in PYL and SnRK2 family genes in S.italica

图3 脱落酸含量测定Fig.3 Determination of abscisic acid content

2.6 谷子脱落酸共表达网络鉴定

通过共表达网络与脱落酸含量进行关联，鉴定到与脱落酸高度正相关的模块为MEpaleturquoise模块，相关性值为 0.56，其次为 MEbrown和 MEdarkolivegreen模块，相关性值均为0.48（图5-A）。其中负相关性中值最高的为 MEpink模块，相关性值为-0.54（图5-B），相关性绝对值最低的为MEdarkorange模块，相关性值为0.024。对脱落酸信号转导途径中的SnRK2和PYL家族基因所在模块的鉴定中发现（图5-C），SnRK2和PYL所在模块数目最多的为MEbrown模块，占总数目的27%，其中该模块中有4个PYL和2个SnRK2，表明该模块也是脱落酸调控及发挥生理效应的关键模块。在MEturquoise模块中有5个SnRK2，但该模块与脱落酸含量的相关性绝对值仅为0.026，不存在相关性。有5个PYL和1个SnRK2不属于任何模块。故MEpaleturquoise和MEbrown为脱落酸及其信号转导响应禾生指梗霉侵染的关键模块。

图4 软阈值确定和共表达网络构建Fig.4 Determination of soft threshold and construction of co-expression network

2.7 模块的功能富集

对筛选的MEpaleturquoise和MEbrown模块进行功能富集分析（表2），其中MEpaleturquoise模块显著富集到 1-3-beta-D-葡聚糖生物合成过程（GO:0006075）、类黄酮代谢过程（GO:0051552）和类黄酮生物合成过程（GO:0051553），还富集到了钙离子平衡（GO:0055074）和磷酸转移酶活性（GO:0016780）等。MEbrown模块富集到响应真菌（GO:0009620）、生物刺激反应（GO:0009607）、免疫系统的过程（GO:0002376）等与胁迫相关的功能；激素响应（GO:0009725）、激素介导的信号通路（GO:0009755）、乙烯激活的信号通路（GO:0009873）和油菜素内酯（GO:0009742）介导的信号通路等与激素相关的一些功能；以及蔗糖响应（GO:0009744）、脂肪酸代谢（GO:0006631）和调节氮化合物代谢过程（GO:0051171）等与基础代谢相关的功能。

对与脱落酸关联的 MEpaleturquoise模块和MEbrown模块中以谷子PYL、SnRK2家族基因构建了共表达网络。其中 MEplaeturquoise中鉴定到 Seita.4G105600、Seita.6G218100和Seita.9G138400为该模块的枢纽基因（图6-A），Seita.4G105600在水稻和拟南芥中都被鉴定为含有 WD结构域的超家族蛋白，Seita.6G218100在水稻中的同源基因为WRKY57转录因子，在拟南芥中同源基因为包含VQ基序的蛋白；Seita.9G138400在水稻和拟南芥中都注释的为植物特有的响应胁迫相关的转录因子TIFY（表3）。在脱落酸信号转导途径中的PYL和SnRK2家族基因有4个处于网络调控的枢纽位置（图6-B），其中PYL家族基因的 Seita.1G013900和SnRK2家族基因的 Seita.2G394500在网络中权重值较高，可能起主要调控作用。

2.8 表达模式分析

在响应禾生指梗霉侵染的表达模式分析中，处理组中的PYL家族基因Seita.1G030500、Seita.3G207900、Seita.9G311900和Seita.9G437300均在前期表达水平较高，SnRK2家族基因Seita.2G394500、Seita.3G003200在TG1时期表达水平较高（图7-A）。对谷子PYL和SnRK2家族基因与脱落酸含量进行相关性分析，鉴定到基因Seita.1G030500、Seita.2G394500和Seita.3G003200显著性最高，分别为0.80、0.65和0.72（图7-B）。

图5 关键模块的鉴定Fig.5 Identification of key modules

图6 MEpaleturquoise模块内基因和MEbrown模块内PYL、SnRK2家族基因共表达网络Fig.6 Gene co-expression network of the MEpaleturquoise module and PYL, SnRK2 family genes in MEbrown module

表2 模块GO富集情况Table 2 GO enrichments of network modules

表3 模块核心基因功能注释Table 3 Functional annotation of modular hub genes

2.9 实时荧光定量PCR

4个基因 Seita.2G394500、Seita.4G105600、Seita.6G218100和Seita.9G138400的荧光定量分析表明，基因Seita.2G394500在TG1、TG2时期的表达量显著上调；基因Seita.4G105600和Seita.6G218100在TG2、TG3、TG4时期的表达量均高于对照组的表达量，且TG2、TG3时期表达量显著上调；基因Seita.9G138400 5个时期表达量显著下调（图8）。

图7 PYL和SnRK2家族基因表达模式及相关性分析Fig.7 Gene expression pattern and correlation analysis of PYL and SnRK2 families

3 讨论

脱落酸可作为抗病性的抑制因子或增强剂，VLEESSCHAUWER等[27]研究发现，通过对水稻叶片施加脱落酸，可增加水稻对褐斑病的抗性，但ULFERTS等[28]发现脱落酸含量的升高导致稻瘟病菌对大麦叶片的危害加重。脱落酸在植物的防御反应中扮演着复杂而矛盾的角色。本研究中，在禾生指梗霉早期侵染阶段，脱落酸含量较对照组显著增加，表明禾生指梗霉的侵染可能引起寄主内源脱落酸含量的积累。鉴定出的 11个谷子SnRK2基因被分为3类，与水稻[29]、拟南芥[30]的研究结果一致，而且与玉米[31]、高粱[32]的研究结果类似，主要包括不依赖于脱落酸、不依赖于激酶或对脱落酸的依赖性弱和对脱落酸具有很强的依赖性 3类。谷子Seita.1G013900、Seita.4G239500基因与水稻中OsPYL3（LOC_Os02G15640）和OsPYL9（LOC_Os06g36670）同源性较高，这两个基因在水稻中可以正向调控脱落酸信号，是真正的水稻脱落酸受体[33]，通过PYL和SnRK2的同源聚类分析验证了其在脱落酸信号转导通路上的重要作用。

图8 不同时期基因表达量变化Fig.8 Changes in gene expression at different stages

在共表达网络分析中，筛选到与脱落酸含量高度相关的MEpaleturquoise模块，通过对该模块的富集分析鉴定到了槲皮素、类黄酮和钙离子平衡相关的功能。贾振华[34]用槲皮素处理野生型拟南芥后，活性氧含量增加，同时还出现活性氧爆发、过敏性细胞坏死、胼胝质沉积等多种抗性反应，进而增强拟南芥对丁香假单胞菌（Pseudomonas syringae）的抗性。PYL和SnRK2在脱落酸信号转导过程中具有重要作用，在顺式作用元件分析中，鉴定到PYL和SnRK2家族基因中存在与防御和胁迫响应、脱落酸、水杨酸、茉莉酸甲酯、类黄酮等响应元件，表明PYL和SnRK2家族基因在参与激素响应胁迫方面可能发挥重要作用。在对PYL和SnRK2高度相关的MEbrown模块富集中，富集到了大量与生物胁迫相关功能的基因，如免疫系统过程、细胞壁大分子生物合成过程、脂肪酸代谢、油菜素内酯、乙烯及生物刺激反应等。脂肪酸代谢物甘油-3-磷酸（G3P）的积累可以提高拟南芥对炭疽病菌（Colletotrichum higginsianum）的抗性[35]，植物角质层中的主要成分也是脂肪酸及其衍生物，对抵御病菌侵染具有重要作用[36]。油菜素内酯可诱导水稻对烟草花叶病毒（Tobacco mosaic virus，TMV）、稻瘟病菌和水稻白叶枯病菌（Xanthomonas oryzae）的抗性[37]。通过对谷子不同侵染、扩展阶段脱落酸含量的测定、生物信息学分析及不同时期表达量的测定，推测脱落酸在禾生指梗霉侵染过程中可能发挥着非常重要的作用。

本研究通过表达模式和相关性分析，鉴定到基因Seita.1G030500、Seita.2G394500和Seita.3G003200在脱落酸响应过程中发挥着重要作用。通过共表达网络分析预测到的3个核心基因中，Seita.4G105600在水稻和拟南芥中的同源基因均为含有WD40重复结构的家族蛋白，WD40家族蛋白在植物生长发育、响应胁迫方面具有重要作用，李宝燕[38]鉴定到含有WD40基序的NtTTG2通过协调生长素和水杨酸提高抗病性。Seita.6G218100在水稻中的同源基因为WRKY57，WRKY转录因子在番茄上被报道为负调控番茄对灰霉病的抗性[39]。Seita.9G138400的同源基因在水稻和拟南芥中均被鉴定为 TIFY家族蛋白基因，TIFY家族蛋白在植物生长发育、激素应答和逆境胁迫方面具有重要作用，过表达VvTIFY9可提高葡萄对白粉病的抗性[40]。笔者推测，谷子内源脱落酸在受到禾生指梗霉胁迫时，可能通过 Seita.6G218100、Seita.4G105600和 Seita.9G138400基因编码的WRKY57、WD40重复结构蛋白以及TIFY蛋白激活抗病响应途径。

4 结论

脱落酸在禾生指梗霉侵染谷子叶片早期大量积累，通过加权基因共表达网络进行脱落酸及其信号转导基因PYL和SnRK2关联分析，共筛选到6个基因参与谷子内源脱落酸响应禾生指梗霉侵染过程，包含 1个PYL家族基因、2个SnRK2家族基因、1个WD40家族蛋白基因、1个WRKY57转录因子和1个TIFY转录因子基因。经qRT-PCR验证发现1个SnRK2家族基因、1个WD40家族蛋白基因和1个WRKY57转录因子基因共3个基因可能在谷子脱落酸响应禾生指梗霉侵染过程中发挥重要作用。