Results
Characteristics of the 492 mother-child dyads are described in Table 1 . Most women in the analytic sample identified as non-Hispanic White (81%), half reported being nulliparous (50%), and few women smoked in the first trimester of pregnancy (5%). Conception season was evenly distributed. The median (25th, 75th percentile) maternal age was 31 years ( 28 , 33 ). The median (25th, 75th percentile) maternal pre-pregnancy BMI was 25 kg/m 2 ( 22 , 29 ). The median AHEI-2010 score across pregnancy (out of 100 points) and early-pregnancy PSS score (out of 40 points) was 52 points and 11 points, respectively. Infant sex was evenly distributed and about half of mothers (46%) arranged childcare at home. Among infants who completed the 5-month visit, 56.0% were exclusively breastfed, 4.2% were exclusively formula fed, and 39.8% received a combination of breastmilk and formula from birth until 5 months of age.
Over 90% of women had nonzero concentrations of all nonpersistent EDC biomarkers, except butylparaben, BPF, and triclocarban (Table S1) ( 63 ). As we previously reported, concentrations of most EDC biomarkers were similar to NHANES, except I-KIDS mothers had lower concentrations of MEP, ethylparaben, methylparaben, and propylparaben, but higher levels of ∑DEHTP, BP-3 and triclosan ( 62 ). Few EDC biomarkers were strongly correlated with each other, except ∑DiNP and MCPP (r = 0.76), methylparaben and propylparaben (r = 0.70), and 2,4-DCP and 2,5-DCP (r = 0.70) ( 61 ).
The distributions of newborn and 5-month AGD and AGI measures are reported in Table S2 ( 63 ). As expected, AGD and AGI at birth and 5 months were larger in males than in females. The % growths in AGD Short and AGI Short from birth to 5 months were only marginally greater in females compared to males ( P = .19 and P = .11, respectively), whereas median % growths in AGD Long and AGI Long were greater in males compared to females (51.56% compared to 32.76%; 19.69% compared to 6.78%; respectively).
In covariate-adjusted linear mixed models, few EDC biomarkers were associated with AGI Short or AGI Long across infancy or at birth in females (Table S3) ( 63 ). At birth, the QGComp phthalate/replacement (but not phenol) mixture was imprecisely associated with 2.03 mm/m (95% CI: −0.31, 4.48) longer AGI Short ( Fig. 1A , Tables S4 and S5) ( 63 ). Similarly, the hierarchical BKMR mixture, driven by phthalates/replacements, was nonlinearly associated with longer AGI Short at birth, with attenuation of the relationship at higher mixture quantiles ( Fig. 2A and Table S6) ( 63 ).
Associations of QGComp mixtures with newborn, 5-month and % AGI growth in A) female infants and B) male infants. QGComp models accounted for race/ethnicity, maternal age, mean pregnancy AHEI-2010, perceived stress, conception season, pre-pregnancy BMI, parity, and childcare arrangements (if applicable). Data are interpreted as the percent change in AGD outcome for each quartile (25%) increase in the QGComp mixture. The phthalate/replacement mixture included MEHP, MEHHP, MEOHP, MECPP, MCOP, MCNP, MCPP, MBzP, MEP, MBP, MHBP, MiBP, MHiBP, MHiNCH, MCOCH, MEHHTP, and MECPTP. The phenol mixture included ethylparaben, methylparaben, propylparaben, BPA, BPS, BP-3, TCS, 2,4-DCP, and 2,5-DCP. 2,4-DCP, 2,4-dichlorophenol; 2,5-dichlorophenol; AGI, anogenital index; AHEI-2010, Alternative Healthy Eating Index 2010; BMI, body mass index; BP-3, benzophenone-3; BPA, bisphenol A; BPS, bisphenol S; MBP, mono-n-butyl phthalate; MBzP, monobenzyl phthalate; MCNP, monocarboxynonyl phthalate; MCOCH, cyclohexane-1,2-dicarboxylic acid-mono(carboxyoctyl) ester; MCOP, monocarboxyoctyl phthalate; MCPP, mono(3-carboxypropyl) phthalate; MECPP, mono(2-ethyl-5-carboxypentyl) phthalate; MECPTP, mono (2-ethyl5-carboxypentyl terephthalate; MEHHP, mono(2-ethyl-5-hydroxyhexyl) phthalate; MEHHTP, mono (2-ethyl-5-hydroxyhexyl terephthalate; MEHP, mono(2-ethylhexyl) phthalate; MEOHP, mono(2-ethyl-5-oxohexyl) phthalate; MEP, monoethyl phthalate; MHBP, mono-hydroxybutyl phthalate; MHiBP, mono-hydroxy-isobutyl phthalate; MHiNCH, cyclohexane-1,2-dicarboxylic acid-monohydroxy isononyl ester; MiBP, mono-isobutyl phthalate; QGComp, quantile g-computation; TCS, triclosan. n = 106–148. * P ≤ .05; # P ≤ .10.
Associations of a nonpersistent EDC biomarker mixture with A) short and B) long newborn, 5-month, and % AGI growth using hierarchical BKMR in female and male infants. Hierarchical BKMR models were fit with 100 000 iterations and accounted for race/ethnicity, maternal age, average pregnancy AHEI-2010, perceived stress, conception season, pre-pregnancy BMI, parity, and childcare arrangements (if applicable). Data are presented as effect estimates and 95% credible intervals, which are interpreted as the association between the mixture at each quantile and AGD compared to when all co-exposures are fixed at the median. The mixture included MEHP, MEHHP, MEOHP, MECPP, MCOP, MCNP, MCPP, MBzP, MEP, MBP, MHBP, MiBP, MHiBP, MHiNCH, MCOCH, MEHHTP, MECPTP, ethylparaben, methylparaben, propylparaben, BPA, BPS, BP-3, TCS, 2,4-DCP, and 2,5-DCP. 2,4-DCP, 2,4-dichlorophenol; 2,5-dichlorophenol; AGI, anogenital index; AHEI-2010, Alternative Healthy Eating Index 2010; BKMR, Bayesian kernel machine regression; BMI, body mass index; BP-3, benzophenone-3; BPA, bisphenol A; BPS, bisphenol S; MBP, mono-n-butyl phthalate; MBzP, monobenzyl phthalate; MCNP, monocarboxynonyl phthalate; MCOCH, cyclohexane-1,2-dicarboxylic acid-mono(carboxyoctyl) ester; MCOP, monocarboxyoctyl phthalate; MCPP, mono(3-carboxypropyl) phthalate; MECPP, mono(2-ethyl-5-carboxypentyl) phthalate; MECPTP, mono (2-ethyl5-carboxypentyl terephthalate; MEHHP, mono(2-ethyl-5-hydroxyhexyl) phthalate; MEHHTP, mono (2-ethyl-5-hydroxyhexyl terephthalate; MEHP, mono(2-ethylhexyl) phthalate; MEOHP, mono(2-ethyl-5-oxohexyl) phthalate; MEP, monoethyl phthalate; MHBP, mono-hydroxybutyl phthalate; MHiBP, mono-hydroxy-isobutyl phthalate; MHiNCH, cyclohexane-1,2-dicarboxylic acid-monohydroxy isononyl ester; MiBP, mono-isobutyl phthalate; TCS, triclosan.
In 5-month-old females, using linear mixed models, most phthalates and replacements were associated with either shorter AGI Short or AGI Long (Table S3) ( 63 ). Consistently, each quartile increase in the QGComp phthalate/replacement (but not phenol) mixture was associated with a 4.13 mm/m (95% CI: −7.22, −1.05) shorter AGI Short , driven by MECPP (16%), MiBP (15%), and MCOCH (15%) ( Fig. 1A , Tables S4 and S5) ( 63 ), and a 3.89 mm/m (95% CI: −7.55, −0.42) shorter AGI Long , driven by MCOCH (26%), MCNP (23%), and MHiBP (22%) ( Fig. 1A , Table S4) ( 63 ). Using hierarchical BKMR, the mixture was linearly associated with shorter AGI Short and was driven by phthalates/replacements (PIP: 0.82), specifically MECPTP (PIP: 0.26) ( Fig. 2A , Table S6) ( 63 ). The BKMR mixture was also linearly associated with shorter AGI Long , but was driven by phenols (PIP: 0.67; primarily methylparaben, PIP: 0.32) ( Fig. 2B , Table S6) ( 63 ).
In multivariable linear regression models, several phthalates/replacements (eg, ∑DiNP and its metabolites; and ΣDEHTP and its metabolites) were associated with less % growth of the female AGI Short or AGI Long from birth to 5 months, whereas the DiNCH metabolite MCOCH was associated with greater % growth of AGI Long (Table S3) ( 63 ). Using QGComp, a quartile increase in the phthalate/replacement mixture (but not the phenol mixture; Fig. 1A , Table S5) ( 63 ) was associated with a 39.49% (95% CI: −62.54, −16.54) less growth of AGI Short from birth to 5 months, driven by MECPTP (25%), MEHHP (17%), and MCOP (13%), ( Fig. 1A , Table S4) ( 63 ) and 8.78% (95%CI: −16.06, −1.51) less growth of AGI Long from birth to 5 months, with MECPTP (23%), MCOCH (21%), and MCNP (12%) driving the observed association ( Fig. 1A , Table S4) ( 63 ). The BKMR mixture was linearly associated with less growth of AGI Short (driven by phthalates/replacements, PIP: 0.94; primarily MCNP, PIP: 0.39) ( Fig. 2A , Table S6) ( 63 ), whereas there was some evidence of a nonlinear association with less growth of AGI Long —with attenuation of the relationship at higher mixture quantiles (driven by phenols, PIP: 0.79; primarily ethylparaben, PIP: 0.32) ( Fig. 2B , Table S6) ( 63 ).
In covariate-adjusted linear mixed models in males, a limited number of EDC biomarkers were associated with AGI across infancy, with inconsistent associations in newborns (Table S7) ( 63 ). Using QGComp, the phthalate/replacement mixture was not associated with newborn AGI, whereas a quartile increase in the phenol mixture was associated with a 4.10 mm/m (95% CI: −7.63, −0.57) shorter newborn AGI Long , with ethylparaben (38%), 2,5-DCP (26%), and BP-3 (24%) contributing to the association ( Fig. 1B , Tables S8 and S9) ( 63 ). Using hierarchical BKMR, there were potentially nonlinear associations of the mixture with shorter AGI Short and AGI Long , driven by phenols (PIPs: 0.59 and 0.69, respectively) ( Fig. 2A and 2B , Table S10) ( 63 ).
In 5-month-old males, using linear mixed models, some nonpersistent EDC biomarkers (eg, ∑DiNCH) were associated with AGI, though few patterns emerged (Table S7) ( 63 ). Neither the QGComp mixtures ( Fig. 1B , Tables S8 and S9) ( 63 ) nor the BKMR mixture ( Fig. 2A and 2B , Table S10) ( 63 ) was associated with AGI in males at 5 months of age.
Using linear regression, few EDC biomarkers were associated with % growth of AGI in males, and findings were inconsistent (Table S7) ( 63 ). Using QGComp, a quartile increase in the phenol mixture (and not phthalate/replacement mixture, Fig. 1B , Table S8) ( 63 ) was imprecisely associated with 1.91% (95% CI: −0.19, 4.01) more growth of AGI Short ( Fig. 1B , Table S9) ( 63 ). This association was potentially strongest (nonlinear) at higher quantiles of the BKMR mixture driven by phenols (PIP: 0.74; primarily 2,5-DCP, PIP: 0.51) ( Fig. 2A , Table S10) ( 63 ).
As reported in Table S11 ( 63 ), our findings from sensitivity analyses focusing on QGComp chemical mixtures did not substantially differ from the main analyses; however, many results from sensitivity analyses had a loss of precision, likely due to decreased sample sizes. For example, the main trending association between the phthalate/replacement QGComp mixture with newborn AGI Short in females was somewhat attenuated in some sensitivity analyses, accompanied by a loss of precision. Similarly, models evaluating associations of the phthalate/replacement QGComp mixture with 5-month AGI Long in females had slight losses in precisions when excluding babies whose mothers took levothyroxine or progesterone and also when excluding babies whose mothers were diagnosed with HDP, likely due to decreased sample sizes in these models. Finally, in males, the previously significant associations between the phenol QGComp mixture and newborn AGI Long and % AGI Short growth also lost some precision in each sensitivity analysis, likely due to decreased sample sizes.
Disclaimer
The findings and conclusions in this report are those of the authors and do not necessarily represent the official position of the CDC. Use of trade names is for identification only and does not imply endorsement by the CDC, the Public Health Service, or the US Department of Health and Human Services.
Discussion
Our findings in a relatively homogeneous sample of healthy pregnant Midwestern U.S. women with high socioeconomic status suggest that in females, in utero exposure to nonpersistent EDC biomarkers, particularly plasticizer phthalates and their replacements, is associated with shorter AGI in mini-puberty and less growth of the AGI from birth to mini-puberty. In male infants, some phenol biomarkers were associated with measures of AGI, but findings were inconsistent and potentially nonlinear. Sensitivity analyses demonstrated the robustness of our main findings; however, additional studies may be needed to further explore if there are differences in associations by other pregnancy complications. Overall, our findings support the need to understand the long-term implications of EDC exposures for reproductive health, particularly in females.
In newborns, we observed few individual phthalate/replacements to be consistently associated with AGI in females or males, but the phthalate/replacement mixture—driven by MCNP and MHiBP—was imprecisely associated with longer AGI Short in females, and the phenol mixture—driven by ethylparaben and 2,5-DCP—was associated with shorter AGI Long in males. To date, several studies have investigated the association of nonpersistent EDC biomarkers with newborn AGD, but many studies only considered males, and those that also evaluated AGD in females have been largely inconsistent. For example, a study in Israel reported that phthalate exposure was not associated with newborn AGD in males, which was similar to our study ( 71 ); however, MBzP was positively correlated with AGD in Israeli newborn females, but MBzP was not associated with female newborn AGD in our study ( 71 ). Similar to our findings, a 2015 study in The Infant Development and the Environment Study (TIDES) cohort also reported no associations of DEHP metabolites with newborn AGD in females ( 78 ). A U.S. study investigating the role of race on prenatal phthalate exposure and newborn AGD reported that a BKMR mixture of MBP, MiBP, MBzP, MEHP, MEOHP, MEHHP, MEP, and monomethyl phthalate was associated with shorter newborn AGD Long in female newborns born to Black women ( 79 ). In our study, the hierarchical BKMR mixture included additional phthalates/replacement and phenols, and therefore is not directly comparable; however, the mixture was nonlinearly associated with newborn AGI Short in females, primarily driven by MCNP and MCOP. Among studies examining environmental phenols, BPA, BPS, and BPF were not associated with newborn AGD in females in Thailand ( 39 ), whereas BPA was associated with shorter AGD Short in newborn females the U.S ( 38 ). Inconsistencies across studies could be due to differing exposure distributions and EDC milieus, timing of EDC exposure assessment, the absence of phthalate or bisphenol replacement chemical metabolites (eg, DEHTP, DiNCH, and BPS) in prior studies, and the limited number of studies considering both male and female newborns.
In our study, many single chemicals, including MCOP, MCNP, MCPP, and metabolites of the phthalate replacement DiNCH, as well as the phthalate/replacements mixture, were negatively associated with female AGI at 5 months. Currently, only a limited number of studies have considered the roles of prenatal nonpersistent EDC exposures with AGD beyond birth, especially during mini-puberty. One notable study in the MIREC cohort reported that maternal MBzP and MEP were associated with shorter and longer AGD Long in 6-month-old female infants, respectively, and MBP was associated with longer AGD Long in 6-month-old male infants ( 68 ). In our study, MBzP was also associated with shorter AGD Long at 5 months in female infants, which is consistent with the MIREC findings. A study in 12-month-old U.S. infants reported a positive, nonlinear association of a phthalate mixture, driven by DEHP and DiNP, with female AGD Short and AGD Long ( 48 ). In our study, the phthalate mixture (driven by the DEHTP metabolite MECPTP) was linearly associated with shorter AGI Short in female infants at 5 months, although we cannot directly compare findings due to the difference in timing of AGD assessment. In the Odense Child Cohort, second-trimester maternal urinary methylparaben was associated with a shorter three-month AGD Short in male infants, but propylparaben was associated with longer AGD Long in both male and female infants ( 41 ), which is unlike our study as we observed no associations between methylparaben or propylparaben with AGI in female infants and a nonstatistically significant association of propylparaben with AGI Long in newborn males. In a separate study in the Odense Child Cohort, prenatal phthalate exposure was not associated with three-month AGD in males ( 73 ), which is similar to the findings in our study. Finally, one study in Chinese infants reported positive associations between prenatal exposure to a BKMR phenol mixture and female AGD at 48 months, but did not observe associations with female AGD at 6 or 12 months of age ( 45 ). This is different from our study in which we reported the BKMR mixture was linearly associated with shorter AGI Long , driven by phenols (primarily methylparaben). It is critical to acknowledge that these studies and ours evaluated infants at different ages, and ours was one of the first to include additional phthalate and phenol replacement like MECPTP and MEHHTP, thereby limiting comparison. Additionally, we used a hierarchical mixture to determine which chemical class was responsible for the observed associations, and therefore included additional chemicals beyond phthalates/replacements, which further complicates comparison between studies. Comparisons were also hindered because we considered and observed possible nonlinear associations of EDCs with AGD outcomes. Due to the inconsistencies in the literature, additional investigations are needed—particularly those measuring AGD during mini-puberty—to better understand the application of AGD as a clinical marker of reproductive system development.
In our study, nonpersistent EDC biomarkers were associated with less % growth of the AGD from birth to 5-months, generally driven by metabolites of the phthalate replacement DEHTP. Several prior epidemiologic studies indicate that there is rapid growth of the AGD between birth and approximately six months of age ( 45-47 ). To that end, three studies have investigated the association of nonpersistent EDC biomarkers with AGD growth trajectories in young children. In one study from the U.K., propylparaben was associated with a slower AGD growth trajectory through 24 months of age in males (but did not consider females) ( 97 ), which differs from our study, as propylparaben was not associated with AGI growth in males. In the other two studies conducted in Chinese infants, BPS and BPF were associated with accelerated AGD growth in females across 48 months, whereas BPA, BPS, BPF, or BPAF were associated with less AGD Long growth from birth to 12 months in males ( 45 , 76 ). In our study, BPA and BPS were not associated with AGI growth in female infants, but BPS was similarly associated with less AGI Long growth in male infants. Although neither study assessed phthalate metabolites, these trajectory studies, like ours, have demonstrated that prenatal exposure to nonpersistent EDC biomarkers is associated with altered AGD growth across infancy. Collectively, these findings suggest that growth of the AGD across infancy, and not just single measurements of AGD at birth or in mini-puberty, may be an important readout of disrupted reproductive system development in utero .
While the precise mechanisms of this observation are not fully elucidated, EDCs can have direct (act on the receptors of developing gonads ( 98 )) and indirect (disrupt placental hormone synthesis ( 99 )) effects on the developing reproductive system. These direct and indirect disruptions in utero can ultimately result in altered levels of sex-steroid hormones and contribute to gynecologic and urologic disorders. For example, several experimental studies have demonstrated that prenatal exposure to parabens and phthalates altered AGD at birth and sex-steroid hormone levels at maturity, and decreased sperm number and motility in male rats ( 100-104 ); findings in female rats are less robust, though paraben and phthalate exposure altered AGD and impaired ovarian function ( 105-107 ). Although there have been no longitudinal studies evaluating associations of prenatal EDC-associated changes in AGD with sex hormone levels or sexual maturity in humans ( 14 ), there are a few studies evaluating associations of prenatal paraben and phthalate exposure with AGD and sex steroid hormones during mini-puberty. A study in the Odense cohort reported that higher prenatal paraben exposure was associated with longer 3-month AGD and lower sex hormone concentrations (follicle stimulating hormone (FSH) and dehydroepiandrosterone sulfate) in females at 3 months of age ( 41 ); conversely prenatal phthalate exposure was associated with a lower testosterone/LH ratio in 3-month-old males ( 108 ). Although studies have not investigated the relationship between AGD and sex-steroid hormones in healthy adult women, a shorter AGD was associated with lower levels of anti-Mullerian hormone and higher levels of FSH among women undergoing in vitro fertilization ( 109 ). Furthermore, while AGD is not standardized for clinical use, several cross-sectional epidemiologic studies have demonstrated that AGD is closely linked to hormone-driven gynecologic disorders: a shorter AGD is associated with a higher risk of being diagnosed with endometriosis ( 16 , 110-113 ), and a longer AGD is associated with a higher risk of a polycystic ovarian syndrome diagnosis ( 16 , 114 ). Though our findings in male infants were less consistent than those in females, it is important to note that cross-sectional studies have reported consistent relationships between shorter AGD and lower testosterone levels among men attending andrology clinics, though there are inconsistencies among healthy men ( 15 , 115-117 ), and studies in adult men have reported associations between longer AGD and better sperm quality, although there was not enough evidence to predict individual fertility ( 118 ). To that end, this body of literature warrants additional studies to further characterize and understand the roles of in utero EDC exposure in AGD growth across infancy, as it may indicate disrupted reproductive system development and predict future reproductive health.
The current study has several limitations, but also many strengths. First, we were unable to assess nonpersistent EDC exposure between birth and the 5-month visit and, therefore, were unable to fully consider the potential role of postnatal EDC exposure (as a precision variable) on AGD growth and development—which may have limited precision of our estimates. However, we accounted for childcare arrangements during the first few months of life, which may account for chemical exposure from the environment and from the diet. Furthermore, our findings are biologically-relevant and AGDs during mini-puberty in our study were generally comparable to, and had overlapping ranges with those previously reported in U.S.- and nonU.S.-based populations during infancy; however, as expected, differences in measurement timing resulted in some inconsistencies ( 41 , 45 , 47 , 65 , 73 , 74 , 119 ). Although we reported our findings as body length-normalized AGI, our QGComp effect estimate for AGD Short (for example) corresponds to approximately 2.6 mm shorter AGD Short , which reflects a quartile decrease in AGD in our population and also in a prior U.K. study reporting trajectories of AGD growth across infancy, thus indicating a potentially biologically-relevant change in AGD ( 47 ). Second, urinary concentrations of a few EDC biomarkers were lower or higher than in prior studies evaluating EDC biomarkers and AGD, making direct comparisons difficult. However, we evaluated an extensive panel of nonpersistent EDC biomarkers from multiple chemical classes, including phthalate alternatives, and quantified nonpersistent EDC biomarkers from three-to-five pooled first-morning urine samples, providing a relatively stable estimate of gestational EDC exposure ( 59 , 120-122 ). Third, our study focused on a relatively small number of nonpersistent EDCs, but pregnant women are exposed to hundreds of unique chemicals daily ( 123 , 124 ). Although some studies have considered other important environmental drivers of AGD, like per- and polyfluoroalkyl substances ( 125-128 ), organophosphate esters ( 129 ), and dioxins ( 130 , 131 ), phthalates and phenols have potentially unique exposure sources, and therefore our findings are likely not confounded by other chemical exposures. Regardless, additional studies that focus on EDCs from other exposure sources are needed to understand the roles of chemical exposures more broadly in reproductive system development. Fourth, we observed that some EDC biomarkers were meaningful contributors to the QGComp mixture but not the BKMR mixture, and vice versa, which requires careful interpretation. However, our findings suggest this is likely because associations of EDC biomarkers with AGI are potentially nonlinear, highlighting an important strength of our study. Fifth, we cannot establish causality due to the observational nature of the study; however, the reported associations are temporal, with our exposures measured before our outcomes, and our hypothesis is biologically plausible. Sixth, although newborn and infant AGDs are difficult to measure and are prone to measurement error, trained researchers measured each outcome in triplicate, thereby improving precision. Additionally, we normalized AGD to child linear growth to improve precision and interpretability, and also considered AGD growth across infancy, which provides a novel readout of prenatal EDC exposure. Lastly, there is the potential for unmeasured confounding in this observational epidemiologic study, and the enrollment of a relatively healthy cohort with very little pathology may reduce generalizability. Additionally, there are other lifestyle factors like maternal stress ( 69 , 132 ) and pre-pregnancy BMI ( 51 ) that have been previously associated with AGD in newborns; unfortunately, in the current study, we were underpowered to further investigate differences in associations by these variables. However, we used a DAG and the current literature to select appropriate covariates for our models (including stress and BMI), and research in a relatively homogeneous sample of healthy women may reduce unmeasured confounding.
Conclusions
The results from our study suggest that exposure to plasticizer chemicals during pregnancy may adversely impact female reproductive programming that may be particularly salient during mini-puberty. Although it is not clear if early-life AGD follows the same trajectory into adulthood, several epidemiologic studies in adult women in France, Spain, and China have reported strong relations between a shorter AGD and higher risk of endometriosis and adenomyosis ( 110-113 ). Because some of our results differ from the prior literature and given the paucity of studies investigating prenatal EDC exposure with AGD and its growth beyond birth, additional studies in diverse populations with higher EDC exposure are needed to corroborate our findings. Additionally, while substantially more research is needed to link our findings to long-term reproductive health, this study further underscores the need for public health policies aimed at limiting environmental EDCs to support reproductive system development and health across the lifespan.
Materials|Methods
The current study leveraged data from pregnant women enrolled in the Illinois Kids Development Study (I-KIDS), an ongoing prospective pregnancy and birth cohort investigating the impact of prenatal environmental chemical exposures on early offspring neurodevelopment. Participant enrollment and recruitment have been extensively described elsewhere ( 18 , 50-52 ). Briefly, pregnant women in Champaign-Urbana, IL were recruited at their first prenatal appointment from two obstetric clinics. Women were ineligible to participate if they participated in I-KIDS with a prior pregnancy, were more than 15 weeks pregnant, were younger than 18 years of age or older than 40 years of age, not fluent in English, in a high-risk (as determined by their doctor) or multiple pregnancy, living farther than a 30-minute drive to the University of Illinois Urbana-Champaign campus, or moving out of the area before the child's first birthday. All women provided written informed consent to participate in I-KIDS, which was approved by the Institutional Review Board at the University of Illinois Urbana-Champaign. Analysis of de-identified specimens at the Centers for Disease Control and Prevention (CDC) Division of Laboratory Sciences was ruled not to constitute human subjects' research.
I-KIDS researchers conducted a home visit at a median of 13 weeks of gestation to collect extensive information on women's sociodemographic, lifestyle, and health characteristics, including their race/ethnicity, annual household income, parity, smoking status in the first trimester of pregnancy, and age via interviewer-administered questionnaires. Participants self-reported height and pre-pregnancy weight, which we used to calculate pre-pregnancy body mass index (BMI; kg/m 2 ). We calculated conception season using reported due dates based on the first day of the last menstrual period and confirmed after the first-trimester ultrasound ( 52 ). Women also reported their early pregnancy perceived stress using the Perceived Stress Scale (PSS), a self-administered ten-item survey validated for use in pregnant populations that assesses thoughts and feelings experienced during the past month ( 53 , 54 ). Participants additionally completed a three-month semi-quantitative food frequency questionnaire (FFQ) adapted for pregnancy from the full-length Block-98 FFQ (NutritionQuest, Berkeley, CA) ( 55-57 ). We calculated the Alternative Healthy Eating Index 2010 (AHEI-2010) in early and late gestation (median 13 and 35 weeks, respectively) using data from the FFQ. The AHEI-2010 is an 11-component diet quality index (maximum of 110 points) based on foods and nutrients predictive of chronic disease risk; a higher score reflects better diet quality ( 58 ). Given that the AHEI-2010 considers moderate alcohol consumption to be beneficial, but guidelines recommend that pregnant women avoid alcohol, we removed the alcohol component to establish a ten-component diet quality index (maximum of 100 points) and calculated the mean of early and late pregnancy AHEI-2010 scores to approximate diet quality across pregnancy.
I-KIDS mothers and their infants completed a postnatal follow-up visit at a median 5 (25th, 75th percentile: 4.7, 5.1) months, where researchers collected updated demographic and pertinent childcare arrangement information. Specifically, women reported childcare arrangements by answering “yes” or “no” to the questions “Since our last interview, has your child been cared for outside the home (other than school)?” and “Since our last interview; has your child been cared for at home by someone other than parents?”.
Due to the relatively short biological half-lives (6-24 hours) and high within-person variability of nonpersistent EDC biomarkers ( 59 ), we quantified EDC biomarkers in three-to-five across-pregnancy first-morning urine samples that were physically pooled before EDC biomarker measurement as previously described ( 60 , 61 ). Briefly, women collected and refrigerated urine samples at home in polypropylene urine cups the morning of their research visit or clinic visit, which we transported to the study laboratory on the same day. Within 24 hours of collection, we layered 900 μL of urine from each timepoint to a frozen sample from previous timepoints to create the urine pool. Specific gravity of pooled samples was measured using a handheld refractometer (TS400; Reichert Technologies, Depew, NY) at the end of pregnancy, when each pooled sample was thawed and vortexed. All urine was stored at −80 °C until it was shipped to the CDC Division of Laboratory Sciences on dry ice in four batches in the order of participant enrollment (batch 1 enrolled December 2013—February 2015; batch 2 enrolled February 2015—July 2016; batch 3 enrolled July 2016—August 2018; batch 4 enrolled September 2018—August 2019). Urinary nonpersistent EDC biomarkers were quantified at the CDC using online solid phase extraction coupled with isotope dilution-high performance liquid chromatography-electrospray ionization-tandem mass spectrometry as previously described and reported ( 62 ).
The following nonpersistent EDC biomarkers were quantified in all batches as reported in Table S1 ( 63 ): mono(2-ethylhexyl) phthalate (MEHP), mono(2-ethyl-5-hydroxyhexyl) phthalate (MEHHP), mono(2-ethyl-5-oxohexyl) phthalate (MEOHP), mono(2-ethyl-5-carboxypentyl) phthalate (MECPP), monocarboxyoctyl phthalate (MCOP), monocarboxynonyl phthalate (MCNP), mono(3-carboxypropyl) phthalate (MCPP), monobenzyl phthalate (MBzP), mono-n-butyl phthalate (MBP), monohydroxybutyl phthalate (MHBP), mono-isobutyl phthalate (MiBP), monohydroxy-isobutyl phthalate (MHiBP), monoethyl phthalate (MEP), cyclohexane-1,2-dicarboxylic acid-mono(carboxyoctyl) ester (MCOCH), and cyclohexane-1,2-dicarboxylic acid-monohydroxy isononyl ester (MHiNCH). In addition, the CDC measured concentrations of triclocarban, butylparaben, ethylparaben, methylparaben, propylparaben, bisphenol A (BPA), bisphenol F (BPF), bisphenol S (BPS), triclosan (TCS), benzophenone-3 (BP-3), 2,4-dichlorphenol (2,4-DCP), 2,5-dichlorophenol (2,5-DCP). One phthalate, mono-isononyl phthalate (MiNP), was removed from the CDC analytical panel for women in batch four and therefore was only quantified in the first three batches. Three additional phthalate and phthalate replacement metabolites were added to the CDC analytical panel for women in batches two, three, and four: monooxononyl phthalate (MONP), mono(2-ethyl-5-hydroxyhexyl) terephthalate (MEHHTP), and mono(2-ethyl-5-carboxypentyl) terephthalate (MECPTP) (Table S1) ( 63 ).
At the hospital research visit within 24 hours of birth and the postnatal follow-up visit at 5 months of age, trained researchers used digital calipers (Mitutoyo Model 500-195-20) to measure AGD in triplicate using previously published protocols ( 64 ). In females, the AGD Short represents the distance from the center of the anus to the clitoral hood, whereas the AGD Long represents the distance from the center of the anus to the posterior fourchette. In males, the AGD Short represents the distance from the center of the anus to the cephalad insertion of the penis, whereas the AGD Long represents the distance from the center of the anus to the base of the scrotum. Additionally, we measured newborn and infant body length (cm) in triplicate using a measuring mat. In the few cases where we were unable to measure newborn body length at the hospital ( n = 35), we used body length information obtained from hospital crib cards ( n = 2) or medical records ( n = 33), if available. In statistical analyses, using the mean of all triplicate AGD and body length measures, we created an anogenital index (AGI, mm/m) by dividing each AGD value (in mm) by the body length (in m) at the corresponding timepoint. We normalized AGD to body length in place of birthweight because prior studies suggest that body length is better correlated with AGD across infancy, and adjusting for body length may make the index less dependent on changes in body size and age ( 47 , 65 ). Finally, we calculated % growth in each AGI Short and AGI Long measure from birth to 5 months of age: 100 × [(5-month AGI − newborn AGI)/newborn AGI)]. We identified and excluded one female newborn's AGD Long value as it was unrealistically short (average of triplicate measures = 13.95 mm).
Of the 563 mother-infant dyads with available maternal nonpersistent EDC data, 505 infants had at least one AGD measure. A total of 492 mother-infant dyads additionally had data on all covariates of interest, with the final sample size for each analysis being dependent on the exposure/outcome relationship, as outlined in Fig. S1 ( 63 ).
We considered many covariates using a directed acyclic graph (DAG) ( 66 , 67 ) and the prior literature (including our own) ( 18 , 38 , 39 , 41 , 45 , 51 , 61 , 62 , 68-80 ). We used the DAG to identify the minimum sufficient adjustment set of covariates ( 67 , 81 ), including those representing sociodemographic, lifestyle, and health factors, as well as several precision variables. We tested for multicollinearity using correlation coefficients, but the selected covariates were weakly correlated with each other (r < 0.35; data not shown ). The final covariate-adjusted statistical models included race/ethnicity, maternal age, parity, conception season, pre-pregnancy BMI, mean pregnancy AHEI-2010, smoking in the first trimester, and early pregnancy perceived stress (operationalized as in Table 1 ). In analyses evaluating 5-month AGD and % growth in AGD from birth to 5 months, we included childcare arrangements categorized as “in-home care” and “out of home care” (which also potentially account for mode of feeding from birth to 5 months, data not shown ), as a precision variable, as it was associated with some 5-month AGD measures. In linear mixed models, we operationalized childcare arrangement as a three-category variable, including the category “did not participate in 5-month visit” to allow inclusion of all infants who completed their visit at birth but did not participate in the 5-month visit, and to improve model fit.
Baseline demographic and lifestyle characteristics of pregnant women in I-KIDS
a
Includes nonHispanic Black, Asian, Native Hawaiian or other Pacific Islander, American Indian or Alaska Native, Multiracial, and Others. b Covariates included in the model. c Missing income and alcohol consumption n = 4, n = 2 female and n = 2 male. d n = 402 women ( n = 208 carrying females, n = 194 carrying males) had medical record data. e n = 10 on progesterone, n = 19 on levothyroxine, n = 1 on both progesterone and levothyroxine. f N = 400 women ( n = 198 females, n = 202 males) did not participate in 5-month visit. g Excluding alcohol.
Abbreviations: AHEI, alternative healthy eating index; BMI, body mass index; AGD, anogenital distance.
We approximated exposure to phthalate/replacement parent compounds that are metabolized and excreted as multiple urinary metabolites using the following molar-sum (in nmol/mL) equations: sum of di(2-ethylhexyl) phthalate metabolites (ΣDEHP) = (MEHP/278) + (MEHHP/294) + (MEOHP/292) + (MECPP/308), sum of di-n-butyl phthalate metabolites (ΣDBP) = (MBP/222) + (MHBP/238), sum of di-iso-butyl phthalate metabolites (ΣDiBP) = (MiBP/222) + (MHiBP/238), sum of di(isononyl) cyclohexane-1,2-dicarboxylate metabolites (ΣDiNCH) = (MHiNCH/314) + (MCOCH/328), and sum of di(2-ethylhexyl) terephthalate metabolites (ΣDEHTP) = (MEHHTP/294) + (MECPTP/308) (Table S1) ( 63 ). In all single chemical analyses, we approximated exposure to di-isononyl phthalate (ΣDiNP) as ΣDiNP = (MCOP/322) + (MONP/306) for n = 327 participants who had a reported value for MONP and ΣDiNP = (MCOP/322) + (MiNP/292) for n = 165 participants for whom MONP was not quantified. In mixtures analyses which do not allow for missing data, we chose to only include MCOP as all participants had these data and as our single pollutant analyses demonstrated that associations of ΣDiNP with our outcomes were consistent with those evaluating MCOP alone. For statistical analyses evaluating parent compounds, molar concentrations of parent compounds were converted back to ng/mL by multiplying ∑DEHP, ∑DiNP, ∑DBP, ∑DiBP, ∑DiNCH, and ∑DEHTP by molecular weights of MECPP, MCOP, MBP, MiBP, MHiNCH, and MECPTP, respectively ( 18 , 61 , 82 , 83 ).
We used all instrument-read values to avoid bias from imputing values below the limit of detection (LOD) ( 84 ). In our statistical analyses, we only included chemical biomarkers with concentrations greater than zero in at least 90% of women (Table S1) ( 63 ). This resulted in butylparaben, BPF, and triclocarban being excluded from further analyses. To prevent undefined estimates for ln-transformed zero concentrations (MEHP = 18; MiNP = 4; MHBP = 6; MHiNCH = 3; MCHOCH = 46; ethylparaben = 6; BPA and TCS = 2; BPS and BP-3 = 1), we used the formula: [ln(chemical biomarker concentration + 0.0001)] ( 85 ). We accounted for urine dilution using specific gravity, which we calculated using the following formula: P c = P [( SG − 1)/( SG i − 1)], where P c is the specific gravity-adjusted chemical metabolite concentration, P is the measured chemical metabolite concentration (ng/mL), SG is the median specific gravity of the pooled urine samples (1.016) and SG i is the specific gravity of each individual urine sample ( 86 ). In all statistical analyses, we ln-transformed nonpersistent EDC biomarkers to improve model fit and interpretation.
To evaluate associations of individual nonpersistent EDC biomarkers with AGD Short and AGD Long from birth to 5-months and at each timepoint, we used covariate-adjusted linear mixed models, including all covariates as fixed effect variables. We specified an unstructured 2 × 2 covariance matrix for each model's residuals. To identify whether associations between individual nonpersistent EDC biomarkers with AGD differed depending on timepoint of AGD measurement, we included a time × EDC biomarker interaction term and extracted timepoint-specific effect estimates and 95% confidence intervals (CI).
At birth and 5 months, we used covariate-adjusted quantile g-computation (QGComp) to model the joint effect of each EDC mixture (in quartiles) and identify the most important contributors within each mixture ( 87 , 88 ). Because MEHHTP and MECPTP (DEHTP metabolites) were not measured at the CDC until our study's second urine shipment, and based on findings from our single-chemical analyses highlighting the importance of DEHTP metabolites, our mixture analyses included women with all chemical data. Based on our prior findings showing little correlation between phthalate and phenol biomarkers in this cohort ( 61 , 62 ), we modeled two separate EDC mixtures in QGComp: the phthalate/replacement mixture composed of 17 individual phthalate metabolites (MEHP, MEHHP, MEOHP, MECPP, MCOP, MCNP, MCPP, MBzP, MEP, MBP, MHBP, MiBP, MHiBP, MHiNCH, MCOCH, MEHHTP, MECPTP) and the phenol mixture composed of nine compounds (methylparaben, ethylparaben, propylparaben, BPA, BPS, BP-3, Triclosan, 2,4-DCP, and 2,5-DCP). We excluded smoking in the first trimester from mixtures models due to poor distribution and model fit. We fit models with 500 bootstraps to obtain effect estimates with more precise confidence intervals. Then, we generated results without bootstrapping to determine the partial negative and positive associations and weights, which denote general importance and direction of each nonpersistent EDC biomarker contributing to the mixture effect. Although QGComp offers a novel extension to accommodate longitudinal data ( 89 ), this extension only reports results from the exposure-outcome association across all timepoints, whereas our single-chemical analyses pointed to distinct timepoint-specific findings. In figures, we reported the most prominent contributors to the observed associations using weight thresholds of 1/n, where n represents the number of EDC biomarkers contributing in the direction of the observed effect estimate ( 90 ). In our summary of results below, for brevity, we reported the top three contributors.
To identify potential nonlinear relationships of joint EDC biomarkers with AGD at birth and 5 months, we conducted hierarchical Bayesian Kernel Machine Regression (BKMR). Briefly, BKMR can flexibly estimate the multivariable exposure-response function, allows for a hierarchical variable selection approach to identify major mixture components, and accommodates highly correlated exposures in the mixture ( 91 , 92 ). Importantly, hierarchical BKMR allows us to identify the relative importance of individual metabolites and chemical classes by treating phthalate/replacements and phenols as two separate groups within one EDC mixture ( 62 ). The BKMR mixture included 26 nonpersistent EDC biomarkers: 17 phthalate/replacement metabolites (MEHP, MEHHP, MEOHP, MECPP, MCOP, MCNP, MCPP, MBzP, MEP, MBP, MHBP, MiBP, MHiBP, MHiNCH, MCOCH, MEHHTP, MECPTP), and nine phenols (methylparaben, ethylparaben, propylparaben, BPA, BPS, BP-3, triclosan, 2,4-DCP, and 2,5-DCP). We ln-transformed, centered, and scaled all EDC biomarkers and continuous covariates; the outcomes were only centered and scaled ( 62 , 91 , 92 ). Categorical covariates were included and not modified; however, smoking was excluded due to poor model fit. We fit hierarchical BKMR models using the Gaussian family and 100 000 iterations. We assessed cumulative nonlinear mixture associations using dose-response curves across quantiles of the full mixture. To identify the important chemicals driving the observed associations of the EDC mixture with AGD measures, we used posterior inclusion probabilities (PIPs).
We used multivariable linear regression models to evaluate associations of nonpersistent EDC biomarkers with % growth in AGD from birth to 5 months. QGComp and BKMR were conducted as described above.
We conducted several sensitivity analyses to evaluate the robustness of our findings, focusing on our QGComp mixtures findings, as this allowed us to understand general trends across our various exposures and outcomes. Maternal and pregnancy health conditions (intrauterine growth restriction (IUGR), progesterone or levothyroxine use, and diagnoses of hypertensive disorders of pregnancy (HDP) and gestational diabetes mellitus (GDM)) for sensitivity analyses were abstracted from electronic medical records, which were available for 82% of the analytic sample. In our first sensitivity analysis, we removed pregnancies diagnosed with IUGR ( n = 8 females). Second, we excluded infants born preterm ( n = 12 females, n = 7 males). Third, we removed infants whose mothers took levothyroxine or progesterone during pregnancy ( n = 18 females, n = 12 males). Fourth, we removed infants whose mothers were diagnosed with any HDP ( n = 16 females, n = 17 males). Fifth, we removed infants whose mothers were diagnosed with GDM ( n = 8 females, n = 8 males).
In linear mixed models and linear regression models, all β-estimates and 95% CIs were estimated as [β × ln(2)] to represent a mm/m change in AGI or % growth in AGI from birth to 5 months for each doubling in individual nonpersistent EDC biomarker. Results from QGComp models are interpreted as the mm/m change in AGI or % growth in AGD from birth to 5 months if all EDC biomarkers in the mixture increase by a quartile (25%).
All linear mixed models and linear regression analyses were conducted in SAS Software, version 9.4 (SAS Institute Inc, Cary, NC) using PROC MIXED and PROC GLM, respectively. We performed QGComp and BKMR analyses in R Software using the R packages “qgcomp: Quantile G-Computation” ( 87 ) and “bkmr: Bayesian Kernel Machine Regression” ( 92 , 93 ), respectively. Relationships were considered statistically significant at P < .05. Given the exploratory nature of our study, and following prior recommendations, we did not adjust for multiple comparisons ( 94-96 ).
Text is read by the "Ask this paper" AI Q&A widget below.
Extraction quality varies by source — PMC NXML preserves structure
cleanly, OA-HTML may include some navigation residue, and OA-PDF can
have broken hyphenation. The publisher copy
(via DOI)
is the canonical version.