Urinary Bisphenol Mixtures at Population-Exposure Levels Are Associated with Diabetes Prevalence: Evidence from Advanced Mixture Modeling

1. Introduction

The contemporary era, often referred to as the “Anthropocene” or the “Plastic Age”, owes its name to the extremely close relationship that humans have developed with plastic polymers [1]. Plastics have become a fundamental element of modern industry because of their countless applications and cost-effectiveness, but they pose a major threat to the environment and to global population health. The centrality of plastic in modern society extends from its use in food packaging to clothing, exposing us directly to xenobiotic compounds [2,3,4]. At the core of this plastic surge is bisphenol A (BPA), a monomer and endocrine disruptor found in epoxy resins and polycarbonates, enhancing the physical properties of diverse polymers [5].

Over time, plastic production has surged dramatically, reaching 368 million tons in 2019 [6], with projections indicating a doubling within the next two decades [7]. This surge, coupled with challenges in waste management and the perpetual cycle of polymer production and recycling, raises concerns about potential far-reaching impacts across various trophic levels [8], posing a substantial threat of chronic human exposure. Recent research unequivocally confirms the widespread presence of plastic monomers, including the emerging bisphenol S and F monomers (BPS and BPF), in water bodies, soil and atmosphere worldwide. [4,9,10].

Parallel to the increased production and use of this class of monomers has been a significant increase in the incidence of diabetes worldwide. Diabetes mellitus (DM) emerges as a prominent global health crisis in the 21st century [11], with a substantial increase in prevalence over recent decades. Back in 1980, approximately 108 million adults aged 20–79 years were affected. Presently, DM affects 10.5% of the global population (536 million individuals), and projections anticipate a potential rise to 12.2% by 2045 (783.2 million) [12]. The risk factors associated with DM encompass a broad spectrum of environmental and genetic elements, including variables such as age, weight, diet, and smoking [13]. This prompts a consideration of the plausible role that environmental pollutants might play in the development or progression of the disease. Existing literature offers insights into this intriguing intersection. Studies suggest a conceivable relationship between DM and environmental pollutants, providing a nuanced perspective on the complex factors contributing to this escalating health challenge [14].

In this regard, a recent meta-analysis demonstrated the need for further analysis of new bisphenols used in industry as substitutes for BPA in view of the increasingly restrictive regulations being promoted in different countries [15]. The results of the meta-analysis showed a statistically significant positive association only with BPS. Nevertheless, it is evident that in the real world, outside the regulated laboratory setting, co-exposure to bisphenols does occur, as demonstrated in human cohorts that quantify the presence of a mixture of bisphenols in urine [16]. However, the paucity of publications studying the possible combined effect of these monomers in the context of DM is striking. In other variations of DM such as gestational DM (GDM), recent research shows that BPS and BPF could be potential risk factors for its development [17,18].

Consequently, the present work aims to analyze the possible combined effect of the mixture of the three most relevant bisphenols in modern industry with the risk of DM in one of the largest global cohorts analyzing urinary bisphenols, the “National Health and Nutrition Examination Survey” (NHANES) cohort. Thus, using correlations, logistic regressions, Weighted quartile sum and quantile G-computational analysis, the possible relationship between combined bisphenol exposure and DM in adults will be explored.

2. Materials and Methods

2.1. Data ExtRaction fRom the NHANES CohoRt

The NHANES cohort databases from all years in which urinary bisphenol levels were quantified were utilized. Currently, the 2013-2014 and 2015-2016 cohorts have data on the quantification of urinary BPA, BPF and BPS. All data are accessible on the official website of the Centers for Disease Control and Prevention [19] (accessed on november 13, 2023). After consolidating and organizing the data by patient code, a total of 20,146 subjects were obtained. The data were then filtered based on subjects aged 18 and above (adults), urinary bisphenol values, and urinary creatinine values, resulting in a total of 3,699 subjects included in the database for analysis. The study population was further subdivided based on the pathology of interest (diabetes), and participants with data on all relevant covariates (to be detailed later) were selected, yielding 3,658 study subjects (2,965 non-diabetics and 693 diabetics) (Figure 1). Individuals with physician-diagnosed diabetes, as well as those with fasting glucose values ≥ 126 mg/dl or haemoglobin A1c ≥ 6.5%, were included in the diabetic group [20].

2.2. Combined URinaRy Bisphenol ExposuRe Study

In the combined study of urinary bisphenols, we first used an approach based on the analysis of the sum of these phenolic compounds. Initially, analyses were conducted based on “Unsupervised summary scores or USS,” as outlined by Li et al. [21]. This involved calculating the molar sum of bisphenols (ΣBPs/Mol) by dividing the concentration of each metabolite by its respective molecular weight and subsequently summing them up:

$$\mathsf{\Sigma }\mathrm{B}\mathrm{P}\mathrm{s}/\mathrm{M}\mathrm{o}\mathrm{l}=({\displaystyle \frac{\left[\mathrm{B}\mathrm{P}\mathrm{A},\frac{\mathrm{n}\mathrm{g}}{\mathrm{m}\mathrm{l}}\right]}{228.29\mathrm{g}\ast \mathrm{m}\mathrm{o}\mathrm{l}}}+{\displaystyle \frac{\left[\mathrm{B}\mathrm{P}\mathrm{S},\frac{\mathrm{n}\mathrm{g}}{\mathrm{m}\mathrm{l}}\right]}{250.27\mathrm{g}\ast \mathrm{m}\mathrm{o}\mathrm{l}}}+{\displaystyle \frac{\left[\mathrm{B}\mathrm{P}\mathrm{F},\frac{\mathrm{n}\mathrm{g}}{\mathrm{m}\mathrm{l}}\right]}{200.23\mathrm{g}\ast \mathrm{m}\mathrm{o}\mathrm{l}}})$$

In the subsequent step, the total concentration of bisphenols, adjusted by the urinary creatinine value (Creat), was employed to mitigate potential result discrepancies stemming from variations in the glomerular filtration capacity among individual study subjects (ΣBPs/Creat).

$${\displaystyle \frac{\mathsf{\Sigma }\mathrm{B}\mathrm{P}\mathrm{s}}{\mathrm{C}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{t}}}(\mathrm{\mu }\mathrm{g}/\mathrm{g})=\left({\displaystyle \frac{\left[\mathrm{B}\mathrm{P}\mathrm{A},\frac{\mathrm{\mu }\mathrm{g}}{\mathrm{L}}\right]}{\mathrm{C}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{t},\frac{\mathrm{m}\mathrm{g}}{\mathrm{d}\mathrm{L}}}}+{\displaystyle \frac{\left[\mathrm{B}\mathrm{P}\mathrm{S},\frac{\mathrm{\mu }\mathrm{g}}{\mathrm{L}}\right]}{\mathrm{C}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{t},\frac{\mathrm{m}\mathrm{g}}{\mathrm{d}\mathrm{L}}}}+{\displaystyle \frac{\left[\mathrm{B}\mathrm{P}\mathrm{F},\frac{\mathrm{\mu }\mathrm{g}}{\mathrm{L}}\right]}{\mathrm{C}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{t},\frac{\mathrm{m}\mathrm{g}}{\mathrm{d}\mathrm{L}}}}\right)\mathrm{x}100$$

The quantification of bisphenols was conducted using an online solid-phase extraction coupled with high-performance liquid chromatography and tandem mass spectrometry (online SPE-HPLC-Isotope dilution-MS/MS) method. The respective detection limits were 0.2, 0.2, and 0.1 for BPA, BPF, and BPS. For more information on the protocols, refer to the NHANES website [22,23].

Basic descriptive statistics were conducted to analyze the variables of interest within each of the pathological subgroups, as well as in the bisphenol’s quartiles. Quantitative variables were reported as the arithmetic mean (standard deviation), AM (SD), or as the geometric mean (95% confidence interval), GM (95% CI), depending on their normality. Descriptive statistical analysis, comparative analysis, and graphical representation were performed using GraphPad Prism 7.0 software (GraphPad Software Inc., San Diego, CA, USA). The distribution of data was assessed using the D’Agostino-Pearson, Shapiro-Wilk, and Kolmogorov-Smirnov normality tests. Subsequently, either the T-student or Mann-Whitney test was employed for the comparative analysis of two variables. For three or more variables, one-way ANOVA or Kruskal-Wallis, followed by a Bonferroni or Dunn’s test, was carried out. T-test and Fisher’s exact test was used for the analysis of dichotomous variables in descriptive statistics. The p-values in the figures and tables correspond to the post hoc test, with a significance level set at p < 0.05.

Next, Pearson correlation analysis was employed to assess potential relationships between urinary bisphenols and quantitative variables such as age, body mass index, cotinine levels, serum glucose, and glycosylated haemoglobin. In the subsequent step, a study model involving binomial and multinomial logistic regression was implemented. The aim of the logistic regression analysis was to identify patterns that could reveal potential associations between the pathological subgroups and the composite of urinary bisphenols, considering the multifactorial origin of the study pathologies. Subsequently, each analysis was conducted in three distinct manners: first, individually (1); second, adjusted for age, gender, and Body Mass Index (BMI) (2); and third, adjusted for (2) + race/ethnicity, poverty-income ratio, hypertension, dyslipidemia, and smoking (3) (inclusion of anthropometric and demographic variables as per relevant literature [21,24,25]). Individuals classified as hypertensive comprised those diagnosed by a healthcare professional and those with systolic pressure ≥ 140 mmHg or diastolic pressure ≥ 90 mmHg. The category of patients with dyslipidemia included those diagnosed with cholesterol disorders or those with fasting total cholesterol levels ≥ 240 mg/dL[26]. Regarding smoking, individuals were included if they answered affirmatively to the question “have you smoked more than 100 cigarettes in your life?” or if they had serum cotinine values exceeding 10 mg/dL [27]. Log-transformation of urinary concentrations of bisphenol metabolites was applied to achieve normalized distributions. The Cochran q test and regression models were performed using IBM SPSS Statistics for Windows software, version 27 (IBM Corp, Armonk, NY, USA).

Following this, the weighted quantile sum (WQS) was constructed utilizing the R package “gWQS” [28] to investigate associations between the exposures and the outcome, concurrently providing a comprehensive summary of the intricate exposure to the particular mixture of interest [29,30]. In essence, WQS regression condenses the overall exposure to the mixture by calculating a singular weighted index, considering the individual contribution of each component through assigned weights [31]. In the combined analysis, the individual value of each of the bisphenols corrected by the urinary creatinine value (BPA/creat, BPS/creat and BPF/creat) was entered.

Subsequently, a quantile G-computation model was executed using the R package “qgcomp” [32] to estimate the collective impact of the bisphenols mixture on Heart Disease. The QG-comp method evaluates the effect of the bisphenol mixture on the outcome by measuring the change in outcome for each quantile increase in the concentration of all bisphenols in the mixture. Furthermore, it calculates the relative contribution of each bisphenol to the overall effect, discerning whether it has a positive or negative direction. The estimated overall mixture effect, denoted as ψ, signifies the alteration in the outcome associated with a quantile increase in the concentration of all bisphenols in the mixture. In essence, ψ captures the joint effect of all bisphenols in the mixture on the outcome, allowing for an assessment of the overall impact of the bisphenol mixture on the outcome of interest.

3. Results

The initial stage in the statistical analysis of subpopulations within the NHANES cohort involved conducting descriptive statistics. As illustrated in Table 1, most of the study covariates employed to adjust the statistical models exhibited noteworthy differences between individuals with and without the pathology of interest. Subjects in the diabetic category were older and had a higher BMI, as is logical in this pathology. In addition, the pathological subgroup presented a higher percentage of subjects with hypertension, dyslipidemia and smokers. With respect to urinary bisphenol levels, there were significant differences between the creatinine-corrected values, while the USS showed no differences between groups. Finally, ethnic differences were also observed, with an increase in the percentage of Mexican American, other Hispanic and non-Hispanic black subjects among healthy and diabetic subjects, as well as a lower percentage of non-Hispanic white diabetics compared to the non-diabetic group.

The subsequent comparative analysis of the BPmix quartiles (Table 2 and Table 3) revealed intriguing ethnic differences, particularly pronounced in values normalized by molecular weight (ΣBPs/Mol). Notably, a negative trend was observed in the ‘Non-Hispanic White’ population, contrasting with a reverse trend in the ‘Non-Hispanic Black’ population. Interestingly, the analysis did not reveal significant results regarding diabetes in either of the two approaches used with corrected bisphenol values. However, ΣBPs/Creat exhibited a significant increase in hypertensive and smoking subjects, and, intriguingly, a higher percentage of women in Q4 compared to Q1. Furthermore, in both approaches, a significant decrease in the Poverty-Income Ratio is observed, and in the case of ΣBPs/Mol, there is an increase in BMI.

The next step involved Pearson correlation analysis between quantitative variables of interest, such as parameters related to diabetes (Glucose, glycated haemoglobin), as well as cholesterol and cotinine. In this case, the results showed a significant and positive relationship with HbA1c but not with glucose, both in the creatinine-corrected and molecular weight-corrected models. It is worth noting that, as indicated in Table 4, the exact number of patients with quantitative values for each of the variables of interest were 3,530 for cotinine, 3,521 for glucose, and 3,546 for HbA1c.

In this case, the results were only statistically significant with the ΣBPs/creat approach. As depicted in Table 5, both individually and corrected for all relevant covariates, including ethnic group (as recommended by NHANES), significant findings were observed. The results suggest that with each one-unit increase in the log-transformed urinary bisphenol mixture corrected for urinary creatinine—accounting for variations in individual glomerular filtration capacity—there is a diabetes risk of 1.103 (1.002-1.214). This risk remains independent of other covariates associated with the condition, including demographic and pathological variables. This result provides new evidence justifying further exploration of potential relationships between the combined action of bisphenols and the risk of diabetes in adults. However, in the case of multinomial logistic regression, no significant results were observed between combined bisphenol quartiles and the risk of diabetes in any of the approaches used in the statistical analysis.

Due to the relatively limited significant results obtained using more conventional statistical methodologies, Weighted Quantile Sum (WQS) and quantile G-computation analyses were employed. In this case, three independent analyses were conducted, utilizing different variables of interest. Firstly, the diabetes parameter (n=3,658) was examined, incorporating the same covariates used in the regression models (age, gender, BMI, ethnicity, poverty-income ratio, hypertension, dyslipidemia, and smoking). In this instance, as depicted in Figure 2A,B, the WQS analysis unveiled that BPS carries significant weight in the diabetes-related effect. Quantile G-computation analysis demonstrates that all factors have a positive impact, yet BPS continues to hold a greater weight, leading to a clear positive association with diabetes when combined, as illustrated in Figure 2D.

Subsequent statistical analyses were conducted on the quantitative variables glucose and Hb1Ac. As evident in Figure 3, in the case of HbA1c, it appears that BPA carries a negative weight, but the predominant influence still lies with BPS, consistently demonstrating a positive association in all analyses with all study variables. In fact, across all results, it is consistently observed that BPS exerts the greatest influence, and the combined effect is consistently positive (Figure 2D, Figure 3D and Figure 4D).

4. Discussion

This manuscript, for the first time, analyzes the potential impact of the three most used bisphenols in the plastic industry on the risk of diabetes in an adult human population cohort. The traditional approach to studying xenobiotic compounds has generally been developed under the premise of “one compound - one disease.” Clearly, in the context of studying multifactorial conditions like diabetes, influenced by lifestyle habits, individual genetics, and even exposure to pollutants, the fundamental axiom used to formulate the methodology must be reconfigured. Therefore, the combined approach to the bisphenol mixture introduces novelty and realism to the findings of this manuscript, which also focuses on a cohort of the general population, the NHANES cohort.

The initial hypothesis, drawn from previous research [33,34,35], suggests a potential additive (or synergistic) influence of BPS and BPF on BPA in the context of diabetes. Considering the substantial structural homology and shared hormonal activity among these three phenolic molecules (evidence of comparable hormonal activity among BPA, BPF, and BPS exists [33]), it is justifiable to treat the urinary bisphenol pool as a unified quantitative variable.

Descriptive statistics have revealed certain differences depending on the approach used to analyze the urinary bisphenol pool. Nevertheless, the most noteworthy findings are the common convergences, specifically the Poverty-Income Ratio. In both models, a significant decrease in the ratio is observed with higher bisphenol exposure (Q4 vs Q1 is significant in both models). This statement suggests that higher bisphenol exposure may be a secondary consequence of socio-economic problems, associated with greater difficulty in obtaining fresh foods [36], prioritizing packaged foods with poorer nutritional quality, as well as limited information about healthy lifestyle habits. This, in turn, has an impact on the healthcare system by increasing non-communicable diseases and the associated morbidity related to conditions such as diabetes mellitus [37,38]. This issue is particularly pertinent in low-income communities where educational resources may be limited. Additionally, heightened consumption of canned foods, frequently lined with bisphenol-containing resins, can contribute to an elevated bisphenol exposure [39]. Similarly, the preparation and heating of food in plastic containers can lead to the leaching of bisphenols into the food, further increasing exposure [40,41].

The use of creatinine normalization is essential for correcting or adjusting the urinary bisphenol value, as standard glomerular filtration estimates may vary among individuals. This is particularly crucial in diabetic patients who frequently experience renal damage, leading to fluctuations in glomerular filtration. Incorporating urinary creatinine is indispensable to ensure clinical relevance in interpreting the urinary bisphenol levels. The results from the binomial logistic regression model also demonstrated a significant and positive odds ratio with the risk of diabetes. However, in the case of the multinomial model, the findings were inconclusive. Hence, more complex statistical analysis models, utilizing WQS and g-computational models, demonstrated that the combined effect is positive not only on diabetes but also on glycemia or HbA1c.

Therefore, the body of evidence in this manuscript suggests the existence of a relationship between exposure to the mixture of bisphenols and the risk of diabetes in the general adult population. At the individual level, the analysis revealed that the predominant impact in the combined effect consistently rests on BPS. This fact is consistent with the evidence described in the manuscript by Moreno-Gómez-Toledano et al. [20], where they analyzed the individual relationship of each of the bisphenols with diabetes, observing a positive and significant relationship only in the case of BPS. Furthermore, a subsequent meta-analysis demonstrated that BPS is the main candidate to consider in future studies on diabetes. In the meta-analysis study, the results showed a combined odds ratio of 1.35 (1.08-1.70), emphasizing the need to delve deeper into the paradigm [15].

The results align with the statistical analysis model conducted by Tang et al. [42] in gestational diabetes. Despite observing a collectively negative impact of BPS, BPA, BPF, TBBPA, and BPB, the most prominent positive effect individually was found with BPS in the quantile-based g-computation model. Three potential factors for comparison arise: First, differences between general population diabetes mellitus and gestational diabetes; second, variations in the analyzed monomers, limited in this study to BPA, BPS, and BPF; and third, the adjustment for creatinine levels could significantly influence the approach to studying the exposure to xenobiotic compounds excreted in urine.

Furthermore, a longitudinal study conducted in 2022 found a positive correlation between the use of BPS and the risk of developing GDM during the first trimester of pregnancy, so perhaps the emerging use of BPS as a substitute for BPA does not guarantee safety and more research is needed to raise awareness of the possible adverse health effects. [43]. Despite the lack of scientific evidence, certain research studies confirm that these alternatives to BPA could affect by causing other types of disorders such as the affectation during embryonic development [44], may be related to other pathologies such as obesity [45] and even interfere in the development of ADHD in childhood [46].

To understand the potential molecular aspects involved in the pathogenesis of the disease, there is evidence developed in experimental cellular and animal models of zebrafish, mouse, and rat. The evidence suggests that BPS may involve the estrogen receptor beta (ERβ), altering key cellular events in pancreatic cells. Marroqui et al. [47] utilized pancreatic β-cells from both wild-type (WT) and ERβ knockout (BERKO) C57BL/6J mice. Following a 48-hour treatment with BPS, they observed enhanced insulin release, reduced ATP-sensitive K+ (KATP) channel activity, and a decreased expression of several ion channel subunits exclusively in β-cells from WT mice, not in those from BERKO mice. This suggests the involvement of an extranuclear-initiated pathway linked to ERβ.

For all these reasons, it is crucial that the relevant institutional authorities in each country take into consideration evidence such as that observed in this manuscript when implementing intervention measures to ensure the safety of citizens. While future studies are needed to conclusively establish the existing causal relationship between the combination of phenolic monomers and diabetes mellitus, the growing body of evidence on the potential deleterious effects of BPA, along with the increasing number of findings related to its substitute molecules in recent decades, underscores the urgency of promoting the application of the precautionary principle. In fact, BPA and other substitutes might remain in sediments of water [48,49] that could pose a risk to humans and wild flora and fauna. Through water flows can reach terrestrial soil, altering ecological processes and reaching agricultural lands[50] that could enter the food chain. There are promising studies focusing on the application of bioremediation to degrade microplastics [51,52] although there is still research to be done in terms of recycling these particles in a large scale. Health policies should be developed with special attention to population groups particularly susceptible to endocrine disruptors but extrapolated to the general population. We must be aware that we live in the “Plastic Age”[53] and act accordingly if we want future generations to lead a life free from endocrine disruptors.

5. Conclusions

The evidence shows the existence of a statistical relationship between the combination of bisphenols in the general population and diabetes mellitus. The results also suggest that BPS could be the monomer with the greatest weight in the context of metabolic disease, indicating that the replacement of BPA by the new BPS and BPF molecules could pose a greater risk to population health, at least in the context of diabetes. Therefore, based on the results obtained and scientific evidence of the last decades regarding BPA, the need to apply the precautionary principle to this class of emerging molecules is evident, since the health of future generations depends on the actions we take today.