This study analyzed data from the 2005–2010 NHANES cycles, spanning three consecutive survey periods. Although BFRs data is available through 2016, a methodology changes in C-reactive protein (CRP) measurement occurred in the 2011–2012 cycle with the introduction of high-sensitivity CRP. Since CRP is a critical marker for calculating phenotypic age, data after 2010 were excluded to ensure methodological consistency.

NHANES, conducted by the CDC, is a nationally representative survey employing a complex, multi-stage sampling design to evaluate the health and nutrition of the U.S. population. Participants undergo interviews, physical exams, and biospecimen collection at mobile examination centers. For this analysis, individuals aged 20 and older were included, providing information on demographics, socio-economic status, dietary intake, chronic diseases, and BFRs exposure from the 2005–2010 cycles. After applying specific inclusion and exclusion criteria, a total of 3,908 participants were selected for the final analysis, as outlined in Figure 1.

BFRs were measured in serum using automated liquid–liquid extraction followed by gas chromatography coupled with high-resolution mass spectrometry. In line with previous studies, serum BFR concentrations were used to reflect individual exposure levels. The analysis included polybrominated biphenyl (PBB153) and eight polybrominated diphenyl ethers (PBDEs) with detection rates above 65% (20). The selected PBDEs were: 2,4,4′-Tribromodiphenyl ether (PBDE28), 2,2′,4,4´-Tetrabromodiphenyl ether (PBDE47), 2,2′,3,4,4´-Tetrabromodiphenyl ether (PBDE85), 2,2′,4,4′,5-Pentabromodiphenyl ether (PBDE99), 2,2′,4,4′,6-Pentabromodiphenyl ether (PBDE100), 2,2′,4,4′,5,5´-Hexabromodiphenyl ether (PBDE153), 2,2′,4,4′,5,6´-Hexabromodiphenyl ether (PBDE154), and Decabromodiphenyl ether (PBDE209). In this study, values below the limit of detection (LOD) were imputed using limit divided by the square root of two, in accordance with the NHANES analytical guidelines (21). This approach is widely adopted in environmental health research. Moreover, previous studies have assessed the sensitivity of various imputation methods for handling values below the detection limit, providing support for LOD/√2 as a reasonable and widely accepted estimation method (22).

Aging acceleration was assessed using phenotypic age, which is calculated based on various biomarkers and a specific algorithm. To control for chronological age, PhenoAgeAccel was determined as the residual from a regression of phenotypic age on chronological age. Participants were categorized into two groups: those with accelerated aging (PhenoAgeAccel ≥ 0) and those with delayed aging (PhenoAgeAccel < 0).

Phenotypic age was calculated using the following formula (14):

Where:

And:

𝑥𝑏 = −19.907–0.0336 × Albumin+0.0095 × Creatinine+0.1953 × Glucose+0.0954 × ln(CRP) − 0.0120 × Lymphocyte Percentage+0.0268 × Mean Cell Volume+0.3306 × Erythrocyte Distribution Width + 0.00188 × Alkaline Phosphatase+0.0554 × Leukocyte Count+0.0804 × Chronological Age. Here, “𝑥𝑏 “represents the linear combination of biomarkers used in the model.

Based on previous literature, covariates that have been proven to be associated with environmental factors and biological aging were included in the analysis. These covariates were: age, sex, race/ethnicity, poverty-to-income ratio (PIR), BMI, marital status, education level, smoking status, alcohol consumption, healthy dietary status (assessed using the Healthy Eating Index-2015 [HEI-2015] (23)), physical activity level, and the presence of chronic diseases, including hypertension, diabetes mellitus (DM), cardiovascular disease (CVD), and cancer, Additionally, the corresponding survey year was included. Hypertension and DM were determined through both self-reports and clinical measurements, while CVD and cancer were self-reported. Estimated glomerular filtration rate (eGFR) was calculated using the Chronic Kidney Disease Epidemiology Collaboration (CKD-EPI) equation. Participants were classified by PIR as low (<1.3), middle (1.3–3.5), and high (>3.5); BMI was categorized based on WHO classifications (normal, overweight, obesity); alcohol consumption was divided into non-drinkers, moderate drinkers (1–3 drinks/day), and heavy drinkers (≥4 drinks/day); and physical activity level was categorized as active, inactive, moderate, or others (24).

Although NHANES employs a complex sampling design, weighting was not applied in this study due to the limitations in incorporating weights into mixture analyses. For baseline characteristics, continuous variables were presented as means with standard deviations (SD), while categorical variables were expressed as counts and percentages. Differences in continuous variables between the delayed and accelerated aging groups were assessed using Student’s t-tests, and chi-square tests were employed for categorical variables.

Generalized Linear Models (GLMs) were used to evaluate the associations between BFRs(log-transformed) and PhenoAgeAccel. Three models were applied: the crude model, which was unadjusted; Model 1 adjusted for age, sex, and race/ethnicity; and Model 2 further adjusted for BMI, PIR, marital status, education, physical activity, smoking, alcohol consumption, hypertension, diabetes, cardiovascular disease, cancer, HEI-2015, and survey year. Restricted cubic spline (RCS) analysis was conducted to investigate the dose–response relationship between BFRs exposure and PhenoAgeAccel after adjusting for all confounders. Subgroup and interaction analyses were also performed based on age and sex.

In the second phase, quantile sum (WQS) regression was used to assess the effect of BFRs mixture exposure on PhenoAgeAccel. The bootstrap method (10,000 iterations) was applied to assign weights to each contaminant, identifying the most influential components in the mixture. Contaminant weights ranged from 0 to 1. This method not only captures real-world co-exposures more effectively but also enhances sensitivity in identifying key predictors compared to univariate analysis. The data were randomly split into 40% training and 60% testing samples, with 10,000 bootstrap iterations.

Additionally, Bayesian Kernel Machine Regression (BKMR) was employed to visualize the joint exposure-response relationship between BFRs mixtures and biological aging risk. BKMR’s strength lies in its flexibility to model exposure-response relationships, accounting for potential nonlinear and non-additive effects often seen in environmental epidemiology. The exposure-response functions were modeled using a Gaussian process, and the analysis was run with 10,000 iterations of the Markov Chain Monte Carlo (MCMC) algorithm. Conditional posterior inclusion probabilities (conPIPs) were calculated to identify the key BFRs components most associated with biological aging risk. All statistical analyses were conducted by R version 4.4.12.

A total of 3,908 U.S. adults were included in the study and categorized into two groups based on PhenoAgeAccel: the delayed aging group (n = 3,004) and the accelerated aging group (n = 904). Participants in the accelerated aging group were significantly older, with a mean age of 54.78 ± 17.46 years, compared to 41.30 ± 17.85 years in the delayed aging group. The accelerated aging group also had a higher proportion of non-Hispanic Black individuals, lower PIR, lower HEI-2015 scores, reduced eGFR, lower education levels, and a lower rate of active physical activity. In contrast, this group exhibited higher BMI, a greater prevalence of widowed/divorced/separated individuals, and an increased incidence of smoking, hypertension, cancer, and diabetes. Furthermore, the accelerated aging group was exposed to higher levels of PBDE28, PBDE47, PBDE85, PBDE99, PBDE100, and PBDE154. Detailed information is provided in Table 1. Besides, the characteristics of study participants by sex in adults are presented in Supplementary Table S1, showing that BFRs exposure levels were significantly higher in males compared to females. Additionally, the characteristics of study participants with accelerated phenotypic age by sex are provided in Supplementary Table S2, where certain BFRs exhibited significantly higher exposure levels in males than in females.

Table 2 presents the results of the GLMs, revealing significant positive associations between several BFRs and PhenoAgeAccel. Specifically, PBDE28 (β = 0.55, 95% CI: 0.13, 0.96), PBDE85 (β = 0.42, 95% CI: 0.11, 0.74), PBDE47 (β = 0.39, 95% CI: 0.04, 0.75), PBDE99 (β = 0.38, 95% CI: 0.07, 0.68), and PBDE154 (β = 0.37, 95% CI: 0.04, 0.70) were all significantly associated with PhenoAgeAccel.

RCS analysis confirmed both linear and nonlinear dose–response relationships for several PBDEs (Figure 2). Specifically, PBDE47 (P for overall = 0.0098, P for nonlinearity = 0.0336) and PBDE99 (P for overall = 0.0041, P for nonlinearity = 0.0233) exhibited significant nonlinear associations with PhenoAgeAccel, indicating that increases in exposure were associated with sharper increases in PhenoAgeAccel at higher exposure levels. In contrast, PBDE28 (P for overall = 0.0225, P for nonlinearity = 0.3143), PBDE85 (P for overall = 0.0118, P for nonlinearity = 0.2902), and PBDE154 (P for overall = 0.0432, P for nonlinearity = 0.2062) demonstrated significant overall associations without evidence of nonlinearity, suggesting that their relationships with PhenoAgeAccel were more likely linear.

The stratified analysis by age and sex revealed differential associations between BFRs exposure and PhenoAgeAccel (Supplementary Table S3). Males were more susceptible to the effects of PBB153, PBDE28, PBDE85, PBDE47, PBDE99, and PBDE154 compared to females, with significant sex interactions (P for interaction <0.0001, 0.002, 0.003, 0.017, 0.003, and 0.025, respectively). Additionally, PBDE99 showed a stronger effect in participants aged 60 and older (P for interaction = 0.027), while PBDE209 had a more pronounced effect in the 40–59 age group (P for interaction = 0.005).

When constraining the analysis in the positive direction, the WQS regression analysis (Figure 3) revealed a significant association between combined exposure to BFRs and an increase in PhenoAgeAccel (β = 0.504, 95% CI: 0.071, 0.937), with this effect being particularly pronounced among males. The chemicals that contributed most to this increase in PhenoAgeAccel were PBB153, PBDE153, PBDE28, and PBDE85. Interaction analyses indicated that males were more strongly affected by this exposure. Conversely, in the negative direction, no significant association was found (β = 0.424, 95% CI: −0.084, 0.932).

The BKMR analysis (Figure 4A) indicated as exposure to BFRs increases from lower to higher quantiles, there is a noticeable trend in the estimated effects, suggesting that higher mixed BFRs exposure may be associated with increased PhenoAgeAccel, particularly in males (red circles). The posterior inclusion probabilities (PIPs) identified PBDE99 as the primary drivers of this effect, Other BFRs, such as PBDE28 and PBDE85, also show significant probabilities, but to a lesser extent than lnPBDE99 (Figure 4B). Males exhibited a stronger response to the BFRs mixture compared to females. Supplementary Figure S2A shows the univariate exposure-response analysis for each BFR, indicating a positive monotonic relationship between PBDE100, PBDE85, and PBDE99 with PhenoAgeAccel. In contrast, PBDE28, PBDE47, and PBDE154 exhibited a negative monotonic trend. As shown in Supplementary Figure S2B, when other chemicals were held constant at the 25th, 50th, and 75th percentiles, a positive association between PBDE100, PBDE85, and PBDE99 and PhenoAgeAccel was observed. Additionally, no statistically significant interactions between BFR components were detected (Supplementary Figure S2C).