The Joint Effect of Prenatal Exposure to Metal Mixtures on Neurodevelopmental Outcomes at 20–40 Months of Age: Evidence from Rural Bangladesh

Background: Exposure to chemical mixtures is recognized as the real-life scenario in all populations, needing new statistical methods that can assess their complex effects. Objectives: We aimed to assess the joint effect of in utero exposure to arsenic, manganese, and lead on children’s neurodevelopment. Methods: We employed a novel statistical approach, Bayesian kernel machine regression (BKMR), to study the joint effect of coexposure to arsenic, manganese, and lead on neurodevelopment using an adapted Bayley Scale of Infant and Toddler Development™. Third Edition, in 825 mother–child pairs recruited into a prospective birth cohort from two clinics in the Pabna and Sirajdikhan districts of Bangladesh. Metals were measured in cord blood using inductively coupled plasma-mass spectrometry. Results: Analyses were stratified by clinic due to differences in exposure profiles. In the Pabna district, which displayed high manganese levels [interquartile range (IQR): 4.8, 18μg/dl], we found a statistically significant negative effect of the mixture of arsenic, lead, and manganese on cognitive score when cord blood metals concentrations were all above the 60th percentile (As≥0.7μg/dl, Mn≥6.6μg/dl, Pb≥4.2μg/dl) compared to the median (As=0.5μg/dl, Mn=5.8μg/dl, Pb=3.1μg/dl). Evidence of a nonlinear effect of manganese was found. A change in log manganese from the 25th to the 75th percentile when arsenic and manganese were at the median was associated with a decrease in cognitive score of −0.3 (−0.5, −0.1) standard deviations. Our study suggests that arsenic might be a potentiator of manganese toxicity. Conclusions: Employing a novel statistical method for the study of the health effects of chemical mixtures, we found evidence of neurotoxicity of the mixture, as well as potential synergism between arsenic and manganese. https://doi.org/10.1289/EHP614


Introduction
Childhood exposure to neurotoxicants is a potential impediment to economic development, as it is most prevalent in developing countries, making this issue particularly poignant in countries such as Bangladesh (Suk et al. 2003;Grandjean et al. 2015). Growing evidence from animal research indicates that the central nervous system is the most vulnerable of all body systems to chemical injury during development (Faustman et al. 2000;Rodier 2004). One of the most widely studied categories of neurotoxicants is metals. Among metals, arsenic, lead, and manganese are prevalent in the environment and have evidence of neurotoxicity. These three metals are thus ideal candidates on which to test new statistical methodologies for mixtures. Arsenic, lead, and manganese exposure is widely prevalent in Bangladesh (Kile et al. 2009), and share the central nervous system as the primary toxicity target in children (Bressler et al. 1999;Clarkson 1987;Pola nska et al 2013;Vahter 2008;Zoni and Lucchini 2013). Exposure to lead even at low levels is commonly accepted as neurotoxic. The neurotoxic effects of arsenic and manganese at levels commonly found in the environment are less well understood, but emerging evidence suggests they too are a concern. Prior to our work, manganese-arsenic interaction studies in the Bangladeshi population were cross-sectional and lacked adequate power to assess interactions among mixture components (Wasserman et al. 2006). An epidemiologic investigation in Mexico found evidence of both an inverted "U" relationship between blood Mn and infant development (i.e., both low and high blood Mn levels were associated with poorer performance) and of a lead-manganese interaction being synergistically more toxic (Claus Henn et al. 2010, 2012. Previous studies on the Bangladeshi population (Wasserman et al. 2004(Wasserman et al. , 2006(Wasserman et al. , 2007(Wasserman et al. , 2008Hamadani et al. 2011) have shown that arsenic exposure during childhood through drinking water is negatively associated with cognition of school-age children. How-ever, this exposure has not been found to be associated with cognitive development at earlier stages in life (Tofail et al. 2009;Hamadani et al. 2010). A recent study of the independent effect of water manganese exposure among school-age Bangladeshi children exposed to low-level arsenic found evidence of manganese neurotoxicity but no evidence of arsenic effects on neurodevelopment (Wasserman et al. 2006). Our group has recently evaluated, using traditional linear regression approaches, the association between postnatal exposure to heavy metals and Bayley neurodevelopment scores measured at 20-40 mo (Bayley 1993). The analyses were conducted in the same population considered in the present study (Rodrigues et al. 2016). The study reported neurotoxic effect of 24-mo exposure to blood lead and water arsenic, as well as an inverted-U dose-response relationship between water manganese and cognitive development. Recently, more attention has been directed towards studying the joint effects of environmental metal mixtures, that is, investigating interactions that may characterize the joint effect of mixtures (Wright et al. 2006;Claus Henn et al. 2012. Traditionally, mixtures have been studied via multivariable parametric regression approaches that concomitantly adjust for the confounding effects of mixture components and estimate the independent effect of each component, adjusting for the others. If multiple metals do act as a mixture, this approach would be limited by both multicollinearity and model misspecification. Moreover, it is challenging to specify a correct parametric model that incorporates the possibility of any type of interaction and nonlinear effects among multiple concurrent exposures; the likelihood that all components of a mixture will always have linear effects seems remote. Statistical models designed to address mixtures are relatively new, and several approaches are now available (Bobb et al. 2015;Feder et al. 2009;Gennings et al. 2013;Zanobetti et al. 2014) that address the extremely complex questions that underlie the relationships between environmental mixed exposures and their health effects.
Our study contributes to the scientific literature by providing new evidence on the effects of metal mixtures on neurodevelopment in Bangladesh, and by employing a novel statistical approach, Bayesian kernel machine regression (BKMR) (Bobb et al. 2015), which flexibly models the joint effect of the mixture components, allowing for potential interactions and nonlinear effects. Specifically, our study provides a prospective analysis (prior studies were cross-sectional) and focuses on child neurodevelopment at 20-40 mo, evaluating its association with prenatal coexposures to lead, arsenic, and manganese. We used umbilical cord blood arsenic, manganese, and lead concentrations as biomarkers of late pregnancy exposure. We focused on prenatal exposure, as there is evidence that this time period is one of heightened susceptibility (Mazumdar et al. 2011). The BKMR approach allows us to study the joint effect of the components of the mixture. Further, we are able to disentangle how the joint effect comes about, allowing for both interactions and nonlinear effects. In particular, we examined a) whether exposure to this mixture jointly is associated with adverse neurodevelopmental effects; b) the dose-response relationships between combinations of metal exposures and cognition; and c) whether the impact of a metal is more pronounced when it occurs as part of a mixture (i.e., whether the components of the mixture interact).

Study Population
The study population has been described previously (Gleason et al. 2014;Kile et al. 2014). Briefly, between 2008 and 2011, we recruited pregnant women to participate in a prospective study to examine the effects of chronic low-level arsenic exposure on reproductive health outcomes. Participants were recruited from two rural health clinics operated by the Dhaka Community Hospital Trust (DCH) in the Sirajdikhan and Pabna Sadar upazilas of Bangladesh. Between 2010 and 2013, healthcare workers invited families to participate in a follow-up study to examine children's neurodevelopment. The study base for this analysis consisted of all children born during the reproductive cohort study who were 20 to 40 mo of age. This study was approved by the Human Research Committees at the Harvard T.H. Chan School of Public Health (Harvard Chan School) and Dhaka Community Hospital (DCH). Boston Children's Hospital formally ceded review of the follow-up study to the Harvard Chan School. Informed consent was obtained from all participants or their parents before enrollment and prior to engaging in any study activities.
A total of 1,613 women were enrolled in the reproductive health study during early pregnancy, of whom 1,458 women had a confirmed singleton pregnancy. Among those with a live birth, a total of 964 children participated in follow-up activities, including 827 who underwent neurodevelopmental assessments. The sample for which both cord blood and neurodevelopmental scores were available was 825, and served as the sample for this analysis (see Table S1 for comparison of baseline maternal and child characteristics in the reproductive health and neurodevelopment studies).

Arsenic, Manganese, and Lead Exposure
We collected umbilical cord venous blood in trace element-free tubes at time of delivery from participants. Samples were kept at 4 C and shipped to the Trace Metals Laboratory at the Harvard Chan School. All samples were processed in a dedicated trace metal clean room outfitted with a Class 100 clean hood and using glassware that was cleaned by soaking in 10% HNO 3 for 24 h and rinsed several times with 18X deionized water. Blood samples were prepared for measurement of arsenic, manganese, and lead concentrations by first weighing ( ∼ 1 g) and then digesting samples in 2 mL concentrated nitric acid for 24 h at room temperature. These samples were then treated overnight with 30% hydrogen peroxide (1 mL per 1 g of blood) and then diluted to 10 mL with deionized water. Acid-digested samples were then analyzed for metal concentrations with a dynamic reaction cellinductively coupled plasma mass spectrometer (Perkin Elmer). The average of five replicate measurements for each individual sample was reported as the final value.

Neurodevelopmental Outcomes
We translated the Bayley Scales of Infant and Toddler Development™, Third Edition (BSID-III™) (Bayley 2006) into Bengali and adapted it for use in rural Bangladesh (M.M.). Two primary outcomes were derived by summing across raw scores of cognitive and language development for each participant: a) raw cognitive development score (CS); and b) raw language development composite score (LCS; sum of the raw scores on the expressive and receptive scales). Trained study personnel who were unaware of participants' umbilical cord blood metal levels administered the tests using standard protocols. An expert child neurologist (M.M.) and neuropsychologist (D.C.B.) oversaw administration of the BSID-III™ and quality control, which included frequent site visits and review of videotaped administration of neurodevelopmental assessments.

Covariates
We collected demographic information using structured questionnaires at three scheduled clinic visits during pregnancy to obtain data on age, education, smoking history, and socioeconomic status. Gestational age at birth was determined by ultrasound measurement taken at time of enrollment (<16 wk gestational age). Infant sex, birth weight, length, and head circumference were recorded at birth and abstracted from medical records. Trained study staff administered questionnaires to the primary caregiver when the child was 12 mo and at the time of the Bayley Scale administration (20-40 mo) to collect medical histories, demographic information, and dietary information. A measure of maternal protein intake during pregnancy was derived from frequency (amount per week) of fish, meat, and egg consumption. The variable used in our analysis was constructed by categorizing the summary measure for nutritional status with the 25th and 75th percentiles as cutoff points. Maternal IQ was assessed using the Raven's Progressive Indices (Raven 1981). Interviewers also administered the HOME instrument that had been previously translated and adapted for use in Bangladesh (Black et al. 2004).

Statistical Analysis
Distributional plots and descriptive statistics were examined for all variables by clinic (Sirajdikhan and Pabna). We modeled manganese, lead, and arsenic as natural log-transformed and centered continuous variables. This transformation achieves a common scale, and accounts for the severe right-skewedness of the metals' concentrations. Neurodevelopment scores were approximately normally distributed and were modeled as z-scored continuous outcomes.
All statistical procedures were performed using R (version 3.1.3; R Foundation for Statistical Computing).
Covariates. The following covariates were included in all models to adjust for confounding: child's sex, child's age at the time of neurodevelopmental testing, mother's age at the time of the child's birth, maternal education (less than high school vs. at least high school), maternal IQ, HOME score, secondhand smoke exposure at baseline (smoking environment vs. nonsmoking environment), and protein intake (low vs. medium vs. high protein intake). Protein effect was found to be linear, and was therefore modeled as ordinal in the analyses. Continuous covariates (age at testing, mother's IQ, mother's age at birth, HOME score) were modeled allowing for quadratic effects. Four percent of eligible children were excluded (n = 33) because of missing covariate data; we assumed these data were missing at random.
Linear regression. We began by conducting multivariable linear regression analysis with models of the form: where Y is the continuous neurodevelopment phenotype (CS, LCS); Mn, As, and Pb are the centered log concentrations of manganese, arsenic, and lead, respectively; and Z = Z 1 , . . . ,Z p are additional p potential confounders. Bayesian kernel machine regression. To allow for potential synergistic and nonlinear effects among mixture components, we implemented BKMR, a new statistical approach for multipollutant mixtures that flexibly models the joint effect of the mixture using a kernel function (Bobb et al. 2015). (Details on the approach are described in Supplemental Material, "Bayesian Kernel Machine Regression" section.). We also provide code for the analyses conducted (Supplemental Material, R Code). This novel method allows estimation of nonlinear and nonadditive dose-response functions for a (potentially high-dimensional) set of correlated exposures accounting for uncertainty. A key feature of BKMR is that the form of the exposure-response function is flexibly modeled, and does not need to be specified a priori (e.g., using quadratic terms). The BKMR model is given below.
The function hðÞ is an exposure-response function that accommodates nonlinearity and/or interaction among the mixture components. In this setting, it can be challenging to specify explicitly a set of basis functions (e.g., spline or polynomial terms) to represent hðÞ. For example, to allow for nonlinearity and interaction by using a spline basis with three degrees of freedom (DF) and including all the interaction terms would require estimating 63 parameters in the case of 3 metals, 255 parameters in the case of 4 metals, and 1,023 parameters in the case of 5 metals [more generally, ð1 + DFÞ M − 1 parameters in the case of M metals]. BKMR handles the potentially high-dimensional parameter space of the exposure-response function by using a kernel exposureresponse machine representation for hðÞ. The kernel machine representation enables formulating model (2) as a linear mixed model , which achieves regularization of the estimated exposure-response function and provides several computational advantages (details in Bobb et al. 2015). Intuitively, the resulting linear mixed model assumes that two individuals who have similar exposure profiles (e.g., similar values of z i = (As i ,Mn i ,Pb i )) will have similar estimated health effects (e. g., h i = hðz i Þ and will be close to h j = hðz j Þ for individuals i and j), where the similarity (or distance) between individuals' exposure profiles is measured using the kernel function.
We next provide details on BKMR specification and model fitting. Note that fitting BKMR depends on the choice of kernel function. We used the Gaussian kernel, which flexibly captures a wide range of underlying functional forms for hðÞ and has been shown to work well in simulation studies based on realistic exposure-response scenarios (Bobb et al. 2015). Fitting the BKMR model yields an estimate of the exposure-response function hðÞ and its uncertainty. From the model fit, we also derived summary statistics that quantify scientifically relevant features of the exposure-response function in order to gain insight on the joint effect of the mixture. Credible intervals obtained from BKMR fit incorporate the additional uncertainty due to estimation of a highdimensional set of exposures, which accounts for multiple-testing penalty (Scott and Berger, 2010).
Generalized additive models. In order to assess the reproducibility of the BKMR results and to test for potential interactions among mixture components, we subsequently fitted multivariable generalized additive models (GAMs). GAMs are an alternative approach to flexibly model dose-response curves, and can be useful to assess the validity of prior assumptions of BKMR in lowerdimensional settings. To test for interaction, we compared two specifications of GAM employing the tensor product smoother, tiðÞ, one that includes the smoothed terms for the metals additively [i.e., tiðAs i ,Mn i ,Pb i Þ = tiðAs i Þ + tiðMn i Þ + tiðPb i Þ] and one that allows for a joint smoothed effect of arsenic and manganese [i.e., tiðAs i ,Mn i ,Pb i Þ].
Sensitivity analyses. Finally, we performed three sensitivity analyses: a) we assessed the influence of extreme observations by fitting models with and without outliers. Outliers were identified using the extreme Studentized deviates statistics with the procedure proposed by Rosner which specifies a cut-off for outlier selection that optimizes the number of points to exclude (Rosner, 1983); b) we also evaluated the influence of prior specification for BKMR by fitting the approach with two alternative prior assumptions allowing for different degrees of smoothness (see Supplemental Material, Sensitivity Analyses for BKMR); and c) finally, we assessed the robustness of our results by further adjusting for child exposure to heavy metals at 20-40 mo, considered in Rodrigues et al. (2016). Arsenic and manganese concentrations in drinking water and blood lead level measures were only available for 524 children. As precision is reduced due to the smaller numbers of observations, we fitted BKMR using grouped variable selection (Bobb et al. 2015):

Study Population Characteristics
Demographic characteristics of participating mother-infant pairs from Pabna and Sirajdikhan clinics in Bangladesh are presented in Table 1. Important differences between these two study sites were observed. Compared to Sirajdikhan, mothers in Pabna were less likely to complete secondary education, more likely to live in a smoking environment, and less likely to report a low-protein intake diet. Compared to Sirajdikhan, children in Pabna had lower birthweight and lower CS. Cord blood manganese and arsenic concentrations were significantly higher in Pabna: manganese [geometric mean ± standard deviation ðSDÞ (range)]: Pabna 9:7 ± 2:6 ð1:7 − 303:1Þ lg=dL, Sirajdikhan 5:4 ± 1:7 ð1:2 − 88:6Þ lg=dL; arsenic: Pabna 1:0 ± 2:2 ð0:1 − 27:7Þ lg=dL, Sirajdikhan 0:4 ± 1:9 ð0:1 − 7:4Þ lg=dL. Cord blood lead concentrations were significantly higher in Sirajdikhan: lead [mean ± SD (range)]: Pabna 1:8 ± 1:9 ð0:3-79:1Þ lg=dL, Sirajdikhan 6:0 ± 1:9ð1:0 − 36:0Þ lg=dL. Metals concentrations were significantly correlated only among participants from the Pabna clinic (Figure 1), indicating that sources of metals might differ across sites. Given the important differences in exposure profiles and baseline characteristics observed across the two study sites, regression models were stratified by clinic in primary analyses and adjusted for clinic in secondary analyses. Table 2 displays associations between metals and neurodevelopment outcomes from multivariable linear models stratified by clinic (see Tables S2-S3 for complete regression output). Among the potential confounders, child's age at testing, clinic, and maternal education were most strongly associated with neurodevelopment outcomes. Protein intake was significantly associated with CS and LCS. HOME inventory scale was significantly associated with CS (Tables S2-S3).

Multivariable Regression Analyses
Linear regression analyses of data from the Pabna clinic revealed a significant negative association of manganese with CS adjusting for arsenic and lead [CS: − 0:2 ð − 0:4, − 0:02Þ]. Given these estimates, we expect CS scores to decline by 0.2 SDs per interquartile range (IQR) change in log manganese concentration (log ðMn lg=dlÞ IQR = 1:3).

Bayesian Kernel Machine Regression Analyses
In order to relax the assumptions of linearity and additive effects imposed in the aforementioned regression models, we implemented BKMR, which yields an estimate of the joint exposureresponse function of the three metals.
Inferences on mixture effects were obtained by computing posterior mean estimates (and 95% posterior credible intervals) of the change in CS associated with changes in the level of each of the mixture components (see Figures 2-3 and Figures S1-S4; Table 3 summarizes these findings). In particular, we estimated a) the cumulative effect of the mixture by estimating the expected change in CS associated with concurrent changes in all of the components of the mixture from their median level; b) the effect of an IQR change of each metal on CS and potential interactions among the metals by estimating the change in the CS for a change in the component of interest from its 25th to 75th percentile, while setting the other metals at the median, 25th, or 75th percentile levels; and c) the dose-response relationship of each mixture component and potential interactions among the metals by estimating the predicted CS for each level of the component of interest, setting the other metals at the median, 25th, or 75th percentile levels. We first display numerical summaries of the overall effect of the mixture, defined as the change in the (z-scored) neurodevelopment scores associated with a simultaneous change in each of the three metals from a particular threshold (25th percentile to 75th percentile) as compared to when those metals are each at their median value (50th percentile) (Figure 2A). We found a joint toxic effect of the mixture in the Pabna clinic. In particular, the joint effect was statistically significant when all metals were at or above their 60th percentile, as compared to when all metals were at their median values, and the toxicity increased at higher levels of the three joint exposures. We then sought to disentangle which of the three metals might be driving this overall association ( Figure 2B) by estimating univariate summaries of the change in the neurodevelopment scores associated with a change in a single pollutant from its 25th percentile to 75th percentile, where all of the other pollutants are fixed at a particular threshold (25th, 50th, or 75th percentile). We found that manganese is the only metal displaying a significant effect in this clinic. Its negative association with CS appears stronger at higher percentiles of other metals. A change in log manganese concentration from the 25th to the 75th percentile is associated with a significant decrease in CS of − 0:2 ð− 0:5, 0:0Þ, − 0:3 ð− 0:5, − 0:1Þ and − 0:4 ð − 0:6, − 0:2Þ SDs when lead and arsenic are set at the 25th, 50th, and 75th percentiles, respectively. To investigate potential nonlinearity of the exposure response function, we then estimated the univariate exposure-response functions ( Figure  2C). The plot shows a suggestion of nonlinear effects of manganese and arsenic. We note that although we observed a U-shaped dose-response curve for manganese, there is considerable variability about the curve for high manganese concentrations. The single-pollutant estimates from Figure 2B suggested possible interaction of the mixture. In particular, we saw that the estimated association between manganese and neurodevelopment scores increased as lead and arsenic both increased from their 25th to their 75th percentiles. To investigate this further, we then plotted bivariate cross-sections of the exposure-response function ( Figure 2B). The figure suggests a potential interaction between arsenic and manganese in their association with CS, whereby the significant negative slope of manganese becomes steeper at higher levels of arsenic concentration. Conversely, the nonsignificant positive slope of arsenic becomes flat at higher levels of manganese concentration. However, the confidence intervals in Figure 2B are highly overlapping, indicating lack of statistical significance. The plots in Figure 3 give summaries of the fit of the model in the Sirajdikhan clinic. Here, we find weaker evidence of a joint toxic effect ( Figure 3A). Although not statistically significant, lead appears to be the most neurotoxic metal in the Sirajdikhan clinic, where lead levels are higher than those found in the Pabna clinic ( Figure 3B) with point estimates close to what was estimated by multivariable linear regression ½ − 0:1 ð − 0:2, 0:0Þ. In contrast with the Pabna population, the Note: Models adjusted for metals listed above, child gender, mother IQ, maternal education, maternal protein intake, smoking environment, age at testing, and maternal age at delivery; models allow for quadratic effect of age at testing, mother IQ and age at birth, HOME score. Environmental Health Perspectives 067015-5 effects appear linear ( Figure 3C) and noninteractive (e.g., slopes for each metal are similar at varying levels of the other metals; Figure 3D). To formally test for the presence of interaction between arsenic and manganese as suggested by the BKMR analyses, we then fitted the multivariable generalized additive regression model allowing for arsenic and manganese interaction, adjusting for the smoothed term for lead and covariates in the Pabna clinic. To test for arsenic and manganese interaction, we compared two specifications of GAM, employing the tensor product smoother, tiðÞ: one that includes the smoothed terms for the metals additively [i.e., tiðAs i ,Mn i ,Pb i Þ = tiðAs i Þ + tiðMn i Þ + tiðPb i Þ] and the other that allows for a joint smoothed effect of arsenic and manganese [i.e., tiðAs i ,Mn i ,Pb i Þ = tiðAs i Þ + tiðMn i Þ + tiðAs i ,Mn i Þ + tiðPb i Þ]. The analyses yielded suggestive evidence of potential interaction (p = 0:1). Note that our stratified analyses might be underpowered to detect interactions. In the next section, we evaluate arsenic-manganese interaction in the whole sample adjusting for clinic. GAM analyses support the results obtained from fitting the BKMR, indicating manganese as the metal most predictive of CS in Pabna. The approach found evidence of a nonlinear U-shaped effect of manganese (results not shown). The same model was fitted for the sample of mother-infant pairs in Sirajdikhan and, again, no evidence of interaction (p = 0:9) was found, and linearity was confirmed.
In models of language composite scores, manganese was also negatively associated among Pabna participants, though associations were weaker than in models of CS ( Figure S3). In Sirajdikhan, arsenic was associated with lower LCS ( Figure S4). Model adjusted for child gender, maternal IQ, maternal education, maternal protein intake, smoking environment, age at testing, and maternal age. (A) Overall effect of the mixture (estimates and 95% credible intervals). This plot compares the neurodevelopment score when all exposures are at a particular quantile to when all are at the 50th percentile. (B) Single pollutant association (estimates and 95% credible intervals, gray dashed line at the null). This plot compares the neurodevelopment score when a single pollutant is at the 75th vs. 25th percentile, when all the other exposures are fixed at either the 25th, 50th, or 75th percentile. (C) Univariate exposure-response functions and 95% confidence bands for each metal with the other pollutants fixed at the median (D) Bivariate exposure-response functions for: arsenic when manganese is fixed at either the 25th, 50th, or 75th percentile and lead is fixed at the median (top left panel); arsenic when lead is fixed at either the 25th, 50th, or 75th percentile and manganese is fixed at the median (top right panel); manganese when arsenic is fixed at either 25th, 50th, or 75th percentile and lead is fixed at the median (bottom left panel); manganese when lead is fixed at either the 25th, 50th, or 75th percentile and arsenic is fixed at the median (bottom right panel). [percentiles: (0.25,0.5,0.75); As = ð0:26,0:40,0:60Þ, Mn = ð3:93,5:14,6:66Þ, Pb = ð3:87,6:12,9:67Þ].

Regression Analyses Adjusted by Clinic
In multivariable linear regression models adjusting for clinic (Tables S2-S3), an IQR change in manganese (log ðMn lg=dlÞ IQR = 0:8) was significantly associated with lower CS [ − 0:1 ð − 0:2, − 0:02Þ]. An IQR change in arsenic [log ðAs lg=dlÞ IQR = 1:1] was marginally associated with higher scores [CS: 0:08 ð − 0:01,0:2Þ]. These results are similar to those observed in the Pabna clinic, suggesting that this clinic may be driving the observed associations. Note that the magnitude of manganese effect appears weaker in the whole sample because the IQR is smaller in the Sirajdikhan clinic than in the Pabna clinic.
BKMR analyses estimated significant toxic effects of manganese and lead. The analyses yielded similar results for manganese as in models with only Pabna participants, and slightly stronger evidence for lead toxicity relative to models with only Sirajdikhan participants ( Figure S1). An IQR change in the log-transformed lead [log ðPb lg=dlÞ IQR = 1:4] was associated with reduced CS by − 0:1 ð − 0:2, − 0:03Þ SDs when arsenic and manganese were at the median level. Compared to lead and arsenic, manganese appears to be the most neurotoxic component of the mixture. Finally, fitting the GAM model for all study participants adjusting for clinic yielded similar results to the GAM model for Pabna participants. Allowing for an interaction between arsenic and manganese significantly improved the fit of the model [p-value = 0:04].

Sensitivity Analyses
When we excluded three outliers of cord blood levels (three highest manganese levels and one highest arsenic level in Pabna clinic), results from BKMR analyses were similar. Generalized additive model analyses yielded weaker evidence of arsenicmanganese interaction in the Pabna clinic (p = 0:3). However, when fitting the GAM model for all study participants adjusting for clinic, arsenic-manganese interaction was still statistically significant (p = 0:01).
BKMR was fitted considering two alternative prior specifications on the smoothness of the dose-response function (higher and lower with respect to the main analyses). Results are robust to changes in prior specification. Sensitivity analyses using BKMR with grouped variable selection that include water arsenic, water manganese, and blood lead at 24 mo, in the smaller sample for which such information was available (n = 524), reveal that cord blood metal exposures are the most important drivers of neurodevelopement at 20-40 mo. BKMR estimates a 99% posterior inclusion probability for the prenatal exposures vs.
28% posterior inclusion probability for the postnatal exposures in the model for CS and 85% posterior inclusion probability for the prenatal exposures vs. 23% posterior inclusion probability for the postnatal exposures in the model for LCS. In stratified analyses by clinic as well as analyses adjusted by clinic, the dose-response relationships for the 24-mo metal exposures are linear, additive, and nonstatistically significant. Even in this smaller sample for which postnatal metal concentrations were available, results for cord blood concentrations of arsenic, manganese, and lead are robust to the adjustment of postnatal exposures.

Discussion
Using BKMR, a new, nonparametric statistical method to capture the complexity of mixed exposures, we evaluated the joint effect of prenatal arsenic, lead, and manganese on neurodevelopmental outcomes among 20-to 40-mo-old children in Bangladesh. This method allows for testing both the overall mixture effect and the effects of each mixture component within the context of the overall joint exposure. Further, the approach can identify the most important windows of vulnerability while accommodating highly correlated exposures using hierarchical variable selection. Finally, employing a Bayesian approach incorporates the uncertainty in the estimation of a high-dimensional set of exposures, accounting for multiple comparisons. Results of the analyses differed by clinic due to differences in exposure profiles. Using BKMR in this study population, we estimated a significant negative joint effect of the mixture of arsenic, manganese, and lead on CS that appears statistically significant when cord blood metal concentrations are all above the 60th percentile in the Pabna clinic (As >0:7 lg=dl, Mn >4:2 lg=dl, Pb >6:6 lg=dl, Figure 2). Further, we reported a) a marginally significant negative and additive linear association of lead with CS in the Sirajdikhan clinic; b) a negative and nonlinear association of manganese with CS within the overall mixture effect; and c) marginally significant interaction between manganese and arsenic on CS within the overall mixture (i.e., multiplicative), such that toxicity of manganese may increase more than additively in the presence of arsenic in the Pabna clinic. Our results suggest that arsenic may be a "potentiator" of manganese, i.e., its toxicity may be dependent to a large extent on the presence of manganese, and alone may not have strong main effects. While some results suggest a positive (although not statistically significant) slope of arsenic, we believe this positive effect may be due to an unmeasured confounder. For example, arsenic exposure sources include food and water, access to which may be a correlate of better nutrition in pregnancy, particularly in a poor country such as Bangladesh. Better nutrition in pregnancy would likely predict higher CS. Other unmeasured Table 3. Descriptive summary of results on the joint association of the metal mixture (arsenic, manganese, lead) and cognitive score. Environmental Health Perspectives 067015-8 confounders are possible. We will conduct exposure studies in future work to try to disentangle arsenic exposure and diet in pregnancy. Our overall finding of a negative joint effect of metal mixtures on neurodevelopment is consistent with studies indicating the importance of considering joint chemical exposure on pediatric health (Carpenter et al. 2002;Hertzberg and Teuschler 2002;Claus Henn et al. 2014). A unique feature of BKMR is that we are able to assess both additive and multiplicative components of a mixture simultaneously.
A strength of BKMR is that it not only addresses the overall mixture effect, but can also tease out the contributions of each component, with the caveat that these contributions are in the context of the joint exposure at the exposure levels seen in the cohort. The BKMR-grouped variable selection analyses indicate that in this sample, prenatal exposure to heavy metals is more predictive of 20-to 40-mo neurodevelopment relative to postnatal exposures measured at 20-40 mo. Our stratified analyses by clinic reveal that metals' relative concentrations and baseline population characteristics influence the metal mixture effect on child neurodevelopment. For example, we observed strong evidence of a manganese effect predicting lower CS with a suggestion of nonlinearity. Previous studies in lower-exposed populations indicated an inverted-U relationship between manganese and neurodevelopment, providing evidence for its beneficial effect at midrange levels and toxic effect at both lower and higher concentrations (Claus Henn et al. 2012). Although our study instead indicates a toxic nonlinear effect of manganese that flattens at higher concentrations ( Figure 2C), this is in the context of arsenic and lead exposures that are much higher than those observed in most populations. The context-dependent nature of a chemical's main effect is perhaps not surprising, given what we know about overlapping biological properties, but is rarely considered when comparing studies of the same chemical or in policy making. Indeed, this issue is likely critical in risk assessment, and may also explain variability in results among studies that address the main effects of neurotoxicants, but have variable effect estimates. The higher manganese levels in our study population, compared to U. S. populations, may mean Bangladesh is at a point in the doseresponse curve where a potential ceiling effect is seen and no further toxicity occurs. Alternatively, they may be due to the different contextual effect of higher population medians for arsenic and lead than that seen in Mexico or U.S. populations. Interestingly, the marginally significant association of lead with CS was found to be linear and additive in this population. The effects are clinically relevant, as the strength of the impact is similar in magnitude for other well-known risk factors adjusted for in the present study (e.g., maternal education, protein intake). These effects are also comparable to lead effects seen in other studies (Budtz-Jørgensen et al. 2013). Blood lead levels in Bangladesh are much higher than in developed populations, and our results are consistent with findings from highly exposed populations from the 1980s and 1990s. This might mean that lead exposure at these higher levels is primarily neurotoxic on its own, and does not interact with other chemicals. More mixtures research on lead is needed in lower-exposed populations in which lead's main effects are less prominent.
The neurotoxicity of manganese exposure was more pronounced in the Pabna population, where concentrations of manganese were highest. In this clinic, we estimated that the change in the log-transformed manganese from the 25th to the 75th percentile is associated with − 0:30 ð − 0:52, − 0:04Þ SDs in CS when lead and arsenic are held constant at their median levels. Note that Pabna displays a larger IQR for manganese than the overall study population. In addition, lead was not associated with scores among children in the Pabna clinic. As previously noted, we believe that in the context of very high Mn exposure, neurotoxicity has plateaued and the addition of other metals, such as arsenic or lead, does not produce any additional neurotoxicity. Consistent with this hypothesis, arsenic and lead appear to be neurotoxic in the Sirajdikhan clinic, where manganese concentrations are lower, although we note their effect estimates are not statistically significant.
In summary, our findings suggest that mixtures may induce different dose-response relationships in communities and that the size and shape of the dose-response curves are dependent on the population's overall metal exposure levels. This perhaps is not surprising, as interactions are likely to be dose-dependent. As an example, one may not expect interaction in the context of a severe poisoning because the toxicity of the high-dose chemical overwhelms the toxic properties of any concurrent levels of toxic chemicals (Könemann and Pieters 1996), and chemical interactions are less likely at high-dose exposures. Several other factors might contribute to differences in exposure profiles and in their impact on neurodevelopment across the two study sites. Heterogeneity might be due to geological differences in the bioavailability of metals at each site, or to cultural/behavioral factors influencing differential exposure to metals. Qualitative studies in the two populations revealed potentially different habits in accessing water sources, with the Sirajdikhan population being more likely to employ pond water to cook and clean, while the Pabna population is more likely to use well water. Differences between the two populations may be due to differences in other contextual factors as well. For example, we observed that participants from Sirajdikhan had lower protein intake, possibly due to their vegetarian habits. Finally, we note that our methodology (BKMR) is robust, and should be able to incorporate data on nonchemical toxicants as well as chemicals. Nutritional intake might also interact with the metals mixture; this will be the subject of future investigation and will be incorporated by BKMR.

Limitations
Umbilical cord blood was used as a biomarker for all metals to represent prenatal exposure. While cord blood will reflect late prenatal exposure to lead, other matrices may better reflect exposure to arsenic or manganese. Manganese levels are tightly regulated via homeostatic mechanisms in the blood; nonetheless, blood manganese is believed to be a reasonable indicator of environmental exposure and body burden in situations of chronic exposure. For arsenic, nail is considered the best biomarker for long-term arsenic exposure, but it was not available for the present analyses. Nonetheless, if our use of blood metals introduced exposure error via misclassification, given the 2 y that elapsed from exposure assessment to Bayley Scale assessment, any error would be expected to be nondifferential and would likely drive effects toward null values. We believe any measurement error might lead to underestimation of the toxic effect of any of the three metals. Loss to follow-up might induce bias in our analyses; however, it is unlikely to be dependent on child neurodevelopment. Potential residual confounding might lead to bias in the reported estimates. However, stratified analyses by clinic reduced heterogeneity in demographics that could act as confounders. Finally, measurement error in confounders might lead to inflated estimates.

Strengths
Metal concentrations display high variability and heterogeneous exposure profiles are found, across clinics, increasing the generalizability of the findings. We employed a novel, flexible statistical method, BKMR, and we were able to quantify and visualize the joint effect of the mixture and potential nonlinearities and nonadditive effects without coarsening the continuous exposures, thus reducing measurement error bias. Finally, by adopting this method, we overcome important limitations of traditional analyses, such as single metal effect estimation, model misspecification, and increased false discovery when fitting many regression models.

Conclusion
Coexposure to arsenic, manganese, and lead during gestation is shown to be jointly neurotoxic and to affect neurodevelopment of Bangladeshi children in early life stages in a complex fashion. Manganese at high concentrations was found to be the most neurotoxic component of the metal mixture. Strategies to monitor and prevent high exposure to manganese and other toxic metals are critical in the Bangladeshi population, especially among pregnant women. Heterogeneous findings across our two study sites that differ in sociodemographic characteristics motivate further investigation of child neurodevelopment, accounting for the complex mixtures of environmental, nutritional, and social risk factors that shape mothers' and children's development. Our study is among the first designed to address the effects of mixtures, and our results suggest that BKMR may be a valuable tool in larger-scale mixture studies. As this research field progresses, we plan to utilize these methods on other chemicals and nonchemical toxicants, and perhaps ultimately employ BKMR for research in the nascent field of "exposomics."