Cell Research | 207份棉花的群体DNA甲基化多态性揭示了表观基因组对复杂性状的贡献

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甲基化图谱，将经典群体遗传学框架扩展到表观遗传学领域，在作物基因组中进行了单核苷酸水平的EWAS分析，首次绘制出纤维发育相关的独立于遗传因素的表观调控网络。通过构建深度学习模型DeepFDML，实现了功能性表观遗传修饰位点的全基因组智能预测，为棉花作物改良提供新途径。

表型变异是遗传、表观遗传变异与环境动态交互作用的结果。虽然GWAS研究已显著推进了对遗传变异的理解，但表观遗传对作物性状多样性的影响仍未被充分探索。DNA甲基化是最广泛研究的表观遗传标记之一，其通过调控基因表达、抑制转座元件活性和维持基因组稳定性，对植物重要农艺性状（如开花时间、产量、抗逆性）产生重要影响。

植物的DNA甲基化可发生在CG、CHG和CHH背景下，且不同背景的甲基化通过特定机制维持。在棉花纤维发育中，DNA甲基化的动态变化为研究表观遗传变异与纤维性状关联提供了机会。本研究整合了207个棉花品种开花后20天纤维的甲基化组、转录组和基因组数据，全面分析了DNA甲基化与纤维品质和产量的关系，鉴定了关键基因和调控位点。这些发现为棉花育种提供了宝贵的基因组和表观基因组资源，有助于通过新型选择策略改良作物性状。

1.  DNA甲基化变异图谱的构建与特征分析

本研究利用207份陆地棉核心种质群体（CUCP1）进行多组学整合研究。实验材料均在2021年于中国杭州种植，研究人员在次生细胞壁增厚阶段采集了开花后20天的棉花纤维，获得了高质量的基因组、甲基化组及转录组数据，并构建了表达数量性状位点（eQTL）和表达数量性状甲基化（eQTM）图谱（图1a）。

研究发现，陆地棉基因组高度甲基化，尤其是在异染色质区域。CG、CHG和CHH背景下的全基因组甲基化水平分别为72%、55%和11%。CG甲基化比例与前人研究一致，显著高于橙子（CG: 41%），但低于大麦（CG: 94.7%）。此外，207个样本中胞嘧啶甲基化水平全基因组范围内存在变异（CG-IQR=8.08%；CHG-IQR=6.75%；CHH-IQR=1.37%）（图1b, c）。全基因组CG DNA甲基化水平与CHG DNA甲基化水平呈显著正相关，而与CHH DNA甲基化水平的相关性较低（图1d）。

研究发现，棉花基因组中常见的SMPs（single methylation polymorphisms）数量（CG为1615万，CHG为3341万，CHH为2.377亿）远高于SNPs（128万）（图1e）。CHH-SMP的MAF值（0.22）显著高于CHG-SMP（0.11）、CG-SMP（0.05）和SNP（0.14）（图1f）。不同基因组特征中，SMP的MAF值差异显著（图1g）。CG-SMP在转座元件（TEs）中的MAF值仅为蛋白编码基因（PCGs，包括外显子和启动子）的一半，而SNP的MAF值在TEs和PCGs中相似。尽管TEs通常高度甲基化，CG-SMP的MAF值在TEs中显著较低。常见的CG-SMP在基因内部区域富集程度显著高于其他SMP类型（CG-SMP：27.53%；CHG-SMP：9.78%；CHH-SMP：10.49%），表明基因内转录区域的CG甲基化更容易发生变异。

为了表征邻近DNA甲基化位点之间的关系，本研究将连锁不平衡（LD）的概念应用于DNA甲基化，即甲基化不平衡（MD）。MD衰减至其最大值一半的平均距离约为50 bp（图1i），衰减速度显著快于LD（300 kb以上）。此外，CHH的MD低于CHG和CG（图1i），表明CG和CHG的甲基化在有丝分裂和减数分裂中可能更容易维持。总体而言，DNA甲基化是基因内部区域重要的变异来源。

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

研究进一步分析了遗传变异对DNA甲基化的影响。随机抽取50,000个CG-SMPs、CHG-SMPs和CHH-SMPs（分别占其类型的0.31%、0.15%和0.021%），重新评估其在顺式（cis）和反式（trans）meQTL中的作用，并定义SNP与相关SMP之间距离在1 Mb以内的为顺式meQTL。结果鉴定出119,685个CG-meQTLs、37,831个CHG-meQTLs和24,683个CHH-meQTLs。尽管通过meQTL分析鉴定出大量反式meQTL（图2a），但顺式meQTL相比反式meQTL表现出更高的显著性水平（图2b）。

为减少假阳性并优化计算，仅选择顺式meQTL进行进一步分析。使用fastQTL分析了2,873亿SMP，鉴定出5,426,782个顺式meQTL，包括940,794个CG-cis-meQTLs、883,280个CHG-cis-meQTLs和3,602,708个CHH-cis-meQTLs。仅有少部分甲基化位点（CG 5.82%、CHG 2.64%、CHH 1.52%）参与了顺式meQTL，一部分顺式meQTL在CG、CHG和CHH中显示出共定位现象（图2c）。此外，CG-cis-meQTL中SNP与SMP之间的距离显著短于CHG和CHH背景（图2d）。顺式meQTL在基因组中分布不均匀，染色体末端密度较高（图2e）。顺式meQTL显著富集于基因内区域（P < 2.2×10^-16），但在转座元件（TEs）中显著缺乏（P <2.2×10^-16）（图2f）。

3. SMPs参与表达调控

研究发现，顺式meQTL在自然群体的蛋白编码基因（PCGs）中富集（图2e, f），因此探索DNA甲基化与基因表达的关系尤为重要。本研究通过对相同组织的转录组进行分析，研究了DNA甲基化对局部基因表达的影响（图3a）。使用TM-1 v2.0参考基因组注释的71,994个PCGs和21,181个长链非编码RNA（lncRNAs），提取了在超过5%群体中表达的41,632个PCGs和5469个lncRNAs用于eQTL和eQTM分析。通过fastQTL软件鉴定出5078个顺式eQTMs，包括3505个PCG-eQTMs和1573个lncRNA-eQTMs（图3b），这些eQTMs共映射到2619个基因，占表达PCGs的5.69%和lncRNAs的29%（图3c）。GO功能分析表明，顺式eQTM基因显著富集于长链脂肪酸代谢、毛状体分枝和葡萄糖稳态等过程，与纤维发育相关。

此外，不同甲基化背景下的顺式eQTM基因存在共定位现象，例如同时与CG、CHG和CHH甲基化相关的顺式eQTM基因占PCGs的30.85%和lncRNAs的60.24%（图3d）。CG甲基化与大多数eQTM基因相关，占PCG和lncRNA eQTM的91%和96%，表明CG甲基化在基因调控中起更重要作用（图3d）。进一步分析发现，90%的顺式eQTMs位于PCGs和lncRNAs的上游（图3e）。此外，位于近端启动子的CG-eQTMs和CHG-eQTMs的甲基化水平与基因表达呈负相关，而远端区域和基因体内的eQTMs影响较小。

通过EMMAX模型结合SNP和表达谱进行eQTL分析，共检测到9157个eQTL，包括5921个eSNP和7398个eGene（PCGs: 5197, lncRNAs: 1014）（图3f）。根据基因位置（1 Mb阈值），这些eQTL分为926个顺式eQTL和8231个反式eQTL（图3g）。

将eQTM基因调控模式分为三类：遗传/顺式表观遗传调控（I型）、遗传/反式表观遗传调控（II型）和仅表观遗传调控（III型）（图3h）。II型eQTM基因占比小于20%，III型eQTM基因占比为33.63%–38.14%，表明DNA甲基化在基因表达调控中具有积极作用。共调控基因显著富集于含氮有机化合物、硫化合物和乙酰辅酶A生物合成过程。

4. 表观基因组广泛关联研究（EWAS）揭示了大量优良的表观等位基因

研究发现，大多数SMPs与遗传变异无关，这与先前的研究一致，表明DNA甲基化变异独立于遗传变异，因此表观遗传学的关联可能无法通过基于SNP的标记捕捉到。研究团队使用全基因组范围内常见的SMPs（MAF≥0.05）针对9个性状进行EWAS分析，共识别出848个CG-EWAS位点、467个CHG-EWAS位点和400个CHH-EWAS位点，总计1715个位点（图4a），其中1010个位点与产量相关，705个位点与纤维品质相关。在不同的甲基化背景下，大多数EWAS位点是独立的，只有22个位点在至少两个甲基化背景下共享（图4b）。

约27.67%的CG-EWAS位点、19.92%的CHG-EWAS位点和16.19%的CHH-EWAS位点位于蛋白质编码基因或lncRNA基因的邻近区域（2kb内）（图4c）。编码核孔复合物相互作用的Nup93基因的启动子区域出现了一个与产量性状（绒花百分比，LP）相关的EWAS信号（图4d, e）。不同的表观等位基因对应不同的LP值（图4f）。

为了分析遗传和表观遗传变异在性状变异中的关系，研究将EWAS位点与GWAS位点结合（图4g）。发现仅有0.93%的EWAS位点位于GWAS位点附近（<20kb），这表明DNA甲基化更多独立于遗传变异位点调控农艺性状。研究还分析了与GWAS位点不重合的EWAS位点，发现992个位点具有转位meQTL效应，其中20个位点与GWAS位点相关（图4h）。最后，通过分析携带多个优良表观等位基因组合的样本，发现携带更多优良等位基因的样本表现出更好的性状。这表明，结合SNPs和SMPs可以提高农艺性状预测的准确性，特别是与纤维产量和品质相关的性状。

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

通过多组学关联分析，研究人员鉴定了与纤维性状相关的基因，包括187个GWAS位点、9157个eQTL、1715个EWAS位点、5078个cis-eQTMs和5426782个cis-meQTLs。随后整合了GWAS位点和eQTLs，基于LD块构建了基因表达调控网络（GRN）（图5a）。发现51个GWAS位点与376个eQTL在同一LD块内共定位（r² > 0.1），包含了397个基因之间的634个连接。

在这些基因网络中，77个（19.4%）eQTL基因也是eQTM基因，说明基因表达受DNA甲基化和遗传变异的共同调控。与四种纤维性状（纤维产量（LP）、强度（FS）、长度（FL）和微纤度（FM））相关的基因网络包含了多种已知与纤维伸长相关的基因，如Expansion A4、纤维素合成类基因（CSL）等。

研究人员通过整合EWAS位点和eQTMs来构建了表观遗传学调控网络（即表观遗传GRN）。此外，研究人员还通过47个与EWAS位点共定位的eQTMs构建了另一个表观遗传GRN。对比这两个网络发现，只有四个基因在两者之间共享，它们分别编码胰蛋白酶、RIBOSOMAL PROTEIN EL8Y、GH_A06G1022和醛脱氢酶，这一结果揭示了纤维性状调控机制的复杂性。

另外，研究人员发现与LP相关的EWAS位点（A03:4217197）位于CIPK10基因的启动子区域，该基因编码一种CBL-交互蛋白激酶，是一个潜在的纤维发育基因。该基因被鉴定为eQTM基因，且其CG-SMP（A03:4217260）的DNA甲基化状态与CIPK10的表达和LP值相关。通过CRISPR/Cas9敲除CIPK10基因，发现纤维长度显著缩短，进一步验证了CIPK10在纤维发育中的重要作用（图5e, f）。

6. 基于DNA序列预测功能性CG甲基化的DeepFDML模型

揭示调控元件的功能性影响是功能基因组学研究中的一项重要挑战，对于推动下一代作物育种策略至关重要。深度学习模型已被应用于通过整合基因组序列和分子特征（如非编码区域转录和启动子中的顺式元件）来揭示遗传元件中的功能模式。然而，针对预测功能性表观修饰位点的方法尚未得到开发。

在本文章中，研究人员开发出了一个深度学习模型DeepFDML，用于预测与基因表达变化相关的功能性SMPs。该模型基于2336个与2423个CG-eQTMs相关的功能性CG位点进行训练，并随机选取了另一组CG-SMPs作为阴性对照。每个CG-SMP位点的邻近序列通过独热编码进行转换（图6a）。

为了评估DNA甲基化对基因表达的影响，研究者首先构建了一个简单的卷积模型，使用五折交叉验证，在受试者工作特征曲线（ROC）和精确率-召回率曲线（PRC）中的准确度均达到了0.65（图6b）。随后，采用更复杂的DeepFDML模型，包含七个卷积池化块和11个transformer编码层，使模型准确度显著提高，ROC为0.82，PRC为0.78（图6b, c）。研究表明，DeepFDML能够通过DNA序列模式有效预测功能性SMPs。