Study population

The NHANES is a cross-sectional study conducted in the United States since the 1960s to assess the health and nutritional status of individuals, including children and adults. Approval was obtained from the Ethical Review Board of the National Center for Health Statistics, and informed consent was secured from all participants. Data from publicly available sources, spanning from 2003 to 2018, were collected for this study. A total of 80,312 participants were initially included, with 63,771 excluded due to missing PFAS level information. Of the remaining 16,541 participants with PFAS measurements, 4,761 were further excluded for missing covariates. Fortunately, no outcome variables were missing. Ultimately, 11,780 participants were included in the final analysis, all of whom had complete information on both PFAS measurements and covariates. Further details on participant identification can be found in Figure 1.

Figure 1

Inclusion and exclusion path diagram

Definition of diabetes

Diabetes diagnoses were made by medical professionals, and information on diabetes status was collected through a questionnaire. Participants aged 20 years or older were asked if they had ever been diagnosed with diabetes by a doctor or healthcare professional. Positive responses to this question were used to identify individuals with diabetes in the survey.

Pollutant measurement

PFASs have become extensively utilized and have been found in various environmental sources, including water and organisms [15–17]. Humans can come into contact with PFASs through inhalation, skin exposure, and consuming contaminated food. Once in the body, PFASs initially circulate in the bloodstream and undergo metabolic changes, leading to the production of metabolites that can accumulate in specific organs and trigger toxic responses, such as endocrine disruption [18, 19]. In this research, blood samples from participants were analyzed for five specific PFASs: perfluorodecanoate (PFDeA), perfluorohexane sulfonate (PFHxS), 2-(N-methyl-perfluorooctane) sulfonamido acetate (Me-PFOSA-AcOH), perfluorononanoate (PFNA), and perfluoroundecanoate (PFUA). These PFASs were identified as persisting in the environment from 2002 to 2018 using high-performance liquid chromatography coupled with tandem mass spectrometry (MS/MS). When the limit of detection (LOD) for a chemical was below 50%, the LOD logarithm/2 method was used to calculate the corresponding data, as 86% of the chemicals met this criterion. To ensure quality, procedural blanks and spiked samples were included in each batch of samples to monitor potential contamination and analytical performance. Detection limits were set for each PFAS analyte, and measurements below these limits were processed using the NHANES LOD/2 method, a standard practice for handling non-detectable values in environmental exposure assessment.

Covariates

A comprehensive review of the relevant literature was conducted to identify potential covariates associated with exposure to perfluoroalkyl and PFAS and the risk of diabetes. The selected covariates were categorized as either continuous or categorical variables. Continuous variables included age, family poverty index ratio (PIR), annual household income, and BMI. Categorical variables encompassed gender (male, female), race (Mexican American, other Hispanic, non-Hispanic white, non-Hispanic black, other race – including multiracial), education (less than ninth grade, grades 9–11 (including grade 12 without a diploma), high school/GED or equivalent, college or AA degree, associate’s degree or above, denial, don’t know), and marital status (married, widowed, divorced, separated, never married, living with a partner).

Statistical analysis

Descriptive statistics were used to summarize the demographic characteristics of the participants and concentrations of biomarkers. Urine creatinine was not needed as a correction factor since it is a laboratory measure of blood and does not require urine dilution. Categorical variable data were presented as the number of cases (percentage) and underwent χ² testing. Non-normal continuous variable data were presented as the median (interquartile range) [M (Q1, Q3)], and analyzed using the Wilcoxon rank sum test. Normal continuous variable data were presented as mean ± standard deviation, and analyzed using the independent sample t-test. A natural log transformation (ln transformation) was applied to PFAS to establish a normal distribution. Spearman’s correlation was used to assess the relationships between chemicals.

A weighted multivariate logistic regression model was utilized to investigate the correlation between natural logarithm (ln) transformed and quartile transformed chemical concentrations and the incidence of diabetes. In order to assess the relationship between combined chemical exposure and diabetes, as well as to pinpoint the key chemicals influencing the occurrence of diabetes, a WQS regression model was employed. The WQS model is a statistical method that assesses the impact of each component of environmental exposure on the overall effect using a weighted index. This model has been widely used in evaluating the effects of exposure to mixtures [20]. The WQS model calculates quantiles for chemicals and assigns individual weights based on their relative importance within the mixture, facilitating the identification of potentially harmful substances.

The study utilized the quantile g-computation model to investigate the association between chemicals and diabetes by evaluating the ‘mixture effect. Quantile g-computation is a unique method for analyzing combined environmental exposures, as it estimates the parameters of marginal structural models and offers causal effect estimates. This model predicts the anticipated change in potential outcomes if a simultaneous intervention is implemented for all exposures, potentially while considering confounding factors. Accurate estimates of the true effect can be achieved through quantile g calculations when the assumptions of exchangeability, causal consistency, positivity, no interference, and correct model specification are met. Statistical analysis was performed using R3.4.3, with a two-sided p-value below 0.05 considered statistically significant [21].

Note: Due to the high proportion of samples with LOD ≤ 50%, they were not excluded from the study population. Instead, these samples were adjusted by dividing them by 2, and then missing values of independent variables were excluded. The NHANES protocol was approved by the NCHS research ethics review board, and all participants provided informed consent.

Baseline characteristics of the study population

A total of 11,780 participants were included in this study for statistical analysis. The general characteristics of the study population can be found in Table I. Among these participants, 1,452 were diagnosed with diabetes, representing a prevalence of 12.32%. The remaining 10,329 participants were diagnosed as non-diabetic, accounting for 87.68% of the sample. Statistical analysis indicated significant differences in age, race, education level, marital status, family income, BMI, and diabetes status between patients with and without diabetes (p < 0.05). However, there was no significant difference in gender between the two groups (p > 0.05).

Distribution of chemical exposure and its association

The nuclear density distribution was analyzed for all five chemicals, and the results are illustrated in Supplementary Figure S1. Supplementary Figure S1 shows that PFUA displays three distinct peaks, while PFNA, monohydroxy polyfluoroalkyl substances (MPAH), PFHxS, and perfluoroalkyl ether acids (PFAeA) each show a single prominent peak.

A Spearman correlation analysis was conducted to investigate the relationship between the measured values of the five chemicals, as shown in Figure 2. The results indicated that the correlation coefficient between PFUA and PFDeA was the highest at 0.76, suggesting a strong positive correlation. Furthermore, the correlation coefficient between PFUA and MPAH was 0.69, and that between PFDeA and MPAH was 0.66, both showing significant positive correlations. On the other hand, the correlation coefficients for the remaining chemicals were below 0.2, indicating either a very weak correlation or no correlation at all. This suggests the possibility of high collinearity among the variables, necessitating the use of multivariate logistic regression to evaluate and confirm the collinearity.

Figure 2

Spearman correlation coefficient between pairwise chemical measurements

Weighted generalized linear regression

Univariate logistic analysis was conducted using three models, each adjusted for specific covariates. The results are presented in Supplementary Tables SI–SIII. Model 1 included adjusted covariates of gender and age. Model 2 expanded on this by including additional covariates such as race, PIR, education level, and BMI. Lastly, Model 3 incorporated sex, age, race, PIR, education level, BMI, marital status, and annual household income as adjusted covariates.

None of the five chemicals showed a significant odds ratio (OR) for diabetes in all three models (p > 0.05). Specifically, in the adjusted Model 3, the highest quartile (Q4) of PFDeA, PFHxS, and PFNA had values above 1, but none were statistically significant (p > 0.05). However, it is important to note that the lack of significance could be due to the small sample size and imbalance in the distribution of disease and non-disease populations in the study. Therefore, the lack of significance in this study should be interpreted with caution, as it may be influenced by these factors.

In the multivariate weighted logistic regression models, all five chemicals showed a variance inflation factor (VIF) of less than 10, suggesting no multicollinearity among them (Table II). The OR of PFNA on the outcome was calculated to be 1.14 (0.91, 1.43) in the multivariate weighted logistic regression model. Nevertheless, this relationship was not deemed statistically significant (p = 0.25).

Weighted quantiles and regression (WQS)

In the positive constraint model of the WQS, the coefficient for the chemical mixture index showed a weak association with diabetes (OR = 1.02), which was not statistically significant. However, strong correlations between the chemical mixture index coefficient and diabetes were observed in Model 1, Model 2, and Model 3 after adjusting for relevant covariates (Model 1: OR = 1.77, Model 2: OR = 1.83, Model 3: OR = 1.83). The results are presented in Supplementary Table SIV. Similarly, in the unadjusted model of the negative constraint in the WQS, the OR was 1.98, but it was not statistically significant (p > 0.05). However, after adjusting for covariates, a robust association between the chemical mixture index and diabetes was evident in Model 1, Model 2, and Model 3 (Model 1: OR = 1.87, Model 2: OR = 1.90, Model 3: OR = 1.90). These findings are summarized in Supplementary Table SV.

The first model (Model 1) considered two covariates: gender and age. The second model (Model 2) included a wider range of covariates, such as gender, age, race, PIR, education level, and BMI. Finally, the third model (Model 3) integrated additional covariates, including sex, age, race, PIR, education, BMI, marital status, and annual household income.

In the positive constraint model of WQS, the substances MPAH, PFHxS, and PFUA are considered to be relatively significant, with weights exceeding 1/18 each. Similarly, within the negative constraint model of WQS, both PFDeA and PFNA are deemed to be relatively important, with weights surpassing 1/18 for each chemical. For more detailed information, see Figure 3.

Figure 3

Weight of regression indicators on diabetes in the WQS model. A – Forward constraints; B – Negative Constraints

Quantile g-computation “mixture effect” evaluation

For further information on the use of binary results in Q g-comp using the Q g-comp.no boot function, please refer to Supplementary Table SVI. The Q g-comp.boot function produces a marginal OR. It is important to understand that introducing nonexposed covariates, such as confounders, into the model will result in estimations of these parameters that are not equivalent. This is because the OR is non-collapsible. Marginal parameters offer estimates of population-average exposure effects, which are typically more interpretable and meaningful than conditional odds.

Further adjustments to the covariates were made using Model I, Model II, and Model III. The results of these adjustments are illustrated in Figure 4, showcasing four different models (1: no adjustment; 2: Model one; 3: Model two; 4: Model three) and their impact on the overall effect for each exposure. It is important to highlight that the constraints imposed on the weights in the WQS model may lead to biased estimations of the effects. In contrast, the Q g-comp model allows for weights to shift in different directions, indicating the potential for some exposures to be beneficial while others may be harmful. This flexibility also addresses sampling variation that can occur with small or medium-sized samples. Within the Q g-comp model, the weights represent the proportion of effects aligning in the same direction, although in certain instances, they may correspond to a smaller or larger proportion compared to the overall ‘mixed’ effect. Notably, Figure 4 illustrates that the combined exposure effect of PFNA on diabetes is positive in all four models, while exposure to the other four chemicals shows negative effects.

Figure 4

Q G-computation weights for the four models. A – Model I; B – Model II; C – Model III; D – Model 4

Nonlinearities in the data were addressed using Qg-comp and Qg-comp. BOOT methods were applied to develop a model that included an interaction term and quadratic terms for each predictor of chemical exposure. By simultaneously accounting for the nonlinear effects of all exposures, a distinct nonlinear trend in the overall exposure effect was observed, as shown in Supplementary Figures S2 and S3. Importantly, the smoothed regression line aligns with the confidence interval of the marginal linear model, indicating consistency. Additionally, it is clear that the OR for the prevalence of PFAS and diabetes increases steadily with each quantile change.