Different epigenetic signatures of newborn telomere length and telomere attrition rate in early life

Study population

Table 1 shows the demographic characteristics of the mother-newborn (mother-child) pairs. The 247 newborns/children had a mean (standard deviation, SD) gestational age of 39.24 (1.37) weeks, and age of 4.58 (0.39) years at follow-up. The study population included 121 (49.0%) girls, and 93.5% of newborns were of European-origin. The mean (SD) cord blood telomere length, as a surrogate for baseline telomere length, was 1.17 (0.27) T/S ratio. The leukocyte telomere length of the children was 1.02 (0.19) T/S and on average 0.14 T/S shorter than cord blood telomere length (p<0.0001). The variations in cord blood telomere lengths and leukocyte telomere lengths were significantly different (p<0.0001). Telomere attrition corrected for regression-to-the-mean (RTM) effect was on average 0.0025 (0.14) T/S ratio, and telomere attrition rate was on average 0.00072 (0.030) T/S per year.

The DNAm-based explanatory model of baseline telomere length

Selection of CpGs

The study design is schematically shown in Figure 1. The epigenome-wide association study (EWAS) of cord blood telomere length, as a pre-screening for candidate CpG predictors from 787,264 CpG probes, identified 244 individual CpGs to be significant (FDR=5%) (Figure 2A). Two CpGs passed the p-value threshold under Bonferroni correction. A loose criterion of p-value<0.01 was applied to allow a sufficiently large number of CpGs in further selection, also considering that CpGs that jointly associate with baseline telomere length might not be the individually significant ones. Therefore, 22,817 CpGs passed the pre-screening and were entered in a subsequent selection in one and same adaptive lasso [31] model which accounts for the between-CpG correlations (Supplementary Table 1). 208 CpGs were jointly selected by adaptive lasso. In order to account for the CpG-CpG interactions, all possible pairwise interactions (21,528 interactions) were added, making it in total 21,736 terms subject to an elastic net selection [32]. In order not to overfit the data, we repeatedly sampled 1000 random subsets from the whole dataset, each containing 85% of the observations (n=210). Variation was observed among the 1000 subsets in terms of the elastic net regularization parameters and the number of selected predictor terms, as shown in Supplementary Figure 1A. Finally, 36 CpG main effects and 7 interactions involving 47 CpGs in total, which were selected in at least 95% of the subsets, were kept as the most relevant. We added the CpGs which were involved in the selected interactions but not in the selected main effects, and thus the final model contained 47 CpG main effects and 7 interactions (54 terms in total).

The final model and interpretations

The Pearson correlation between baseline telomere length and the weighted sum of the 54 terms was 0.89 and the R² was 0.76. Figure 3A and Supplementary Table 2 show the contributing weights (final model coefficients) and the EWAS coefficients of the CpGs and their interactions, as well as the genomic locations. 8 out of the 47 CpG main effects in the final model had been identified significant in EWAS. In the final model, 26 main effect CpGs had a positive coefficient (being positively associated with baseline telomere length) and 21 had a negative coefficient (being inversely associated with baseline telomere length), consistent with their EWAS coefficients. 15 CpGs were promoter-associated and 1 was related to differentially methylated regions (DMRs). The genomic locations of the involved 47 CpGs were distributed across autosomes except chromosomes 5, 9, 18, 20 and 21. All the 7 interaction terms involved CpGs from different chromosomes, among which 3 had positive coefficients indicating that the interactors are mutually promoting their associations with baseline telomere length. The other 4 pairs of interactions with negative coefficients indicated mutual suppression of the associations between their methylation level and baseline telomere length.

40 of the 47 CpGs had gene name annotations, from which the pathway analysis suggested one significantly enriched biological pathway “translocation of SLC2A4 (GLUT4) to the plasma membrane” (Figure 4, left panel). Correlating the selected CpGs to RNA probes from the global gene expression microarray suggested 75 significantly correlated RNAs (p<0.05), which had no biological pathways significantly enriched.

The DNAm-based explanatory model of telomere attrition rate in early childhood

Selection of CpGs

Three CpGs were identified to associate with telomere attrition rate based on an FDR of 5% and are highlighted in Figure 2B. Only cg26557617 passed Bonferroni correction. 9,862 CpGs (p<0.01) were used for subsequent selection. 247 CpGs were then selected by an adaptive lasso model, based on which there were 30,381 pairwise interactions, summing up to 30,628 terms. The most relevant terms selected using an elastic net in at least 95% of the subsets were 25 CpGs and 5 interactions. Together with 6 CpGs which were involved in the selected interactions but not in the selected main effects, the final model contained 31 CpG main effects and 5 interactions. No overlap was observed between these terms and the model terms of baseline telomere length.

The final model and interpretations

The Pearson correlation between telomere attrition rate and the weighted sum of the 36 terms was 0.87 and the R² was 0.72. As is shown in Figure 3B and Supplementary Table 3, in the final model 18 of the 31 CpG main effects had telomere attrition rate-related increase in DNA methylation and 13 had inverse associations, also consistent with their EWAS associations. 7 CpGs were promoter-associated and 1 was related to a reprogramming-specific DMR (RDMR). 4 out of the 5 interactions had negative coefficients and the other had a positive coefficient.

25 CpGs in the final model were annotated with gene names but no pathway was found significantly enriched. Pathway analysis of the 54 RNA transcripts correlated to the selected 31 CpGs resulted in 30 enriched pathways (Figure 4, right panel) which are mostly related to signaling pathways in immune response.

Study population

The ongoing prospective ENVIRONAGE birth cohort [30] was initiated in 2010. This cohort study has been approved by the ethical committee of Hasselt University and East-Limburg Hospital (Genk, Belgium) (EudraCT B37120107805), and is carried out according to the Declaration of Helsinki.

The mother-newborn pairs were recruited at delivery in the East-Limburg Hospital and were invited for a follow-up visit when the children were around 4 years old. Cord blood samples were collected at birth and blood samples of the children were collected at the follow-up visit. Genome-wide DNA methylation profiled in cord blood was available in 372 samples. Among the 423 children attending the follow-up examination between October 2014 and October 2019, 247 children had matching cord blood (baseline) telomere length, leukocyte telomere length (follow-up) and cord blood DNA methylation. These samples were used to evaluate the relationship between DNA methylation and baseline telomere length or longitudinal change in telomere length from birth to childhood. The selected 247 participants are representative for the total ENVIRONAGE population in terms of most characteristics as described in Table 1, but slightly differed in maternal education and newborn ethnicity (Supplementary Table 4).

Written informed consent was obtained from the participating mothers and questionnaires were filled out by the mothers to collect information on lifestyles and socioeconomic status. The information on newborns’ sex, birth weight and gestational age, and maternal age and parity was collected from medical records in the hospital. The date of conception was estimated based on the first day of the mother’s last menstrual period in combination with the first ultrasonographic examination. Maternal body mass index (BMI) was determined as the ratio of the maternal weight to the squared maternal height measured at the first antennal visit (weeks 7-9 of pregnancy). Ethnicity of a newborn was categorized as European if at least 2 grandparents were Europeans and classified as non-European otherwise. Educational level of the mothers was coded as 0 if they did not obtain any diploma, as 1 if they obtained at highest a high school diploma, as 2 if the highest diploma obtained was from a 3-year college and as 3 if obtained a 4-year college or university degree. Maternal smoking status was classified into “never smoker”, “former smoker” (when having quit smoking before pregnancy), or “smoker” (in case of continuing smoking during pregnancy). Mothers were regarded as with pregnancy complications if they had any of gestational diabetes, hypertension, hyper- or hypothyroidism, infectious disease, preeclampsia, vaginal bleeding, phenylketonuria and allergy or asthma during pregnancy.

Blood sample collection and DNA extraction

Cord blood was collected immediately after delivery in BD Vacutainer® plastic whole blood tubes with spray-coated K2EDTA (BD, Franklin Lakes, NJ, USA). Samples were centrifuged at 3,200rpm for 15 minutes. Plasma was separated and the remainder with the buffy coats was stored at -80° C until further analysis. At the follow-up visit, blood was collected from the participating children by venipuncture using a winged steel needle in spray-coated K2EDTA tubes. After the sample process as described for cord blood, samples were stored at -80° C until analysis.

Cord blood and child blood DNA was extracted using the QIAamp DNA mini kit (Qiagen, Inc., Venlo, The Netherlands). DNA purity and concentration were verified by Nanodrop 1000 spectrophotometer (Isogen, Life Science, Belgium). DNA was considered pure when the A260/280 was greater than 1.80 and A260/230 greater than 2.0. DNA integrity was assessed with agarose gel electrophoresis.

Average relative telomere length measurement and data processing

Average relative telomere length was measured in triplicates using a previously described quantitative, real-time polymerase chain reaction (qPCR) protocol [14]. Detailed assay specifications are provided in the Supplementary Method. Telomere length were expressed as the ratio of telomere copy number to single-copy gene number (T/S) relative to the average T/S ratio of the entire sample set. Assay-reliability was assessed using an intra-class correlation coefficient (ICC) [56]. The inter-assay ICC was 0.936 (95%CI: 0.808 to 0.969) and the intra-assay ICC was 0.952 (95%CI: 0.947 to 0.956).

The cord blood telomere lengths were measured along with the leukocyte telomere lengths at follow-up in order to ensure that matching samples of individuals were assayed under the same condition: samples from the same individuals were on the same qPCR plate. Having been measured in a single batch, the telomere length data did not undergo batch effect removal. Telomere attrition from birth to follow-up was defined as the difference between telomere length at birth and at follow-up, the latter subtracted from the former. The difference between telomere lengths at birth and at follow-up was assessed using a paired t-test. Whether the children had differential telomere length change (homogeneity) was tested using the test-statistic suggested by Berry et al [54]. To account for the RTM effect which might be a result of the variation in measurement, we adopted a previously published method [54] to obtain the RTM-corrected telomere attrition:

$\Delta TL=\rho \theta (X_1-{\mu }_1)-(X_2-{\mu }_2)$ (1)

where X₁ and X₂ are the notations for telomere length at birth and at follow-up, respectively, $\widehat{\rho }=r$is the Pearson correlation between the telomere lengths measured at the two time points, and $$\widehat{\theta }=\frac{S_2}{S_1}$$is the ratio of the standard deviation of telomere lengths at follow-up to the standard deviation of baseline telomere lengths. ${\widehat{\mu }}_1=\overline{X_1}$ and ${\widehat{\mu }}_2=\overline{X_2}$ denote the mean baseline telomere length and mean leukocyte telomere length at follow-up, respectively. Telomere attrition rate was calculated as telomere attrition divided by age (years) at follow-up.

DNA methylation measurement and data processing

Cord blood DNA samples were bisulphite-converted, amplified and hybridized to the Illumina HumanMethylationEPIC Bead-Chip array (Illumina, San Diego, CA, USA), at the service lab GenomeScan (Leiden, The Netherlands). The array measurements were scanned using an Illumina iScan and the data quality was assessed using the R script MethylAid. DNA methylation data were preprocessed using the minfi package in R [57]. Briefly, probes with a call rate lower than 95% based on detection p-value > 10e-16 [58], samples with a call rate lower than 98% and samples with discordant sex, predicted using shinyMethyl [59], were removed. Methylation data was normalized using functional normalization [60]. For each CpG locus the methylation level was expressed as beta value calculated by M/(M+U) where M stands for the signal intensity from methylated probes and U stands for the signal intensity of unmethylated probes. The preprocessing resulted in a DNA methylation data with 857,898 CpGs. The missingness in the matrix of beta values was imputed by K-nearest neighbor (KNN) imputation and technical confounding effects (batch and position) were removed from the beta matrix using an empirical Bayes method [61]. Subsequently, we trimmed the data CpG-wisely for outliers which were lower than 3 inter-quarter-range below the 1^st quartile, or higher than 3 inter-quarter-range above the 3^rd quartile. CpG probes were filtered to exclude CpGs on X and Y chromosomes, those known to be common SNPs and those having cross-reactivity with multiple genomic locations [51]. 787,264 CpGs remained in the DNA methylation data were used in the EWAS of telomere lengths.

Based on the DNA methylation data, blood cell proportions (nucleated red blood cells, granulocytes, monocytes, natural killer cells, B cells, CD4⁺ T cells, and CD8⁺ T cells) in the cord blood samples were estimated using Bakulski algorithm [62].

Pre-screening by EWAS

The design of the present study is illustrated in Figure 1. Since variable selection from the whole EPIC array is computationally expensive, we performed EWAS as a pre-screening for candidate CpGs. For the j-th CpG (j=1, 2, …, 787264), a multiple linear regression model was fitted to test the association between baseline telomere length (Y_i, i=1, 2, …, 247) and the DNA methylation level at the CpG (X_ij), with robust covariance matrix estimated using the R package sandwich [63, 64]:

$Y_i={\beta }_0+{\beta }_m{X}_{ij}+C_i^T{\beta }_c$ (2)

CpGs entered models in their original scale without extra centering or standardization. Models were adjusted for covariates selected a priori (C_i): newborn gestational age, sex, ethnicity and birth weight, maternal pre-pregnancy BMI, age, parity, medical condition during pregnancy and education, paternal age, and cell type heterogeneity. Similarly, an EWAS was performed for telomere attrition rate by fitting the same linear model as equation (2) for each CpG adjusting for the same set of covariates. The test-statistics of β_m were corrected for an inflation factor and bias estimated based on their empirical null distribution [65], hence giving the corrected p-values. The Benjamini-Hochberg adjusted p-values were calculated based on the corrected p-values to identify significant EWAS probes (FDR=5%). In order to select a sufficiently large number of CpGs for the subsequent selection procedure, CpGs whose inflation- and bias-corrected p-value was smaller than a pre-determined threshold of 0.01 were the candidates that become the inputs of the next selection procedure.

Variable selection by statistical learning

The variable selection in this study was within the framework of elastic net where CpG selection and estimation of the contributing weight of each predictor take place simultaneously. We aimed to build explanatory models that express the baseline telomere length or telomere attrition rate as a weighted sum of the methylation level of a set of CpGs. The variable selection method used for both outcomes were identical, therefore we only describe the method in detail for baseline telomere length here.

The aim of this study was to select a number of CpGs (X_k) and CpG-CpG interactions (X_sX_t) that can express the i-th newborn’s baseline (cord blood) telomere length (Y_i) as:

$Y_i=w_0+{\displaystyle \sum }_k{w}_k{X}_{ik}+{\displaystyle \sum }_{s,t}{v}_{s,t}{X}_{is}{X}_{it}$(3)

where w₀ is the model intercept, being the mean of Y_i due to the standardization of CpGs. Parameters w_k and v_s,t denote the contributing weights of the CpGs and interaction terms in the model. The model was trained with the matrix of candidate CpGs as the inputs and the vector of baseline telomere lengths as the output. In order to reduce the influence by the difference between CpG’s magnitude and possible collinearity, the CpGs were standardized to have a mean of 0 and standard deviation of 1 before entering the model. When interactions are subject to selection, the pairwise interaction of CpGs were made by multiplying pairwisely the standardized CpGs.

The elastic net regression minimizes the objective function (simplified to include only main effects)

$\begin{array}{l}\frac{1}{2N}{\displaystyle \sum _{i=1}^N{(Y_i-\underset{k}{{{\displaystyle \sum }}^{\text{}}}{w}_k{X}_{ik})}^2}\hfill \\ +\left[\lambda \frac{1}{2}(1-\alpha )\underset{k}{{{\displaystyle \sum }}^{\text{}}}{w}_k{}^2+\alpha \underset{k}{{{\displaystyle \sum }}^{\text{}}}\left|w_k\right|\right]\hfill \end{array}$ (4)

where N is the number of observations. The second term is the penalty on coefficient estimates of variables (w_k). The regularization parameter λ was determined by 10-fold cross validation based on model deviance. In order to balance between model predictability and model overfitting, the λ value corresponding to a 1σ deviation (λ_1σ) from λ=λ_min was chosen, which led to a slightly smaller number of predictors than with the λ_min. The tuning parameter α between lasso (α=1) and ridge (α=0) regression was determined by computation over a sequence of α value ranging from 0 to 1 with step of 0.05 and the α with the smallest deviance was selected. The models were fit in R with package glmnet [66], where coordinate descent algorithm was used and convergence threshold was set to 10^-7.

Adaptive lasso, a consistent penalized model, variable selection, was used to select candidate CpGs in a multiple linear model including all 22,817 CpG’s selected from pre-screening by EWAS. The adaptive lasso performs lasso regression (α=1) by shrinking coefficients of less relevant variables to zero, while at the same time assigning unequal penalties to each variable [31]. The penalty term shown in equation (4) for adaptive lasso was λ∑_k a_k |w_k|, where a_k was the penalty weight for the k-th CpG and was determined in a data-driven way as the inverse of the absolute value of the ridge regression coefficient of the corresponding CpG, as the ridge regression accounts for correlations among variables. By doing so, the CpGs given more shrinkage in ridge regression were assigned even higher penalty in the adaptive lasso variable selection. 208 CpGs were selected by adaptive lasso, which forms 21,528 pairwise interactions. Thus, a total number of 21,736 terms underwent the elastic net selection.

To encounter the issue that one of the components of elastic net, the lasso penalty, is not statistically consistent, we used 1000 repeatedly randomly sampled subsets from the whole data set for model selection. During each round, 85% of the samples (n=210) were taken and used to fit an elastic net, for which the parameters α, λ_1σ and the number of predictors corresponding to the selected model were recorded. Among all the CpG and interactions, those has been selected in no less than 95% of the 1000 subsets were regarded as active terms and put in the final model. The weights w_k in equation (3) were estimated by fitting a ridge regression model on all the selected CpGs and interactions to the whole data set. Since statistical inference was not available within the ridge regression framework, the p-values of the ridge coefficients were simulated via bootstrap with 10,000 replicates.

Functional annotation of the selected CpGs

The selected CpGs were annotated with genomic regions (promoter region or DMRs) and UCSC reference gene names. Chord diagram visualization was performed using R package circlize [67].

To assist the translation of the identified CpGs into biological pathways, the DNA methylation level of the CpGs were correlated to the gene expression data available in the ENVIRONAGE birth cohort. Global gene expression profiling by microarrays was performed in cord blood samples using a method described previously [68]. 44 samples which were available in both DNA methylation and global gene expression profiles were used in searching for RNA probes that are correlated with the selected CpGs. The correlation was assessed by computing the Pearson correlation coefficient within each CpG-RNA pair, which were defined based on gene name annotation or genomic location. Among the 47 CpGs explaining baseline telomere length, 36 CpGs found RNA probes with the same gene name annotated. Additionally, each of the 47 CpGs was correlated to the RNA transcripts on the same chromosome whose starting sites were within 1Mb from the CpG locus. In total, 1,311 correlations were calculated and 75 correlations had p-value<0.05. The gene names corresponding to the 75 correlated RNA transcripts were used for a pathway analysis. Similarly, there were 676 RNA transcripts correlated to the 31 CpGs explaining telomere attrition rate: 22 were found based on gene name and 654 were found within genomic neighborhood. 53 correlations with p-value < 0.05 were used for a pathway analysis.

The gene names of selected CpGs and the correlated transcripts, respectively, were mapped to the Reactome pathway database [69] using the R package ReactomePA [70] for an over-representation analysis (ORA). A pathway was considered significantly enriched if its q-value was smaller than FDR of 5% and at least 3 genes were involved in the pathway.

Ethics approval and consent to participate

The study protocol was approved by the Ethical Committee of Hasselt University and East-Limburg Hospital in Genk (Belgium) and has been carried out according to the declaration of Helsinki. Written informed consent was obtained from all participating mothers.

Data availability

Data are available on request.