Associations of urinary phthalate metabolites with DNA methylation algorithms of aging among U.S. adults: 1999-2002 NHANES data

Study participants

NHANES is an annual, nationally representative survey of non-institutionalized civilians in the U.S., conducted by the Centers for Disease Control and Prevention (CDC) through the National Center for Health Statistics (NCHS)^[33]. Utilizing a complex sampling design that includes oversampling of specific subgroups^[34], NHANES gathers demographic and health data from participants who provide blood and urine samples at mobile examination centers and give written consent^[35]. For our study, we combined data from two cycles (1999-2002) in accordance with NCHS recommendations.

Quantification of urinary PAE metabolite concentrations

Using data from the publicly available NHANES database, we analyzed urinary concentrations of several PAE metabolites: mono-ethyl phthalate (MEP), mono-n-butyl phthalate (MBP), mono-cyclohexyl phthalate (MCHP), mono-(2-ethyl) hexyl phthalate (MEHP), mono-n-octyl phthalate (MOP), mono-benzyl phthalate (MBzP), mono-isononyl phthalate (MiNP), mono-(2-ethyl-5-hydroxyhexyl) phthalate (MEHHP), mono-(2-ethyl-5-oxohexyl) phthalate (MEOHP), mono-(3-carboxypropyl) phthalate (MCPP), and mono-isobutyl phthalate (MiBP). Quantification of urinary PAE metabolite concentrations followed NHANES guidelines^[36,37]. In summary, urine samples underwent enzymatic deconjugation of glucuronides, followed by solid-phase extraction. Chromatographic separation of PAE metabolites was achieved using reversed-phase high-performance liquid chromatography (HPLC), with detection performed using atmospheric pressure chemical ionization tandem mass spectrometry (APCI-MS/MS). Metabolites were quantified through isotope dilution, and assay precision was enhanced by incorporating ¹³C₄-labeled internal standards for each analyte, along with a conjugated internal standard to assess deconjugation efficiency. Urinary concentrations of the PAEs below the level of detection (LOD) were assigned the limit of detection divided by the square root of two as recommended by NHANES^[38]. This standardized approach, consistently adopted in peer-reviewed analyses of NHANES 1999-2002 biomarker data^[39], ensures methodological alignment with established protocols for handling left-censored observations in this cohort. We included only PAE metabolites with detection frequencies ≥ 75%. This detection frequency (75%) was selected following the established methodology for urinary PAE analyses in the NHANES 1999-2002 dataset, as implemented in the study of PAEs and leukocyte telomere length^[39].

DNAm aging algorithm measurements

Following NHANES DNA Methylation Array and Epigenetic Biomarkers Data Documentation^[40], we assessed DNAm using blood samples from adults aged ≥ 50 years. These samples were collected during the 1999-2000 and 2001-2002 survey cycles. Methylation profiling employed the Illumina EPIC array (San Diego, CA, USA), with preprocessing and normalization conducted per manufacturer protocols in Dr. Yongmei Liu’s laboratory at Duke University.

We computed five established epigenetic aging biomarkers: HannumAge^[41], SkinBloodAge^[42], PhenoAge^[43], GrimAge^[44], and GrimAge2^[45]. For each algorithm, we calculated age-adjusted residuals (epigenetic age acceleration metrics) to quantify deviations from chronological age. Designated as HorvathAgeacc, SkinBloodAgeacc, PhenoAgeacc, GrimAgeMortacc, and GrimAge2Mortacc, these metrics indicate accelerated or decelerated aging, providing insight into biological aging relative to chronological age.

Covariates

To control for potential confounders, we included the following covariates: demographic characteristics (age in years, continuous; sex), lifestyle factors (smoking pack-years and alcohol consumption), leukocyte composition, urine creatinine, and C-reactive protein levels, along with survey years. As previously described^[46,47], age, sex, and alcohol consumption data were collected through self-report questionnaires. We categorized alcohol consumption as Never, Ever, Current, and Don’t Know based on the question: “In any one year, have you had at least 12 drinks of any type of alcoholic beverage?” Smoking pack-years were estimated using a DNAm-based method^[44], while leukocyte composition was determined via the Houseman algorithm^[48]. Urinary creatinine (ng/mL) was measured via a Jaffe rate reaction on a CX3 analyzer and entered into models as a natural log-transformed variable.

Statistical analysis

We summarized baseline demographic characteristics using standard descriptive methods, with analyses weighted for urinary PAE subsamples per NCHS recommendations. Using Spearman correlation coefficients, we assessed correlations among age, urinary PAE metabolites, and the five DNAm aging algorithms, presenting results in a correlation matrix plot.

We examined associations between urinary PAE metabolites and DNAm aging algorithms using linear regressions. Model 1 was adjusted for age, sex, leukocyte composition, and survey years, while Model 2 included additional adjustments for smoking pack-years, alcohol consumption, natural log-transformed urine creatinine, and natural log-transformed C-reactive protein. We also stratified analyses by sex and age to explore demographic-specific associations specific to demographic groups, which could further elucidate how PAEs impact biological aging. To rigorously evaluate the assumption of linear dose-response relationships between PAEs and DNAm aging algorithms, we employed generalized additive models (GAMs) using the “mgcv” package (v1.8-42) in R.

All analyses were performed using SAS, version 9.4 (SAS Institute, Inc., Cary, NC) and R, version 4.3.1. Statistical significance was defined as a two-sided P-value < 0.05.

Population characteristics

Table 1 presents the characteristics of the study population, which consists of 831 participants with an equal distribution of males and females. The average age was 66.66 years in the 1999-2000 cycle and 65.13 years in the 2001-2002 cycle, with no significant difference between the two cycles. Additionally, alcohol consumption, body mass index, smoking pack-years, and C-reactive protein levels were comparable across the two cycles.

Profiles of PAE metabolites in urine

Table 2 summarizes the distribution of urinary PAE metabolites in the study population and lists the LOD values provided by NHANES^[49]. MEP, MEHP, MBzP, and MBP were highly detectable across all participants, with detection frequencies of 99%, 78%, 96%, and 99%, respectively. The average concentrations of MEP, MEHP, MBzP, and MBP were 551.86, 6.31, 16.26, and 40.39 ng/mL, respectively. MEP and MBzP concentrations were significantly higher in males than in females. MEP levels also differed significantly between participants aged 50-64 years (average: 669.55 ng/mL) and those aged ≥ 65 years (average: 443.97 ng/mL). This finding aligns with a notable negative correlation between age and MEP (P-value = 0.0005, r = -0.12), suggesting that MEP exposure declines with increasing age [Figure 1]. These results indicate that MEP is the predominant PAE metabolite in humans, and males aged 50-64 years may face greater health risks from MEP exposure.

Figure 1. Spearman correlation among age, urinary phthalate metabolites, and the five DNAm aging algorithms in the entire participant population. ^***P < 0.001; ^**P < 0.01; and ^*P < 0.05. DNAm: DNA methylation.

Correlations of age, PAE metabolites, and DNAm algorithms

Assessments of acceleration values from the five DNAm algorithms showed 0.89 for HannumAgeacc, -1.50 for SkinBloodAgeacc, -10.47 for PhenoAgeacc, 0.31 for GrimAgeMortacc, and 6.15 for GrimAge2Mortacc [Table 3]. Significant differences were observed: HannumAgeacc, GrimAgeMortacc, and GrimAge2Mortacc were higher in males than in females, and individuals aged 50-64 years exhibited greater accelerations compared to those aged ≥ 65 years. MBzP showed notable correlations with GrimAgeMortacc (P-values = 0.0001, r = 0.13) and GrimAge2Mortacc (P-values < 0.0001, r = 0.14) [Figure 1].

Associations between PAE metabolites and acceleration of DNAm algorithms

Table 4 presents the associations between MEP, MEHP, MBzP, and MBP with the five DNAm algorithms. After multivariable adjustments, MBzP exhibited significant associations with all DNAm algorithms, with P-values ranging from 0.0001 to 0.0393 and β coefficients varying from 0.22 (95%CI: 0.07-0.38) for GrimAgeMortacc to 0.68 (95%CI: 0.18-1.18) for PhenoAgeacc. In contrast, MEP, MEHP, and MBP did not demonstrate significant associations with any DNAm algorithm.

Sex-stratified analyses [Table 5] revealed that MBzP was significantly associated with PhenoAgeacc, GrimAgeMortacc, and GrimAge2Mortacc in males, with the strongest effect observed for PhenoAgeacc (β = 0.75, 95%CI: 0.12-1.38, P-value = 0.0205). Among females, a significant association was identified only between MBzP and GrimAge2Mortacc (β = 0.36, 95%CI: 0.05-0.66, P-value = 0.0250). These findings are corroborated by Supplementary Figure 1, which illustrates that MBzP and other metabolites correlate with GrimAgeMortacc and GrimAge2Mortacc in males [Supplementary Figure 1A], whereas no associations were observed in females [Supplementary Figure 1B].

Age-stratified analyses [Table 6] indicated significant associations between MBzP and all five DNAm algorithms in participants aged 50-64 years (P-values = 0.0012-0.0155), with β values ranging from 0.30 (95%CI: 0.10-0.51) for GrimAgeMortacc to 0.93 (95%CI: 0.34-1.52) for PhenoAgeacc. In the ≥ 65 age group, significant associations were observed between MEP and GrimAgeMortacc (β = -0.17, 95%CI: -0.31 to -0.03, P-values = 0.0212) and between MBzP and GrimAge2Mortacc (β = 0.33, 95%CI: 0.05-0.62, P-values = 0.0239). These associations are emphasized in Supplementary Figure 2, which shows correlations for MBzP across all five DNAm algorithms in the 50-64 age group [Supplementary Figure 2A]. Overall, significant associations between GrimAge2Mortacc and MBzP were observed in the general sample, with particularly stronger correlations between PhenoAgeacc and MBzP among males (P-interaction < 0.0001) and those aged 50-64 (P-interaction < 0.0001). GAMs detected no significant nonlinearity in the associations between PAEs and epigenetic aging (all nonlinearity P > 0.05), validating the use of linear regression in this analysis.

Discussion

Our study demonstrates that urinary MBzP is significantly associated with accelerated epigenetic aging, particularly affecting GrimAge-based algorithms in the overall population and PhenoAgeacc in males aged 50-64 years. This association persisted after multivariable adjustment and remained robust in age- and sex-stratified analyses. Notably, MBzP exerted stronger effects in males and middle-aged participants (50-64 years), suggesting demographic-specific susceptibility. In contrast, despite their high detection rates, MEP, MEHP, and MBP showed no significant associations with DNAm aging biomarkers.

This demographic-specific association points to intersecting biological vulnerabilities. First, MBzP exposure is linked to anti-androgenic effects in humans, as evidenced by its correlation with delayed pubarche in girls^[50]. This property may synergize with age-related testosterone decline in males, creating a unique window of susceptibility^[51]. Second, MBzP directly induces oxidative stress in human erythrocytes by increasing reactive oxygen species (ROS) and hydroxyl radicals at concentrations ≥ 5 μg/mL, while simultaneously altering superoxide dismutase and glutathione peroxidase activity^[52]. Middle-aged males demonstrate reduced plasma antioxidant capacity, which could exacerbate MBzP-induced DNA damage^[53]. This is supported by epidemiological evidence showing that the association between MBzP and uterine fibroid risk is mediated by oxidatively generated DNA damage (8-OHdG)^[54].

These findings provide context for prior research on PAE. While Oluwayiose et al. reported broad links between PAEs and epigenetic age in sperm, with MBzP showing only a non-significant positive trend (β = 0.24, P = 0.10)^[32], our study identifies MBzP as a significant driver of systemic aging in blood-derived DNAm algorithms among vulnerable demographics. This divergence highlights tissue-specific effects and underscores MBzP’s distinct role in systemic biological aging. The association between MBzP and aging measures appears to vary across populations. In elderly Koreans (samples 2012-2014), MBzP had a lower detection frequency (731 out of 1,190) but a higher average urinary concentration (2,190 ng/mL)^[54]. In this elderly group, linear correlations were observed between urinary MBzP quartiles and measures of physical performance^[54]. In contrast, a recent study of children (2017-2019) reported a very high MBzP detection frequency (99.5%) with a lower average concentration (11.8 ng/mL)^[55]. Furthermore, compounds with lower detection frequencies (< 75%) in our study (MEHHP, MEOHP, MCPP, MiBP) were highly detectable in children’s urine (100%, 100%, 99.4%, 100%) with average concentrations of 28.5, 6.42, 1.78, and 20.1 ng/mL, respectively^[55]. The attenuated MBzP-aging association in elderly Koreans contrasts with heightened susceptibility in children, further emphasizing population-specific vulnerabilities, where age-dependent exposure sources and metabolic shifts modulate risks^[55,56]. This aligns with observations of progressive declines in antioxidant capacity during aging, which could amplify susceptibility to MBzP-induced oxidative stress in midlife^[53].

The mechanistic basis for MBzP’s potency involves its direct biological activities. First, MBzP induces oxidative stress by elevating ROS and hydroxyl radicals in human erythrocytes, leading to DNA damage (8-OHdG) that mediates disease pathogenesis^[52,54]. Second, MBzP exposure elevates systemic inflammation markers (C-reactive protein and interleukin-6) and promotes hyperuricemia through immune inflammation mediation^[57,58]. Third, MBzP correlates with altered DNAm patterns in human sperm and urine (this study), confirming its potential to disrupt epigenetic processes^[31]. These pathways converge on epigenetic aging mechanisms, through which MBzP-induced inflammation promotes GrimAge acceleration, a recognized mediator of cognitive impairment^[59]. The stronger association in males aged 50 to 64 likely reflects three synergistic factors: (1) Interaction between MBzP’s anti-androgenic effects and age-related hormonal changes^[50,53]; (2) Compromised antioxidant defenses in midlife males^[53]; and (3) Enhanced sensitivity to inflammaging, as GrimAge incorporates inflammatory mortality predictors^[59]. These findings shift focus from general PAE exposure to specific risk scenarios: MBzP exposure during midlife in males represents a critical window for the acceleration of biological aging. Future studies should validate these associations longitudinally while measuring oxidative stress and inflammation biomarkers to confirm mechanistic pathways.

Several limitations should be considered when interpreting the findings. First, as an observational study, it is susceptible to potential biases, including reverse causality and the influence of unmeasured or residual confounding factors (such as diet, occupation, or other environmental/lifestyle factors). These limitations restrict the ability to establish causal relationships between urinary PAE metabolites and DNAm-based aging algorithms. Second, the study sample comprised only U.S. participants aged 50 years and older. While the sample shares key sociodemographic characteristics with the broader U.S. population, it is not a random sample. This constrains the generalizability of the findings to younger populations, other geographical regions, or more diverse cohorts. Third, DNAm signatures were assessed at a single time point, limiting the ability to track their stability or changes over time. Longitudinal data would be necessary to better understand the temporal dynamics of these signatures and their associations with health outcomes. Finally, further validation in independent studies involving more diverse populations with varying age ranges and risk profiles is essential to confirm the robustness and generalizability of these results.