2024年10月18日，浙江大学棉花精准育种团队的方磊教授和张天真教授，与中国农科院生物技术所的谷晓峰研究员及阿里巴巴达摩院（湖畔实验室）顾斐博士团队，在国际学术期刊Cell Research上发表了 “Population-wide DNA Methylation Polymorphisms at Single-nucleotide Resolution in 207 Cotton Accessions Reveal Epigenomic Contributions to Complex Traits”研究论文。本研究利用遗传学、生物大数据和人工智能技术，构建了迄今为止作物中最大的全基因组DNA甲基化图谱，首次绘制出纤维发育相关的独立于遗传因素的表观调控网络。通过构建深度学习模型DeepFDML，实现了功能性表观遗传修饰位点的全基因组智能预测。
棉花是全球最大的天然纺织纤维来源，同时也是研究细胞分化、伸长、细胞壁发育调控和作物多倍体化的重要模式植物。解析棉花复杂农艺性状的遗传和生物过程的分子机制，对棉花的分子育种具有重要意义。全基因组关联研究（GWAS）常用于解析作物的复杂农艺性状，但无法揭示表观遗传层面的群体多样性与性状变异的相关性。
DNA甲基化是表观遗传的重要组成部分，可以通过自然进化和人工选择传递，但其群体多样性和在育种中的作用仍然未知。为了探究DNA甲基化对农艺性状的贡献，本研究利用207个棉花种质资源的全基因组亚硫酸盐测序数据、基因组测序数据、转录组测序数据，结合农艺性状的表型数据，构建了群体尺度的棉花全基因组DNA甲基化图谱，并将经典群体遗传学分析框架扩展到了表观遗传学标记中，开展了亿级多组学的关联分析。
棉花的基因组较大（~ 2.5 G），通过15倍覆盖度的亚硫酸盐测序（Whole-Genome Bisulfite Sequencing，WGBS），鉴定了2.87 亿个单甲基化多态性（Single Methylation Polymorphism, SMP）位点，其数量是单核苷酸多态性（SNP）位点的100倍。不同甲基化类型中，全基因组CG与CHG的DNA甲基化高度相关（PCC = 0.95），而与CHH甲基化的相关性较低（PCC = 0.22）。不同基因组区域的DNA甲基化多态性受到不同选择压力。转座子通常高度甲基化，其CG-SMP的次要等位基因频率（Minor Allele Frequency , MAF）只有蛋白质编码基因MAF值的一半，表明甲基化在抑制转座子中发挥重要作用。为表征相邻DNA甲基化位点间的关系，研究将连锁不平衡（Linkage Disequilibrium, LD）概念扩展到DNA甲基化标记，并称之为甲基化不平衡（Methylation Disequilibrium , MD）。MD衰减到其最大值一半的平均距离约为50 bp，这与人类和拟南芥的先前报道的DNA甲基化的MD一致。MD的衰减速度显著快于LD。此外，CHH的MD低于CHG和CG的MD。这表明CG甲基化可能相比非CG甲基化更易于遗传。
通过多组学联合分析，识别出5,426,782个顺式甲基化数量性状基因座（cis-meQTL）、5,078个顺式表达数量性状甲基化（cis-eQTM）以及9,157个表达数量性状基因座（eQTL）。
meQTL（Methylation Quantitative Trait Loci）是指通过将DNA甲基化水平作为分子表型，与遗传变异进行关联分析而获得的QTL。研究鉴定出5,426,782个meQTL，其中包括940,794个CG、883,280个CHG和3,602,708个CHH Cis-meQTL，分别占CG位点总数的5.82%、CHG位点的2.64%和CHH位点的1.52%。尽管上述结果依赖于计算参数，但在相同参数下，CG类型的DNA甲基化与遗传变异的关系更为紧密。
eQTM（Expression Quantitative Trait Methylation）指的是与基因表达相关的甲基化数量性状位点。研究鉴定出5,078个cis-eQTM，包括3,505个基因-eQTM和1,573个长链非编码RNA-eQTM，占表达的蛋白质编码基因的5.69%和lncRNA的29%。通过分析发现，一个基因可同时与不同类型的DNA甲基化相关。从eQTM的数量来说，与CG甲基化相关的cis-eQTM基因在所有cis-eQTM中占较大比例（蛋白质编码中为30.85%，lncRNA中为60.24%）。从基因个数上说，大多数eQTM基因与CG甲基化有关，分别占基因和lncRNA 的91%和96%。这表明CG甲基化在基因调控中比其他类型甲基化更关键。
在群体水平上，90%的cis-eQTM基因的SMP倾向于位于PCG和lncRNA的上游，暗示其可能通过改变启动子的甲基化状态影响基因表达。有趣的是，36.39%的cis-eQTM基因表达与遗传变异无关，这表明大量与基因表达差异相关的SMP独立于SNP。研究讨论了表观与遗传共同或独立调控基因表达的三种模式。
由于大量的甲基化位点并不与SNP连锁，因此表观关联分析能够发现SNP关联分析无法检测到的位点，可以作为GWAS的重要补充。通过表观全基因组关联分析（EWAS），鉴定了1,715个与棉花农艺性状相关的表观遗传位点，将具有育种价值的遗传标记拓展到表观遗传标记。其中有36个（2.10%）与全基因组关联研究（GWAS）位点重合。多组学调控网络的构建揭示了43个潜在参与纤维发育的eQTM基因，为棉花功能基因组提供了重要的候选基因列表。
解析调控元件的功能影响是功能基因组学研究中的关键挑战，对推动下一代作物育种策略至关重要。深度学习模型已被应用于揭示遗传元件的功能模式。然而，目前尚未开发出用于预测功能性表观遗传修饰位点的方法。本研究发现，影响基因表达的eQTM位点通常位于启动子上游，其DNA甲基化水平和所在的DNA序列位置都可能具有重要作用。据此我们研发了深度学习模型DeepFDML，基于DNA甲基化位点附近的序列，预测影响基因表达的甲基化位点。
本研究系统分析了群体水平DNA甲基化多态性对基因表达的调控模式和对农艺性状改良的贡献，揭示了表观与遗传协同调控农艺性状的复杂网络，为实现表观遗传设计育种提供了重要的育种资源与指导。
浙江大学棉花精准育种团队赵汀特聘研究员、关雪莹研究员、胡艳教授、硕士研究生张子迁和阿里巴巴算法专家杨晗硕士为论文共同第一作者，浙江大学方磊教授、中国农科院生物技术所谷晓峰研究员、阿里巴巴达摩院（湖畔实验室）顾斐博士、浙江大学张天真教授为论文共同通讯作者。本研究工作获得国家科技重大专项、国家自然科学基金、浙江大学-阿里巴巴AZFT等项目资助。
图：群体单碱基DNA甲基化多态性的
精准鉴定、调控网络解析与预测模型构建
