The synergistic coordination of DNA 5-hydroxylmethylcytosine and RNA 5-methycytosine regulates human foetal development

After implantation, complex and highly specialized molecular events render functionally distinct organ formation, whereas how the epigenome shapes organ-specic development remains to be fully elucidated. Here, Nano-hmC-Seal and RNA bisulte sequencing (RNA-BisSeq) were performed and the rst multilayer landscapes of both DNA 5hmC and RNA m 5 C epigenome, as well as the accompanying transcriptome were obtained, in the heart, kidney, liver and lung of the human foetuses at 13 − 28 weeks with 123 samples in total. DNA 5hmC plays a consistently leading role, while RNA m 5 C acts as a dynamic downstream regulatory valve in regulating gene expression with a turning point of 18–19 weeks. Although, these modications appear to coordinate in general, they also play different roles with regard to specic functions. Our integrated studies illustrate the epigenetic maps during different stages of human foetal organogenesis, which provide a foundation for understanding the in-depth epigenetic mechanisms for early development and birth defects. Tissue collection and ethics approval of this study. All of the human foetal organs at different developmental stages were obtained from the donors voluntarily at the First Aliated Hospital of Zhengzhou University with signed informed consent and the approval of the Ethics Committee (license no. 2020-KY-261). Foetuses were at 13 to 28 weeks in the second trimester of pregnancy. Foetal heart, liver, kidney, and lung were collected from selective induced labor fetus of donors with normal chromosome karyotypes and no family history of infectious diseases, hereditary diseases and smoking history. The foetal organs were isolated after a mechanical dissection, serially washed with medical saline for three times, stored in RNAlater Stabilization Reagent (AM7024, Invitrogen) and immediately frozen in liquid nitrogen. The details of this research, including the research procedure, the benets and risks of donation, the condentiality of the donor’s information, the right of withdrawal, the generation and analysis of sequencing data, and the publication of results, were explained to the donors before they gave consent. All the experiments were also under the regulations of International Ethical Guidelines for Biomedical Research Involving Human Subjects and in compliance with the policies and laws of China. in at least one of the four stages. To obtain the enrichment scores of binding sites for CEBPB and TRPS1 with these Alu DhMRs, we rst calculated the density of the intersection of either TF binding regions with these Alu DhMRs. This density was then divided by either the density of the intersection of these TF binding sites with all Alu elements in the human genome, or by the density of these TF binding sites in the whole genome. The enrichment scores of global TF binding sites in all Alu elements were calculated by the ratio of the density of the intersection of corresponding TF binding sites with all Alu elements in genome divided by the density of corresponding TF binding sites in the whole genome.


Introduction
Epigenetic modi cations, including DNA methylation 1-3 , histone modi cations [4][5][6] , chromatin accessibility 7-10 and RNA modi cations 11,12 , not only act as heritable and stable indications for the speci cation of chromatin organization and structure, but also participating in regulation of transcriptional states during mammalian early development. The epigenomes of human embryos and the adult organs in both physiological and pathological conditions have been investigated extensively 3,6,11−15 ; however, the epigenetic network implicated in the normal development of human foetal organs remains unclear.
As an essential intermediate product of targeted, active DNA demethylation, 5-hydroxymethylcytosine (5hmC) is oxidized from 5-methylcytosine (5mC) by the ten eleven translocation (TET) family of enzymes (TET1/2/3) 16 . This modi cation is widely distributed in the locations of promoters, enhancers and gene body, regulating target gene expression 17 . Now accumulating evidence showed that abundant 5hmC and TET enzymes are involved in mediating developmental progression through regulation of the chromatin state and/or transcription in various embryonic and adult cell types, including zygotes, primordial germ cells, Purkinje neurons, and embryonic stem cell (ESCs) [18][19][20][21][22] . In fertilized oocytes, paternal 5mC is converted to 5hmC allowing drastic reprogramming of sperm DNA in mice 23 . However, 5hmC accumulation is not a requirement for the initial loss of paternal 5mC. This phenomenon depends on TET3 activity and the de ciency of zygotic DNMT3A and DNMT1 24 . Although the regulatory power and plasticity of 5hmC have been revealed in human adults 20,25 , its function and reprogramming are still unclear during human foetal organ development after implantation.
In addition to DNA modi cations, RNA methylations, especially N 6 -methyladenosine (m 6 A) and 5methylcytosine (m 5 C), contribute to mammalian development at post-transcriptional regulation level [26][27][28][29] . The RNA m 6 A pro le in human adults 15 and mixed foetal samples from 18-25 weeks 12 has been investigated, and the tissue-speci c differences were found to be highly correlated with developmental processes 12 . Another prevalent modi cation, m 5 C, is reported to regulate mRNA nuclear export 28,30 , mRNA stability 31,32 and disease pathogenesis 33 . The conserved, organ-speci c and dynamic features across mammalian transcriptomes of m 5 C indicate its potential importance during mammalian development 28 .
It has been hypothesizsed that the modi cations on both DNA and RNA may coordinately work together in an intricate regulatory network. As histone H3 trimethylation at Lys36 (H3K36me3) and H3 trimethylation at Lys27 (H3K27me3) could guide or impede RNA m 6 A deposition in mice and humans 34,35 , this issue has been a subject of intense investigation in aim to elucidate the cross-talk between epigenetic factors in physiological and pathological conditions such as development and diseases. However, it is still elusive how the epigenome spatiotemporally shapes organ-speci c development in a cooperative manner and further affects normal organ development in humans.
To address these questions, we have established the rst pro les of genome-wide DNA 5hmC and single base resolution RNA m 5 C of various samples including heart, kidney, liver and lung from human foetuses at 13-28 weeks. Importantly, we found that both DNA 5hmC and RNA m 5 C not only show organ-and stage-speci c characteristics in human fetuses, but also exhibit differential and synergistic roles during human foetal development.

Results
Genome-wide pro ling of DNA 5hmC in human foetal organs. We collected samples of 4 major organs including, heart, kidney, liver, and lung, from 8 human foetuses at 13-15, 18-19, 21-23 and 25-28 weeks with 2-4 replicates for each stage ( Fig. 1a and Supplementary Table 1). Nano-5hmC-Seal 36 was rst conducted for DNA 5hmC pro ling and displayed high reproducibility between replicates in different periods (Fig. S1a). To investigate genome-wide 5hmC features, we rst examined the distributions of normalized reads of several development-related genes (Fig. 1b). The 5hmC enrichment of these genes prompted us to consider its organ-speci c features in human foetuses at different stages.
To evaluate 5hmC global changes, we then illustrated 5hmC distributions along gene bodies (± 3 kb), which were similar to those of human adults 17,25 (Fig. 1c). Most of the 5hmC peaks were localized on vital functional elements, such as promoters (Fig. S1b). These results showed that the distribution of DNA 5hmC was conserved in human foetuses and adults, similar to other mammals 17,37 . Moreover, the distinct signal intensity of 5hmC among stages in each organ indicated that 5hmC may play different roles during organ development (Fig. 1c). For instance, in the foetal heart, the 5hmC signal was relatively more stable than that of other three organs probably because it was almost formed at 8 weeks and displayed fundamental biological functions (Fig. 1a, c). In the liver, the 5hmC signal was higher when haematopoiesis began, and even higher when the liver showed a strong haematopoietic ability (Fig. 1a,  c). Together, these results indicate that DNA 5hmC may dynamically regulate organ development by altering the activity of functional elements in an organ-speci c manner.
We next identi ed differentially 5hmC-modi ed bins among four organs (ANOVA, fold change ≥ 1.5, P value < 0.05) (Fig. S1c), and the majority of these bins showed higher 5hmC levels than those of other organs (84.71-99.80%). Then adjacent bins were merged into one region called differentially hydroxymethylated regions (DhMRs) and 173,330-485,961 organ-speci c DhMRs were obtained with around100-200 nt in length (Fig. S1d, S1e). No obvious preference of 5hmC modi cation was detected between gene body and intergenic regions, which was the same as for different chromosomes (Fig. S1f, g). Consistent with the annotation results of the 5hmC peaks in each organ, most of the organ-speci c DhMRs were localized in functional elements (Fig. 1d). The merged top 1% least signi cant bins were de ned as common hydroxymethylated regions (ChMRs) acting as controls. As expected, 5hmC enrichment in DhMRs was different across the 4 types of organs. For example, in the liver, 5hmC enrichment in promoters was higher than that of other organs (Fig. 1d). Moreover, most of the organspeci c DhMRs were located near transcription start sites (TSSs), with prominent in the liver (Fig. 1e). These results indicated that mammalian-conserved 5hmC modi cations are widespread but show signi cant organ-speci c feature in human foetuses.
We further identi ed 10,658-20,419 organ-stage-speci c DhMRs (ANOVA, fold change ≥ 1.5, P value < 0.05) and found that samples from the same organ clustered more closely, whereas that samples from the same stage were scattered on the periphery ( Fig. 1f and Supplementary Table 2). The results indicated that 5hmC could serve as a potential marker speci c for human organs with different developmental stages. Moreover, cluster analysis showed that both 5hmC signal intensity and the corresponding functions of the DhMRs were both organ-and stage-speci c (Fig. 1g, h). Taken together, these results suggest that organ-and stage-speci c 5hmC modi cation is associated with morphogenesis and organ formation during early human development.
The reprogramming of DNA 5hmC during organ development. To further explore the dynamic sequential changes of 5hmC, we rst classi ed the DhMRs into 6 clusters according to the overall changes of 5hmC signal in different stages, considering that it contains all types of nonrepetitive trends for 6 clusters (Figs. 2a and S2a). Most of the DhMR-containing genes belong to at least 2 clusters, suggesting that the same gene may be regulated through multiple pathways in a 5hmC-dependent manner (Fig. S2b). Therefore, we separated these DhMRs from different clusters into 2 groups: a single-regulation group (gene sets with DhMRs only belonged to 1 cluster) and a mixed-regulation group (gene sets with DhMRs belonged to at least 2 clusters) (Fig. 2b). Then, TF prediction was performed on the DhMRs from both groups, and the top 15 TFs signi cantly related to development and three most commonly enriched TFs are shown (Fig. 2b). To explore how TFs are involved in 5hmC-dependent organ developmental processes, we performed RNA-Seq on samples used for the Nano-hmC-Seal data (Supplementary Table 3). There was a strong correlation between replicates (Fig. S2c). In both the single-and mixedregulation groups, development-related TFs, including T-box transcription factor 20 (TBX20), paired box protein Pax-8 (PAX8), Krueppel-like factor 1 (KLF1), and transcription factor 21 (TCF21), were highly enriched on DhMRs with organ-speci c expression feature (Fig. 2b). Notably, a high 5hmC signal was detected on the DhMRs from cluster 5 in the heart at 21-28 weeks, cluster 3 in the kidney at 18-23 weeks, cluster 1 in the liver at 13-15 weeks, and cluster 3 in the lung at 18-23 weeks, which were speci cally enriched by TBX20, PAX8, KLF1 and TCF21, respectively. For example, TBX20 is a transcriptional activator/repressor required for cardiac development 38 . The signi cant enrichment of this TF suggests its key role in the later stage of development, conferring full maturity of foetal heart. PAX8 was reported to be highly associated with the mesenchymal to epithelial transition involving metanephros morphogenesis 39 . From 18-23 weeks, corticomedullary differentiation is ongoing, and the nephron frequently increases simultaneously. PAX8 may mainly function during this period. A signi cantly higher binding of most TFs was also found in human adults, but the expression level of those TFs in foetuses was obviously higher than that in adults (Fig. S2d, e). Collectively, these results illustrate that TFs may participate in speci c organ development with the assistance of sequential changes in 5hmC in foetuses.
DNA methylation is one of the major mechanisms adopted by the host to suppress the expression of transposable elements (TEs) 40 , some of which, however, could be used by their host genome as regulatory elements in certain organs during development [41][42][43] . Therefore, we investigated whether 5hmC modi cations are preferentially installed on certain families of TEs in an organ-and stage-speci c manner, which potentially facilitates the switch of the repressive state to the accessible state of the chromatin for TF embedding and subsequent regulation of gene expression in cis. We found that members of the Alu family of SINE elements show signi cant 5hmC enrichment on the DhMRs within cluster 3 and 4 in the foetal kidney (Figs. 2c and S2g-i). In particular, the binding sites of two TFs, namely CEBPB and TRPS1, are highly enriched in TSS-proximal Alu elements within these DhMRs (Figs. 2d and S2j), indicating that certain TFs preferentially bind to these Alu DhMRs, and their binding and cis-regulatory effect on target genes may be affected by differential 5hmC modi cations of these Alu elements. The pattern of expression level changes of CEBPB target genes was consistent with its role as an activator in transcriptional regulation (Fig. 3e). Indeed, the expression levels of almost all CEBPBtargeted Alu-proximal genes (12/13) followed the 5hmC signals of the corresponding Alu DhMRs (Fig. 3e). Interestingly, approximately half of these CEBPB-targeted Alu-proximal genes (7/13) have been previously reported to be tightly associated with kidney development or responses to injury (Supplementary Table 4). This result indicates that 5hmC can epigenetically mark a transition state of Alu elements, through which to affect the binding of TFs and further downstream gene expression.
DNA 5hmC is correlated with gene expression homeostasis across human organs. Most DhMRs were located around the TSSs in all foetal organs (Fig. 1e) and the sequential changes of 5hmC signal on DhMRs showed a strong association with RNA abundance during organ development (Fig. 2b, e). We then asked whether the dynamic reprogramming of DNA 5hmC on promoters regulates gene expression at developmental stages. In support, a positive correlation between transcriptional activity and the 5hmC signal on promoters was observed in foetal organ samples (Fig. 3a), and most of the promoters contained 1-2 DhMRs (Fig. S3a). Importantly, their association was not set in stone during the development of each organ, and higher correlation coe cients appeared at the later stages of almost all four organs (Fig. 3a). Then we identi ed genes with signi cant 5hmC alterations in promoters across adjacent stages (changes ≥ 25%, fold change ≥ 2, Student's t-test, P value < 1.03 × 10 − 6 ) (Figs. 3b and S3b, c), and most of them were organ-speci c promoters with obvious 5hmC changes (Fig. 3c). In particular, more than half of differentially expressed genes, such as those in the heart at 18-23 weeks (77.06%), liver at 18-23 weeks (74.09%) and 21-28 weeks (74.70%) (Figs. 3d and S3d), were subjected to 5hmC changes, and genes which were highly related to development were shown in Fig. 3d.
Additionally, we downloaded public DNase sequencing (DNase-Seq) data from the foetal heart, kidney and lung at 13-23 weeks 44 and found that there was a high correlation among the 5hmC signal, chromatin openness on promoters and gene expression level, for which genes could be detected from all three omics datasets (Fig. 3e). This nding suggests that 5hmC overlaps with open chromatin and is positively associated with transcriptional activity. However, we did not observe the above similar results between adjacent stages (Fig. S3e). Then, we wondered whether TFs regulate the process across different organs through speci c binding to 5hmC-enriched promoters. To discover corresponding 5hmC peaks based on known TF motifs, we performed pattern matching to reveal their one-to-one potential correspondence. To increase the reliability and speci city of the binding sequence, we used 30% of the single base proportion to obtain the optimal combinatorial sequences ( Fig. S3f, g). Four kinds of random peaks/fragments were used to eliminate the effects of background noise (paired Student's t-test, P value < 0.05) (Fig. S3h). According to the results, the number of 5hmC-modi ed target genes is different among distinct TF families. Cys2His2 zinc nger (C2H2 ZF), basic helix-loop-helix (bHLH), basic leucine zipper (bZIP), nuclear receptor, etc, are the main TF families in eukaryotes 45 , with more target genes than others in foetal organs (Fig. S3i). We further calculated the correlation of the 5hmC signal and the gene expression level of target genes from different TF families (Fig. 3f). The regulatory relationship of DNA 5hmC and transcription remained almost stable in C2H2 ZF, bHLH, bZIP, etc., whereas there were different regulatory roles in speci c TF families, such as adaptor protein 2 (AP-2), early 2 factor (E2F), and homeodomain ( Fig. 3f). All these ndings indicate that 5hmC-modi ed promoters are highly associated with gene expression homeostasis during foetal organ development.
RNA m 5 C acts as a downstream regulator during foetal organ development. In the setting of organspeci c DNA 5hmC in the regulation of upstream gene expression, we wondered how RNA m 5 C drives organ-speci c post-transcriptional regulatory networks, since it is one of the prevalent RNA posttranscriptional modi cations associated with RNA stability 31,32 . RNA bisul te sequencing (RNA-BisSeq) was performed on four foetal organs at three periods from 13 to 28 weeks. The samples had a dependable C-to-T conversion rate with good repeatability (Fig. S4a, b and Supplementary Table 5).
Consistent with the ndings in other vertebrates 28,29,31,32 , the majority of m 5 C sites were enriched near translation initiation sites (Fig. 4a). A total of 2,679-15,655 m 5 C sites within 316-1,622 mRNAs were identi ed at each stage during organ development (Fig. 4b). Eliminating the effect of continuous Cs, most m 5 C sites were strongly enriched in CG-rich and CDS regions (Fig. S4d, e). Considering that both the m 5 C methylation level and the number of m 5 C sites may affect the functions of corresponding mRNAs, we calculated the total methylation level across different organs (Figs. 4c and S4c). Globally, m 5 C levels dynamically change during the development of the kidney, liver and lung.
The global distributive features of m 5 C in human foetal organs prompted us to explore the potential role of m 5 C during development. Upon comparing RNA abundance changes between m 5 C-modi ed and unmodi ed mRNAs across adjacent stages, we found that the m 5 C-modi ed mRNAs were unstable before 18-19 weeks, whereas they became stabilized afterwards (P value < 5.40×10 − 7 ), indicating that m 5 C modi cation has a dual role in the regulation of organ development (Fig. 4d). This association was independent of coverage changes (Fig. S4f).
Additionally, consistent with that of DNA 5hmC, organ-speci c m 5 C-modi ed mRNAs were signi cantly enriched in the pathways related to cell differentiation, proliferation, and adhesion, and participated in the regulation of gene expression and translation during embryogenesis (Fig. S4g). Importantly, the organstage-speci c m 5 C-modi ed mRNAs were closely related to RNA abundance in each organ ( Fig. 4e and Supplementary Table 6). These genes revealed a high level of enrichment for organ physiology-related processes. Thus, these results indicate that as a post-transcriptional regulator, RNA m 5 C is indispensable during human foetal organ development.
DNA 5hmC and RNA m 5 C synergistically regulate organ-speci c developmental genes which are critical in maintaining normal foetal growth. We next proposed that DNA 5hmC and RNA m 5 C may cooperate to spatiotemporally modulate the regulatory networks in foetal organs. Overall, promoters of genes with high transcriptional activity presented relatively higher 5hmC signals than the m 5 C methylation level on their corresponding mRNAs, suggesting a dominant role of DNA 5hmC at the transcriptomic level during development (Fig. S5a). To further investigate the respective functions of DNA 5hmC and RNA m 5 C, we divided the organ-speci c modi ed genes into 4 groups: differentially expressed genes speci cally regulated by both DNA 5hmC and RNA m 5 C (group A), differentially expressed genes regulated only by DNA 5hmC (group B), differentially expressed genes regulated only by RNA m 5 C (group C), and commonly expressed genes speci cally regulated by DNA 5hmC and RNA m 5 C (group D) (Figs. 5a, S5b, c and Supplementary Table 7). As group D only accounts for 0.84% of the genes, these two types of modi cations indeed affect the expression of the vast majority of organ-speci c genes (Fig. S5b).
Speci cally, in group A, there were close associations among DNA 5hmC signal intensity, RNA m 5 C methylation level, and gene expression level (Fig. 5a).
Importantly, genes in all 3 groups (A, B and C) were involved in the process of organ formation and functional maturity ( Fig. 5b and Supplementary Table 8). In the foetal heart, the genes from all three groups were signi cantly enriched in heart reconstruction, cardiomyocyte differentiation, and heart raterelated cardiomyocytes. The co-regulatory genes were mainly enriched in angiogenesis. In contrast to DNA 5hmC mainly regulating cardiac myo bril assembly, related ion transportation, and cell communication through electrical coupling, m 5 C modi cation is involved in atrial action potential and action potential depolarization, especially for SA and AV node cells. Nevertheless, both ventricular action potential and action potential repolarization need the involvement of both modi cations. Although the foetal heart performs basic functions earlier than other organs, its gradual functional maturity and maintenance of normal development also require these two modi cations. Additionally, RNA m 5 C was even more critical in cell differentiation of kidney, such as the regulation of mesenchymal epithelial cell transformation during metanephric morphogenesis, glomerular development, and renal vesicle morphogenesis. DNA 5hmC is the key driver of organic acid metabolic process in the liver. In the lung, both DNA 5hmC and RNA m 5 C synergistically regulate respiratory gaseous exchange, and they separately modi ed different gene sets to modulate lung morphogenesis. These results suggest that DNA 5hmC and RNA m 5 C generally play a synergistic role in the process of human foetal development, but in the foetal kidney and liver they exhibit differential and speci c functions.
Interestingly, we found that 5hmC enrichment on promoters with m 5 C modi cation on their corresponding mRNAs was much higher than that without m 5 C (Fig. 5c), indicating a coordinated role of 5hmC and m 5 C in foetal organs. TET3 was reported to be the only existed DNA 5mC oxidase after fertilization, mediating the oxidation of most paternal nuclear cytosines and a large-scale removal of 5mC 23,46 . To determine the potential mechanisms underlying for the cross-talk between the two types of modi cations, we calculated the correlation coe cients between the expression of TET3 and RNA m 5 C regulators. Apparently, the correlation coe cients of the expression levels of TET3 and NOP2/Sun RNA methyltransferase 2 (NSUN2)/Aly/REF export factor (ALYREF) were high (0.7652 and 0.7664) (Fig. 5d). However, no signi cant relationship was found between TET3 and Y-box binding protein 1 (YBX1). These results suggest that DNA 5hmC and RNA m 5 C might drive organ-speci c regulatory networks through pre-and post-transcriptional mechanisms, respectively, which are mediated by TET3 and NSUN2/ALYREF. However, how they mechanistically contribute to developmental processes warrants further detailed investigation.

Discussion
The importance of epigenetic modi cations in embryogenesis is well known, but the global function and co-regulation by DNA 5hmC and RNA m 5 C at pre-and post-transcriptional levels during foetal organ development are currently poorly understood. Our study provides the rst overview of DNA 5hmC and RNA m 5 C during human foetal development, and reveals their cooperation as key molecular events in the development of foetal organs, including heart, kidney, liver and lung. We found that the key TFs, TBX20, PAX8, KLF1, TCF21 and CEBPB, speci cally contribute to the formation of distinct organs at different stages. Additionally, 5hmC-enriched Alu elements may participate in the regulation of expression of TF targeted genes. Both DNA 5hmC and RNA m 5 C dynamically modulate the changes in RNA abundance and are highly associated with the developmental functions of the corresponding organs. Moreover, this synergistic effect is more obvious at the later stages. The high correlation between TET3 and NSUN2/ALYREF strongly suggests the potential cross-talk between DNA 5hmC and RNA m 5 C. Collectively, the modi cations in the layers of DNA and RNA perform their functions and serve synergistic roles during foetal organ development (Fig. 6).
A previous study showed that the correlation coe cient between 5hmC and gene expression 2 kb upstream of the TSS is approximately 0.21 in the adult liver and lung 20 , suggesting its positive regulatory roles in human organs. In foetuses, we found that the relationship between the 5hmC signal on promoters and transcriptional activity was highly dynamic and that the 5hmC modi cation localized on open chromatin was strongly associated with the gene expression level (Fig. 3a, e). Chromatin accessibility of cis-regulatory sequences drives cell differentiation and zygotic genome activation (ZGA) by regulating gene expression in early mammalian embryos [7][8][9]13 . Studies have shown that post-translational modi cations, such as histone acetylation and ubiquitination, enhance chromatin accessibility by disrupting chromatin folding and globule formation 47,48 . Additionally, as a 5hmC binding protein, methyl-CpG binding protein 2 (MeCP2) competes with H1 for nucleosomal binding sites leading to chromatin higher-order structure changes 48-50 . According to the existing theoretical basis, we reasonably speculate that chromatin accessibility could be modulated by speci c binding of 5hmC-related regulators in a histone modi cation-dependent manner, resulting in high transcriptional activity.
RNA m 5 C modi cation has been recently demonstrated to be an important post-transcriptional regulator of mRNA stability during maternal-to-zygotic transition (MZT) in zebra sh and the pathogenesis of bladder cancer in humans 31,32 . In human foetuses, the involvement of RNA m 5 C was inferred to be divided into two periods with 18-19 weeks as the boundary. Compared with those in the former stages, m 5 C-modi ed mRNAs were more stable than unmodi ed mRNAs after 18-19 weeks in all four organs. Generally, the fate of modi ed-transcripts is determined by their regulators: "writers", "readers" and "erasers" 26 . In some cases, RNA methylation or its regulators play dual roles in cancer progression by affecting the stability of different target genes 33,51 , such as the m 6 A "writer" protein methyltransferase like 3 (METTL3) in colorectal cancer (CRC). We hypothesize that the differential roles of RNA m 5 C may result from the abundance/binding a nity changes of m 5 C regulators on different target genes to maintain normal development. In addition, the high correlation between the expression level of TET3 and NSUN2/ALYREF suggests that DNA methylase and RNA m 5 C methylase/binding proteins may synergistically regulate the early foetal organ development. A previous study demonstrated that H3K36me3 facilitates METTL3/METTL14/WTAP-dependent m 6 A deposition at the 3'end of coding sequences 34 . We speculate that TET3 not only mediates 5hmC formation, but also might interact with NSUN2/ALYREF, inducing mRNAs to obtain m 5 C modi cation and thus maintaining stability. However, the exact molecular mechanisms need to be further studied.
Taken together, this study not only uncovered dynamic DNA 5hmC and RNA m 5 C reprogramming during organogenesis in post-implanted embryos, but also revealed a putative cross-talk between DNA modi cation and RNA methylation, which adds another layer by which the transcriptome can be spatiotemporally co-regulated to form a well-coordinated network at speci c developmental stages ensuring normal foetal development. Furthermore, a comprehensive description of the characteristics and functional importance of modi cations at both the DNA and RNA levels provides a foundation for understanding the in-depth epigenetic mechanisms for early human development.

Methods
Tissue collection and ethics approval of this study. All of the human foetal organs at different developmental stages were obtained from the donors voluntarily at the First A liated Hospital of Zhengzhou University with signed informed consent and the approval of the Ethics Committee (license no. 2020-KY-261). Foetuses were at 13 to 28 weeks in the second trimester of pregnancy. Foetal heart, liver, kidney, and lung were collected from selective induced labor fetus of donors with normal chromosome karyotypes and no family history of infectious diseases, hereditary diseases and smoking history. The foetal organs were isolated after a mechanical dissection, serially washed with medical saline for three times, stored in RNAlater Stabilization Reagent (AM7024, Invitrogen) and immediately frozen in liquid nitrogen. The details of this research, including the research procedure, the bene ts and risks of donation, the con dentiality of the donor's information, the right of withdrawal, the generation and analysis of sequencing data, and the publication of results, were explained to the donors before they gave consent. All the experiments were also under the regulations of International Ethical Guidelines for Biomedical Research Involving Human Subjects and in compliance with the policies and laws of China.
Genomic DNA and total RNA preparation. A maximum amount of 20 mg of tissue ground in a 100:1 mixture of 600 µl Buffer RLT Plus and 6 µl 14.3M β-mercaptoethanol for genomic DNA and total RNA extraction by user operation manual of AllPrep DNA/RNA/miRNA Universal Kit (80224, Qiagen). For rRNA depletion, the puri ed total RNA was further treated by a 3:1 mixture of probe and total RNA followed by RNase H (EN0202, Thermo Fisher) treatment at 37°C. TURBO DNase (AM2238, Invitrogen) treatment was used to eliminate DNA contamination. Then, the rRNA-depleted RNAs were puri ed by ethanol precipitation and used for construction of the RNA-seq and RNA-BisSeq libraries. The concentration and quality of puri ed genomic DNA and rRNA-depleted RNA were detected by NanoDrop (Thermo Fisher) and Qubit (Life Technologies).
RNA-seq library generation and sequencing. The RNA-seq libraries of human foetal organs were constructed by using user manual of KAPA RNA HyperPrep Kit (KK8541, KAPA). In brief, the high-quality, rRNA-depleted RNAs of foetal organs were fragmented for 6 min at 85℃ (about 300 nucleotides) by 2× Fragment, Prime and Elute Buffer in different periods. The DNA concentration and quality of constructed libraries were detected by Qubit (Life Technologies) and 2100 Bioanalyzer (Agilent) before sequencing.
Bisul te conversion of RNA. RNA fragmentation and bisul te conversion were consistent with previous description 28 , but had some modi cations. Brie y, in vitro transcribed mouse Dhfr mRNA, as a methylation conversion reference, was mixed with puri ed rRNA removal RNA at a ratio of 1:200. The mixture was fragmented into ~ 200 nucleotides by 10× RNA Fragmentation Reagent (AM8740, Thermo Fisher) for 22s at 88℃ which was terminated by 10× RNA stop solution. Then fragmented mRNAs were obtained by ethanol precipitation. After washing precipitation with 100% ethanol, the RNA precipitates were dissolved in 100 µl bisul te solution that mixed 40% sodium bisul te (243973,Sigma) and 600 µM hydroquinone (H9003 Sigma) in a ratio of 100:1. The synthesis was incubated at 75 ℃ for 4 hours. Nanosep columns with 3K Omega membranes (OD003C35, Pall Corporation) were used to desalt the reaction mixture with centrifugation. The RNA pellet was washed with 1 M Tris-HCl (pH 9.0) followed by centrifugation for ve times. Finally, the RNA was re-suspended in 75µl of nuclease-free water and incubated at 75°C for 1 hour with equal volume of 1 M Tris-HCl (pH 9.0) for desulfonation. After ethanol precipitation The bisul te-converted RNA was dissolved in RNase-free water and prepared for library construction of RNA-BisSEq. After reverse transcription with SuperScript II Reverse Transcriptase (18064014, Invitrogen) and ACT random hexamers to form cDNA, KAPA RNA HyperPrep Kit (KK8541, KAPA) was used to perform the subsequent procedures as described in the manufacturer's instructions. Paired-end sequencing was performed on the Illumina HiSeq2500 instrument with 125 bp read length.
Nano-hmC-Seal in human tissues. Genomic DNA was prepared from the above steps of genomic DNA and total RNA preparation. The input DNA ranged from 5 ng to 50 ng. Nano-hmC-Seal library construction was similar to the description of previous study 36  The mapped reads of public DNase-Seq data from foetal heart, kidney and lung were downloaded and used to do peak calling by DFilter (version 1.6) 63 with the parameters: -ks = 50 -lpval = 2 -f = bam.
Differential 5hmC regions (DhMRs) identi cation. A non-overlapping 100bp window was applied to perform DhMRs identi cation across the samples from different organs. Brie y, the genome was rst binned into consecutive 100bp windows, and the signal intensity of each window was calculated based on the normalized RPM. Then the average 5hmC signal of replicates were calculated. ANOVA (Analysis of Variance) was used to compute the magnitude of the difference for each bin in R (version 3.4.1) (https://www.r-project.org/). The P values were calculated using Benjamini-Hochberg method. The bins with fold change > 1.5 and P value < 0.05 were considered statistically signi cant. The 5% least signi cant bins were considered as common bins (control). The signi cantly differential bins and common bins were merged respectively according to their positions on the chromosome and were Identi cation of clusters with temporal changes. The clusters with sequential changes were obtained by K-means with the method Pearson correlation (MEV, version 4.9.0) (http://mev.tm4.org). When the number of clusters was greater than 6, similar 5hmC trends appeared. Thus, the DhMRs was divided into 6 clusters nally. HOMER (version 4.9.1) 56 was used to performed TF prediction in DhMRs of each clusters with the following command: ndMotifsGenome.pl input.bed hg19 -size given -len 6 -norevoppcache 1000. Genes with TF-potentially-binding DhMRs only from one clusters are considered as singleregulation groups, whereas genes with TF-potentially-binding DhMRs from two or more clusters are considered as mixed-regulation groups.
Differential The results were further ltered to remove binding sites of TFs that were scarcely expressed at all 4 stages in foetal kidney (RPKM < 0.1). In order to eliminate ambiguity in the correlation between TF binding sites, differential 5hmC modi cations and downstream gene expression, Alu elements that were bound by only one TF were kept for further analysis. The resulting Alu elements were then ltered by their proximity to downstream genes, keeping those within 1000 bp upstream of a gene with proper orientation and observable expression (RPKM > 0.1) in at least one of the four stages. To obtain the enrichment scores of binding sites for CEBPB and TRPS1 with these Alu DhMRs, we rst calculated the density of the intersection of either TF binding regions with these Alu DhMRs. This density was then divided by either the density of the intersection of these TF binding sites with all Alu elements in the human genome, or by the density of these TF binding sites in the whole genome. The enrichment scores of global TF binding sites in all Alu elements were calculated by the ratio of the density of the intersection of corresponding TF binding sites with all Alu elements in genome divided by the density of corresponding TF binding sites in the whole genome.
Dynamic changes of DNA 5hmC between adjacent stages. The promoters for which the fold change of 5hmC signal between samples from two adjacent stages exhibited a greater than 25% were de ned as signi cant different promoters (Student's t-test with P value < 0.05). Thus, the promoters were separated into two group sets: promoters with decreased 5hmC signal and promoters with increased 5hmC signal.
To explore the relationship between 5hmC signals in promoters and gene expression, the gene expression levels of corresponding group sets were compared and the P values were calculated using two-sided Wilcoxon and Mann-Whitney test. The genes with fold change over 1.2 were de ned as differential expressed genes in 5hmC-decreased and -increased groups.
Corresponding 5hmC peak discovery based on known TF motifs. To further explore how different TFs regulate gene expression via binding to 5hmC-enriched promoters, TFs-binded 5hmC peaks were identi ed based on the known motifs of TFs by pattern matching. The motifs of TFs were downloaded from HumanTFs database 45 (http://humantfs.ccbr.utoronto.ca/) which contains ~ 1600 TFs' information (~ 5000 possible motifs). Different TFs could bind to the sequence with the same motif and one TF could recognize the sequence with different motifs. In addition, there are different kinds of combinatorial sequences for each motif. Thus, we exhausted all the possible combinatorial sequences of 10 nt motifs for each TF. To increase the reliability of the sequence, 0.1-0.9 gradients of threshold for single base proportion were set to obtain the optimal combinatorial sequences. The single base proportion represents the probability of a speci c base (A, T, C and G) at a speci c position. 5hmC peaks on promoters were subjected to cyclic matching to those combinatorial sequences of motifs. Four kind of backgrounds were used to eliminate the effects of background noise: 1) random 5hmC peaks; 2) random 5hmC peaks in distal region; 3) random genome fragments in the same length of that of the average of 5hmC peaks; 4) random genome fragments in distal region in the same length of that of the average of 5hmC peaks. The fragments were randomly selected for 20 times per background. The motif sequence remained if the number of actual peaks divided by the number of fragments in the background > 2 for at least one random fragments group. The P values were calculated using Paired Student t-test. The remaining motifs that could be detected based on all the four kinds of backgrounds were considered credible.
Dynamic changes of RNA m 5 C. All expressed genes were divided into two groups based on with or without m 5 C modi cation at adjacent time points. To further evaluate the effects of m 5 C dynamic changes on mRNA abundance change, mRNAs with m 5 C modi cation at last stage but without m 5 C modi cation at next stage were de ned as mRNAs with m 5 C loss; mRNAs without m 5 C modi cation at last stage but with m 5 C modi cation at next stage were de ned as mRNAs with m 5 C gain. The mRNA abundance changes of these two groups were compared. The P values were calculated using two-sided Wilcoxon and Mann-Whitney test. The organ-speci c m 5 C-modi ed genes were identi ed according to the DhMRs methods.
Association analysis with different organ-speci c groups. To analyze the relationship between DNA 5hmC, RNA m 5 C and transcriptome, the 5hmC signal of promoter, m 5 C methylation level and gene expression level (RPKM) were subjected to log 2 transformation and deviation standardization [(xmin)/(max-min)] (x, input speci c value; min, the minimum value of one speci c omics data; max, the maximum value of one speci c omics data). The organ-speci c expressed genes were identi ed using ANOVA in R (version 3.4.1). To further explore how epigenetic marks differentially regulate foetal organ development, organ-speci c genes were separated into 4 groups: differentially expressed genes that DNA 5hmC and RNA m 5 C speci cally regulated, differentially expressed genes that were only speci cally regulated by DNA 5hmC, differentially expressed genes that were only speci cally regulated by RNA m 5 C, and commonly expressed genes that DNA 5hmC and RNA m 5 C speci cally regulated.

Figure 2
Page 23/29 DNA 5hmC dynamically shapes foetal organ development by transcription factors (TFs) as well as transposable elements (TEs). a, K-means clustering analysis of organ-speci c DhMRs from heart, kidney, liver and lung. Each organ was separated into 6 clusters. b, Heatmap showing the TF motifs identi ed by DhMRs. The TF-binding signi cance of single-and mis-regulation gene groups, and the expression levels of each TF are shown from left to right. c, Enrichment of speci c TE families in organ-stage-speci c DhMRs of foetal kidney. d, Bar plots showing the enrichment score of CEBPB binding sites in Alu elements in cluster 3 and cluster 4 compared with their overall enrichment in all Alu elements (red) and their genome-wide distribution (blue). The overall enrichment of all CEBPB binding sites in all Alu elements compared with their genome-wide distribution is also shown (gray). e. The z-score normalized 5hmC signal of CEBPB-targeted Alu elements as in d and the expression level of downstream genes during human foetal kidney development expression level of 18 organ-speci c genes from 5hmC decrease and increase groups. e, Biplots comparing changes in 5hmC signals (x axis), changes in DNase I hypersensitive signals (y axis) with the RNA abundance changes (color) between the adjacent stages. Each dot represents a single gene. Public DNase-Seq data were from GSE18927. f, Dynamic of pearson's correlation coe cient of different TF family groups from heart, kidney, liver and lung. Pearson's correlation coe cients were calculated based on 5hmC signals and gene expression level in corresponding gene sets.

Figure 4
Page 26/29 The RNA m5C pro le during foetal organ development. a, The m5C distributions within different regions in different samples from heart, kidney, liver and lung. b, Bar charts show the numbers of m5C sites and m5C-modi ed mRNAs. c, Overall total RNA m5C methylation level across different developmental stages in heart, kidney, liver and lung. d, Cumulative distribution displaying the expression level change of m5Cmodi ed and unmodi ed mRNAs from RNA-Seq data comparing samples from adjacent stages. The P values in c and d were calculated using two-sided Wilcoxon and Mann-Whitney tests. e, Heatmap showing the strong positive correlation between the m5C methylation level of organ-stage-speci c genes and the expression level of corresponding genes during foetal organ development. GO biological processes for each organ-stage-speci c gene set are listed on the right. biological process. c, The 5hmC signal on promoters of corresponding genes with or without m5C modi cation. The P values were calculated using two-sided Wilcoxon and Mann-Whitney tests. d, Scatterplots showing the correlations between the expression levels of 5hmC methyltransferase TET3 and m5C "writer" NSUN2, "reader" ALYREF and YBX1 in foetal organs. Figure 6