Pre-adolescence DNA methylation is associated with lung function trajectories from pre-adolescence to adulthood

The pattern of lung function development from pre-adolescence to adulthood plays a significant role in the pathogenesis of respiratory diseases. Inconsistent findings in genetic studies on lung function trajectories, the importance of DNA methylation (DNA-M), and the critical role of adolescence in lung function development motivated the present study of pre-adolescent DNA-M with lung function trajectories. This study investigated epigenome-wide associations of DNA-M at cytosine-phosphate-guanine dinucleotide sites (CpGs) at childhood with lung function trajectories from childhood to young adulthood. DNA-M was measured in peripheral blood at age 10 years in the Isle of Wight (IOW) birth cohort. Spirometry was conducted at ages 10, 18, and 26 years. A training/testing-based method was used to screen CpGs. Multivariable logistic regressions were applied to assess the association of DNA-M with lung function trajectories from pre-adolescence to adulthood. To detect differentially methylated regions (DMRs) among CpGs, DMR enrichment analysis was conducted. Findings were further tested in the Avon Longitudinal Study of Parents and Children (ALSPAC) cohort. Pathway analyses were performed on the mapped genes of the identified CpGs and DMRs. Biological relevance of the identified CpGs was assessed with gene expression. All analyses were stratified by sex. High and low trajectories of FVC, FEV1, and FEV1/FVC in each sex were identified. At PBonferroni < 0.05, DNA-M at 96 distinct CpGs (41 in males) showed associations with FVC, FEV1, and FEV1/FVC trajectories in IOW cohort. These 95 CpGs (cg24000797 was disqualified) were further tested in ALSPAC; 44 CpGs (19 in males) of these 95 showed the same directions of association as in the IOW cohort; and three CpGs (two in males) were replicated. DNA-M at two and four CpGs showed significant associations with the corresponding gene expression in males and females, respectively. At PFDR < 0.05, 23 and 10 DMRs were identified in males and females, respectively. Pathways were identified; some of those were linked to lung function and chronic obstructive lung diseases. The identified CpGs at pre-adolescence have the potential to serve as candidate markers for lung function trajectory prediction and chronic lung diseases.


Introduction
The patterns of lung function development, from preadolescence to adulthood, play a major role in the pathogenesis of respiratory diseases. Recent studies have highlighted that reduced lung function development in Open Access *Correspondence: hzhang6@memphis.edu 1 Division of Epidemiology, Biostatistics, and Environmental Health Sciences, School of Public Health, University of Memphis, Memphis, TN 38152, USA Full list of author information is available at the end of the article young adulthood predisposes to respiratory and other chronic diseases in later life and is also associated with early mortality [1,2]. Lung function grows dramatically throughout childhood and reaches its peak in adolescence or early adulthood. After a brief period of stable lung function in early adulthood, a gradual decline ensues with aging [3][4][5]. Previous studies have demonstrated that early decline of lung function, and/or failure to reach maximal level lung function (even with a normal rate of decline), is associated with the development of chronic obstructive pulmonary disease (COPD) in later life [3][4][5][6], suggesting that the origins of COPD lie, in part, in early life [6,7]. COPD is projected to become the third leading cause of death worldwide by 2030 [8,9], highlighting how insights into the trajectories of lung function development from childhood-to-young adulthood would be beneficial for COPD prediction, prevention, and management.
Encouraged by the significance and advantage of longitudinal designs, we and others examined the temporal trend of lung function growth and decline through multiple important stages of life: childhood, adolescence, and adulthood [4,[10][11][12]. These studies demonstrated that there are distinct groups of individuals with a persistently low lung function trajectory from childhood-toadulthood, suggested a potential connection with COPD in later life. One study showed weak evidence that persistently low FEV 1 trajectory is associated with genetic factors in addition to multiple childhood exposures [10]. A recent study based on repeated measurement of lung function in adults reported that genetic variants associated with cross-sectional lung function measurements were not associated with a longitudinal decline of lung function [13]. These inconsistent findings in genetic studies and the clear impact of environmental factors on lung function motivated the investigation of the role of epigenetic factors such as DNA methylation (DNA-M) in determining variation in lung function between people and over time.
DNA-M represents an epigenetic mechanism that regulates gene expression, which consequently influences disease risk [14,15]. Growing evidence indicates that DNA-M in whole blood is associated with lung function and its related diseases such as asthma and COPD [15][16][17][18]. Pre-adolescence adverse exposure is shown to be associated with adulthood chronic lung diseases [19]. As an epigenetic memory of past exposures, the role of preadolescence DNA-M on lung function trajectories from pre-adolescence to young adulthood is unknown [19]. We hypothesized that differential methylation at certain cytosine-phosphate-guanine dinucleotide sites (CpGs) in childhood is associated with the trajectories of lung function. Given that lung function growth and decline is sex-dependent and such dependence is attributable to multiple biological determinants, including dimensional/ anatomical, immunological, and hormonal determinants [20][21][22][23], we examined the hypothesis in male and female participants, separately [12,24]. The study was carried out in the birth cohort located on the Isle of Wight (IOW) in the UK. To assess the potential of generalizability, an independent UK birth cohort, the Avon Longitudinal Study of Parents and Children (ALSPAC) cohort, was used for replication.

Results
In the complete IOW cohort (n = 1456), lung function measurements at ages 10, 18, and 26 years were available for 980 (67.3%), 838 (57.6%), and 547 (37.6%) participants, respectively. A total of 377 male and 432 female participants were included for trajectory analyses, and each of the participants had spirometry tests at two or more of the three ages (Fig. 1). The analysed sub-sample (n = 809) was not statistically different from the complete cohort (n = 1456) regarding FVC, FEV 1 , and FEV 1 /FVC at the corresponding ages (Table 1).

Lung function trajectories
In trajectory analyses, two distinct lung function trajectories from pre-adolescence to adulthood (10 to 26 years of age) were identified in the IOW cohort, labelled as 'low trajectory' and 'high trajectory' for FVC, FEV 1 , and FEV 1 /FVC in both male and female participants (Fig. 2). Among the 377 male participants, 199 (52.8%), 204 (54.1%), and 96 (25.5%) were assigned to low FVC, FEV 1 , and FEV 1 /FVC trajectories ( Fig. 2: gray "dashed lines"), respectively, using probability > 0.5 to define class membership. Among these male participants, at least 82% and 92% of them had a trajectory assignment probability ≥ 0.7 for low and high trajectories, respectively.

Pre-adolescence DNA-M and lung function trajectories
In total, 176 of the 377 male and 136 of the 432 female participants included in the analyses had DNA-M data available at age 10 years (Fig. 1). In screening, 119 distinct CpGs for males ( associations with lung function trajectories, with the effects of confounders adjustment. There was no overlap between the 119 and 56 CpGs identified in males and females, respectively. Using multivariable logistic regression models, DNA-M levels at age 10 years of 11, 13, and 17 CpGs in males, and 21, 21, and 16 CpGs in females were statistically significantly associated with FVC, FEV 1 , and FEV 1 / FVC trajectories from pre-adolescence to adulthood, respectively, after correcting for multiple testing using the Bonferroni approach. Among the 96 distinct CpGs identified in the IOW cohort, 95 were further examined in the ALSPAC (cg24000797 was excluded during quality control in ALSPAC).

Testing IOW cohort findings in the ALSPAC
In total, 4,861 participants (males = 2216) in ALSPAC had FVC, FEV 1 , and FEV 1 /FVC measurements for more than a single time point at ages 8-, 15-, and 24-year follow-up and were included in the trajectory analyses. Of these participants, 691 had DNA-M data at age 7 years.
We identified two trajectories, low and high, for FVC, FEV 1 , and FEV 1 /FVC across 8, 15, and 24 years, comparable to those from the IOW cohort (Fig. 2). Next, for the 95 CpGs identified in the IOW cohort, we tested the association of DNA-M at age 7 years with lung function trajectories. Among the 95 CpGs, DNA-M at 44 distinct CpGs (2 CpGs overlapped between FEV 1 and FVC in females) showed the consistent association with lung function trajectories with those in the IOW cohort in terms of the direction of association. These include 19 CpGs (6 CpGs for FVC, 7 for FEV 1 , and 6 for FEV 1 / FVC) in males, and 25 CpGs (11 CpGs for FVC, 9 for FEV 1 , and 7 for FEV 1 /FVC) in females ( Table 2, Fig. 4, Additional file 1: Table S1). These 44 CpGs were noted as IOW-ALSPAC consistent CpGs. Among these CpGs, cg14669749 and cg21131402 showed statistically significant associations with FEV 1 /FVC trajectories in males and cg23987789 with FEV 1 trajectories in females. DNA-M at three CpGs showed marginal statistical significance in females, two for FEV 1 /FVC trajectories (cg23190164, p = 0.08 and cg24479027, p = 0.09) and one for FVC (cg05597624, p = 0.088). At 74% of the 19 identified CpGs in males, a higher DNA-M was associated with a higher odds of being in the low lung function trajectories, while in females, the percentage was 44%.

Detection of differentially methylated regions
For DMR enrichment analysis, a frequency of 20 and above was focused in screening to secure a sufficient number of CpGs. In females, 486, 518, and 461 CpGs and in males, 419, 559, and 842 CpGs for FVC, FEV 1 , and FEV 1 /FVC trajectory, respectively, were selected after screening and included in DMR analyses for each trajectory. After controlling FDR of 0.05, 23 statistically significant DMRs in males and 10 in females were identified (Table 4 with breakdown for each lung function trajectory). In total, 78 CpGs were in the 33 (23 + 10) identified DMRs (Additional file 1: Table S2), of which two CpGs (cg09707262 and cg02304879) were also among the 44 Fig. 2 Distinct lung function trajectories from childhood-to-adulthood following comparable patterns in the IOW and ALSPAC cohorts. Note: Among the subjects assigned to each trajectory with a probability > 0.5, most assignments were with a probability ≥ 0.7, much higher than 0.5. In the following, we provide, for each sex and lung function parameter, the percentages of assignments with an assignment probability ≥ 0.7: (1) Among the male participants who were assigned to persistently low lung function trajectories with a probability > 0. 5 CpGs individually identified CpGs. The CpGs in the identified DMRs with the mapped genes and chromosomes, and the corresponding p values of DMRs for each sex are presented in Additional file 1: Figure S1.

Biological pathways of mapped genes
The 44 IOW-ALSPAC consistent CpGs (19 CpGs for males and 25 for females) and the 78 CpGs in the 33 identified DMRs (23 DMRs in males and 10 in females) were mapped to 42 ( Table 2) and 33 genes (Table 4), respectively. The distinct 73 genes were included in pathway enrichment analysis (focusing on pathways with at most 2000 genes) via bioinformatics tool ToppFun. After controlling the FDR of 0.05, six and 12 pathways were identified (Table 5) in males and females, respectively.

Discussion
Two major distinct lung function trajectories from pre-adolescence to adulthood in each sex were identified using latent class trajectory analyses in two  population-based birth cohort. We showed that preadolescence DNA-M at 44 CpGs was associated with the trajectories. These CpGs mapped to 42 genes, which were enriched in 18 KEGG and REACTOME pathways. We identified 23 and 10 DMRs associated with the lung function trajectories in males and females, respectively. We further evaluated the functional effects of the identified CpGs by integrating gene expression. DNA-M at two CpGs in males and four in females at age 10 years was longitudinally associated with gene expression at age 26 year among the distinct set of 15 CpGs in each sex. Since in our study lung function trajectories cover from pre-adolescence transition till adulthood, and DNA-M was measured at an age close to the transition rather than a number of years before adolescence, the identified 44 CpG sites have a strong potential of high sensitivity to predict an individual's lung function development from pre-adolescence to young adulthood. The high and low trajectories of lung function identified in this study were the same as in our previous study [12] except that in this study subjects who had lung function measurement only at one-time point were excluded to improve the accuracy of trajectory assignment. The identified trajectories were also consistent with the main two trajectories of the previous reports from populationbased studies including ours [4,10,12]. In a large-scale study, Belgrave et al. [10] (age 5-24 years) identified 4 distinct trajectories of FEV 1 ; persistently high; normal; below average; and persistently low. In our study, on the other hand, two (high and low) trajectories of FVC, FEV 1 , and FEV 1 /FVC represented the data the best, which was likely due to the relatively smaller sample size. Moreover, like Belgrave et al. study, in this study, participants in the low FVC and FEV 1 trajectory group did not achieve maximally attainable FVC and FEV 1 and in the low FEV 1 / FVC trajectory group showed an accelerated decline from age 10 to 26 years (15% and 11% decline in males and females, respectively), compared to the declines in the high trajectory (9% and 10%, respectively) ( Fig. 2), suggesting a risk of future COPD. The observations in this study also support the findings in the previous longitudinal studies by Lange et al. and Bui et al. [6,11] that the persistently low lung function trajectory is associated with the risk of COPD in adults.
To our knowledge, this is the first study examining the association of pre-adolescence DNA-M with lung function trajectories ranging from pre-adolescence to post-adolescence. The identified CpGs and DMRs at childhood may provide insight into the pathogenesis of variations in lung function growth in adolescence. In addition, the associations of methylation at some identified CpGs with gene expression, such as cg16709691 (LMF1), cg12655437 (SMAD2), cg16049690 (BTNL9), cg07562175 (FBRSL1), cg13168117 (KLHL30), and cg23987789 (VAMP3) manifest the functional importance of the CpGs as biomarkers. Among these genes, SMAD2, FBRSL1, and VAMP3 were associated with lung function, its related pathway, and COPD in previous studies [25][26][27][28]. Although most individually identified CpGs through logistic regressions were different from those in the DMRs due to different assumptions and statistical approaches between these two analyses, their mapped genes jointly involved at the biological pathways (Table 5).
Among the listed biological pathways linked to the mapped genes (Table 4), several pathways play a significant role in lung function and/or COPD, including downregulation of SMAD2/3: SMAD4 transcriptional log odds ratios (ORs) of IOW-ALSPAC consistent CpGs for the association of DNA-M at childhood with lung function trajectories from childhood-to-adulthood in males and females, separately The CpGs showed the same direction of associations and were significant at 0.05 or 0.1 level and were in bold font Chr. no chromosome number, ORs odds ratios, CIs confidence intervals a Genes located at intergenic location were not found in Illumina annotation file and were identified using online tool SNIPPER activity, circadian entrainment, GABA B receptor activation, and activation of G protein-gated potassium channels [25,[29][30][31][32]. For example, downregulation of SMAD2/3: SMAD4 transcriptional activity plays a role in the regulation of TGF-β1-induced collagen expression in lung. Excessive collagen deposition is one of the characteristics of idiopathic pulmonary fibrosis that lead to impaired lung function later in life. The association of cg12655437 with SMAD2 expression in this study revealed the pathway as functionally meaningful. Another pathway, the circadian rhythm regulates physiological diurnal variation of lung function through the autonomous peripheral circadian clock mechanisms. Clara cells in the bronchioles play a major role in such variations of lung function. These physiological oscillations are driven by transcriptional factors and genes such as PER3.
The CpGs showing consistent direction of associations with statistical significance at 0.05 or < 0.1 in both cohorts included cg14669749 (SKI), cg21131402 (C12orf50), cg23987789 (VAMP3), cg23190164 (LGR5), cg24479027 (ABR), and cg05597624 (RNF220). Some of these genes, such as VAMP3, LGR5, PER3, and SDC1, were found to be involved in the different physiological functions of lung and chronic lung disease [27,28,[33][34][35][36]. Among these genes, VAMP3 is found as one of the soluble N-ethylmaleimide-sensitive factor attachment protein receptors regulating mucin granule exocytosis. Mucin secretion is an innate immunity mechanism, which is harmfully upregulated in obstructive lung diseases including COPD [27,28]. In addition, being in the intergenic region, the significant positive association of methylation at cg23987789 with the expression of VAMP3 revealed a potential of this CpGs' functionally regulatory role.
LGR5 is related to the WNT signalling cascades, which are the critical regulators of different developmental and pathophysiological processes in lung. Dysregulated LGR5 expression influences to reduced WNT-β catenine signalling cascades, which is further linked to chronic lung disease including COPD [33,34]. Among the annotated genes of the DMRs and also identified in circadian entrainment pathway, in females, PER3 has been previously associated with childhood and adolescence lung function (FEV 1 ) [35]; in males, SDC1 was found as a differentially expressed genes in COPD development by robust rank aggregation method and in KEGG pathway in the previous study [36].
Although overall patterns of lung function trajectories in males and females were similar, for each identified trajectory, there existed large differences in volumes and flows between the two sex (Fig. 2). Such differences were expected and acknowledged in the literature [20][21][22][23] and were the major driving factors for the stratified study design. The uniqueness of identified CpGs for each sex led us to postulate the possibility of either different underlying epigenetic mechanisms in males and females in the regulation of gene activity and may act as the biomarkers of physiology and/or exposures that influence lung function trajectory. Another strength of this study was the time-lagged study design. Because of this, CpGs identified in this study have the potential to serve as candidate predictors for future lung function trajectories and will be beneficial to the detection of early lung diseases

Table 3 Association of DNA-M at IOW-ALSPAC consistent CpGs with gene expression
The CpGs which showed a significant association with gene expression at 0.05 level were in bold font a Genes located at intergenic location were not found in Illumina annotation file and were identified using online tool SNIPPER and subjects with a higher risk of developing those diseases.
Some issues related to the study designs and data analyses are worth discussing. In this study, participants with lung function measurements available only at one-time point were excluded from the analyses. The exclusion of subjects with missing data plus stratification by sex made the sample size smaller for each group, especially

Table 4 DMRs for lung function trajectory in relation to childhood DNA-M identified by DMRcate method
DMRcate annotates to UCSC RefGene from the Illumina annotation file a Genes were not found in Illumina annotation file and were identified using online tool SNIPPER b CpGs names are provided in Additional file 1: for females (n = 136). However, including participants with lung function measurement available at two or more time points ensured a high probability of trajectory group assignment rather random group assignment with a probability around 0.50. Also, based on the comparison with the whole study cohort, the study samples represented the whole cohort indicating that such restriction (≥ 2 repeated lung function measurements) did not bring statistically significant selection bias into the study samples. The use of a screening process, together with stringent control of multiple testing via the Bonferroni approach instead of controlling FDR, and the utilization of the replication cohort ensured that the findings from our study are robust with the potential of being generalized. Besides, in the IOW cohort, age 10 was treated as the pre-adolescence age, since almost all children (males 98% and females 92%) included in this study had not entered any phase of puberty. The children with minimal pubertal events were not excluded from the present study, since excluding them was not expected to alter the findings and conclusions but might have decreased testing power. Another perspective related to age is that, in the IOW cohort, the analyses were based on data collected at ages 10, 18, and 26 years representing pre-and post-adolescence. In the ALSPAC, the corresponding ages were 7/8, 15, and 24 years. The decline phase of lung function might have not started yet at age 24 years in the ALSPAC for some subjects (Fig. 3). This inconsistency between the two cohorts might be the cause of a lack of replication for some CpGs. Finally, the identified CpGs had minimal overlapping among FVC, FEV 1 , and FEV 1 / FVC trajectories, although in females, some mapped genes of identified CpGs associated with FVC, FEV 1 , and/or FEV 1 /FVC trajectories were involved in common pathways across these three different lung function parameters. Since DNA-M was measured in whole blood rather than in airways, although several studies support this non-invasive sampling approach of assessing DNA-M, the relevance of epigenetic changes measured in leukocytes in whole blood to gene expression in the lung remains unanswered and deserves further investigations of their biological evidence.

Study subjects and design The IOW cohort-Discovery cohort
The Isle of Wight (IOW) birth cohort is a populationbased birth cohort established in 1989, UK. The study was originally approved by the IOW Local Research Ethics Committee at recruitment, and further assessments of this cohort are approved by the National Research Ethics Service, Committee South Central-Southampton B (06/Q1701/34). Informed written consent was obtained from participants or their parents before participating.

Lung function
Pre-bronchodilator spirometric measurements, including forced vital capacity (FVC), forced expiratory volume in one second (FEV 1 ), and the ratio of FEV 1 over FVC (FEV 1 /FVC), were conducted at ages 10 (n = 980), 18 (n = 838), and 26 (n = 546) years and included in the study. FVC and FEV 1 were measured using a Koko spirometer and software with a portable desktop device (both PDS Instrumentation, Louisville, KY, USA), and the FEV 1 /FVC ratio was calculated. Spirometry was conducted and evaluated according to the American Thoracic Society (ATS) guidelines [38,39]. To conduct spirometry, participants were required to be free of respiratory infection for two weeks, not taking oral steroids, not taking any β-agonist for six hours and caffeine for at least 4 h.

Measurement of DNA methylation (DNA-M)
Peripheral blood samples collected at ages 10 years from n = 330 randomly selected subjects were used for DNA extraction via a standard salting out procedure [40]. DNA concentration was estimated by Qubit quantitation.
For each sample, one-microgram DNA was bisulphitetreated for cytosine to thymine conversion using the EZ 96-DNA methylation kit (Zymo Research, Irvine, CA, USA), following the manufacturer's protocol. DNA samples were randomly distributed on microarrays to control against batch effects. DNA-M was measured using HumanMethylation450K and HumanMethylationEPIC BeadChips (Illumina, Inc., SanDiego, CA, USA). Arrays were processed using a standard protocol as described elsewhere [41], with multiple identical control samples assigned to each bisulphite conversion batch to assess assay variability.

Preprocessing of DNA-M
Probes not reaching a detection p value of 10 -16 in at least 95% of samples were excluded. The same criterion was applied to exclude samples, i.e. a sample with detection p value of > 10 -16 in more than 5% of CpGs was excluded. CpGs on sex chromosomes were also excluded to avoid bias. DNA-M was then pre-processed using the "CPACOR" pipeline for data from both platforms [42]. DNA-M intensities were quantile normalized using the R computing package, minfi [43]. DNA-M β values for each CpG were then calculated as a ratio of methylated (M) over the sum of methylated and unmethylated (U) probes (β = M/[c + M + U]) interpreted as the percentage of methylation [44], where c was used as a constant to prevent zero in the denominator. Principal components (PCs) inferred based on control probes were used to represent latent variables explaining chip-to-chip and technical (batch) effects on DNA-M variations. Since DNA-M data were from two different platforms (450 k and EPIC), we determined the PCs based on DNA-M at shared control probes between the two platforms. In total, 195 control probes were overlapped between 220 control probes from the 450 K BeadChips and 204 from the EPIC array. These 195 shared probes were then used to calculate the control probe PCs, top 15 of which were used to represent latent batch factors [42]. After pre-processing, a total of 473,864 and 847,155 CpGs were available in the 450 K and EPIC methylation array data, respectively, and 439,635 overlapping CpGs were identified between the two platforms. CpGs with a single nucleotide polymorphisms (SNP) overlapping the detection probe with minor allele frequency ≥ 0.7% (corresponding to at least 10 subjects in the IOW cohort with n = 1456) within 10 base pairs of the targeted CpG were excluded due to potential bias that those SNPs brought to the measurement of DNA-M. After excluding probe SNPs, 402,714 CpGs were included in the statistical analyses.

Gene expression data
RNA-seq gene expression data for subjects at age 26 years were available in IOW cohort, which was used to evaluate biological relevance of CpGs shown to have time-lagged associations with lung function. We used paired-end (2 × 75 bp) RNA sequencing with the Illumina Tru-Seq Stranded mRNA Library Preparation Kit with IDT for Illumina Unique Dual Index (UDI) barcode primers following manufacturer's recommendations. RNA samples were extracted from whole blood of the IOW cohort participants at age 26 years. All samples were sequenced second time using the identical protocol, and for each sample, the output from both runs was combined. FASTQC were run to assess the quality of the FASTQ files [45]. Reads were mapped against Human Genome (GRch37 version 75) using HISAT2 (v2.1.0) aligner [46]. The alignment files, produced in the Sequence Alignment Map (SAM) format, were converted into the Binary Alignment Map (BAM) format using SAMtools (v1.3.1) [47]. HTseq (v0.11.1) was used to count the number of reads mapped to each gene in the same reference genome used for alignment [48]. Normalized read count FPKM (Fragments Per Kilobase of transcript per Million mapped reads) were calculated using the countToFPKM package (https ://githu b.com/AAlhe ndi17 07/count ToFPK M) and were included for subsequent data for analysis.

Confounders
Based on prior knowledge in the published literature, variables potentially associated with lung function trajectories in addition to DNA-M were included in the model as confounders. The potential confounders were birth weight, gestational age, duration of breastfeeding, maternal smoking exposure during pregnancy, recurrent chest infection collected at ages 1/2 years, second-hand smoking exposure at age 10 years (childhood), height and body mass index (BMI) at age 10 years, exposure to pets at age 10 years, age of puberty onset, and socioeconomic status (SES) [4,11,12].
Gestational age information was recorded during delivery in the hospital. Birth weight was measured immediately after birth. Heights and weights at age 10 were measured before spirometric measurements, and BMI was calculated from height and weight accordingly. The age of puberty onset was estimated based on self-reports about age of initiation of five pubertal markers for each sex, a growth spurt, body hair growth, skin changes, deepening voice of male, facial hair of male, breast development of female, and initiation of menstruation of female. The National Institute of Child and Human Development questionnaire from the Study of Early Child Care and Youth Development was used to identify pubertal stages. Information on second-hand smoke exposure in childhood was collected from parents. SES was classified using the composite "SES-cluster" method based on the following three variables: (a) the British socioeconomic classes [1][2][3][4][5][6] derived from parental occupation reported at birth; (b) the number of children in the index child's bedroom (collected at age 4 years); and (c) family income at age 10 years [49]. This composite variable captures the family social class across the entire study period. Pet exposure information was collected at age 10 years through questionnaires.

The ALSPAC cohort-Replication cohort
ALSPAC is a population-based birth cohort study established in 1991 in Avon, UK, approximately 75 miles from the IOW [50,51]. All pregnant women who were expecting to deliver between 1 April 1991 and 31 December 1992, and residing in the Avon region of the South West of England were eligible to be recruited. In total, 14,541 pregnant women were recruited for the study, of those 13,761 were eligible with 10,321 having DNA sampled. Information on environment, lifestyle, and health of the child and family was collected through annual questionnaires since the child's birth. At age 7 an additional 913 children were enrolled. The total sample size for analyses using any data collected after the age of 7 is therefore 15,454 pregnancies, resulting in 15,589 foetuses. Of these 14,901 were alive at 1 year of age. From age 7 years, the participants were invited to an annual research clinic to obtain the exposure and other demographic data annually. Spirometry (Vitalograph 2120; Vitalograph, Maids Moreton, UK) was performed at 8, 15, and 24 years of age according to ATS standards [39,52,53], the same method as that applied in the IOW cohort. The study website contains details of all the data that are available through a fully searchable data dictionary and variable search tool (http://www.brist ol.ac.uk/alspa c/resea rcher s/our-data/).
DNA-M data of children at ages 7 (n = 968) years were included in the study. DNA-M in peripheral blood was assessed using the Infinium HumanMethylation450K BeadChip [54]. The procedure for DNA sample preparation was comparable to that applied in the IOW cohort. The pre-processing of DNA-M was performed by adjusting batch effect, excluding CpGs with detection p value ≥ 0.01, and excluding samples that were flagged a sex-mismatch based on X-chromosome methylation. Details of pre-processing, quality control, and quantile normalization of DNA-M data have been described elsewhere [54,55].

Statistical analyses Descriptive analyses
To evaluate whether subjects included in the study reasonably represented those in the complete study cohort, we compared lung function tests at ages 10, 18, and 26 years in the study samples with those of the complete cohort using one-sample t tests.

Determining distinctive lung function trajectories
Our previous publication [12] of lung function trajectory was based on at least single spirometry test to attain a maximum sample size. In this study, subjects with at least two-time point tests were included for trajectory analyses to improve the average posterior probability and to avoid the random assignment of the subjects into a trajectory. An unsupervised semi-parametric mixture modelling implemented in the SAS procedure PROC TRAJ [56] was applied to identify developmental lung function trajectories of FVC, FEV 1 , and FEV 1 /FVC over time (10,18, and 26 years) for males and females separately [57], the same approach applied in our previous study [12]. This method combines the latent growth curve and mixture modelling approaches to detect distinct groups of trajectories [56]. All possible models were evaluated each with different numbers of groups (i.e. 2, 3, and 4) and different shapes of the trajectories (linear, quadratic, and cubic) for each group. Trajectory parameters were estimated using the maximum likelihood approach [58,59]. The best model was selected based on two criteria, being as parsimonious as possible to summarize the distinctive features and with high Bayesian information criterion (BIC) [57,60,61]. To improve the quality of identified trajectories, in addition to BIC, probability of trajectory assignment as well as sample sizes in each group was further incorporated; the average posterior probabilities of assignment to a group were set at least 0.7, and the sample size of each group was required to be at least 5% of the total sample size [60]. Individuals were assigned to one of the trajectories/groups based on their highest estimated group-membership probabilities. The assigned group (categorical variable) of distinct lung function trajectories was used in subsequent analyses.

Association analyses
To assess the association of DNA-M at an earlier age with lung function trajectories at later ages in the IOW cohort, we followed a two-step analytical approach. In the first step, CpGs were screened to exclude CpGs potentially not associated with lung function trajectories using ttScreening (R package 3.5.1 version) [62,63]. This method utilizes training and testing data in robust linear regressions with surrogate variables included in the regressions to adjust for unknown effects. The training and testing steps were repeated 100 times. A CpG that was statistically significant in both training and testing steps at least 50 times was included in the final set for subsequent regression analyses. The screening was performed for each lung function parameter, stratified by sex.
In second step, CpGs that passed screening were further assessed in logistic regression models in SAS 9.4 for the trajectories of each lung function parameter stratified by sex and adjusting for the above-mentioned confounders. Lung function trajectory was treated as the outcome variable, and the DNA-M at each CpG that passed screening was used as an independent variable. Multiple testing was corrected using the Bonferroni method with an experiment-wise significance level of 0.05. In all analyses, DNA-M adjusted for cell types, principle components, and batch effects at each CpG was used.

Replication analysis-in the ALSPAC cohort
The identified CpGs from the IOW cohort were further examined in the ALSPAC. Comparable analytical methods as those implemented in the IOW cohort were applied except for the exclusion of two covariates in regression analyses, recurrent chest infection, and pet exposure at childhood, which were unavailable in the ALSPAC.

Gene expression analysis
To assess the potential biological relevance of the identified CpGs, we examined the time-lagged association of DNA-M at those CpGs with the expression of their mapped genes. Linear regressions were applied to DNA-M at age 10 years with gene expression at age 26 years to each CpG which showed the consistent direction of association between the IOW cohort and ALSPAC.

Analyses of differentially methylated regions (DMRs)
To detect regional differential methylation signals among the CpGs that passed screening, an R package DMRcate was used [64]. The default settings in DMRcate include having at least two CpGs in a region and a minimum length of 1000 nucleotides. In our study, a DMR was considered to be statistically significant if the false discovery rate (FDR)-adjusted p value was < 0.05 [64]. Since DMR analyses focus on contribution of a region as a whole unit, a significant DMR can be detected even if there is no genome-wide significant individual CpGs in the region.

Pathway analyses
Genes annotated to the CpGs explored in ALSPAC with respect to the direction of odds ratios (ORs in the log scale) and DMRs were identified based on Illumina's manifestation file and SNIPPER (https ://csg.sph. umich .edu/boehn ke/snipp er/) version 1.2. Bioinformatic assessment of the genes was conducted using the online bioinformatics tool ToppFun, available in the ToppGene Suite [65]. Multiple testing was adjusted by controlling FDR of 0.05.