研究背景和待解决的科学问题

在发育过程中，DNA甲基化和去甲基化的动态变化会对基因表达产生影响，指导细胞特异性的分化方向。5mC修饰主要由脱氧核糖核酸甲基转移酶（DNA methyltransferases，DNMTs) 催化，通常是转录抑制的信号。去甲基化主要由TET蛋白家族介导。通过DNMT和TET蛋白的协同作用，DNA甲基化稳态得到了很好的维持。
人类胚胎干细胞（human embryonic stem cells，hESCs)的定向分化为理解细胞如何对调节因子做出反应提供了体外研究模型。例如，胚胎干细胞向胰腺内胚层的逐步分化可以通过将胚胎干细胞暴露于不同的外源性信号分子来控制。最近有研究表明，与polycomb蛋白或增强子活性相关的关键组蛋白修饰对于hESC-PE分化至关重要。然而，在hESCs向分化为胰腺的过程中，甲基化和去甲基化的变化知之甚少，甲基化和去甲基化如何影响染色质的可及性和转录因子的结合也仍未明确。
将hESC暴露于不同的外部信号分子，hESC逐步向胰腺内胚层分化（pancreatic endoderm, PE）的这一过程模拟了体内胰腺发育的关键步骤。胰腺的体外分化步骤主要是： 因此，本文需要解决的问题是：hESC诱导分化成胰腺细胞的过程中，DNA甲基化水平、染色质可及性，组蛋白修饰发生了哪些动态变化，对分化过程有什么影响？因此，作者开展了以下研究。


研究内容 作者首先按照高效分化的方法将H1-hESC分化成胰腺祖细胞，其间的过渡细胞包括定形内胚层（DE）、原始肠管（GT）和后前肠（FG）细胞。分化的主要阶段如下图所示。


通过RNA-seq，可以得到hESC向胰腺谱系分化的证据：各类细胞的marker基因，下图可以清晰看到各阶段marker基因的表达变化。



为了进行表观层面的研究，通过WGBS-seq、CMS-IP-seq和ATAC-seq，作者获得了不同分化阶段细胞的DNA甲基化位点、DNA羟基化位点以及染色质可及性区域。作者首先分析了DNA甲基化情况，其中5hmC甲基化程度由dot-blot assay完成（下图左）。结果显示，在ES到PE分化过程中，5hmC修饰水平的全局变化呈U型分布（下图右），而5mC修饰水平在整个分化过程中只有2%左右的起伏。这说明某些关键基因座的5mC修饰水平改变决定分化方向，而不像早期胚胎发育过程中DNA甲基化全局变化。



ES分化到DE阶段时，5hmC修饰水平降低了接近70%，可能的原因是TET1和DNMT表达减少（下图左）。而后续分化阶段（GT-FG-PE）5hmC修饰水平逐渐回升，这可能是因为TET2和/或TET3表达量的增加弥补了TET1表达量的减少（下图左）。



作者同时对5个分化阶段分别进行了ATAC-seq，发现染色质可及性变化模式和5hmC修饰水平变化相似（上图右）。
第一部分作者主要总结了DNA甲基化修饰水平和染色质可及性在分化过程中的全局变化，下面作者将开展更加细致的分析和研究。
  为了进一步探索分化过程中各阶段DNA甲基化动态变化情况，作者将所有CpG位点按照每个阶段﹥20%变化进行总结分类，最终分为3类：a. 低甲基化区域（LMS，0-20%）；b. 中甲基化区域（IMS，20-80%）；c. 高甲基化区域（HMS，80%-100%）。在所有CpG区域中，IMS在分化过程中动态变化是最多的（下图），LMS和HMS则保持相对稳定。



与5mC动态变化一致，作者在5hmC甲基化富集区域（HERG）中检测到IMS而不是HMS区域的富集。进一步分析发现IMS区域主要定位于远端调控区域（与TSS平均距离约32 kB），而且染色质状态是高度开放的（下图右）。



同样，按照之前的筛选标准进行汇总和分类，作者将DMR/DHMR区域分成甲基化水平降低（hypo-DMR/DHMR）和甲基化水平升高（hyper-DMR/DHMR）的区域，各分化阶段得到的DMR/DHMR区域数量如下图所示。



为了探究这些DMR和DHMR区域是否和分化阶段具有关联，作者对分化第三阶段的区域所关联的基因进行富集分析，结果显示这些基因主要是胰腺和消化道发育相关的基因（如下图）。



在这些基因中，作者发现了一个典型基因——PTF1A基因座在分化过程中5hmC修饰水平升高，5mC修饰水平降低，染色质可及性增加，这与PTF1A表达水平升高103倍相吻合（下图）。



进一步系统性比较hypo/hyper-DMR/DHMR区域中5mC和5hmC修饰水平发现：1) hypo-DMR中，5mC与5hmC修饰呈负相关（下图左），具体表现在分化过程中，5mC修饰水平的减少伴随5hmC修饰水平的升高，而在hyper-DMR中，这两种甲基化修饰水平没有明显的相关性；2) hyper-DHMR区域中，5mC修饰水平在分化的3个阶段均有所减少（下图右上）；而hypo-DHMR区域5mC修饰水平并没有明显的变化（下图右下）。



这些结果共同说明在分化后期阶段，5hmC修饰占据了5mC去甲基化后的位点；而在整个分化过程中，5hmC修饰水平降低并没有影响5mC修饰水平。
  前面提到IMS区域具有很高的5hmC信号，定位于远端调控元件并且伴随很高的染色质开放程度，由此，作者假设5hmC修饰与增强子活性的调控有关。为了验证这一假设，作者将5hmC富集区域与ATAC-seq数据以及之前发表的H3K4me1、H3K4me3、H3K27ac和H3K27me3的ChIP-seq数据进行联合分析，结果显示5hmC富集区域与H3K4me1信号以及染色质开放区域重叠，而H3K27ac信号关联程度弱于H3K4me1（下图）。另外两种组蛋白修饰H3K4me3和H3K27me3的信号与5hmC修饰在关联程度不高。



为了进一步研究5hmC与其它表观遗传修饰在hESC分化过程中的相关性，作者将5hmC图谱进行分类来鉴定分化阶段特异性的5hmC富集区域，并将其与H3K4me1、H3K4me3、H3K27ac、H3K27me3以及5mC富集模式进行关联，结果显示5hmC与H3K4me1、H3K27ac修饰水平正相关，而与H3K4me3、H3K27me3修饰关联程度较弱（下图）。





一般来说，增强子区域都会伴随H3K4me1和H3K27ac的修饰，作者在研究PDX1（胰腺发育的关键转录因子）的远端调控区域时发现，这个区域在DE分化至PE时，5hmC修饰水平逐渐升高，同时H3K4me1和H3K27ac修饰水平也逐渐升高（下图）。



根据H3K4me1和/或H3K27ac富集情况，增强子在胰腺分化过程中分为三类：a. 阶段性具有活性；b. 一直静止；c. 一直具有活性。作者根据5个分化阶段的H3K4me1/ H3K27ac富集情况进一步将这三类增强子的5hmC修饰特征展现了出来（下图）。结果显示5hmC和H3K4me1在三种类型的增强子中都具有强烈的关联性（H3K27ac稍弱）。



与之前的分析一致，5hmC高度富集区域的5mC修饰水平明显降低，而且H3K4me1/H3K27ac与5mC修饰负相关。这些结果共同说明5hmC修饰占据了增强子中5mC甲基化的位点，且5hmC与决定胰腺细胞谱系分化的增强子有关。
  为了研究分化过程中5hmC如何调控增强子活性，作者调取了静止和具有活性的增强子中5hmC的分布模式。5hmC在静止增强子和活性增强子中都有富集，但活性增强子的富集更加明显（下图上）。另一方面，作者观察到静止增强子和活性增强子区域的甲基化修饰水平整体较低，其中活性增强子的甲基化水平低于静止增强子（下图下）。



为确认以上结果，作者将5hmC的数据与数据库中GRO-seq（GRO-seq结果可以反映增强子活性）数据进行联合分析，发现在分化的各个阶段，5hmC修饰的增强子比5hmC缺失的增强子表现出更高的转录活性（通过更高的GRO-seq信号反映），如下图。



为了进一步揭示5hmC修饰在静止增强子和活性增强子分化过程中的动态变化，作者选择增强子活性发生变化的调控基因组区域（从活性到静止状态；或反之），并检测其5hmC修饰富集情况，结果显示增强子活性的动态变化与5hmC修饰富集的变化之间有很强的相关性，具体表现在在活性增强子变成静止时，5hmC修饰富集水平明显降低（下图）。



同时，作者还分析了5hmC修饰富集与H3K4me1修饰富集重叠的区域，并和通过GRO-seq鉴定到的RNA Poll II介入的基因转录起始区域进行整合，结果同样显示在分化各个阶段，5hmC修饰富集与H3K4me1修饰水平的升高或者GRO-
seq信号呈正相关（下图）。



那么，5hmC与染色质可及性是否有相关性呢？通过ATAC-seq，作者发现5hmC信号与染色质可及性呈正相关，具体表现在基因组5hmC修饰水平升高的区域的染色质可及性也同样升高，反之亦然（下图）。



最后，为了验证5hmC修饰和染色质可及性的变化是否和增强子的活性相关，作者将GRO-seq信号与同时具有H3K4me1 peak、5hmC/ATAC-seq peak的区域关联起来，结果显示5hmC修饰富集、转录起始活性或增强子染色质可及性之间存在明显的正相关（下图）。



上述结果充分证明，在hESC分化至胰腺过程中，DNA羟甲基化与染色质可及性还有增强子活性密切相关。
  考虑到与发育相关的关键转录因子通常与染色质可及的静止/活性增强子区域结合，作者进一步研究了两个先驱转录因子FOXA1和FOXA2的5hmC修饰富集情况，这两个转录因子都属于相同的转录因子家族，并且对内胚层的分化至关重要。虽然5hmC修饰在FOXA1和FOXA2结合区域均有富集（下图），且伴随着开放的染色质状态，但FOXA1结合区域的5hmC修饰富集在DE、GT和PE阶段更为显著和一致。



作者将数据库中已有的FOXA1和FOXA2 ChIP-seq数据和自己的5hmC和5mC数据进行联合分析，结果显示增强子5mC/5hmC修饰的动态变化促使FOXA1或FOXA2结合增强（下图）。



作者对DNA甲基化与FOXA1/FOXA2结合的相关性进行分析，并根据GT、DE和PE不同阶段的DNA甲基化水平划分了它们的结合区域。分析结果揭示了：a. 相比FOXA2结合位点，FOXA1的结合位点甲基化修饰水平更低；b. DE到PE阶段，FOXA1结合位点甲基化修饰减少，而FOXA2没有变化（下图）。



作者又分析了主导胰腺谱系分化的转录因子PDX1在PE阶段时5hmC修饰情况，结果显示PDX1结合位点富集5hmC（下图）。



除了PDX1，作者还将5hmC富集区域进行motif扫描，发现除了FOXA1的motif，CUX1、ONECUT1、ATF1和GMEB1的motif也显著富集，其中CUX1和ONECUT1是胰腺功能相关的关键转录因子。



以上研究表明增强子中5hmC和转录因子的结合位点的DNA甲基化状态改变有关，而且会促进潜在DNA甲基化敏感转录因子的结合，从而调控胰腺谱系特异性基因的表达。
  根据文献报道，5hmC修饰富集在DNA甲基化水平较低的大型基因组区域的边缘（定义为甲基化峡谷，甲基化率<20%，长度>3.5 kb），而且5hmC富集程度与峡谷区域大小的动态变化密切相关。为了进一步了解分化过程中5hmC修饰在大型基因组功能区域富集的情况，作者检测了典型大型基因组功能区域的5hmC修饰的分布情况，包括超级增强子、H3K4me3 broad peak以及甲基化峡谷，结果显示5hmC富集在大型基因组功能区域的边界（定义为5hmC环），而这些区域的5mC修饰水平较低（下图）。



由于H3K4me3修饰的基因组区域的长度与转录延伸和增强子活性相关，故作者检测了H3K4me3 broad peak里的5hmC的修饰情况，结果发现5hmC修饰在边界的富集与所有五个分化阶段的H3K4me3 broad peak的长度正相关（下图上）。这一结果与之前文献报道一致，H3K4me3的peak的长度与相关基因表达量正相关（下图下）。



得到这个规律之后，作者进一步分析了H3K4me3 broad peak的大小动态变化与5hmC修饰水平的关系，结果显示H3K4me3 broad peak的长度收缩与5hmC修饰的减少有关，反之亦然（下图）。



具体到个别基因，作者发现在GT到PE分化阶段，H3K4me3宽峰边界处5hmC水平显著增加（74.36%），同时ONECUT1启动子处peak范围扩大了3.6倍，而且ONECUT1表达上调（下图左）。与ONECUT1相似，RFX6（与胰岛素生成相关的关键转录因子）的启动子区也出现相同的景象（下图右）。



这些结果共同说明5hmC修饰在基因启动子区内H3K4me3 broad peak边界的富集与组蛋白修饰peak的范围和相应的基因表达水平正相关。


总结

1. 通过WGBS、CMS-IP-seq、ATAC-seq获得hESC向胰腺谱系分化各阶段的DNA甲基化/羟基化、染色质可及性的动态变化数据；
2. 通过表观和转录层面的综合性生物信息学分析，发现DNA甲基化/羟甲基化与胰腺谱系决定有关；
3. 在分化过程中，5hmC标记的DNA区域在后续分化阶段的细胞中去甲基化，这与增强子活性和染色质可及性呈正相关，促进谱系决定和先锋转录因子结合，从而诱导相关基因的表达；
4. 5hmC修饰环（5hmC-rim）定位于大型功能基因组区域的边界，包括增强子、DNA甲基化峡谷以及H3K4me3的宽peak中，可能是为了维持该区域中低水平的DNA甲基化。