Individual variation and longitudinal pattern of genome-wide DNA methylation from birth to the first two years of life

Abstract

Prenatal development and early childhood are critical periods for establishing the tissue-specific epigenome, and may have a profound impact on health and disease in later life. However, epigenomic profiles at birth and in early childhood remain largely unexplored. The focus of this report is to examine the individual variation and longitudinal pattern of genome-wide DNA methylation levels from birth through the first two years of life in 105 Black children (59 males and 46 females) enrolled at the Boston Medical Center. We performed epigenomic mapping of cord blood at birth and venous blood samples from the same set of children within the first two years of life using Illumina Infinium Humanmethylation27 BeadChip. We observed a wide range of inter-individual variations in genome-wide methylation at each time point including lower levels at CpG islands, TSS200, 5′UTR and 1st Exon locations, but significantly higher levels in CpG shores, shelves, TSS1500, gene body and 3′UTR. We identified CpG sites with significant intra-individual longitudinal changes in the first two years of life throughout the genome. Specifically, we identified 159 CpG sites in males and 149 CpG sites in females with significant longitudinal changes defined by both statistical significance and magnitude of changes. These significant CpG sites appeared to be located within genes with important biological functions including immunity and inflammation. Further studies are needed to replicate our findings, including analysis by specific cell types, and link those individual variations and longitudinal changes with specific health outcomes in early childhood and later life.

Keywords: :

• CpG
• DIP test
• DNA Methylation
• empirical bayes
• genome-wide
• normal mixture

Introduction

DNA methylation on the carbon-5 position of the pyrimidine ring of cytosine in a CpG dinucleotide is one of the epigenetic marks that influences gene expression.¹ In contrast to relatively stable genetic sequences, DNA methylation patterns are dynamic in response to internal and external environmental stimuli throughout the lifecycle.^2,3 For example, demographic factors (e.g., age and gender), environmental exposures (e.g., cigarette smoking and folate intake) and genetic variations (e.g., SNPs and CNVs) all can influence DNA methylation patterns.^4-8 Epigenetic variations are associated with normal growth and development and with health outcomes^2,3 and diseases including obesity, diabetes, cardiovascular disease, cancer and allergic disease.^1,9-14 There is a growing recognition that epigenetic mechanisms are likely the interface that connects the influence of genetic and environmental factors on gene expression, cell/tissue function and the development of disease.³

This report attempts to fill several gaps in the field. Prenatal development and early childhood are critical periods for the establishment of the epigenome.^3,15-18 Epigenomic profiles at birth reflect a net result of underlying genetic variations, in utero environmental exposures and stochastic epigenetic changes.¹⁶ Animal studies have demonstrated that early maternal care could influence the epigenetic processes of the offspring and subsequently contribute to inter-individual variations in future behavior.^15,17,18 Despite the critical role of the epigenome at birth and during early childhood in the health outcomes of later life, there is a particular lack of epigenomic data during this time period in humans. Several studies have examined DNA methylation in cord blood,^19-26 but these studies were either limited to global methylation (e.g., Alu and LINE-1 repeated sequences throughout the genome) or, in the case of genome-wide DNA methylation at a single nucleotide resolution, limited to a small sample size (i.e., 12 cord blood samples²⁰) or limited to the study of DNA methylation determinants^22,23,25 rather than looking into DNA methylation profiles at birth and longitudinal patterns from birth onward. While several studies have evaluated the longitudinal changes of genome-wide DNA methylation in peripheral blood cells in adults,^27-29 longitudinal epigenomic profiles from birth to early childhood remain unexplored. Previous studies have suggested that DNA methylation levels vary by specific locations across CpG and gene structures in the genome in both the normal condition and in the disease state. For example, when considering CpG islands relative to the RefSeq genes, the percentage of methylated CpG islands was highest in gene body, followed by 3′ end of a gene, intergenic, and 5′ end of a gene in human B cells.³⁰ In general, CpG islands have significantly lower methylation levels than non-CpG island regions across different tissues;^20,31-33 and most methylation alterations, for instance in colon cancer, occur at CpG island shores (i.e., < 2,000 bp from the borders of a CpG island) rather than at CpG islands.³⁴ However, no cross-sectional or longitudinal data are available to address whether DNA methylation levels vary by specific locations across CpG and gene structures on a genome-wide scale at birth and within the first two years of life.

Using a prospective birth cohort design, we sought to address the following questions: (1) what is the genome-wide DNA methylation profile at birth and in the first two years of life and does it vary by gender and among individuals? (2) What is the longitudinal pattern of DNA methylation from birth to the first two years of life and does it vary by gender and among individuals? And (3) do cross-sectional and longitudinal patterns of DNA methylation vary across CpG and gene structures throughout the genome?

A total of 105 Black children (59 males and 46 females) enrolled at Boston Medical Center were included in this study. We performed epigenomic mapping for paired cord blood samples at birth and venous blood samples from the same set of children within the first two years of life using Illumina Infinium Humanmethylation27 BeadChip. To our knowledge, this is the first and the largest epigenomic study of its kind in a prospective US birth cohort.

Results

Characteristics of study population

Our study included 105 pairs of cord blood samples at birth and venous blood samples within the first two years of life. Blood samples were obtained from 46 (43.81%) females and 59 (56.19%) males with a mean follow-up age of 0.99 (± 0.33) years. The study children were all Black, singleton, full-term births, with an average African ancestry proportion of 81% (± 15%). The average birth weight was 3301 (± 490) grams and average gestational age was 40 (± 1.2) weeks.

Distributions and variations of methylation levels

Figure 1 shows the distributions of methylation levels of 27K CpG sites within and across each study sample. A bimodal pattern was observed for the methylation-level distribution of all CpG sites within each sample. Approximately 45% of the CpG sites had methylation levels less than 0.1 in both male and female samples at each time point. In addition, less than 20% of the CpG sites had methylation levels greater than 0.8. The inter-individual variability of CpG sites’ proportions with methylation levels < 0.2 appeared higher than those with methylation levels > 0.2. We also found a wide range of inter-individual variations in genome-wide methylation levels (Figs. S1andS2), although the majority of CpG sites (~90%) had a median absolute deviation (MAD) of less than 0.04 for both cord and postnatal samples.

Figure 1. Gender specific distributions [(A) male, (B) female] of methylation levels of 27K CpG sites at two time points. CB stands for cord blood and PB stands for postnatal blood. The y-axis is the proportions of CpG sites corresponding to methylation categories specified by the x-axis. A bimodal distribution pattern was observed for both cord and postnatal blood samples. Comparisons of proportions of CpG sites in different methylation categories were conducted between males and females, and between cord and postnatal blood samples using Wilcoxon Rank Sum test or Wilcoxon Signed-Rank test. The FDR adjusted p-values are presented in the table below the figure.

Six major types of methylation-level distributions of the 27K probes were identified across the 105 samples (Fig. 2). The majority of CpG sites had positively skewed (~60%) or negatively skewed distributions (~30%). About 2% of the CpG sites had bimodal distributions, and less than 1% was normally distributed. The number of CpG sites across six types of distributions was significantly different between males and females both for cord and postnatal samples. For example, the percentage of CpG sites with a negatively skewed distribution was 41% in males compared with 35% in females in cord blood samples. We further explored how bimodal CpG sites were distributed across chromosomes. In all 105 samples, we found that 538 (92%) bimodal CpG sites in cord blood samples and 521 (93%) sites in postnatal samples were on the X chromosome. However, in male cord blood samples, 52 (53%) bimodal sites were on the X chromosome, but in female cord blood samples only 2 (3%) sites were on the X chromosome. Similar results were observed in postnatal samples. It appears that CpG sites on the X chromosome account for the majority of bimodal distributions in combined samples and in male samples but not in female samples. We further investigated whether CpG sites with bimodal distribution across samples were associated with SNPs or frequent structural variations. Interestingly, the majority of these CpG sites on the X chromosome were away from the known SNPs.

Figure 2. Six major types of DNA methylation distributions across 105 samples. The number of CpG sites (and the proportion) with each specific distribution is presented in the table below the figure. Bimodal (A) and unimodal distributions (B to F) were identified using the DIP test (the threshold value for combined samples is 0.051, for male samples is 0.074, and for female samples is 0.067). Kurtosis (K) and Skewness (S) parameters were used to further classify unimodal distributions into five categories: (B) S < -0.03, (C) S > 0.03, (D) -0.03 ≤ S ≤ 0.03 and K > 0.05, (E) -0.03 ≤ S ≤ 0.03 and -0.05 ≤ K ≤ 0.05, and (F): -0.03 ≤ S ≤ 0.03 and K < -0.05. The kurtosis was estimated using the moment method and then subtract 3.

Longitudinal changes in methylation from birth to the first two years of life

DNA methylation changes over time were evaluated based on continuous β values and/or changes in methylation categorical status from birth to the first two years of life. Tables S1 and S2 listed CpG sites with significant genome-wide methylation level changes in the first two years of life for males and females, separately. Some gender differences were noticed at each time point, especially on the X chromosome (Fig. S3 and S4). However, changes over time were not significantly different between males and females in all chromosomes except for the X chromosome (Figs. S5 and S6). The inter-individual variability of methylation changes ranged from 0.003 to 0.1, which was similar within each gender. Categorical classification of CpG sites showed that about 95% of CpG sites stayed in the same category from birth through the first two years of life, and only 5% of CpG sites had methylation categorical changes ranging from hypo to intermediate, intermediate to hyper, intermediate to hypo, or hyper to intermediate (Table 1). When we combined the results of continuous methylation level changes and categorical status changes according to stringent criteria as described in the Statistical Analysis Methods section, we identified a total of 149 significant CpG sites in female samples and 159 significant CpG sites in male samples, and of those 98 CpG sites were shared by both males and females (Table S3). Among the 98 sites, 2% changed from intermediate to hypo-methylation, 27% changed from intermediate to hyper-methylation and 71% changed from hyper-methylation to the intermediate methylation category.

Table 1. Dynamic changes of 27K probes from birth to the first 2 years of life

Postnatal Cord		Hypo-methylation (% ± SD)	Intermediate methylation (% ± SD)	Hyper-methylation (% ± SD)
Overall	Hypo-methylation	55 ± 1.2	1 ± 0.6	0 ± 0.04
	Intermediate methylation	2 ± 0.5	16 ± 1.8	1 ± 1.1
	Hyper-methylation	0 ± 0.08	1 ± 0.7	24 ± 0.9
Male	Hypo-methylation	^a57 ± 1.1	1.5 ± 0.7	0 ± 0.04
	Intermediate methylation	1.8 ± 0.6	^b15 ± 1.3	1 ± 1.1
	Hyper-methylation	0 ± 0.03	1 ± 0.8	^c24 ± 0.9
Female	Hypo-methylation	^a54 ± 0.5	1 ± 0.4	0 ± 0.05
	Intermediate methylation	2 ± 0.4	^b17 ± 1.0	1 ± 1.1
	Hyper-methylation	0 ± 0.1	1 ± 0.6	^c23 ± 0.7

The comparisons were conducted between males and females for proportions of CpG sites classified in each cell. A three-component normal mixture model was used to classify CpG sites within each sample into three categories (Fig. S14). ^a,b,cWilcoxon-Rank Sum test p < 0.0001. All other cells between males and females were not significantly different.

Using the significant CpG sites identified above, we further conducted network analyses using IPA software and presented five networks showing the top scores for males and females, separately (Tables S1 and S2). The top functions of networks were largely overlapped across gender, and were mainly associated with hematological system development and function, or immune and inflammatory response. We present the top network with the highest score in Figures S7 (males) and S8 (females) and their detailed biological functions in Tables S7 and S8. Specifically, for males the methylation levels of the tested probes in six genes increased while ten decreased during the first two years of life. For females, the number of CpG site annotated genes with increased and decreased methylation levels during the first two years of life was seven and nine, respectively. We also conducted IPA analysis for 98 common CpG sites across gender (Table S6) and detailed results are presented in Table S3.

Methylation patterns within CpG island groups and gene structures

We investigated methylation patterns within CpG island groups (Figs. 3 and S9) and different gene structures (Figs. 4 and S10) by gender at each time point. Among 27K probes, 43% were within CpG islands, 17% distributed in the north shore,³⁴ 15% were in the south shore, about 1% was within the north shelf (a gene feature beyond the north shore³⁵) and the south shelf (a gene feature beyond the south shore³⁵) each, while the rest were located at non-CpG islands and their nearby features. We found a “U” shape pattern of methylation levels across CpG island groups, i.e., CpG islands had the lowest methylation levels, and the further away from the CpG islands the higher the methylation levels became (Fig. 3). This “U” shape pattern was similar for both males and females at each time point. In parallel with the “U” shape cross-sectional pattern, longitudinal changes were similar between north and south shelves and between north and south shores, though they all differed from that found in CpG islands (all p < 0.0001). This observation was true for both males and females (Fig. 3). We did find a significant gender difference in the distributions of methylation levels within CpG islands for both cord (p < 0.0001) and postnatal samples (p = 0.0084), but this was not the case in other CpG island groups.

Figure 3. Distributions of mean methylation levels are presented at different CpG island groups for cord blood samples at birth (A) and postnatal venous samples within 2 y of life (B). For each boxplot, the dot in the box represents the mean methylation level, and the horizontal bar in the box is the median methylation level. The two horizontal bars below the 25th percentile bars (the bottom side of the box) reflect the minimum and 5th percentile, respectively. Similarly, the top two bars above the 75th percentile (the top side of the box) reflect the 95th percentile and the maximum, respectively. The lower two panels are density plots of changes in methylation levels (postnatal minus cord) for males (C) and females (D), respectively.

Figure 4. Distributions of mean methylation levels are presented at different locations of the gene structures at two time points. For each boxplot, the dot in the box represents the mean methylation level, and the horizontal bar in the box is the median methylation level. Two horizontal bars below the 25th percentile bar (the bottom side of the box) reflect the minimum and 5th percentile, respectively. Similarly, the top two bars above the 75th percentile bar (the top side of the box) reflect the 95th percentile and the maximum, respectively. The lower two panels are probability density plots of the mean methylation level changes (postnatal minus cord) for males and females, respectively.

Mean methylation levels were comparable among TSS200 [left 200 bp from the transcription start sites (TSS)], 5′UTR, and 1stExon, but were higher in the other three locations including gene body, TSS1500 (left 1500 bp from the transcription start sites), and 3′UTR (Fig. 4). The probability densities of mean methylation levels were similar among TSS200, 5′UTR, and 1st Exon, and between TSS1500 and gene body, but they were quite different compared with the levels at 3′UTR (Fig. S10). Compared with TSS200, the distributions of mean methylation levels were significantly different from other gene structures for both cord and postnatal samples (all p < 0.0001) in males and females. Mean longitudinal changes were similar at TSS200, 5′UTR, and 1st Exon, but they shaped differently at 3′UTR (Fig. 4). Compared with that at TSS200, the distributions of mean methylation longitudinal changes were significantly different at all other gene structures (p < 0.0001) except 5′UTR for both genders.

Inter-individual variability of methylation levels was the smallest in the CpG islands and was significantly higher in other island groups (all p < 0.0001). This pattern was observed at each time point and in both genders (Fig. S11A, B and C). Inter-individual variability at TSS200 was similar to that at 5′UTR and in the 1stExon regions (p > 0.05), but was significantly lower than that at TSS1500, gene body and 3′UTR (all p < 0.0001). This observation was identified in both genders at each time point, and also was true across two time points (Fig. S11D, E and F). There was no gender difference for inter-individual variations at different CpG island groups or at different gene locations at each time point, except for at birth when males had a higher inter-individual variability at CpG islands than did females (p < 0.0001).

Among 98 common significant CpG sites identified through longitudinal analysis, 53% were not annotated on the known CpG islands, 9% were in the CpG islands, 2% were in the north shelf, 5% were in the south shelf, 13% were in the north shore and 18% were in the south shore. In terms of their locations within gene structures, 35% were in TSS1500, 13% in TSS200, 21% in gene body, 14% in 5′UTR, and 17% in 1stExon (Table S3). It appears that methylation levels at CpG shores and/or distant gene structures tend to experience more active methylation changes over time.

Discussion

This is the first and the largest epigenetic study to investigate genome-wide variation of DNA methylation levels and longitudinal changes from birth to the first two years of life in a prospective US birth cohort. This study contributes new information as discussed below.

Distribution and individual variability of DNA methylation levels

At single-nucleotide resolution on a genome-wide scale at birth and in the first two years of life, we demonstrated bimodal distributions for methylation levels of all CpG sites within each study sample. Our findings are consistent with previous reports.^36,37 We also observed that a similar bimodal pattern is sustained across the structural features of the genes and the CpG islands, but undergoes a shifted distribution with regards to methylation levels: low methylation levels (e.g., < 0.2) is prominent in CpG islands and the 5′ region of the gene; while high methylation levels (e.g., > 0.8) is pronounced in distant regions from CpG islands and the gene body (Figs. S9 and S10). These patterns were similar across gender and across two time points. Our findings corroborated previous findings showing that DNA methylation distributions were different between non-CpG island regions and CpG islands, as was found in a UK epigenomic study of 12 fetal cord blood samples.²⁰ More importantly, we provided a comprehensive depiction of distributions for the detailed classification of CpG islands and gene structures in a large prospective US birth cohort.

We also observed that different CpG sites had different methylation distributions across the sample population, including skewed distributions in the majority and a small percentage of bimodal distributions. This fact suggests the importance of choosing appropriate analytical methods that consider the diversity of methylation distributions, for example, non-parametric tests (e.g., Wilcoxon Rank Sum test) represent better choices than parametric tests (e.g., t-test) for methylation data analysis.

In addition, we found that overall inter-individual methylation variations were small (Figs. S1 and S2). However, the variations of methylation levels had a large range and were different with regards to CpG islands groups (Figs. 3 and S9) and gene structures (Figs. 4 and S10): they were smallest in CpG islands, followed by shores and shelves; and they were smaller in nearby TSS regions (i.e., TSS200, 5′UTR, 1stExon) than in distant TSS regions (i.e., TSS1500, gene body, and 3′UTR). The largest methylation variation observed was at 3′UTR, however this finding should be interpreted with caution given that there were only a small number of probes covering this region (e.g., only 0.5% of the probes were located at 3′UTR). Still, we found that neither gender nor different time points had an influence on these variation patterns (Fig. S11). These findings suggest that the magnitude of inter-individual methylation variations observed for a given locus is associated with its biological function, i.e., the CpG islands under strongest functional constraint have the smallest inter-individual methylation variations. Similar conclusions apply to gene regions close to TSS vs. those further away from TSS. Consistent with a recent study on DNA methylation entropy,³⁸ our results indicated that methylation variation could be loci-specific and associated with the functional aspects of the corresponding genomic sequences. Higher inter-individual variations in the CpG shores and shelves and in TSS1500 indicated that these features were possibly more likely to be influenced by environmental factor exposure in utero and during early life, and thus could be strong candidates for disease differential methylation regions (DMR).²⁹

Longitudinal changes of methylation levels

We demonstrated that DNA methylation levels in general are also stable within the first two years of life, which is similar to a previous report in adults in which there was 8 to 11 y separating the sample collection.^3,28 Our findings support the belief that epigenomes originate mainly in utero and then are maintained through mitotic inheritance. On the other hand, we did identify approximately 0.5% CpG sites with substantial meaningful changes in DNA methylation levels in the first two years of life using stringent criteria as described in the Statistical Analysis Methods section. Interestingly, these longitudinal changes were mainly located in the CpG and gene features with higher variability (e.g., TSS1500 or CpG shores). We speculate that early environmental factors such as breastfeeding, infection and immunization could contribute to dynamic epigenetic changes, which in turn affect individual risk for the development of disease. Data from our IPA indicate that many of the significant CpG sites are associated with genes that are known to be involved in the development and function of various systems including the immune system, which begins in utero and matures during the first few years of life.^39,40

Our study has the following merits. We have presented the first genome-wide DNA methylation profiles at birth and within the first two years of life in a prospective US birth cohort. This study is also the first attempt to evaluate the longitudinal changes of DNA methylation levels within the first two years of life. To minimize technique variability, each pair of samples was placed on the same chip. The high correlations and high reproducibility (Fig. S12) of the duplicated samples, which were placed on different chips, indicate the high quality of the methylation data. Additionally, in the analytical phase we removed potential batch effects in the experiments using an empirical Bayes method implemented in program ComBat,⁴¹ which in general outperforms the other five batch adjustment methods commonly used in gene expression studies.⁴² We performed in-depth analysis by considering gender, two time points, CpG islands and gene structures. Finally, we ranked top probes/genes based on the combined criteria by not only considering the changes in methylation levels and methylation categories, but also the number of subjects within the category change. This method is especially useful to reduce false positive findings due to small variations in methylation data.

Limitations and Future Recommendations

When interpreting our findings, the following limitations need to be considered. We used Illumina HumanMethylation27 BeadChip for this longitudinal study. It is well known that the BeadChip has a limited coverage of the full range of genomic CpG sites, however it has high quantitative accuracy because of the high correlation coefficients between the Infinium assay and the GoldenGate and MethyLight assays (r² > 0.8)⁴³ (www.illumina.com/support/literature.ilmn). A recent study used the HumanMethylation27 assay as a common reference to evaluate the accuracy of sequencing-based methods for methylation.⁴⁴ Consequently, the quantitative accuracy rather than the coverage of this assay may be more critical to the present study, which focuses on the variability of DNA methylation levels and longitudinal changes. Another limitation is that this study assessed DNA methylation levels in whole blood samples, which reflect a net result from a mixture of various types of blood cells. As such, the observed variations can be affected by the proportions of various cell types. The degree to which our observed patterns and variability can reflect DNA methylation in each specific type of blood cell remains to be determined. Nevertheless, DNA samples extracted from whole blood have been commonly used in methylation studies, including two longitudinal studies in adults conducted by Feinberg et al.^27,28 Previous studies also have demonstrated that the cellular heterogeneity in leukocytes did not contribute to the variation of DNA methylation for the majority of the candidate loci.⁴⁵ Whether this is the case for a genome-wide scan needs to be investigated in the future. From a clinical and public health perspective, blood samples are the most easily accessible biospecimens, and the variability of genome-wide DNA methylation in peripheral blood, if linked to health and disease, should be of great value for disease diagnosis, prevention, and therapy. In addition, follow-up time of these 105 children (59 males and 46 females) varied from 6 mo to 2 y of age. The median follow-up is 10 mo for males and 11 mo for females. Due to sample size limitation, we were unable to further stratify by the follow-up length in the analyses since the estimates will be unstable. As such, we could not eliminate any possible influence of varied length of follow-up age on DNA methylation patterns. Finally, since the Illumina BeadChip does not provide a direct measure of DNA methylation, our findings need to be replicated in independent samples and validated by pyrosequencing. Validation using orthogonal techniques will provide meaningful biological interpretations if methylation findings and phenotypes are linked together.

Despite the limitations, our study represents the first step in understanding the patterns, variations, and longitudinal changes of the epigenome in early life, which may be relevant to health and disease in later life. The Boston Medical Center cohort is an ongoing birth cohort, and the continued follow-up of this cohort will enable us to distinguish the consistently stable CpG sites from the dynamic ones. With this knowledge, stable CpG sites could contribute to personalized epigenotypes and become the focus of future studies to identify the earliest biomarkers for complex diseases, while dynamic CpG sites could be considered environmentally inducible epigenetic markers and become the focus of future studies to understand the mechanisms underlying the associations between postnatal exposure and the development of disease.

Materials and Methods

Study subjects

The study samples are part of an ongoing birth cohort study in which mother-infant pairs are enrolled at Boston Medical Center at birth and infants are followed from birth onward. The initial recruitment and follow-up have been previously described.^46-48 This report included 105 Black children (59 males and 46 females), who were singleton, full-term births and whose cord blood and venous blood were collected at birth and within the first 2 years of life, respectively. Written informed consents were obtained from all participants in this study. The study protocols have been approved by the Institutional Review Boards (IRBs) of Children’s Memorial Hospital in Chicago and the Boston University Medical Center.

Epigenetic mapping and quality control

Genomic DNA was extracted from whole blood leukocytes and quantified with SpectraMax M2 (www.moleculardevices.com). Genome-wide DNA methylation was assessed using Illumina HumanMethylation27 BeadChip at the Genomics Core Facility, Northwestern University. Each pair of study samples was placed on the same chip. We also included six duplicates, which were randomly selected from the 105 paired DNA samples. A β value signifying the percentage of methylation ranging from 0 to 1.0 was reported for each of the 27,578 probes from Illumina’s GenomeStudio. Sample-independent (e.g., hybridization and extension) and sample-dependent (e.g., bisulfite conversion and mismatch) control analyses illustrated the high quality of the microarray process and high performance of the samples, respectively. Pearson correlation coefficients for the six duplicates were > 0.99, and 99% of CpG sites had an absolute methylation-level difference of less than 0.05 between two duplicates (Fig. S12). All of the samples had greater than 98% CpG call rates (i.e., > 27,000 probes).

Statistical analyses method

The 105 paired samples were hybridized on 18 chips. To eliminate potential batch effect, an empirical Bayes methodology implemented in ComBat software⁴¹ was applied to all data from each of the different chips to effectively remove batch effects in the experiment.^42,49 Analysis of variance (ANOVA) was conducted to evaluate the batch effects for normalized β values and showed that about 96% of CpG sites had significant batch effects before normalization, but this proportion dropped to 0.12% after normalization using ComBat. All of the analyses were based on methylation levels after performing ComBat and were stratified by gender due to potential gender differences in DNA methylation.^6,50

We first assessed the overall distributions of methylation levels across all 105 samples and within each sample. The DIP test⁵¹ was used to identify bimodal distributions across samples. We used a graphic interpolation method to identify appropriate significant cut-off values for DIP statistics specific to a sample size (Fig. S13). The cut-off value for the DIP statistics at a sample size of n = 105 (combined samples) was 0.051, and was 0.074 for n = 59 (males) and 0.076 for n = 46 (females). Kurtosis (K) and skewness (S) parameters were used to classify unimodal methylation distributions in the sample population. Specifically, a normal distribution had -0.03 ≤ S ≤ 0.03 and -0.05 ≤ K ≤ 0.05 (K was estimated using the moment method and subtract 3); a negative skewed distribution had S < -0.03; a positive skewed distribution had S > 0.03; a platykurtic distribution had -0.03 ≤ S ≤ 0.03 and K < -0.05; and a leptokurtic distribution had -0.03 ≤ S ≤ 0.03 and K > 0.05. We used 0.03 and 0.05 as cutoff values of skewness and kurtosis instead of their theoretic values to consider the small sample variability of estimation.

Second, we investigated distributions and variations of methylation levels for different CpG island groups (i.e., CpG islands, CpG shores and CpG shelves) and gene structures (i.e., transcription start sites 200 (TSS200), TSS1500, 5′UTR, gene body, and 3′UTR) defined by Illumina at each time point for males and females, separately. The Kolmogorov-Smirnov test⁵² was utilized to compare empirical distributions of methylation levels between different CpG island groups and gene structures. Median absolute deviation (MAD) was calculated to measure inter-individual variations of methylation levels, and the Wilcoxon Rank Sum test was used to compare MAD shift between different groups (e.g., between different CpG island groups) and to compare differences in proportions of CpG sites and methylation level changes between male and female samples at each time point.

Finally, we investigated longitudinal changes in methylation levels within the first two years of life and conducted pathway analysis using Ingenuity Pathway Analysis (IPA), as described below. We used the Wilcoxon Signed-Rank test to compare methylation changes for all CpG sites over two time points, and to compare differences of proportions of CpG sites for each methylation level category over time for each gender. A three-component normal mixture model^53,54 was applied to classify CpG sites into different methylation categories (i.e., hypo, intermediate or hyper) (Fig. S14). This approach is less arbitrary compared with that used in previous studies, which defined hypo-methylation or hyper-methylation using a cutoff of 20% or 80%.^36,55 The Stuart-Maxewell^56,57 test was used to test methylation categorical status changes over time. We selected CpG sites with significant longitudinal changes by merging results derived from the Wilcoxon Signed-Rank test and Stuart-Maxewell test. The combined results were sorted ascendingly by two FDR adjusted p-values, and then were sorted descendingly by both absolute changes in methylation levels and total number of discordant pairs (i.e., the number of subjects with discordant methylation categories). For each gender group, we included CpG sites for the network analyses by IPA only when the following criteria were met: methylation level changes greater than 0.1 in β values and multiple comparison corrected p values less than 0.05; and methylation categorical status changes in at least half of the samples. The Benjamini and Hochberg false discovery rate⁵⁸ was used to adjust p values for multiple testing. A sample size of 46 has 91% power to detect an effect size of 0.8 for an individual CpG site assuming the standard deviation is 0.1 (i.e., the difference is 0.08). We also assume that 150 CpG sites among 27K probes are truly significantly different over time, and the false discovery rate is controlled to be 5% for the power analysis. Thus, we had appropriate power to detect methylation difference of 0.1, which was used as one criterion for selecting CpG sites for IPA analysis.

Network analyses

The CpG sites to be included in the IPA were annotated to the corresponding UCSC genes. For these selected genes, we conducted network analyses using IPA (www.ingenuity.com/index.html). Networks were generated based on ‘seeds’ (i.e., network eligible molecules), which are encoded by the selected genes and also interact with other molecules in the Ingenuity Pathways Knowledge Base. A high score for a network indicates a more approximate fit between network eligible molecules and the molecules that constitute the network. For the top-ranked network, IPA also identified their top biological functions.

Abbreviations:
DIP	=	a name of statistical test to identify bimodal distribution
TSS	=	the transcription start sites
MAD	=	median absolute deviation
CpG	=	genomic regions that contain a high frequency of CpG sites
IPA	=	ingenuity pathway analysis
DMR	=	differential methylation regions

Acknowledgments

The parent study is supported in part by the March of Dimes PERI grants (PI: Wang, 20-FY02–56), the NIEHS (PI: Wang, R21 ES011666), and the NICHD (PI: Wang, R01 HD041702). The study is supported in part by the Food Allergy Initiative, the March of Dimes Birth Defects Foundation (PI: Tsai, 21-FY07-605) the NIAID (PI: Wang, R21AI079872; R21AI088609; U01AI090727), and the NICHD (PI: Wang, R21HD066471). This study also is supported in part by the Department of Defense (PI: Wang, W81XWH-10–1-0123). Dr. Liu is supported by the NIAID (PI: Liu, R21AI087888). None of the authors have a conflict of interest pertaining to this work.

Disclosure of Potential Conflicts of Interest

No potential conflicts of interest were disclosed.

Author Contributions

Designed the study and provided direction: XW XL DW. Conducted data analyses and interpretation of results: DW XL HX YZ XW. Performed the lab analysis: XH GW RL. Draft and revision of the manuscript: DW XL HX HJT XH GW RL XW. Provided financial support: XW XL HJT.