{'fraction_non_alphanumeric': 0.029189222440944882, 'fraction_numerical': 0.01879306102362205, 'mean_word_length': 5.128746465598492, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The promoter methylomes of monochorionic twin placentas reveal intrauterine growth restriction-specific variations in the methylation patterns

Intrauterine growth restriction (IUGR) affects the foetus and has a number of pathological consequences throughout life. Recent work has indicated that variations in DNA methylation might cause placental dysfunction, which may be associated with adverse pregnancy complications. Here, we investigated the promoter methylomes of placental shares from seven monochorionic (MC) twins with selective intrauterine growth restriction (sIUGR) using the healthy twin as an ideal control. Our work demonstrated that the IUGR placental shares harboured a distinct DNA hypomethylation pattern and that the methylation variations preferentially occurred in CpG island shores or non-CpG island promoters. The differentially methylated promoters could significantly separate the IUGR placental shares from the healthy ones. Ultra‐performance liquid chromatography/tandem mass spectrometry (UPLC‐MS/MS) further confirmed the genome‐wide DNA hypomethylation and the lower level of hydroxymethylation statuses in the IUGR placental shares. The methylation variations of the LRAT and SLC19A1 promoters, which are involved in vitamin A metabolism and folate transportation, respectively, and the EFS promoter were further validated in an additional 12 pairs of MC twins with sIUGR. Although the expressions of LRAT, SLC19A1 and EFS were not affected, we still speculated that DNA methylation and hydroxymethylation might serve a functional role during in utero foetal development.


Results
Sample characteristics and promoter LHC-BS sequencing. In this study, seven pairs of placental shares from MC twins, with one that was diagnosed as sIUGR, were enrolled for LHC-BS sequencing and mass spectrometry analysis. In addition, another twelve pairs of placental samples from MC twins with sIUGR were employed to validate the differentially methylated promoters. The subjects' clinical characteristics were summarized in Table S1.
The promoter methylomes of the healthy and IUGR placental shares were obtained by LHC-BS sequencing, respectively 26,27 . Promoter-targeted LHC-BS could interrogate 1,860,932 CpG sites, of which 138,171 and 5,387 were located in the X and Y chromosomes, respectively, by designing probes covering the upstream 2,200 bp and downstream 500 bp across the transcriptional start sites (TSS) of the Crick strand 27 . A total of 66.6% of the targeted promoters contain CGI(s) and 69.0% of the associated CGI shores are covered as well. Thus, 47.7% of the total CGIs and 35.8% CGI shores in the genome are included in the probes.
From the sequencing results, we generated an average of 5.73 Gb of clean sequencing data for each sample and reached an average of 21.39X coverage depth for the CpGs (Table S2). A total of 86.26% of the generated reads were mapped to at least one genomic position and 74.51% were uniquely mapped. To evaluate the capture efficiency, we found that approximately 97% of the target regions were detected. Furthermore, on average, 91.4% of the designed CpG sites in the Watson strand were detected by at least one depth and 70.3% were detected by at least five depths (Table S2). Although the capture specificity was just approximately 40%, most of the reads that missed the target in the Watson strand were mapped to their complementary strand (the Crick strand) and some reads were flanked on the designed promoter regions. CpG sites with less than 5X sequencing depth and reads mapped to the Crick strand were filtered out, and, as a result, we obtained an average of approximately 1.31 million CpG dinucleotides for each sample for the subsequent analyses (Table S2).
Global DNA hypomethylation occurred in the placental shares of the IUGR foetus. First, we evaluated the whole promoter methylation levels of the seven pairs of samples and found that the methylation levels of the IUGR placental shares were significantly lower compared to those of the matched healthy twin (Fig. 1A). To further estimate the distribution of promoter DNA methylation between the IUGR twins' and co-twins' placental shares, we divided the promoters into three categories based on the CpG ratio and the length of the CpG-rich regions 29 , namely high-CpG promoters (HCPs), intermediate-CpG promoters (ICPs) and low-CpG promoters (LCPs). As a result, we found that all three kinds of promoters showed reduced DNA methylation in the IUGR placental shares and the variation was more likely to have occurred in the ICP and LCP promoters ( Figure S1). In agreement with our result of the tendency for whole promoter DNA hypomethylation in the IUGR placental shares, wide-spread gene enhancer demethylation has been recently reported in placental tissues of early-onset pre-eclampsia 30 , which is a common cause of pre-term birth and foetal growth restriction. 5-hydroxymethylcytosine (5 hmC) is a recently discovered DNA modification in mammalian cells that is catalysed by ten-eleven translocation (TET) dioxygenases using the 5-methylcytosine (5 mC) substrate 31,32 . However, bisulfite sequencing could not distinguish 5 hmC from 5 mC 33 . Genomic profiling has revealed that 5 hmC was particularly enriched in the 5 mC-depleted regions, which was thus believed to be involved in the positive DNA demethylation process. 5 hmC is abundant in embryonic stem cells and the brain. However, it is reduced in cancer tissue and is frequently associated with gene and gene regulatory elements 32,[34][35][36] , which further indicated that it may serve as more than simply a transient intermediate during the DNA demethylation process and may function in human development and disease.
To further investigate if there is a genome-wide DNA demethylation event in the placental shares of the IUGR foetus and to estimate if 5 hmC plays a functional role in the placental tissue during foetal development in utero, we conducted UPLC-MS/MS experiments on the seven MC twin pairs' placentas. Consistent with our promoter methylome results, we confirmed that the IUGR placental shares were significantly hypomethylated across the genome (Fig. 1B). Interestingly, in contrast with a previous report that 5 hmC was abundant in 5 mC-depleted regions, the UPLC-MS/MS results showed that the 5 hmC levels were also decreased in the IUGR placental shares (Fig. 1C). Pair-wise comparisons between the placental shares of the IUGR and normal foetus within one MC twin further confirmed this conclusion ( Figure S2). Based on these findings, we hypothesized that in addition to placental DNA methylation, the hydroxymethylation modification of the genomic DNA might also be important for in utero foetal development. Thus, profiling their genome-wide distribution in future studies might further assist us in understanding the involvement of epigenomics in foetal development and associated diseases.
The promoter methylomes separated the placental shares of the IUGR foetuses from their co-twins. MC twins are derived from one zygote by splitting into two separate embryos, and they share the same placenta. The DNA methylation statuses were thus expected to be similar between the two individuals or their shared placental tissues 37 . To estimate the DNA methylation variations between the placental shares of the healthy and IUGR foetuses, we selected the commonly covered CpG sites (approximately 1.02 M) in all seven pairs of samples and performed principal component analysis (PCA) based on the methylation level of each site. Interestingly, we found that the IUGR placental tissues were significantly separated from the normal tissues ( Fig. 2A). Because the variations may have occurred in the ICP and LCP regions ( Figure S1), we further performed a hierarchical clustering analysis on the average methylation levels of the CGI (Fig. 2B) and non-CGI promoters (Fig. 2C), respectively, using the "pvclust" algorithm 38 . In agreement with the PCA results, the methylation levels in both the CGI and non-CGI promoters showed a similar tendency in differentiating the IUGR placental shares from the healthy placental shares, which indicated that the IUGR placental shares harboured a distinct promoter methylation pattern.
To further evaluate the distinct pattern of promoter methylomes between the IUGR and healthy placental shares, we conducted the Mann and Whitney test and identified 4,605 significantly differentially methylated sites (DMSs) with a P-value < 0.05 and a methylation difference greater than 20%. Unsupervised clustering of the placentas using DMSs could perfectly discriminate the tissues according to their pathological characteristics (Fig. 2D). In agreement with the finding that the CGI regions are more refractory to methylation variations 37,39 , we confirmed that the IUGR placental DNA demethylation mostly occurred in the non-CGI promoters and CGI shores (Fig. 2E). Although a fraction of CpG sites in the CGI promoter displayed variations in methylation, most of the sites became hypermethylated, which is in contrast to the genome-wide demethylation in the IUGR placental shares (Fig. 2E). CGI promoter hypermethylation is frequently associated with transcriptional repression in cancer 22 , which sheds light on its functional roles in IUGR pathogenesis. Taken together, these findings indicated that the IUGR placental shares harboured a distinct promoter methylation characteristic, which might serve as a novel biomarker for monitoring the foetal intrauterine growth status.
The differentially methylated genes are involved in multiple pathways in the maternal placentas. To minimize the inter-group methylation variations, we performed pair-wise comparisons on the promoter methylation patterns to identify the differentially methylated regions (DMRs). Applying a sliding window strategy with Fisher's exact test, we identified a total of 2,930 promoters that contained at least one DMR, based on the criteria described in the Methods (Table S3). To estimate the inter-group methylation variations, we further calculated the DMR distributions within each pair of samples. As a result, the data showed that with similar DNA methylation variations, the seven pairs of placentas contained nearly 600 DMRs, respectively ( Figure S3). LOC647288 is a pseudogene that is hypomethylated in all seven IUGR placental shares. Moreover, lecithin retinol acyltransferase (LRAT) is a major enzyme involved in vitamin A metabolism by catalysing the esterification of all-trans-retinol into all-trans-retinyl ester, which regulates various physiological processes. However, it is frequently inhibited by promoter hypermethylation in numerous cancers [40][41][42] . In this study, we showed that the LRAT promoter (GRCh37/Hg19 chr4: 155663541-155664999) was hypermethylated in six IUGR placental shares and the inter-group comparisons displayed a similar result (Fig. 3A), which implied that aberrant DNA methylation might affect in utero foetal development through vitamin A metabolism.
To further estimate the functional effects of the DNA methylation variations, we screened 128 promoters that were differentially methylated in no less than three pairs of IUGR placental shares and performed KEGG pathway enrichment and gene ontology (GO) analysis using WebGeStalt 43 . Notably, eight KEGG pathways were significantly enriched with an adjusted p-value less than 0.05 (Table S4). In addition to LRAT, SLC19A1, whose promoter was hypomethylated in four pairs of IUGR placental shares, was also involved in vitamin metabolism pathways. An inter-group comparison confirmed the differential methylation of SLC19A1 (GRCh37/Hg19 chr21: 46963151-46964484) (Fig. 3B). It has been reported that disturbances in maternal vitamin absorption are connected with restrictions of foetal growth and development and could increase the risk of chronic disease in later life 44 . Furthermore, multiple genes that were involved in the protein digestion and absorption pathway, including KCNN4, COL12A1 and PRSS3, were also aberrantly methylated in most IUGR placentas. Other enriched pathways included Dorso-ventral axis formation, Axon guidance, Maturity onset diabetes of the young, Amyotrophic lateral sclerosis (ALS), Olfactory transduction and Neuroactive ligand-receptor interaction. Gene ontology analysis based on the molecular function and biological processes revealed that the differentially methylated genes were mainly enriched in DNA binding activity and the organismal development process (Fig. 3C, Table S5).
Comparisons of the promoter methylation of imprinted genes between the IUGR and healthy placental shares. Imprinted genes, which have been thought to function in the regulation of embryonic and placental development as well as in intrauterine foetal growth conditions, are theoretically controlled by DNA methylation in imprinting control regions (ICRs) 2,13 . However, their promoter methylation characteristics in IUGR and normal placental tissues have not been studied. Furthermore, it has been reported that some of the allele-specific expressed imprinted genes were even controlled by the promoter DNA methylation statuses, including Gsα, XLsα, A/B, NESP55 and GNAS-AS1 45 . To evaluate if the promoters of imprinted genes are susceptible to methylation changes in response to IUGR, we collected the currently known human imprinted genes from the geneimprint database (www.geneimprint.com) and divided them into two groups based on their paternal allele-specific or maternal allele-specific expression patterns. As a result, we found that the mean methylation levels of the promoters for both maternally and paternally imprinted genes attended to decrease in the IUGR placental shares (Fig. 3D).
The methylation variations of the EFS, SLC19A1 and LRAT promoters were validated in MC twins with sIUGR. To further identify the aberrantly methylated genes in the IUGR placental shares, we chose four genes, LRAT, SLC19A1, EFS and SR140, and validated their promoter methylation variations in an additional twelve pairs of MC twins with sIUGR (Table S1). The regional information for validation and the associated primers were summarized in Table S6. Bisulfite genomic sequencing with the multiple clones approach (BGS) was employed, and twelve to twenty subclones for candidate genes were sequenced on a 3730 Genetic Analyzer (Applied Biosystems), respectively. No less than five subclones of a candidate gene were successfully sequenced in one sample. In agreement with the LHC-BS sequencing results, the inter-group comparisons showed that the methylation levels of the EFS, SLC19A and LRAT promoters were significantly changed in the IUGR placental shares (Fig. 4A). Consistently, pair-wise analyses of the methylation statuses of the CpG sites for each validation region in the IUGR and normal placental shares also showed similar results (Fig. 4B,D).
The protein encoded by EFS plays an important role in intracellular signal transduction 46 . Aberrant promoter methylation could significantly inhibit EFS expression and is correlated with poor prognosis in cancer 47 . To further explore the correlations between the promoter methylation variations and gene expression in the IUGR placental shares, we concomitantly tested the relative expressions of EFS, LRAT and SLC19A1 in eleven pairs of placental shares using RT-QPCR (Table S1). However, they were not significantly changed by normalizing to the endogenous housekeeping gene of GAPDH, based on the collected tissues (Fig. 5). GAPDH has been widely used as a standard during testing gene expressions 14,15,48 , but was reported to be not stably expressed in placental tissues 49 . Thus, we selected another four genes, SDHA, TBP, YWHAZ and RPL19, whose expression levels were confirmed to be stable in placental tissues 49,50 , as four internal controls to validate the relative expressions. The geometric mean of the four carefully selected housekeeping genes were used for normalization as previously described 51 . However, we still did not find a significant correlation between the promoter methylation variations and gene expression ( Figure S4). Further studies to test the relationship between promoter DNA methylation and the gene expression of these candidate genes in homogeneous cells from placentas and further expanding the number of samples in future studies might help us to understand their contributions to IUGR pathogenesis.

Discussion
Numerous studies have begun to evaluate the correlation between IUGR and the epigenomic modifications in foetal tissues, including the placenta 2,8,12,13 . The placenta plays a central role in foetal development by providing the foetus with nutrients and removing waste from the foetal circulation throughout pregnancy. However, it is subject to different environmental factors, each with the potential to alter the placental epigenomic profiles, which might then contribute to placental dysfunction and adverse pregnancy complications 12 . Many of the current studies were focused on genomic imprinting or identifying the imprinted genes that were speculated to be involved with the occurrence of IUGR [14][15][16][17] . The genome-wide DNA methylation profiles of human placentas were recently investigated using an Illumina Infinium HumanMethylation27 BeadChip array, which confirmed that placental DNA methylation was significantly associated with foetal growth and could function as a marker for the intrauterine environment 18 . This finding provided new insight into epigenomic modifications other than genomic imprinting that are involved in the occurrence of foetal growth restriction in utero, as presented here.
We hypothesized that alterations of DNA methylation in promoters, which are recognized as an essential regulatory element for gene expression, might be found in the placentas of neonates and significantly associated with foetal growth in utero. LHC-BS is our recently developed method, with single base-pair resolution analysis of the targeted regional DNA methylation status 26,27 . Here, we applied the newly developed promoter-targeted LHC-BS approach to investigate the promoter methylome of seven pairs of MC twins' placentas and further identify the genes that were regulated by promoter methylation but not genomic imprinting throughout the genome. In agreement with Banister's finding that the IUGR-associated placentas harboured distinct DNA methylation patterns 18 , our result showed that the promoter methylome could significantly distinguish the IUGR placental shares from the normal placental shares (Fig. 2). Furthermore, we found that most of the variations occurred in non-CGI promoters or CGI shores (Fig. 2D,E), which were also preferentially associated with tissue-specific differential methylation regions (T-DMRs) and cancer-specific differential methylation regions (C-DMRs) 52,53 .
Previous research demonstrated that the levels of human placental DNA methylation were approximately 2.5-3.5%, which was lower than most somatic tissue (3.5-5%) and human tumours (2-4%) 8,54 . However, here, we performed UPLC-MS/MS analyses on the seven pairs of MC twin placentas with sIUGR and showed that the genomic DNA methylation levels in the placental shares from the normal foetuses were lower than 2% (Fig. 1B). This discordance may primarily be due to differences in the techniques used and the regions selected for detection in the previous study. Most interestingly, we simultaneously found that the IUGR placental shares were methylated at much lower levels compared to the normal shares (Fig. 1B). However, the hydroxymethylation level was not increased as anticipated, but was also decreased (Fig. 1C). There is a step-wise genomic demethylation event from fertilization until the morula stage in mammals. After that, DNA methylation is then re-established in a tissue-specific manner during differentiation, which is important for placental morphology and physiology 9 . The inner cell mass, which gives rise to all of the tissues of the adult and some extraembryonic tissues, is hypermethylated, while the trophectoderm, which forms the placenta, becomes hypomethylated 2,8 . It is not clear whether the much lower methylation levels in the IUGR-associated placentas were formed during the de novo DNA methylation process or affected by environmental factors, but they are very likely to affect placental functions, as revealed by our KEGG pathway enrichment and gene ontology analyses.
One limitation of this study is that the revealed promoter methylome was commixed with the DNA hydroxymethylation modification. Similar to other bisulfite conversion-based methods, LHC-BS could not distinguish 5-hydroxymethylcytosine (5 hmC) from 5-methylcytosine (5 mC) 27,33 . Previous studies demonstrated that the 5 hmC levels varied in tissue-specific and cell-type-specific manners 32,34,36 and were significantly associated with human disease, including cancer 34,35 . Integrated, single base-pair resolution analyses of the distributions of DNA methylation and hydroxymethylation in IUGR and matched placentas will be needed to fully identify the epigenomic contributions to in utero foetal development.
The advent of examining the circulating foetal DNA in the maternal peripheral blood provided an optimal approach for non-invasive prenatal diagnosis. In maternal plasma, 3 to 6% of the cell-free DNA during pregnancy is derived from the foetus, and the placenta is the predominant source 55 . Furthermore, the placenta contains different types of cells in both the trophoblast layer and stromal bed 8 . Placental inhomogeneity increases the difficulties of sample collection and might affect the detection accuracy in related studies. Directly identifying the methylation variations in placental tissues could help us to understand the underlying mechanisms, but is hard to apply to prenatal diagnosis. Therefore, it will be important to test if there are methylation variations in the maternal peripheral blood and validate the identified promoters in the cell-free circulating foetal DNA in future studies for applications in non-invasive prenatal diagnosis and foetal development monitoring.
Although we confirmed the methylation variations in the EFS, LRAT and SLC19A1 promoters by both LHC-BS and GBS approaches, their expression levels were not significantly changed between the IUGR and healthy placental shares. One possible reason for this observation may be related to the stability of the mRNA by post-transcriptional mechanism to control gene expression 25 , and the placental inhomogeneity may also affect the detection result. Furthermore, numerous studies have indicated that variation in DNA methylation may not necessarily relate to variation in gene expression, and their relationship is more complex than anticipated 56 . The reason for this is unclear, but it might be related to multiple additional layers of epigenetic modifications cooperating with DNA methylation to regulate gene expression.
In summary, we explored the promoter methylomes and the genome-wide DNA methylation/hydroxymethylation levels in the placentas of MC twins with sIUGR to study the contribution of placental DNA methylation towards foetal development. We revealed that the placental shares of the IUGR foetuses were significantly hypomethylated in both the promoter and throughout the genome and that the methylation profile of the promoters was significantly connected with the foetal IUGR status. In addition, we found that some of the aberrantly methylated genes were involved in multiple metabolic pathways, such as vitamin/protein digestion and absorption during foetal development. Therefore, further estimations of these variations in a homogeneous type of cells from placental tissue with an expanded sample number might help us to understand the contributions of the placenta to IUGR pathogenesis.

Materials and Methods
Collection of placenta samples. We collected the placental tissues of 19 pairs of MC twins with sIUGR between May 2012 and April 2013. The study was performed with the approval of the ethics committee of the First-Affiliated Hospital of Sun Yat-sen University and the ethics committee follows the World Medical Association's (WMA) Declaration of Helsinki-Ethical Principles for Medical Research Involving Human Subjects. The collection of placental tissues was performed with the patients' informed consent. The chorionicity was determined using an ultrasound at 11-14 weeks as previously described 57 , and further confirmed by placental examination after delivery. The sIUGR is used to define cases when one twin's birth weight is less than the 10th percentile for the gestational age, as judged by the previous reference for the MC twins birth weight 4 . The pregnancies with severe maternal complications, TTTS and foetal death were excluded. The sIUGR is classified by the Doppler pattern of umbilical artery diastolic flows of the growth restricted twin 58 . sIUGR type I has a normal diastolic flow. Type II is defined by persistently absent/reverse end-diastolic flow, while type III is defined by the presence of intermittent absent/reverse end-diastolic flow (iAREDF).
The placentas were obtained after vaginal delivery or caesarean section. The placental tissue around the individual insertion region of the umbilical cord was collected within 30 minutes after delivery. The placental tissue was excised from inside the placental lobules, avoiding both the maternal surface and the amniotic membrane. The tissue samples were excised and washed 3 times in sterilized, ice-cold PBS to eliminate any blood. The tissues were placed in freezing tubes, deep frozen in liquid nitrogen overnight, and then transferred to a − 80 °C freezer for storage.

DNA extraction and construction of the LHC-BS libraries.
The genomic DNA was extracted using the QIAamp DNA Blood Mini Kit (Qiagen) according to the manufacturer's instructions. After purification, the concentration and integrity of the isolated genomic DNA were analysed on a Qubit fluorometer using a Quant-iT dsDNA HS Assay kit (Invitrogen) and 1% agarose gel electrophoresis, respectively. The promoter-targeted bisulfite conversion libraries were constructed as previously described 27 , with some modifications. Briefly, 1 μ g of the genomic DNA was fragmented by an ultrasonic instrument (Covarias) to a mean size of 250 bp, followed by end-repairing, dA addition to the 3′ end and barcode adapter ligation. Four adapter-ligated samples were equally pooled together for the subsequent hybridization step using our previously designed probes 27 ; each sample was 250 ng. The specifically captured DNA ligated with adapters was further treated using an EZ DNA Methylation-Gold Kit (ZYMO), followed by purification using a MiniElute PCR purification Kit (Qiagen). After PCR amplification, the products were purified using a 0.8-fold dilution of AMPure beads (Agencourt), and then quantified using the Bioanalyzer analysis system (Agilent) and real time PCR assay. The qualified libraries were finally analysed using an Illumina HiSeq2000.

UPLC-MS/MS analysis of the global DNA methylation and hydroxymethylation.
One microgram of genomic DNA was incubated at 37 °C for 6 h in 150 μ l of digestion buffer (6.7 mM MgCl2, 33 mM NaCl and 6.7 mM Tris-HCl pH = 7.9) containing 2.5 U benzonase nuclease, 3 mU phosphodiesterase I from crotalus adamanteus venom and 2 U alkaline phosphatase, calf intestinal (CIP). A microcon centrifugal filter device with a 3,000 D cut-off membrane was used to the remove protein from the digested DNA samples by centrifuging them at 12,000 rpm for 60 min. The mobile phase consisted of 0.1% formic acid (solvent A) and methanol containing 0.1% formic acid (solvent B). The flow rate was set to 500 μ L/min. The enzymatically digested DNA samples (5 μ L each) were injected for UPLC-MS/MS analysis with the following LC gradient: 0.0-4.0 min, 0 to 50% solvent B; 4.0-6.0 min, 50% solvent B; 6.0-6.1 min, 50% to 5% solvent B; and 6.1-15 min, 5% solvent B. The separated analytes were detected using a 5500 Qtrap linear ion trap quadrupole mass spectrometer equipped with a Turbo V ion source operated in the ESI mode (AB Sciex) with Analyst software (Version 1.5). The Source and Gas were set as follows: gas 1, nitrogen (45 psi); gas 2, nitrogen (40 psi); ion spray voltage, 5,500 V; ion source temperature, 400 °C; and curtain gas, nitrogen (30 psi). The mass spectrometer was operated in the multi-reaction monitoring mode (MRM).
Bisulfite genomic sequencing. The PCR primers were designed using the online MethPrimer software (www.urogene.org/methprimer) and were listed in Table S6. Approximately 500 ng of the genomic DNA were subjected to the sodium bisulfite treatment and purification using an EZ DNA Methylation-Gold Kit (ZYMO), and a one-fifth volume of the eluted products was used as the template. PCR amplification was performed with a thermal cycling program of 95 °C for 1 min; 38 cycles of 95 °C for 30 s, 58 °C for 30 s, and 72 °C for 30 s; and then a final 5 min extension at 72 °C. The products were purified using the QIAquick Gel Extraction Kit (Qiagen) and sub-cloned. Twenty colonies for each PCR product were selected for sequencing using the 3730 Genetic Analyzer (Applied Biosystems) to evaluate the methylated cytosine levels.
RNA isolation and quantitative RT-PCR. RNA isolation and reverse-transcription were performed as previously described 59 . The primers used for real-time PCR were shown in Table S7. Real-time PCR was performed with the SYBR ® Premix Ex Taq TM (TaKaRa) on an ABI 7500 Real-Time PCR System using the SYBR Green detection protocol. Briefly, the amplification mixture was 0.5 μ M primers, 12.5 μ L of SYBR ® Premix Ex Taq TM , and 2 μ L of template cDNAs in a total volume of 25 μ L. The samples were amplified with the following program: initial denaturation at 95 °C for 30 sec, followed by 45 cycles of denaturation for 5 s at 95 °C, and annealing/elongation for 20 s at 60 °C. All qPCR reactions were run in triplicate. The qPCR results were exported as cycle thresholds (Ct). The relative level of expression of each gene was normalized to the expression level of the GAPDH transcript.
Computational processing of the LHC-BS data. After removing the adapter sequences and filtering out the low quality reads, the bisulfite sequencing reads were directly aligned to the human reference genome (UCSC hg19) using BSMAP 2.73 with the default parameters 60 , which combines genome hashing and bitwise masking to achieve fast and accurate bisulfite mapping. The ratio of the DNA methylation measurement/level of a specific cytosine was then calculated as the number of reads supporting methylation divided by the total number of reads covering that cytosine. CGIs were defined as regions greater than 200 bp with a GC fraction greater than 0.5 and an observed-to-expected ratio for CpG greater than 0.65, as annotated in the UCSC genome browser. Further explanations of the assessments are detailed in the corresponding illustrations of the figures or tables.
DMRs were identified using a sliding window strategy; commonly covered CpG sites with sequencing depth ≥ 5X between MC-twin samples were selected for the analysis. Fisher's exact test was first performed to evaluate the significance of the differences in the methylation levels between two samples. CpG sites with significantly different methylation levels (P < 0.05) were identified as the initial DMR loci, and their subsequent sites were gradually merged into the candidate DMRs according to following criteria: first, the distance between two neighbouring CpG sites ≤ 300 bp; second, all CpG sites in one candidate DMR maintain the same methylation tendency (hyper or hypo); and third, one candidate DMR contains no less than 5 CpG sites. Then, for each of the above candidate DMRs, Fisher's exact test was performed again to filter out those regions with a P-value > 0.05 and a difference in mean methylation levels between the two sample < 20%. Finally, the Benjamini-Hochberg procedure (False Discovery Rate, FDR) was used to adjust the P-values for multiple comparisons with 0.05 thresholds. Data Access. The sequencing and processed data have been deposited in the Gene Expression Omnibus (GEO) with the accession number GSE67377. All of the raw UPLC-MS/MS data were deposited in the figshare database (http://dx.doi.org/10.6084/m9.figshare.1368311).