The study population was selected from the Led-Fertyl cross-sectional study “Lifestyle and Environmental Determinants of Seminogram and other Male Fertility-related Parameters,” a cohort of healthy men aged 18 - 40 years established between February 2021 and April 2023. Recruitment details have been described in detail elsewhere.^18,19 Briefly, the exclusion criteria accounted for severe chronic illnesses, reproductive disorders, major organ transplants, cardiovascular disease, HIV or hepatitis B/C, active or recent cancer, severe psychiatric or endocrine disorders, liver failure, use of certain medications (like antidepressants, corticosteroids, or immunosuppressants), significant recent weight loss, and any condition that might hinder adherence to the study protocol. The present study included 48 participants randomly selected from the pool of 200 participants with available biological samples and measured outcomes. All procedures followed were in accordance with the ethical standards of the Helsinki Declaration. The project protocol received approval from the Ethics Committee of the Pere Virgili Health Research Institute (Ref. CEIM: 181/2019). An online and written informed consent was provided from all participants.

Blood, urine and semen samples were collected during the same morning visit at the Hospital Universitari Sant Joan de Reus (Reus, Tarragona, Spain). The details for the biocollection have been described elsewhere.^17,19 Fasting spot urine samples were stored in aliquots in polypropylene tubes and kept at −80°C until the laboratory analysis. Semen samples were collected following a minimum of 3 days of sexual abstinence, with the P25-P75 range being 3–5 days. Samples were obtained through masturbation and deposited into sterile, standard polypropylene containers. After a 20-min liquefaction period at 37°C, semen quality parameters and chemical exposome were analysed. Venous blood samples were collected after a minimum fasting period of 8 h into 10 mL Vacutainer tubes containing lithium heparin as the anticoagulant, and were also kept at −80°C until the laboratory analysis.

Macroscopic parameters, including semen volume and pH, were assessed. Microscopic examination involved the use of a phase-contrast microscope and a computer-assisted sperm analysis (CASA) system SCA, Microptic, version 6.5.0.67 (Microptic, Barcelona, Spain). This comprehensive analysis covered conventional factors such as sperm count and concentration, sperm motility, sperm vitality, and sperm morphology. Collection and examination procedures adhered to the World Health Organization standards (WHO, 2021).²⁰ Briefly, sperm count, and concentration were measured with the 10× phase contrast objective and expressed as millions of spermatozoa per ejaculate or millions of sperm cells per mL, respectively. Sperm motility was evaluated in 200 spermatozoa by analysing various real-time images captured by the CASA system, with each sperm cell categorized as progressive motile, non-progressive motile, or immotile. Motility was further quantified as a percentage of the total motility observed, encompassing both progressive and non-progressive motility. Sperm vitality was assessed using the hypoosmotic swelling test (HOS test) at 60× magnification and analysing 200 sperm cells. Additionally, sperm morphology was examined utilizing the Hemacolor (Millipore, Sigma-Aldrich, Darmstadt, Germany) staining protocol, observing 200 sperm cells under 60× brightfield optics. Morphology assessment involved quantifying the percentage of normal forms or abnormalities in the head, midpiece, terminal piece or combined abnormalities.

Analytical-grade solvents including methanol (MeOH, HPLC-grade), acetonitrile (ACN, HPLC-grade), water (HPLC-grade), formic acid (>99% purity), ammonium acetate, and ammonium formate (≥99.0% purity) were obtained from Merck (Darmstadt, Germany). Ultrapure water was produced using a Milli-Q purification system (Aurium, PRO-VFT, Sartorius, Göttingen, Germany). Captiva 3 mL Non-Drip filter cartridges were obtained from Agilent (Madrid, Spain). Analytical standards and isotopically labeled internal standards (ISs), were purchased from LGC Standards (Teddington, UK), Merck (Darmstadt, Germany), and Sigma-Aldrich (St Louis, MO, USA).

A comprehensive chemical database was employed for wide-scope targeted screening (Table S1), previously applied in human biomonitoring studies in urine, blood plasma, and semen.^15-17. The database was obtained from the Bruker TargetScreener method (Bruker Daltonics, Bremen, Germany) and was complemented with additional in-house compounds, including MS/MS information and retention times (RT). The final database comprises 2029 individual compounds covering the following major chemical classes(subclasses): artificial sweeteners, industrial chemicals & transformation products (TPs) (eg, PFAS, flame retardants, plasticizers, benzothiazoles, benzotriazoles, surfactants, UV stabilizers or tyre additives), naturally occurring substances & stimulants & TPs (eg, food-derived, or tobacco-derived), personal care products & TPs (eg, UV-filters), pharmaceuticals (PhACs) & TPs, and plant protection products & TPs (eg, biocides, or preservatives). Analytes were selected based on their potential presence in the aforementioned biological matrices.

Previously validated methods have been used to conduct a chemical exposome profiling of more than 2000 compounds after minimal sample preparation, ensuring the preservation of a maximal number of chemicals of interest.^15,17 Briefly, simplified sample preparation was performed for the three matrices as follows: semen samples were thawed at room temperature, incubated at 37°C for 15 min to achieve liquefaction, and centrifuged (10 000 g, 10 min) to obtain seminal plasma. An aliquot of the supernatant (100 μL) was mixed with 300 μL acetonitrile for protein precipitation and centrifuged again (10 000 g, 10 min) before transfer to a chromatographic vial. Blood serum/plasma samples were thawed, 200 μL aliquots were mixed with 600 μL acetonitrile, vortexed, centrifuged (10 000 g, 10 min), and the supernatant transferred to chromatographic vials. Urine samples were centrifuged (3500 rpm, 5 min), filtered through Captiva cartridges (Agilent, Madrid, Spain), and 475 μL of filtrate mixed with 25 μL methanol in a vial. Sample extracts were injected in an ultrahigh-performance liquid chromatography (UHPLC) system with a Bruker Elute Pump HPG 1300 coupled to a QTOF Impact II (Bruker, Bremen, Germany). The chromatographic separation was performed on an Intensity Solo column (2.1 × 100 mm, 1.8 μm) from Bruker, preceded by a guard column, CORTECTS C18, 1.7 μm (2.1 × 5) mm from Waters (Milford, USA), thermostated in an air oven at 40°C. Instrument was operated in broadband collision-induced dissociation (bbCID), a data-independent acquisition (DIA) mode, at 3 scans per second. The aqueous mobile phase consisted of H2O: MeOH (99:1) with solvent modifiers consisting in 5 mM ammonium formate and 0.01% formic acid for ESI(+), and 5 mM ammonium acetate for ESI(−). The organic phase was MeOH with identical additives for each ionization mode, respectively.

Concerning quantification, matrix-matched calibration curves were used for chemicals with available standards. Matrix-matched calibration curves were prepared for each matrix under identical experimental conditions as the samples to ensure reproducibility and comparability of the results. For compounds without available matrix-matched calibration curves, semi-quantification was performed using a quantitative structure–ionization relationship (QSIR) model based on ionization efficiency (IE).^15,21 This approach uses a log-transformed ratio of calibration curve slopes relative to reference compounds [bisphenol G for ESI(−), O-desmethyl venlafaxine for ESI(+)], normalized to reflect ionization potential across compounds. Table S3 presents the chemicals used to construct the matrix-matched calibration curve employed in the semi-quantification model. To account for matrix effects, a limited number of calibrant compounds were projected onto the biological matrix. The pseudo-molecular ions ([M − H]⁻ and [M + H]⁺) were used as diagnostic signals to calculate peak areas and logIE values. For these semi-quantified compounds, limits of detection (LOD) were estimated based on a signal-to-noise ratio approach, defined as three times the noise level (LOD = 3×S/N), using the signal in sample.

Strict QA/QC protocols were implemented to ensure data integrity and minimize contamination or instrumental bias throughout sample processing and analysis. All glassware was cleaned with Milli-Q water and acetone, and baked at 450 °C for 6 h prior to use. Work surfaces were decontaminated daily using water and acetone. A total of 26 procedural blanks—accounting for ∼18% of the total sample number—were processed using HPLC-grade water instead of biological matrix, and subjected to identical extraction and analysis procedures. Blank signals plus three times the standard deviation were subtracted from sample peak areas prior to quantification to correct for background contamination. Instrumental performance was monitored by injecting methanol and a 50 ng/mL standard mixture every 20 injections. Samples were injected in a randomized sequence, including first blood plasma, then seminal plasma, and finally urine samples. Clothianidin-d₃ was added before extraction to monitor potential sample preparation losses, while an internal standard (IS) mixture was added prior to LC-HRMS analysis to correct for matrix effects and instrumental fluctuations. IS list is provided in Table S3. As the coefficient of variation (CV%) for all internal standards (IS) remained low along the batch (average <12%), the analyte areas were corrected using the mean IS area per sample, providing a robust general correction in the absence of compound-specific IS (Table S3). To ensure high mass accuracy, external calibration using sodium formate was performed before each analytical batch. Additionally, internal calibration was carried out in the first 15s of each run using sodium formate/acetate in a 2-propanol/water mixture (50:50, v/v), allowing real-time mass correction.

Descriptive analysis was conducted for demographics, individual and semen parameter data. Partition coefficients between biological matrices (seminal plasma/blood plasma and seminal plasma/urine) were calculated as the ratio of median concentrations for each detected compound. Chemicals with a detection frequency above 75% were used as continuous variables in regression analysis substituting the non-detected observations by the limit of detection (LOD)/√2.²² Chemicals detected between 10% and 75% of samples were modeled in regression models as binary variables (detected Versus non-detected). Spearman correlation analysis was conducted to characterize the bivariate associations between chemicals in the same biological matrix and between matrices. Multiple linear regression was conducted to characterize the associations between semen quality parameters (log-transformed) and chemical concentrations (log-transformed) or presence (binary) as independent variables. The models were adjusted by confounding variables selected a priori using a minimally adjusted causal structure considering age (continuous, years), abstinence (continuous, days) and smoking (binary). A set of sensitivity analysis were conducted to assess the influence of physical activity (Metabolic Equivalent of Task (MET), continuous, MET min/week), body mass index (continuous, kg/m²), professional activity or residence in the proximity of a petrochemical industry in the results (binary). Due to the limited sample size and the exploratory nature of the study, the results from regression analysis were not adjusted by false discovery proportion. In order to discern potential co-exposure patterns or interactions, multipollutant models were built using gradient boosting trees allowing to combine continuous and binary exposure variables. We implemented the method using the R packages “gbm” and “dismo.” The gbm.step function was computed using the input parameters setup 0.001 for the learning rate and 4 for the tree complexity, the 0.8 bag fraction at 0.8, and a 10‐fold cross‐validation. We computed a two-step approach considering first a full model with all variables followed up by reduced model with a subset of variables with contributions above 2%. The relative contribution of biomarkers to the overall model fit was used as a variable importance metric and the marginal effects were inspected with the partial dependency plots. All statistical analysis were conducted with R software v. 4.2.3.²³

The wide-scope target screening method enabled the identification of 61 chemicals in at least one sample across the three analyzed matrices including a variety of chemical families such as artificial sweeteners, biocides, flame retardants, food-derived compounds, PFAS, pharmaceuticals or tobacco smoke chemicals. The detection frequencies (%DF), limits of quantification (LOQ), median concentrations and whether concentrations were determined by quantification (n_serum= 12; n_{seminal plasma}= 15; n_urine= 15) or semi-quantification (n_serum= 30; n_{seminal plasma}= 30; n_urine= 34) of all detected chemical are summarized in Figure 1 and supplemental Table S5. The average blanks concentration was also expressed in Table S5. The number of chemicals frequently detected (eg, above 75% of samples) ranged between 9 in seminal plasma or urine and 15 in plasma. The majority of chemicals were detected in less than 25% of analyzed samples.

Twenty-six chemicals were found in the three matrices (Figure S2), and the number of chemicals identified in just one matrix ranged between 4 (blood plasma) and 10 (urine). Strong positive correlations between matrices were observed for most chemicals, particularly for tobacco smoke compounds (See diagonal in the correlation heatmaps of supplemental Figures S3-S5). In turn, some chemicals showed no association or even mild negative correlations, highlighting the specificity of the matrix. That is the case for instance of mono(2-ethylhexyl) phthalate (MEHP) or mono n/isobutyl phthalate (MxBP).

The average partition coefficient of chemicals between matrices showed an overall trend to concentrate in seminal plasma over urine, mostly in the range 10 - 100 fold, but in some cases, such as quinmerac and nicotinamide, that exceeded a 100-fold concentration in seminal plasma (supplemental Figure S6). In case of the ratio seminal plasma/plasma we observed a large proportion of chemicals with a ratio below the unit (ie, indicating higher concentrations in plasma) including perfluorooctanesulfonic acid (PFOS), MEHP or hydroxycotinine.

The multivariate regression analysis adjusting for covariates (abstinence, age and smoking) showed a wide variety of chemical families systematically associated with multiple parameters of poor sperm quality. These included the artificial sweetener acesulfame, biocides (nitenpyram or 1,2-Benzisothiazol-3(2H)-one [BIT]), personal care products (propyl- and methylparabens), the plasticizer bisphenol S (BPS), the surfactant N, N-Dimethyldodecylamine (DMDDA), the tire additive 1,3-dicyclohexylurea (DCU) and a list of tobacco smoke compounds (Figure 2). While tobacco smoke compounds showed strong correlations among matrices, the associations with semen quality parameters were substantially strengthened when exposures were assessed in seminal plasma. Some chemicals showed isolated statistical associations with individual semen parameters or in specific matrices. For instance, the presence of PFOS in blood plasma was negatively associated with seminal volume or pentachlorophenol (PCP) showed a negative association with the normal morphology. We also identified some food-derived chemicals positively associated with multiple semen quality parameters across different matrices including nicotinamide (vitamin B3 derivative) or theophylline.

Multipollutant gradient-boosting tree (GBT) models showed complex associations between chemicals with multiple semen quality parameters. The analysis reveals that chemical mixtures exhibit both individual and interactive effects on sperm function, with variable importance and directionality differing across endpoints (Figure 3 and supplemental Figures S10-S12). Several interactions substantially contributed to the models for sperm concentration, total motility, count and vitality, indicating that chemical co-exposures may be more important than individual chemical effects for these endpoints. Nicotinamide appeared to be one of the most influential factors for most models, interacting with multiple chemicals. Pentachlorophenol (PCP) appeared to be negatively associated to total motility, normal form and vitality, whereas negative contributions of PFOA or PFOS were revealed for total motility, normal form, volume or vitality. DCU appeared to be the top ranked chemical negatively associated to vitality, interacting with paraxanthine and caffeine. Overall, these results highlight the importance of considering both individual chemical toxicity and mixture interactions when assessing reproductive health impacts of environmental chemical exposures.