Full length article

# Using mixture and exposome methods to assess the associations and interactions of chemical biomarkers with semen parameters

<figure><img src="image_1.png" /></figure>

Maximilien Génard-Walton<sup>a, *, e</sup>, Charline Warembourg<sup>b, e</sup>, Paige L. Williams<sup>c, d, e</sup>, Irene Souter<sup>g</sup>, Jennifer B. Ford<sup>e, f</sup>, Jorge E. Chavarro<sup>c, f, h</sup>, Russ Hauser<sup>c, d, e, i</sup>, Lidia Mínguez-Alarcón<sup>e, h</sup>, for the Earth Study Team

## ARTICLE INFO

Handling Editor: Marti Nadal

## ABSTRACT

**Background:** Although epidemiologic studies have found associations between several chemical biomarkers and semen parameters, only a few studies have evaluated their combined effects and interactions.

**Methods:** We included 420 men aged 18–55 years who participated in the Environment and Reproductive Health (EARTH) study between 2005 and 2019 at the Massachusetts General Hospital Fertility Center. Four mixture/ exposome models, including exposome-wide association studies (ExWAS), Deletion/Substitution/Addition algorithm (DSA), Bayesian Kernel Machine Regression (BKMR), and the non-linear interaction model (NLInteraction), were used to evaluate associations between 24 chemical biomarkers (urinary bisphenol A, four parabens, triclosan, 17 phthalate metabolites, and hair mercury [Hg]) and semen parameters (sperm count, concentration, morphology, motility, and semen volume).

**Results:** Overall, Hg was positively, and triclosan was negatively significantly associated with all semen outcomes except semen volume across multiple models. In BKMR, when increasing the mixture from the 10th to 90th percentile, the chemical mixture was negatively associated with percent normal sperm morphology (−1.42 %, 95 % Credibility Interval [CRI] −12.36; 9.52), motility (−0.10 %, 95 % CRI −1.78; 1.58), and volume (−0.26 mL, 95 % CRI −0.79; 0.28), and positively associated with sperm count (0.21 mill, 95 % CRI −0.23; 0.65) and concentration (0.33 mill/mL, 95 % CRI −0.14; 0.80). We detected one interaction between urinary triclosan and butylparaben in relation to normal sperm morphology using NLInteraction.

**Conclusion:** Among 24 biomarkers tested, Hg and triclosan showed consistent associations with semen parameters while other biomarkers were less consistent depending on the model used. Similar research studies should be encouraged to use multiple models to incorporate interaction assessment and assess the robustness of results.

## 1. Introduction

Poor semen quality, and especially low sperm count, is known to significantly impact male fertility but is also associated with increased all-cause morbidity and mortality (Eisenberg et al., 2014; Latif et al., 2017). Furthermore, studies have observed overall declines in sperm

count and concentration over the last few decades in multiple western countries (Levine et al., 2017; Mínguez-Alarcón et al., 2018a). Exposure to endocrine disrupting chemicals (EDCs) has been brought forward to partly explain these observations of reduced semen quality, with some families of EDCs having received particular attention because their presence in everyday products leads to widespread human exposure.

<footnote><sup>a</sup> Université Grenoble Alpes, Inserm, CNRS, IAB, Site Santé, Allée des Alpes, 38700 La Tronche, France</footnote>

<footnote><sup>b</sup> Univ Rennes, Inserm, EHESP, Irset (Institut de recherche en santé, environnement et travail) - UMR_S 1085, F-35000 Rennes, France</footnote>

<footnote><sup>c</sup> Departments of Epidemiology, Harvard T.H. Chan School of Public Health, Boston, USA</footnote>

<footnote><sup>d</sup> Department of Biostatistics, Harvard T.H. Chan School of Public Health, Boston, USA</footnote>

<footnote><sup>g</sup> Vincent Obstetrics and Gynecology, Massachusetts General Hospital, Harvard Medical School, Boston, MA, USA</footnote>

<footnote><sup>e</sup> Department of Environmental Health, Harvard T.H. Chan School of Public Health, Boston, USA</footnote>

<footnote><sup>f</sup> Department of Nutrition, Harvard T.H. Chan School of Public Health, Boston, USA</footnote>

<footnote><sup>h</sup> Channing Division of Network Medicine, Harvard Medical School & Brigham and Women's Hospital, USA</footnote>

<footnote><sup>i</sup> Department of Obstetrics, Gynecology and Reproductive Biology, Harvard Medical School, Boston, USA</footnote>

<footnote>* Corresponding author at: Institut pour l'Avancée des Biosciences, Site Santé, Allée des Alpes, 38700 La Tronche, France.
E-mail address: maximilien.genard-walton@inserm.fr (M. Génard-Walton).</footnote>
EDCs of interest include chemicals in personal care products (e.g., parabens, phthalates), and food (e.g., metals, phthalates) or food packaging (e.g., bisphenols, phthalates). Recent reviews have highlighted the presence of strong evidence for an association between phthalate metabolites, including those of di(2-ethylhexyl)phthalate (DEHP), and reduced sperm concentration and motility (Wang et al., 2023). Parabens such as methylparaben and butylparaben have been associated with lower sperm count, lower motility, and lower morphology (Smarr et al., 2018), and bisphenol A (BPA) and triclosan with poor semen quality (Daza-Rodríguez et al., 2023; Presunto et al., 2023). Additionally, metals such as lead, and mercury (Hg) have been linked to various semen parameters (Giulioni et al., 2023; Henriques et al., 2019). Among men attending a fertility center, we have previously observed associations between certain phthalate biomarkers as well as Hg (Mínguez-Alarcón et al., 2018a,b), and higher sperm quality (Hauser et al., 2006), but not for parabens (Meeker et al., 2011).

Whilst the literature has reported evidence regarding the associations between select chemicals individually and semen parameters, the literature on their combined effects and interactions is rather limited. This is a critical aspect because studying the impact of correlated chemicals as single exposures without also assessing their mixture effects could lead to spurious associations due to confounding by co-exposure (Weisskopf et al., 2018). The multiplicity of chemicals assessed in a single paper or across several articles from a given cohort can also lead to false positive risk inflation (Braun et al., 2016). In addition, as the number of chemicals assessed increases, the more complex it becomes to study chemical interactions but not evaluating interactions could lead to an underestimation of effects when interaction between chemicals is present (Sun et al., 2013). Multiple methods have been developed in epidemiology to explore associations between an ensemble of chemical exposures and an outcome and avoid issues such as multicollinearity that would arise from including all exposure biomarkers in a simple regression model. However, no single method is able to answer all questions that are raised when working with multiple chemical exposures (Gibson et al., 2019). In the field of male fertility, some studies have previously tried to address this issue. Chung et al. (2019) conducted exposome-wide association studies (ExWAS) to explore associations between 128 chemical biomarkers and semen quality but found no significant associations after false positive correction. Chang et al. (2024) used elastic net and mixture methods including weighted quantile sum regression (WQS) and Bayesian kernel machine regression (BKMR) to investigate the relationship between 21 pesticide biomarkers and semen parameters, and reported associations between most exposure biomarkers as well as the mixture and sperm motility parameters in all models. Lastly, Gao et al., (2024) found, using BKMR and quantile g-computation, significant associations between a mixture of six non-persistent chemical biomarkers and sperm count and motility. However, these studies were mainly focused on assessing an overall mixture effect, rather than identifying interactions, except in Chang et al.'s study (2024) in which BKMR was used to explore the possibility of interactions between chemical biomarkers.

The objective of our analysis was to assess associations of chemical biomarkers with semen parameters in the context of co-exposures and their interactions. However, because no superior method to answer these objectives exists, and in order to assess the robustness of model results, we opted to use methods from both the mixture and exposome fields to assess the relationship between multiple chemical biomarkers and semen parameters among men attending a fertility center who participated in the EARTH cohort.

# 2. Methods

## 2.1. Population

Participants were recruited into the EARTH Study, a prospective cohort study aimed at assessing the impact of environmental and dietary

factors on reproductive health in both women and men at the Massachusetts General Hospital (MGH) Fertility Center (Boston, United States of America). Full details for the inclusion of participants and their follow-ups have previously been reported (Messerlian et al., 2018). In brief, both partners of couples seeking fertility care at MGH were invited to participate in the study – which they could do either as couples or individually – and followed until they were pregnant through birth (for those achieving pregnancy), or until they discontinued treatment or withdrew from the study. Between 2005 and 2019, 511 men aged 18–55 years who had not had a vasectomy and who were not taking hormones at the time of recruitment, were enrolled and provided at least one semen sample. At enrolment, research staff administered sociodemographic, lifestyle, and medical history questionnaires to participants. Weight and height were measured by trained study staff. Study participants also completed comprehensive questionnaires on family, medical, reproductive and occupational history, perceived stress, consumer products use, smoking history, diet, and physical activity.

At enrolment and for each IVF or IUI cycle (on the day their female partner underwent their scheduled fertility treatment procedure), male participants provided a spot urine sample, non-fasting blood sample and semen sample along with the abstinence time questionnaire. Hair samples were collected after recruitment into the study by instructing participants to collect hair when going for a haircut to the barber/ hairdresser (excluding men with very short hair or bald men).

The study was approved by the Human Subject Committees of the Harvard T.H. Chan School of Public Health (HSPH), the Institutional Review Boards of MGH, and the Centers for Disease Control and Prevention (CDC). Informed consent was obtained from all participants.

## 2.2. Semen assessment

Semen was collected on site at the MGH via masturbation in a sterile plastic specimen cup following a recommended 48-hour abstinence period. Some men (N = 235) provided multiple semen samples for further evaluation of their semen quality and/or because their female partners underwent more than one cycle of fertility treatment. Before analysis, the sample was liquefied at 37 °C for 20 min after collection and analyses were conducted 10–20 min afterwards. Semen volume (ml) was measured by an andrologist using a graduated serological pipet. Sperm concentration (million/ml) and motility (% motile) were assessed using a computer-aided semen analyzer (10HTM-IVOS, Hamilton-Thorne Research, Beverly, MA, USA). To measure sperm concentration and motility, 5 ml of semen from each sample was placed into a pre-warmed (37 °C) Makler counting chamber (Sefi-Medical Instruments, Haifa, Israel). A minimum of 200 sperm cells from at least four different fields were analyzed from each specimen. Motile sperm were defined according to the WHO classification (World Health Organization, 2010). Total sperm count (million/ejaculate) was calculated by multiplying sperm concentration by semen volume. Sperm morphology (% normal) was measured on two slides per specimen (with a minimum of 200 cells assessed per slide) with a microscope using an oil-immersion 100× objective (Nikon, Tokyo, Japan). Strict Kruger scoring criteria were used to classify men as having normal or below normal morphology (Kruger et al., 1988). For quality assurance/quality control purposes, the laboratory conducted weekly monitoring of sperm morphology smears. In addition, the laboratory performed quarterly competence evaluations of all technicians and proficiency testing by an outside evaluator occurred every 6 months. Semen analysis was not available for 45 men, and a further 8 azoospermic men were excluded from this study analysis (N = 458 men).

## 2.3. Chemical exposure assessment

Urine samples were collected in a sterile polypropylene specimen cup. Similar to semen samples, some men (N = 235) provided multiple urine samples because their female partners underwent more than one
cycle of fertility treatment. Specific gravity was measured within hours of collection at the HSPH at room temperature using a handheld refractometer (National Instrument Company, Inc., Baltimore, MD, USA) calibrated with deionized water before each measurement and was used as a covariate for adjustment in the statistical models. The urine was stored at -80 °C and frozen samples were shipped on dry ice overnight to the CDC for analysis. As previously described (Silva et al., 2007; Silva et al., 2013; Zhou et al., 2014), we used online solid-phase extraction coupled with isotope dilution-high-performance liquid chromatography-tandem mass spectrometry to quantify the urinary concentrations of BPA, four parabens (methylparaben, propylparaben, butylparaben, ethylparaben), triclosan, 17 phthalate metabolites: mono-methyl phthalate (MMP), mono-ethyl phthalate (MEP), mono-n-butyl phthalate (MnBP), mono-3-hydroxybutyl phthalate (MHBP), mono-isobutyl phthalate (MiBP), mono-2-methyl-2-hydroxypropyl phthalate (MHiBP), monobenzyl phthalate (MBzP), mono(3-carboxypropyl) phthalate (MCPP), mono(2-ethylhexyl) phthalate (MEHP), mono(2-ethyl-5-oxohexyl) phthalate (MEOHP), mono(2-ethyl-5-hydroxyhexyl) phthalate (MEHHP), mono(2-ethyl-5-carboxypentyl) phthalate (MECPP), mono-isononyl phthalate (MNP), mono-carboxyoctyl phthalate (MCOP), monocarboxy-isononyl phthalate (MCNP), and two metabolites of the phthalate substitute 1,2-cyclohexane-1,2-dicarboxylic acid diisobonyl ester (DINCH): cyclohexane-1,2-dicarboxylic acid-mono (hydroxy-isononyl) ester (MHiNCH), cyclohexane-1,2-dicarboxylic acid-mono(carboxyoctyl) ester (MCOCH). Limits of detection (LOD) of urinary biomarkers ranged from 0.1 to 1.0 µg/L. Along with study samples, each analytical run included a set of calibrators, reagent blanks, and high- and low-concentration quality control (QC) materials. Concentrations of the QCs were evaluated using standard statistical probability rules (Caudill et al., 2008). Precision over three months, calculated from the coefficient of variation of repeated analysis of QC samples, was 2.7 %–14 %, depending on the biomarker.

Men who could not provide a hair sample on the day of semen collection (N = 133) sent their hair sample by mail, loose in an envelope. These men were briefed and given written instructions to collect 2 cm of hair close near the scalp in the occipital area. Upon reception, hair was cleaned before analysis at the Institute for Exposomic Research at the Icahn School of Medicine at Mount Sinai (New York, USA) to remove extraneous contaminants by sonication for 15 min in a 1 % Triton X-100 solution. After sonication, samples were rinsed with distilled deionized water and dried 5 times at 60 °C for 48 h. Total Hg in parts per million (ppm) was measured using 0.02 g of the 2 cm of proximal hair using a Direct Mercury Analyzer 80 (Milestone Inc, Monroe, CT) with a matrix matched calibration curve. Certified reference material GBW 07601 (human hair; Institute of Geophysical and Geochemical Exploration, China) containing 360 ppm Hg was used as the quality control standard. The LOD for Hg was 0.01 ppm with the percentage recovery for quality control standards ranging from 90 to 110 percent. Percentage difference in duplicate sample (same time point) analysis was <10 %.

Overall, chemical biomarker assessment was available for 420 men.

## 2.4. Statistical analyses

The clinical and demographic characteristics of the study population were examined using means ± standard deviations (SDs) for continuous variables, medians and inter-quartile ranges for discrete variables, and counts and proportions for categorical variables. We used Spearman correlation coefficients to summarize the correlations between urinary biomarker concentrations. Chemical biomarkers with concentrations detected in less than 60 % of samples were excluded from the analyses. Concentrations below the LOD for the remaining chemical biomarkers were imputed using multiple imputation through chained equations. Missing data for covariates were also imputed within this same process.

In modelling relationships of chemical biomarkers with semen parameters, all chemical biomarker concentrations were first natural log-transformed. For the purpose of reducing measurement error inherent

to chemical biomarkers with relatively short biological half-lives (Roggeman et al., 2022), and retaining maximum information despite model incompatibility with repeated measures, sample concentrations of both chemical biomarkers and semen were averaged (arithmetic mean) over the whole study period. A sensitivity analysis was also performed using information from the first cycle only.

ExWAS, a method widely used in exposome studies, was used as a first step to obtain single-exposure associations, and allow comparability with other studies (Juarez and Matthews-Juarez, 2018). ExWAS results are reported using volcano plots representing significance versus fold-change of the outcome. In order to best answer the study questions, model choice for the statistical analyses was based on the capacity of the model to detect interactions between chemical biomarkers, and to model the large number of chemical biomarkers included, as well as being easily implementable through the use of available R packages. Models included were the Deletion/Substitution/Addition algorithm (DSA), BKMR, Group Lasso Interaction Network (GLINTERNET), and the non-linear interaction model (NLInteraction).

The DSA model aims to select, through an iterative process, a final model that minimizes the value of the root mean squared error of predictions using cross-validated data (Vrijheid et al., 2020). This is achieved by starting with an empty model and adding, deleting or substituting (replacing selected variables by unselected ones) variables in the pending model until it reaches a final model. This process is then repeated multiple times using a different subset of the data every time as the validation set. Finally, the results from each validation step are averaged to produce a more robust estimate of the associations. We allowed for quadratic relationships and interaction terms with a maximum order of two, and made no restriction on the number of predictors. In order to stabilize the results, we ran 50 DSA models on the same data for each outcome using different seeds, and retained chemical biomarker concentrations that were selected in at least 10 % of the DSA runs (Kusters et al., 2022). These variables were then used in a linear regression model with results reported as mean differences of the outcomes and their corresponding 95 % confidence intervals (CI).

BKMR is a non-parametric Bayesian method which uses kernel functions to handle the high-dimensional parameter space of the exposure-response function in a flexible way that allows for nonlinear and non-additive effects (Bobb et al., 2015). It can provide an overall mixture effect when increasing mixture exposure biomarker concentrations from one percentile to another (thereby reducing the number of dimensions), and can also perform variable selection. Moreover, it can assess the presence of interactions between mixture components by assessing the effect of one exposure biomarker on the outcome whilst one or all others are fixed at various quantiles. Model estimations are summarized over multiple iterations (20,000 in this study) with different sets of parameter values to reduce model uncertainty. Mixture associations are presented as mean differences of the outcomes corresponding to an increase in all urinary concentrations of the chemical biomarker mixture from the 10th to the 90th percentiles, along with their corresponding 95 % credible intervals (CRI). We report (i) the overall association of the mixture corresponding to a simultaneous quantile increase in all mixture components (versus median concentrations), (ii) exposure-response relationships for each individual component while holding all other components at their median levels, (iii) and single-exposure biomarker associations (mean difference comparing the 75th to 25th percentile) when all the other exposure biomarker concentrations are fixed at the 25th, 50th, and 75th percentiles.

NLInteraction is a flexible Bayesian model that allows exposure biomarkers to interact with each other and have nonlinear relationships with the outcome using multivariate spike and slab priors (Antonelli et al., 2020). This model also uses iterations (20,000 in this study) to stabilize results. The interactions are interpretable through the use of posterior probabilities of inclusion (PIPs) into the model to identify components that interact with each other.

We also intended to use a tree-based method, GLINTERNET, a model
based on the overlapped group-lasso method that fits linear pairwise-interactions in its variable selection algorithm (Barrera-Gómez et al., 2017), but the models were not stable enough to provide estimations.

ExWAS consists in applying independent linear regression models and correcting association p-values for multiple testing. The correction used was an adaptation of the Bonferroni procedure that estimates the number of truly independent tests according to the correlation structure (Li et al., 2011) and then divides the nominal significance by these calculated effective number of tests.

When possible (ExWAS, DSA, and BKMR), multiple imputations were used. A set of covariates was selected using a directed acyclic graph (DAG) based on prior knowledge (Supplementary Fig. S1) and included age (continuous), body mass index (BMI, categorized as ≤20, between 20 and 25, and between 25 and 55 kg/m²), and education level (no graduate degree and graduate degree). In addition, all models were adjusted on abstinence time, as well as specific gravity to account for urine dilution. All models assessed the associations between sperm count (mill), concentration (mill/mL), morphology (% normal), motility (% motile), and volume (mL). Sperm count and concentration were natural log-transformed to ensure linear regression assumptions were respected. All analyses were performed using R (version 4.4.1). The threshold for significant results for all our analyses except ExWAS was 5 %.

# 3. Results

## 3.1. Population description

In total, 420 men provided 917 semen (median = 2, interquartile range [IQR] = 2) and urine samples (median = 2, IQR = 2) and one hair sample each. These men were included in the EARTH cohort between

<caption>**Table 1**</caption>

Demographic and reproductive characteristics among 420 men enrolled in the EARTH Study between 2005 and 2019.

<table><thead><tr><th>Characteristic<sup>a</sup></th><th>N = 420</th></tr></thead><tbody><tr><td>Age, mean ± SD</td><td>36.4 ± 5.2</td></tr><tr><td>Race</td><td></td></tr><tr><td>White/Causasian</td><td>367 (87)</td></tr><tr><td>Black</td><td>12 (3)</td></tr><tr><td>Asian</td><td>26 (6)</td></tr><tr><td>Other</td><td>15 (4)</td></tr><tr><td>Education level</td><td></td></tr><tr><td>No graduate degree</td><td>167 (40)</td></tr><tr><td>Graduate degree</td><td>253 (60)</td></tr><tr><td>Body Mass Index (kg/m<sup>2</sup>), mean ± SD</td><td>27.7 ± 5.4</td></tr><tr><td>Smoking status</td><td></td></tr><tr><td>Never smoker</td><td>282 (67)</td></tr><tr><td>Past smoker</td><td>115 (27)</td></tr><tr><td>Current smoker</td><td>23 (6)</td></tr><tr><td>Alcohol units per day, median (IQR)</td><td>0.3 (0.55)</td></tr><tr><td>Missing data</td><td>139</td></tr><tr><td>Coffee consumption per day, median (IQR)</td><td>1.0 (0.5)</td></tr><tr><td>Missing data</td><td>139</td></tr><tr><td>Tea consumption per day, median (IQR)</td><td>0.1 (0.5)</td></tr><tr><td>Missing data</td><td>139</td></tr><tr><td>Total physical activity (hours per week), mean (± SD)</td><td>7 (6)</td></tr><tr><td>Previous partner pregnancy (yes v. no)</td><td>158 (39)</td></tr><tr><td>Missing data</td><td>13</td></tr><tr><td>Previous varicocele (yes v. no)</td><td>39 (9)</td></tr><tr><td>Missing data</td><td>1</td></tr><tr><td>Sperm count (mil), mean ± SD</td><td>160 ± 166</td></tr><tr><td>Sperm concentration (mil/mL), mean ± SD</td><td>64 ± 72</td></tr><tr><td>Sperm morphology (% normal), mean ± SD</td><td>6.3 ± 3.7</td></tr><tr><td>Sperm motility (% motile), mean ± SD</td><td>43 ± 23</td></tr><tr><td>Semen volume (mL), mean ± SD</td><td>2.8 ± 1.4</td></tr><tr><td>Urine specific gravity, mean ± SD</td><td>1.017 ± 0.008</td></tr><tr><td>Number of urine samples, median ± IQR</td><td>2 ± 2</td></tr></tbody></table>

Abbreviations: N, number; IQR, interquartile range; SD, standard deviation.

2005 and 2017. They had an average age at enrolment of 36.4 (± 5.2) years (Table 1). Most men were Caucasian (87 %), whilst 3 % were Black, 6 % were Asian and 4 % self-identified with another racial group. The majority of men had a college and graduate degree (60 %) and had never smoked (67 %). The average hours per week of total physical activity was 7 h (±6). Less than half (39 %) of their partners had a previous pregnancy and 9 % of men had a history of varicocele. The yearly distribution of inclusion of men from the EARTH cohort included in our study population is available in Sup. Table S1. Men showed a trend (unadjusted) of increased sperm count and concentration over calendar year (Sup. Table S2) but no significant trends were identified for normal sperm morphology, sperm motility, and semen volume.

Among all 420 men, chemical biomarker detection frequencies ranged from 16.7 % (MMP) to 100 % (Table 2). Seventeen chemical biomarkers were detected in at least 60 % of samples. Median urinary chemical biomarker concentrations were similar to those reported among all male participants in the National Health and Nutrition Examination Survey (NHANES) (CDC, 2024). Urinary DEHP metabolite concentrations were highly correlated (r > 0.83) as were urinary methylparaben and propylparaben concentrations (r = 0.79). Except for MEP and MHiNCH, the other DINCH and phthalate biomarkers correlation coefficients were above r = 0.34. Correlations between families of chemical biomarkers were weak to moderate (−0.01 < r < 0.39).

## 3.2. Single-exposure biomarker associations

Using ExWAS, Hg was positively associated with sperm concentration (0.10 mill/mL, 95 % CI 0.01; 0.20), and close to positively associated with sperm count (0.08 mill, 95 % CI 0.00; 0.17) for each one log unit increase (Table 3). Triclosan was negatively associated with percent normal sperm morphology (−0.19 %, 95 % CI −0.34; −0.04), and mobility (−1.01 %, 95 % CI −1.97; −0.05). No significant associations were found for sperm count or volume though almost all chemical biomarkers had negative effects on semen volume. None of the associations remained significant when corrected for multiple testing (corrected p-value threshold = 0.006). These results were consistent with BKMR exposure-response relationships for each individual biomarker component while holding all other components at their median levels. BKMR also detected a significant positive association between Hg and sperm count and motility. In DSA (Table 4), Hg was positively associated with sperm count (0.10 mill, 95 % CI 0.06; 0.13), concentration (0.14 mill, 95 % CI 0.10; 0.17), and percent normal morphology (0.30 %, 95 % CI 0.18; 0.42). These relationships were and in the opposite direction for triclosan (−0.04 mill, 95 % CI −0.06; −0.02, −0.05 mill/mL, 95 % CI −0.12; −0.02, and −0.22 %, 95 % CI −0.28; −0.15, respectively). MnBP and MBzP were positively associated with sperm concentration (0.11 mill/mL, 95 % CI 0.05; 0.17, and 0.09 mill/mL, 95 % CI 0.03; 0.11, respectively) while MiBP was negatively associated with sperm concentration (−0.19 mill/mL, 95 % CI −0.25; −0.14).

## 3.3. Mixture effect

BKMR was the only model with the capacity of assessing overall mixture effects. When increasing the overall mixture of biomarkers concentration from the 10th percentile to the 90th percentile, the mixture was negatively associated with sperm morphology (−1.42 %, 95 % CRI −12.36; 9.52), motility (−0.10 %, 95 % CRI −1.78; 1.58), and volume (−0.26 mL, 95 % CRI −0.79; 0.28), and positively associated with sperm count (0.21 mill, 95 % CRI −0.23; 0.65) and concentration (0.33 mill/mL, 95 % CRI −0.14; 0.80), although these CRIs all included 0. The positive associations with the latter two outcomes were primarily driven by Hg (Fig. 1).

## 3.4. Interactions

Regarding interactions, the DSA model did not detect any

<footnote>a Values are presented as N (%) unless otherwise noted.</footnote>
<caption>**Table 2**</caption>

Unadjusted urinary bisphenol, triclosan, parabens, DINCH and phthalates biomarker (µg/L) and hair mercury concentrations (mg/L) among 420 men contributing 917 urine and 420 hair samples.

<table><thead><tr><th>Biomarker</th><th>Minimum LOD</th><th>Detection frequency (%)</th><th>q5</th><th>q25</th><th>Median</th><th>q75</th><th>q95</th></tr></thead><tbody><tr><td>MMP</td><td>0.3</td><td>16.7</td><td>&lt;LOD</td><td>&lt;LOD</td><td>&lt;LOD</td><td>&lt;LOD</td><td>2.3</td></tr><tr><td>MEP</td><td>1.2</td><td>100</td><td>3.1</td><td>15.4</td><td>39.2</td><td>126</td><td>740.5</td></tr><tr><td>MnBP</td><td>0.4</td><td>100</td><td>0.6</td><td>4.3</td><td>10.9</td><td>23.7</td><td>59.9</td></tr><tr><td>MHBP</td><td>0.3</td><td>54.3</td><td>&lt;LOD</td><td>&lt;LOD</td><td>0.5</td><td>1.1</td><td>9.8</td></tr><tr><td>MiBP</td><td>0.2</td><td>100</td><td>0.6</td><td>3.1</td><td>7.7</td><td>15.5</td><td>39.6</td></tr><tr><td>MHiBP</td><td>0.3</td><td>31.8</td><td>&lt;LOD</td><td>&lt;LOD</td><td>&lt;LOD</td><td>1.1</td><td>7.3</td></tr><tr><td>MBzP</td><td>0.2</td><td>100</td><td>0.2</td><td>1.2</td><td>3.4</td><td>8.0</td><td>27.3</td></tr><tr><td>MCPP</td><td>0.18</td><td>100</td><td>0.3</td><td>1.3</td><td>3.2</td><td>9.0</td><td>39.1</td></tr><tr><td>MEHP</td><td>0.2</td><td>100</td><td>0.3</td><td>0.8</td><td>2.1</td><td>5.8</td><td>38.5</td></tr><tr><td>MEOHP</td><td>0.2</td><td>100</td><td>0.5</td><td>2.5</td><td>7.6</td><td>19.0</td><td>122.2</td></tr><tr><td>MEHHP</td><td>0.2</td><td>100</td><td>0.9</td><td>4.4</td><td>13.0</td><td>31.7</td><td>208.4</td></tr><tr><td>MECPP</td><td>0.4</td><td>100</td><td>1.8</td><td>6.8</td><td>19.9</td><td>49.4</td><td>313.2</td></tr><tr><td>MNP</td><td>0.3</td><td>46.5</td><td>&lt;LOD</td><td>&lt;LOD</td><td>&lt;LOD</td><td>3.0</td><td>29.9</td></tr><tr><td>MCOP</td><td>0.3</td><td>94.4</td><td>&lt;LOD</td><td>5.6</td><td>18.5</td><td>67.3</td><td>327.2</td></tr><tr><td>MCNP</td><td>0.2</td><td>94.4</td><td>&lt;LOD</td><td>1.3</td><td>3.0</td><td>7.0</td><td>26.4</td></tr><tr><td>MHiNCH</td><td>0.3</td><td>43.9</td><td>&lt;LOD</td><td>&lt;LOD</td><td>&lt;LOD</td><td>0.3</td><td>1.4</td></tr><tr><td>MCOCH</td><td>0.3</td><td>20.5</td><td>&lt;LOD</td><td>&lt;LOD</td><td>&lt;LOD</td><td>&lt;LOD</td><td>0.7</td></tr><tr><td>BPA</td><td>0.1</td><td>100</td><td>0.2</td><td>0.6</td><td>1.3</td><td>2.5</td><td>6.9</td></tr><tr><td>Butylparaben</td><td>0.1</td><td>95.4</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td><td>7.3</td></tr><tr><td>Methylparaben</td><td>0.1</td><td>95.4</td><td>0.7</td><td>7.5</td><td>18.9</td><td>69.7</td><td>458.2</td></tr><tr><td>Propylparaben</td><td>0.1</td><td>95.4</td><td>0.1</td><td>0.4</td><td>1.6</td><td>9.5</td><td>115.0</td></tr><tr><td>Ethylparaben</td><td>0.1</td><td>35.7</td><td>&lt;LOD</td><td>&lt;LOD</td><td>&lt;LOD</td><td>0.7</td><td>5.4</td></tr><tr><td>Triclosan</td><td>1.0</td><td>65.2</td><td>&lt;LOD</td><td>&lt;LOD</td><td>1.6</td><td>12.7</td><td>275.8</td></tr><tr><td>Hg</td><td>0.01</td><td>100</td><td>45.2</td><td>339.6</td><td>625.1</td><td>1153.3</td><td>2551.0</td></tr></tbody></table>

N: number of participants; LOD: limit of detection; MMP: monomethyl phthalate; MEP: monoethyl phthalate; MnBP: mono-n-butyl phthalate; MHBP: mono-3-hydroxybutyl phthalate; MiBP: mono-isobutyl phthalate; MHiBP: mono-2-methyl-2-hydroxypropyl phthalate; MBzP: monobenzyl phthalate; MCPP: mono(3-carboxypropyl) phthalate; MEHP: mono(2-ethylhexyl) phthalate; MEOHP: mono(2-ethyl-5-oxohexyl) phthalate; MEHHP: mono(2-ethyl-5-hydroxyhexyl) phthalate; MECPP: mono(2-ethyl-5-carboxypentyl) phthalate; MNP: monononyl phthalate; MCOP: monocarboxyoctyl phthalate; MCNP: monocarboxy-isononyl phthalate; MHiNCH: cyclohexane-1;2-dicarboxylic acid-mono(hydroxy-isononyl) ester; MCOCH: cyclohexane-1;2-dicarboxylic acid-mono(carboxyoctyl) ester; BPA: bisphenol A; Hg: mercury.

<caption>**Table 3**</caption>

Single-chemical biomarker associations with semen parameters using ExWAS (N = 420).

<table><thead><tr><th rowspan="2">Chemical biomarker</th><th colspan="2">Count (mill)</th><th colspan="2">Concentration (mill/mL)</th><th colspan="2">Normal morphology (%)</th><th colspan="2">Motility (%)</th><th colspan="2">Volume (mL)</th></tr><tr><th>Beta</th><th>95 % CI</th><th>Beta</th><th>95 % CI</th><th>Beta</th><th>95 % CI</th><th>Beta</th><th>95 % CI</th><th>Beta</th><th>95 % CI</th></tr></thead><tbody><tr><td>MEP</td><td>0.02</td><td>-0.05; 0.09</td><td>0.04</td><td>-0.03; 0.11</td><td>-0.01</td><td>-0.26; 0.24</td><td>-0.02</td><td>-1.59; 1.54</td><td>-0.06</td><td>-0.15; 0.03</td></tr><tr><td>MnBP</td><td>0.01</td><td>-0.08; 0.09</td><td>0.04</td><td>-0.04; 0.13</td><td>-0.18</td><td>-0.51; 0.15</td><td>-0.41</td><td>-2.47; 1.64</td><td>-0.09</td><td>-0.21; 0.03</td></tr><tr><td>MiBP</td><td>-0.03</td><td>-0.12; 0.06</td><td>-0.04</td><td>-0.13; 0.06</td><td>-0.15</td><td>-0.48; 0.18</td><td>-1.01</td><td>-3.14; 1.13</td><td>0.02</td><td>-0.10; 0.15</td></tr><tr><td>MBzP</td><td>0.03</td><td>-0.05; 0.11</td><td>0.05</td><td>-0.03; 0.14</td><td>0.13</td><td>-0.17; 0.43</td><td>0.01</td><td>-1.93; 1.96</td><td>-0.02</td><td>-0.13; 0.09</td></tr><tr><td>MCPP</td><td>-0.01</td><td>-0.08; 0.07</td><td>0.02</td><td>-0.05; 0.10</td><td>-0.04</td><td>-0.31; 0.23</td><td>-0.82</td><td>-2.55; 0.91</td><td>-0.06</td><td>-0.16; 0.04</td></tr><tr><td>MEHP</td><td>0.01</td><td>-0.06; 0.08</td><td>0.02</td><td>-0.05; 0.09</td><td>-0.20</td><td>-0.44; 0.04</td><td>-0.02</td><td>-1.59; 1.55</td><td>-0.03</td><td>-0.12; 0.06</td></tr><tr><td>MEOHP</td><td>0.01</td><td>-0.06; 0.09</td><td>0.02</td><td>-0.05; 0.09</td><td>-0.19</td><td>-0.45; 0.08</td><td>-0.25</td><td>-1.90; 1.40</td><td>-0.02</td><td>-0.12; 0.07</td></tr><tr><td>MEHHP</td><td>0.02</td><td>-0.06; 0.09</td><td>0.02</td><td>-0.05; 0.10</td><td>-0.20</td><td>-0.46; 0.07</td><td>-0.26</td><td>-1.93; 1.40</td><td>-0.03</td><td>-0.12; 0.07</td></tr><tr><td>MECPP</td><td>0.04</td><td>-0.04; 0.11</td><td>0.03</td><td>-0.04; 0.11</td><td>-0.21</td><td>-0.48; 0.07</td><td>0.13</td><td>-1.59; 1.86</td><td>-0.01</td><td>-0.11; 0.09</td></tr><tr><td>MCOP</td><td>-0.04</td><td>-0.11; 0.03</td><td>-0.03</td><td>-0.09; 0.04</td><td>-0.05</td><td>-0.31; 0.20</td><td>-1.17</td><td>-2.80; 0.46</td><td>-0.03</td><td>-0.12; 0.07</td></tr><tr><td>MCNP</td><td>0.02</td><td>-0.07; 0.11</td><td>0.02</td><td>-0.07; 0.11</td><td>-0.23</td><td>-0.56; 0.10</td><td>-0.14</td><td>-2.27; 1.98</td><td>-0.01</td><td>-0.13; 0.11</td></tr><tr><td>BPA</td><td>-0.04</td><td>-0.14; 0.07</td><td>-0.04</td><td>-0.14; 0.07</td><td>-0.10</td><td>-0.54; 0.33</td><td>-1.75</td><td>-4.23; 0.74</td><td>0.00</td><td>-0.14; 0.14</td></tr><tr><td>Butylparaben</td><td>-0.01</td><td>-0.05; 0.03</td><td>-0.01</td><td>-0.04; 0.03</td><td>-0.10</td><td>-0.21; 0.02</td><td>-0.38</td><td>-1.11; 0.36</td><td>0.01</td><td>-0.03; 0.05</td></tr><tr><td>Methylparaben</td><td>-0.01</td><td>-0.10; 0.07</td><td>-0.01</td><td>-0.10; 0.08</td><td>0.00</td><td>-0.24; 0.24</td><td>0.30</td><td>-1.42; 2.02</td><td>-0.02</td><td>-0.11; 0.07</td></tr><tr><td>Propylparaben</td><td>0.00</td><td>-0.07; 0.06</td><td>0.01</td><td>-0.05; 0.08</td><td>0.03</td><td>-0.15; 0.20</td><td>-0.02</td><td>-1.30; 1.26</td><td>-0.04</td><td>-0.10; 0.02</td></tr><tr><td>Triclosan</td><td>-0.03</td><td>-0.08; 0.02</td><td>-0.03</td><td>-0.07; 0.02</td><td><strong>-0.19</strong></td><td><strong>-0.34; -0.04</strong></td><td><strong>-1.01</strong></td><td><strong>-1.97; -0.05</strong></td><td>-0.01</td><td>-0.07; 0.04</td></tr><tr><td>Hg</td><td>0.08</td><td>0.00; 0.17</td><td><strong>0.10</strong></td><td><strong>0.01; 0.20</strong></td><td>0.23</td><td>-0.12; 0.57</td><td>1.24</td><td>-0.58; 3.07</td><td>-0.08</td><td>-0.20; 0.04</td></tr></tbody></table>

N: number of participants; CI: confidence interval; MEP: monoethyl phthalate; MnBP: mono-n-butyl phthalate; MiBP: mono-isobutyl phthalate; MBzP: monobenzyl phthalate; MCPP: mono(3-carboxypropyl) phthalate; MEHP: mono(2-ethylhexyl) phthalate; MEOHP: mono(2-ethyl-5-oxohexyl) phthalate; MEHHP: mono(2-ethyl-5-hydroxyhexyl) phthalate; MECPP: mono(2-ethyl-5-carboxypentyl) phthalate; MCOP: monocarboxyoctyl phthalate; MCNP: monocarboxy-isononyl phthalate; BPA: bisphenol A; Hg: mercury.

All models were adjusted for age, body mass index, education level, specific gravity, and days since last ejaculate. All chemical biomarker concentrations were natural log-transformed.

interactions between chemical biomarkers. Similarly, no interactions were detected in BKMR models when inspecting plots of single-exposure associations when all the other exposure biomarkers concentrations were fixed at the 25th, 50th, and 75th percentiles. PIPs for interactions in NLInteraction models were all below 20 % excepted for a triclosan-butylparaben negative interaction in the sperm normal morphology model (PIP = 35.8 %) (Fig. 2), but this interaction was found to be not significant when tested in a linear regression model (data not shown).

### 3.5. Sensitivity analyses

In our sensitivity analyses using only the semen sample from the first IVF cycle, some differences in results were observed. MnBP was selected in DSA instead of MBzP in the sperm count model (0.09 mill, 95 % CI 0.03; 0.12) (Sup. Table S4), and was consistently associated in ExWAS (0.14 mill/mL, 95 % CI 0.04; 0.24) (Sup. Table S3), DSA (-0.06 mill/mL, 95 % CI -0.20; 0.07) and BKMR in the sperm concentration models
<table><caption>**Table 4**
Multiple chemical biomarker associations with semen parameters using DSA (N = 420).</caption><thead><tr><th rowspan="2">Biomarker*</th><th colspan="2">Count (mill)</th><th colspan="2">Concentration (mill/mL)</th><th colspan="2">Morphology (%)</th><th colspan="2">Motility (%)</th><th colspan="2">Volume (mL)</th></tr><tr><th>Beta</th><th>95 % CI</th><th>Beta</th><th>95 % CI</th><th>Beta</th><th>95 % CI</th><th>Beta</th><th>95 % CI</th><th>Beta</th><th>95 % CI</th></tr></thead><tbody><tr><td>Hg</td><td>0.10</td><td>0.06; 0.13</td><td>0.14</td><td>0.10; 0.17</td><td>0.30</td><td>0.18; 0.42</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>MnBP</td><td>-</td><td>-</td><td>0.11</td><td>0.05; 0.17</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>MiBP</td><td>-</td><td>-</td><td>-0.19</td><td>-0.25; -0.14</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>MBzP</td><td>-</td><td>-</td><td>0.09</td><td>0.03; 0.11</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Triclosan</td><td>-0.04</td><td>-0.06; -0.02</td><td>-0.05</td><td>-0.12; -0.02</td><td>-0.22</td><td>-0.28; -0.15</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

Bold = significant at the 5 % level.

All chemical biomarker concentrations were natural log-transformed.

CI: confidence interval; MnBP: mono-n-butyl phthalate; MiBP: monoisobutyl phthalate; MBzP: monobenzyl phthalate; Hg: mercury.

* Only effects for chemical biomarker exposures or variable selected in an interaction with a chemical biomarker in at least one of the models are shown. All models were adjusted for age, body mass index, education level, specific gravity, and days since last ejaculate.

(Supplementary Figure S3). An interaction term between Hg and MnBP was also selected in DSA in the sperm concentration model (0.03 mill/mL, 95 % CI 0.01; 0.05). MEP was consistently negatively associated with semen volume in all models. Lastly, no interaction with a PIP above 20 % was found in the NLInteraction models (Supplementary Fig. S4).

## 4. Discussion

In this study, our aim was to assess associations and interactions between 24 chemical biomarkers and semen parameters by using multiple models used in both mixture and exposome contexts and by evaluating the model robustness in the context of men consulting for couple infertility. Overall, some single-chemical biomarker associations were more consistently associated with several semen outcomes including Hg (positively associated with all semen parameters except for volume), and triclosan (negatively associated with all semen outcomes except for volume). Other chemical biomarkers showed sporadic associations and were significantly associated with sperm concentration (MBzP, MiBP and MnBP) and suggestive associations with sperm normal morphology (MBzP) or volume (MiBP and MnBP), but not for other outcomes. While these associations were not significant across all models, these trends were found in all of them. In our main models, the absence of interactions between chemical biomarkers was consistent across all models apart from one possible interaction between triclosan and butylparaben detected by NLInteraction for the sperm normal morphology outcome. These differences in results highlight the importance of considering several models when evaluating chemical mixtures.

Non-persistent chemicals such as phthalates, bisphenols and parabens have been shown to reduce the Leydig cells' ability to produce steroidogenic enzymes by increasing reactive oxygen species (ROS) and induce germ cell apoptosis in the testis (Presunto et al., 2023). In turn, ROS production causes DNA fragmentation. Parabens can also generate mitochondrial ROS in spermatozoa, leading to reduced motility (Boguenet et al., 2021). Moreover, these chemicals can bind to receptors on spermatozoa and change sperm functions, change the Y:X sperm chromosome ratio, and modify sperm protein composition through the deregulation of proteins; resulting in reduced semen parameters (Virant-Klun et al., 2022).

The positive associations we found between hair Hg and semen parameters were consistent with our previous analysis in the EARTH cohort (Minguez-Alarcón et al., 2018a,b). In this earlier study, adjustment for total fish (known to be the primary source of Hg exposure and whose intake improves semen parameters possibly due to the beneficial effects of the long chain omega-3 fatty acids fish may contain) and total calorie intakes, attenuated these associations. Whilst negative associations between triclosan and semen parameters are consistent with the literature (Adegbola et al., 2024), the positive associations found for MnBP and MBzP run contrary to our hypothesis. However, even though these biomarkers have generally been negatively associated with sperm concentration in epidemiological studies, it is not always the case (Cai

et al., 2015; Wang et al., 2023). The opposite results found between MnBP and MiBP in our DSA model despite high correlation between the metabolites could be the consequence of differences in metabolic pathways or biological activities (Hoffman et al., 2025). Some toxicological studies have shown that phthalates do not follow a linear dose-response function, which could lead to variations in associations depending on the level of exposure (Tian et al., 2019). Moreover, these associations could be due to false positive inflation, as evidenced by the lack of significance when corrected for multiple testing in our ExWAS model.

Nevertheless, differences in results between single-chemical biomarker associations in the literature and mixture or exposome associations found in this study are to be expected, provided the exposure biomarkers do not exhibit completely independent relationships with the outcome from the other exposure biomarkers. If interactions between chemical biomarkers exist and are captured by the models, differences in single-exposure and mixture effects will likely be observed. Moreover, a model such as DSA will only select exposure biomarkers and covariates that minimized the root mean squared error of predictions which could lead to variations in effect directions when accounting for different confounders. ExWAS simply corrects single-pollutant regression results for multiple testing and does not account for other exposure biomarkers. This model also does not allow different sets of confounders for different exposure biomarkers, which could also lead to differences in results if the causal structure surrounding exposures and outcomes are different across exposure biomarker concentrations.

In their exposome-wide association study of 128 chemical biomarkers, Chung et al. (2019) did not find any significant associations between phthalates, BPA, or Hg (measured in serum) and semen quality, but the study was performed in a general population sample (N = 473). Gao et al. (2024) found a negative association between overall mixture (six chemical biomarkers: BPA, ethylparaben, methylparaben, MCPP, benzophenone 1 and benzophenone 3) and sperm count, concentration, and motility, but not volume, in 155 men from a reproductive health center in China, using BKMR and quantile g-computation, and found no evidence of interactions between chemical biomarkers in BKMR models. In this study, MBzP was also non-statistically positively associated with semen parameters for medium concentrations (but not high concentrations) in single-exposure models.

ExWAS is a method borrowed from the genetic field which performs a high number of regression analyses while correcting for multiple testing, but does not allow to detect interactions unless pre-specified as a covariate. This model has a very high sensitivity (detecting true associations) but also a high false discovery rate despite its p-value correction feature (Agier et al., 2016). It is often used as a two-step approach to select exposure variables prior to incorporating them in another multiple exposure model as a way to increase its performance. However, we were not interested in doing so in this study but instead used ExWAS to obtain a baseline to compare the other model results too. The high threshold proposed by the multiple correction also makes it challenging
<figure><img src="image_2.png" /><figcaption>B</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>C</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>D</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>E</figcaption></figure>

<figure><img src="image_6.png" /></figure>

(caption on next page)
**Fig. 1.** BKMR associations of chemical biomarkers concentrations with semen parameters (N = 420) Outcomes: A, sperm count (mill); B, sperm concentration (mill/mL); C, sperm normal morphology (%); D, sperm motility (%); E, semen volume (mL). From left to right: 1. Overall association of the chemical biomarker mixture at different concentration percentiles compared to the 10th percentile and their 95 % credible intervals. 2. Exposure-response relationships for each chemical biomarker while holding all other biomarkers at their median concentrations. 3. Mean difference comparing the 75th to 25th percentile of each exposure biomarker (estimates and 95 % credible intervals) when all the other biomarker concentrations are fixed at the 25th, 50th, and 75th percentiles. All models were also adjusted on age, body mass index, education level, abstinence time and specific gravity. N: number of participants; MEP: monoethyl phthalate; MnBP: mono-n-butyl phthalate; MiBP: mono-isobutyl phthalate; MBzP: monobenzyl phthalate; MCPP: mono(3-carboxypropyl) phthalate; MEHP: mono(2-ethylhexyl) phthalate; MEOHP: mono(2-ethyl-5-oxohexyl) phthalate; MEHHP: mono(2-ethyl-5-hydroxyhexyl) phthalate; MECPP: mono(2-ethyl-5-carboxypentyl) phthalate; MCOP: monocarboxyoctyl phthalate; MCNP: monocarboxy-isononyl phthalate; BPA: bisphenol A; MP: methylparaben; PP: propylparaben; BP: butylparaben; Hg: mercury.

to appropriately select exposure biomarkers associated with a given outcome. Santos et al. (2020) pointed out that to achieve power of 80 % to detect the 95th effect sizes would require a sample size of 1000–2000 subjects to deal with 100 exposure biomarkers, which could also explain the lack of associations found in Chung et al. (2019). DSA is a selection method that offers a fairly good compromise between sensitivity and false discovery rate, that can be programmed to perform two-way interactions, including between potential covariates and exposure biomarkers, and can test quadratic terms. Nonetheless, its results require to be stabilized through iteratively launching computing DSAs a number of times using different seeds, thereby increasing computational time, and selecting variables that have been identified in a pre-specified number of iterations. This process ensures minimization of prediction errors but renders reproducibility difficult. It also means the selected exposure biomarker is not necessarily causally related to the outcome, but simply well correlated with one that may be. BKMR is a Bayesian method that uses iterations to select the most relevant variables and includes, as opposed to the other methods employed, the possibility of measuring the overall mixture effect, accounting for non-linearity and interactions. However, these are only interpretable using the visual inspection of plots, an issue that NLinteraction, another Bayesian method, tries to solve by probabilities of how much the model favors the inclusion of a given interaction. This model induces sparsity using multivariate spike and slab priors to determine which exposure biomarkers are associated with the outcome, and which exposure biomarkers interact with each other. Nonetheless, the advantages provided by BKMR have to be weighted with the necessity for large sample sizes given the effect, which impacts computational time.

Forest tree-based methods, such as GLINTERNET, seem like a promising yet underused option to perform variable selection and detect pairwise interactions. Unfortunately, all our models had large estimates of cross-validation error with all values of the regularisation parameter, and the instability of the variable selection over different values led us to remove the model from our analysis. However, other studies that have managed to fit the model, and doing so whilst using bootstrap samples to increase the stability of their findings, have emphasized the computational efficiency of this method but also the difficulty to interpret results (Cano-Sancho et al., 2023).

In addition to the ExWAS method that was performed to provide a baseline for single-exposure biomarker associations, the other methods were chosen for their capacity to allow confounder specification, multiple imputation (except GLINTERNET), integration of non-linearity, and the reduction of potential confounding by co-exposures. Nonetheless, despite these methods ability to detect interactions, it can be challenging to do so when the interaction is weak, especially in a context where power can be lacking, which is often the case with chemical exposures due to the cost of measuring these chemical biomarkers. This might explain the difficulties faced in this study to detect any interaction, and to confirm them across models when one is detected.

This study has several limitations. First, we only included Hg measured in hair and select non-persistent chemical biomarkers into our analyses whose concentrations can be reliably measured in urine because other families of pollutants such as persistent fluorinated or brominated chemical biomarkers are measured in blood and not urine. In addition, we did not include other metals like cadmium and lead which have previously been associated with semen parameters (Gao

et al., 2023; Giulioni et al., 2023) because these were not available in the cohort. Nonetheless, the correlation structure of the over 100 environmental exposure biomarkers measured in the European HELIX multi-cohort revealed that the correlation coefficients were stronger within than between families of exposure biomarkers, which may suggest that concentrating on a single family of exposure biomarkers at a time may not be confounded to a great extent by unmeasured exposures from other chemical families (Tamayo-Uria et al., 2019). Second, exposure measurement error is possible given the likely episodic nature of exposures and short biological half-lives (<24 h) of most of the chemicals evaluated (Koch et al., 2004; Shin et al., 2023; Völkel et al., 2002). Nevertheless, we previously demonstrated in the EARTH cohort that a single urine sample may reasonably represent an individual's exposure to short-lived chemicals such as BPA and phthalates over several months (Braun et al., 2012). Also, we employed several methods that have been shown to partially reduce variability in exposure biomarker concentrations including the collection of several urine samples to calculate an average concentration, and the adjustment of urinary concentrations using specific gravity (Bastiaensen et al., 2020; Hoppin et al., 2002; Roggeman et al., 2022). Third, men in this study were mostly White, with a high socioeconomic status and seeking fertility treatment, which prevents us from generalizing the associations to other populations and could make our sample not representative of the population attending fertility centers. However, semen parameters were not significantly different according to self-reported races (data not shown). Fourth, our analysis did not include sperm kinetic parameters such as straight-line velocity (VSL) or straightness (STR) which could show different patterns of associations with non-persistent chemicals because they were not included in the study protocol and were therefore not measured. We also did not assess the impact of chemical biomarkers on semen parameters using thresholds proposed by the WHO (World Health Organization, 2021) because we believed our sample size was too low to provide adequate power for dichotomous outcomes, and because these thresholds do not necessarily predict IVF success (Lamb and Marinaro, 2023). Lastly, reverse causation is of possible concern given the cross-sectional nature of the study. The main strength of this study is the use of diverse and modern methods that account for multiple imputation, non-linearity of effects, and known confounders to analyze mixtures and interactions between a large number of chemical exposure biomarkers in relation to several semen parameters. Another main strength is the collection of multiple urine and semen samples which reduces exposure and outcome misclassification.

Despite the increasing trend of multipollutant models, little attention has been given to evaluating interactions. Those who have attempted the latter have mostly used BKMR, which is known to require both large samples and computational time. Based on our results showing some variations in estimated associations between models, research studies including multiple chemical biomarkers should be encouraged to diversify their model portfolio to incorporate interaction assessment, assess the robustness of their results, and highlight reasons for potential differences to better inform the reader on model limitations. The selection of these models should be guided by the research objective of the study. Additionally, forest tree-based models should be further considered and be made more easily implementable through coding packages.
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

**Fig. 2.** NLinteraction posterior inclusion probabilities in association with semen parameters (N = 420). A, sperm count (mill); B, sperm concentration (mill/mL); C, sperm normal morphology (%); D, sperm motility (%); E, semen volume (mL). All models were also adjusted on age, body mass index, education level, abstinence time and specific gravity. N: number of participants; MEP: monoethyl phthalate; MnBP: mono-n-butyl phthalate; MiBP: mono-isobutyl phthalate; MBzP: monobenzyl phthalate; MCP: mono(3-carboxypropyl) phthalate; MEHP: mono(2-ethylhexyl) phthalate; MEOHP: mono(2-ethyl-5-oxohexyl) phthalate; MEHHP: mono(2-ethyl-5-hydroxyhexyl) phthalate; MECPP: mono(2-ethyl-5-carboxypentyl) phthalate; MCOP: monocarboxyoctyl phthalate; MCNP: monocarboxy-isononyl phthalate; BPA: bisphenol A; MP: methylparaben; PP: propylparaben; BP: butylparaben; Hg: mercury.
## CRediT authorship contribution statement

Maximilien Génard-Walton: Writing – review & editing, Writing – original draft, Methodology, Formal analysis. **Charline Warembourg:** Writing – review & editing, Methodology. **Paige L. Williams:** Writing – review & editing, Conceptualization. **Irene Souter:** Writing – review & editing, Data curation. **Jennifer B. Ford:** Writing – review & editing, Data curation. **Jorge E. Chavarro:** Writing – review & editing, Conceptualization. **Russ Hauser:** Writing – review & editing, Conceptualization. **Lidia Mínguez-Alarcón:** Writing – review & editing, Methodology, Conceptualization.

## Funding

The project was financed by Grants (R01ES022955, R01ES033651, R01ES009718 and P30ES000002) from the National Institutes of Health.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

The authors gratefully acknowledge all members of the EARTH study team, specifically the Harvard T. H. Chan School of Public Health research staff Myra Keller, Ramace Dadd and Alex Azevedo, physicians and staff at Massachusetts General Hospital Fertility Center as well as CDC lab personnel. We thank the National Institutes of Health for their funding. A special thank you to all of the study participants.

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.envint.2025.109568.

## Data availability

The authors do not have permission to share data.

## References

* Adegbola, C.A., Akhigbe, T.M., Adeogun, A.E., Tvrdá, E., Pizent, A., Akhigbe, R.E., 2024. A systematic review and meta-analysis of the impact of triclosan exposure on human semen quality. *Front. Toxicol* 6, 1469340.

* Agier, L., Portengen, L., Chadeau-Hyam, M., Basagana, X., Giorgis-Allemand, L., Siroux, V., Robinson, O., Vlaanderen, J., González, J.R., Nieuwenhuijsen, M.J., et al., 2016. A systematic Comparison of Linear regression-based statistical methods to assess exposome-health associations. *Environ. Health Perspect* 124, 1848–1856.

* Antonelli, J., Mazumdar, M., Bellinger, D., Christiani, D., Wright, R., Coull, B., 2020. Estimating the health effects of environmental mixtures using bayesian Semiparametric regression and Sparsity inducing priors. *Ann Appl Stat* 14, 257–275.

* Barrera-Gómez, J., Agier, L., Portengen, L., Chadeau-Hyam, M., Giorgis-Allemand, L., Siroux, V., Robinson, O., Vlaanderen, J., González, J.R., Nieuwenhuijsen, M., et al., 2017. A systematic comparison of statistical methods to detect interactions in exposome-health associations. *Environ. Health* 16, 74.

* Bastiaensen, M., Malarvannan, G., Gys, C., Ait Bamai, Y., Araki, A., Covaci, A., 2020. Between- and within-individual variability of urinary phthalate and alternative plasticizer metabolites in spot, morning void and 24-h pooled urine samples. *Environ. Res.* 191, 110248.

* Bobb, J.F., Valeri, L., Claus Henn, B., Christiani, D.C., Wright, R.O., Mazumdar, M., Godleski, J.J., Coull, B.A., 2015. Bayesian kernel machine regression for estimating the health effects of multi-pollutant mixtures. *Biostat Oxf Engl* 16, 493–508.

* Boguenet, M., Bouet, P.-E., Spiers, A., Reynier, P., May-Panloup, P., 2021. Mitochondria: their role in spermatozoa and in male infertility. *Hum Reprod Update* 27, 697–719.

* Braun, J.M., Gennings, C., Hauser, R., Webster, T.F., 2016. What can epidemiological studies tell us about the impact of chemical mixtures on human health? *Environ. Health Perspect* 124, A6–A9.

* Braun, J.M., Smith, K.W., Williams, P.L., Calafat, A.M., Berry, K., Ehrlich, S., Hauser, R., 2012. Variability of Urinary phthalate metabolite and bisphenol a concentrations before and during pregnancy. *Environ. Health Perspect* 120, 739–745.

* Cai, H., Zheng, W., Zheng, P., Wang, S., Tan, H., He, G., Qu, W., 2015. Human urinary/ seminal phalates or their metabolite levels and semen quality: a meta-analysis. *Environ. Res.* 142, 486–494.

* Cano-Sanche, G., Warembourg, C., Güil, N., Stratakis, N., Lertxundi, A., Irizar, A., Llop, S., Lopez-Espinosa, M.-J., Basagana, X., González, J.R., et al., 2023. Nutritional modulation of associations between prenatal exposure to persistent organic pollutants and childhood obesity: a prospective cohort study. *Environ. Health Perspect* 131, 037011.

* Caudill, S.P., Schleicher, R.L., Pirkle, J.L., 2008. Multi-rule quality control for the age-related eye disease study. *Stat. Med.* 27, 4094–4106.

* CDC, 2024. National Exposure Report Home Page. Natl Rep Hum Expo Environ Chem. Available at: https://www.cdc.gov/environmental-exposure-report/index.html. Accessed September 11, 2024.

* Chang, C., Dai, Y., Zhang, J., Wu, Z., Li, S., Zhou, Z., 2024. Associations between exposure to pesticides mixture and semen quality among the non-occupationally exposed males: four statistical models. *Environ. Res.* 257, 119400.

* Chung, M.K., Buck Louis, G.M., Kannan, K., Patel, C.J., 2019. Exposome-wide association study of semen quality: systematic discovery of endocrine disrupting chemical biomarkers in fertility require large sample sizes. *Environ. Int* 125, 505–514.

* Daza-Rodríguez, B., Aparicio-Marenco, D., Márquez-Lázaro, J., 2023. Association of triclosan and human infertility: a systematic review. *Environ Anal Health Toxicol* 38, e2023015.

* Eisenberg, M.L., Li, S., Behr, B., Cullen, M.R., Galusha, D., Lamb, D.J., Lipshultz, L.I., 2014. Semen quality, infertility and mortality in the USA. *Hum Reprod Oxf Engl* 29, 1567–1574.

* Gao, C., Sun, N., Xie, J., Li, J., Tao, L., Guo, L., Shi, L., He, X., Shen, X., Wang, H., et al., 2024. Co-exposure to 55 endocrine-disrupting chemicals linking diminished sperm quality: mixture effect, and the role of seminal plasma docosapentaenoic acid. *Environ. Int* 185, 108571.

* Gao, X., Li, G., Pan, X., Xia, J., Yan, D., Xu, Y., Ruan, X., He, H., Wei, Y., Zhai, J., 2023. Environmental and occupational exposure to cadmium associated with male reproductive health risk: a systematic review and meta-analysis based on epidemiological evidence. *Environ. Geochem. Health* 45, 7491–7517.

* Gibson E, Nunez Y, Abuawad A, R. Zota A, Renzetti S, L. Devick K, Gennings C, Goldsmith J, A. Coull B, Kioumourtzoglou M-A., 2019. An overview of methods to address distinct research questions on environmental mixtures: an application to persistent organic pollutants and leukocyte telomere length. *Environ Health 2019*; 18. Available at: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6714427/. Accessed April 12, 2021.

* Giulioni, C., Maurizi, V., De Stefano, V., Polisini, G., Teoh, J.-Y.-C., Milanese, G., Galosi, A.B., Castellani, D., 2023. The influence of lead exposure on male semen parameters: a systematic review and meta-analysis. *Reprod. Toxicol.* 118, 108387.

* Hauser, R., Meeker, J.D., Duty, S., Silva, M.J., Calafat, A.M., 2006. Altered semen quality in relation to Urinary concentrations of phthalate monoester and oxidative metabolites. *Epidemiology* 17, 682.

* Henriques, M.C., Loureiro, S., Fardilha, M., Herdeiro, M.T., 2019. Exposure to mercury and human reproductive health: a systematic review. *Reprod. Toxicol.* 85, 93–103.

* Hoffman, S.S., Tang, Z., Dunlop, A., Brennan, P.A., Huynh, T., Eick, S.M., Barr, D.B., Rushing, B., McRitchie, S.L., Sumner, S., et al., 2025. Impact of prenatal phthalate exposure on newborn metabolome and infant neurodevelopment. *Nat Commun* 16, 2539.

* Hoppin, J.A., Brock, J.W., Davis, B.J., Baird, D.D., 2002. Reproducibility of urinary phthalate metabolites in first morning urine samples. *Environ. Health Perspect* 110, 515.

* Juarez, P.D., Matthews-Juarez, P., 2018. Applying an exposome-wide (ExWAS) approach to cancer Research. *Front Oncol* 8, 313.

* Koch, H.M., Bolt, H.M., Angerer, J., 2004. Di(2-ethylhexyl)phthalate (DEHP) metabolites in human urine and serum after a single oral dose of deuterium-labelled DEHP. *Arch Toxicol* 78, 123–130.

* Kruger, T.F., Acosta, A.A., Simmons, K.F., Swanson, R.J., Matta, J.F., Oehninger, S., 1988. Predictive value of abnormal sperm morphology in in vitro fertilization. *Fertil. Steril.* 49, 112–117.

* Kusters, M.S.W., Essers, E., Muetzel, R., Ambrós, A., Tiemeier, H., Guxens, M., 2022. Air pollution exposure during pregnancy and childhood, cognitive function, and emotional and behavioral problems in adolescents. *Environ. Res.* 214, 113891.

* Lamb, D.J., Marinaro, J.A., 2023. Can semen parameters predict pregnancy outcomes? *Fertil. Steril.* 120, 709–714.

* Latif, T., Kold Jensen, T., Mehlsen, J., Holmboe, S.A., Brinth, L., Pors, K., Skouby, S.O., Jørgensen, N., Lindahl-Jacobsen, R., 2017. Semen quality as a predictor of subsequent morbidity: a danish cohort study of 4,712 men with long-term follow-up. *Am J Epidemiol* 186, 910–917.

* Levine, H., Jørgensen, N., Martino-Andrade, A., Mendiola, J., Weksler-Derri, D., Mindlis, I., Pinotti, R., Swan, S.H., 2017. Temporal trends in sperm count: a systematic review and meta-regression analysis. *Hum Reprod Update* 23, 646–659.

* Li, M.-X., Yeung, J.M.Y., Cherny, S.S., Sham, P.C., 2011. Evaluating the effective numbers of independent tests and significant p-value thresholds in commercial genotyping arrays and public imputation reference datasets. *Hum. Genet.* 131, 747.

* Meeker, J.D., Yang, T., Ye, X., Calafat, A.M., Hauser, R., 2011. Urinary concentrations of Parabens and serum hormone levels, semen quality Parameters, and sperm DNA damage. *Environ. Health Perspect* 119, 252–257.

* Messerlian, C., Williams, P.L., Ford, J.B., Chavarro, J.E., Mínguez-Alarcón, L., Dadd, R., Braun, J.M., Gaskins, A.J., Meeker, J.D., James-Todd, T., et al., 2018. The
Environment and Reproductive Health (EARTH) Study: a Prospective Preconception Cohort. Hum. Reprod. Open 2018, hoy001.

Mínguez-Alarcón, L., Afeiche, M.C., Williams, P.L., Arvizu, M., Tanrikut, C., Amarasiriwardena, C.J., Ford, J.B., Hauser, R., Chavarro, J.E., 2018a. Hair mercury (Hg) levels, fish consumption and semen parameters among men attending a fertility center. *Int J Hyg Environ Health* 221, 174–182.

Mínguez-Alarcón, L., Williams, P.L., Chiu, Y.-H., Gaskins, A.J., Nassan, F.L., Dadd, R., Petrozza, J., Hauser, R., Chavarro, J.E., 2018b. Secular trends in semen parameters among men attending a fertility center between 2000 and 2017: identifying potential predictors. *Environ. Int* 121, 1297–1303.

Presunto, M., Mariana, M., Lorigo, M., Cairrao, E., 2023. The effects of bisphenol a on human male infertility: a review of current epidemiological studies. *Int J Mol Sci* 24, 12417.

Roggeman, M., Gys, C., Klimowska, A., Bastiaensen, M., Wielgomas, B., Ait Bamai, Y., Covaci, A., 2022. Reviewing the variability in urinary concentrations of non-persistent organic chemicals: evaluation across classes, sampling strategies and dilution corrections. *Environ. Res.* 215, 114332.

Santos, S., Maitre, L., Warembourg, C., Agier, L., Richiardi, L., Basagaña, X., Vrijheid, M., 2020. Applying the exposome concept in birth cohort research: a review of statistical approaches. *Eur J Epidemiol* 35, 193–204.

Shin, M.-Y., Choi, J.W., Lee, S., Kim, S., Kho, Y., Choi, K., Kim, S., 2023. Pharmacokinetics of transdermal methyl-, ethyl-, and propylparaben in humans following single dermal administration. *Chemosphere* 310, 136689.

Smarr, M.M., Honda, M., Kannan, K., Chen, Z., Kim, S., Buck Louis, G.M., 2018. Male Urinary Biomarkers of antimicrobial exposure and bi-directional associations with semen quality Parameters. *Reprod. Toxicol. Elmsford N* 77, 103–108.

Sun, Z., Tao, Y., Li, S., Ferguson, K.K., Meeker, J.D., Park, S.K., Batterman, S.A., Mukherjee, B., 2013. Statistical strategies for constructing health risk models with

multiple pollutants and their interactions: possible choices and comparisons. *Environ. Health* 12, 85.

Tamayo-Uria, I., Maitre, L., Thomsen, C., Nieuwenhuijsen, M.J., Chatzi, L., Siroux, V., Aasvang, G.M., Agier, L., Andrusaityte, S., Casas, M., et al., 2019. The early-life exposome: description and patterns in six european countries. *Environ. Int* 123, 189.

Tian, M., Liu, L., Zhang, J., Huang, Q., Shen, H., 2019. Positive association of low-level environmental phthalate exposure with sperm motility was mediated by DNA methylation: a pilot study. *Chemosphere* 220, 459–467.

Virant-Klun, I., Imamovic-Kumalic, S., Pinter, B., 2022. From oxidative stress to male infertility: review of the associations of endocrine-disrupting chemicals (bisphenols, phthalates, and Parabens) with human semen quality. *Antioxidants* 11, 1617.

Völkel, W., Colnot, T., Csanády, G.A., Filser, J.G., Dekant, W., 2002. Metabolism and kinetics of bisphenol a in humans at low doses following Oral administration. *Chem. Res Toxicol* 15, 1281–1287.

Vrijheid, M., Fossati, S., Maitre, L., Márquez, S., Roumeliotaki, T., Agier, L., Andrusaityte, S., Cadiou, S., Casas, M., de Castro, M., et al., 2020. Early-life environmental exposures and childhood obesity: an exposome-wide approach. *Environ. Health Perspect* 128, 067009.

Wang, H., He, H., Wei, Y., Gao, X., Zhang, T., Zhai, J., 2023. Do phthalates and their metabolites cause poor semen quality? a systematic review and meta-analysis of epidemiological studies on risk of decline in sperm quality. *Environ. Sci Pollut Res* 30, 34214–34228.

Weisskopf, M.G., Seals, R.M., Webster, T.F., 2018. Bias amplification in epidemiologic analysis of exposure to mixtures. *Environ. Health Perspect* 126, 047003.

World Health Organization, 2025. WHO laboratory manual for the examination and processing of human semen, 6th ed. 2021. Available at: https://www.who.int/publi cations/i/item/9789240030787. Accessed April 18, 2025.