Research Paper Volume 12, Issue 12 pp 11942—11966
Biomarkers of aging and lung function in the normative aging study
- 1 Department of Environmental Health, Harvard T.H. Chan School of Public Health, Boston, MA 02115, USA
- 2 Department of Environmental Medicine and Public Health, Icahn School of Medicine at Mount Sinai, New York, NY 10029, USA
- 3 Department of Biostatistics, Harvard T.H. Chan School of Public Health, Boston, MA 02115, USA
- 4 Department of Preventive Medicine, Feinberg School of Medicine, Northwestern University, Chicago, IL 60611, USA
- 5 VA Normative Aging Study, VA Boston Healthcare System, Boston, MA 02130, USA
- 6 Department of Medicine, Boston University School of Medicine, Boston, MA 02118, USA
- 7 Department of Epidemiology and Environmental Health Sciences, Columbia University, New York, NY 10027, USA
Received: February 20, 2019 Accepted: May 20, 2020 Published: June 19, 2020https://doi.org/10.18632/aging.103363
How to Cite
Copyright © 2020 Wang et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Elderly individuals who are never smokers but have the same height and chronological age can have substantial differences in lung function. The underlying biological mechanisms are unclear. To evaluate the associations of different biomarkers of aging (BoA) and lung function, we performed a repeated-measures analysis in the Normative Aging Study using linear mixed-effect models. We generated GrimAgeAccel, PhenoAgeAccel, extrinsic and intrinsic epigenetic age acceleration using a publically available online calculator. We calculated Zhang’s DNAmRiskScore based on 10 CpGs. We measured telomere length (TL) and mitochondrial DNA copy number (mtDNA-CN) using quantitative real-time polymerase chain reaction. A pulmonary function test was performed measuring forced expiratory volume in 1 second / forced vital capacity (FEV1/FVC), FEV1, and maximum mid-expiratory flow (MMEF). Epigenetic-based BoA were associated with lower lung function. For example, a one-year increase in GrimAgeAccel was associated with a 13.64 mL [95% confidence interval (CI), 5.11 to 22.16] decline in FEV1; a 0.2 increase in Zhang’s DNAmRiskScore was associated with a 0.009 L/s (0.005 to 0.013) reduction in MMEF. No association was found between TL/mtDNA-CN and lung function. Overall, this paper shows that epigenetics might be a potential mechanism underlying pulmonary dysfunction in the elderly.
In 2050, there will be 84.7 million people ≥ 65 years of age in the U.S., almost double its estimated population of 43.1 million in 2012 . The aging process reduces physiological capacity, which can result in functional impairment (e.g., lower lung function), chronic disease, and mortality. Chronological age is undoubtedly a major risk factor for aging-related diseases and death . However, there is still great heterogeneity in the health outcomes of older individuals who have the same chronological age , especially for lung function, which has substantial heterogeneity among elderly individuals who are never smokers but have the same height and chronological age. These different vulnerabilities to age-related diseases and death are likely reflective of differences in their underlying biological aging processes . The prevalence of pulmonary dysfunction increases sharply with age among elderly people. Moreover, poor lung function is related to impaired cognition [5, 6], adverse findings on brain imaging , and increased mortality [8–11]. However, the underlying mechanisms of lung function decline with aging are not fully characterized.
In the past decade, a number of biomarkers of aging (BoA) have emerged. For example, telomere length (TL) [12–14] and mitochondrial DNA copy number (mtDNA-CN) [15, 16] are well-known BoA and are related to aging-related diseases. As increasing body of evidence has suggested associations between aging and the epigenome , with DNA methylation (DNAm) levels being used to generate composite BoA [18–21]. The first generation of epigenetic BoA was calculated from DNAm levels alone, such as the Hannum epigenetic clock based on 71 5'-C-phosphate-G-3' sites (CpGs) in leukocytes , and the Horvath epigenetic clock based on 353 CpGs in multiple tissues . Recently, the second generation of epigenetic BoA has been developed incorporating additional markers. For example, Levine et al.  created DNAm PhenoAge, based on 513 CpGs selected by predicting phenotypic age that were constructed using 9 clinical biomarkers (e.g., albumin). Lu et al.  generated DNAm GrimAge from 1,030 CpG sites that were highly predictable of 7 plasma proteins (e.g., DNA leptin) and smoking pack-years. Further, refined measures of BoA-epigenetic age acceleration were shown to be associated with age-related diseases and mortality [4, 20, 22]. The widely used BoA on epigenetic age accelerations include intrinsic and extrinsic epigenetic age acceleration (IEAA, EEAA) , PhenoAgeAccel , and GrimAgeAccel . Additionally, Zhang’s DNAmRiskScore, which was calculated based on 10 CpG sites , was shown to be strongly associated with all-cause mortality . The relationship between TL and lung function has been reported to be positive or null [23, 24]. Meanwhile, mtDNA-CN reduction is associated with chronic obstructive pulmonary disease . Even epigenetic modifications have been shown to be associated with lung function in the elderly [26–28], however, no studies have investigated the associations between epigenetic aging biomarkers and lung function.
In this present study, we hypothesized that some of these BoA are associated with lower lung function. To examine this hypothesis, we tested for associations between each of seven BoA and three measures of lung function [forced expiratory volume in 1 second (FEV1), forced expiratory volume in 1 second / forced vital capacity (FEV1/FVC), maximum mid-expiratory flow (MMEF)] in the Normative Aging Study. Lower lung function in this paper is defined as a decrease in any of these three studied pulmonary function parameters.
The present study included 696 elderly men with 1,070 visits during years of 1999-2013. Of the 696 subjects, 345 (50%) had two follow-up visits, and 29 (4%) completed three follow-up visits. Table 1 shows the personal characteristics of the participants across visits. The mean chronological age was 72.5, 75.2, and 76.3 in the first, second, and third visit, respectively. The cohort was well educated and made up of one third never smokers. Summary statistics of seven different BoA are presented in Table 2. The means of four age accelerations were close to zero, with GrimAgeAccel had the smallest standard deviation (4 years) (see Table 2).
Table 1. Personal characteristics of 696 (N=1,070) elderly white men from the Normative Aging Study, 1999-2013.
|Variable||First visit (N=696)||Second visit (N=345)||Third visit (N=29)||All visits (N=1,070)|
|Chronological age (years)|
|Mean ± SD||72.51 ± 6.74||75.20 ± 6.39||76.34 ± 5.94||73.48 ± 6.74|
|BMI (kg/m2), mean ± SD||28.09 ± 4.10||27.81 ± 4.13||27.96 ± 4.32||28.00 ± 4.11|
|Height (m), mean ± SD||1.74 ± 0.07||1.74 ± 0.07||1.74 ± 0.07||1.74 ± 0.07|
|Smoking status, n (%)|
|Never smokers||216 (31.03)||112 (32.46)||13 (44.83)||341 (31.87)|
|Current or former smokers||480 (68.97)||233 (67.54)||16 (55.17)||729 (68.13)|
|Pack-year smoked (years), mean ± SD||21.33 ± 25.33||19.07 ± 23.60||14.88 ± 20.18||20.43 ± 24.68|
|Years of education, mean ± SD||14.98 ± 2.97||15.17 ± 3.18||15.24 ± 3.45||15.05 ± 3.05|
|Corticosteroids, n (%)||51 (7.33)||30 (8.70)||1 (3.45)||82 (7.66)|
|Estimated cell types (%) |
|Monocytes||10.7 (1.7-25.1)||10.1 (4.0-18.7)||9.5 (6.4-12.6)||10.5 (1.7-25.1)|
|B cells||1.5 (0.0-31.4)||1.2 (0.0-17.5)||0.6 (0.0-5.3)||1.4 (0.0-31.4)|
|CD4+ T lymphocytes||8.7 (0.0-22.6)||8.3 (0.0-24.1)||8.5 (0.9-20.1)||8.6 (0.0-24.1)|
|CD8+ T lymphocytes||8.7 (0.0-20.2)||8.5 (0.0-15.5)||8.0 (3.0-13.8)||8.6 (0.0-20.2)|
|Natural killer cells||7.3 (0.0-22.1)||7.5 (0.0-25.6)||7.8 (1.6-21.6)||7.4 (0.0-25.6)|
|FEV1 (mL 1st sec), mean ± SD||2500 ± 580||2540 ± 610||2680 ± 440||2520 ± 586|
|FEV1/FVC, mean ± SD||0.75 ± 0.08||0.74 ± 0.07||0.73 ± 0.08||0.75 ± 0.08|
|MMEF (L/s), mean ± SD||0.24 ± 0.11||0.23 ± 0.10||0.22 ± 0.08||0.24 ± 0.11|
|Abbreviations: SD = standard deviation; BMI = body mass index; FEV1 = forced expiratory volume in 1 second; FVC = forced vital capacity; MMEF = maximum mid-expiratory flow.|
Table 2. Summary statistics of different biomarkers of aging in 1,070 visits among 696 elderly white men from the Normative Aging Study, 1999-2013.
|Variables||mean ± SD||Percentile|
|GrimAgeAccel (years)||0.01 ± 4.07||-5.52||-0.63||7.82|
|PhenoAgeAccel (years)||-0.17 ± 5.99||-9.31||-0.56||10.09|
|IEAA (years)||-0.16 ± 5.08||-7.87||-0.62||8.26|
|EEAA (years)||-0.10 ± 5.81||-9.41||-0.24||9.61|
|Zhang’s DNAmRiskScore||-1.87 ± 0.44||-2.54||-1.90||-1.07|
|TL||1.25 ± 0.49||0.58||1.19||2.15|
|mtDNA-CN||1.03 ± 0.22||0.73||1.01||1.38|
|Abbreviations: SD = standard deviation; IEAA = intrinsic epigenetic age acceleration; EEAA = extrinsic epigenetic age acceleration; TL = Telomere length; mtDNA-CN = mitochondrial DNA copy number.|
The trends of lung function and methylation based BoA among 29 men who had three visits are presented in supplementary material (Supplementary Figures 1 and 2). The supplementary material also contains trends of four DNAm ages (i.e. DNAm GrimAge, DNAm PhenoAge, the Horvath’s Clock, and the Hannum’s Clock), which increased consistently across visits (Supplementary Figure 3). We also generated four DNAm ages (i.e., DNAm GrimAge, DNAm PhenoAge, the Horvath’s clock, and the Hannum’s clock) on the Age Calculator website. Summary statistics of chronological age and four epigenetic ages are presented in Supplementary Table 1. Supplementary Figure 4 shows the distribution of chronological age and four DNAm ages, among which the distribution of DNAm GrimAge was closest to the one of chronological age (Supplementary Figure 4).
We calculated pairwise correlations between chronological age and the seven BoA (Figure 1). According to Hung et al. , only correlation coefficient (r) 0.30 should be considered as "clinically significant" when P < 0.05. In this study, chronological age was not clinically associated with seven BoAs, while it was statistically significantly correlated with Zhang’s DNAmRiskScore (correlation coefficient ® = 0.23; P ≤ 0.001), TL (r = -0.09; P ≤ 0.01), and mtDNA-CN (r = -0.13; P ≤ 0.01). GrimAgeAccel was clinically associated with PhenoAgeAccel (r = 0.33; P ≤ 0.001) and Zhang’s DNAmRiskScore (r = 0.54; P ≤ 0.001). IEAA and EEAA were also highly correlated with each other (r = 0.45; P ≤ 0.001). TL and MtDNA-CN was not clinically associated with other BoAs. We also computed pairwise correlations between chronological age and four types of epigenetic ages. All of them were highly correlated with each other (Supplementary Figure 5). For epidemiologic purposes, non-clinically significant associations may still represent important associations.
BoA and lung function
Figure 2 shows the results of the primary analyses. FEV1, FEV1/FVC and MMEF was significantly negatively associated with each of GrimAgeAccel and Zhang’s DNAmRiskScore. For example, each increment of one year in GrimAgeAccel was associated with a decrease in FEV1 of 13.64 mL (95% CI: -22.16 to -5.11; FDRB-H = 0.006), in FEV1/FVC of 0.002 (95% CI: -0.003 to -0.001; FDRB-H = 0.022), and in MMEF of 0.004 L/s (95% CI: -0.005 to -0.002; FDRB-H < 0.001). In addition, a 0.2 increase in Zhang’s DNAmRiskScore was associated with a decrease in FEV1 of 32.01 mL (95% CI: -51.26 to -12.75; FDRB-H = 0.006), in FEV1/FVC of 0.004 (95% CI: -0.007 to -0.001; FDRB-H = 0.037), and in MMEF of 0.009 L/s (95% CI: -0.013 to -0.005; FDRB-H < 0.001). IEAA was related to lower FEV1/FVC. No association was found between lung function and non-epigenetic BoA.
Figure 2. Associations between Seven BoAs and lung function for 696 men (1,070 visits), the NAS, 1999-2013. For GrimAgeAccel, PhenoAgeAccel, TL, and mtDNA-CN, we adjusted for chronological age, BMI, height, smoking status, cigarette pack-years, years of education, corticosteroid use, estimated cell types, and batch effects. For EEAA and IEAA, we adjusted for chronological age, BMI, height, smoking status, cigarette pack-years, years of education, corticosteroid use, and batch effects. For Zhang’s DNAmRiskScore, we adjusted for chronological age, BMI, height, smoking status, cigarette pack-years, years of education, corticosteroid use, estimated cell types, batch effects, and technical covariates: Non polymorphic Red, Specificity I Red, Bisulfite Conversion I Red, Bisulfite Conversion II, Extension Red. Abbreviations: IEAA = intrinsic epigenetic age acceleration; EEAA = extrinsic epigenetic age acceleration; TL = Telomere length; mtDNA-CN = mitochondrial DNA copy number; BoA = biomarkers of aging; BMI = body mass index; FDRB-H = Benjamin-Hochberg false discovery rate.
Sensitivity analyses further validated the reliability of our primary findings. Results were remarkably consistent after excluding visits at which corticosteroid use was reported (N=82) (Supplementary Figure 6) or adjusting for coronary heart disease (31% of visits), stroke (8%), and diabetes (15%) (Supplementary Figure 7). When we used inverse probability weighting to account for potential selection bias, the associations were similar to those from the primary analysis (Supplementary Figure 8). Subset analyses suggested strong negative associations of lung function and BoA (e.g., GrimAgeAccel, and Zhang’s DNAmRiskScore) among smokers, whereas almost no associations between BoA and lung function among never smokers were present (Supplementary Figure 9A and 9B).
In this longitudinal cohort of 696 elderly males, we found that GrimAgeAccel and Zhang’s DNAmRiskScore were associated with lower lung function, including FEV1, FEV1/FVC, and MMEF. These effects were stable in three sensitivity analyses. These results suggest that epigenomic variation might shed new insights into the pathogenesis of lower lung function with age. To the best of our knowledge, this is the first study to examine the associations between epigenetic BoA and lung function. Non-epigenetic aging biomarkers, including TL and mtDNA-CN were not related to pulmonary function.
According to the American Lung Association, the lung matures by age 20-25 years, and its function declines gradually after the age of 35 . Poor lung function is often indicative of risk of diseases [5–7] and mortality [8–11], and hence a public health priority. In this study, we focused on three spirometric indices that provide slightly different information on lung function. FEV1 is a measure of overall lung function. FEV1/FVC often decreases due to changes in elasticity, airspace enlargement, and other physiologic change that occur with age . In this study, we found that mean FEV1/FVC was high at 0.75 in the first visit and decreased slightly in the next two visits, which was consistent with our previous study . The potential reasons were: 1). the subjects in NAS were free of known chronic medical conditions at enrollment; 2). they have high educational attainment and relative healthy lifestyle. MMEF, expressed as forced expiratory flow between 25-75% of FVC, is determined by the size of the mid-airways. A low flow of rate of MMEF reflects airway narrowing or obstruction. In the present study, TL and mtDNA-CN were not related to lung function, whereas GrimAgeAccel and Zhang’s DNAmRiskScore were significantly and stably associated with lower lung function, including FEV1, FEV1/FVC, and MMEF. IEAA was significantly related to lower FEV1/FVC.
In keeping with the unprecedented growth rate of the world’s aging population, there is a clear need to better understand the biological aging process. TL shortens with each cell division , and thus serves as a measure of biological aging. Although a number of studies showed TL was associated with several aging-related diseases [14, 15], we didn’t find the relationship between TL and lower lung function in the main analyses. Consistent with our finding, Andujar et al.  found no association between TL at baseline and FEV1 decline among 448 middle-aged adults after 11 years of follow-up in Europe. It could be because telomere attrition does not have marked effects on cell physiology until a critical TL is reached , at which point the cell becomes senescent , or this result may be specific to lung function.
Mitochondria are double-membrane-bound organelles, present in almost all mammalian cells. Although most of a cell’s DNA is contained in the cell nucleus, the mitochondrion has its own independent genome, which is more susceptible to oxidative attack due to lack of introns and protective histone proteins as well as limited capacity to repair . The most important roles of mitochondria are the generation of adenosine triphosphate and regulation of oxidative stress . Decrease in mtDNA-CN, which serves as a surrogate marker of mitochondrial function, has been linked with diseases [15, 16]. However, there are differences among tissues in mitochondrial number, function, protein composition , which might partly explain why the present study did not find an association between leucocyte mtDNA-CN and lung function.
Several studies have suggested that epigenetic mechanisms like DNAm may provide explanation for lower lung function [26–28]. For example, Carmona et al.  found a positive association between DNAm of aryl-hydrocarbon receptor repressor (AHRR) gene and FEV1, and MMEF in the NAS and the Cooperative Health Research in the Region of Augsburg cohort. Lepeule et al.  showed that decreased methylation of genes carnitine O-acetyltransferase, coagulation factor-3, and toll-like receptor 2 was associated with lower lung function in the NAS.
Over the last six years, several epigenetic BoA have been shown associated with health [35, 36, 38, 39]. The first generation of epigenetic BoA was developed to predict chronological age (correlation coefficients > 0.9), such as the Hannum’s clock, which was based on 71 CpG sites from adults’ blood DNA . Instead of a single-tissue DNAm aging biomarker, Horvath et al.  created the Horvath’s clock based on 353 CpGs using over 30 different tissue and cell types collected from children and adults. Nevertheless, there are two main limitations of the first generation of BoA. First, while they exhibit statistically significant associations with many age-related diseases [38, 39], the effect sizes are typically small or moderate . Second, only weak associations with clinical measures of physiological dysregulation were observed [35, 36].
In the interest of obtaining more powerful DNAm-based estimators of epigenetic BoA, Levine et al.  developed DNAm PhenoAge based on 513 CpGs, which derived from clinical biomarkers. While the first generation selects CpGs to optimize prediction of chronological age, the CpGs in DNAm PhenoAge is optimized to predict a multi-system proxy of physiological dysregulation (i.e., phenotypic aging). In doing so, the investigators were able to not only capture CpGs that exhibited strong correlations with chronological age, but those that capture variations in the risk of death and diseases among same aged individuals .
Recently, Lu et al.  employed a novel two-stage procedures and calculated DNAm GrimAge, which is a linear combination of chronological age, sex, DNAm-based surrogate biomarkers, and smoking pack-years. Compared with DNAm PhenoAge, DNAm GrimAge includes more clinical biomarkers, such as cardiovascular disease related plasma proteins (e.g. C-reactive protein, and growth differentiation factor-15). Hence, DNAm GrimAge outperforms all other DNAm-based biomarkers, on a variety of health-related metrics .
Since chronological age was highly correlated with epigenetic ages (Supplementary Figure 5 in the supplementary material), we used the corresponding age accelerations, which are new indices of accelerated biological aging and display a high risk of various health conditions (e.g., lung cancer, all-cause mortality, Parkinson’s disease, and Alzheimer’s disease) [22, 36, 39, 40] to investigate the relationship between epigenetic BoA and lung function in this study. IEAA is defined as the residual resulting from regressing the Horvath’s clock on chronological age and blood cell count estimates. By definition, IEAA does not depend on chronological age or on blood cell counts . EEAA is calculated in three steps: i) calculation of the Hannum’s clock; ii) inclusion the blood cell types to generate the age estimate (called BioAge4 in epigenetic clock software); iii) calculation the residual from a univariate model regressing BioAge4 on chronological age . EEAA captures age-related changes in blood cell types. Moreover, IEAA and EEAA, by definition, are not correlated (r<0.01) with chronological age.
GrimAgeAccel is the corresponding raw residual in the linear regression model with DNAm GrimAge regressed on chronological age (i.e. the difference between the observed value of DNAm GrimAge minus its expected value) . Similarly, PhenoAgeAccel is the residual calculated by a linear regression model in which DNAm PhenoAge is the outcome, and chronological age is the independent variable . A positive value of GrimAgeAccel/PhenoAgeAccel indicates that the biological age of the observation is higher than expected based on chronological age. GrimAgeAccel and PhenoAgeAccel, which by definition, are not correlated (r<0.01) with chronological age. To date, only limited studies have investigated the relationship between epigenetic age acceleration and pulmonary health. For instance, Levine et al.  found that IEAA was significantly associated with lung cancer incidence (hazard ratio: 1.5, P = 0.003) among 2,029 females in the Women’s Health Initiative.
Unlike the above epigenetic BoA, which are in units of years, Zhang’s DNAmRiskScore is a continuous risk score ranging from -4 to 0. Zhang et al.  developed this risk score based on ten CpG sites. The DNAmRiskScore is a predictor of all-cause mortality. Even though it includes only 10 CpGs, Zhang’s DNAmRiskScore had significantly and stable negative associations with lung function in our analyses. To understand this, we investigated the common CpG sites among different epigenetic BoA. Since DNAm GrimAge has been patented, Liu et al. could not release the 1,030 CpG identifiers. We found that DNAm PhonoAge had 41 CpGs in common with the Horvath’s clock, and 6 CpGs in common with the Hannum’s clock. The Horvath’s clock shared 6 CpG sites with the Hannum’s clock. However, Zhang’s DNAmRiskScore didn’t share any CpGs with the other three BoA [4, 18, 19, 21].
To further explore Zhang’s DNAmRiskScore, we investigated the genes where these CpG sites are located. Interestingly, one of the 10 CpGs (cg05575921) maps to the AHRR gene. Recently, studies have shown that hypo-methylation in AHRR at cg05575921 not only strongly reflects smoking history, but relates to lower lung function , Additionally, it predicts future smoking-related mortality . Another CpG site (cg19572487) maps to the retinoic acid receptor alpha gene, which plays an important role in various human cancers such as lung cancer [35, 39]. Moreover, both cg05575921 and cg19572487 are not contained in the other three epigenetic BoA in the present study (We were not able to investigate the 1,030 CpGs involved in DNAm GrimAge due to the patent issue as mentioned above).
For smoking-stratified analysis, some epigenetic BoA were associated with lower lung function in ever smokers (current and former). The subgroup analyses indicate that the significant associations between lower lung function and BoA were evident only among smokers. With aging, lung function experiences a progressive decrease due to changes in physiology and structure in the lung. In the present study, among never smokers, the lung function decline was mainly associated with chronological age, BMI, and height; whereas among smokers, the impaired pulmonary function was related to pack-years and epigenetic BoA. These findings indicate that smoking status may result in significant abnormal alteration in age-related genes. Epigenomic variation may help to explain features of lower lung function and related pathophysiology in smokers. Additionally, the findings support the importance of smoking cessation, especially for elderly people.
Several study limitations should be noted. First, our findings are based on an elderly cohort of white males. Hence, additional studies involving diverse demographic groups (younger or women) will be helpful to complement our conclusions. Second, we used an old spirometer to measure pulmonary health (it was not replaced to ensure continuity of the measurements). However, the acceptability of spirograms was judged according to American Thoracic Society/European Respiratory Society criteria  and the spirometric values were adjusted by body temperature and pressure . Third, we cannot rule out unmeasured confounding. However, we used existing literature and a priori knowledge of clinical relevance to adjust for potential confounders.
The major strengths in the study are the following: First, we investigated different types of BoA, including non-epigenetic- and epigenetic-based measures. Therefore, our work is based on multifaceted analyses of the associations between age estimators and lung function. Second, this evidence not only sheds new light into pathogenesis and susceptibility to lower lung function, but since DNAm changes are reversible, epigenetic BoA might be useful for identifying anti-aging interventions in the lung, especially for elderly smokers.
In conclusion, DNAm GrimAge, and Zhang’s DNAmRiskScore were associated with lower lung function. TL and mtDNA-CN were not related to pulmonary impairment. Epigenetic mechanisms such as DNAm may provide further explanation for decreases in lung function as individual age.
Materials and Methods
The Normative Aging Study is a closed cohort established by the U.S. Veterans Administration in 1963 . It enrolled 2,280 men volunteers living in the Boston area, aged between 21-80 years, who were free of known chronic medical conditions at enrollment. The study population had undergone clinical examinations every 3 to 5 years. These examinations took place in the morning after an overnight fast. Data on demography, health status, mediation use, and blood sample were collected at each visit.
We narrowed down our analyses to visits with DNA samples, which were collected from 1999 to 2013 . To reduce study heterogeneity, those who were non-white (3%) were excluded. Subjects with a diagnosis of leukemia (n=11) were also be removed, due to its potential influence on DNAm . This study was approved by the Institutional Review boards of the Department of Veterans Affairs and the Harvard TH Chan School of Public Health. All participants provided their written informed consent. The flowchart of our study participants is shown in Figure 3.
Pulmonary health measures
Spirometry was assessed in the standing position with a nose clip using a Collins Survey Spirometer with Eagle II Microprocessor (Warren E. Collins, Braintree, Massachusetts), and a protocol that adhere to American Thoracic Society standards . Acceptability of spirograms was judged according to American Thoracic Society/European Respiratory Society criteria . A minimum of three acceptable spirograms was obtained, of which at least two were reproducible within 5% for FVC, and FEV1. Spirometric values were adjusted by body temperature and pressure .
DNA samples were extracted from the buffy coat of whole blood collected from each visit using the QIAamp DNA Blood Kit (Qiagen, CA, USA) between 1999 and 2013. To ensure a similar age distribution and avoid confounding across plates, we randomized samples based on a two-stage age stratified algorithm . We measured DNAm using the Illumina Infinium HumanMethylation450 BeadChip (Infinium HD Methylation protocol, Illumina, San Diego, CA, USA), which includes ~ 485, 000 CpGs at a single nucleotide resolution. In the quality control step, we removed samples with a detection p-value < 0.01 based on an estimation of the background distribution using non-specific fluorescence , and corrected for dye-bias .
In the website of New Methylation Age Calculator (https://dnamage.genetics.ucla.edu/new), we uploaded our methylation data file and sample annotation file. We selected “Normalize Data”, and “Advanced Analysis” before submitting our data. We acquired four measures of BoA - epigenetic age accelerations (i.e. GrimAgeAccel, PhenoAgeAccel, EEAA, IEAA), and four DNAm ages (i.e. DNAm GrimAge, DNAm PhenoAge, the Hannum’s clock, the Horvath’s clock) via email, which was used to register in the website.
Zhang’s DNAm RiskScore was calculated based on 10 selected CpGs. The formula  is: cg01612140*(-0.38253) + cg05575921*(-0.92224) + cg06126421*(-1.70129) + cg08362785*(2.71749) + cg10321156*(-0.02073) + cg14975410*(-0.04156) + cg19572487*(-0.28069) + cg23665802*(-0.89440) + cg24704287*(-2.98637) + cg25983901*(-1.80325).
We measured leukocyte TL by quantitative real-time polymerase chain reaction . Relative leukocyte TL was calculated using a method previously described [12, 13]. Briefly, it is calculated as the ratio of the telomere (T) repeat copy number to a single-copy gene (S) copy number (T:S ratio) in a given sample, and reported as relative units expressing the ratio between TL in the test DNA and TL in a DNA pool used to generate a stand curve in each PCR run. We ran all samples in triplicates, and the average of the three T measurements was divided by the average of the three S measurements to calculate the average T:S ratio.
For mitochondrial DNA copy number (mtDNA-CN), we also adapted a real-time polymerase chain reaction and quantified it using mtDNA 12S ribosomal RNA TaqMan probe, as described in detail previously . We adapted a multiplex RT-PCR to measure mtDNA content . The mtDNA 12 S ribosomal ribonucleic acid TaqMan (Applied Biosystems, Waltham, Massachusetts) probe (6FAM-5′ TGCCAGCCACCGCG 3′-MGB) was used to measure mtDNA-CN. The sequences of primers used for amplification in mtDNA were mtF805 (5′CCACGGGAAACAGCAGTGATT3′) and mtR927 (5′CTATTGACTTGGGTTAATCGTGTGA3′). All samples were run in triplicate. The mean of the three measurements was used for statistical analyses.
Separate linear mixed-effect models were applied to estimate the associations between each BoA and lung function, with subject-specific intercept to account for repeated estimation of lung function. For outcome Y ij representing either FEV1, FEV1/FVC or MMEF for subject i on occasion j, the model is
where β0 is the intercept, μi is the random intercept for participant i; β1 is the association between a chosen BoA measure and lung function for that subject on that occasion; X2ij to Xpij represent the p-1 covariates, the selection of which differed slightly depending on the measure of BoA used in the analysis (see below). Effect estimates are expressed per one year increase in GrimAgeAccel / PhenoAgeAccel / IEAA / EEAA, or 0.2 increase in Zhang DNAmRiskScore / mtDNA-CN, or 1 increase in TL.
We selected the fundamental covariates a priori: chronological age, body mass index (BMI), height, smoking status [never smokers who never smoked before versus former or current smokers)], cigarette pack-years, years of education, corticosteroid use, estimated cell types (CD4+ T lymphocytes, CD8+ T lymphocytes, natural killer cells, B cells, and monocytes) derived by the method of Houseman et al. , batch effects for all BoA except for IEAA which measures “pure” epigenetic aging effects that are not confounded by differences in blood cell counts, and EEAA which is accounted for cell types in the calculation steps .
For Zhang’s DNAmRiskScore, which was calculated from the formula (1), we additionally adjusted for important technical covariates because normalization on the methylation data might diminish biologic signals. Therefore, we included important technical covariates from the control metrics monitoring the performance of various experimental steps . The potential technical covariates included time when DNA methylation measures were performed and experimental parameters including Staining Green, Staining Red, Extension Green, Extension Red, Hybridization High Medium, Hybridization Medium Low, Target Removal 1, Target Removal 2, Bisulfite Conversion I Green, Bisulfite Conversion I Red, Bisulfite Conversion II, Specificity I Green, Specificity I Red, Specificity II, Non Polymorphic Green, and Non Polymorphic Red. We used elastic-net regularized generalized linear model and selected five important technical covariates: Non polymorphic Red, Specificity I Red, Bisulfite Conversion I Red, Bisulfite Conversion II, Extension Red.
Statistical analysis were performed using R (3.5.1) with “lme4” package (linear mixed-effect model) and “glmnet” package (lasso and elastic-net regularized generalized linear model). To adjust for multiple comparisons, we applied Benjamin-Hochberg false discovery rate (FDRB-H) and set the false positive threshold as 0.05 .
We repeated the models in several different ways to conduct sensitivity analyses. First, we excluded visits at which corticosteroid use was reported (N=82); second, we adjusted models for chronic diseases (coronary heart diseases, stroke, and diabetes) when their DNA samples were collected; third, to adjust for the possibility that healthier men are more likely to return for subsequent exams, we applied inverse probability weighting  using logistic regression to calculate the probability of having a subsequent visit given chronological age, education, BMI, blood pressure, smoking status, cigarette pack years, alcohol consumption, C-reactive protein, asthma, chronic bronchitis, and emphysema at previous visit. Finally, we divided the subjects into two groups based on their smoking status, and tested the BoA-lung function associations separately. One groups were subjects who never smoked before (never smokers), and the other one were subjects who were former or current smokers.
SD: standard deviation; BMI: body mass index; FEV1:: forced expiratory volume in 1 second; FVC: forced vital capacity; MMEF: maximum mid-expiratory flow; TL: Telomere length; mtDNA-CN: mitochondrial DNA copy number; FDR B-H: Benjamin-Hochberg false discovery rate.
C.W. takes full responsibility for the content of the manuscript, including data and analysis. C.W., A.B., J.S. conceived of the idea; P.S.V., A.B., J.S., D.S. contributed data; C.W., A.J., J.H., L.H., Y.Z. interpreted data; C.W., J.S drafted the manuscript; J.H., L.H., Y.Z., B. C., D.S., A.B., J.S revised the manuscript for important intellectual content; A.B., J.S. obtained funding.
We thank all those who participated in the Normative Aging Study, including the subjects as well as the technicians for lung function testing, blood sampling, methylation measures, sample storage, data management, and all those who supervised the study.
Conflicts of Interest
No financial interests or potential conflicts of interest.
This work was supported by the National Institute of Environmental Health Sciences (grants ES000002, R01ES015172, 5R01ES015172-09, 5R01ES027747-02, P30ES009089, R01ES021733, R01ES025225). The VA Normative Aging Study is supported by the Cooperative Studies Program/Epidemiology Research and Information Center of the U.S. Department of Veterans Affairs and is a component of the Massachusetts Veterans Epidemiology Research and Information Center in Boston, MA.
- 1. Ortman JM, Velkoff VA, Hogan H. An Aging Nation: The Older Population in the United States Washington D.C.: U.S. Department of Commerce May 2014. P25–1140.
- 2. Jaul E, Barron J. Age-related diseases and clinical and public health implications for the 85 years old and over population. Front Public Health. 2017; 5:335. https://doi.org/10.3389/fpubh.2017.00335 [PubMed]
- 3. Lowsky DJ, Olshansky SJ, Bhattacharya J, Goldman DP. Heterogeneity in healthy aging. J Gerontol A Biol Sci Med Sci. 2014; 69:640–49. https://doi.org/10.1093/gerona/glt162 [PubMed]
- 4. Levine ME, Lu AT, Quach A, Chen BH, Assimes TL, Bandinelli S, Hou L, Baccarelli AA, Stewart JD, Li Y, Whitsel EA, Wilson JG, Reiner AP, et al. An epigenetic biomarker of aging for lifespan and healthspan. Aging (Albany NY). 2018; 10:573–91. https://doi.org/10.18632/aging.101414 [PubMed]
- 5. Anstey KJ, Windsor TD, Jorm AF, Christensen H, Rodgers B. Association of pulmonary function with cognitive performance in early, middle and late adulthood. Gerontology. 2004; 50:230–34. https://doi.org/10.1159/000078352 [PubMed]
- 6. Singh-Manoux A, Dugravot A, Kauffmann F, Elbaz A, Ankri J, Nabi H, Kivimaki M, Sabia S. Association of lung function with physical, mental and cognitive function in early old age. Age (Dordr). 2011; 33:385–92. https://doi.org/10.1007/s11357-010-9189-x [PubMed]
- 7. Liao D, Higgins M, Bryan NR, Eigenbrodt ML, Chambless LE, Lamar V, Burke GL, Heiss G. Lower pulmonary function and cerebral subclinical abnormalities detected by MRI: the atherosclerosis risk in communities study. Chest. 1999; 116:150–56. https://doi.org/10.1378/chest.116.1.150 [PubMed]
- 8. Drummond MB, Hansel NN, Connett JE, Scanlon PD, Tashkin DP, Wise RA. Spirometric predictors of lung function decline and mortality in early chronic obstructive pulmonary disease. Am J Respir Crit Care Med. 2012; 185:1301–06. https://doi.org/10.1164/rccm.201202-0223OC [PubMed]
- 9. Lee HM, Le H, Lee BT, Lopez VA, Wong ND. Forced vital capacity paired with framingham risk score for prediction of all-cause mortality. Eur Respir J. 2010; 36:1002–06. https://doi.org/10.1183/09031936.00042410 [PubMed]
- 10. Mannino DM, Buist AS, Petty TL, Enright PL, Redd SC. Lung function and mortality in the united states: data from the first national health and nutrition examination survey follow up study. Thorax. 2003; 58:388–93. https://doi.org/10.1136/thorax.58.5.388 [PubMed]
- 11. Neas LM, Schwartz J. Pulmonary function levels as predictors of mortality in a national sample of US adults. Am J Epidemiol. 1998; 147:1011–18. https://doi.org/10.1093/oxfordjournals.aje.a009394 [PubMed]
- 12. McCracken J, Baccarelli A, Hoxha M, Dioni L, Melly S, Coull B, Suh H, Vokonas P, Schwartz J. Annual ambient black carbon associated with shorter telomeres in elderly men: veterans affairs normative aging study. Environ Health Perspect. 2010; 118:1564–70. https://doi.org/10.1289/ehp.0901831 [PubMed]
- 13. Colicino E, Wilson A, Frisardi MC, Prada D, Power MC, Hoxha M, Dioni L, Spiro A, Vokonas PS, Weisskopf MG, Schwartz JD, Baccarelli AA. Telomere length, long-term black carbon exposure, and cognitive function in a cohort of older men: the VA normative aging study. Environ Health Perspect. 2017; 125:76–81. https://doi.org/10.1289/EHP241 [PubMed]
- 14. Blasco MA. Telomeres and human disease: ageing, cancer and beyond. Nat Rev Genet. 2005; 6:611–22. https://doi.org/10.1038/nrg1656 [PubMed]
- 15. Morten KJ, Ashley N, Wijburg F, Hadzic N, Parr J, Jayawant S, Adams S, Bindoff L, Bakker HD, Mieli-Vergani G, Zeviani M, Poulton J. Liver mtDNA content increases during development: a comparison of methods and the importance of age- and tissue-specific controls for the diagnosis of mtDNA depletion. Mitochondrion. 2007; 7:386–95. https://doi.org/10.1016/j.mito.2007.09.001 [PubMed]
- 16. Alvarez-Mora MI, Podlesniy P, Gelpi E, Hukema R, Madrigal I, Pagonabarraga J, Trullas R, Mila M, Rodriguez-Revenga L. Fragile x-associated tremor/ataxia syndrome: regional decrease of mitochondrial DNA copy number relates to clinical manifestations. Genes Brain Behav. 2019; 18:e12565. https://doi.org/10.1111/gbb.12565 [PubMed]
- 17. Oberdoerffer P, Sinclair DA. The role of nuclear architecture in genomic instability and ageing. Nat Rev Mol Cell Biol. 2007; 8:692–702. https://doi.org/10.1038/nrm2238 [PubMed]
- 18. Hannum G, Guinney J, Zhao L, Zhang L, Hughes G, Sadda S, Klotzle B, Bibikova M, Fan JB, Gao Y, Deconde R, Chen M, Rajapakse I, et al. Genome-wide methylation profiles reveal quantitative views of human aging rates. Mol Cell. 2013; 49:359–67. https://doi.org/10.1016/j.molcel.2012.10.016 [PubMed]
- 19. Horvath S. DNA methylation age of human tissues and cell types. Genome Biol. 2013; 14:R115. https://doi.org/10.1186/gb-2013-14-10-r115 [PubMed]
- 20. Lu AT, Quach A, Wilson JG, Reiner AP, Aviv A, Raj K, Hou L, Baccarelli AA, Li Y, Stewart JD, Whitsel EA, Assimes TL, Ferrucci L, Horvath S. DNA methylation GrimAge strongly predicts lifespan and healthspan. Aging (Albany NY). 2019; 11:303–27. https://doi.org/10.18632/aging.101684 [PubMed]
- 21. Zhang Y, Wilson R, Heiss J, Breitling LP, Saum KU, Schöttker B, Holleczek B, Waldenberger M, Peters A, Brenner H. DNA methylation signatures in peripheral blood strongly predict all-cause mortality. Nat Commun. 2017; 8:14617. https://doi.org/10.1038/ncomms14617 [PubMed]
- 22. Horvath S, Ritz BR. Increased epigenetic age and granulocyte counts in the blood of parkinson’s disease patients. Aging (Albany NY). 2015; 7:1130–42. https://doi.org/10.18632/aging.100859 [PubMed]
- 23. Elina Sillanpää, Sarianna Sipilä, Timo Törmäkangas, Jaakko Kaprio, Taina Rantanen. Genetic and environmental effects on telomere length and lung function: a twin study. The journals of gerontology. Series A., Biological sciences and medical sciences. 2017; 72:1561–1568. https://doi.org/10.1093/gerona/glw178
- 24. Minh Thien Nguyen, Richard Saffery, David Burgner, Kate Lycett, Regan Vryer, Anneke Grobler, Terence Dwyer, Sarath Ranganathan, Melissa Wake. Telomere length and lung function in a population-based cohort of children and mid-life adults. Pediatric pulmonology. 2019; 54:2044–2052. https://doi.org/10.1002/ppul.24489 [PubMed]
- 25. Liu SF, Kuo HC, Tseng CW, Huang HT, Chen YC, Tseng CC, Lin MC. Leukocyte mitochondrial DNA copy number is associated with chronic obstructive pulmonary disease. PLoS One. 2015; 10:e0138716. https://doi.org/10.1371/journal.pone.0138716 [PubMed]
- 26. Lepeule J, Baccarelli A, Motta V, Cantone L, Litonjua AA, Sparrow D, Vokonas PS, Schwartz J. Gene promoter methylation is associated with lung function in the elderly: the normative aging study. Epigenetics. 2012; 7:261–69. https://doi.org/10.4161/epi.7.3.19216 [PubMed]
- 27. Bolund AC, Starnawska A, Miller MR, Schlünssen V, Backer V, Børglum AD, Christensen K, Tan Q, Christiansen L, Sigsgaard T. Lung function discordance in monozygotic twins and associated differences in blood DNA methylation. Clin Epigenetics. 2017; 9:132. https://doi.org/10.1186/s13148-017-0427-2 [PubMed]
- 28. Carmona JJ, Barfield RT, Panni T, Nwanaji-Enwerem JC, Just AC, Hutchinson JN, Colicino E, Karrasch S, Wahl S, Kunze S, Jafari N, Zheng Y, Hou L, et al. Metastable DNA methylation sites associated with longitudinal lung function decline and aging in humans: an epigenome-wide study in the NAS and KORA cohorts. Epigenetics. 2018; 13:1039–55. https://doi.org/10.1080/15592294.2018.1529849 [PubMed]
- 29. Hung M, Bounsanga J, Voss MW. Interpretation of correlations in clinical research. Postgrad Med. 2017; 129:902–06. https://doi.org/10.1080/00325481.2017.1383820 [PubMed]
- 30. Lung Capacity and Aging. American Lung Association. https://www.lung.org/lung-health-and-diseases/how-lungs-work/lung-capacity-and-aging.html. Updated March 11, 2020. Accessed April 28, 2020.
- 31. García-Río F, Pino JM, Dorgham A, Alonso A, Villamor J. Spirometric reference equations for european females and males aged 65-85 yrs. Eur Respir J. 2004; 24:397–405. https://doi.org/10.1183/09031936.04.00088403 [PubMed]
- 32. Greider CW. Telomere length regulation. Annu Rev Biochem. 1996; 65:337–65. https://doi.org/10.1146/annurev.bi.65.070196.002005 [PubMed]
- 33. Sanders JL, Newman AB. Telomere length in epidemiology: a biomarker of aging, age-related disease, both, or neither? Epidemiol Rev. 2013; 35:112–31. https://doi.org/10.1093/epirev/mxs008 [PubMed]
- 34. Quach A, Levine ME, Tanaka T, Lu AT, Chen BH, Ferrucci L, Ritz B, Bandinelli S, Neuhouser ML, Beasley JM, Snetselaar L, Wallace RB, Tsao PS, et al. Epigenetic clock analysis of diet, exercise, education, and lifestyle factors. Aging (Albany NY). 2017; 9:419–46. https://doi.org/10.18632/aging.101168 [PubMed]
- 35. Horvath S, Raj K. DNA methylation-based biomarkers and the epigenetic clock theory of ageing. Nat Rev Genet. 2018; 19:371–84. https://doi.org/10.1038/s41576-018-0004-3 [PubMed]
- 36. Levine ME, Hosgood HD, Chen B, Absher D, Assimes T, Horvath S. DNA methylation age of blood predicts future onset of lung cancer in the women’s health initiative. Aging (Albany NY). 2015; 7:690–700. https://doi.org/10.18632/aging.100809 [PubMed]
- 37. Fernández-Vizarra E, Enríquez JA, Pérez-Martos A, Montoya J, Fernández-Silva P. Tissue-specific differences in mitochondrial activity and biogenesis. Mitochondrion. 2011; 11:207–13. https://doi.org/10.1016/j.mito.2010.09.011 [PubMed]
- 38. Marioni RE, Shah S, McRae AF, Ritchie SJ, Muniz-Terrera G, Harris SE, Gibson J, Redmond P, Cox SR, Pattie A, Corley J, Taylor A, Murphy L, et al. The epigenetic clock is correlated with physical and cognitive fitness in the lothian birth cohort 1936. Int J Epidemiol. 2015; 44:1388–96. https://doi.org/10.1093/ije/dyu277 [PubMed]
- 39. Levine ME, Lu AT, Bennett DA, Horvath S. Epigenetic age of the pre-frontal cortex is associated with neuritic plaques, amyloid load, and alzheimer’s disease related cognitive functioning. Aging (Albany NY). 2015; 7:1198–211. https://doi.org/10.18632/aging.100864 [PubMed]
- 40. Perna L, Zhang Y, Mons U, Holleczek B, Saum KU, Brenner H. Epigenetic age acceleration predicts cancer, cardiovascular, and all-cause mortality in a german case cohort. Clin Epigenetics. 2016; 8:64. https://doi.org/10.1186/s13148-016-0228-z [PubMed]
- 41. Kodal JB, Kobylecki CJ, Vedel-Krogh S, Nordestgaard BG, Bojesen SE. AHRR hypomethylation, lung function, lung function decline and respiratory symptoms. Eur Respir J. 2018; 51:1701512. https://doi.org/10.1183/13993003.01512-2017 [PubMed]
- 42. Bojesen SE, Timpson N, Relton C, Davey Smith G, Nordestgaard BG. AHRR (cg05575921) hypomethy-lation marks smoking behaviour, morbidity and mortality. Thorax. 2017; 72:646–53. https://doi.org/10.1136/thoraxjnl-2016-208789 [PubMed]
- 43. Miller MR, Hankinson J, Brusasco V, Burgos F, Casaburi R, Coates A, Crapo R, Enright P, van der Grinten CP, Gustafsson P, Jensen R, Johnson DC, MacIntyre N, et al, and ATS/ERS Task Force. Standardisation of spirometry. Eur Respir J. 2005; 26:319–38. https://doi.org/10.1183/09031936.05.00034805 [PubMed]
- 44. Standardization of spirometry—1987 update. Statement of the american thoracic society. Am Rev Respir Dis. 1987; 136:1285–98. https://doi.org/10.1164/ajrccm/136.5.1285 [PubMed]
- 45. Bell B, Rose CL, Damon A. The veterans administration longitudinal study of healthy aging. Gerontologist. 1966; 6:179–84. https://doi.org/10.1093/geront/6.4.179 [PubMed]
- 46. Wang C, Koutrakis P, Gao X, Baccarelli A, Schwartz J. Associations of annual ambient PM 2.5 components with DNAm PhenoAge acceleration in elderly men: the normative aging study. Environ Pollut. 2020; 258:113690. https://doi.org/10.1016/j.envpol.2019.113690 [PubMed]
- 47. Dai L, Mehta A, Mordukhovich I, Just AC, Shen J, Hou L, Koutrakis P, Sparrow D, Vokonas PS, Baccarelli AA, Schwartz JD. Differential DNA methylation and PM 2.5 species in a 450K epigenome-wide association study. Epigenetics. 2017; 12:139–48. https://doi.org/10.1080/15592294.2016.1271853 [PubMed]
- 48. Heiss JA, Just AC. Improved filtering of DNA methylation microarray data by detection P values and its impact on downstream analyses. Clin Epigenetics. 2019; 11:15. https://doi.org/10.1186/s13148-019-0615-3 [PubMed]
- 49. Xu Z, Langie SA, De Boever P, Taylor JA, Niu L. RELIC: a novel dye-bias correction method for illumina methylation BeadChip. BMC Genomics. 2017; 18:4. https://doi.org/10.1186/s12864-016-3426-3 [PubMed]
- 50. Cawthon RM. Telomere measurement by quantitative PCR. Nucleic Acids Res. 2002; 30:e47. https://doi.org/10.1093/nar/30.10.e47 [PubMed]
- 51. Zhong J, Cayir A, Trevisi L, Sanchez-Guerra M, Lin X, Peng C, Bind MA, Prada D, Laue H, Brennan KJ, Dereix A, Sparrow D, Vokonas P, et al. Traffic-related air pollution, blood pressure, and adaptive response of mitochondrial abundance. Circulation. 2016; 133:378–87. https://doi.org/10.1161/CIRCULATIONAHA.115.018802 [PubMed]
- 52. Andreu AL, Martinez R, Marti R, García-Arumí E. Quantification of mitochondrial DNA copy number: pre-analytical factors. Mitochondrion. 2009; 9:242–46. https://doi.org/10.1016/j.mito.2009.02.006 [PubMed]
- 53. Houseman EA, Accomando WP, Koestler DC, Christensen BC, Marsit CJ, Nelson HH, Wiencke JK, Kelsey KT. DNA methylation arrays as surrogate measures of cell mixture distribution. BMC Bioinformatics. 2012; 13:86. https://doi.org/10.1186/1471-2105-13-86 [PubMed]
- 54. BeadArray Controls Reporter. San Diego, CA: Illumina; October 2015. 1000000004009 v00.
- 55. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal statistical society: series B (Methodological). 1995; 57:289–300. https://doi.org/10.1111/j.2517-6161.1995.tb02031.x