Surrounding Greenness and Biological Aging Based on DNA Methylation: A Twin and Family Study in Australia

Abstract

Background: High surrounding greenness has many health benefits and might contribute to slower biological aging. However, very few studies have evaluated this from the perspective of epigenetics.

Objectives: We aimed to evaluate the association between surrounding greenness and biological aging based on DNA methylation.

Methods: We derived Horvath’s DNA methylation age (DNAmAge), Hannum’s DNAmAge, PhenoAge, and GrimAge based on DNA methylation measured in peripheral blood samples from 479 Australian women in 130 families. Measures of DNAmAge acceleration (DNAmAgeAC) were derived from the residuals after regressing each DNAmAge metric on chronological age. Greenness was represented by satellite-derived Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI) metrics within 300-, 500-, 1,000-, and 2,000-m buffers surrounding participant addresses. Greenness-DNAmAgeAC associations were estimated using a within-sibship design fitted by linear mixed effect models, adjusting for familial clustering and important covariates.

Results: Greenness metrics were associated with significantly lower DNAmAgeAC based on GrimAge acceleration, suggesting slower biological aging with higher greenness based on both NDVI and EVI in 300–2,000m buffer areas. For example, each interquartile range increase in NDVI within 1,000m was associated with a 0.59 (95% CI: 0.18, 1.01)–year decrease in GrimAge acceleration. Greenness was also inversely associated with three of the eight components of GrimAge, specifically, DNA methylation-based surrogates of serum cystatin-C, serum growth differentiation factor 15, and smoking pack years. Associations between greenness and biological aging measured by Horvath’s and Hannum’s DNAmAgeAC were less consistent, and depended on neighborhood socioeconomic status. No significant associations were estimated for PhenoAge acceleration.

Discussion: Higher surrounding greenness was associated with slower biological aging, as indicated by GrimAge age acceleration, in Australian women. Associations were also evident for three individual components of GrimAge, but were inconsistent for other measures of biological aging. Additional studies are needed to confirm our results.

Introduction

As a result of decreased fertility and improved survival, the population is aging rapidly across the world (United Nations 2017). The global population aged ≥60y has increased from 382 million in 1980 to 962 million in 2017 and is predicted to double again by 2050, reaching nearly 2.1 billion (United Nations 2017). With this rapid demographic shift, it is important to better understand the determinants of the biological aging process and how to achieve healthy aging (WHO 2015). Several biomarkers of biological aging have been proposed to predict the functional capability of a person, and the DNA methylation age (DNAmAge) has been suggested as the most robust (Ferrucci et al. 2020; Jylhävä et al. 2017).

DNAmAge is a measurement of biological aging based on DNA methylation, a biological process that regulates gene expression by adding methyl groups to the cytosine bases of the DNA (Horvath and Raj 2018). Four main algorithms are used to calculate DNAmAge: Horvath’s age (Horvath 2013), Hannum’s age (Hannum et al. 2013), PhenoAge (Levine et al. 2018), and GrimAge (Lu et al. 2019). The DNAmAge acceleration (DNAmAgeAC) can be estimated as the residuals after regressing individual measures of DNAmAge on chronological age (the calendar time that has passed since birth). Increased DNAmAgeAC has been positively associated with mortality and conditions related to aging, such as cancer and cardiovascular diseases (Horvath and Raj 2018; Levine et al. 2018; Lu et al. 2019). DNAmAgeAC calculated from the most recently developed estimator (GrimAge) showed the strongest association (GrimAge>PhenoAge>Hannum>Horvath ) with deaths, cancers, coronary heart disease, and early age at menopause in an analysis of longitudinal data from 6,395 adults (mean age: 63.0 y) in four cohorts (Framingham Heart Study, Women’s Health Initiative, Jackson Heart Study, and InCHIANTI Study) followed for 13.7 y on average (Lu et al. 2019).

An important feature of DNAmAge is that DNA methylation is modifiable by interventions and environmental factors, making it a useful tool for identifying or validating healthy aging determinants (Horvath and Raj 2018). One of the potentially important determinants is surrounding greenness, also called residential green space or surrounding vegetation (Reid et al. 2018), referring to a measurement of the density of vegetated land (e.g., gardens, parks, grassland) in surrounding areas (Taylor and Hochuli 2017). Greenness is a critical component of a healthy built environment. A high level of surrounding greenness may decrease mental stress, facilitate social interaction, boost physical activity, and reduce harm from noise, air pollution, and heat exposure (James et al. 2015; Markevych et al. 2017). Numerous epidemiological studies have reported evidence supporting health benefits of high surrounding greenness, including slower cognitive decline and reduced risks of mental illness, obesity, cardiovascular diseases, type 2 diabetes, and death (de Keijzer et al. 2018, 2020; Fong et al. 2018; James et al. 2015; Jia et al. 2018; Lachowycz and Jones 2011; Rojas-Rueda et al. 2019; Sarkar 2017).

To the best of our knowledge, only one study of infants, which was reported as a conference abstract, has evaluated the relationship between greenness and DNA methylation-based measures of biological aging (Sbihi et al. 2018). The study estimated gestational age (GA) based on DNA methylation measured in cord blood samples and reported a borderline association between greenness surrounding the maternal residence and DNA methylation-based GA acceleration (i.e., older DNA methylation GA compared with chronological GA) {0.6 d, 95% [confidence interval (CI)]: −0.6 to 1.8, for one interquartile range (IQR) exposure increase}. This suggested that high greenness may accelerate development in utero (Sbihi et al. 2018).

Given the evidence supporting the benefits of high surrounding greenness on health, we aimed to evaluate whether high surrounding greenness can slow biological aging in adults. This can be indicated by inverse associations between greenness exposures and DNA methylation-based measures of biological aging.

Methods

Study Design and Population

We analyzed data from the Australian Mammographic Density Twins and Sisters Study (AMDTSS), a twin and family study population (Stone et al. 2007). The AMDTSS was originally designed to investigate genetic, environmental, and lifestyle factors associated with mammographic density, a major risk factor for breast cancer (Boyd et al. 2002). In 1995–1999, female twins in three Australian cities (Melbourne, Sydney, and Perth) were recruited through the Australian Twin Registry. In 2004–2009, the twins were invited to participate further, and their sisters were invited to participate as well. The study included participants 40–70 years of age and excluded pregnant women or nursing mothers and women with a history of breast cancer. Written informed consent was obtained from all participants. Participants were asked to complete a telephone-administered questionnaire and, if they consented to provide a blood sample, were sent a blood collection kit with a shipping container, tubes, and detailed instructions and asked to visit a nearby contracted pathology laboratory where a 27-mL peripheral whole-blood sample was collected. In the case a participant was unable to travel to a pathology laboratory, a trained phlebotomist visited their home to perform the blood collection. All blood samples were delivered to the Genetic Epidemiology Laboratory at The University of Melbourne within 48 h after collection. These samples were then processed to generate Guthrie cards, which were used to measure DNA methylation.

The AMDTSS recruited over 2,700 participants during 2004–2009. All AMDTSS participants with DNA methylation measurements were eligible for the present study, which comprised 479 women from 130 families, including 66 monozygotic twin pairs, 66 dizygotic twin pairs, and 215 sisters of these twins. Each family had 3–6 participants. Each eligible participant had one DNA methylation measurement. The AMDTSS was approved by the Human Research Ethics Committee of The University of Melbourne. The present analysis was approved by the Monash University Human Research Ethics Committee.

DNA Methylation Data

DNA methylation was measured in peripheral blood samples using the Illumina Infinium HumanMethylation450 BeadChip array, as detailed in our previous study (Li et al. 2015). We processed the raw methylation data with the Bioconductor minfi package (Aryee et al. 2014). This process included normalization of data using Illumina’s reference factor-based normalization methods (preprocessIllumina) and correction of type I and II probe bias using subset-quantile within array normalization (SWAN) (Maksimovic et al. 2012). We used an empirical Bayes batch-effects removal method, ComBat (Johnson et al. 2007), to minimize the technical variation across batches. After excluding 65 probes corresponding to known single nucleotide polymorphisms, and probes with detection p>0.01 in one or more samples, 479,957 probes were left for further analyses (Li et al. 2015). The methylation level of each cytosine–phosphate–guanine site (CpG) was represented by a probe’s beta value ranging from 0 (not methylated) to 1 (fully methylated).

The processed DNA methylation data were used to calculate DNAmAge and cell-type proportions (i.e., proportions of leukocyte cells) in peripheral blood. Four measures of biological age based on methylation of specific CpGs were derived using an online calculator ( http://dnamage.genetics.ucla.edu/new): Horvath’s age, based on 353 CpGs (Horvath 2013); Hannum’s age, based on 71 CpGs (Hannum et al. 2013); Pheno age, based on 513 CpGs (Levine et al. 2018); and Grim age, based on 1,030 CpGs (Lu et al. 2019). In our processed DNA methylation data set, methylation data were missing for 1 CpGs used to calculate Horvath’s age, 2 CpGs used to calculate PhenoAge, and 21 CpGs used to calculate GrimAge, whereas methylation data were complete for all CpGs used to derive Hannum’s age. The methylation levels for the missing CpGs were imputed by the online calculator using the default k nearest neighbor (KNN) method through the impute.knn function, with default parameters embedded in the impute R package (Di Lena et al. 2019; Hastie et al. 2018; Horvath 2013; Troyanskaya et al. 2001). In a simulation study using data from 1,495 healthy control samples in the Gene Expression Omnibus, values of Horvath’s DNAmAge derived using KNN to impute methylation data that were missing completely at random for 10% of the CpGs used to derive the metric were consistent with Horvath’s DNAmAge estimates based on the complete data (Pearson correlation coefficient=0.99 , mean absolute error=1.8y ) (Di Lena et al. 2019). Proportions of CD8+ T cells , CD4+ T cells , natural killer (NK) cells, B cells, monocytes, and granulocytes were estimated using a widely used published algorithm (Houseman et al. 2012). The proportions of exhausted CD8 T cells, naïve CD8+ T cells , naïve CD4+ T cells , and plasma cells (effector B cells) were estimated using the algorithm embedded in the online DNAmAge calculator (Horvath 2013). The online calculator also provides separate estimates for the eight individual components that contribute to GrimAge, including smoking pack years and the concentration of seven plasma proteins [adrenomedullin, beta-2-microglobulin, cystatin-C, growth differentiation factor 15 (GDF-15), leptin, plasminogen activator inhibitor 1 (PAI-1), and tissue inhibitor metalloproteinases 1 (TIMP-1)], which were all estimated by DNA methylation data based on the model trained in Framingham Heart Study samples (Lu et al. 2019).

Exposure Assessment

Participants who agreed to donate a blood sample provided residential address information (state, city, suburb, street name, and house number) for shipment of the blood collection kit. We transformed participant addresses into longitudes and latitudes using the Google Map Application Programming Interface (API) through the geocode function embedded in ggmap R package (version 3.0.0; R Development Core Team). These longitudes and latitudes were then used for exposure assessment. According to the 2011 Australian Statistical Geography Standard, 84% (i.e., 403/479) of the residential addresses were in urban areas (ABS 2011, 2020). Because residential histories were not collected, we assumed that the participants lived at the same residence during the 12 months before their blood draw. This assumption is supported by a national survey conducted in 2007–2008 that reported that nearly 80% of Australians 40–79 years of age had not changed residences in the last 5 y (ABS 2010). At the time of blood draw, 12 (2.5%) participants (two monozygotic twin pairs, two dizygotic twin pairs, and two non-twin sister pairs) lived at the same residential address with her twin or sister.

Residential surrounding greenness was measured using the Normalized Difference Vegetation Index (NDVI) and the Enhanced Vegetation Index (EVI), which were derived from Moderate resolution Imaging Spectroradiometer (MODIS) images collected from the National Aeronautics and Space Administration’s Terra satellite (Didan 2015). Green vegetation strongly absorbs visible red light (wavelengths from 0.6 to 0.7μm ) for photosynthesis, whereas it reflects most near-infrared light (wavelengths from 0.7 to 1.1μm ). According to this characteristic, NDVI is calculated as NDVI=(NIR− Red)/(NIR + Red) , where Red and NIR represent land surface reflectance of red light and near-infrared light measured by satellite sensors, respectively. The EVI is calculated similarly to the NDVI, but it corrects for distortions in the reflected light due to atmospheric aerosols and the land surface background beneath the vegetation. The EVI has a higher sensitivity to the vegetation difference in high biomass regions (e.g., forests) than the NDVI. The MODIS provides gridded NDVI and EVI values at a 250×250m spatial resolution for every 16 d since 2000. The 16-d composites were created on a per-pixel basis, using the Constrained Maximum Value Composite algorithm. Briefly, the algorithm selects the best observations from all images within 16 d, minimizing the effects of image quality, cloud, and view angle. Both NDVI and EVI values range from −1 to 1, with negative values representing ice, water, and nonvegetated soil and higher positive values corresponding to higher greenness.

Based on the 16-d gridded NDVI/EVI data set, we calculated the residential surrounding greenness in several steps. First, we selected all 16-d NDVI/EVI gridded observations at a 250m resolution during 2003–2009 across Australia. Second, we converted negative NDVI/EVI values in all pixels to zero so that negative values would not offset positive values in calculating the areal average (Reid et al. 2018). Third, we calculated the surrounding NDVIs/EVIs within different radial buffers (300, 500, 1,000, and 2,000m ) centered on each participant’s home address for every 16 d, by averaging the values of all pixels weighted by area size covered by the buffer areas. We chose 300-, 500-, 1,000-, and 2,000-m buffers, because the health impacts of surrounding greenness are most likely to be evident within these buffers (Browning and Lee 2017; Reid et al. 2018). Finally, to reduce the influence of seasonal variations across different locations, we classified exposures for each index and radial buffer using the highest 16-d NDVI (or EVI) during the year (365 d) before each participant’s blood draw.

Covariate Data

Through a telephone-administered questionnaire, we collected data on birth date, educational level, marital status, smoking behavior, and alcohol use behavior (see the section “Variables measured by telephone administrated questionnaire” in the Supplemental Material for details). For 95% of the participants (457/479), the telephone interview was conducted within 1 y before the blood draw, and the other participants’ telephone interviews were conducted within 1–2 y before the blood draw or up to 3 months after the blood draw (Figure S1). Chronological age was defined as the calendar time period between birth date and the day when the blood sample was collected. Area-level socioeconomic status (SES) was represented by the percentile of Index of Relative Socioeconomic Advantage and Disadvantage (IRSAD) in 2006 for the population census collection district (CD) in which the participant lived. The IRSAD is a linear combination of multiple indicators of community-level SES (e.g., household income, educational level, occupation, housing conditions) collected from the population census in Australia (ABS 2006). A higher IRSAD represents a higher area-level SES. The CD was the smallest geographic area for which the IRSAD data were available. In our study, the participants were distributed across 366 CDs, with the usual population size ranging from 82 to 1,686.

We collected daily ambient temperature and relative humidity from a gridded climate data set at 0.05°×0.05° (about 5×5km ) spatial resolution (Jeffrey et al. 2001). Each participant was assigned the values of the pixel where her home address was located. For each participant, we used the average ambient temperature and relative humidity within 1 y (365 d) prior to the day blood was taken in order to represent the participant’s exposure to meteorological factors.

Statistical Analysis

Main model analyses.

We used the same procedure to estimate associations between each DNAmAgeAC metric and surrounding greenness metric. First, we calculated the DNAmAgeAC as residuals after regressing DNAmAge on chronological age (Horvath and Raj 2018; Levine et al. 2018; Lu et al. 2019). Then, we estimated associations between greenness and DNAmAgeAC using a within-sibship design (Li et al. 2017; Stone et al. 2010) fitted by a linear mixed-effect model, with random family-specific intercepts to capture familial clustering:

E(DNAmAgeACij)=β0+βwithin×(Greenij−¯Greeni)+βbetweenׯGreeni+β×Xij+μi, [1]

where DNAmAgeACij is the DNAmAge acceleration of subject j from family i ; Greenij j is the surrounding greenness index for subject j from family i (e.g., NDVI within 1,000m ); ¯Greeni is the average surrounding greenness of all subjects in family i , with each subject’s surrounding greenness defined according to their residential address and blood draw date; and μi is the random intercept for family i . We modeled greenness as a linear term because the linear model had a lower Bayesian information criterion (BIC) value than models of greenness as a natural cubic spline with 2–4 degrees of freedom (df) (Table S1).

βbetween and βwithin represent effect estimates (difference in mean DNAmAgeAC per IQR increase in greenness) for between- and within-sibship associations, respectively. The between-sibship association represents a cross-sectional comparison between families that may be biased by familial factors and genetic background that vary between families. By contrast, the within-sibship association represents intra-family differences between twins and sisters, thus reducing or eliminating confounding by genetic background and shared familial factors (Li et al. 2017; Stone et al. 2010). Therefore, we report βwithin and its 95% CI throughout this article.

Xij represents a set of covariates that might bias greenness–DNAmAgeAC associations (James et al. 2015; Markevych et al. 2017), including educational level (secondary or below; vocational training; university), marital status (married or de facto; never married; widowed or separated or divorced), smoking behavior (current smoker, former smoker, never smoked), alcohol use (current drinker, former drinker, never drank), and area-level SES (four categories according to quartiles in Australia); plus continuous variables for chronological age, annual mean temperature, annual mean relative humidity, and the proportions of seven types of leukocyte cells (naïve CD8+ T cells , CD8+ T cells , plasma cells, CD4+ T cells , NK cells, monocytes, and granulocytes), as suggested by (Lu et al. 2019). Given adjustment for cell-type proportions, effect estimates from our main model should be interpreted as associations between greenness and differences in intrinsic DNAmAgeAC (independent of differences in leukocyte proportions) (Horvath and Raj 2018).

Missing values for educational level (two participants) and area-level SES (one participant), were addressed by five multiple imputations using the default method embedded in the mice R package (van Buuren and Groothuis-Oudshoorn 2011), with the four DNAmAgeAC metrics, the eight greenness metrics, and the main model covariates (other than the target variable of the imputation) used as predictors.

Subgroup analyses.

We performed subgroup analyses by educational level (secondary or below vs. higher than secondary, including vocational training and university), area-level SES (lowest 50% vs. highest 50%), age group (<60 vs.≥60y ), smoking behavior (never smoked vs. current or former smoker), DNA methylation-based smoking index (described below, ≤median vs.>median ) (Gao et al. 2017), and alcohol use (never drank vs. current or former drinker) to evaluate the potential modification of the greenness-DNAmAgeAC association by these factors. We tested the statistical difference in effect estimates between subgroups by adding an interaction term (e.g., education–greenness) to the main model.

Sensitivity analyses.

We performed several sensitivity analyses. First, we further adjusted for season (four categories: spring, summer, fall, winter) and survey year (continuous variable) at the time of the blood draw to account for potential seasonal variation and long-term trends in DNA methylation. Second, we modeled annual mean temperature and relative humidity as nonlinear (vs. linear) terms using natural cubic splines with 3 df (Xu et al. 2020). Third, we repeated models without adjusting for smoking and alcohol use because there was no clear evidence that these two variables were related to greenness. Fourth, we repeated the models without adjusting for cell-type proportions in order to generate effect estimates that reflected extrinsic differences in DNAmAgeAC associated with greenness. The extrinsic differences reflect both differences in DNAmAgeAC among different leukocyte subpopulations and intrinsic differences in DNAmAgeAC independent of leukocyte type (Horvath and Raj 2018). Finally, instead of adjusting for self-reported smoking as a categorical variable (i.e., never, former, or current smoker) we adjusted for a continuous DNA methylation-based smoking index based on nine CpGs, using the algorithm provided by Gao et al. (2017), which the authors reported to be a stronger predictor of frailty in a cross-sectional sample of German adults than self-reported smoking. In our study population, the smoking index tended to be higher in current smokers than in former or never smokers, although smoking index distributions overlapped among the three smoking groups based on self-report (Figure S2). We used fixed-effects meta-regressions to test whether effect estimates from sensitivity analyses differed significantly from corresponding estimates based on the main model (Xu et al. 2019).

Associations with components of GrimAge.

To evaluate the association between surrounding greenness and individual components of GrimAge, we fit separate models using Z -scores for each component as the outcome instead of GrimAge. Z -score was calculated as (actual value – mean value)/ standard deviation.

All data analyses were performed using R software (version 3.5.1; R Development Core Team) with packages nlme (linear mixed effect model), splines (nature cubic splines), mvmeta (meta-regression), and mice (multiple imputation by chained equations). In this analysis, a two-sided p<0.05 was considered as statistically significant.

Results

We included 479 participants who resided in many different locations across Australia (Figure 1). All of them were female, including 132 monozygotic twins, 132 dizygotic twins, and 215 sisters of the twins from 130 families (Table 1). Their chronological ages ranged from 39.7 to 77.5 y (median= 55.7 y; IQR: 11.75 y). Among the four DNAmAgeAC estimates, the GrimAge acceleration showed the smallest variation [standard deviation (SD)= 3.12 y; minimum to maximum: −7.10 to 13.74 y; median= −0.44y ; IQR: 3.91 y], whereas the PhenoAge acceleration showed the largest variation (SD= 5.78 y; minimum to maximum: −15.45 to 21.25 y; median= 0.31 y; IQR: 8.17 y) (Table 2). The median NDVI and EVI values within 1,000m were 0.55 (IQR: 0.164) and 0.35 (IQR: 0.119), respectively. The distributions of NDVI and EVI were similar across the four buffer zones for each metric (Table 2), and the greenness indicators were highly correlated (Pearson correlation coefficients=0.85 –0.98 and 0.82–0.97 among NDVI and EVI measures, respectively, and 0.77–0.92 between NDVI and EVI measures) (Figure S3). Among the four DNAmAge estimates, the GrimAge showed the strongest association with chronological age (Pearson correlation coefficients=0.88 for GrimAge, 0.72 for Hannum’s age, 0.66 for Horvath’s age, 0.62 for PhenoAge) (Figures S3–S4). The four DNAmAge and DNAmAgeAC estimates were also correlated with each other (Pearson correlation coefficients=0.67 –0.80 between DNAmAge estimates, and 0.26–0.64 between DNAmAgeAC estimates). GrimAge acceleration was positively associated with the eight individual components of GrimAge, with Pearson correlation coefficients ranging from 0.19 for DNA methylation-based serum TIMP-1 to 0.80 for DNA methylation-based smoking pack years (Figure S5).

[Image omitted - see PDF]

Table 1 Basic characteristics of the 479 participants in Australia during 2004–2009.

Table 1, in two columns, lists Categories and Lowercase italic n (percentage).

Characteristics	n (%)
Women	479 (100.0)
Age group (y)
<60	323 (67.4)
≥60	156 (32.6)
Twin pairs
Dizygotic twins	132 (27.6)
Monozygotic twins	132 (27.6)
Sisters of twins	215 (44.9)
Education
Secondary or below	198 (41.3)
Vocational training	135 (28.2)
University	144 (30.1)
Missing	2 (0.4)
Marital status
Married or de facto	369 (77.0)
Never married	28 (5.8)
Widowed/separated/divorced	82 (17.1)
Area-level SES percentile in Australia
0, 25	88 (18.4)
25, 50	103 (21.5)
50, 75	136 (28.4)
75, 100	151 (31.5)
Missing	1 (0.2)
Smoking
Current	41 (8.6)
Former	147 (30.7)
Never	291 (60.8)
Alcohol use
Current	235 (49.1)
Former	52 (10.9)
Never	192 (40.1)

Note: SES, socioeconomic status.

Table 2 Descriptive statistics of continuous variables the 479 participants in Australia during 2004–2009.

Table 2, in six main columns, lists Categories, Mean plus or minus Standard Deviation, Minimum, Percentiles, Maximum, and Interquartile Range. The Percentiles column is subdivided into three columns, namely, Twenty-fifth, Fiftieth, and Seventy-fifth.

Variables	Mean±SD	Min	Percentiles			Max	IQR
			25th	50th	75th
Chronological age (y)	56.41±7.90	39.69	50.16	55.66	61.91	77.53	11.75
DNAmAge (y)
Horvath’s age	55.54±6.46	38.62	50.89	55.19	60.02	72.89	9.13
Hannum’s age	57.28±6.37	40.89	52.91	57.11	61.50	74.88	8.59
PhenoAge	53.04±7.37	34.92	47.96	52.75	57.73	72.79	9.78
GrimAge	54.20±6.46	39.64	49.51	53.86	58.54	72.76	9.04
DNAmAgeAC (y)
Horvath’s age acceleration	0.00±4.88	−16.33	−3.23	−0.15	3.02	15.40	6.26
Hannum’s age acceleration	0.00±4.41	−14.98	−2.75	0.21	2.66	18.42	5.40
PhenoAge acceleration	0.00±5.78	−15.45	−4.20	0.31	3.96	21.25	8.17
GrimAge acceleration	0.00±3.12	−7.10	−2.17	−0.44	1.74	13.74	3.91
Surrounding greenness
NDVI (buffer)
300m	0.56±0.12	0.25	0.47	0.55	0.64	0.87	0.175
500m	0.56±0.12	0.26	0.47	0.55	0.64	0.87	0.165
1,000m	0.56±0.12	0.26	0.48	0.55	0.64	0.87	0.164
2,000m	0.57±0.12	0.27	0.48	0.56	0.66	0.85	0.178
EVI (buffer)
300m	0.37±0.10	0.15	0.30	0.36	0.42	0.76	0.122
500m	0.36±0.10	0.14	0.30	0.35	0.42	0.77	0.121
1,000m	0.36±0.10	0.14	0.30	0.35	0.42	0.74	0.119
2,000m	0.37±0.10	0.13	0.29	0.35	0.43	0.70	0.134
Annual mean temperature (°C)	17.68±2.81	11.35	15.26	17.76	20.26	25.07	5.00
Annual mean relative humidity (%)	69.81±3.76	52.97	67.92	70.13	72.13	80.80	4.21

Notes: The chronological age was calculated as the difference between the date blood was taken and birth date. 2-y average, the average of the current year and previous year prior to blood being taken; DNAmAge, DNA methylation age; DNAmAgeAC, DNA methylation age acceleration; EVI, Enhanced Vegetation Index. IQR, interquartile range; max, maximum; min, minimum; NDVI, Normalized Difference Vegetation Index; SD, standard deviation.

GrimAge acceleration showed consistent significant (all p<0.05 ) inverse associations with NDVI and EVI metrics across all buffers (Figure 2). Each IQR increase in NDVI value within 1,000m was associated with a 0.59-y decrease (95% CI: −1.01 , −0.18 ; p=0.005 ) in GrimAge acceleration, whereas the effect estimates were slightly weaker for other greenness metrics. However, none of the greenness exposure metrics were significantly associated with DNAmAgeAC based on Horvath’s age, Hannum’s age, or PhenoAge (all p>0.10 ) in the study sample as a whole.

[Image omitted - see PDF]

There were significant interactions between surrounding greenness and area-level SES on Horvath’s and Hannum’s DNAmAgeAC. In general, IQR increases in surrounding greenness were associated with increased Horvath’s and Hannum’s DNAmAgeAC among participants whose area-level SES was below the median of Australia [e.g., for an IQR increase in EVI within 300m , 1.21 y (95% CI: 0.14, 2.28) and 0.87 y (95% CI: 0.05, 1.70), respectively], whereas corresponding estimates suggested slower age acceleration with increased greenness among those with a higher area-level SES [i.e., −0.43y (95% CI: −1.35 , 0.50; pinteraction=0.03 ) for Horvath’s DNAmAgeAC and −0.86y (95% CI: −1.64 , −0.08 ; pinteraction=0.004 ) for Hannum’s DNAmAgeAC] (Figure 3; Table S2). For GrimAge acceleration, there was a consistent pattern of inverse associations with greenness among those with a higher area-level SES, and null or weaker inverse associations among those with a lower area-level SES, although the differences were not statistically significant (pinteraction=0.10–0.66 ). Associations with PhenoAge acceleration varied, without consistent patterns by greenness metric or area-level SES (pinteraction≥0.28 ).

[Image omitted - see PDF]

Differences in associations with Horvath’s and Hannum’s DNAmAgeAC by education were somewhat consistent with differences by area-level SES, with increased greenness tending to have positive or null associations with age acceleration among those with lower education, and inverse or null associations among those with a higher education, although the differences were not statistically significant (pinteraction≥0.08 and ≥0.23 , respectively) (Figure S6, Table S2). For GrimAge acceleration, associations were very similar between those with lower and higher education (pinteraction≥0.26 ), whereas associations with PhenoAge acceleration varied, without clear differences by education (pinteraction≥0.17 ).

Differences in associations between greenness and DNAmAgeAC by smoking varied depending on the DNAmAgeAC metric and the variable used to classify smoking. For Horvath’s and Hannum’s DNAmAgeAC, associations tended to be positive or null among current and former smokers (combined) and inverse or null among never smokers (pinteraction= 0.04–0.13 and 0.03–0.24 for Horvath’s and Hannum’s metric, respectively) (Figure S7, Table S2). However, associations followed an opposite pattern for the DNA methylation-based smoking index, such that associations tended to be positive or null among those with a higher smoking index (above the median), and inverse or null for those with a lower smoking index (pinteraction= 0.07–0.39 and 0.08–0.26, respectively) (Figure S8). In contrast, associations with GrimAge DNAmAgeAC were inverse and similar between current or former smokers and never smokers (pinteraction≥0.28 ) and between those with a high or low smoking index (pinteraction≥0.28 ) (Figures S7–S8, Table S2.) Associations with PhenoAge DNAmAgeAC varied by greenness metric but were similar by smoking status (pinteraction≥0.19 ) and high or low smoking index (pinteraction≥0.50 ).

Associations between greenness metrics and DNAmAgeAC varied somewhat by age (<60 or≥60y ), but the differences were not statistically significant (all pinteraction≥0.15 ) and there were no consistent patterns by age among the DNAmAgeAC outcomes (Figure S9, Table S2). In general, associations between greenness and DNAmAgeAC were similar between those who never used alcohol and those who were current or former users, regardless of the DNAmAgeAC metric used (all pinteraction≥0.18 ) (Figure S10, Table S2).

Z -scores for five of the individual components used to derive the overall GrimAge metric (DNA methylation-based predictors of smoking pack years, serum cystatin-C, serum GDF-15, serum beta-2-microglobulin, and serum PAI-1) were inversely associated with all of the greenness metrics (Figure 4; Table S3). The associations were strongest for DNA methylation-based surrogates of smoking pack years (significant for all exposures), followed by serum cystatin-C and serum GDF-15. Of the remaining components, associations with exposures were mixed (inverse, null, and positive) for methylation-based estimators of serum adrenomedullin, close to the null for serum TIMP-1, and null or positive for serum leptin.

[Image omitted - see PDF]

Effect estimates were similar to the main model (pdifferences≥0.45 ) in sensitivity analyses of models adjusted for season and survey year, for annual mean temperature and relative humidity as nonlinear terms, and for a DNA methylation-based smoking index instead of self-reported smoking, as well as for models that were not adjusted for self-reported smoking and alcohol use or for cell-type proportions (Figure S11, Table S4).

Discussion

To the best of our knowledge, this is the first study to evaluate the association between surrounding greenness and measures of biological aging based on DNA methylation. Higher surrounding greenness was associated with slower biological aging, as measured by GrimAge acceleration. Of the eight components used to derive GrimAge, DNA methylation-based estimators of smoking pack years, serum GDF-15, and serum cystatin-C showed the strongest associations across both greenness metrics in all buffers when analyzed as separate outcomes. The association between higher surrounding greenness and slower biological aging measured by Horvath’s, Hannum’s, and PhenoAge DNAmAgeAC were close to the null in the population as a whole. Some significant positive associations were estimated for Horvath’s and Hannum’s age acceleration metrics among those with lower area-level SES, and some significant associations were estimated for Hannum’s metric among those with a higher area-level SES, whereas no significant associations were estimated for PhenoAge acceleration.

Given the well-documented associations between higher surrounding greenness and a lower risk of many aging related conditions (e.g., cognitive decline, cardiovascular diseases, type 2 diabetes) (de Keijzer et al. 2018, 2020; Fong et al. 2018; James et al. 2015; Jia et al. 2018; Lachowycz and Jones 2011; Rojas-Rueda et al. 2019; Sarkar 2017), we hypothesized that high surrounding greenness would slow the rate of biological aging. Our findings for GrimAge support this hypothesis in our study population. This is generally consistent with a cross-sectional study from Hong Kong, China, which reported that average telomere length was longer in residents of a 69-km2 region with more green space than in residents of four highly populated regions (area size: 7–10 km2 ) with little green space after adjusting for age, smoking, SES, and physical activity level (Woo et al. 2009). However, although telomere length has been interpreted as a marker of biological aging, DNAmAge is more strongly associated with mortality and morbidity (Jylhävä et al. 2017) and is thought to provide a more accurate measure of age acceleration, particularly in differentiated cells (Horvath and Raj 2018).

Based on meta-analyses of longitudinal data from four cohorts (6,395 adults followed for 13.7 y on average), Lu et al.(2019) estimated that each 1-y increase in GrimAge acceleration was associated with a 10% [hazard ratio (HR)=1.10 (95% CI: 1.09, 1.12)] increase in all-cause mortality. Using this information and the estimated difference in GrimAge acceleration with an IQR increase in NDVI within 500m in our study population (–0.51y per 0.165-unit increase in NDVI), we estimate that a 0.1-unit increase in NDVI within 500m would be associated with a 3% (i.e., HR=1.10−0.51×0.1/0.165=0.97 ) decrease in all-cause mortality. This estimate is consistent with a recent meta-analysis that reported a pooled HR of 0.96 (95% CI: 0.94, 0.97) for all-cause mortality associated with every 0.1-unit increase in NDVI within 500m (Rojas-Rueda et al. 2019).

Although all measures of greenness exposure were associated with significantly slower biological aging based on the GrimAge DNAmAgeAC metric, associations with the other measures of biological age acceleration were variable (depending on the greenness measure or study subpopulation) or null. These apparent inconsistencies were not unexpected given recent evidence suggesting that GrimAge is a better indicator of biological aging. For example, among the four DNAmAge estimators, DNAmAgeAC estimated by GrimAge showed the strongest association (GrimAge>PhenoAge>Hannum>Horvath ) with aging-related conditions (death, cancer, coronary heart disease, and age at menopause) based on an analysis of longitudinal data from four cohorts (6,395 participants followed for 13.7 y on average) (Lu et al. 2019). Similarly, a longitudinal analysis of data from two cohorts of older individuals (n=2,462 , mean follow-up of ∼9y ) reported that, of the four DNAmAgeAC metrics assessed, GrimAge acceleration had the strongest and most consistent associations with mortality, myocardial infarction, and stroke (Wang et al. 2021). GrimAge acceleration also showed the strongest association with fatty liver and excess visceral fat in 1,177 adults from the Framingham Heart Study (Lu et al. 2019). A recent analysis of data from three cross-sectional surveys of UK adults 45–87 years of age (n=1,388–1,685 , depending on the outcome) also reported that GrimAge acceleration was the strongest predictor of physical and cognitive functioning, as assessed by grip strength, forced expiratory volume in 1 s (FEV1 ), mental speed, and episodic memory (Maddock et al. 2020). For example, a 1-y increase in GrimAge acceleration was associated with a 0.03-mL (95% CI: –0.05 , –0.01 ) lower mean FEV1 , whereas corresponding estimates for Horvath’s, Hannum’s, and PhenoAge DNAmAgeAC metrics were 0.00 (95% CI: –0.01 , 0.01), 0.00 (95% CI: –0.02 , 0.01), and –0.01 (95% CI: –0.02 , –0.01 ), respectively (Maddock et al. 2020). GrimAge acceleration has also shown stronger associations with potential determinants of aging than other metrics, including smoking and alcohol use in a cross-sectional analysis of 1,100 adults in the United States (Zhao et al. 2019) and physical activity and neighborhood SES in cross-sectional analyses of 2,700 women in the United States (Kresovich et al. 2020; Lawrence et al. 2020). Taken together, these findings suggest that GrimAge DNAmAgeAC may be a better measure of biological aging than the other DNAmAgeAC metrics examined in our analysis.

Of the eight individual DNA methylation-based components of GrimAge, three had the strongest and most consistent associations with greenness when modeled as individual outcomes: a) DNA methylation-based markers of smoking pack years; b) GDF-15; and c) cystatin-C. DNA methylation-based smoking pack years represents a DNA methylation signature resulting from prolonged exposure to cigarette smoke (Lu et al. 2019), and this biomarker can be reversed following smoking cessation (McCartney et al. 2018). A plausible explanation for our results is that high greenness may help reverse the DNA methylation alterations resulting from exposure to cigarette smoke, given that high greenness is related to reduced risk of many adverse outcomes (e.g., cardiovascular diseases, cognitive decline, mental illness) that can be caused by smoking (Banks et al. 2019; Deal et al. 2020; Fong et al. 2018; Wootton et al. 2020).

GDF-15 is a promising biomarker for disease prognosis and an emerging target for cancer immunotherapy because of its important immune-regulatory function (Wischhusen et al. 2020). Serum levels of GDF-15 often increase with age, and highly elevated serum levels of GDF-15 are mostly linked to inflammation, myocardial ischemia, and cancer (Wischhusen et al. 2020). GDF-15 is also a potential target in the regulation of body weight and energy expenditure because high levels of GDF-15 contributes to appetite loss (Tsai et al. 2016). The DNA methylation-based GDF-15 was built to predict serum levels of GDF-15 (Lu et al. 2019). In our study, its inverse association with greenness may reflect the documented associations between high greenness and improved immune function (Kuo 2015) or improved metabolic health (e.g., lower adiposity) (Sarkar 2017).

High serum levels of cystatin-C is a widely used clinical indicator of reduced kidney function or declined glomerular filtration rate (Di Somma and Marino 2019). In our study, the DNA methylation-based cystatin-C, a predictor of serum level of cystatin-C (Lu et al. 2019), was inversely associated with greenness, although the association was only statistically significant for EVI within 300m . This may suggest that high greenness may be associated with improved kidney function. However, current evidence about the relationship between greenness and kidney diseases remains scarce (Kondo et al. 2018). A prospective cohort study in 108,630 women from the U.S.-based Nurses’ Health Study followed up between 2000 and 2008 found a relationship between high greenness and reduce mortality risk from kidney diseases, but the association was not statistically significant [HR=0.63 (95% CI: 0.38, 1.04) for each 0.1-unit increase in NDVI within 250m ], possibly due to the small number of cases (139 cases) (James et al. 2016). The other potential explanation for our results for cystatin-C might be that high greenness is associated with lower adiposity (Sarkar 2017). Adipose tissue was found to be among the tissues expressing the highest level of cystatin-C mRNA, which contribute to the higher serum levels cystatin-C seen in obese people compared with nonobese people (Naour et al. 2009).

Patterns of associations between greenness and age acceleration according to neighborhood SES varied among the DNAmAgeAC metrics, with some significant differences for Horvath’s and Hannum’s metrics resulting from positive associations among those with lower neighborhood SES and inverse associations among those with a higher neighborhood SES. There were no significant differences by neighborhood SES for the GrimAge metric, although associations were consistently inverse for those with higher neighborhood SES and closer to the null for those with lower SES. Associations between greenness and Hannum’s and GrimAge DNAmAgeAC were similar between those with higher and lower education, a proxy measure of individual-level SES, whereas associations with Horvath’s DNAmAgeAC tended to be positive for those with lower education and null for those with higher education (with a significant difference for only one exposure metric.) In a Swiss cohort, the association between surrounding greenness and lower mortality was stronger among participants living in communities with a higher area-level SES (Vienneau et al. 2017), whereas in residents of Rome, Italy, residential greenness was associated lower mortality in the population as a whole but with higher mortality among residents in the lowest quintile of area-level SES (Orioli et al. 2019). Another study reported that a positive association between neighborhood socioeconomic disadvantage and a DNA methylation-based mortality risk index was attenuated in areas with high neighborhood greenness (Ward-Caviness et al. 2020). Differences in associations between greenness and biological age acceleration or other health outcomes could be related to differences in greenspace quality or accessibility between communities with higher SES or lower SES (Vienneau et al. 2017). Higher neighborhood SES might also be associated with a stronger social cohesion (Feldman and Steptoe 2004) that could boost the health benefits of green space (Jennings and Bamkole 2019).

The twin and family design is a strength of our study. This design compares each participant with her siblings, thus reducing potential confounding by genetic background and shared familial factors (e.g., childhood living environment, familial dietary behavior) (Li et al. 2017; Stone et al. 2010). This could be an important advantage given that familial factors explain a substantial proportion of variation in DNAmAge (Horvath and Raj 2018; Li et al. 2020; Lu et al. 2019).

Several limitations also should be acknowledged. This is a cross-sectional analysis, and further longitudinal studies are warranted to verify our findings. Our sample size may have limited statistical power to estimate precise and stable greenness–DNAmAgeAC associations and potential effect modifiers. In addition, the study included only women. Previous studies reported that associations between surrounding greenness and decreased mortality and body fat were weaker in men than women (Ji et al. 2019; Sarkar 2017). Associations between greenness and DNAmAgeAC also might be weaker in men.

In conclusion, higher surrounding greenness was associated with slower biological aging measured by GrimAge acceleration in Australian women. The association was also evident for three individual components of GrimAge and seemed to be stronger among those with higher vs. lower neighborhood SES, although the differences were not statistically significant. Associations between greenness and biological aging measured by Horvath’s and Hannum’s DNAmAgeAC were less consistent, and varied depending on neighborhood SES. Our study and previous epidemiological evidence suggests that GrimAge may be a more reliable measure of biological aging than other DNAmAge metrics, but additional studies are needed to confirm our findings in other populations.

Acknowledgments

R.X. is supported by China Scholarship Council (201806010405). Shuai Li is supported by an Early Career Fellowship of the Australian National Health and Medical Research Council (NHMRC; APP1109193). Y.G. is supported by a Career Development Fellowship of the Australian National Health and Medical Research Council (APP1163693). Shuai Li is supported by an Early Career Research Fellowship of the Victorian Cancer Agency (ECRF 19020). M.C.S. is a NHMRC Senior Research Fellow (APP1155163). J.L.H. is a NHMRC Senior Principal Research Fellow. The Australian Mammographic Density Twins and Sisters Study (AMDTSS) was facilitated through access to Twins Research Australia, a national resource supported by a Centre of Research Excellence Grant (1079102) from the NHMRC. The AMDTSS was supported by NHMRC (1050561 and 1079102), Cancer Australia and National Breast Cancer Foundation (509307). The raw and processed DNA methylation data set are open accessed or free on Gene Expression Omnibus (accession number GSE100227). As required by the ethics approval, the authors are not allowed to open other data (e.g., data on covariates) used in this study in order to protect the privacy of participants. If anyone wants to use the data to repeat our analyses or to perform other research, please contact J.L.H. ([email protected]), who is happy to engage with external collaborators. This person would have to be added to the ethics approval to access the data. Working with the AMDTSS group as collaborators would likely result in better science given this is family data and the nuances of this, as well as the sampling issues, need to be understood to guard against making false conclusions from a naïve use of data without such knowledge.