农学TOP期刊 ▏90个样本的DNA甲基化谱揭示了胡杨不同异形叶片发育的动态表观遗传调控

• 中文标题：DNA甲基化谱揭示了胡杨不同异形叶片发育的动态表观遗传调控

• 期    刊：Industrial Crops & Products

• 影响因子：5.9

• 发表时间：2024年5月14日

• 主要作者单位：塔里木大学、浙江师范大学、中南民族大学

• 通讯作者：李志军教授、吴智华副教授

• 研究背景

胡杨（Populus euphratica）是干旱地区的重要树种和沙漠环境的重要生态组成部分，也是研究木本树木环境适应性的典型模型。不同树冠上生出异形叶是胡杨引人注目的一个特点，异形叶按Li（条形）、La（披针形）、Ov（卵形）和Bo（阔卵形）从上至下连续性分布，四种典型异形叶片以叶长（LL）和叶宽（LW）的比值划分，这被认为是杨树在长期进化过程中对干旱环境的一种适应性特征。

最近的研究主要集中在胡杨异形叶之间的功能和结构差异上。与窄叶（Li和La）相比，宽叶（Ov和Bo）具有增强抗旱能力的特征，包括更厚的栅栏组织、更大的气孔密度和尺寸以及更稳定的持水能力。胡杨异形叶的这些差异特征在多个分子水平上受到调控，如转录、翻译和群体遗传学。然而，随着树木从幼年发育到成年，异形叶的整体表观遗传调控仍然未知。研究异形叶的DNA甲基化将有助于加深对其结构和功能特征以及生态适应性分子机制的理解。

2024年5月14日，塔里木大学李志军教授团队及其合作者团队联合在农学TOP期刊Industrial Crops & Products上发表题为“DNA methylation profile revealed the dynamically epigenetic regulation of  the distinct heteromorphic leaf development in Populus euphratica”的研究性论文。研究通过分析90个异形叶样本的DNA甲基组和转录组，揭示了DNA甲基化与基因表达之间的动态调控景观，以及异形叶的发生和发育，并在树龄和叶期的异形叶中鉴定出一组由DNA甲基化特异性调控的基因。研究结果为胡杨异形叶形态形成的表观遗传调控机制提供了一个新的视角。

• 研究路线

 研究结果 

1. 胡杨异形叶发育过程中DNA甲基化的全基因组分析

研究对总共90个样本进行了WGBS分析以探索胡杨叶片发育过程中的全局DNA甲基化景观（图1a）。在每个叶阶段，异形叶的19条染色体的基因密度和甲基化水平之间存在显著的负相关（图1b）。与几乎所有植物物种一样，胡杨的DNA甲基化模式也分为CG、CHG和CHH三种类型。其中，CG甲基化水平最高，其次是CHG，最后是CHH，与它们的甲基化比例完全相反。三种甲基化背景的分布表明，P1的甲基化水平低于P2和P3，Ov和Bo叶的甲基化水平低于Li和La叶，但在不同的树龄中没有观察到显著的甲基化模式（图1c）。

为了更清楚地了解异形叶跨叶期的甲基化模式，将异形叶分为窄叶（Li和La）和宽叶（Ov和Bo）两类。值得注意的是，随着叶期的增加，窄叶和宽叶甲基化水平的差异逐渐减小。P1阶段，窄叶和宽叶之间的甲基化差异最大，尤其是CG和CHG甲基化模式，分别达到5.51%和5.54%（图1d）。这些结果表明，DNA甲基化可能在胡杨早期叶阶段异形叶的调控中起着至关重要的作用。

图1. 胡杨异形叶片发育的DNA甲基化景观。（a）胡杨4个树龄在3个发育叶期的不同叶形。（b）P1、P2和P3时期胡杨异形叶中各样品基因密度和甲基化水平的Circos图。a代表基因密度；b-j分别表示A_Li、B_Li、C_Li、D_Li、B_La、C_La、D_La、C_Ov、D_Ov和D_Bo。（c）不同比较组的CG、CHG和CHH背景的DNA甲基化水平。（d）P1、P2和P3期窄叶和宽叶三种背景DNA甲基化水平的差异。

2. 基因和TE区域的DNA甲基化模式

基因的DNA甲基化对于转录的表观遗传调控非常重要，而转座因子（TEs）的超甲基化限制了它们的流动性，增强了基因组的稳定性。分析所有样本的编码区域和TE区域的DNA甲基化水平，观察到三种背景的甲基化都在基因区表现出一致的趋势，即朝着转录起始位点（TSS）下降，以钟形方式覆盖基因体，并在转录终止位点（TSS）下游增加（图2a）。在TE区域，与上游TSS和下游TTS相比，TE体内所有三种背景的甲基化水平都有所升高（图2b）。

有趣的是，与CG和CHG甲基化相比，CHH甲基化在异形叶的发育过程中表现出更明显的基因和TEs差异（图2a-b）。长末端重复反转座子（LTR）是植物基因组中数量最多的TEs，研究对不同类型的LTRs的DNA甲基化模式进行了比较分析，重点关注Gypsy和Copia LTRs。值得注意的是，Gypsy和Copia LTRs的DNA甲基化谱表现出惊人的相似性。然而，对于所有三种背景，Copia LTRs在TSS上游和TTS下游的甲基化水平低于Gypsy LTRs，可能表明Copia在胡杨基因组中具有更高的转录活性。

据报道，基因启动子区及基因体区的DNA甲基化与基因表达密切相关。此研究发现窄叶中的DNA甲基化率高于宽叶，除了P3中CHG背景的（图2a）；而在基因体和TE体区域的所有甲基化类型中，无论叶的发育阶段如何，窄叶始终表现出比宽叶更高水平的甲基化（图2a-b）。从P1到P3，窄叶和宽叶在所有基因和TE背景下的甲基化水平差异逐渐减小，直到在P3差异不明显（图2c）。因此，窄叶和宽叶甲基化水平的显性差异主要发生在早期叶片阶段（P1和P2），这是异形叶形成和扩展的关键阶段。

图2. 胡杨异形叶的DNA甲基化模式。（a）编码基因的全局DNA甲基化景观。柱状图显示了上游TSS、基因体区域和下游TTS内的平均甲基化水平。（b）TEs的全局DNA甲基化景观。柱状图显示了TE体内的平均甲基化水平。（c）不同背景下宽叶和窄叶在不同叶期基因和TE甲基化水平的差异。

3. 叶形比较组中关键甲基化基因的鉴定

为了更好地阐明窄叶和宽叶DNA表观遗传学的动态差异，将30个对照组分为组间比较[(Ov or Bo) vs . (Li or La)]和组内比较[(La vs . Li)和(Bo vs . Ov)]两类。结果显示，组间差异甲基化区域（DMRs）和差异甲基化基因（DMGs）的数量大于组内差异的，尤其是P1和P2阶段，数量差异显著（图3a-b）。此外，组间DMRs和DMGs均值数从P1-P3逐渐减少，而组内的DMRs和DMGs均值数几乎不变，尽管在各比较组中存在动态差异（图3b）。进一步的分析确定了P1和P2时期D树龄中DMRs的定位，揭示了TSS上游和TTS下游DMRs数量的显著增加。这表明DNA甲基化主要通过转录调节元件的调节影响基因表达（图3c），强调了其在叶形态调节中的关键作用。

P1和P2期的C树龄和D树龄的DMGs存在显著重叠。在P1，分别有2203和1406个DMGs在树木C期和树木D期中重叠（图3d）；而在P2，相应的数量变为1712和525个。GO富集分析显示，P1期的重叠DMGs显著富集“微管功能”（图3e），一些微管运动激酶基因如PeKINESIN-12B在阔叶中上调，它在细胞分裂过程中对维持细胞骨架的完整性起着关键作用。此外，增强叶片抗旱性结构的基因，如PeCUS2和PeRab-H1b也在阔叶中上调表达。值得注意的是，术语“叶片发育”显著丰富（图3e），其中参与多种胁迫的PeDGK2在窄叶中上调。此外，大多数参与吲哚、萜类和类固醇次生代谢生物合成过程的DMGs在窄叶中表达上调（图3e-f）。这些结果表明，窄叶通过次生代谢对逆境胁迫做出响应，而宽叶通过发育功能结构对逆境胁迫做出适应。

图3. DMRs和DMGs的表征。（a）成对比较鉴定DMRs和DMGs。（b）组间比较和组内比较的dmr和dmg数量。*p < 0.05;** p < 0.01;*** p < 0.001。（c）DMRs在组间比较沿TSS、基因体和TTS的分布。（d）P1期的C树龄和D树龄DMGs重叠维恩图。P1中重叠DMGs的GO（e）和KEGG（f）富集分析。

4.叶期和树龄比较组中的DMRs和DMGs鉴定

统计胡杨不同叶期和树龄比较组中的DMRs和DMGs数量，发现在叶期比较组中，P2 vs. P1中DMRs和DMGs的数量高于P3 vs. P2（图4a），表明甲基化变化在叶片发育初期更为明显。引人注目的是，D树龄的宽叶比窄叶显示出明显更高的DMRs丰度（图4a）。此外，各叶期的比较组中，hyper-DMRs占主导地位（图4b）。调节叶面积的PeGRF2的基因表达在Ov成熟过程中降低，并与启动子区域的DNA甲基化呈负相关（图4c）。相反，该基因的表达和甲基化水平在La叶片成熟过程中保持相对稳定（图4c）。这些结果表明，DNA甲基化可以影响PeGRF2的表达，从而调节胡杨叶片的发育。然而，DMRs和DMGs的数量在不同树龄没有显著变化，hyper-DMRs和hypo-DMRs的分布相似，表明树龄对DNA甲基化的影响很小。

图4. 叶期比较组的差异甲基化区域（DMRs）。（a）DMRs数目。（b）P2 vs. P1和P3 vs. P2中hyper-DMRs和hypo-DMRs的比例。（c）6个样本中PeGRF2的DNA甲基化基因组浏览图。（d）PeGRF2在6个样本中的表达情况。

5. DNA甲基化与基因表达的联合分析

为了鉴定受甲基化影响的调节叶片形状的潜在基因，对DNA甲基化和转录组数据进行了联合分析。首先，研究鉴定了叶形组间和组内比较的差异表达基因（DEGs）。与图3a所示DMRs的分布模式相似，组间比较DEGs数量明显高于组内比较（图5a-c）。鉴于C树龄和D树龄树木都有宽叶和窄叶，这项研究集中于此，并分别在P1、P2和P3三个叶期中鉴定到C树龄和D树龄的871、2222、609、1373、205和651个叶形DMGs和DEGs重叠基因（图5d-f）。将C树龄和D树龄下的各重叠基因合并后，P1、P2和P3分别有2627、1784和805个基因重叠（图5d-f），所有这些基因被进一步用作功能富集分析。它们中的大部分在“蛋白质泛素化”、“DNA复制”、“对激素的反应”、“脂肪酸生物合成过程”等GO功能类别上显著富集（图5g），对“对激素的反应”的11个基因在所有样本中的表达量被详细分析。值得注意的是，作为生长素激素信号通路核心的6个PeARF基因在Ov和Bo叶片中显著上调（图5h），反映出生长素可能参与了异形叶片的生长过程。

图5. 基因表达与甲基化之间的关系。P1（a）、P2（b）和P3（c）中上调和下调的DEGs数量。维恩图显示了P1（d）、P2（e）、P3（f）叶形比较中DMGs和DEGs的重叠基因。（g）三个叶期所有重叠基因的GO富集。（h）所有样品中激素相关基因的表达。

6. 生长素相关基因调控胡杨异形叶片发育

在上述研究结果的基础上，利用BLASTp从胡杨基因组中鉴定了40个拟南芥ARF同源基因。其中，有24个基因在宽叶中的表达量高于窄叶，而大约一半的PeARF基因在甲基化和基因表达之间表现出负相关（图6a），这意味着在生长素信号转导中涉及多种调节机制。先前已有研究表明，在叶片发育中起关键作用的ARF5直接与WOX1启动子结合并激活其表达，WOX1在叶片边缘区域表达，使叶片变平形成宽叶片。此研究发现，两个PeARF5基因在阔叶中唯一高表达，其下游的PeWOX1基因也表现出类似的表达模式（图6a-c）。这些结果表明，与窄叶相比，宽叶的叶面积更大可能是由于其生长素信号通路更活跃。

为了进一步研究PeARF表达上调的原因，还深入研究了AUX/IAA-TIR1核信号通路，还分别鉴定出33个PeAUX/IAA、3个PeAFB和4个PeTIR1基因。先前的研究表明，TIR1/AFB感知生长素导致AUX/IAA蛋白的降解，这反过来释放ARF蛋白，用于生长素调控基因的转录（图6d）。在这里，PeARF5的表达水平上调，而抑制PeARF5表达的上游基因PeAUX/IAA的表达也上调，说明PeARF5的上调可能受到AUX/IAA蛋白泛素化的影响（图6d）。此外，PeTIR1在阔叶中不高表达，而PeAFB在阔叶中高表达（图6a）。因此，研究提出了PeAFB升高加剧了PeAUX/IAA泛素化降解，导致AUX/IAA含量降低后不能抑制PeARF5的表达，从而促进了下游PeWOX1基因的表达，最终导致叶片变宽的假设。值得注意的是，PeARF5的甲基化水平和基因表达在窄叶和宽叶之间表现出显著差异（图5g-f）。此外，当比较不同的叶形状时，鉴定了PeARF5、PeIAA16和PeAFB2基因上游的稳定DMR，表明DNA甲基化在胡杨异形叶的发育过程中起关键作用。

图6. 生长素信号通路相关基因的表达模式。（a）P1中ARF、AUX/IAA、TIR1和AFB基因的表达水平。“★”和“💧”分别表示DNA甲基化与基因表达呈正相关和负相关。（b）P1时WOX1基因在所有样品中的表达水平。（c）D树龄P1 期WOX1基因的qRT- PCR和RNA-seq结果。（d）提出的胡杨异型叶片发育生长素信号通路模型。

 结  论 

总之，本研究提供了胡杨异形叶生长过程中DNA甲基化模式的综合分析。结果表明，不同叶形和不同时期的DNA甲基化水平不同。基因和TE区的DNA甲基化分析显示出不同的模式，其中CHH甲基化表现出更明显的差异。此外，研究还揭示了生长素激素信号通路参与异形叶片的形成，DNA甲基化密切调节生长素反应因子的表达。总的来说，本研究扩大了我们对表观遗传机制参与异形叶片形成的理解，并为胡杨的适应进化提供了新的见解。