2025, Vol. 55
2. 中国海洋大学海洋生物多样性与进化研究所，山东 青岛 266003
DNA 5mC甲基化一般是指以S-腺苷甲硫氨酸(SAM) 为甲基供体，在DNA甲基转移酶的催化下，将甲基基团转移到胞嘧啶的5号碳原子的过程[1]，是目前为止研究最为广泛的表观遗传学修饰之一。动物DNA 5mC甲基化的目标序列主要限于回文CpG二核苷酸，通常以对称方式甲基化，而非CpG甲基化在植物中普遍存在[2]。已有研究表明，基因组中DNA 5mC甲基化的调控对基因组稳定性、基因表达和胚胎发育具有重要影响[3-4]。
在人类早期胚胎发育过程中，第一次全基因组去甲基化发生在受精卵至2细胞阶段，并维持到桑葚胚阶段，第二次全基因组去甲基化发生在桑葚胚至囊胚阶段[5-7]。小鼠早期胚胎发育过程中同样发生明显的去甲基化现象。这一过程始于1细胞阶段，表现为全基因组的去甲基化。这种早期的全局低甲基化水平保持了幼体的多能性，并保证了未来准确的分化调节[8-10]，之后在谱系特异性区域中重新建立甲基化[11]。无脊椎动物海葵早期胚胎发育过程中不存在明显的DNA 5mC甲基化动态变化，其精子和卵细胞的整体DNA 5mC甲基化水平相似，各种基因组元件的DNA 5mC甲基化水平也没有检测到显著差异[12]。在海胆中，从配子时期到早期胚胎都观察到DNA 5mC甲基化水平的动态变化，各种基因组元件的DNA 5mC甲基化水平也变化剧烈[12]。值得注意的是，海葵和海胆的基因组中均会显示出马赛克的DNA 5mC甲基化修饰模式[12-13]，即基因组由高度甲基化区域和非甲基化区域穿插组成。在模式动物果蝇和线虫的基因组中，5mC甲基化水平均非常低，果蝇胚胎基因组DNA 5mC甲基化的缺失对胚胎发育也不会产生重大影响[14-15]。总而言之，DNA 5mC甲基化在不同物种中的分布存在一定的差异，因此需要在更多的物种中研究其对早期胚胎发育的作用和意义。
海鞘属于脊索动物门尾索动物亚门[16]，在进化上被认为是脊椎动物的近亲。玻璃海鞘(Ciona robusta)和萨氏海鞘(Ciona savignyi)都拥有较小的基因组且早期胚胎发育速度极快，当胚胎发育至8细胞时期时，整个发育命运就已经被决定，属于经典的镶嵌型发育式生物。玻璃海鞘和萨氏海鞘所拥有的独特的进化地位、快速的早期胚胎发育、较小的基因组以及成体易获取等特性，使得它们成为研究早期胚胎发育的理想模型。Suzuki等[17]曾以玻璃海鞘为研究对象，用质谱法检测了玻璃海鞘早期胚胎中的8细胞期、32细胞期、112细胞期、幼虫时期以及成体的精子和肌肉组织的DNA 5mC甲基化整体水平，结果显示，DNA 5mC甲基化在早期发育胚胎以及成体组织中虽有变化，但变化微小，且证实在差异甲基化位点相关基因中不包含合子基因。因此本研究认为在玻璃海鞘早期胚胎发育中，DNA 5mC甲基化不调控合子基因的激活。Xu等[12]通过对萨氏海鞘的精子、卵母细胞、卵裂时期以及原肠胚进行全基因组亚硫酸氢盐测序(Whole genome bisulfite sequencing, WGBS)发现: 在配子时期以及早期胚胎发育中，萨氏海鞘全基因组DNA 5mC甲基化存在动态变化，且全基因组甲基化是由高度甲基化区域和非甲基化区域穿插组成，具有马赛克分布的特征; 在萨氏海鞘中，甲基化会优先发生在基因体区域而非转座子区域，并且启动子在发育过程中始终保持低甲基化状态。以往并未对DNA 5mC甲基化的分布特征和功能进行深入分析和预测。
本研究以萨氏海鞘为研究对象，通过全基因组亚硫酸氢盐测序技术构建了萨氏海鞘早期胚胎发育过程中32细胞期、112细胞期、原肠胚期和尾芽中期胚胎的全基因组DNA 5mC甲基化图谱，预测了DNA 5mC甲基化在萨氏海鞘早期胚胎发育中的调控功能，为处于无脊椎动物和脊椎动物之间的进化节点生物DNA 5mC甲基化研究奠定了基础，同时也为脊椎动物甲基化调控的起源研究提供了参考。
1 材料与方法 1.1 海鞘胚胎不同发育阶段样本的收集本研究中所用的萨氏海鞘采自山东省荣成市寻山集团养殖区和山东省青岛市黄岛区附近海域。将采集获得的活体动物暂养在实验室内部养殖系统中，温度依据季节设置在(18±2) ℃，维持海水中NaCl浓度在(30±1)‰，养殖过程中持续供氧，并根据生活习性设置光照时间并定期进行投喂，之后选取处于最佳生长状态且精子或卵子处于成熟期的个体进行后续实验。取性腺发育成熟的海鞘成体数只，分别吸取足量的精子或卵子进行人工体外受精，保证胚胎同步发育。将受精卵置于培养皿中，然后放于16 ℃恒温培养箱中进行培养。收集萨氏海鞘32细胞期(受精后3 h)、112细胞期(受精后4.5 h)、原肠胚期(受精后6 h)和尾芽中期(受精后10 h)共4个发育时期的胚胎，装入1.5 mL RNase-free的离心管中，用液氮保存。每个发育时期设置3组生物学重复，共计12个样品，总计大约32 000只胚胎用于后续测序分析。
1.2 DNA提取及WGBS建库和测序将胚胎捣碎并加入600 μL组织裂解液，之后加入8 μL浓度为20 mg/mL的蛋白酶K溶液，充分混匀至液体呈现透明状，再加入200 μL乙酸铵溶液，混匀后冰上放置5~10 min，4 ℃下12 000 r/min离心15 min。吸取约600 μL上清液于另一干净的离心管中，加入等体积(约600 μL)的异丙醇，充分转摇，然后4 ℃下12 000 r/min离心15 min，最后用无水乙醇清洗DNA沉淀，晾干后将DNA溶解在20 μL双蒸水中，置于-20 ℃冰箱中保存备用。
WGBS建库和测序由北京贝瑞和康生物技术有限公司完成。具体流程如下: 样本DNA质检合格后，将掺有0.5% 未甲基化Lambda DNA(Promega)的基因组DNA超声打断至各片段长度100~500 bp，然后将片段DNA用末端修复酶混合物(End repair enzyme mix, NEB)、Klenow 3′-5′exo-(NEB)和T4 DNA连接酶(NEB)分别进行末端修复、dA拖尾和连接。使用EZ DNA methylation-Gold kit试剂盒进行亚硫酸氢盐转化，最后用KAPA HiFi HotStart Uracil+ ReadyMix高保真DNA聚合酶预混液扩增亚硫酸氢盐处理的DNA。测序在Illumina Novaseq 6000平台上进行，测序模式为PE150。
1.3 萨氏海鞘胚胎的5-氮胞苷处理用二甲基亚砜(Dimethyl Sulfoxide, DMSO)将5-氮胞苷粉末配置成100 mmol/L的抑制剂母液。将萨氏海鞘受精卵置于小培养皿中，然后放在16 ℃恒温培养箱中，每个皿中约50只胚胎，设置不同浓度(20、40、100、150和200 μmol/L)的5-氮胞苷处理萨氏海鞘胚胎，对照组为加入等量的DMSO，每隔约1 h，将胚胎从16 ℃恒温培养箱中取出，观察其发育状态。实验组和对照组都设置3个生物学重复。
1.4 点杂交实验将DNA原液经热变性处理后点在尼龙膜上，孵育5mC一抗(全称: 5-Methylcytosine (5-mC) antibody (mAb); 品牌: Active Motif; 货号: 39649;用1×PBST溶液按照1∶4 000稀释)和二抗(全称: ProteinFindⓇ Goat Anti-Mouse IgG (H+L), HRP Conjugate; 品牌: 全式金; 货号: HS201-01;用5%脱脂奶粉液体按照1∶8 000稀释)，然后用ECL显色液(品牌: PerkinElmer; 货号: NEL105001EA)显色，用化学发光成像系统(品牌: 上海天能)曝光后进行拍摄，之后清洗尼龙膜上的显色液，取5 mL甲基蓝溶液(品牌: 生工; 货号: A610622-0025)覆盖尼龙膜孵育5 min后可显现出信号，信号的深浅代表DNA上样量。
1.5 RNA提取及转录组建库和测序将胚胎捣碎并加入1 mL RNAsio Plus (品牌: Takara，货号: 9109)充分裂解，之后加入200 μL的氯仿，涡旋震荡40 s，在4 ℃下12 000 r/min离心15 min; 吸取上清液约500 μL至新的离心管中，加入500 μL异丙醇(与上清液体积1∶1)充分混匀后，于4 ℃下12 000 r/min离心15 min; 用无水乙醇清洗RNA沉淀，之后将RNA溶解在DEPC处理后的双蒸水中，检测RNA的完整性，以挑选完整性高的RNA样品，放于-20 ℃冰箱中保存待用。本研究中将RNA样本交给北京百迈客生物科技有限公司进行转录组建库和测序。首先由公司对RNA样品进行完整性检测，质检合格后进行文库构建，在上机进行高通量测序前，对文库进行芯片制备，将文库模板固定到芯片上，每个RNA分子会形成一个簇，一个簇就是一个测序位点，测序读长为PE150。
1.6 甲基化数据生物信息学分析WGBS原始下机数据为经过一级质控的原始数据(Raw data)，将Raw data使用fastp软件[18]进行二级质控检查后，使用Trimmomatic软件[19]进一步过滤掉接头和低质量碱基序列，参数设置为LEADING=3，TRAILING=3，SLIDINGWINDOW=4∶15，MINLEN=36。之后使用Bismark软件[20](内部默认调用Bowite2软件[21])进行基因组比对分析，获得甲基化位点数据。在计算甲基化水平及后续分析时，仅保留reads数不小于5的位点。基因组平均甲基化水平计算方法为基因组所有甲基化的胞嘧啶位点数除以基因总的胞嘧啶位点数; 特定区域内平均甲基化水平为该区域内所有甲基化的胞嘧啶位点数除以位于该区域胞嘧啶位点数。使用DSS软件包[22]鉴定差异甲基化位点，对鉴定出的差异甲基化位点进行Wald检验，获得p值，之后根据p值的严格筛选和差异甲基化区域的鉴定条件来确定差异甲基化区域。在使用DDS软件的smoothing模块过程中鉴定参数设置如下: p.threshold=1×e-5，delta=0.1，minlen=50(差异甲基化区域最短长度)，minCG=3(差异甲基化区域所包含的最低CpG位点数)，pct.sig=0.5，dis.merge=100。使用Chipseeker软件[23]注释差异甲基化区域，定义被差异甲基化区域覆盖的基因为差异甲基化基因。
甲基化图谱是使用Tbtools-Ⅱ软件[24]中的Advanced Circos功能模块绘制而成的，每圈以10 kb为滑动窗格，以柱状图的形式来展示该区域内甲基化水平的均值，柱子高度代表该区域的平均甲基化水平。
1.7 转录组数据生物信息学分析原始下机结果Raw data用fastp软件[18]进行质量检测，获得干净的reads，使用HISAT2(v2.0.5)软件[25]将测序数据与参考基因组进行精确比对后，计算每个样本中的基因表达矩阵。本研究采用TPM值表示基因的表达量。基因本体(Gene ontology，简称GO)富集分析使用clusterProfiler软件[26]，同时调用ggplot2软件[27]进行可视化和作图，以研究输入基因的富集结果。利用DESeq2软件[28]进行差异表达基因的鉴定，首先对基因表达矩阵进行标准化，并使用二项分布检验去计算不同样本比较之间的p值和差异倍数(Fold change，FC)，差异基因筛选标准如下: 符合log2(Fold change)≥1且p≤0.05的基因为表达量显著上调基因; 符合log2 (Fold change)≤-1且p≤0.05的基因为表达量显著下调基因。
2 结果 2.1 萨氏海鞘早期胚胎发育不同阶段的DNA 5mC甲基化图谱收集萨氏海鞘早期胚胎发育中的32细胞期、112细胞期、原肠胚期和尾芽中期四个阶段的胚胎，各时期胚胎形态如图 1A所示。对四个发育阶段共12个胚胎样本(每个阶段有3个重复样本)进行WGBS测序后，使用Bismark软件鉴定到甲基化C位点的不同背景类型，发现CG序列在各样本中占比最大(84.7%~86.4%)，而CHG和CHH序列在各样本中占比分别为2.0%~3.0%和10.7%~12.4%(见表 1)，所以萨氏海鞘DNA 5mC甲基化以mCpG为主要模式，后续分析仅针对mCpG甲基化进行。
|
(A: 萨氏海鞘早期胚胎发育各时期的形态; B: 从里层到外层依次为萨氏海鞘染色体，以及32细胞期的3个样本、112细胞期的3个样本、原肠胚期的3个样本和尾芽中期的3个样本的甲基化位点分布; C: 萨氏海鞘不同发育阶段甲基化水平。代表 0.001< p < 0.01，代表 0.01< p < 0.05。A: Morphology at various stages of early embryonic development in C. savignyi. B: From the inner layer to the outer layer, there are 13 chromosomes of C. savignyi, 3 samples each at the 32-cell, 112-cell, gastrula and mid-tailbud stage. C: Methylation levels at different developmental stages in C. savignyi. represents 0.001< p < 0.01, represents 0.01< p < 0.05.) 图 1 萨氏海鞘早期胚胎发育四个时期形态及全基因组DNA 5mC甲基化图谱构建 Fig. 1 Construction of morphology and whole genome DNA 5mC methylation atlas for four stages of early embryonic development in C. savignyi |
|
表 1 不同序列背景下甲基化C位点的比例 Table 1 Proportion of methylated C sites in different sequence backgrounds |
在确定萨氏海鞘基因组DNA 5mC甲基化序列类型后，在全基因组水平构建萨氏海鞘早期胚胎发育中四个时期的DNA 5mC甲基化图谱(见图 1B)。结果显示，萨氏海鞘早期胚胎四个发育阶段的各染色体中广泛存在DNA 5mC甲基化的修饰。
接下来对四个发育时期的全基因组甲基化水平进行计算，从而发现32细胞期全基因组DNA 5mC甲基化水平最高(约22%)，112细胞期、原肠胚期和尾芽中期基因组DNA 5mC甲基化水平较低(约18%)。对四个发育时期的DNA 5mC甲基化水平进行显著性检验发现，32细胞期基因组DNA 5mC甲基化水平显著高于112细胞期、原肠胚期和尾芽中期，但112细胞期、原肠胚期、尾芽中期之间无显著性差异(见图 1C)。
2.2 萨氏海鞘早期不同胚胎发育时期的差异甲基化分析用DSS软件鉴定各时期间的差异甲基化区域，结果显示，各时期间均存在差异甲基化区域，其中在32细胞期至112细胞期检测到的差异甲基化区域最多，有1 285个区域属于甲基化水平升高区域(Hyper-differentially methylated region, Hyper-DMR)，有7 720个区域属于甲基化水平降低区域(Hypo-differentially methylated region, Hypo-DMR); 在112细胞期至原肠胚期，有2 505个区域属于Hyper-DMR，有3 199个区域属于Hypo-DMR; 在原肠胚时期至尾芽中期，有5 059个区域属于Hyper-DMR，有2 501个区域属于Hypo-DMR(见图 2A)。
|
(A: 萨氏海鞘早期胚胎不同发育时期间差异甲基化区域数目柱状图; B: 差异甲基化基因GO富集柱状图。A: Bar graph of the number of differentially methylated regions during different developmental stages of early embryos in C. savignyi; B: Bar graph illustrating GO enrichment analysis of differentially methylated genes.) 图 2 萨氏海鞘早期胚胎发育不同时期间差异甲基化分析 Fig. 2 Differential methylation analysis during different stages of early embryonic development in C. savignyi |
将差异甲基化区域注释到基因上可以得到差异甲基化基因，各时期间鉴定到的差异甲基化基因数目如表 2所示。为了解差异甲基化基因的功能，揭示萨氏海鞘早期胚胎发育过程中DNA 5mC甲基化的潜在功能调控机制，对以上注释到的差异甲基化基因进行GO富集分析，挑选各时期p.adjust值小于0.05的条目进行展示。GO富集结果显示: 32细胞期与112细胞期之间的差异甲基化基因会显著富集在泛素蛋白转移酶活性、丝氨酸型内肽酶活性这两个分子功能中; 112细胞期与原肠胚时期之间的差异甲基化基因显著富集在泛素蛋白转移酶活性和ATP依赖性微管运动活性这两个分子功能中; 原肠胚时期与尾芽中期之间的差异甲基化基因在分子功能中显著富集于泛素蛋白转移酶活性、蛋白激酶活性、ATP依赖性微管运动活性和蛋白质丝氨酸/苏氨酸激酶活性，在生物过程中显著富集于蛋白质磷酸化和GTP酶活性的正向调节(见图 2B)。对以上各时期间差异甲基化基因富集到的功能条目进行梳理发现，泛素蛋白转移酶活性通路在不同发育时期均会被富集，这表明在萨氏海鞘早期胚胎发育过程中DNA 5mC甲基化修饰与蛋白泛素化修饰存在紧密联系。
|
|
表 2 萨氏海鞘早期胚胎不同发育时期的差异甲基化基因数目 Table 2 Differentially methylated genes at different developmental stages |
5-氮胞苷(5-Azacytidine)是DNA甲基转移酶(DNA methyltransferases, DNMT)的抑制剂，可以特异性作用于DNA甲基转移酶，致使细胞中DNA甲基转移酶活性的缺失和DNA去甲基化[29]。本研究中通过5-氮胞苷抑制萨氏海鞘全基因组DNA 5mC甲基化水平，以探讨其胚胎发育对表型及相关功能通路的影响。
实验结果显示，200 μmol/L 5-氮胞苷处理组的胚胎在112细胞时期出现了显著畸形，畸形的胚胎细胞排列紊乱，不会呈现左右对称的状态(见图 3A)，这种畸形胚胎数目约占胚胎总数的17%，而对照组和其他浓度处理组的胚胎未观察到显著畸形。点杂交实验检测发现，与对照组相比，200 μmol/L 5-氮胞苷处理组的DNA 5mC甲基化水平显著降低(见图 3B)。
|
(A: 5-氮胞苷处理组在112细胞期畸形胚胎和对照组正常胚胎的形态比较; B: 使用5mC特异性抗体通过点杂交分析5-氮胞苷处理组和对照组的全基因组DNA甲基化水平。“5-AzaC”表示5-氮胞苷处理组; “DMSO”表示加入DMSO的对照组; “Anti-5mC”表示使用了5mC特异性抗体的样品; “Input”表示在尼龙膜上点入相同量的基因组DNA作为上样量对照。A: Comparison of aberrant embryonic morphology at the 112-cell stage between the 5-azacytidine treated group and the control group. B: Analysis of global genomic DNA methylation levels in the 5-azacytidine treatment group and the control group using dot blot with a 5mC-specific antibody. "5-AzaC" represents the 5-azacytidine treatment group; "DMSO" represents the addition of DMSO to the control group; "Anti-5mC" represents samples using a 5mC-specific antibody; "Input" represents the loading of equal amounts of genomic DNA onto the nylon membrane as a loading control.) 图 3 5-氮胞苷抑制萨氏海鞘全基因组DNA 5mC甲基化表型分析结果 Fig. 3 Phenotypic analytic result of 5-azacytidine inhibiting whole genome DNA 5mC methylation in C. savignyi |
为研究全基因组DNA 5mC甲基化水平被抑制后受影响的基因及其通路，对5-氮胞苷处理后112细胞期的胚胎和112细胞期对照组的胚胎进行二代转录组测序(每组3个生物学重复)。测序结束后，计算每个基因在不同样本的TPM值，利用DESeq软件分析5-氮胞苷处理组和对照组的差异表达基因。结果发现，与对照组相比，5-氮胞苷处理后显著影响到925个基因的差异表达，其中有192个基因表达显著升高，733个基因表达显著降低(见图 4A)。对基因表达量显著降低的733个基因进行GO富集分析发现，其主要影响的通路是G蛋白偶联受体、DNA结合转录因子活性(见图 4B)。对基因表达量显著升高的192个基因进行GO富集分析发现，其主要影响的通路是乙酰转移酶活性(见图 4C)。
|
(A: 5-氮胞苷处理组与对照组差异表达基因火山图; B: 5-氮胞苷处理后显著下调基因GO富集点图; C: 5-氮胞苷处理后显著上调基因GO富集点图。A: Volcano map of differentially expressed genes; B: Dot plot depicting GO enrichment of significantly downregulated genes; C: Dot plot depicting GO enrichment of significantly upregulated genes.) 图 4 萨氏海鞘早期胚胎5-氮胞苷处理组与对照组转录组差异分析结果 Fig. 4 Analytic result of transcriptome differences between the 5-azacytidine treated group and the control group of early embryos in C. savignyi |
为验证RNA-seq测序和分析结果的准确性，本研究从显著富集的G蛋白偶联受体信号通路、G蛋白偶联受体活性、DNA结合转录因子活性和乙酰转移酶活性通路中随机选取12个显著差异表达基因，通过荧光定量PCR技术进行验证。结果显示，这些基因在5-氮胞苷处理组和对照组中的表达变化趋势与RNA-seq数据一致，表明本研究中RNA-seq数据分析结果具有较高准确性和可靠性(见图 5)。
|
(A.Alpha-1A肾上腺素受体Alpha-1 A adrenergic receptor; B.转录因子AP-1 Transcription factor AP-1; C.黏附型G蛋白偶联受体G7 Adhesion G-protein coupled receptor G7; D.D(2)多巴胺受体A D(2) dopamine receptor A; E.G蛋白偶联受体84 G-protein coupled receptor 84 evm.model.5.873; F.视黄酸受体alpha Retinoic acid receptor alpha; G.含ETS结构域的ERF样转录因子ETS domain-containing transcription factor ERF-like; H.大麻素受体1 Cannabinoid receptor 1; I.G蛋白偶联受体84 G-protein coupled receptor 84 evm.model.10.865; J.芳胺N-乙酰转移酶1 Arylamine N-acetyltransferase 1; K.芳胺N-乙酰转移酶Arylamine N-acetyltransferase; L.芳胺N-乙酰转移酶松果体同工酶NAT-3 Arylamine N-acetyltransferase, pineal gland isozyme NAT-3. 代表p < 0.001，代表 0.001< p < 0.01，代表 0.01< p < 0.05。represents p < 0.001, represents 0.001< p < 0.01, represents 0.01< p < 0.05.) 图 5 萨氏海鞘早期胚胎5-氮胞苷处理组(5-AzaC)与对照组(CT)差异表达基因的表达模式验证 Fig. 5 Verification of the expression patterns of differentially expressed genes between the 5-azacytidine treated group (5-AzaC) and the control group (CT) in C. savignyi |
本文通过WGBS测序手段构建了萨氏海鞘早期胚胎发育四个时期的甲基化图谱，揭示了mCpG是萨氏海鞘基因组DNA 5mC甲基化的主要模式，这与大部分动物基因组中DNA 5mC甲基化模式相似，即DNA 5mC甲基化主要存在于回文CpG二核苷酸上，通常是以对称的方式甲基化，而非CpG甲基化在植物中更为普遍[2]。对四个时期全基因组DNA 5mC甲基化水平统计分析发现，32细胞期全基因组DNA 5mC甲基化水平最高(22%)，112细胞期、原肠胚期和尾芽中期全基因组DNA 5mC甲基化水平较低(均约18%)，这与Greer等[15]发现的萨氏海鞘配子时期及早期胚胎发育中全基因组DNA 5mC甲基化存在动态变化的结论一致。到目前为止，动物配子时期及早期胚胎发育时期全基因组DNA 5mC甲基化水平在脊椎动物和无脊椎动物中存在一定的差异[30]。在脊椎动物如小鼠和人类中，DNA 5mC甲基化几乎发生在整个基因组中[31-33](甲基化水平>0.75);在斑马鱼中，当受精卵经历细胞分裂时，DNA 5mC甲基化减少，并在64细胞阶段达到最低。随后的再甲基化是渐进的，到合子基因组激活时，卵母细胞甲基化被重新调整以匹配精子[34]; 在两栖动物非洲爪蟾(Xenopus laevis)中，父本基因组在受精后不会立即发生去甲基化，卵母细胞和精子基因组中的高水平DNA 5mC甲基化在早期胚胎中保持不变，但在卵裂阶段逐渐下降，在中囊胚转变期和随后的原肠胚形成期间具有最低的DNA 5mC甲基化水平[35]; 而无脊椎动物具有低甲基化模式，例如海葵(Nematostella vectensis)全基因组DNA 5mC甲基化水平约为0.11[12]，栉水母(Pleurobrachia bachei)全基因组DNA 5mC甲基化水平约为0.08[36]，大理石花纹小龙虾(Procambarus virginalis)和牡蛎(Crassostrea gigas)全基因组DNA 5mC甲基化水平约为0.15[37-38, 12]，蜜蜂(Apis mellifera)、蚂蚁(Camponotus floridanus)、黄蜂(Nasonia vitripennis)、蚕(Bombyx mori)、甲虫(Tribolium castaneum)和水蚤(Daphnia pulex)全基因组DNA 5mC甲基化水平约为0.01[39-42]。这些数据表明，不同动物之间全基因组DNA 5mC甲基化水平存在差异，且脊椎动物DNA 5mC甲基化水平普遍高于无脊椎动物。
为探究萨氏海鞘早期胚胎发育过程中的差异甲基化，本研究对不同时期之间差异甲基化基因进行GO富集分析。结果显示，差异甲基化基因在泛素蛋白转移酶活性通路中富集最为显著，这表明在萨氏海鞘早期胚胎发育过程中，DNA 5mC甲基化修饰与蛋白泛素化修饰之间存在紧密联系。蛋白泛素化修饰系统在以往的研究中已被报道可以在配子阶段及早期胚胎发育的关键进程中发挥重要作用[43]。例如，在人类胚胎中，基因组激活后，蛋白泛素化修饰系统能够确保母体蛋白被降解，从而使胚胎基因组接管胚胎后续发育进程[44]。此外，在胚胎着床之后的发育过程中，蛋白泛素化修饰系统还会建立胚胎发育的中线以及关闭神经管，使胚胎能够正常生长发育[43-44]。DNA 5mC甲基化修饰和蛋白泛素化修饰在细胞中都是重要的表观遗传修饰方式，但在海鞘中，两者之间的内在联系及分子机制尚未被揭示，未来深入研究这两者关系可能对揭示海鞘细胞调控网络和早期胚胎发育命运决定机制具有重要意义。
为研究DNA 5mC甲基化在萨氏海鞘早期胚胎发育中的功能，本研究使用DNA 5mC甲基化抑制剂5-氮胞苷成功抑制了萨氏海鞘胚胎全基因组DNA 5mC甲基化水平。结果显示，全基因组DNA 5mC甲基化水平被抑制后会导致胚胎发育畸形，这表明维持一定水平的DNA 5mC甲基化对萨氏海鞘胚胎发育至关重要(但也要注意防止过高浓度的5-氮胞苷处理所导致的非特异性细胞毒性)。在高等脊椎动物如小鼠中，通过敲除或抑制DNA甲基转移酶使得全基因组DNA 5mC甲基化被异常抑制，进而会导致胚胎发育的严重缺陷和死亡[45-46]。在斑马鱼中，将斑马鱼胚胎从受精卵暴露于10 μmol/L 5-氮胞苷溶液中，并一直追踪斑马鱼的胚胎发育状态发现，与对照组相比，5-氮胞苷处理后的第15天，鱼鳔发生异常膨胀，肠道也发育异常，并且在F1代会导致斑马鱼性别发生显著转变，这说明生命体的运转和正常生长发育都需要维持一定水平的DNA 5mC甲基化。为研究全基因组DNA 5mC甲基化水平被抑制后受影响的基因及其通路，本研究还对5-氮胞苷处理后112细胞时期的畸形胚胎和对照组发育正常的胚胎进行二代转录组测序。结果显示，与对照组相比，5-氮胞苷处理显著影响了925个基因的表达，这些基因的主要通路为G蛋白偶联受体、DNA结合转录因子活性和乙酰转移酶活性等。G蛋白偶联受体是一类重要的跨膜蛋白受体，在细胞信号传导过程中发挥着关键作用[47]，因此本研究推测在萨氏海鞘中DNA 5mC甲基化可能通过某些转录因子调节G蛋白偶联受体相关基因的表观遗传修饰(例如通过影响组蛋白修饰、染色质重塑的方式)，进而影响G蛋白偶联受体的功能。另外，虽然甲基化抑制剂抑制DNA 5mC甲基化，但并不一定会使基因表达水平上调，这表明海鞘中DNA 5mC甲基化在胚胎发育中调控基因表达的机制是一个高度复杂的调控网络。
总而言之，DNA 5mC甲基化在早期胚胎发育过程中的调控机制复杂而精密。本研究为处于无脊椎动物和脊椎动物之间的进化节点生物的表观基因组学研究提供参考数据，同时也为脊椎动物DNA 5mC甲基化的起源演化探索提供新的线索。
| [1] |
Sulewska A, Niklinska W, Kozlowski M, et al. DNA methylation in states of cell physiology and pathology[J]. Folia Histochemica et Cytobiologica, 2007, 45(3): 149-158. (  0) |
| [2] |
Law J A, Jacobsen S E. Establishing, maintaining and modifying DNA methylation patterns in plants and animals[J]. Nature Reviews Genetics, 2010, 11(3): 204-220. DOI:10.1038/nrg2719 ( 0) |
| [3] |
Jaenisch R, Bird A. Epigenetic regulation of gene expression: How the genome integrates intrinsic and environmental signals[J]. Nature Genetics, 2003, 33(S3): 245-254. DOI:10.1038/ng1089 ( 0) |
| [4] |
Smith Z D, Meissner A. DNA methylation: Roles in mammalian development[J]. Nature Reviews Genetics, 2013, 14(3): 204-220. DOI:10.1038/nrg3354 ( 0) |
| [5] |
Guo H, Zhu P, Yan L, et al. The DNA methylation landscape of human early embryos[J]. Nature, 2014, 511(7511): 606-610. DOI:10.1038/nature13544 ( 0) |
| [6] |
Fulka H, Mrazek M, Tepla O, et al. DNA methylation pattern in human zygotes and developing embryos[J]. Reproduction (Cambridge, England), 2004, 128(6): 703-708. DOI:10.1530/rep.1.00217 ( 0) |
| [7] |
Okae H, Chiba H, Hiura H, et al. Genome-wide analysis of DNA methylation dynamics during early human development[J]. PLoS Genetics, 2014, 10(12): e1004868. DOI:10.1371/journal.pgen.1004868 ( 0) |
| [8] |
Nichols J, Smith A. Naive and primed pluripotent states[J]. Cell Stem Cell, 2009, 4(6): 487-492. DOI:10.1016/j.stem.2009.05.015 ( 0) |
| [9] |
Theunissen T W, Friedli M, He Y, et al. Molecular criteria for defining the naive human pluripotent state[J]. Cell Stem Cell, 2016, 19(4): 502-515. DOI:10.1016/j.stem.2016.06.011 ( 0) |
| [10] |
Peng G, Suo S, Cui G, et al. Molecular architecture of lineage allocation and tissue organization in early mouse embryo[J]. Nature, 2019, 572(7770): 528-532. DOI:10.1038/s41586-019-1469-8 ( 0) |
| [11] |
Zhang Y, Xiang Y, Yin Q, et al. Dynamic epigenomic landscapes during early lineage specification in mouse embryos[J]. Nature Genetics, 2018, 50(1): 96-105. DOI:10.1038/s41588-017-0003-x ( 0) |
| [12] |
Xu X, Li G, Li C, et al. Evolutionary transition between invertebrates and vertebrates via methylation reprogramming in embryogenesis[J]. National Science Review, 2019, 6(5): 993-1003. DOI:10.1093/nsr/nwz064 ( 0) |
| [13] |
Bird A P, Taggart M H, Smith B A. Methylated and unmethylated DNA compartments in the sea urchin genome[J]. Cell, 1979, 17(4): 889-901. DOI:10.1016/0092-8674(79)90329-5 ( 0) |
| [14] |
Kunert N, Marhold J, Stanke J, et al. A Dnmt2-like protein mediates DNA methylation in Drosophila[J]. Development (Cambridge, England), 2003, 130(21): 5083-5090. DOI:10.1242/dev.00716 ( 0) |
| [15] |
Greer E L, Blanco M A, Gu L, et al. DNA methylation on N6-adenine in C. elegans[J]. Cell, 2015, 161(4): 868-878. DOI:10.1016/j.cell.2015.04.005 ( 0) |
| [16] |
Delsuc F, Brinkmann H, Chourrout D, et al. Tunicates and not cephalochordates are the closest living relatives of vertebrates[J]. Nature, 2006, 439(7079): 965-968. DOI:10.1038/nature04336 ( 0) |
| [17] |
Suzuki M M, Mori T, Satoh N. The Ciona intestinalis cleavage clock is independent of DNA methylation[J]. Genomics, 2016, 108(3-4): 168-176. DOI:10.1016/j.ygeno.2016.10.001 ( 0) |
| [18] |
Chen S, Zhou Y, Chen Y, et al. Fastp: An ultra-fast all-in-one FASTQ preprocessor[J]. Bioinformatics, 2018, 34(17): 884-890. DOI:10.1093/bioinformatics/bty560 ( 0) |
| [19] |
Bolger A M, Lohse M, Usadel B. Trimmomatic: A flexible trimmer for illumina sequence data[J]. Bioinformatics, 2014, 30(15): 2114-2120. DOI:10.1093/bioinformatics/btu170 ( 0) |
| [20] |
Krueger F, Andrews S R. Bismark: A flexible aligner and methylation caller for bisulfite-seq applications[J]. Bioinformatics, 2011, 27(11): 1571-1572. DOI:10.1093/bioinformatics/btr167 ( 0) |
| [21] |
Langmead B, Salzberg S L. Fast gapped-read alignment with Bowtie 2[J]. Nature Methods, 2012, 9(4): 357-359. DOI:10.1038/nmeth.1923 ( 0) |
| [22] |
Feng H, Conneely K N, Wu H. A Bayesian hierarchical model to detect differentially methylated loci from single nucleotide resolution sequencing data[J]. Nucleic Acids Research, 2014, 42(8): e69. DOI:10.1093/nar/gku154 ( 0) |
| [23] |
Yu G, Wang L G, He Q Y. ChIPseeker: An R/Bioconductor package for ChIP peak annotation, comparison and visualization[J]. Bioinformatics, 2015, 31(14): 2382-2383. DOI:10.1093/bioinformatics/btv145 ( 0) |
| [24] |
Chen C, Chen H, Zhang Y, et al. TBtools: An integrative toolkit developed for interactive analyses of big biological data[J]. Molecular Plant, 2020, 13(8): 1194-1202. DOI:10.1016/j.molp.2020.06.009 ( 0) |
| [25] |
Kim D, Paggi J M, Park C, et al. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype[J]. Nature Biotechnology, 2019, 37(8): 907-915. DOI:10.1038/s41587-019-0201-4 ( 0) |
| [26] |
Wu T, Hu E, Xu S, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data[J]. The Innovation, 2021, 2(3): 100141. DOI:10.1016/j.xinn.2021.100141 ( 0) |
| [27] |
Ito K, Murphy D. Application of ggplot2 to pharmacometric graphics[J]. CPT: Pharmacometrics & Systems Pharmacology, 2013, 2(10): 1-16. ( 0) |
| [28] |
Love M I, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2[J]. Genome Biology, 2014, 15(12): 550. DOI:10.1186/s13059-014-0550-8 ( 0) |
| [29] |
Christman J K. 5-Azacytidine and 5-aza-2'-deoxycytidine as inhibitors of DNA methylation: Mechanistic studies and their implications for cancer therapy[J]. Oncogene, 2002, 21(35): 5483-5495. DOI:10.1038/sj.onc.1205699 ( 0) |
| [30] |
Schübeler D. Function and information content of DNA methylation[J]. Nature, 2015, 517(7534): 321-326. DOI:10.1038/nature14192 ( 0) |
| [31] |
Lister R, Pelizzola M, Dowen R H, et al. Human DNA methylomes at base resolution show widespread epigenomic differences[J]. Nature, 2009, 462(7271): 315-322. DOI:10.1038/nature08514 ( 0) |
| [32] |
Feng S, Cokus S J, Zhang X, et al. Conservation and divergence of methylation patterning in plants and animals[J]. Proceedings of the National Academy of Sciences, 2010, 107(19): 8689-8694. DOI:10.1073/pnas.1002720107 ( 0) |
| [33] |
Hon G C, Rajagopal N, Shen Y, et al. Epigenetic memory at embryonic enhancers identified in DNA methylation maps from adult mouse tissues[J]. Nature Genetics, 2013, 45(10): 1198-1206. DOI:10.1038/ng.2746 ( 0) |
| [34] |
Potok M E, Nix D A, Parnell T J, et al. Reprogramming the maternal zebrafish genome after fertilization to match the paternal methylation pattern[J]. Cell, 2013, 153(4): 759-772. DOI:10.1016/j.cell.2013.04.030 ( 0) |
| [35] |
Stancheva I, El-Maarri O, Walter J, et al. DNA methylation at promoter regions regulates the timing of gene activation in Xenopus laevis embryos[J]. Developmental Biology, 2002, 243(1): 155-165. DOI:10.1006/dbio.2001.0560 ( 0) |
| [36] |
Dabe E C, Sanford R S, Kohn A B, et al. DNA methylation in basal metazoans: Insights from Ctenophores[J]. Integrative and Comparative Biology, 2015, 55(6): 1096-1110. DOI:10.1093/icb/icv086 ( 0) |
| [37] |
Gatzmann F, Falckenhayn C, Gutekunst J, et al. The methylome of the marbled crayfish links gene body methylation to stable expression of poorly accessible genes[J]. Epigenetics & Chromatin, 2018, 11(1): 57. ( 0) |
| [38] |
Wang X, Li Q, Lian J, et al. Genome-wide and single-base resolution DNA methylomes of the Pacific oyster Crassostrea gigas provide insight into the evolution of invertebrate CpG methylation[J]. BMC Genomics, 2014, 15(1): 1119. DOI:10.1186/1471-2164-15-1119 ( 0) |
| [39] |
Bonasio R, Li Q, Lian J, et al. Genome-wide and caste-specific DNA methylomes of the ants Camponotus floridanus and Harpegnathos saltator[J]. Current Biology, 2012, 22(19): 1755-1764. DOI:10.1016/j.cub.2012.07.042 ( 0) |
| [40] |
Wang X, Wheeler D, Avery A, et al. Function and evolution of DNA methylation in Nasonia vitripennis[J]. PLoS Genetics, 2013, 9(10): e1003872. DOI:10.1371/journal.pgen.1003872 ( 0) |
| [41] |
Xiang H, Zhu J, Chen Q, et al. Single base-resolution methylome of the silkworm reveals a sparse epigenomic map[J]. Nature Biotechnology, 2010, 28(5): 516-520. DOI:10.1038/nbt.1626 ( 0) |
| [42] |
Song X, Huang F, Liu J, et al. Genome-wide DNA methylomes from discrete developmental stages reveal the predominance of non-CpG methylation in Tribolium castaneum[J]. DNA Research, 2017, 24(5): 445-457. DOI:10.1093/dnares/dsx016 ( 0) |
| [43] |
Cruz Walma D A, Chen Z, Bullock A N, et al. Ubiquitin Ligases: Guardians of mammalian development[J]. Nature Reviews Molecular Cell Biology, 2022, 23(5): 350-367. DOI:10.1038/s41580-021-00448-5 ( 0) |
| [44] |
Zhang Y, Zhao L, Zhang J, et al. DCAF 13 promotes pluripotency by negatively regulating SUV 39H1 stability during early embryonic development[J]. The EMBO Journal, 2018, 37(18): e98981. DOI:10.15252/embj.201898981 ( 0) |
| [45] |
Bourc'his D, Xu G L, Lin C S, et al. Dnmt3L and the establishment of maternal genomic imprints[J]. Science, 2001, 294(5551): 2536-2539. DOI:10.1126/science.1065848 ( 0) |
| [46] |
Bourc'his D, Bestor T H. Meiotic catastrophe and retrotransposon reactivation in male germ cells lacking Dnmt3L[J]. Nature, 2004, 431(7004): 96-99. DOI:10.1038/nature02886 ( 0) |
| [47] |
Dogra S, Yadav P N. Biased agonism at kappa opioid receptors: Implication in pain and mood disorders[J]. European Journal of Pharmacology, 2015, 763: 184-190. DOI:10.1016/j.ejphar.2015.07.018 ( 0) |
2. Institute of Evolution and Marine Biodiversity, Ocean University of China, Qingdao 266003, China