中国科学家方磊/谷晓峰团队通过群体DNA甲基化多态性分析 揭示表观遗传育种在作物改良中的巨大潜力(国人佳作) ：Cell Res/IF25.9

近日，浙江大学农业与生物技术学院方磊教授和张天真教授、中国农业科学院生物技术研究所谷晓峰研究员、阿里巴巴达摩院（湖畔实验室）顾斐博士团队合作，通过对207份棉花材料开花后20天（DPA）纤维进行单碱基DNA甲基化测序（WGBS）、转录组测序（RNA-seq）和全基因组测序（WGS），利用GWAS、eQTL、EWAS、eQTM、深度学习模型DeepFDML等算法，将群体遗传学的经典框架扩展到表观遗传学领域。研究揭示了DNA甲基化多态性（SMPs）与基因表达和纤维性状的关联，并证明了DNA甲基化数据在作物育种中的重要应用价值。相关研究成果以“Population-wide DNA Methylation Polymorphisms at Single-nucleotide Resolution in 207 Cotton Accessions Reveal Epigenomic Contributions to Complex Traits”为题发表在国际学术期刊《Cell Research》。

标题：Population-wide DNA Methylation Polymorphisms at Single-nucleotide Resolution in 207 Cotton

Accessions Reveal Epigenomic Contributions to Complex Traits（207份棉花的群体DNA甲基化多态性揭示了表观基因组对复杂性状的贡献）

发表时间：2024年10月17日

发表期刊：Cell Research

影响因子：IF 25.9/Q1

作者单位：浙江大学农业与生物技术学院

组学技术：WGBS、RNA-seq、WGS（易基因金牌技术）

分析算法：eQTL、eQTM、DeepFDML等

DOI: 10.1038/s41422-024-01027-x

DNA甲基化在作物发育中发挥多种调控作用，但自然作物群体甲基化多态性与遗传多态性、基因表达和表型变异之间的关系仍不清楚。本研究对207份棉花材料开花后20天（DPA）纤维的单碱基DNA甲基化组（methylomes）、转录组（transcriptomes）和基因组（genomes）进行了整合分析，并将群体遗传学的经典框架扩展到表观遗传学领域。研究共鉴定出超过2.87亿个单甲基化多态性（SMPs），这一数量是单核苷酸多态性（SNPs）的100倍。这些SMPs显著富集在基因内区域，而在转座元件中则相对缺失。关联分析进一步鉴定出总计5,426,782个顺式甲基化数量性状位点（cis-meQTLs）、5,078个顺式表达数量性状甲基化（cis-eQTMs）和9,157个表达数量性状位点（eQTLs）。其中36.39%的顺式eQTM基因与遗传变异无关，表明许多与基因表达变异相关的SMPs不依赖于SNPs。

在与产量和纤维品质性状相关的1,715个表观遗传位点中，仅有36个（2.10%）与全基因组关联研究（GWAS）位点一致。多组学调控网络的构建揭示了可能43个参与纤维发育的cis-eQTMs基因，这些基因无法仅靠GWAS分析鉴定。在这些基因中，通过基因编辑成功验证了编码CBL互作蛋白激酶10的基因在纤维长度调控中的作用。总之，本研究结果表明，DNA甲基化数据可以作为育种的额外资源，并为加强和加速作物改良进程提供机会。

研究方法

植物材料：选取207份棉花材料，在开花后20天（20-DPA）时采集纤维样本，这一时期是纤维发育的关键阶段，特别是在次生细胞壁加厚过程中。

组学测序：采集20-DPA纤维样本进行WGBS、RNA-seq和WGS。WGBS计算甲基化水平和SMPs注释。

eQTM分析：使用fastQTL软件进行eQTM分析，研究DNA甲基化变异与基因表达的关系。

cis-meQTLs分析：使用fastQTL软件进行cis-meQTLs分析，研究遗传变异与DNA甲基化的关系。

eQTLs分析：通过连锁不平衡（LD）分析减少meQTLs的冗余性，并将关联SNPs合并为一个meQTL。最终鉴定出的meQTLs进一步分类为cis-eQTLs和trans-eQTLs。

]EWAS分析：鉴定与纤维性状相关的表观遗传位点。

植物材料、载体构建与遗传转化：使用CRISPR-Cas9技术进行基因编辑，验证关键基因的功能。

表型预测：使用G2Pdeep和GBLUP方法进行表型预测，评估模型性能。

基于深度学习的功能性CG甲基化位点预测：开发DeepFDML模型，基于DNA序列预测功能性CG甲基化位点。

结果图形

（1）DNA甲基化变异图谱的构建与表征

研究人员对207份核心种质陆地棉群体（CUCP1）的20天开花后（20-DPA）的纤维进行了全基因组亚硫酸盐测序（WGBS）和转录组测序（RNA-seq），以构建单碱基DNA甲基化变异图谱。所有样本的甲基化组数据均比对到棉花参考基因组TM-1 v2.1，平均比对率为74.90%±3.55%。每个测序的甲基化组平均覆盖深度＞15X。经过严格的数据处理和质控，研究人员在CG、CHG和CHH三种甲基化背景下分别定量了62.32M、66.06M和4.3301M甲基化胞嘧啶。

棉花基因组的DNA甲基化水平较高，尤其是在异染色质区域。全基因组范围内，CG、CHG和CHH位点的DNA甲基化水平分别为72%、55%和11%。研究还发现，DNA甲基化变异在基因内区域显著富集，而在转座元件（TEs）中则相对缺失。此外，研究人员定义了单甲基化多态性（SMPs），并基于SMPs进行了系统发育分析，将207份材料分为四个分支。这些结果表明，DNA甲基化是棉花基因组中一个重要的变异来源，且在基因内区域具有较高的变异频率。

图1：自然群体中DNA甲基化变异的广泛模式。 a.多组学关联分析的工作流程。 b.不同材料间DNA甲基化多样性的基因组区域示例。 c.三维图展示不同材料的DNA甲基化多样性。 d.三种不同DNA甲基化背景之间的相关性。 e.条形图显示MAF大于0.5的SMP的数量和比例。 f.密度图展示CG、CHG、CHH-SMPs和SNPs的MAF分布。 g-h. 盒须图展示不同基因组特征中CG位点的SMP（g）和SNP（h）的MAF分布。 i.   不同DNA甲基化背景之间的连锁不平衡（LD）衰减比较（垂直轴：LD水平；水平轴：成对距离）。

（2）基因富集区域的遗传变异对甲基化组的显著影响

研究人员进一步分析了遗传变异对DNA甲基化的作用。通过全基因组随机抽样，研究人员在CG、CHG和CHH三种甲基化背景下分别鉴定出119,685、37,831和24,683个甲基化数量性状位点（meQTLs）。其中，顺式meQTLs（cis-meQTLs）表现出更强的显著性。研究人员对所有SMPs进行了顺式meQTL分析，共鉴定出5,426,782个顺式meQTLs，包括940,794个CG-cis-meQTLs、883,280个CHG-cis-meQTLs和3,602,708个CHH-cis-meQTLs。这些顺式meQTLs在基因组中的分布不均匀，特别是在染色体末端密度较高。此外，顺式meQTLs在基因内区域显著富集，而在TEs中则显著缺失。这些结果表明，基因富集区域的遗传变异对DNA甲基化具有重要影响。

图2：三种DNA甲基化背景的遗传基础。 a.SMP及其相关SNPs的基因组分布。x轴表示显著SNPs的基因组位置，y轴表示相应SMPs的基因组位置。选择50,000个CG、CHG和CHH的SMP进行全基因组meQTL分析。 b.盒须图显示顺式（cis-）和反式（trans-）meQTL的-log10(P) 分布。显示中位数和四分位间距（IQR）。 c.UpsetR图展示不同DNA甲基化背景下共有顺式meQTLs比例。 d.DMR与显著SNP之间的距离。 e.顺式meQTL在全基因组范围内的分布。 f.顺式meQTL在不同基因组特征中的富集和缺失情况。

（3）SMPs在表达调控中的作用

研究人员进一步探讨了DNA甲基化与基因表达之间的关系。通过eQTM分析，研究人员共鉴定出5078个顺式eQTMs，涉及3505个蛋白编码基因（PCGs）和1573个长链非编码RNA（lncRNAs）。这些eQTM基因在长链脂肪酸代谢、毛状体分支和葡萄糖稳态等过程中富集，可能与纤维发育相关。此外，研究人员发现，CG甲基化在基因表达调控中起着更为重要的作用。在顺式eQTMs中，90%的eQTM基因与CG甲基化相关。这些结果表明，DNA甲基化在基因表达调控中具有重要作用，尤其是CG甲基化。

图3：受DNA甲基化调控的基因表达变异。 a.eQTL分析的工作流程。 b.在PCGs和lncRNAs中鉴定出的顺式eQTMs数量。 c.受DNA甲基化调控的PCGs和lncRNAs比例。 d.重叠和特定的顺式eQTM基因的UpsetR图。右侧：PCGs；左侧：lncRNAs。 e.靠前SMPs与相关转录起始位点的距离。左侧：PCGs；右侧：lncRNAs。 f.高置信度eSNP-表达关联散点图。每个点代表一个检测到的eQTL。 g.顺式和反式eQTLs数量的饼图。 h.eQTM和meQTL分析中鉴定出的eQTM基因表征。这些位点被分为三组。遗传和顺式表观遗传调控（类型I）、遗传和反式表观遗传调控（类型II）以及仅表观遗传调控（类型III）。

（4）表观基因组关联研究（EWAS）揭示了大量与农艺性状相关的优异表观等位基因

研究人员利用SMPs进行了表观基因组关联研究（EWAS），鉴定出1715个与纤维产量和品质性状相关的表观遗传位点。其中，1010个位点与产量相关性状相关，705个位点与纤维品质性状相关。这些表观遗传位点中，只有16个（0.93%）与GWAS位点重叠。例如，研究人员在A11染色体上发现了一个与衣分（lint percentage, LP）相关的表观遗传位点，该位点位于一个编码核孔复合体相互作用组分（Nup93）的基因启动子区域。不同表观等位基因对应的LP值存在显著差异。这些结果表明，DNA甲基化为农艺性状提供了额外的调控层，并且大多数表观遗传位点与遗传变异独立存在。

图4：EWAS位点分布及其对农艺性状的累积效应。 a.与农艺性状相关的EWAS位点分布。纤维产量性状包括衣分（LP）；纤维品质：纤维长度（FL）、强度（FS）、伸长率（FE）、细度（FM）和均匀度（FU）。与每个性状相关的位点在染色体图上以黑色垂直线表示。 b.UpsetR图展示CG-EWAS、CHG-EWAS和CHH-EWAS之间的重叠。 c.与小于2-kb和5-kb区域内有侧翼基因的EWAS位点比例。 d.EWAS分析LP性状的曼哈顿图。红色箭头表示在A11染色体上的信号。 e.放大图显示靠前SMP代表A11染色体上LP的EWAS位点，且信号坐标位于同一甲基化不平衡块中。 f.不同表观等位基因的LP值，用于e中所示的位点。 g.根据CG-、CHG-和CHH-EWAS位点以及GWAS位点显示自然种群中单倍型分布的热图。优异等位基因用红色表示。每一列代表一个材料，每一行指基因组中的一个位点。 h.对在EWAS和GWAS中鉴定出的位点表征描述。这些位点被分为四组。仅表观遗传调控（类型I）、仅遗传调控（类型II）、遗传/顺式表观遗传调控（类型III）以及遗传/反式表观遗传调控（类型IV）。

（5）通过多组学关联分析鉴定纤维相关基因

研究人员通过整合GWAS、eQTL和EWAS结果，构建了纤维性状的基因调控网络。研究人员发现，51个GWAS位点与376个eQTL在同一个连锁不平衡（LD）块内共定位。基于LD块分析，研究人员构建了一个包含397个基因的基因调控网络。此外，研究人员还构建了一个基于EWAS和eQTM的表观遗传调控网络。通过比较这两个网络，研究人员发现它们之间只有四个共有基因。这些结果表明，纤维性状的调控机制非常复杂，涉及遗传和表观遗传的多重调控。

图5：与纤维发育相关的遗传和表观遗传调控网络。 a.功能性基因调控网络（GRN）构建的分析流程。分别进行eQTM和eQTL分析，以获得EWAS和GWAS位点中的因果位点。将同一连锁不平衡（LD）块内的位点合并为一个主要SNP，LD块内的eGenes被聚类成一个GRN。对EWAS位点也进行相同步骤。 b.调控棉花纤维性状的基因网络。右侧：通过整合GWAS和eQTL构建的遗传变异依赖网络；左侧：通过整合EWAS和eQTM构建的表观遗传调控网络。 c.通过共定位分析确定的候选基因热图。 d.不同表观等位基因中CIPK10的表达水平和LP值。 e.基因编辑（CRISPR敲除，CR-KO）对调控纤维性状的eQTM基因GhCIPK10的性能的图像。 f.两条CIPK10CR-KO品系的纤维长度。

（6）基于DNA序列使用DeepFDML预测功能性CG甲基化位点

研究人员开发了一个名为DeepFDML的深度学习模型，用于预测功能性CG甲基化位点。该模型基于2336个与基因表达变异相关的CG位点进行训练。研究人员使用one-hot encoding将每个CG位点的侧翼序列转换为模型input。通过五折交叉验证，该模型的准确率达到了0.65。研究人员进一步构建了一个更复杂的DeepFDML模型，该模型采用了预训练的Enformer模型作为其主干网络。该模型的接收者操作特征曲线（ROC）和精确召回曲线（PRC）分别达到了0.82和0.78，显著优于卷积模型。这些结果表明，基于DNA序列的深度学习模型可以有效预测功能性CG甲基化位点。

图6：用于功能性CG位点预测的卷积神经网络。a.所提出的深度学习方法的流程示意图。它主要包括四个部分：input序列、主干网络、特征选择和output层。每个input是一个以CG位点为中心的8192bp的DNA序列，经过one-hotencoding处理。主干网络采用预训练的Enformer模型。在特征选择阶段，利用了中间八个位置的特征。output层是一个作为二元分类器的全连接层。b.在整个数据集上测量的接收者操作特征（ROC）曲线。c.在整个数据集上测量的精确率-召回率曲线（PRC）。

讨论和启示

本研究证明了DNA甲基化数据在作物育种中的重要应用价值。通过整合表观组（DNA甲基化）、转录组、基因组的多组学数据，研究人员不仅发现了大量与纤维性状相关的表观遗传位点，还通过基因编辑技术验证了关键基因的功能。这些发现为棉花育种提供了新的策略，有助于加速作物改良进程。

此外，研究人员开发的DeepFDML模型为预测功能性CG甲基化位点提供了一种新的方法，这将有助于未来在其他物种中进行类似的研究。本研究强调了DNA甲基化在基因表达调控和性状改良中的重要作用，为表观遗传学在作物改良中的应用提供了新的思路。

参考文献：

Zhao, T., Guan, X., Hu, Y.et al.Population-wideDNA methylation polymorphisms at single-nucleotide resolution in 207 cottonaccessions reveal epigenomic contributions to complex traits.Cell Res34,859–872 (2024). Doi:10.1038/s41422-024-01027-x