Epigenome-wide association study of leukocyte telomere length

Telomere length is associated with age-related diseases and is highly heritable. It is unclear, however, to what extent epigenetic modifications are associated with leukocyte telomere length (LTL). In this study, we conducted a large-scale epigenome-wide association study (EWAS) of LTL using seven large cohorts (n=5,713) – the Framingham Heart Study, the Jackson Heart Study, the Women’s Health Initiative, the Bogalusa Heart Study, the Lothian Birth Cohorts of 1921 and 1936, and the Longitudinal Study of Aging Danish Twins. Our stratified analysis suggests that EWAS findings for women of African ancestry may be distinct from those of three other groups: males of African ancestry, and males and females of European ancestry. Using a meta-analysis framework, we identified DNA methylation (DNAm) levels at 823 CpG sites to be significantly associated (P<1E-7) with LTL after adjusting for age, sex, ethnicity, and imputed white blood cell counts. Functional enrichment analyses revealed that these CpG sites are near genes that play a role in circadian rhythm, blood coagulation, and wound healing. Weighted correlation network analysis identified four co-methylation modules associated with LTL, age, and blood cell counts. Overall, this study reveals highly significant relationships between two hallmarks of aging: telomere biology and epigenetic changes.

Shown are results of the region-based (binomial) and gene-based (hypergeometric) functional enrichment test for all coagulation and wound-healing related annotations tested using the GREAT method. Differentially methylated probes were used as input and show a consistent signal of enrichment for LTL-associated CpG sites.

LTL-DNAm correlation in subtelomeric regions
Here, subtelomeric regions included both ends of chromosomes. The ends were proportional (5% for each) to the length of the corresponding chromosome (Supplementary Figure 2).

Supplementary Figure 2. Definition of subtelomeric regions.
We focused on the 823 significant CpG sites that were associated with the fully adjusted LTL. We counted the number of positive and negative Z-scores in nonsubtelomeric and subtelomeric regions. The proportion of the positive LTL-DNAm correlations was 9.9% in non-subtelomeric bodies, whereas the corresponding proportion was 17.1% in subtelomeric regions.

Non-subtelomeric
The chance of having significant CpGs was slightly higher in non-subtelomeric bodies than in subtelomeric regions (P=0.0427).

Supplementary Figure 3. Comparison of the BICOR and OLR method.
We conducted a sensitivity analysis to compare the biweight midcorrelation (BICOR) and ordinary linear regression (OLR) method for the EWAS of LTL. The panel above displays the Z scores generated by the BICOR and by the OLR method in the six strata (Supplementary Figure 3). The LTL was adjusted for age and the blood cell counts in each sex and ethnicity specific stratum as it was in the original analysis using the BICOR method. We replaced the BICOR with the OLR. For a clear presentation of results, we focused on the 823 significant CpG sites (fully adjusted) in this sensitivity analysis.
www.aging-us. We conducted a sensitivity analysis with and without additional adjustment for BMI using the three cohorts (FHS, JHS and WHI). Supplementary Figure 4 reveals that the Z scores adjusted for the existing variables (age, sex, ethnicity and blood cell counts) were almost same as the Z score adjusted for the existing variable and BMI. We did not observe any significant change.

Supplementary Figure 5. Sensitivity analysis with or without additional adjustment for BMI and education.
We conducted another sensitivity analysis with and without additional adjustment for BMI and education using FHS and WHI. The education attainment of the JHS subjects was missing. Again, we did not observe any change by the additional adjustment for BMI and education. We have also conducted a cohort-specific metaanalysis. We found the rough but consistent LTL-DNAm associations at the 823 CpG sites across the seven cohorts. Supplementary Figure 6 displays the Z scores from two cohorts using scatter plots (21 possible pairs/panels). The Stouffer"s method ( , where is the square root of the sample size in the th stratum) was used to combine corresponding strata for each cohort.

The Framingham Heart Study (FHS)
FHS started in 1948 to investigate common risk factors for cardiovascular disease (CVD) [2]. FHS recruited 5,209 subjects who lived in Framingham, Massachusetts, USA, and who were free from symptoms of CVD, heart attack or stroke at enrollment. The FHS Offspring Cohort started in 1971 and enrolled the original participants" grown-up children and the children"s spouses (n=5,124) [3]. Our study included 874 participants from the FHS Offspring Cohort who attended the sixth and eighth examination and consented to the use of their biospecimens for research purposes. Data from FHS can be retrieved from dbGaP (under accession numbers phs000363.v16.p10 and phs000724.v2.p9).
Telomere length measurement: DNA was extracted from leukocytes collected from the sixth examination cycle. Participants who had enough buffy coat available were selected for LTL measurement. The Southern blot method was used to obtain the mean length of the terminal restriction fragment (TRF) as previously described [4]. The coefficient of variation was 2.4% for the LTL measurement of duplicate and triplicate DNA samples.

The Jackson Heart Study (JHS)
JHS recruited 5,306 African Americans to investigate risk factors for cardiovascular disease in the Jackson metropolitan area, Mississippi, USA [7]. Participants provided medical and social records, physical and biochemical measurements, information on diagnostic procedures, and DNA samples during a baseline examination (2000)(2001)(2002)(2003)(2004), and two follow-up examinations (2005-2008 and 2009-2012). JHS follows up the participants every year and maintains cohort surveillance. Our study included the participants who visited at the baseline examination as part of the ancillary study ASN0104 in JHS.
Telomere length measurement: The Puregene kit (Gentra System, Minneapolis, MN, USA) [8] was used to extract DNA from whole blood. LTL (in kilobases) was measured using Southern blot [9]. The inter-assay coefficient of variation was 2.0%.
The interclass correlation coefficient was 0.95 for individual measures of LTL.

The Women's Health Initiative (WHI)
WHI started in 1992 and enrolled 64,500 postmenopausal women aged between 50 and 79 years into either clinical trials or observational studies [11]. Our study included WHI participants with available phenotype and DNA methylation array data referred to as the WHI "Broad Agency Award 23" (WHI BA23). WHI BA23 aimed to identify miRNA and genomic biomarkers of coronary heart disease.
Telomere length measurement: DNA was extracted from blood samples collected at the time of the 2012-2013 visit, using the 5-prime method (5 PRIME, Inc.; Gaithersburg, MD). Prior to LTL measurement, DNA integrity was assessed visually after ethidium bromide-stained 1% agarose gel electrophoresis (200 V for 2 hours). The Southern blot method was used to measure the average length of the terminal restriction fragments (in kilobases) [12]. Individual samples were measured in duplicate on different gels. The average inter-assay coefficient of variation was 2.0%.
DNA methylation: Among many sub-studies, only WHI BA23 provided both blood-based LTL and DNAm array data. WHI BA23 used the background correction method from Illumina"s proprietary GenomeStudio software.

The Bogalusa Heart Study (BHS)
BHS started in 1972 and recruited multiple waves of participants from childhood, adolescent and adulthood in a biracial community in Bogalusa, Louisiana, USA, comprising 65% whites and 35% African Americans [13]. The longitudinal cardiovascular risk factor phenotype and genotype data of the BHS cohort are available via application through the NHLBI Biologic Specimen and Data Repository Information Coordinating Center website (https://biolincc.nhlbi.nih.gov/studies/bhs).
Telomere length measurement: LTL was measured by Southern blot, as in the above studies. DNA was hybridized to a digoxigenin 3"-end labeled 5"-(CCCTAA)3 telomeric probe after overnight DNA digestion with 10 U Hinf I and 10 U Rsa I restriction enzymes as previously described [14]. Digitized autoradiograms of LTL measurement were analyzed for each sample resolved in duplicate on AGING different gels, and the coefficient of variation for the duplicate samples was 1.4% [14].
DNA methylation: Genomic DNA was isolated from whole blood samples in BHS using the FlexiGene DNA extraction kit (Qiagen, Hilden, Germany). The Infinium HumanMethylation450K BeadChip (Illumina, San Diego, California, USA) was used for whole-genome DNAm analysis. Samples were processed at the Microarray Core Facility Lab, University of Texas Southwestern Medical Center, Dallas, Texas, USA. For each subject, 750 ng of genomic DNA was bisulfiteconverted using the 96-well EZ DNAm kit (Zymo Research, Irvine, California, USA) according to the manufacturer"s instructions. The efficiency of the bisulfite conversion was confirmed by in-built controls on the 450K array. The methylation profile of each participant was measured by processing 4 μl of bisulfiteconverted DNA, at a concentration of 50 ng/μl, on an Illumina 450K array. The bisulfite-converted DNA was amplified, fragmented and hybridized to the array following the protocol. We scanned the arrays by using an Illumina iScan scanner, and then the raw methylation data was extracted using Illumina"s GenomeStudio Methylation (M) Module. BHS used the data-driven separate normalization method (dasen, [15]) from the R package wateRmelon [15]. The probe exclusion criteria for filtering samples and probes were: 1) samples having 1% of CpG sites with a detection p-value greater than 0.05; 2) probes having 5% of samples with a detection pvalue greater than 0.05; 3), and probes with bead count less than 3 in 5% of the samples.

The Lothian Birth Cohorts (LBC)
The Since then, extensive phenotypic data have been collected roughly every three years in four further waves of testing. The data collection includes detailed physical, cognitive, psychosocial and lifestyle measures. In addition, genetic and epigenetic data are available in both LBC21 and LBC36. More details on recruitment and testing can be found elsewhere [17,18]. Our study here uses data obtained in the first wave of testing.
Telomere length measurement: Telomere length in the LBC21 and LBC36 was measured in wave 1 using a quantitative real-time polymerase chain reaction (qPCR) assay [19]. DNA was extracted from whole blood at the Wellcome Trust Clinical Research Facility Genetics Core at the Western General Hospital in Edinburgh using standard procedures. A 7900HT Fast Real Time PCR machine with 384-well plate capacity (Applied Biosystems; Pleasanton, California, USA) was used to perform the PCRs. Telomere length was measured as the ratio of telomeric template to glyceraldehyde 3-phosphate dehydrogenase. Four internal control DNA samples derived from cell lines of known absolute telomere length were included on each plate to correct for plate-to-plate variation, and measurements were performed in quadruplicate and the mean was used in further assessments.
DNA methylation: LBC used internal controls from the R package minfi to correct for background noise. Following this, samples of low quality, e.g. those with bisulfite conversion, staining signal, inadequate hybridization or nucleotide extension, were excluded. In addition, probes with a detection rate <95% at p<0.01 and samples with a low call rate (<450,000 probes detected at p<0.01) were removed. Finally, samples for which DNA-methylation predicted sex did not match reported sex and samples that showed a poor match between SNP control probes and genotype were removed.

The Longitudinal Study of Aging Danish Twins (LSADT)
LSADT was initiated in 1995 and recruited all Danish twins aged 70 years or more [9,10]. Surviving twins were surveyed every other year until 2007. In 1997, whole-blood samples were collected from 689 same-sex twins. For 310 of these individuals, genome-wide DNAm data was available. The present study includes all twin pairs who participated in the 1997 wave and for whom genome-wide DNA methylation data and LTL measurements were available.
Telomere length measurement: LTL was measured using Southern blot. The LTL measurement referred to the average terminal restriction fragments after digestion with HinfI and RsaI restriction enzymes as previously described [9]. The two LTL measures presented a high correlation (r=0.88, p=0.000). Each sample was used in duplicate on different gels. The inter-assay coefficient of variation was 2.5% (for the HinfI/RsaI digest).

DNA methylation:
The EZ Methylation Gold kit (Zymo Research, Orange County, California, USA) was used to isolate DNA