The NHANES is a series of surveys conducted by the National Center for Health Statistics (NCHS), employing a sophisticated multistage probability cluster design to assess the health and nutritional well-being of a representative sample from the non-institutionalized US population. Comprehensive details are accessible in the NHANES survey methods and analytic guidelines (23). Extracted data originated from the 2017–2018 NHANES cycles, specifically focusing on participants aged 20 years and above. The data underwent analysis between September and December 2023. Ethical approval for NHANES procedures and protocols was obtained from the NCHS Research Ethics Review Board, and all participants provided written informed consent. The study followed the reporting guidelines set forth by the Strengthening the Reporting of Observational Studies in Epidemiology.

The current study included 2,986 participants from the NHANES 2017–2018 cycle who underwent urinary phthalate metabolite measurements. Exclusion criteria were applied post-data collection, removing participants <20 years old (n = 1,221), those with incomplete data on urinary phthalates metabolite and gallstones (n = 65), and individuals lacking demographic information (n = 230), including age, sex, race/ethnicity, marital status, poverty income ratio (PIR), and educational level. Additional exclusions were made for those without data on physical activity (PA), drinking status, smoking status, body mass index (BMI), and urinary creatinine (n = 86). The final dataset for analysis comprised 1,384 participants (Supplementary Figure S1).

In NHANES, the presence of gallstones was evaluated through the question by asking people over the age of 20, “Has a doctor or other health professional ever told you that you had gallstones?.” A response of “yes” to this question was considered indicative of gallstones, while an exact response of “no” was categorized as the absence of gallstones, as previously mentioned in the literature (24). The occurrence of gallstones served as the outcome variable for analysis.

This study investigated 19 urinary phthalate metabolites, as outlined in Supplementary Table S1, including Mono (3-carboxypropyl) phthalate (MCPP), Mono-ethyl phthalate (MEP), Mono-isobutyl phthalate (MiBP), Mono-2-hydroxy-iso-butyl phthalate (MHiBP), Mono-n-butyl phthalate (MBP), Mono-3-hydroxybutyl phthalate (MHBP), Monobenzyl phthalate (MBzP), Mono(2-ethylhexyl) phthalate (MEHP), Mono(2-ethyl-5-carboxypentyl) phthalate (MECPP), Mono(2-ethyl-5-hydroxyhexyl) phthalate (MEHHP), Mono(2-ethyl-5-oxohexyl) phthalate (MEOHP), Mono-isononyl phthalate (MNP), Monocarboxyisooctyl phthalate (MCOP), Mono-oxo-isononyl phthalate (MONP), Monocarboxy-isononyl phthalate (MCNP), Cyclohexane-1,2-dicarboxylic acid mono(hydroxy-isononyl) ester (MHINCH), Cyclohexane-1,2-dicarboxylic acid mono(carboxyoctyl) ester (MCOCH), Mono(2-ethyl-5-hydroxyhexyl) terephthalate (MEHHTP), and Mono(2-ethyl-5-carboxypentyl) terephthalate (MECPTP). Among them, MNP is excluded from the analysis due to over 80.07% of participants having concentrations below the limit of detection (LOD) to mitigate potential bias. The LOD and molecular weight (MW) of each metabolite and its corresponding parent compound are shown in Supplementary Table S1.

The quantification of phthalate metabolites in urine employed a high-performance liquid chromatography-electrospray ionization-tandem mass spectrometry (HPLC-ESI-MS/MS) method (25). This process involves enzymatic deconjugation of glucuronidated phthalate monoesters in urine samples, followed by on-line solid phase extraction (SPE) coupled with reversed phase HPLC-ESI-MS/MS. The use of isotopically-labeled internal standards for the phthalate metabolites improves assay precision. In addition, 4-methyl umbelliferone glucuro- nide is employed to monitor deconjugation efficiency. This method enables the rapid and selective detection of monoester metabolites of commonly used phthalate diesters in human urine, with LOD in the low ng/mL range. Concentrations below the LOD were replaced using the LOD/√2 method (26).

Sum phthalate metabolites were calculated to assess the influence of the parent chemical, including Di(2-ethylhexyl) phthalate (ΣDEHP, expressed as MEHP, MW = 278.2 g/mol), Di-isobutyl phthalate (ΣDiBP, expressed as MiBP, MW = 222.2 g/mol), Di-n-butyl phthalate (ΣDBP, expressed as MBP, MW = 222.2 g/mol), Di-isononyl phthalate (ΣDNP, expressed as MONP, MW = 306.4 g/mol), 1,2-Cyclohexane dicarboxylic acid diisononyl ester (ΣDINCH, expressed as MHINCH, MW = 314.4 g/mol), and Di(2-ethylhexyl) terephthalate (ΣDEHTP, expressed as MEHP, MW = 294.3 g/mol). High MW (> 250) phthalate (high-MWP) is the molar sum of MCPP, MBzP, MHiBP, MEHP, MECPP, MEHHP, MEOHP, MCOP, MONP, MCNP, MHINCH, MCOCH, MEHHTP, and MECPTP (expressed as MCPP, MW = 252.2). Low MW (< 250) phthalate (low-MWP) is the molar sum of MEP, MiBP, MHiBP, MBP, and MHBP (expressed as MEP, MW = 194.2) (27).

Several potential confounding variables were considered based on published research and clinical judgment, including age, sex, marital status, race/ethnicity, education level, PIR, BMI, PA, smoking status, drinking status, and creatinine (28). Age, PIR, BMI, PA, and creatinine were included in the adjustment model as continuous variables. In the subgroup analysis, age was categorized as <40, 40–59, and ≥ 60 when utilized as an exposure variable or considered in subgroup analyses. Sex was categorized into female and male. Marital status was grouped into married, never married, living with a partner, and other (e.g., widowed, divorced, or separated). Self-reported race/ethnicity was categorized as Mexican American, non-Hispanic White, non-Hispanic Black, other Hispanic, and others (including multi-racial participants). Educational level was classified as less than high school, high school or equivalent, and above high school. BMI was calculated as the weight (in kilograms) divided by the square of the height (in meters). PA was defined as the time individuals reported spending during the week on activities such as walking or biking, tasks around the home or yard, work, and recreational activity (29). Smoking status was characterized as never (<100 cigarettes smoked in life), former (>100 cigarettes smoked but currently quit), and now (>100 cigarettes smoked and currently smoking) (30). Drinking status was categorized as never (<12 drinks lifetime), former (<12 drinks lifetime or none in the past year), mild (≤1 drink per day for females and ≤ 2 drinks per day for males in the past year), moderate (≤2 drinks per day for females and ≤ 3 drinks per day for males in the past year), and heavy (≥3 drinks per day for females and ≥ 4 drinks per day for males in the past year) (31). Creatinine was measured from the laboratory.

Complex sampling design and environmental subsample weights were considered in our analyses. Participant characteristics were calculated based on the presence or absence of gallstone, with categorical variables presented as numbers and percentages, and continuous variables as means and standard error (SE). Differences across groups were analyzed using the chi-squared test with Rao & Scott’s second-order correction (for categorical variables) and the Wilcoxon rank-sum test (for continuous variables).

The concentrations of individual and sums of phthalate metabolite were naturally log-transformed to achieve a normal distribution, and they were divided into three groups (Tertile1-Tertile3) according to percentile. Survey-weighted multivariable logistic regression was performed to assess associations of individuals and sums phthalate metabolite (treated as both a categorical and continuous variable) with gallstone. The crude model did not account for covariates, while the adjusted model accounted for age, sex, race/ethnicity, PIR, marital status, education level, BMI, PA, smoking status, drinking status, and creatinine. Possible nonlinear effects were modeled using restricted cubic spline (RCS) models with 3 knots at 10, 50, and 90%.

Further to calculating the sum of phthalate metabolites based on MW, we employed both weighted quantile sum (WQS) regression and Bayesian kernel machine regression (BKMR) to analyze the joint effects of phthalate metabolites. The WQS regression model built a weighted index to estimate the combined effects of all predictive factors related to the outcomes, demonstrating high sensitivity and specificity in identifying significant exposures (32). Bootstrapping with 1,000 iterations was employed to construct the WQS index as a representative measure of mixed exposure levels to phthalates, and the contribution weights of each compound to the overall impact of the mixture were calculated. The dataset was randomly divided, with 40% allocated to the training set and the remaining 60% used as the validation set. The WQS model estimated the incremental risk of gallstones corresponding to a 25% (one-quarter) increase in the WQS index. BKMR is an approach that models the health effects of complex chemical mixtures, using flexible functions and Bayesian methods to identify important mixture components and account for correlations, particularly in high-dimensional settings (33). After adjusting for all covariates, this model underwent 10,000 iterations using the Markov Chain Monte Carlo algorithm. Posterior incorporation probabilities (PIPs), spanning from 0 to 1, were derived from the BKMR model to ascertain the relative significance of each exposure concerning the outcome, with a defined threshold of 0.5.

In additionally analysis, we further explored potential mediators of the association between phthalate exposure and gallstone, including oxidative stress markers (including bilirubin, uric acid, and gamma glutamyl transferase) (34), inflammatory markers (including C-reactive protein, alkaline phosphatase, neutrophils, and ferritin) (35), metabolic syndrome (36), body composition [including BMI, waist circumference, visceral adiposity index (37), lean/fat mass measured using dual-energy X-ray absorptiometry (38)], diabetes and insulin resistance (including fasting blood glucose [FBG], fasting serum insulin [FINS], triglyceride glucose index, and insulin resistance [HOMA-IR], insulin sensitivity [HOMA-IS], and β-cell function [HOMA-β] assessed through the homeostasis model assessment [HOMA]) (39, 40). Mediation analyses were conducted by using the Sobel test, Bootstrap, and the quasi-Bayesian Monte Carlo method with 1,000 simulations based on normal approximation (41, 42). Stratified analyses were performed based on age (<40, 40–59, or ≥ 60 years), sex (female or male), race/ethnicity (Mexican American, non-Hispanic Black, non-Hispanic White, other Hispanic, or other), and BMI (<25, 25–29, ≥30). The interaction was assessed by adding a cross-product term to the regression model and calculating the p value using the likelihood ratio test, comparing models with and without the interaction term.

All analyses were conducted with R (4.2.3) and Free Software Foundation statistics software (version 1.9.2). Statistical significance was determined by two-sided p values below 0.05.

Table 1 presents the characteristics of a sample representing 200.6 million U.S. adults for analysis. The weighted mean age was 47.45 years, with females comprising 50.94%. Among this population, 21.20 million had gallstones, more prevalent among those aged ≥60 years (weighted percentage, 50.45%) and women (weighted percentage, 75.76%). Compared to those without gallstones, individuals with gallstones had higher BMI, HOMA-IR, HOMA-β, fasting insulin, MECPP, MEOHP, MCOP, and ΣDEHP levels, alongside lower HOMA-IS and MCPP concentrations.

Figure 1 illustrated the association between individual phthalate metabolites and gallstones. In the adjusted multivariable logistic regression, each one-unit increase in the natural log-transformed concentrations of MECPP, MHINCH, and MCOCH was associated with a higher incidence of gallstones by 37, 41, and 54%, respectively. Individuals in the third tertile exhibited a lower incidence of gallstones compared to those with MHiBP concentrations in the first tertile. Higher concentrations of MEOHP, MBzP, MCOP, and MCPP were associated with a higher incidence of gallstones in the crude model, but this association was no longer observed after accounting for confounding factors.

Figure 2 depicted the association between the sums of phthalate metabolites and gallstone occurrence. Following adjustments for confounding factors, logistic regression indicated a 31% increase in gallstone risk per unit rise in ΣDEHP and an 87% elevation for ΣDINCH. WQS regression revealed a 91% escalation in gallstone risk per interquartile range increase for High-MWP, 87% for ΣDINCH, and 57% for ΣDEHTP mixtures. The major contributors in High-MWP were MCOCH, MECPP, and MBzP, while in ΣDINCH, it was MCOCH, and in ΣDEHTP, it was MECPTP (Supplementary Figure S2). No significant associations were observed for Low-MWP, DBP, DiBP, and DNP with gallstones.

RCS regression indicated a significant dose–response association with gallstones for the individual phthalate metabolites, such as MECPP, MEOHP, MHINCH, MCOCH, and MEHHTP, as well as for the sum phthalate metabolites, such as ΣDEHP and ΣDINCH (Supplementary Figures S3, S4; Figure 3). In BKMR analysis, the overall effect of mixtures was associated with a lower risk of gallstone when all DEHP metabolites were below the 25th percentile compared to their median (Figure 3), with the PIP for both MECPP and MEHHP exceeding the 0.5 threshold.

Given the significance of ΣDEHP and ΣDINCH in the adjusted multivariate logistic regression and RCS, our study has centered on exploring their potential mediating roles in relation to gallstones. Our study reveals that ΣDEHP is linked to higher levels of HOMA-IS, HOMA-β, and FINS. Notably, lower HOMA-IS and higher HOMA-β and FINS levels are significantly associated with a higher risk of gallstones (Figure 4). Multiple testing methods consistently highlight the significant mediating effects of these three indicators, suggesting that ΣDEHP may elevate FINS levels by reducing HOMA-IS and increasing HOMA-β, ultimately contributing to the occurrence of gallstones. Furthermore, we observed a significant mediating effect of HOMA-IS in the association between ΣDINCH and gallstones (Figure 4). However, no notable mediating effects were observed for indicators related to oxidative stress, inflammation, body composition, or metabolic syndrome.

No significant interactions were observed, indicating a consistent relationship between phthalate metabolites and gallstones across these subgroups (Supplementary Tables S2, S3).