茶树中的动态DNA甲基化及其与盐碱胁迫下基因表达变化的相关性

《Molecular Horticulture》：Dynamic DNA methylation in tea plants and its association with changes in gene expression under salt and alkali stress

茶树在盐碱胁迫下基因表达动态及表观遗传调控机制研究。通过RNA测序和全基因组bisulfite测序，发现胁迫延长导致差异表达基因（DEGs）数量显著增加，KEGG通路分析显示与黄酮类合成、光合作用及MAPK信号通路相关。甲基化水平变化与基因表达呈负相关，特定区域（如基因体、2Kb上下游）的甲基化修饰介导了胁迫响应基因的调控，并鉴定CsDRM2等关键基因参与甲基化水平升高。

　　

盐分和碱分胁迫会诱导茶树（Camellia sinensis (L.) O. Kuntze）的基因表达变化，进而影响植物的生长发育（Wan等人 2024；Zhang等人 2017, 2023）。近期在表型组学和转录组学方面的进展揭示了茶树幼苗在盐胁迫下的生理和基因表达变化（Wan等人 2018, 2020, 2024；Zhang等人 2017, 2023）。然而，这些基因表达变化背后的表观遗传调控机制仍不清楚。因此，本研究探讨了盐分和碱分胁迫下转录和5mC甲基化修饰的动态变化。

在RNA测序（RNA-seq）实验中，一年生的茶树幼苗分别用盐溶液（200 mM NaCl）或碱溶液（150 mM NaHCO₃）处理0天（处理前）、3天和7天。盐分和碱分胁迫处理后，幼叶出现明显萎蔫现象，且随着时间的推移，萎蔫程度、过氧化物酶（POD）和超氧化物歧化酶（SOD）的含量逐渐增加（图 1A，S1）。该研究从15个样本中获得了99.79 Gb的干净RNA-seq数据（每个处理组有三个生物学重复），平均读段映射率为93.59%（表S1和图S2、S3）。RNA-seq分析显示，与处理3天的样本相比，处理7天的样本中差异表达基因（DEGs）的数量增加了4倍以上（NC-7d和NHC-7d vs CK-0d）。这一发现表明，较长的胁迫时间会导致更广泛的基因表达变化（图 1B）。相应地，处理后第7天的植物表现出比第3天更严重的萎蔫现象（图 1A）。京都基因与基因组百科全书（KEGG）对DEGs的功能注释显示，所有组中下调的DEGs在苯丙素类和黄酮类生物合成途径（图 1C）以及光合作用和DNA复制途径（图 1C和D）中显著富集，表明盐分和碱分胁迫影响了茶树的黄酮类和儿茶素合成及光合作用系统。相反，上调的DEGs在二萜类和油菜素内酯生物合成途径（图 1C）以及几个已知的非生物胁迫相关途径中显著富集，包括MAPK信号传导、植物激素信号传导和环境适应途径（图 1D）。

图1

A 茶树幼苗在盐分和碱分胁迫下的形态变化（CK、NC和NHC分别代表对照、NaCl和NaHCO₃处理组）。每个处理组包含三个独立的生物学重复。 B NC-3d与CK-0d（NC-3d）、NC-7d与CK-0d（NC-7d）、NHC-3d与CK-0d（NHC-3d）以及NHC-7d与CK-0d（NHC-7d）之间的差异表达基因（DEGs）数量。Up和Down分别代表上调和下调的DEGs。 C和D KEGG注释结果（P值<0.05）显示的差异DEGs（差异比较的缩写与图1B相同）。 E和F 盐分和碱分胁迫下基因及转座子（TEs）本体及其上下游2 Kb区域的甲基化水平变化（处理组的缩写与图1A相同）。 G和H NC-3d、NC-7d、NHC-3d和NHC-7d中由五种不同DMR介导的下调DEGs（Down-DEGs）和上调DEGs（Up-DEGs）的数量（差异比较的缩写与图1B相同）。 I 五种不同DMR介导的下调DEGs或上调DEGs的KEGG注释（差异比较的缩写与图1B相同）。 J 黄酮类合成途径和MAPK信号途径中几种DMR介导的DEGs的表达倍数变化（Log₂处理后）。 K CsTGY05G0002034（CHS）、CsTGY06G0001806（PYL）和CsTGY11G0001022（PP2C）的相对基因表达结果（qRT-PCR分析）。 *** 表示

值<0.001（配对t检验）（处理组的缩写与图1A相同）。 L 与茶树基因组中的甲基化和去甲基化相关的基因数量。 M 盐分和碱分胁迫下与甲基化和去甲基化相关基因的表达变化（处理组的缩写与图1A相同）

全尺寸图像

维恩图分析将DEGs分为两类：核心DEGs，指在处理后第3天和第7天相对于第0天持续上调或下调的基因；以及时期特异性DEGs，指仅在处理后第3天或第7天相对于第0天差异表达的基因。在碱分胁迫下，该分析识别出6610个时期特异性DEGs和1098个核心DEGs（上调：545 + 下调：553）（图S4A）。同样，在盐分胁迫下，检测出3837个时期特异性DEGs和331个核心DEGs（上调：166 + 下调：165）（图S4B）。进一步分析这四组核心DEGs发现，在碱分胁迫下有413个下调核心DEGs和483个上调核心DEGs，在盐分胁迫下有25个下调核心DEGs和104个上调核心DEGs，在盐分和碱分胁迫条件下共有的有140个下调核心DEGs和62个上调核心DEGs（图S4C）。KEGG结果显示，这些核心DEGs参与了不同的途径，DEGs的数量与途径参与度之间存在相关性（图S5）。四个已知的非生物胁迫相关途径——信号传导、植物激素信号传导、碳水化合物代谢和淀粉及蔗糖代谢——在所有六组核心DEGs中均显著富集，尽管涉及的基因数量有所不同（图S5）。

此外，我们对这六组核心DEGs进行了转录因子（TF）鉴定分析，分别在NC.Core.Down、NC.Core.Up、NHC.Core.Down、NHC.Core.Up以及NC和NHC.Core.Up DEG组中识别出10、56、14、5和7个TF编码的核心DEGs（图S4C和图S6）。在先前的研究中，其中一些TF已在盐分和碱分胁迫下通过实验验证了其功能效应（图S4D），包括WRKY33和ERF98（Jiang和Deyholos 2009；Zhang等人 2012）。此外，我们还在ABA信号传导途径中识别出5个核心DEGs，在MAPK级联途径中识别出2个核心DEGs，在Ca²⁺信号传导途径中识别出9个核心DEGs，在14个核心DEGs中识别出编码不同转运蛋白的基因（图S4E）。这些发现表明，维持盐分和碱分胁迫下的稳态涉及多种转运蛋白。

考虑到茶树基因组的高度杂合性质，我们开发了一种结合全基因组亚硫酸盐测序（WGBS-seq）和重测序的新方法（图S7），能够在之前通过RNA-seq分析的相同样本上进行DNA甲基化测序（表S2）。WGBS-seq分析显示，CG型、CHG型和CHH型的甲基化水平分别为89.51%–90.88%、69.37%–70.65%和7.80%–8.58%（图S8），这与之前的研究结果不同（CG、CHG和CHH甲基化水平分别为84.88%、68.41%和12.60%（Kong等人 2023）。这种差异主要源于我们改进的方法，该方法能够纠正由全基因组C/T杂合位点引起的假阳性甲基化量化。分析显示，处理样本中的所有类型甲基化水平均显著升高（图S8），表明盐分和碱分胁迫导致了甲基化修饰的显著变化。进一步比较基因及转座子（TEs）本体、上游和下游区域的甲基化水平发现，基因及其上游和下游区域以及TEs的本体和侧翼区域的甲基化水平发生了显著变化（图 1E和F）。这些结果表明，茶树可能通过调节基因和TEs的甲基化水平来应对盐分和碱分胁迫。差异甲基化分析和注释结果显示，高甲基化差异位点（hyper-DMPs）和高甲基化差异区域（hyper-DMRs）的数量显著多于低甲基化DMPs（hypo-DMPs）和低甲基化DMRs（hypo-DMRs）（图S9A）。DMRs主要位于基因间的远端和启动子区域，仅有少量位于基因编码区域（图S9B）。

为了研究DNA甲基化对基因表达的影响，本研究分析了基因区域及其上游和下游2 Kb区域内mCG、mCHG和mCHH水平与基因表达水平之间的关系。分析显示，mCG在基因下游2 Kb区域、mCHG在基因本体、mCHH在基因下游2 Kb区域以及mCHH在基因本体中的表达水平与甲基化水平呈负相关（图S10）。值得注意的是，mCHH在基因上游2 Kb区域的甲基化水平与基因表达水平呈正相关（图S10），这与moso竹（Ding等人 2022）和苹果（Xu等人 2018）的研究结果一致。基于这些相关性，研究识别出特定区域内由DMR介导的DEGs：mCG在基因下游2 Kb区域、mCHG在基因本体、mCHH在基因上游2 Kb区域以及mCHH在基因本体。分析分别识别出NC-3d vs CK-0d（NC-3d）、NC-7d vs CK-0d（NC-7d）、NHC-3d vs CK-0d（NHC-3d）和NHC-7d vs CK-0d（NHC-7d）中由这五种不同DMR介导的22、206、23和95个非冗余Down-DEGs（图 1G和表S3）。相应地，NC-3d、NC-7d、NHC-3d和NHC-7d中分别由五种不同DMR介导的33、187、22和158个非冗余Up-DEGs（图 1H）。KEGG注释突出了差异甲基化修饰对黄酮类合成和胁迫响应途径相关基因表达变化的影响（图 1I），包括FLS（黄酮醇合成酶）和CHS（查尔酮合成酶）（与黄酮类合成相关）在Down-DEGs中，以及PYL（脱落酸受体PYR/PYL家族）（Cutler等人 2010；Ren等人 2022）、MAPKKK17/18（丝裂原活化蛋白激酶激酶17/18）（Chardin等人 2017；Zhang等人 2018）和PP2C（蛋白磷酸酶2C）（Cutler等人 2010）（与MAPK信号途径相关）在Up-DEGs中（图 1J、S11、K）。此外，三种核心DEGs，即WRKY33、MPK3和ABC转运蛋白，同时受到盐分和碱分胁迫的影响（表S4）。

通过对单倍体TGY基因组进行BLASTP和HMM算法分析，识别出了编码去甲基化酶、甲基转移酶和组蛋白修饰因子的基因（图 1L和M）。RNA-seq和qRT-PCR分析（引物见表S5）显示，一个甲基化相关基因（CsDRM2）的表达上调有效解释了胁迫下甲基化水平的增加（图 1M、S12），表明CsDRM2的上调可能是胁迫下甲基化水平增加的主要原因。然而，考虑到某些时间点甲基化相关基因的表达变化也可以解释甲基化水平的变化，因此全基因组加权甲基化的增加也可能是由多个甲基化相关基因的共同作用所致。

总之，这些发现增强了我们对非生物胁迫下茶树表观遗传机制的理解，并确定了通过表观遗传变异改良未来育种的重要基因靶标。