Exposome Profiling of Environmental Pollutants in Seminal Plasma and Novel Associations with Semen Parameters.

From 2020 to 2021, 158 couples undergoing a fresh IVF cycle at Sheba Medical Center, a tertiary medical center in Tel HaShomer, Israel, were approached, of which 100 consented. Couples were enrolled during ovarian stimulation and those with severe oligoasthenoteratozoospermia were excluded at the time of enrollment. Demographic, lifestyle, anthropometric, and medical history data were collected by trained clinic personnel. Standard semen parameters, including volume, percent motility, and total sperm count data were collected as part of routine IVF protocol following WHO protocol on semen sample treatment and semen parameter assessment. 23 A composite semen index was derived using principal component analysis (PCA). Semen samples were collected in a sterile plastic specimen cup after a recommended 2–7 day abstinence period, per standard IVF protocol. Following liquefaction, sperm concentration and motility were evaluated by Makler counting cell (Sefi-Makler). The liquefied semen was gently added to the density gradient medium (80%/40% PureCeption, Quinn’s, SAGE) and centrifuged at 600 g for 20 min to separate the seminal plasma. Samples were stored at −80 °C, defrosted in a 4 °C refrigerator, and then brought to room temperature prior to extraction for NTA. Extensive details of the targeted and nontargeted high resolution platform are provided in the Supporting Information , and are based on previously described methods. 16 In brief, seminal plasma was extracted using a modified QuEChERS extraction. A 500 μL volume of seminal plasma was mixed with 500 μL hexane/acetone:dichloromethane, vortexed for 30 s, and transferred to a QuEChERS tube (150 mg MgSO 4 + 50 mg C18). The tubes were centrifuged and the supernatant was transferred to a glass vial for collection. This was repeated two more times with 250 μL hexane/acetone/dichloromethane. The final extract was reduced to 150 μL under nitrogen, transferred to an amber autosampler vial, spiked with 10 μL of an internal standard solution containing 62.5 μg/L of phenanthrene d -10 and chrysene d -12 (AccuStandard) and with 10 μL of diluted retention time marker (AccuStandard DRH-TX-003-CNM), and sealed with a cap. The extracts were analyzed using a using a high-resolution Thermo Q Exactive Orbitrap MS equipped with a Thermo Trace 1300 GC and a TriPlus RSH Autosampler. Limits of detection were determined by injecting a calibration standard seven times, and calculated as [ s × t (d f , 1 – α = 0.99)]/ m where s is the standard deviation, t is the student’s t -value, d f is the degrees freedom, α is the significance level ( n = 7, α = 0.01, t = 3.14), and m is the slope of the calibration curve. Full details of the analytical methods, including the QC, analytical sequence, chromatography and mass spectrometry settings, and data processing, are provided in the Supporting Methods . Table S1 shows the recovery spikes of the contaminant mixture in fetal bovine serum. Table S2 shows the extraction recovery of the contaminants from NIST SRM. For all statistical analyses, we restricted the analysis of targeted measurements to chemicals with ≥40% detection rate, resulting in 43 chemicals. Any measurements below the limit of detection were replaced with half the value of limit of quantitation or the minimum detected, whichever was smaller. The nontargeted chemical data analysis was restricted to peaks that were detected in all samples, with a minimum abundance >50,000 in order to ensure peak quality, resulting in 814 features for analysis. An overview of the statistical approach is shown in Figure S1 . Chemicals measured via the targeted approach were categorized based on entries for each chemical on the EPA CompTox Chemicals Dashboard. 24 They were broadly categorized as PCB ( n = 46), PAH ( n = 11), pesticides were subclassified as chlorobenzene ( n = 7), organochlorine ( n = 18), triazine and triazole ( n = 4), organophosphate ( n = 6), and other pesticides ( n = 8), and all other chemicals were categorized as “other” ( n = 18). Further, we examined the correlation between detection and physicochemical properties obtained from the CompTox Dashboard, namely ExpoCast predicted median exposure, predicted octanol–water partition coefficient, predicted bioconcentration factor, and predicted biodegradation half-life. To characterize the nontargeted data, we used chemoinformatic tools to generate simplified molecular-input line-entry system notation by querying the chemical identifier resolver using chemical names, annotated by compound discoverer, using the webChem R package 25 (version 1.1.3). We determined the elemental composition and functional groups present using the rcdk package. 26 For the NTA, features were detected with 10 ppm mass tolerance, 10,000 total ion chromatogram threshold, signal-to-noise ratio of 3, and 99% allowable ion overlap. Each chromatogram was retention time-aligned using the carbon distribution marker (contains 9 alkanes; only compounds containing greater than 8 carbons were used since the compounds smaller than this eluted during the solvent delay) spiked into each sample and retention indices (RIs) were calculated for each peak detected. For suspect screening analysis, the RI of each peak was used to limit suspects during identification; the allowed maximum RI difference was 300. Compounds were identified by searching their mass spectra in the NIST Mass Spectra Library (NIST/EPA/NIH EI and NIST Tandem Mass Spectral Library Version 2.3) and a high-resolution library developed in-house using certified standards containing 354 unique compounds. A minimum match factor and reverse match factor score of 500 was used for assigning library matches. Peaks with scores <500 were not assigned an identification. Chemicals that matched to our in-house library were assigned level 1 annotation if they were also detected in our standard mixture, based on previously described scoring scheme. 27 For all models, we decided to include age, body mass index (BMI), smoking and cannabis use, and infertility diagnoses (male, female, unexplained) as covariates. These were selected due to their known association with male fertility. Education and whether someone lived or worked near agriculture were not included in the model as it was not associated with any semen parameter outcomes in our population. For the targeted organic pollutants, we applied an Exposome Wide Association Study (ExWAS) approach where we tested each chemical individually with the semen parameters as outcomes, analogous to GWAS studies. 28 Using linear regression models, chemicals with 40–70% detection rate were modeled as binary exposures (above and below detection) and those with >70% detection rate were modeled as continuous variables. This was done to balance the statistical power gained by modeling variables as continuous when possible and the overall model fit as variables with high nondetects may be less appropriate for linear regression. Models were adjusted for age, BMI, smoking and cannabis use, and infertility diagnoses (male, female, unexplained). The model estimates were corrected for false discovery rate. 29 As a sensitivity analysis, we excluded 13 individuals who were diagnosed with male factor infertility and/or had total sperm count <39 million, which corresponds to fifth percentile of men whose partners became pregnant within 12 months. 23 Varying this cutoff made no material difference to our model estimates. For the 814 highly abundant nontargeted features, we used linear regression as a verification of our PCP–PCA-BKMR results (see below) because potential bias introduced by coexposure adjustment may lead to spurious findings in statistical mixture models under specific scenarios. 30 , 31 As all features were in arbitrary units that scale with concentration, we decided to model the feature variables after log 2-transformation to allow the resulting model estimates to be interpreted as the difference in outcome per doubling of exposure. This analysis was restricted to exposures identified by the PCP–PCA-BKMR pipeline. Models were adjusted for age, BMI, smoking and cannabis use, and infertility diagnoses. PCP is a robust method for dimensionality reduction and pattern recognition. 32 , 33 Its theory and application to environmental health research has been described previously. 19 , 20 In brief, we can consider any set of exposure data to comprise two underlying matrices—the low-rank L-matrix and a sparse S-matrix. The L-matrix represents the latent underlying patterns that can be estimated but are not directly observed. The S-matrix represents the unusual or uniquely rare values that cannot be explained by the latent patterns observed by the L-matrix. By partitioning the observed data with PCP, we can effectively separate the underlying latent patterns from the unusual and rare events. PCP offers two unique advantages. First, by removing the influence of unusual and rare events, we can remove outliers and statistical noise that are often products of random variability. Second, while the resulting L-matrix has the same dimensions as the original exposure matrix, it will have a low rank, effectively reducing the dimensionality of the data. Any subsequent matrix factorization techniques, such as PCA, will identify only a few factors/components that explain a large proportion of variability. This stands in contrast with PCA on typical observed—and noisy—exposure data matrices (i.e., raw exposure data) where a large number of factors/components are required to explain the majority of the variability in the data. For our presented analysis, we used non-negative nonconvex square root PCP with a noise-independent universal choice of regularization parameters. 19 , 34 This approach requires estimation of two parameters, λ and μ. We tested both default universal parameters (λ = 1/√ n and μ = p /√2 where n is sample size and p is number of features) 32 , 34 as well as parameter values obtained via a grid search. In each case, the resulting L matrices were very similar both in rank and subsequent regression results. Thus, we chose to present the results using the default universal parameters which have been shown to have theoretical guarantees. PCP was applied separately for the targeted chemicals and on the nontargeted spectra. The input for the targeted chemicals included 43 pollutants with ≥40% detection rate. The input for the nontargeted chemicals included 814 high abundance features that were detected in all samples. All inputs were scaled to 0–1 for PCP. We applied PCA on the L-matrices and for downstream regression analyses, we chose the top 5 PCs from the priority common organic pollutants (99.1% cumulative variance) and the top 6 PCs from the nontargeted spectra (99.4% cumulative variance). We extracted the eigenvectors as well as the loadings of each chemical in the data set. Linear regression was then used to model each PC as the predictor and semen parameters (total motile sperm, concentration, and percent motility) as outcomes. All models were adjusted for age, BMI, smoking and cannabis use, and infertility diagnoses (male, female, unexplained). For PCs that were found to be associated with outcomes, we identified the top chemicals or features with the highest loadings as potential drivers and used this to inform our mixtures models. BKMR 21 is a flexible statistical mixtures method that can model multiple exposures to estimate the dose–response relationship of each mixture component while holding other components constant. Bayesian kernel machine regression (BKMR) makes no assumptions regarding the shape of the relationships and can be used to capture nonlinear curves and statistical interactions between the mixture components. For all BKMR models, we scaled observed exposures, outcomes, and continuous covariates, and specified 50,000 iterations with default priors. All models were adjusted for age, BMI, smoking and cannabis use, and infertility diagnoses. From PC regression, we identified candidate chemicals or features that may be associated with semen parameters. We then modeled the observed continuous exposure values (i.e., observed concentrations) of these chemicals or features in BKMR models with our outcomes. We reasoned that real associations should be detectable using the observed values and that using the latent representations (i.e., L-matrix) may result in spurious findings induced by the data handling procedure. Furthermore, models using observed exposure values are more easily verified and replicated by others. Also, to ensure that our observations regarding individual chemicals were not products of arbitrary decisions for the number of mixture components, we tested BKMR models with top 5, 10, and 15 chemicals or features to ensure that the estimated dose–response relationships were a faithful reflection of the underlying data. All analyses were conducted in R v4.0.1. We used pcpr package ( https://github.com/Columbia-PRIME/pcpr ) for PCP and bkmr package for BKMR.

Population characteristics are shown in Table 1 . Among the male partners of the 100 couples seeking assisted in vitro fertilization, five reported use of cannabis and 19 reported current smoking. Their mean age was 38.1 (SD = 5.7) years old and mean BMI was 26.3 (SD = 4.6). A substantial portion of the study population lived or worked near agricultural sites. The population semen parameter was similar to the WHO population parameters 23 and the majority of participants had normal volume and concentration according to WHO standards. Only 9 couples sought IVF treatment due to male-factor infertility. Our study design, methods, and key findings are shown in (Graphical Abstract). Missing data on 2 individuals. Range is 8–21 years of formal education. Median is 15 years. Starting with 119 commonly measured persistent or trace organic pollutants (all with level 1 identification confidence, i.e. confirmed chemical structures), we found that nearly all tested organic pollutant exposures were detectable in at least one sample, including all measured classes of pesticides, polyaromatic hydrocarbons (PAH), polychlorinated biphenyl (PCBs), and other common environmental contaminants such as dioxins and dioxin-like compounds, phthalates, and solvents ( Figure 1 ). Cyanazine was the only chemical not detected in any of the samples. We observed predominantly positive correlations among the targeted chemicals, which is consistent with the fact that we expect chemical exposures to arise from complex chemical mixtures and not be mutually exclusive. Summary statistics for individual compounds can be found in Table S3 . High-resolution gas chromatography mass spectrometry efficiently detects widespread detection of organic pollutants in seminal plasma. Left side—the detection rate of 118 targeted organic pollutants. Right—correlation plot of chemicals with ≥40% detection rate. As expected we found a positive correlation between detection rate and predicted median exposure and (rho = 0.57, p < 0.001). We found negative correlations between detection rate and predicted octanol–water partition coefficient (rho = −0.56, p < 0.001), predicted bioconcentration factor (rho = −0.52, p < 0.001), and predicted biodegradation half-life (rho = −0.19, p < 0.05). Several specific exposure patterns can also be observed in our data ( Figure 1 ). There were expected correlations between different isomers and related metabolites such as lindane (γ-hexachlorohexane [HCH]), α-lindane, δ-lindane, and β-HCH. There was also a cluster that comprised numerous pesticides, such as triadimefon, N , N -diethyl-3-methylbenzamide (DEET), lindane (and its related isomers and metabolites), 2,4,5-trichlorophenol, etridiazole, 1,4-dichlorobenzene, and chlopyrifos, and reagents that are commonly used for production of pesticides, including diethyl phthalate and isosafrole. PCB-49 is also part of this cluster of correlated chemicals, which may reflect the usage of PCBs in certain pesticides similar to diethyl phthalate and isosafrole. Other specific coexposure patterns and exposure profiles were observed among dibenzofuran, safrole, and two PAHs, fluorene and acenaphthene. Although it is not clear which products or behaviors underlie this profile, these four chemicals also had positive correlations with many pesticides in our data. We used two parallel approaches to quantify the associations between the targeted organic pollutants and semen parameters, including motility, concentration, total motile sperm, and a combined index that captures overall poor semen parameters ( Figure S2 ). First, we testing each exposure individually using an ExWAS approach, adjusting for age, BMI, smoking and cannabis use, and infertility diagnoses and correcting for multiple comparisons. We found that etridiazole was associated with lower total motile sperm (FDR q = 0.01), concentration (FDR q = 0.07), and overall index of semen parameters (FDR q = 0.004) ( Figure 2 A,B, Table S4 ). These results were robust when we excluded those with male factor infertility or very low total motile sperm count ( Table S5 ). Exposome-wide association study using standard single-exposure regression models and machine learning pattern recognition based approach both identified etridiazole as negatively associated with worse semen parameters including lower total motile sperm, percent motility, and concentration. (A) Volcano plot of exposome-wide association study results, with the dashed line representing Bonferroni-corrected p-value cutoff threshold. (B) Boxplot of semen parameters among those with detectable levels of etridiazole compared to those with nondetectable levels of Etridiazole. (C) Heat plot of results PC regression. Green indicates negative association and the asterisk denotes statistical significance after Bonferroni correction. (D) Bar graph of posterior inclusion probability calculated from BKMR. (E) Estimated dose-response plots and associated 95% credible intervals from BKMR. Green line and credible interval indicates that the credible interval deviates from the null. We then sought to use our novel analysis approach, which combines a machine learning pattern recognition approach with modern mixture methods, to identify latent patterns in the exposure data and evaluate these exposure patterns, rather than individual exposures, with the semen parameter outcomes. In this approach, we first used PCP, a robust pattern recognition and data dimension reduction technique, to derive a low-rank exposure matrix that represents the latent underlying patterns without outliers and rare events ( Figure S3 ). Despite the apparent complexity in the raw observed data, ∼5 PCs (PCP–PCs) explained >99% of the variability in the low-rank exposure matrix ( Figure S4 ). Of these 5 PCP–PCs, we found that PCP–PC2 was associated with total motile sperm after adjusting for age, BMI, smoking and cannabis use, and infertility diagnoses and correcting for multiple comparisons ( Figure 2 C). The loadings on PCP–PC2 showed high positive loadings for PAHs and triazine/triazole pesticides and highest negative loadings for etridiazole and 2,3,7,8-tetrachlorodibenzofuran ( Figure S5 ). To formally investigate which exposures contributed to the association between PCP–PC2 and total motile sperm, the original observed exposure concentrations for the top 10 chemicals with highest absolute loadings were fit as a mixture in a BKMR model adjusting for age, BMI, smoking and cannabis use, and infertility diagnoses. The BKMR model showed that within this mixture, etridiazole had the highest posterior inclusion probability ( Figure 2 D) and was the only one that showed a negative association with any tested semen parameter ( Figure 2 E), which is evidence that etridiazole was the sole driver of this relationship. Indeed, we also saw that etridiazole was clearly associated with lower sperm motility in the BKMR model, which was also evidenced in the ExWAS models, but did not pass multiple testing corrections. To ensure that our results were robust and not specific to model parameters, we conducted sensitivity analyses where we specified alternative model parameters for PCP and the number of mixture components for BKMR. There were no appreciable changes in the results as etridiazole was consistently associated with poor semen parameters ( Figures S6–S7 ). In addition to the set of targeted organic chemicals, we conducted a NTA of all high-quality spectra tentatively identified in seminal plasma with the goal of applying and leveraging the exposomic approach to identify potential male reproductive toxicants among less-investigated pollutants. We characterized the tentatively identified spectra for elemental composition and presence of functional groups ( Figure S8 ). We found that the 39% of the tentatively identified features contained only hydrogen and carbon and ∼45% did not carry any functional groups. Starting from 814 spectral peaks with abundances detected in all study samples, we again applied our PCP-based approach. In these nontargeted data, PCP derived a low-rank exposure matrix with 6 components that explained >99% of the variance ( Figure S9 ), one of which was associated with total motile sperm ( Figure 3 A). Using BKMR model with the top 10 loading peaks of this component, we found that one peak had the highest posterior inclusion probability for total motile sperm, concentration, and combined index ( Figure 3 B). This peak had a retention time of 5.734 min and the mass-to-charge ( m / z ) ratio was 57.03352 and was the only one in this mixture that showed negative association with semen parameters ( Figure 3 C). These results did not differ when we considered alternative BKMR inputs ( Figure S10 ). Machine learning pattern recognition based approach identified a spectral peak to be negatively associated with worse semen parameters. (A) Heat plot of results PC regression. Green indicates negative association and the asterisk denotes statistical significance after Bonferroni correction. (B) Bar graph of posterior inclusion probability calculated from BKMR. (C) Estimated dose-response plots and associated 95% credible intervals from BKMR. Green line and credible interval indicates that the credible interval deviates from the null. (D) Forest plot of estimates and associated 95% confidence intervals from standard single-exposure linear regression models. To verify the observations from our PCP- and BKMR-based approach, we also modeled each of the top 10 peaks individually with linear regression, adjusting for age, BMI, smoking and cannabis use, and infertility diagnoses. These models also showed that the same peak (retention time = 5.734 min, m / z = 57.03352) was negatively associated with total motile sperm ( p = 0.01) and combined index ( p = 0.01) and marginally associated with lower sperm concentration ( p = 0.07) ( Figure 3 D). These results persisted after excluding those with male factor infertility or very low total motile sperm count ( p = 0.001 for total motile sperm, p = 0.02 for concentration, p = 0.003 for combined index). In order to determine the chemical identity of the peak, we used a suspect screening analysis pipeline described previously 16 and identified NDEA as a high quality match for this m / z feature with exact mass match, reverse match factor >600, and near exact retention time match (within 0.0008%). We compared this feature in our samples to the NDEA reference standard and found very similar chromatograph profile ( Figure S11 ) and ion ratios ( Table S6 ). Together, this provides level 1 identification confidence for this m / z feature as NDEA.

In this study, we show that numerous organic pollutants can be detected in seminal plasma and can be observed as several pollutants frequently co-occur with a range of other pollutants, forming distinct exposure profiles. Chemicals with low octanol–water partition coefficient, bioconcentration factor, and biodegradation half-life had higher detection rates among the participants. Additionally, we applied a novel pattern recognition approach to the exposure data to reduce dimensionality and subsequently identified etridiazole and NDEA as two potential male reproductive toxicants negatively associated with semen parameters. These discoveries were supported by the classical ExWAS approach and linear regression analyses, demonstrating the effectiveness and validity of our novel data analysis approach, and persisted when men with male factor infertility and low sperm count were removed. On a broader level, our study shows that we can adapt high dimensional approaches to exposomic studies of environmental determinants of reproductive health, which may be necessary to better understand the environmental contributions of the global male fertility decline. Pesticides have been previously link to semen parameters. 35 , 36 To our knowledge, this is the first study to report the potential association of etridiazole, a pesticide commonly used to control rot due to fungi and oomycetes, 37 and male fertility parameters. Etridiazole is sold as the active ingredient in numerous commercial pesticides such as terrazole, Truban, and Banrot and has been used on golf courses, cotton ornaments, lawns, and agricultural seed products. 37 In addition to occupational exposures, 38 − 40 etridiazole can be present in water and air near sites where it is frequently applied, 37 , 41 which may lead to low levels of intermittent exposure in the general population. Etridiazole is not typically found on food, 37 , 42 which suggests that the observed exposure is unlikely to be a result of confounding via dietary patterns and associated habits. Ultimately, there is likely substantial variation in detection rate of etridiazole across different populations, depending on underlying exposure levels and the analytical methods. It is notable that a substantial portion of our population live or work near agricultural sites, which likely led to higher detection rates than other populations. Etridiazole is classified as a probable human carcinogen, 43 , 44 in part due to evidence of testes tumor and testicular interstitial cell hyperplasia in rats, 37 which suggests potential reproductive toxicity. Mechanistically, etridiazole inhibits fungal growth by hydrolysis of phospholipid membrane of fungal mitochondria via increased phospholipase A activity. 45 The effect of etridiazole on mammalian cells is not well characterized, but there is evidence that it causes hemolysis and lipid peroxidation of human erythrocyte cell membranes by free radicals. 45 Thus, the biological mechanisms underlying the associations of etridiazole with semen parameters should be investigated in future studies. NDEA is a nitrosamine and a probable human carcinogen 43 , 44 with well characterized hepatotoxic, carcinogenic, and mutagenic properties. 46 NDEA is used in a wide range of industrial applications, including as a gasoline and lubricant additive, antioxidant, plastics stabilizer, and can be found in tobacco 47 , 48 and food products. 49 It can also be found as a water treatment disinfection byproduct. 50 , 51 A well-characterized mechanism of carcinogenesis by NDEA is through lipid peroxidation and generation of free radicals. However, lipid peroxidation and free radicals may result in a variety of other adverse health effects, including impaired spermatogenesis and poor semen parameters. 52 , 53 Although there was no prior human evidence, there is compelling evidence demonstrating reproductive toxicity of NDEA in animal models. NDEA increased abnormal sperm and markedly decreased sperm count, sperm motility, and male reproductive organ weight in chronically exposed rats in two different experiments. 54 , 55 These changes were complemented by depletion of antioxidants, increase in malondialdehyde and other markers of oxidative stress, elevated indicators of lipid peroxidation in the testes, pathological changes in seminiferous tubules, and changes in levels of key sex hormones. 54 , 55 Similar induction of oxidative stress, changes in sex hormones, and pathological changes in seminiferous tubules were observed in rabbits. 56 Lycopene, an antioxidant, appeared to rescue the NDEA-induced effects, suggesting that oxidative stress is the underlying mechanism. 54 Thus, our study is consistent with the emergent experimental evidence and there is strong evidence that NDEA can lead to impaired spermatogenesis and reduced fertility in humans. In addition to the discovery of two new potential human reproductive toxicants, our study also introduces a new workflow for the analysis of high dimensional exposome data. Our approach leverages a popular machine learning pattern recognition approach to remove outliers and rare events from noisy exposure data to uncover latent patterns. The resulting latent patterns can then use other machine learning and statistical techniques such as factor analysis (e.g., PCA) and mixtures analysis (e.g., BKMR) to screen hundreds of exposures with a smaller number of statistical tests while maintaining the ability to assess potential mixture effects and interactions. As proof of concept, we show that this workflow was able to recapitulate associations that were detectable in classical ExWAS analyses while identifying associations present in our data that were not detectable via ExWAS analyses due to lack of power. In effect, our proposed workflow can maximize statistical efficiency and overcome a key challenge in exposome studies 15 to efficiently screen hundreds of exposures using realistic study sample sizes in reproductive health studies to identify plausible reproductive toxicants. There are additional considerations when applying this approach. First, we were able to identify specific patterns and chemicals related to our outcomes, but it is possible that PCs with a highly toxic chemical will not be associated with the outcome of interest if the rest of the chemicals loading highly on the pattern were null or protective, thereby diluting the effect of the highly toxic chemical. It may be advisible to conduct parallel ExWAS analyses, but whether ExWAS can identify such associations depend on the statistical power. Second, a collection of weakly toxic chemicals may have a significant joint effect without a single primary suspect. Although we did not observe this in our study, this is another strength of BKMR that should be kept in mind during the interpretation. Lastly, it is important to note that the individual chemical estimates are the product of a prior selection step, which means that the estimates do not reflect uncertainty in dimension reduction. While our results are promising and consistent with experimental evidence, it raises some questions for future work. Given the absence of recent biomonitoring data for these chemicals, a comparison of exposure values in our study population to other populations cannot be conducted. For example, we detected relatively high levels of dioxin-like chemicals and etridiazole in our study compared to previous studies, 57 , 58 but there are major differences between populations, biomatrices, and assay methods, as well as potential differences across time. Ultimately, this highlights the urgent need for better human exposure biomonitoring efforts for known environmental contaminants. 59 While there are several advantages to measuring parent compounds directly, future efforts may also consider measuring known human metabolites. For example, etridiazole 60 and NDEA 61 have metabolites that may be detected at higher frequencies and quantities, making those metabolites potentially more attractive for targeted studies. Low detection may lead to nondifferential measurement error, which on average biases the results toward the null. Thus, while this likely did not result in spurious findings, in the case of etridiazole, additional efforts involving etridiazole metabolites and other populations would be helpful because the mechanism of action is unclear and complementary experimental evidence is necessary to understand how these putative male reproductive toxicants act mechanistically. Our study also has some limitations. First, our study population was recruited from couples seeking IVF in Israel and the associations observed with etridiazole and NDEA need to be replicated in other populations. However, there is compelling support for the validity of the observed associations from prior literature and experimental evidence. Although recruitment of men from infertility center might not be representative of the general population, most of those included in our cohort did not have male factor infertility: 42% of the participants underwent IVF for preimplantation genetic diagnosis with normal sperm parameters, 9% underwent planned embryo cryopreservation and 32% due to female infertility. Only a few study participants were diagnosed with male factor infertility and had obvious poor semen parameters. Furthermore, our results were stronger when we excluded individuals with male factor infertility or very low total motile sperm, which gives us confidence that these associations are likely also present in the general population. Additional limitation was that we assessed only three sperm parameters and not the full semen quality parameters according to the WHO. Upon the first visit in our unit, male partners were asked to provide a full semen analysis, including sperm morphology. However, these tests were performed in different laboratories and variable timing and only sperm volume, percent motility, and concentration were collected in our clinical IVF lab on the day of sample collection. Another limitation was that we only assessed semen parameters as the evidence suggested a global decline in sperm count. However, the relationship between semen parameters and fecundity is nonlinear 7 , 62 and it is unknown whether etridiazole or NDEA are associated with measures of fecundity such as time-to-pregnancy. Lastly, we were unable to adjust for other potentially relevant factors such as diet, which could influence both our observed exposures and outcomes.

Infertility, defined as the inability to achieve a pregnancy after ≥12 months of trying, is a global public health concern that affects 15% of all reproductive age couples in the world. 1 , 2 Infertility has significant negative physical, emotional, and social impacts on the affected individuals; addressing infertility will help realize the fundamental right for individuals and couples to have children and to improve the health of the affected individuals. Male infertility contributes to 40–50% of overall cases. 3 While changes in gross rates of infertility have not been observed, 2 many alarming trends have been observed in male fertility, including notable declines in testosterone and sperm count, with concurrent rises in male genital anomalies and testicular cancer. 4 Recent data show that sperm count declined by 51.6% between 1973 and 2018 and this decline has accelerated since 2000. 5 Thus, if this trend persists, a greater proportion of the male population could drop below the threshold 6 , 7 where lower sperm count could lead to a rise in male infertility, reflecting not only difficulty in conceiving but also a signal of poorer general health. The cause of this decline is likely multifactorial, including both prenatal and adult factors, 8 − 10 with environmental factors playing a major role. 4 However, our current knowledge of the determinants of male fertility does not sufficiently explain this constant and global decline. Over the past decades, several environmental pollutants have been found to be male reproductive toxicants, such as airborne pollutants, certain types of persistent organic pollutants, and chemicals commonly found in plastics. However, these known environmental reproductive toxicants might not be the primary drivers of this alarming trend because the overall trends of these specific pollutants are declining across the same time frame as the observed decline in sperm count and many observed effect sizes are relatively small. Thus, it is likely that there are other unidentified reproductive toxicants or uncharacterized interactions between chemicals that may impact male fertility. The emerging field of exposomics offers a comprehensive approach for environmental health research. The concept of the “exposome” was first introduced by C.P. Wild in 2005 as a way to represent the totality of environmental (i.e., nongenetic) drivers of health and disease. 11 The exposome is a function of external forces and internal biological processes. 12 In practice, exposome studies attempt to capture a large set of environmental exposures simultaneously and have the potential to address limitations of existing studies by investigating the impact of “real-life” exposures and their combinations and interactions between exposures. Historically, there have been several major challenges to identifying previously unknown reproductive toxicants. First, populations around the world are exposed to thousands of chemicals and other environmental exposures with large temporal and geospatial variations where most chemicals are introduced without rigorous reproductive toxicity data. This presents an obvious biomonitoring challenge as it is difficult to assess thousands of toxicant exposures and their metabolites simultaneously, particularly chemicals without reproductive toxicity data or those that have never been profiled in humans. Second, not every environmental pollutant has accessible assays, creating additional barriers to investigators and biomonitoring efforts. For pollutants that can be measured, the analytical platforms used are usually highly specific and lack the ability for high-throughput simultaneous assessment of many exposures from a reasonable amount of biospecimen. This restricts investigators to only a small subset of putative or known toxicants with accessible assays, impeding efforts for nontargeted discovery. This is compounded by a third challenge; human exposures are typically measured using available and minimally invasive biospecimen such as urine, blood, and hair. This can create exposure misclassification and spurious findings as male reproductive organs and processes are protected by the blood-testis-barrier, which can lead to different exposure profiles between systemic circulation and the male reproductive system. 13 , 14 Fourth, the standard environment-wide (i.e., exposome) association studies (ExWAS) approach can be limited by multiple testing penalties, leading to high sample size requirements. 15 While ExWAS remains useful, additional strategies are necessary to analyze high-dimensional exposure data. Together, these constraints collectively hinder investigators from more efficiently and comprehensively identifying environmental pollutant influences on male reproductive health. The objectives of this study were to (1) profile a large set of organic pollutants relevant to male reproductive tissues using a targeted approach and a discovery-driven nontargeted analysis (NTA) and (2) assess associations between the detected pollutants and percent motility, concentration, and total motile sperm to discover previously uncharacterized male reproductive environmental toxicants. We leveraged state-of-the-art analytical advances in high-resolution mass spectrometry 16 to simultaneously profile thousands of potential environmental chemicals in seminal plasma, which is more proximal and relevant for male reproductive health compared to measures of chemicals in systemic circulation. 17 , 18 We then combined a novel machine learning pattern recognition approach, principal component pursuit (PCP), 19 , 20 with modern statistical mixtures analyses 21 to efficiently detect associations of environmental chemicals with male reproductive health. Typical studies model one feature (e.g., genetic polymorphism or environmental exposure) at a time, 22 repeated through all features, which incurs severe multiple testing penalties on statistical power. In addition, environmental exposures do not occur in isolation as represented in these models. In real life, individuals are subjected to complex sets of environmental exposures that jointly act to induce endogenous responses. Our novel approach helps address both limitations of traditional statistical approaches by removing noise, reducing the number of statistical comparisons, and allowing for discovery of interactions and exposure patterns.