Cross-Sectional Estimation of Endogenous Biomarker Associations with Prenatal Phenols, Phthalates, Metals, and Polycyclic Aromatic Hydrocarbons in Single-Pollutant and Mixtures Analysis Approaches

Background: Humans are exposed to mixtures of toxicants that can impact several biological pathways. We investigated the associations between multiple classes of toxicants and an extensive panel of biomarkers indicative of lipid metabolism, inflammation, oxidative stress, and angiogenesis. Methods: We conducted a cross-sectional study of 173 participants (median 26 wk gestation) from the LIFECODES birth cohort. We measured exposure analytes of multiple toxicant classes [metals, phthalates, phenols, and polycyclic aromatic hydrocarbons (PAHs)] in urine samples. We also measured endogenous biomarkers (eicosanoids, cytokines, angiogenic markers, and oxidative stress markers) in either plasma or urine. We estimated pair-wise associations between exposure analytes and endogenous biomarkers using multiple linear regression after adjusting for covariates. We used adaptive elastic net regression, hierarchical Bayesian kernel machine regression, and sparse-group LASSO regression to evaluate toxicant mixtures associated with individual endogenous biomarkers. Results: After false-discovery adjustment (q<0.2), single-pollutant models yielded 19 endogenous biomarker signals associated with phthalates, 13 with phenols, 17 with PAHs, and 18 with trace metals. Notably, adaptive elastic net revealed that phthalate metabolites were selected for several positive signals with the cyclooxygenase (n=7), cytochrome p450 (n=7), and lipoxygenase (n=8) pathways. Conversely, the toxicant classes that exhibited the greatest number of negative signals overall in adaptive elastic net were phenols (n=20) and metals (n=21). Discussion: This study characterizes cross-sectional endogenous biomarker signatures associated with individual and mixtures of prenatal toxicant exposures. These results can help inform the prioritization of specific pairs or clusters of endogenous biomarkers and exposure analytes for investigating health outcomes. https://doi.org/10.1289/EHP7396


Introduction
Observational studies are leveraging increasingly high-dimensional data sets, both in terms of data collected on exogenous environmental toxicants and endogenous biomarkers. Historically, studies have characterized single-pollutant associations between individual chemicals or select classes of toxicants [e.g., phthalates, phenols, parabens, polycyclic aromatic hydrocarbons, or metals], health outcomes, and biomarkers of intermediate effects (e.g., biomarkers of inflammation, oxidative stress, or lipid metabolism) (Dominici et al. 2010;Patel 2017). The advantages of single-pollutant analyses are that results are directly applicable for policy intervention through environmental regulations. This is particularly relevant for toxicants that do not share similar point sources of pollution. For example, phthalate exposure occurs commonly through the manufacture of consumer products (Schettler 2006), and although they may have biologically interactive effects with human exposure to PAHs, the latter class of chemicals enter our environment predominantly through combustion and air contamination (Alegbeleye et al. 2017).
Although there is utility in evaluating individual associations and interpreting their relationship within a narrow physiological scope, there is also a critical need to conduct exhaustive estimation of exposure-wide associations with endogenous biomarkers simultaneously. From a policy perspective, this approach can inform risk assessment by allowing robust and efficient identification of multiple biological targets associated with exposures to various classes of chemicals (Drakvik et al. 2020). Multipollutant approaches also allow for the evaluation of the potential additive or synergistic effects that result from interactions between chemicals (Stafoggia et al. 2017;Sun et al. 2013). Furthermore, in a complex mixture of toxicants, regularized regression approaches such as adaptive elastic net can be used to identify the most predictive exposure analytes in association with specific endogenous biomarkers (Zou and Zhang 2009). This is particularly useful for regulatory prioritization and precisely characterizing high-dimensional biomarkers signals associated with environmental contaminants.
Inflammation and oxidative stress are highly discussed mechanisms of action for several classes of environmental toxicants. On a cellular level, toxicants can impact systemic inflammation and oxidative stress through interference with transcriptional regulation, enzymatic activity, intracellular balance of reactive oxygen species, or membrane permeability (Ferguson and Chin 2017;Thompson et al. 2015). For example, exposures to phthalates and parabens have been shown in vitro to inhibit cytochrome p450 activity (Ozaki et al. 2016). In vitro studies of PAHs have shown that exposure can lead to up-regulation of cyclooxygenase activity (Bauer et al. 2018;Siegrist et al. 2019). In addition, several mechanistic studies have reported the elevation of inflammatory biomarkers, such as cytokines, in response to trace metals exposures (Milnerowicz et al. 2015). Determining associations between multiple environmental toxicant classes and biomarkers of various biological pathways in human studies is critical for developing precise risk estimation for health outcomes.
The developmental origins of health and disease hypothesis proposes that early life environmental factors can impact health status throughout the life span (Preston et al. 2018). Therefore, the characterization of associations between environmental toxicants and endogenous biomarkers during pregnancy has wide applications for future studies seeking to link toxicological pathways as mediators for several maternal and child health outcomes. Systemic perturbations in maternal inflammation, oxidative stress, and lipid metabolism during pregnancy can result in tissue damage, altered vascularization, and interference of fetal nutrient supply (Colella et al. 2018;Ferguson and Chin 2017). The consequences of these toxicological mechanisms include altered fetal development and increased risk for neonatal mortality and morbidities (Colella et al. 2018;Jiang et al. 2018).
In the present study, our primary goal was to characterize an exposure-wide association assessment of a large panel of endogenous biomarkers. By doing so, we estimated single-pollutant associations between individual toxicants and endogenous biomarkers. Building upon this, our secondary goal was to evaluate the exposure variable importance for each endogenous biomarker in the context of a) mixtures of individual toxicants and b) whole classes of toxicants within a grouped multipollutant mixtures framework. We hypothesized that toxicants would have unique profiles of associations with endogenous biomarkers and that toxicant class grouping structure would drive these profiles.

Study Population
Participants in this cross-sectional study were a subset of the LIFECODES pregnancy cohort recruited at the Brigham andWomen's Hospital in Boston, Massachusetts, between 2006 and2008. Recruitment occurred at <15 wk gestation, and study participants were >18 years of age. Questionnaires were administered at recruitment to collect for demographic and health information. Urine and blood specimens were collected at a clinic visit occurring between 23.1 and 28.9 wk (median 26 wk) gestation. From the overall cohort of 1,600 pregnant women, a sample of 173 participants was available for the present study. The sample was selected based on prioritizing participants with availability of biological samples (plasma and urine) where the greatest number of exposure analytes were measured. Secondary to this was prioritizing the proportion of cases and controls (1:2) as close as possible to that in the overall LIFECODES cohort. The larger LIFECODES cohort had 27% cases, whereas our study sample had 33.5%, which included 58 women who delivered preterm (<37 wk gestation) and 115 randomly selected women who delivered after 37 wk gestation (Table S1). The Brigham and Women's Hospital administered institutional review board approval for this study. Detailed descriptions regarding recruitment and study design for the LIFECODES cohort have been published previously (Ferguson et al. 2014;McElrath et al. 2012).

Exposure Biomarkers
Aliquots of urine samples were analyzed for exposure biomarkers at NSF International (Ann Arbor, MI), and extensive details on protocols have been previously described (Aung et al. 2019a;Ferguson et al. 2015bKim et al. 2018). Quantification of exposure analytes were based on protocols and methods developed by the Centers for Disease Control and Prevention for use in the National Health and Nutrition Examination Survey. We measured 17 trace metals using a Thermo Fisher iCAP RQ inductively coupled plasma mass spectrometer with a Teledyne CETAC Technologies ASX-520 autosampler. We used isotope dilution liquid chromatography with tandem mass spectrometry (ID-LC-MS/MS) to quantify 8 PAH metabolites, 7 phenolderived compounds, and 4 parabens. Finally, we quantified 9 phthalate metabolites using high-performance LC-electrospray ionization-MS/MS (HPLC-ESI-MS/MS). Figure 1 shows each exposure analyte according to its toxicant subclass and indicates abbreviations for each analyte. All exposure analytes that were below the limit of detection (LOD) were imputed by dividing the LOD value by the square root of 2.

Endogenous Biomarkers
Both urine and plasma samples were used to measure endogenous biomarkers (see Figure 1 for details and abbreviations for each analyte and biomarker). We measured inflammatory markers in plasma samples at the University of Michigan Cancer Center Immunology Core, including four cytokines using the Milliplex Multiplex Assay Simultaneous High Sensitivity Human Cytokine Magnetic Bead Panel (EMD Millipore Corp.) and C-reactive protein using a DuoSet enzyme-linked immunosorbent assay (R&D Systems). We quantified plasma concentrations of the angiogenic biomarkers PGF and sFlt-1 using the ARCHITECT immunoassay (Abbott Laboratories). In addition, we measured a panel of 53 eicosanoids in plasma using a 6490 triple quadrupole mass spectrometer (Agilent). Three unique protein damage markers NY, DY, and CY) were measured in plasma samples using ESI-MS/MS. Finally, two oxidative stress markers were measured in urine samples at Cayman Chemical: 8-IP, which was quantified via affinity column chromatography and enzyme immunoassay, and 8-OHdG, which was quantified with direct dilution and enzyme immunoassay. Endogenous biomarkers that were below the LOD were imputed with the LOD value divided by the square root of 2. More extensive details on analysis and measurement of endogenous biomarkers have been previously described (Aung et al. 2019b;).

Descriptive Statistics
We tabulated all potential covariates for the study sample, including maternal age, race (White, African American, or other), education level (high school graduate, technical school, some college, or college graduate), health insurance provider (public or private), body mass index (BMI) at the initial study visit (<25 kg=m 2 , 25-29:9 kg=m 2 , or ≥30 kg=m 2 ), alcohol use during pregnancy (yes or no), tobacco use during pregnancy (yes or no), and fetal sex (male or female). Next, we estimated the distributions of individual exposure analytes and endogenous biomarkers. Specifically, for descriptive univariate statistics, distribution of exposure analytes (phthalates, phenols, parabens, PAHs, and metals) and endogenous biomarkers (8-IP and 8-OHdG) measured in urine samples were corrected using specific gravity as follows: where UB SG is the specific gravity-adjusted urinary biomarker concentration, and UB is the observed uncorrected urinary biomarker concentration. The constant SG Median (1.015) is the specific gravity population median, and SG is the observed specific gravity of the individual urine sample (Meeker et al. 2009).

Single-Pollutant Regression
Our first study aim was to estimate individual pair-wise relationships between exposure analytes and endogenous biomarkers using single-pollutant analysis. We implemented multiple linear regression to test for these associations, and the regression model was specified as follows: where M denotes the endogenous biomarker (q = 65 biomarkers), and A represents the exposure analyte (p = 38 analytes), both of which were transformed using the natural log. To account for urinary dilution for all urinary exposure and endogenous biomarkers in regression models, specific gravity was modeled as a covariate, rather than using specific gravity-corrected exposure concentrations. From the potential aforementioned covariates, we selected a subset based on a priori assessment from previous studies in the LIFECODES cohort (Aung et al. 2019a(Aung et al. , 2019cFerguson et al. 2014). The final covariate matrix Z contains specific gravity, maternal age, race, education level, health insurance provider, and BMI at initial study visit. Maternal age and specific gravity variables were modeled as continuous and the rest of the covariates were modeled as categorical. Models were fit on complete cases of exposure and endogenous biomarker pairs; therefore any missing values in the exposure, endogenous biomarker, or covariates were excluded. Figure S1 illustrates a directed acyclic graph of the statistical model. To account for the nested case-control sampling, we calculated and applied inverse probability weighting to the model in Equation 2, to weight the cases and controls proportional to the rates in the parent LIFECODES cohort. As a sensitivity analysis, we also modeled preterm birth case-control status as a covariate rather than using inverse probability weights. Based on all possible pairs of comparisons, there were a total of 2,535 unique linear regression models. To account for multiple testing and to control the false-discovery rate (FDR), we applied the Benjamini-Hochberg procedure to the set of p-values obtained from Equation 2, respectively, with FDR being 0.2.

Multipollutant Mixtures Analysis
Our study consisted of two mixtures analysis goals: a) to conduct an exposure-wide mixtures analysis to identify the most predictive exposures among individual analytes and b) evaluate the variable importance of whole toxicant classes. For the first goal, we fit a multiple exposure model of each endogenous biomarker, specified as follows: where M and Z are defined in the same way as in Equation 2. M is now an n × p matrix, where n = 160 and p = 38. We applied adaptive elastic net (Yang and Zou 2013;Zou and Zhang 2009), a regularized regression approach, to each of the endogenous biomarkers separately using the R package gcdnet (version 1.0.5; https://www.rdocumentation.org/packages/gcdnet/versions/1.0.5). The estimates of the vector b = ðb 0 , b A , b z Þ T were obtained by solving where hðm, a, z; bÞ is an assumed loss function for Equation 2, b w = ðb w 1 , Á Á Á , b w p Þ T is a vector of adaptive weights obtained from a regular elastic net, and j Á j takes the absolute value of each component in b A . In addition, the coefficients of the covariates are not penalized in this setup. The elastic net uses both least absolute shrinkage and selection operator (LASSO) (Tibshirani 1996) and ridge (Hoerl and Kennard 1970) type penalties to perform variable selection with a correlated set of exposure analytes. We optimized for the ' 2 penalty (k 2 ) parameter by considering potential values on a grid from 0 to 2, in increments of 0.1. We forced the covariates Z to always be included in the model. For each potential k 2 value, the elastic net function computes the solutions for a fine grid of k 1 values, which is the ' 1 regularization parameter and assumes a Huberized squared hinge loss. From the potential k 1 and k 2 values, we selected the combination of values that yielded the smallest 10-fold cross-validated misclassification error and used their corresponding coefficient estimates as adaptive weights.
In our second mixtures analysis goal to evaluate the variable importance of the whole toxicant classes, we applied sparsegroup LASSO. This method assumes a linear predictor matrix framework and provides evidence of toxicant class importance in association with endogenous biomarkers. Coefficients within the sparse-group LASSO framework were estimated using where LðbÞ represents the likelihood function defined in  (Figure 1). Among the groups, p k represents the number of exposure analytes in the kth group. We applied the mixing parameter a and tuning parameter k, which collectively controls both group-wise sparsity and within-group sparsity of non-null predictor variables. Specifically, we designated a as 0.95, a value that assumes strong overall sparsity while encouraging grouping (Simon et al. 2013). As a sensitivity analysis, we also explored a more conservative a value of 0.99. Based on this implementation, predictor groups with the greatest importance relative to individual endogenous biomarkers are contained in the model as the tuning parameter k increases. Therefore, the groups exiting the model last are the most important.
We also conducted hierarchical Bayesian kernel machine regression (BKMR) to assess class-level importance and selection using the bkmr package (version 0.2.0) (Bobb et al. 2015). In contrast to sparse-group LASSO, hierarchical BKMR accounts for nonlinearity in the exposure matrix. The exposure analytes are partitioned into groups S g , g = 1, Á Á Á , G, corresponding to G = 5 toxicant classes (metals, PAHs, parabens, phenols, and phthalates) (Figure 1). Group-level posterior inclusion probabilities were estimated from hierarchical BKMR to signify the relative importance of the toxicant classes. For each endogenous biomarker, the hierarchical BKMR framework was modeled as where hðÁÞ is a kernel function that may incorporate nonlinear associations and/or exposure interactions, and e ∼ MVNð0, r 2 IÞ. Hierarchical BKMR incrementally allows individual group-level predictors to enter the kernel function at a given time. The relationship between hðS 1 , Á Á Á ,S G Þ and each outcome variable M in the model is further adjusted by a covariate matrix Z. Implementation of hierarchical BKMR relied on a Gaussian kernel distribution, flat priors, and 10,000 iterations for the Markov chain Monte Carlo.
The PAH metabolites had 17 FDR-adjusted associations with endogenous biomarkers, including 13 positive associations (Figure 2A and Excel Table S3). Four PAH metabolites (2-FLU, 2-and 3-PHE, 1-PHE, and 1-PYR) were positively associated with the oxidative stress marker 8-IP. Figure 2. (A) Heat map of pair-wise associations between exposure analytes and endogenous biomarkers estimated using multiple linear regression and inverse probability weights. Models adjusted for maternal age, race, education, health insurance provider, body mass index at first visit, and specific gravity. The sample size for most models was N = 160, except for the following biomarkers: BCPGE1 (N = 159), LTE4 (N = 132), sFlt-1 (N = 156), and PGF (N = 157). Black and blue grids indicate positive and negative associations, respectively. Color intensities are representative of p-values, that is, darker grids indicate smaller p-values. Pair-wise associations that remain significant after controlling the false-discovery rate at 0.2 are labeled by white X symbols. Dashed lines delineate biomarker subgroups and toxicant classes. See Excel Table S3 for the corresponding numeric data. (B) Heat map of multipollutant associations between exposure analytes and endogenous biomarkers selected by adaptive elastic net. Models adjusted for maternal age, race, education, health insurance provider, body mass index at first visit, specific gravity, and preterm birth case status. The sample size for most models was N = 160, except for the following biomarkers: BCPGE1 (N = 159), LTE4 (N = 132), sFlt-1 (N = 156), and PGF (N = 157). Black and blue grids indicate positive and negative signals, respectively. Dashed lines delineate biomarker subgroups and toxicant classes. See Excel Table S4 for the corresponding numeric data. Abbreviations and subclasses of toxicants are defined in Figure 1.

A
Trace metals had 18 FDR-adjusted associations with endogenous biomarkers, and half were positive (n = 9) (Figure 2A and Excel Table S3). However, in the lipoxygenase pathway, all three FDR-adjusted associations were negative (Ba and RVD1; Ba and 13-oxoODE; and Ni and 13-oxoODE).
The directions of association were largely consistent when we adjusted for preterm birth case-control status as a covariate rather than using inverse probability weights ( Figure S2). However, precision was reduced relative to the primary models, and fewer associations met the FDR threshold.

Multipollutant Associations
When adaptive elastic net was used as a multipollutant penalized regression approach, phthalates were selected as predominantly positive signals among eicosanoids in the cyclooxygenase (n = 7), cytochrome p450 (n = 7), and lipoxygenase pathway (n = 8) groups, as well as parent lipid compounds (n = 3) ( Figure 2B and Excel Table S4). PAHs were also predominantly selected as positive signals among eicosanoids in the cytochrome p450 pathway (n = 8). The chlorinated phenols (2,4-DCP, 2,5-DCP, and TCS) were predominantly selected as negative signals among eicosanoids belonging to each of the enzymatic pathways: cyclooxygenase (n = 8), cytochrome p450 (n = 5), and lipoxygenase (n = 4). The selection trend for metals was largely negative among eicosanoids in the lipoxygenase pathway (n = 5) while exhibiting both negative (n = 6) and positive (n = 7) selection among eicosanoids in the cytochrome p450 pathway. Among the endogenous biomarkers, the oxidative stress marker 8-OHdG exhibited the greatest consistency of positive selection across multiple toxicant classes, including phthalates, phenols, PAHs, and metals.
Adaptive elastic net coefficient estimation accounts for confounding and co-adjustment of multiple pollutants; therefore, estimates differ from single-pollutant regression coefficients. For example, for 9,10-DiHOME, one of the most highly selected cytochrome p450 metabolites, adaptive elastic net regression with adjustment for all exposures selected MEOHP   Tables  S3 and S4).
Sparse-group LASSO was used to rank the importance of each toxicant class based on the number of times it appeared with increasing k values. Using an a value of 0.95, phthalates were . Heat map of posterior inclusion probabilities (PIPs) for toxicant classes estimated using Bayesian kernel machine regression (BKMR). Models adjusted for maternal age, race, education, health insurance provider, body mass index at first visit, specific gravity, and preterm birth case status. The sample size for most models was N = 160, except for the following biomarkers: BCPGE1 (N = 159), LTE4 (N = 132), sFlt-1 (N = 156), and PGF (N = 157). Directions of association are estimated directly from joint-wise associations in the BKMR models. Black cells designate pairs that were positively associated, and blue cells designate pairs that were negatively associated. The color intensities are indicative of PIP, that is, darker colors represent higher PIPs. Dashed lines delineate biomarker subgroups and toxicant classes. See Excel Table S5 for the corresponding numeric data. Abbreviations and subclasses of toxicants are defined in Figure 1. selected as the most important class more often than other toxicant classes (47% of the time) (Excel Table S6). Across all of the endogenous biomarkers, other toxicant classes were selected as the most important (or were tied with another class for the most important rank) in the following proportions: phenols (23%), parabens (12%), PAHs (26%), and metals (23%). Phthalate selection as the most important toxicant class was also evident in the context of eicosanoid enzymatic pathways: cyclooxygenase products (47% of the time), cytochrome p450 products (61%), and lipoxygenase products (44% of the time). Sensitivity analysis with a stricter a value of 0.99 yielded similar results, with phthalates being selected as the most important toxicant class 54% of the time (Excel Table S6).

Discussion
The present study characterizes potential intermediate toxicological pathways by implementing an exposure-wide association assessment of multiple classes of toxicants using a unique panel of biomarkers indicative of lipid metabolism, inflammation, oxidative stress, and angiogenesis. We applied a comprehensive analysis of pair-wise single-pollutant relationships and performed shrinkage estimation and variable selection in a multipollutant framework using adaptive elastic net, hierarchical BKMR, and sparse-group LASSO. Overall, this study is a robust assessment of associations across various biological pathways and can be used to develop models of exposure-response relationships and inform toxicological mechanisms of health outcomes.
Across multiple mixtures analysis methods there were consistencies and differences of note. One of the major consistencies was that phthalates exhibited overwhelmingly positive signals with cytochrome p450 and lipoxygenase pathway products based on adaptive elastic net, sparse-group LASSO, and hierarchical BKMR. However, although adaptive elastic net and sparse-group LASSO continued to select phthalates as highly important for the cyclooxygenase pathway, the importance of phthalates in hierarchical BKMR was diminished. Findings also differed with regard to phthalates and the cyclooxygenase metabolite 9-oxoODE, where phthalates were selected as the most important toxicant class by sparse-group LASSO but were not identified as important by adaptive elastic net and hierarchical BKMR. Each of the mixtures analysis methods indicated largely negative signals for phenols with products from the cyclooxygenase, cytochrome p450, and lipoxygenase pathways. However, for the lipoxygenase metabolite 13-oxoODE, the sparse-group LASSO selected phenols as the most important, whereas the hierarchical BKMR identified phthalates as the most important. Both the adaptive elastic net and sparse-group LASSO methods assume linear relationships between exposures and outcomes, whereas BKMR can fit nonlinear exposure-response functions. That may be one factor influencing the differences in findings across methods. In addition, each method has differential tolerance for two important factors in toxicant effects: sparsity of signals and multicollinearity. These nuances in methodological differences underline the importance of interpreting the mixtures results in the context of single-pollutant results, which can help substantiate findings. Ultimately, differences in findings across mixtures methods underscore the need to transparently apply multiple mixtures methods within both frequentist and Bayesian frameworks.
Focusing on the eicosanoid biomarkers, phthalate metabolites were selected as positive signals for multiple products from each of the enzymatic pathways (cytochrome p450, lipoxygenase, and cyclooxygenase). Phthalate metabolites were also positively associated with the parent lipid compounds a-LA and AA, which can be metabolized by cytochrome p450 enzymes. The cytochrome p450 enzymes belong to a superfamily of intrinsic membrane-bound enzymes that are dispersed throughout the endoplasmic reticulum and mitochondria (Guengerich et al. 2016). These enzymes metabolize parent fatty acids into metabolites such as 9,10-DiHOME and 8(9)-EET. Cytochrome p450 enzymes are also involved in the metabolism of xenobiotics (e.g., steroids and environmental chemicals) (Guengerich et al. 2016). In vitro studies have reported evidence that exposure to parent phthalate compounds (dibutyl phthalate and dimethyl phthalate) and metabolites (MEHP and MBP) can inhibit the functions of cytochrome p450 enzymes (Ozaki et al. 2016;Peng et al. 2019) . Our findings illustrate potential biochemical processes that could be sensitive to phthalate exposure, including the metabolism of aLA and other parent compounds into 9,10-DiHOME and 8(9)-EET.
One of the cytochrome p450 metabolites, 9,10-DiHOME, is a leukotoxin and ligand of the peroxisome proliferator-activated receptor gamma (PPAR-c) (Lecka-Czernik et al. 2002). Mechanistic insights from in vitro studies suggest that 9,10-DiHOME can directly alter mitochondrial membrane permeability (Sisemore et al. 2001) and interfere with neutrophil mediated clearance of pathogens (Thompson and Hammock 2007). In addition, 8(9)-EET is also a ligand of PPAR-c and has broad ranging functions including anti-inflammatory and anti-hypertensive effects (Barnych et al. 2017;Bystrom et al. 2011). Interestingly, mechanistic studies have reported evidence that phthalate metabolites are capable of stimulating PPAR-c (Engel et al. 2017;Schlezinger et al. 2004). Inhibition of cytochrome p450 activity may potentially result in increased substrate and metabolite concentrations. Furthermore, phthalate stimulation of PPAR-c can also create competition for cytochrome p450 metabolites, which could partly explain the positive associations we observed with products from this pathway. The conceptual hypothesis of phthalate-induced inhibition of cytochrome p450 and interference of PPAR-c should be tested further in additional replication studies.
PAHs were predominantly selected as positive signals for cytochrome p450 metabolites and negative signals for cyclooxygenase metabolites. PAH exposure can lead to reduced gap junction activity (Bauer et al. 2018), which is integral for epithelial cell-to-cell signaling and immune cell adhesion and migration (Okamoto and Suzuki 2017). Furthermore, PAHs have proinflammatory capacity upon ligand binding of aryl hydrocarbon receptors and interactions with the nuclear factor kappa-lightchain-enhancer of activated B cells (NF-jB), which is a transcription factor that regulates the expression of inflammatory signaling molecules (Gdula-Argasi nska et al. 2016). Mechanistic studies have demonstrated that PAHs can form DNA adducts and influence gene expression of several enzymes, including cytochrome p450 and cyclooxygenases (Bauer et al. 2018;Gdula-Argasi nska et al. 2016). PAHs may potentiate an increase in lipoxygenase products as an indirect result of stimulating enhanced expression of cytochrome p450 and cyclooxygenases.
The selection of phenols and parabens as negative signals was most evident in the lipoxygenase and cyclooxygenase pathways and, to a lesser extent, in the cytochrome p450 pathway. Mechanistic in vitro models have shown that chlorinated phenols can alter cytosolic calcium ion levels (Michałowicz et al. 2018), which can influence calcium-dependent metabolism of fatty acids (Schlottmann et al. 2014). In another example, the chlorinated phenol TCS has been shown to inhibit the aromatase enzyme CYP19A1 in vitro, which can alter production of estradiol (Li et al. 2017). Inhibition of estradiol production may result in decreased lipid metabolism (Naukam and Curtis 2019), which may partially explain the negative signals that we observed for metabolites across the eicosanoid enzymatic groups. Additional replication studies are needed to inform potential mechanisms by which phenols and parabens impact lipid metabolism.
Lipoxygenases, which are expressed throughout epithelial tissue and by circulating immune cells, catalyze the formation of inflammatory mediators from parent fatty acid compounds such as LA and AA (Mashima and Okuyama 2015). The catalytic activity of lipoxygenases is modified by intracellular Ca 2+ concentrations and the process by which lipoxygenases translocate from the cytosol to biomembranes is Ca 2+ -dependent (Rouzer and Kargman 1988;Walther et al. 2004). In addition to phenols, metals (i.e., As, Ba, and Pb) were also selected as negative signals for lipoxygenase products by elastic net regression. Several trace metals have cationic states and can, therefore, potentially alter lipoxygenase activity by competing with Ca 2+ . For example, Pb exposure in animal models has demonstrated inhibition of Ca 2+ -dependent enzymes such as protein kinases, resulting in altered cytosolic concentrations (Kasten-Jolly and Lawrence 2018). Essential metals such as Mn are also potent ligands of lipoxygenases (Zhu and Richards 2017). As such, increasing intracellular concentrations of multiple trace metals may alter the catalytic activity of lipoxygenases.
Oxidative stress in damaged tissues and blood lymphocytes can lead to the oxidation of macromolecules such as mitochondrial and nuclear DNA, proteins, and lipids. These oxidative stress signals can be quantified by measuring oxidation products such as 8-OHdG, NY, DY, CY, and 8-IP. Notably, there were several positive associations with 8-OHdG, representing multiple classes (trace metals, PAHs, phenols, phthalates). After adjusting for FDR, this included MBP, BPA, 2,5-DCP, 2-FLU, 1-PHE 1-PYR, Se, and Zn. On the other hand, we observed much more subtle patterns with the protein damage markers. Similarly, there were also less consistent signals for the lipid oxidative stress marker 8-IP compared with 8-OHdG. These findings suggest that 8-OHdG may be a common biomarker of oxidative stress across several toxicant classes, although the utility of 8-OHdG measured via immunoassay has been challenged in validation studies (Barregard et al. 2013;Cooke et al. 2009;Møller et al. 2012;Rossner et al. 2016). Furthermore, 8-OHdG and 8-IP are urinary outcome biomarkers, and even with adjustment of specific gravity, residual bias due to hydration status may have differentially influenced the signals we observed (O'Brien et al. 2016). Therefore, additional biomarkers of oxidative stress should be measured to complement inferences drawn from 8-OHdG.
In a separate prospective birth cohort of singleton pregnancies based in Michigan (n = 56), specific gravity-adjusted phenols, phthalates, and metals were measured in maternal urine samples in the first trimester (8-14 wk gestation) for cross-sectional and prospective analyses with oxidative stress biomarkers measured in maternal urine samples and neonatal cord blood samples (Puttabyatappa et al. 2020). Consistent with our findings, adjusted linear regression models in the prospective analyses of that study reported a negative association between first-trimester maternal urinary Pb and neonatal cord blood CY (Puttabyatappa et al. 2020). Exposure levels of Pb in this sample [geometric mean ðGMÞ = 0:22 lg=L] were comparable, but slightly lower than levels observed in the present LIFECODES subsample (GM = 0:37 lg=L) (Goodrich et al. 2019). However, their findings for first-trimester urinary phthalate metabolites in adjusted linear regression models were less consistent, where inverse cross-sectional associations were observed between two metabolites (MCPP and MCINP) and the first-trimester maternal plasma CY level (Puttabyatappa et al. 2020). In addition, in prospective analyses, firsttrimester maternal urinary MBzP was inversely associated with neonatal cord blood NY (Puttabyatappa et al. 2020). Puttabyatappa et al. (2020) analyzed phthalates with lower GMs of first-trimester maternal urinary MCPP (1:38 lg=L) and MBzP (4:2 lg=L) compared with what we observed in our LIFECODES subsample for MCPP (GM = 2:64 lg=L) and MBzP (GM = 8:22 lg=L) (Goodrich et al. 2019). The GM for first-trimester MCINP in that comparison study was 1:29 lg=L (Goodrich et al. 2019), although this metabolite was not measured in our LIFECODES subsample. Another prospective study based in Michigan (n = 24) measured plasma BPA and oxidative stress markers (NY, DY, and CY normalized for total tyrosine) in women (8-14 wk gestation) with singleton term pregnancies and in cord blood (Veiga-Lopez et al. 2015). Plasma BPA was operationalized as low (range: 0:03-0:14 ng=mL) and high (range: 4:1-96:4 ng=mL) (Veiga-Lopez et al. 2015).
Cross-sectional crude analyses in that study reported a positive correlation between maternal plasma concentrations of tyrosine normalized NY and unconjugated BPA (Veiga-Lopez et al. 2015). Differences in biological matrices for BPA exposure assessment limit comparisons in this context, but they circumstantially support evidence of associations across tissue types. Inconsistent associations observed across independent studies may be due to several factors, including exposure distributions, sample size differences, and the timing of sample collection during pregnancy. For example, samples collected in the present LIFECODES subsample were between 23.1 and 28.9 wk of gestation, whereas the Michigan based study samples measured biomarkers in biological media collected between 8 and 14 wk of gestation. Beyond these factors, inconsistencies in associations may be present due to spurious findings within our study and comparison studies.
In the larger LIFECODES cohort (N = 464), estimates from linear mixed effects models used to estimate associations between maternal urine 8-OHdG and maternal urinary exposures measured in samples collected at multiple time points (5-38 wk gestation) from each participant (n = 1,555 samples) indicated positive associations for phenols and parabens with 8-IP and 8-OHdG (Ferguson et al. 2019), consistent with the present analysis. A separate cross-sectional study of a subset of 200 LIFECODES participants ) analyzed urine samples collected at median 26 wk gestation and reported associations between maternal urinary PAH metabolites, 8-OHdG, and 8-IP that were in the same direction as in the present study. In a crosssectional study of women (n = 230) who delivered singleton term infants based in Taiwan, Huang et al. (2017) measured urinary BPA, nonylphenol, and 8-OHdG from samples collected between 27 and 38 wk of gestation. That study did not report notable associations between maternal urinary BPA and 8-OHdG concentrations in adjusted linear regression models (n = 200) (Huang et al. 2017). However, their study did find that higher levels of nonylphenol was associated with higher 8-OHdG (Huang et al. 2017). In adjusted linear regression models of another analysis with a subset of the sample (n = 202), Huang et al. (2018) found that pregnant women with urinary BPA and nonylphenol both above the median was associated with higher 8-OHdG. The creatininecorrected GM BPA concentration in third-trimester urine samples from the Taiwanese study population was 2:24 lg=g creatinine, whereas in the present LIFECODES subset sample, the specific gravity-corrected GM urine BPA concentration in samples collected 23-29 wk of gestation was 1:50 lg=L (Excel Table S1).
Associations between phthalate metabolites and urinary 8-OHdG and 8-IP were consistent between the repeated measures analysis of multiple samples from the larger LIFECODES study sample (Ferguson et al. 2015a) and the present subset analysis. In a prospective study (n = 46) based in South Korea, paraben exposure and oxidative stress biomarkers were measured in urine samples collected from mothers a day before delivery (one participant was collected after delivery) (Kang et al. 2013). That study identified suggestive (p = 0:07) positive cross-sectional associations in adjusted linear regression models between maternal urinary EPB (specific gravity-corrected median = 44:6 lg=L) and 8-OHdG (Kang et al. 2013). Although we did not observe similar associations in the present LIFECODES subset analysis, this may partially be explained by exposure levels of EPB, which were substantially lower in the LIFECODES cohort subset (specific gravity-corrected median = 4:34 lg=L).
Among the trace metals, Al-Saleh et al. (2017) conducted a large cross-sectional study (n = 944) of mothers without chronic conditions (e.g., hypertension, renal and cardiac conditions, diabetes) in Saudi Arabia and measured trace metals and 8-OHdG in maternal urine samples collected 3-12 months after their recent pregnancy. Adjusted linear regression models in that study (n Hg = 881, n Cd = 881, n Pb = 875) indicated that maternal urinary Hg (GM = 0:68 lg=g creatinine), Cd (GM = 0:31 lg=g creatinine), and Pb (GM = 3:77 lg=g creatinine) were associated with higher maternal urinary 8-OHdG (Al-Saleh et al. 2017). In correlation estimations from another subset analysis in this sample (n = 316), Al-Saleh et al. (2016) also observed a positive correlation between total mercury in breast milk (GM = 0:77 lg=L) and maternal urinary 8-OHdG. The inconsistent relationships may be partly due to geographic differences or, in the case of Pb, a significant difference in exposure levels (LIFECODES subsample GM = 0:37 lg=L), albeit the use of creatinine in the Saudi Arabian sample limits the ability to make direct exposure distribution comparisons. Additional replication analyses will further inform oxidative stress associations in observational studies.
Collectively, the associations reported in the present study underscore the potential for enzymatic perturbations through exposure to environmental toxicants. The clinical implications can be postulated based on knowledge of the pathophysiology linked to cytochrome p450, lipoxygenases, and cyclooxygenases. Given the immunoresponsive nature of eicosanoid products, alterations in these signaling molecules may influence the risk for cardiovascular conditions such as hypertension and myocardial infarction (Fan et al. 2015;Parente and Perretti 2003;Rowland and Mangoni 2014;Singh and Rao 2019). During pregnancy, cardiovascular disease risk also affects maternal and neonatal mortality and morbidity (Leon et al. 2019;Roos-Hesselink et al. 2019). Changes in these endogenous eicosanoid biomarkers also pose altered risk for chronic kidney and liver disease (Afshinnia et al. 2020;Fan et al. 2015;Rowland and Mangoni 2014). Altogether, multiple health end points are susceptible to environmental perturbations in enzymatic pathways that govern the concentrations of circulating eicosanoids and biomarkers of inflammation and oxidative stress.
There are some limitations in our study that are important to underline in order to improve future studies. First, the study design was cross-sectional, and our findings are therefore susceptible to reverse causation. For example, disease states caused by other stimuli (e.g., infection or injury) can impact systemic inflammation and oxidative stress, which in turn could alter renal clearance of toxicants. In addition, cross-sectional associations are unable to confirm that the exposures preceded the outcome variables. We recommend that future studies explore repeated measurements of the panel of endogenous biomarkers that we focused on in order to assess windows of susceptibility. The associations reported in this study also have the potential for being noncausal due to uncontrolled confounding, selection bias, exposure misclassification, and other sources of bias that cannot be ruled out. We also conducted multiple comparisons; therefore, there is risk for false discoveries. This was partly ameliorated by accounting for FDR in our interpretations. FDR adjustment does not, however, reduce the risk of false discoveries due to bias. For example, the risk of amplification bias in multipollutant models may reduce accuracy in comparability between single and multipollutant models (Weisskopf et al. 2018). In terms of effect estimation, an additional limitation to consider is the potential for differences in the amount of exposure measurement error for various biomarkers of exposure, which underscores the possibility that larger magnitudes of associations do not necessarily confirm stronger true associations. Another limitation that should be highlighted is our use of circulating endogenous biomarkers in plasma, which limited our inferences to systemic perturbations. We may be missing important tissue-specific relationships associated with toxicant exposures, and the relevance of findings to subsequent health outcomes is uncertain. Last, the findings of this study are not necessarily generalizable to other populations.
Our study also contains several strengths that should be considered. First, we conducted our analyses in a well-characterized birth cohort. Drawing from a case-control sample of preterm births, we partly reduced selection bias by constructing inverse probability weights, which enhances the quality of inference in our findings. Furthermore, we leveraged high-sensitivity instrumentation to measure a robust panel of endogenous biomarkers. Similarly, our exposure assessment was also highly comprehensive, covering a wide range of toxicant classes. Because of these strengths, future studies can build upon our model parameterization to evaluate replication across cohorts. In addition, the statistical approach in our study is innovative because it used penalized regression, hierarchical BKMR, and sparse-group LASSO to disentangle pathway specific associations for toxicants within a multipollutant framework. This approach is particularly useful for policy applications in order to help inform prioritization of toxicants based on biological mechanisms.
In conclusion, this study leveraged multiple statistical methods to characterize signals of endogenous biomarkers associated with individual toxicants and whole toxicant classes. When we applied adaptive elastic net, BKMR, and sparse-group LASSO, phthalates were consistently linked to positive signals with plasma eicosanoids from the cytochrome p450 and lipoxygenase pathway. Phenols were also consistently linked to negative signals with all eicosanoid enzymatic pathways. This study informs the prioritization of specific pairs or clusters of endogenous biomarkers and exposure analytes for future studies of health outcomes.