两种诱变处理下紫花苜蓿全基因组DNA甲基化图谱分析

申晓慧

(宜春学院， 生命科学与资源环境学院， 江西 宜春 336000)

紫花苜蓿(MedicagosativaL.)具有产量高、适口性好、利用年限长、营养价值丰富和广适应性等特点，随着我国畜牧业的快速发展，家畜数量急剧上升，对饲草数量和品质要求越来越高[1]，为此，紫花苜蓿在畜牧业发展中起着举足轻重的作用。然而我国一些极端地区如干旱、高寒、盐碱地等特殊的环境条件给紫花苜蓿丰产性带来了严重影响，这就导致生产中迫切需要更新紫花苜蓿种质资源来应对这种特殊环境条件。因此筛选、培育抗性品种成为育种工作中的关键。实践证明，利用人工诱变遗传变异是丰富作物种质资源，创制物种新材料及新品种选育的重要手段之一[2]。然而，诱变后代能否保证优良性状能够稳定遗传对新品种选育至关重要。表观遗传可以通过调控基因表达进而影响植物生长发育，而DNA甲基化作为表观遗传学研究的一种，其在调节植物生长发育、植物应对生物或非生物胁迫环境以及维持基因组稳定等方面发挥着至关重要的作用[3-5]。

目前，随着基因二代测序技术的不断向前发展，带动着全基因组的DNA甲基化的飞快进步，其研究成果在多学科、多领域中都发挥着重要的积极作用。如薛倩[6]利用BS-seq技术构建了京海黄鸡腿肌全基因组DNA甲基化图谱，揭示了其DNA甲基化特征。苏畅[7]利用重亚硫酸盐测序技术构建了白桦基因组甲基化图谱，综合分析了DNA甲基化对于基因表达的作用。赵倩等[8]研究表明，水稻在高剂量重离子辐射下，DNA去甲基化的多态性水平明显高于DNA甲基化的多态性水平。还有一些在小麦[9]、中草药[10]以及绿化植物长春花[11]等植物方面也都有诸多研究。在我国诱变技术在紫花苜蓿中已有广泛应用，如韩微波等[12]通过零磁空间处理肇东紫花苜蓿种子，从中选育出了新品种‘农菁10号’。然而，在众多研究中通过诱变手段导致紫花苜蓿全基因组DNA甲基化水平的变化研究却鲜有报道。基于此，作者在已有诱变对苜蓿生长、生理变化的影响[13-14]的研究基础上，利用BS-seq技术进一步深入分析诱变处理前后紫花苜蓿全基因组甲基化图谱变化情况，为紫花苜蓿突变体筛选及诱变条件下苜蓿表观遗传变化的分子机制研究提供科学依据。

1 材料与方法

1.1 试验材料

选取紫花苜蓿‘龙牧806’，该品种由黑龙江省畜牧研究所提供。

1.2 试验方法

1.2.1试验处理 甲基磺酸乙酯(EMS)处理[15]：首先将成熟饱满的苜蓿种子用浓硫酸处理5 min，之后用蒸馏水反复冲洗。然后在浓度为100 mmol·L-1,pH值 7.0的4℃磷酸缓冲液中浸泡12 h。最后用磷酸缓冲液配制成浓度为0.4%的EMS溶液，将浸泡过的种子置于室温黑暗条件下处理15 h，而后用蒸馏水反复冲洗，去除残留种子表面的EMS溶液。每处理保留100粒种子，3次重复。

60CO-γ射线处理[15]：将成熟饱满的苜蓿种子进行150 Gy剂量的60CO-γ射线辐照处理，处理周期为7 d。以未处理种子作为对照，每处理保留100粒种子，3次试验重复。

1.2.2DNA提取、建库及测序 将上述各重复处理的种子分别进行混合，每100粒作为一个样本，每个样品3个生物学重复，进行DNA提取，DNA的提取按照张延召等[16]的方法进行。DNA样品经武汉康测科技有限公司质检合格后进行重亚硫酸盐(BS-seq)文库构建，主要包括基因组DNA提取(CTAB法)，片段化，利用EZ DNA Methylation-Gold kit进行Bisulfite处理及基于单链的DNA文库构建等4个步骤。文库质量检测后通过Illumina XTen (Illumina，San Diego，CA，USA)测序仪进行双末端测序，具体方法参照文献[17]。试验选取样本材料分别为对照组(A34-1A)，EMS处理组(A34-2A)及60CO-γ射线处理组(A34-3A)。

表1 单链DNA浓度展示表Table 1 Single-stranded DNA concentration display table

1.2.3测序数据的过滤、比对分析 测序结果经过识别碱基后转化为原始测序序列(Raw reads)，然后除去Raw reads中含接头序列、低质量、含N比例>10%的Reads后得到高质量序列(Clean Reads)。使用Bismark(v0.17.0)软件对Clean reads进行比对。通过将reads和参考基因组分别进行C>T和G>A的转换之后，分别进行比对，筛选出比对结果最好的组合方式[18]。参考基因组选用Medicagotruncatula，源自ftp://ftp.ensemblgenomes.org/pub/release-33/plants/fasta/medicago_truncatula/dna/Medicago_truncatula.MedtrA17_4.0.dna.toplevel.fa.gz。

1.2.4苜蓿组织差异甲基化分析 使用swDMR软件(version 1.0.0)检测DMR[19]。swDMR采用滑动窗口(sliding window size:1 000 bp;step length:100 bp)检测DMR。大于2个样本时，采用ANOVA方法。

1.2.5DMR(Differentially methylated region)基因KEGG通路富集分析 通路富集分析以KEGG Pathway为单位[20]，应用超几何检验，找出与整个苜蓿基因组背景相比，在DMR相关基因中显著富集的通路。采用KOBAS[21]对DMR所在基因列表进行基因通路富集分析，筛选出通路富集P-value(或校正后P-value)小于0.05以及查询基因在此通路中的基因数目大于2的通路。

2 结果与分析

2.1 紫花苜蓿全基因组DNA甲基化测序数据分析

使用Illumina XTen测序平台对供试样品进行测序，共得到58.11 G的Clean reads (表2)。使用软件SOAPnuke(version 1.6.0)将A34-1A，A34-2A 及A34-2A样品得到的Clean reads比对到苜蓿基因组序列上，比对率分别为 90.40%，88.72%和90.60%。

表2 不同处理苜蓿基因测序数据统计及比对到参考基因组结果Table 2 Statistical results of sequencing data and clean reads map to reference genome results

2.2 紫花苜蓿全基因组DNA甲基化C特征分析

C位点的DNA甲基化主要有3种类型：mCG，mCHG和mCHH。这 3种类型的数量及构成比例反映了特定物种的全基因组甲基化的特征。本研究发现，A34_1A，A34_2A，A34_3A中分别有32.00%，35.79%和30.82%的C位点发生了甲基化，且DNA甲基化变化主要有mCG，mCHH和mCHG 3种类型，其中A34_2A的各位点甲基化比例高于A34_1A，A34_3A的各位点甲基化比例低于A34_1A (表3)。对甲基化C碱基中CG，CHG和CHH的分布比例结果表明，苜蓿诱变前后其mCHH位点所占比例最高，变化趋势为先升后降，而mCG和mCHG位点出现频率表现为先降后升的变化趋势(图1)。研究结果表明，不同诱变处理紫花苜蓿后，C位点发生甲基化呈现出不同变化趋势，不同类型的mC模式均发生了变化，说明DNA甲基化水平变化在一定范围内是保持相对稳定的。

表3 苜蓿不同分布类型甲基化C的数量及比例Table 3 The number and proportion methylated C in alfalfa

图1 甲基化C碱基中CG，CHG和CHH的分布比例Fig.1 The distribution of mCG,mCHG and mCHH in all methylcytosine

2.1.2紫花苜蓿全基因组DNA甲基化水平分析 由供试苜蓿全基因组C位点DNA甲基化水平分析结果(表4)可知，EMS处理后苜蓿DNA甲基化水平由15.67% 提高到17.43%，达到差异显著水平(P<0.05)，而60CO-γ射线辐射处理后苜蓿DNA甲基化水平由15.67% 下降到12.85%(P<0.05)。此外，除EMS处理下mCG水平由52.45%上升到54.89%(P<0.05)，mCHG和mCHH水平及60CO-γ射线辐射处理下mCG，mCHG和mCHH水平均呈现下降趋势，由此表明，DNA甲基化主要发生在CG序列，CHH序列在mC中所占比例最小，从而说明不同诱变处理甲基化的变化与mCG的类型变化关系更密切。

表4 不同诱变处理前后C，CG，CHG，CHH甲基化水平Table 4 Methylation levels of C，CG，CHG and CHH in different mutagenesis treatments

2.2 紫花苜蓿全基因组甲基化在基因区的分布规律

不同基因组元件甲基化水平可以反映出基因组的甲基化模式，同时也有利于充分理解与认识甲基化功能机制。图2是紫花苜蓿诱变前后基因不同区域上发生甲基化水平情况，由图可知，在各基因区域上，CG甲基化水平最高，CHG甲基化水平次之，CHH甲基化水平最低。此外还可以发现，EMS处理在紫花苜蓿许多基因区域的CG甲基化水平都有所上升，说明EMS处理在诱变过程中存在多个位点的超甲基化过程；而60CO-γ射线处理在紫花苜蓿许多基因区域的CG甲基化水平都有所下降，说明60CO-γ射线处理在诱变过程中存在多个位点的去甲基化过程。研究结果表明，紫花苜蓿诱变过程中内含子区域发生甲基化水平均较高，其次是外显子区和启动子区域，5′UTR和3′UTR甲基化程度最低。这可能是因为外显子及内含子参与基因表达调控及防止了其序列的异常转录所导致的。

图2 不同基因区域上各类型C的甲基化水平Fig.2 Methylation level of various types of C in different gene regions

2.3 紫花苜蓿DNA甲基化组的差异甲基化及KEGG通路富集分析

2.3.1DNA甲基化组的DMR分析 将苜蓿对照组与EMS处理组(A34-1A VS.A34-2A)及对照组与60CO-γ辐射处理组(A34-1A VS. A34-3A)比较中分别检测出的3 066和1 964个DMRs。将其比对到苜蓿参考基因组中，发现分别有1 212和517个差异甲基化基因被DMR覆盖到，这些基因被称为DMR相关基因或差异甲基化基因(DMG)，其中有105个基因在两组比较中均存在差异甲基化区域(见图3)。将各组比较检测出的DMR比对到不同的基因组元件上(图4)，发现DMR主要分布于启动子区、内含子区和外显子区，而在5′UTR和3′UTR区的分布较少。

图3 紫花苜蓿不同诱变处理差异甲基化基因Fig.3 Differentially methylated genes in each alfalfa mutagenesis comparison

图4 DMR在不同基因组元件中的分布Fig.4 Distribution of DMRs in different genomic elements

2.3.2差异甲基化基因的KEGG通路富集分析 对所得到的DMR相关基因即差异甲基化基因进行KEGG功能富集分析，挖掘出与苜蓿生长及生理变化相关的生物学过程和pathway调控通路，进而探讨差异甲基化基因在紫花苜蓿生长发育中的调控作用。对各组比较的差异甲基化基因进行KEGG通路富集分析，结果表明，对照组和EMS处理组之间差异表达基因显著富集于的Pathway(P<0.1)包括淀粉和蔗糖代谢途径(Starch and sucrose metabolism)、氨基酸糖和核苷酸糖代谢(Amino sugar and nucleotide sugar metabolism)、氮代谢通路(Nitrogen metabolism)嘧啶代谢途径(Pyrimidine metabolism)、精氨酸和脯氨酸代谢(Arginine and proline metabolism)等生理变化相关代谢及通路(图5-1)，涉及到的基因共33个(表5)，其中氮代谢通路、淀粉和蔗糖代谢通路以及脯氨酸代谢通路是与苜蓿逆境生理息息相关的3个主要通路；对照组与60CO-γ辐射组的差异甲基化基因显著富集于植物病原体互作(Plant-pathogen interaction)、鞘脂代谢通路(Sphingolipid metabolism)、代谢通路(Metabolic pathways)、次生代谢产物的生物合成(Biosynthesis of secondary metabolites)、生物合成植物激素信号转导(Plant hormone signal transduction)、及氨基酸类生理代谢变化及通路(图5-2)，涉及到的基因共9个(见表5)，这些通路参与调节细胞的生长、分化、衰老和细胞程序性死亡等许多重要的信号转导过程，使细胞产生各种不同的生物学功能。综合通路富集分析结果可推测紫花苜蓿在受到诱变胁迫条件下，这些通路及相关基因在生理代谢过程中可能发挥着重要的调控作用。

图5-1 A34-1A VS.A34-2A差异表达基因pathway富集分析Fig.5-1 The enrichment analysis of DGEs in the comparison of A34-1A VS.A34-2A

图5-2 A34-1A VS.A34-3A差异表达基因pathway富集分析Fig.5-2 The enrichment analysis of DGEs in the comparison of A34-1A VS.A34-3A

表5 差异甲基化基因富集的与生理代谢相关通路及基因Table 5 The pathways and genes involved in physiological metabolism

3 讨论

诱变能够引起生物体基因组[22]、表达组[23]、蛋白质组[24]和代谢组[25-26]的变化。其中DNA甲基化也被称为基因组中的“第五碱基”，它是维持基因组稳定性和调控基因表达的另一种方式，在植物各种胁迫中以及生长发育过程中扮演着至关重要的角色[27]。有研究报道，DNA甲基化水平往往与物种基因组大小表现出一定的正相关关系[28-29]。前人研究CG位点的甲基化水平变化范围从30.5%的拟南芥(Arabidopsis)到92.5%甜菜(Beet)；CHG甲基化水平的变化范围从9.3%的盐水水芹(Saltwater cress)到81.2%甜菜(Beet)，CHH甲基化水平变化范围从1.1%葡萄(Grape)到18.8%甜菜(Beet)[30]。在本研究中，诱变前后紫花苜蓿基因组在C，CG，CHG和CHH位点的甲基化水平范围分别为12.85%～15.67%，46.62%～54.89%，18.51%～22.68%和5.94%～8.35%。显然与其他物种相比，紫花苜蓿DNA甲基化水平在变化范围之内，且紫花苜蓿基因组在CG和CHG序列环境下也存在相对较高的甲基化水平，且甲基化水平随着不同诱变胁迫条件改变而变化。

此外，本研究发现，C位点的DNA甲基化主要为mCG，mCHH和mCHG 3种类型，其中EMS诱变后苜蓿甲基化水平上升，而60CO-γ辐射处理后苜蓿的甲基化水平下降，且甲基化C碱基中CG，CHG和CHH的分布比例因诱变处理不同而不同，其中CG类型甲基化胞嘧啶所占比例最高，这与Shi等[31]研究发现，成熟期水稻DNA甲基化在重离子辐射处理下发生在CG位点多于CNG位点的结论相一致。此外，也有研究表明CNG位点甲基化变化在植物生长发育和抵御逆境过程中发挥着重要作用[32-34]。进一步表明研究诱变后甲基化水平变化对突变体筛选工作具有重要作用。

不同基因组元件甲基化水平能够反应基因组的甲基化模式，有利于分析甲基化功能机制研究。本实验通过计算各元件内胞嘧啶位点的平均甲基化水平，分析各功能元件的甲基化水平差异。结果发现CG，CHG，CHH各序列环境下基因组各元件甲基化水平均发生了不同程度的变化，其中外显子及内含子区域发生甲基化水平均很高。这可能是因为外显子及内含子在参与基因表达调控为了阻止其序列的异常转录所导致的。这与赵振利研究的白花泡桐丛枝病发生过程中全基因组DNA甲基化差异分析的研究结果相似[35]。

本研究进行对照与EMS处理及对照与60Co-γ射线处理(A34-1A VS.A34-2A和A34-1A VS. A34-3A)的DNA甲基化进行比较，分别检测出的3066和1964个DMR。将其比对到参考基因组中，发现分别有1212和517个差异甲基化基因被DMR覆盖到，且主要分布在启动子区，说明大部分的甲基化差异主要发生在启动子区。将这些基因进行KEGG功能富集分析，发现这些基因显著富集于氮代谢途径、嘧啶代谢途径、多种氨基酸代谢途径及糖代谢等生理变化相关的信号通路中(P< 0.05)。推测这些通路及基因在紫花苜蓿生理代谢过程中发挥着重要的调控作用。进一步说明，DNA甲基化在一定程度上调控基因的表达，是受外界因素调控的众多基因的公共开关，通过改变生理代谢变化，从而促进或者抑制作物生长过程中某些能量物质的代谢与转换，进而影响植物的生活力，协助植物渡过逆境胁迫[36]。

4 结论

本研究利用BS-seq技术对EMS及60CO-γ两种诱变处理下紫花苜蓿全基因组DNA甲基化的变化进行了研究，分析了DNA甲基化的水平和模式变化情况，发现EMS诱变后苜蓿甲基化水平上升，60CO-γ辐射处理后苜蓿的甲基化水平下降，对两组处理进行甲基化区域的差异研究及KEGG通路富集分析，发现了3个与逆境生理相关的信号通路，涉及到的相关基因共有42个。