Associations of a Metal Mixture Measured in Multiple Biomarkers with IQ: Evidence from Italian Adolescents Living near Ferroalloy Industry

Background: Research on the health effects of chemical mixtures has focused mainly on early life rather than adolescence, a potentially important developmental life stage. Objectives: We examined associations of a metal mixture with general cognition in a cross-sectional study of adolescents residing near ferromanganese industry, a source of airborne metals emissions. Methods: We measured manganese (Mn), lead (Pb), copper (Cu), and chromium (Cr) in hair, blood, urine, nails, and saliva from 635 Italian adolescents 10–14 years of age. Full-scale, verbal, and performance intelligence quotient (FSIQ, VIQ, PIQ) scores were assessed using the Wechsler Intelligence Scale for Children-III. Multivariable linear regression and Bayesian kernel machine regression (BKMR) were used to estimate associations of the metal mixture with IQ. In secondary analyses, we used BKMR’s hierarchical variable selection option to inform biomarker selection for Mn, Cu, and Cr. Results: Median metal concentrations were as follows: hair Mn, 0.08μg/g; hair Cu, 9.6μg/g; hair Cr, 0.05μg/g; and blood Pb, 1.3μg/dL. Adjusted models revealed an inverted U-shaped association between hair Cu and VIQ, consistent with Cu as an essential nutrient that is neurotoxic in excess. At low levels of hair Cu (10th percentile, 5.4μg/g), higher concentrations (90th percentiles) of the mixture of Mn, Pb, and Cr (0.3μg/g, 2.6μg/dL, and 0.1μg/g, respectively) were associated with a 2.9 (95% CI: −5.2, −0.5)–point decrease in VIQ score, compared with median concentrations of the mixture. There was suggestive evidence of interaction between Mn and Cu. In secondary analyses, saliva Mn, hair Cu, and saliva Cr were selected as the biomarkers most strongly associated with VIQ score. Discussion: Higher adolescent levels of Mn, Pb, and Cr were associated with lower IQ scores, especially at low Cu levels. Findings also support further investigation into Cu as both beneficial and toxic for neurobehavioral outcomes. https://doi.org/10.1289/EHP6803


Introduction
Recent progress in understanding how chemical mixtures affect health can be largely attributed to collaboration between scientific fields, innovation of statistical methods, and prioritization from funding agencies (Carlin et al. 2013;Taylor et al. 2016). Metals are an important class of chemicals in which to conduct mixtures research given that they are ubiquitous, commonly co-occur in the environment, and have interactive potential (Claus Henn et al. 2014;Valeri et al. 2017;Wright et al. 2006). Communities are exposed to mixtures of metals that may be both naturally occurring and from anthropogenic sources. Anthropogenic activities such as gasoline combustion, waste incineration and recycling, mining, fungicide application, and steel industry emissions, including those from ferroalloy plants, are known to emit metals into local environments (ATSDR 2012).
Manganese ferroalloys are metal mixtures primarily composed of iron (Fe) in combination with manganese (Mn) and other metals such as chromium (Cr), lead (Pb), copper (Cu), nickel (Ni), zinc (Zn), cadmium (Cd), and aluminum (Al) (Olsen et al. 2007). Most environmental epidemiology studies in communities living near ferromanganese industries have focused on Mn only, reporting associations with parkinsonism in adults and neurobehavioral decrements in children and the elderly, such as worse overall cognition, behavior, attention, learning and memory, olfactory function, and motor abilities (Bauer et al. 2017;Chiu et al. 2017;Haynes et al. 2015;Lucchini et al. 2007Lucchini et al. , 2014Menezes-Filho et al. 2014;Rodrigues et al. 2018). Yet other metals in the ferroalloy mixture, such as Pb and Cu, are also known to be neurotoxic (Bulcke et al. 2017;Caito and Aschner 2017). Furthermore, prior studies have shown interactions between some of these metals on neurodevelopment (Claus Henn et al. 2012;Gorell et al. 1997;Lin et al. 2013;Liu et al. 2018;Menezes-Filho et al. 2018). Despite the presence of multiple metals in ferroalloy emissions that may have additive or synergistic associations with neurobehavioral outcomes, few prior studies have considered the health effects of exposures to multiple metals (Lucchini et al. 2012b;Menezes-Filho et al. 2018).
Neurobehavioral effects from exposure to environmental contaminants, such as metals, are understudied in adolescence and early adulthood although these are potentially important developmental life stages (Wasserman et al. 2018). Neurodevelopment occurs as a dynamic process from the in utero period to adulthood (e.g., ∼ 30 years of age) (Blakemore and Choudhury 2006;Sowell et al. 2001). In addition to the prenatal and early childhood periods, preadolescence and adolescence are periods of rapid growth and development during which the brain is sensitive to neurotoxicants. During this time frame, the brain changes dramatically via synaptic pruning, myelination, neuronal transmission, and neural circuitry restructuring (Blakemore and Choudhury 2006;Giedd et al. 1999;Wahlstrom et al. 2010). These processes lead to refinement in brain areas, particularly the prefrontal cortex and hippocampus, which are important for cognitive abilities such as higher-order executive functioning (Sander et al. 2012).
In this article, we examine the association of a mixture of Mn, Pb, Cu, and Cr with cognition in adolescents living near varied ferromanganese industrial activity that utilizes and emits these metals. Our study investigated this mixture association during adolescence, a potentially critical window for neurodevelopment. Data analysis included both traditional regression and Bayesian kernel machine regression (BKMR)-a statistical approach developed to assess the health effects of environmental mixtures (Bobb et al. 2015)-to estimate individual, joint, and interactive associations. As a secondary analysis, given the lack of scientific consensus on the best exposure biomarker for Mn, Cu, and Cr, the hierarchical variable selection option in BKMR was used to help inform the choice of metal biomarker.

Study Participants
Participants for this analysis were part of the ongoing Public Health Impact of Mixed Element exposure (PHIME) study, designed to investigate associations between multiple metal exposures from ferroalloy emissions and neurobehavior in adolescents 10-14 years of age. We have previously described the study design in detail (Lucchini et al. 2012a). Briefly, 720 participants were recruited from three demographically similar, but geographically distinct, sites in the province of Brescia, Italy, a region with varied ferroalloy activity: a) Bagnolo Mella, an area with currently active ferroalloy industry since 1974; b) Valcamonica, an area with historical ferroalloy production for over a century that ended in 2001; and c) Garda Lake, a tourist region with no history of ferroalloy activity (see Lucchini et al. 2007 for a map).
Participants were eligible for the study if a) they were living in the study area since birth; b) their families had lived in the study area since 1970; and c) their children were 10-14 years of age at the time of enrollment. Participants were excluded from the study if they a) had been diagnosed with neurological, hepatic, metabolic or endocrine diseases, or clinically evident hand/ finger motor deficits; b) were currently prescribed psychoactive drugs or had psychiatric disturbances; c) had inadequately corrected visual deficits; or d) had ever received parenteral nutrition.
We enrolled 720 participants into the study: 312 participants in the first phase (2007)(2008)(2009)(2010) and 408 in the second phase (2010)(2011)(2012)(2013)(2014). The two phases reflect two waves of funding, but they were conducted by the same researchers using identical questionnaires and study protocols. The second phase included a third site (Bagnolo Mella), as well as measurement of the Home Observation Measurement of the Environment (HOME) Short Form, a measure of cognitive stimulation and emotional support in the home (National Longitudinal Surveys 1979). The second phase also included measurement of metals concentrations in additional biomarkers (saliva, urine, fingernails). Eligible children received a detailed description of the study procedures before consenting to participate. The institutional review boards at the Ethical Committee of the Public Health Agency of Brescia, the Icahn School of Medicine at Mount Sinai, and the University of California, Santa Cruz approved all PHIME study protocols.

Metals Concentrations in Biological Samples
Blood, hair, fingernails, urine, and saliva were collected from participants at 10-14 years of age, concurrent with neurobehavioral assessment. Methods of collection and analysis of biological samples have been previously described (Eastman et al. 2013;Lucchini et al. 2012a;Smith et al. 2007). Briefly, we collected venous whole blood samples using a 19-gauge butterfly catheter into lithium-heparin Sarstedt Monovette ® Vacutainers; hair samples using stainless steel scissors (2-3 cm section of hair from the nape of the neck, proximal to the scalp); fingernails using stainless steel nail clippers; spot urine samples in polyethylene containers; and passive saliva samples through a plastic straw into trace metal-free microfuge tubes after rinsing the mouth three times with ultrapure water and waiting 10 min. Hair and fingernail samples were cleaned of exogenous metal contamination using a validated method (Eastman et al. 2013;Lucas et al. 2015).
Concentrations of Mn, Pb, Cu, and Cr were measured in all samples using magnetic sector inductively coupled plasma mass spectrometry (Eastman et al. 2013;Lucchini et al. 2012aLucchini et al. , 2012bSmith et al. 2007). The analytical detection limits for all biomarkers were based on repeated measurements of procedural blanks on four separate analysis days (for details see Butler et al. 2019). Most of the biological measurements were above the limits of detection (LODs). Samples below the LOD were assigned a value of one-half the LOD (hair Mn, n = 1; Mn urine, n = 5; Mn nails, n = 25; hair Cr, n = 1; saliva Cr, n = 1; Cr nails, n = 31; Cu nails, n = 2). education and occupation. We categorized participants into low, medium, or high socioeconomic status (SES) based on methodology developed in Italy that combines parent education and occupation (Cesana et al. 1995;Lucchini et al. 2012b). Hemoglobin (in grams per deciliter) was measured in the same blood samples collected for metals analysis and was considered a proxy for iron status, which may influence Mn and Pb absorption and neurodevelopment (Park et al. 2013;Smith et al. 2013).

Statistical Analysis
We performed univariate and bivariate analyses for each variable. Metals concentrations were natural log (ln) transformed to satisfy model assumptions of normality of residuals and then standardized. Bivariate correlations between metals concentrations were examined using Spearman's rank correlation coefficients. We used age-adjusted WISC IQ scores (FSIQ, VIQ, and PIQ) and modeled them as continuous outcome variables.
There is a lack of consensus about which biological sample is the most consistent and valid biomarker of exposure for Mn, Cu, and Cr (Bertinato and Zouzoulas 2009;Coetzee et al. 2016;Jursa et al. 2018;Lukaski 1999;Viana et al. 2014). In primary analyses, we considered concentrations of these three metals in hair because hair has been used for these metals in previous environmental epidemiologic studies (Bertinato and Zouzoulas 2009;Bouchard et al. 2011;Haynes et al. 2015;Randall and Gibson 1989;Riojas-Rodríguez et al. 2010) and had the largest sample size in our data. In secondary analyses, we used a data-driven approach (described below) to select the biological sample to use in models with IQ. Concentrations of Pb measured in blood were used for all analyses, given its wide acceptance as a Pb biomarker.
Confounder selection. We used directed acyclic graphs to examine potential confounders that might be associated with metals concentrations and cognition but are not causal intermediates or colliders (see Figure S1). Continuous covariates and exposures were modeled in a way that best represented the shape of the association with IQ (e.g., HOME score was modeled with an additional quadratic term) based on visual assessment of penalized splines (constrained to 4 knots) in adjusted generalized additive models (GAMs). All final models were adjusted for sex, SES (high, medium, low), hemoglobin (in grams per deciliter), HOME score, and HOME score squared (due to a nonlinear relationship with IQ).
Multiple imputation. Some participants were missing data on the HOME score and some biomarkers. We used multiple imputation to impute missing values using chained equations with the MICE package in R (van Buuren and Groothuis-Oudshoorn 2011).
We imputed missing values for all participants (n = 709) (White et al. 2011) and included in the imputation process variables thought to be related to the missing data, including outcomes, metals concentrations in biological as well as environmental samples (Lucchini et al. 2012a), and potential confounders (see Table  S1 for list of variables used in imputation). We assumed data were missing at random, that is, that the missingness did not depend on unobserved data. We generated 40 imputed data sets. In our analytic data set, we included only participants with measured IQ and hair metals concentrations (n = 635).
GAMs and multivariable linear regression models. We first examined the shape of the ln-metals (hair Mn, hair Cu, hair Cr, blood Pb) and IQ associations using GAMs with penalized splines (knots = 4), adjusting for previously selected confounders. In GAMs, to visualize the association with the outcome, we used standard multiple imputation methods to combine the point estimates of coefficients parameterizing the nonlinear terms from multiple GAM fits across imputed data sets, averaging results from the 40 imputations. We calculated standard errors using Rubin's rule that combines the within-and between-imputation variances of these coefficient estimates (Rubin 2004). To plot the association between each ln-metal and IQ, we estimated nonlinear effect estimates across each metal's range of values in 0.1 increments while setting all other continuous variables in the model at their median and as linear combinations of these averaged coefficient estimates and computed point-wise 95% confidence intervals (CIs) of these linear combinations accordingly (see Figure S2).
Based on adjusted GAMs, all ln-metals appeared to be linearly associated with IQ except Cu, which appeared to have an inverted U-shaped association. To test nonlinearity, we compared model fit using a likelihood ratio test from the MICE R package (van Buuren and Groothuis-Oudshoorn 2011) for multivariable linear models. We compared model fit between models with an additional quadratic term for ln-Cu to a nested model with a linear ln-Cu term only.
To obtain effect estimates, we subsequently used multivariable linear regression, adjusting for selected confounders with Mn, Pb, and Cr modeled as continuous ln-transformed linear terms and ln-Cu modeled as a categorical variable (tertiles: high, low versus middle). We initially included all pairwise interactions between metals as cross-product terms in this model, but herein we report estimates from the final model only, which included only statistically significant (p < 0:10) pairwise interactions. Our final multivariable linear regression model was specified as follows: We obtained pooled estimates and CIs by combining information from the 40 mean and variance estimates using Rubin's method (Rubin 2004).
This approach of multivariable regression modeling, however, limits the investigation of the mixture to pairwise interactions. We therefore considered BKMR models to further examine this potentially complex relationship.
Bayesian kernel machine regression. To allow for potential nonlinear associations between ln-metals and IQ as well as higher-order interactions between metals, we used Bayesian kernel machine regression (BKMR) to flexibly model the adjusted associations of the metal mixture with IQ. This method has been described previously in detail (Bobb et al. 2015Valeri et al. 2017). Briefly, this approach estimates the joint association of mixture components and a health outcome such as IQ, allowing for nonlinear and nonadditive associations and potentially higher-order interactions among correlated mixture components, while adjusting for confounders. We used the component-wise variable selection option of BKMR, which controls for multiple testing and provides a relative measure of variable importance for each individual component of the mixture. The BKMR model for our analyses is given as follows: The model consists of a health outcome Y i ; a smooth function hðÞ that is modeled by a kernel function K that links the exposure contribution to the outcome for different subjects by the distance between individuals in the exposure space; a vector X i that allows for adjustment of confounders; and a random error term e i . All metals were ln-transformed to remain consistent with the multivariable linear regression models. The choice of a Gaussian kernel function permits hðÞ to potentially contain nonlinear and interactive effects among mixture components, allowing the data to dictate the shape of this exposure-response surface. For primary analyses, we obtained 10,000 posterior samples of all model parameters using a Markov chain Monte Carlo sampler and default noninformative priors specified in the R package. Credible intervals for estimates, which account for uncertainty both due to estimating the exposure-response function and due to variable selection from among a high-dimensional set of exposures, were calculated using the bkmr R package . The bkmr package also calculates posterior inclusion probabilities (PIPs) to quantify the relative importance of each metal to the outcome.
We fit a separate BKMR model for each of the 40 imputed data sets. Using each BKMR fit, we estimated the posterior mean and variance for every contrast of interest. For each contrast, we obtained an overall estimate and credible interval by combining information from the 40 posterior mean and variance estimates using Rubin's method. We developed publicly available R functions to combine information from multiple BKMR fits and to provide estimates of the overall association and credible intervals (Devick 2019) (see code file in supplemental material).
In order to describe associations from BKMR's highdimensional parameter space, we computed summary measures to quantify and visualize the exposure-response surface . We estimated the following: a) the exposureresponse relationship for each metal (e.g., Mn) on IQ when all other metals (i.e., Pb, Cu, and Cr) are fixed at their median; b) the association of an interquartile range (IQR; i.e., 25th to 75th percentile) increase in a particular metal (e.g., Mn) on IQ, at varying levels of the mixture (i.e., at the 25th, 50th, and 75th concurrent percentiles of Pb, Cu, and Cr); c) the pairwise dose-response relationship for each metal (e.g., Mn) on IQ at varying levels of a second metal (e.g., Cu at 25th, 50th, and 75th percentiles), when remaining metals (i.e., Pb and Cr) are fixed at their median; and d) the joint association of an incremental increase in all metals on IQ, compared with when all metals are their median.
Hierarchical variable selection option in BKMR. Given the lack of consensus about which biological sample is the most consistent and valid biomarker of exposure for Mn, Cu, and Cr (Bertinato and Zouzoulas 2009;Coetzee et al. 2016;Jursa et al. 2018;Lukaski 1999;de Sousa Viana et al. 2014), we explored the use of hierarchical, or grouped, variable selection in BKMR as a secondary analysis to identify the biological sample for each metal that was most strongly associated with IQ. Details on the hierarchical variable selection approach of BKMR have been published previously (Bobb et al. 2015). Briefly, this function allows for the selection of one component within a prespecified group of correlated components. These secondary analyses proceeded in two steps. In the first step, we applied the hierarchical variable selection function to identify the most predictive biological sample type for each of the three metals (Mn, Cu, and Cr) when considering IQ as the outcome. In the second step, we fit new BKMR models with these newly selected biomarkers (i.e., one biological sample type per metal) along with blood Pb to evaluate associations with IQ. Specifically, we first entered the suite of biomarkers into the model as a group for each of the three metals (e.g., for Mn: blood Mn, hair Mn, saliva Mn, urine Mn, nail Mn). The hierarchical variable selection process quantifies the relative importance of each group (i.e., each metal) as well as each biological sample type within each group to the outcome using a PIP. For example, for Mn, conditional PIPs are estimated for blood Mn, hair Mn, saliva Mn, urine Mn, and nail Mn, given that the group of Mn biomarkers is an important group of exposures in the model. We selected the biological sample type with the highest conditional PIP to use, as a second step, in a standard BKMR model. This final model included only the selected biomarkers for Mn, Cu, and Cr as well as for blood Pb.
Sensitivity analyses. In sensitivity analyses we examined the robustness of our findings in three ways: a) We used two alternative prior specifications for BKMR models to allow for varying degrees of smoothness of the function (Valeri et al. 2017); b) we examined the potential confounding effect of study site by adding this variable to models; and c) we performed analyses using complete data (i.e., without multiple imputation). All statistical analyses were conducted in R (version 3.5.1; R Development Core Team). BKMR was run using the bkmr package .

GAMs and Multivariable Linear Regression
We began by fitting GAMs and linear regression models for VIQ. We observed a nonlinear, inverted U-shaped association between ln-hair Cu and VIQ (see Figure S2). Using multivariable linear regression, the pooled estimate from 40 models with an additional quadratic term for ln-Cu was a better fit (likelihood ratio test p = 0:03) compared with the pooled estimate of nested models with ln-Cu specified only as a linear term, supporting the observed nonlinear Cu-VIQ association. We then fit multivariable linear regression models with ln-Cu categorized into tertiles to obtain more interpretable estimates while maintaining the inverted U-shaped association. Compared with the middle tertile, low and high Cu were associated with 1.2 (−3:7, 1.2) and 2.2 (−4:7, 0.3)-point decreases in VIQ after adjusting for sex, age, SES, hemoglobin, HOME score, and HOME score squared as well as Mn, Pb, and Cr (Table 3). There was evidence of interaction between ln-Mn and ln-Cu (interaction term p-value between ln-Mn and low ln-Cu: p interaction = 0:04): a 1-SD increase in Mn was associated with decreased VIQ scores among individuals with low or mid-levels of Cu (b 1 + b 6 = − 2:0 for low Cu; b 1 = − 1:4 for mid-level Cu), whereas Mn was associated with a 1.1-point increase in VIQ among individuals with high Cu levels (b 1 + b 7 ). Similarly, the adverse ln-Pb association was associated with a 2.1-point decrease in VIQ among individuals with low ln-Cu levels (b 2 + b 8 = − 2:14), whereas the Pb association was close to null among those with mid-and high-level Cu. For PIQ, the association with Cu did not appear to be nonlinear and pairwise interactions with Mn or Pb were not evident (p interaction = 0:45-0:93) (Table 3). Cu was nonlinearly associated as an inverted U with FSIQ, but there was no evidence of pairwise interaction with Mn or Pb (p interaction = 0:23-0:77).

Bayesian Kernel Machine Regression Analyses
We implemented BKMR to obtain an estimate of the joint exposure-response function of hair Mn, hair Cu, hair Cr, and blood Pb on IQ in adolescents. We began by investigating the doseresponse relationships of each individual (ln-transformed) metal in the mixture with VIQ. Figure 1 displays the exposure-response functions for each metal, when all other metals are set at their medians, after adjusting for sex, age, SES, hemoglobin, HOME score, and HOME score squared. Similar to findings from GAMs and multivariable linear regression, findings from the BKMR analyses suggested a nonlinear, inverted U-shaped association between ln-Cu and VIQ, where both low and high concentrations of Cu were associated with lower VIQ scores when other metals were set to their medians. A decrease in Cu from the 40th to 10th percentile (8.5 to 5:4 lg=g; Table 2) was associated with a 0.8 (−1:9, 0.3)-point reduction in VIQ score, whereas an increase in Cu from the 60th to 90th percentile (11.8 to 24:9 lg=g) was associated with a 1.3 (−2:7, 0.2)-point reduction in VIQ. The association between Cr and VIQ was null when other metals were set to their medians. To more fully investigate possible effect modification by Cu in light of its inverted U-shaped association, we estimated the association of an IQR increase in an individual metal on VIQ score, at varying levels (25th, 50th or 75th percentiles) of all three other metals and stratified by level of Cu (at the 10th, 50th, and 90th percentiles) (Figure 2). Interaction was suggestive if there was variation in the effect estimate of a given metal at different percentiles of the other metals. There was variation between estimates for each metal across the three different plots of Cu percentiles, although CIs were relatively wide. For example, when Cu was set at its 10th, 50th, or 90th percentile (5.5, 9.6, or 24:9 lg=g, respectively), an IQR increase in Mn (0:10 lg=g) was associated with respective reductions in VIQ of approximately 1.3 (−3:0, 0.3), 1.0 (−2:5, 0.4) and 0.5 (−2:2, 1.1) points, regardless of Pb and Cr variation. This suggests that the adverse associations of Mn are more apparent at lower Cu levels. This finding is also reflected in Figure 3, showing  Table 3. Adjusted associations of verbal, performance, and full-scale intelligence quotient (VIQ, PIQ, and FSIQ) and metal biomarkers from multivariable linear regression models.  and VIQ, when other metals are set to their respective medians. Model was adjusted for age, sex, SES, HOME score, quadratic HOME score, and hemoglobin. Note: HOME, Home Observation Measurement of the Environment; SES, socioeconomic status; VIQ, verbal intelligence quotient.
pairwise interactions, where the negative slope of ln-Mn with VIQ is slightly steeper at lower (e.g., 25th percentile) Cu levels than at higher Cu levels. Associations of ln-Cr and VIQ were null (omitted from figure). Finally, we investigated the overall joint association of the metal mixture on VIQ score, stratified by Cu level (Figure 4). Each posterior mean estimate and the corresponding CI represented the change in VIQ when each metal in the mixture was concurrently set to a particular percentile (10th to 85th percentile) as compared with when each metal was set to its median value (50th percentile). Overall, increasing levels of Mn, Pb, and Cr were associated with decreases in VIQ score, especially at low levels of Cu. For example, at low Cu levels (10th percentile), metals set at their respective 70th, 75th, and 80th percentiles were associated with 1.1 (−1:9, −0:2), 1.4 (−2:6, −0:2), and 1.8 (−3:2, −0:3)-point decreases in VIQ, respectively, compared with when all metals were set at their 50th percentiles.
In models of FSIQ and PIQ (see Figures S3 and S4), Mn and Pb associations were also negative and the Cu association also appeared to be nonlinear, but associations were smaller in magnitude and less precise (see Figures S3A and S4A). Associations of the overall mixture with FSIQ and PIQ were negative but weaker than associations with VIQ (see Figures S3D and S4D).

Hierarchical Variable Selection Using BKMR
We implemented hierarchical variable selection within BKMR to allow the model to select, for each metal, the biological sample that was most predictive of the VIQ score. We used 10 randomly selected data sets from the 40 imputed data sets. For (ln-transformed) Mn, Cu, and Cr (blood was forced into the model as the Pb biomarker), the highest average conditional PIPs (cPIPs) across the 10 data sets were estimated for saliva for Mn (cPIP = 0:99), hair for Cu (cPIP = 0:60), and saliva for Cr (cPIP = 0:74). In BKMR models using these selected biomarkers plus blood Pb, saliva Mn was negatively associated with VIQ ( Figure 5A,B) at all levels of the other metals: An IQR increase in saliva Mn (71:4 lg=L) was associated with reductions of 5.1 (95% CI: −8:6, −1:7), 4.8 (−7:9, −1:7), and 4.2 (−7:7, −0:8) points in VIQ score, when other metals were at their 25th, 50th, and 75th percentiles, respectively ( Figure 5B; Cu results were omitted from Figure 5B due to their nonlinear shape). There was also a suggested interaction between Mn and at least one other metal: Mn effect estimates varied slightly with increasing levels of the mixture, with the strongest association estimated when all other metals were set at their 25th percentile ( Figure 5B). Jointly increasing levels of all metals were associated with decreases in VIQ score, which was stronger than the joint association observed previously in models using hair as the biomarker of exposure for Mn, Cu, and Cr (cf. Figure 5D versus Figure 4).

Sensitivity Analyses
Results from complete case analyses for both primary analyses (using hair Mn, hair Cu, hair Cr, and blood Pb; n = 338) and secondary analyses (using saliva Mn, hair Cu, saliva Cr, and blood Pb; n = 330) were similar to models using multiple imputation, although estimates were attenuated (see Figures S5 and S6). Adjusting for study site as an additional confounder in models did not appreciably change results, nor did changing the prior specifications to allow for varying degrees of smoothness for the parameter of the h function in the BKMR model (see Figures S7-S9). Overall, our conclusions from primary analyses were unchanged after performing sensitivity analyses.
We found that joint increases in metal concentrations were associated with lower IQ scores but particularly at lower Cu levels (e.g., fixed at its 10th percentile) and when other metals were above their respective 50th percentiles. Associations with VIQ were stronger than with PIQ. Associations of (ln-transformed) Mn, measured either in hair or saliva, with VIQ were linear and negative, suggesting adverse effects on aspects of cognition that involve language comprehension, verbal expression, attention, working memory, and executive function. Although other studies of children 7-12 years of age with industrial airborne Mn exposure did not estimate the association of Mn on VIQ in the context of a metals mixture, they have similarly reported adverse, linear associations between hair Mn and VIQ (Menezes-Filho et al. 2011;Riojas-Rodríguez et al. 2010;Wright et al. 2006) as well as with domains that are influenced by verbal abilities, including verbal learning, memory Torres-Agustín et al. 2013;van Wendel de Joode et al. 2016;Wright et al. 2006), and verbal working memory Haynes et al. 2018;van Wendel de Joode et al. 2016). These findings are also consistent with some (Bouchard et al. 2011;Wasserman et al. 2011), but not all , pediatric studies of environmental Mn exposure via drinking water and VIQ. Inconsistencies among the drinking-water studies might be due in part to reduced variability of Mn exposure between subjects and differences in outcome measurement given that some of the studies used a condensed neurobehavioral assessment battery to calculate VIQ ).
An inverted U-shaped association was observed between ln-hair Cu and VIQ, suggesting that both low and high levels of Cu may Figure 5. The adjusted association of the metal mixture [manganese (Mn), lead (Pb), copper (Cu) and chromium (Cr)] on VIQ after using hierarchical variable selection to identify most associated biomarkers of exposure (saliva Mn, blood Pb, hair Cu, saliva Cr). Models were adjusted for age, sex, SES, HOME score, quadratic HOME score, and hemoglobin. (A) Exposure-response functions for each metal (Mn, Pb, Cu, or Cr) and VIQ, when other metals are set to their respective medians. (B) Associations (estimates and 95% credible intervals) for an IQR increase in each metal with VIQ, when other metals are fixed at either the 25th, 50th, or 75th percentile. IQRs for metals are: saliva Mn (57:5 lg=L); saliva Cr (2:2 lg=L); blood Pb (0:9 lg=dL). The 25th, 50th, and 75th percentiles for metals: saliva Mn (3.5, 16.5, 75:0 lg=L); blood Pb (1.0, 1.3, 1:9 lg=dL); saliva Cr (0.28, 0.82, 2:43 lg=L). Cu results in B omitted due to nonlinear shape. (C) Exposure-response functions for Cu, Mn, and Pb at varying levels (25th, 50th, 75th percentile) of another metal, when other metals are set at their median [hair Cu (7.05,9.56,15:37 lg=g)]. (D) The overall association of the mixture of Mn, Pb, Cu, and Cr on VIQ (estimates and 95% credible intervals). The figure plots the estimated change in VIQ when all metals are at a given percentile (x-axis) compared with when all metals are at their 50th percentile. Note: expos, exposure; HOME, Home Observation Measurement of the Environment; IQR, interquartile range; SES, socioeconomic status; VIQ, verbal intelligence quotient. have adverse effects on neurobehavior. This nonlinear association is consistent with the properties of Cu as both an essential micronutrient and, in excess, a neurotoxicant (Otten et al. 2006). There are few environmental studies of the association of Cu and neurobehavior, although Cu contamination in the environment is common: Cu is mined extensively in the United States to produce a variety of products (i.e., electrical wiring, sheet metal, preservatives), is widely used as a pesticide, and is a common drinking-water contaminant (ATSDR 2004;U.S. EPA 2008). In addition, environmentally elevated Cu levels have been identified at more than half of 1,647 current or former U.S. Superfund sites (ATSDR 2004).
Animal behavioral research has demonstrated that increasing Cu dose by oral gavage or drinking water reduces spatial memory, attention, learning, motor coordination, and strength (Behzadfar et al. 2017;Kalita et al. 2018;Kumar et al. 2015). In humans, cross-sectional studies have reported higher Cu accumulation in hair and nails of children with autism spectrum disorder (ASD) compared with typically developing children (Lakshmi Priya and Geetha 2011;Obrenovich et al. 2011). Yet few epidemiologic studies have investigated preclinical decrements of neurobehavior in relation to Cu concentrations. Among 606 Belgian adolescents living near industrial areas, whole blood Cu was inversely related to sustained attention and working memory (Kicinski et al. 2015). Likewise, a cross-sectional study of 826 Chinese children 10-14 years of age reported associations of high total serum Cu (≤110 lg=dL) and reduced working memory among boys only (Zhou et al. 2015). In an ecological study of fourth-graders in New Orleans, Louisiana, soil Cu was associated with lower academic performance (Zahran et al. 2012).
On the other hand, Cu is biologically essential, and there is evidence that Cu supports brain development and possible neuroprotection (Bica et al. 2014;Hung et al. 2012;Kaler 2011). In vitro studies have demonstrated neuroprotective effects of Cu, in particular that Cu enhances neurite elongation (Bica et al. 2014) and improves Parkinson's disease outcomes (Hung et al. 2012). At nutritionally insufficient low levels, Cu has been linked with adverse health outcomes in humans, including Menkes disease, a rare X-linked recessive disorder where the enzyme ATP7A is absent, leading to severe muscle dystonia, gray matter neurodegeneration, and infantile death (Kaler 2011). Collectively, the aforementioned research supports both a neuroprotective and neurotoxic role for Cu, consistent with our finding of an inverted U-shaped association between hair ln-Cu and VIQ in GAMs, multivariable linear regressions, and BKMR models.
When examined on its own, Cu appears to be optimal for VIQ at midrange concentrations in our data. However, when examined in the context of a mixture, Cu may mitigate some of the neurotoxicity of Mn and Pb. In our study, we found suggestive interactions between Mn and Cu and between Pb and Cu, where the adverse associations of Mn with VIQ were strongest at lower Cu levels (i.e., 10th percentile, or 5:5 lg Cu/g hair). This suggests an opposing effect of Cu: Mn and Pb toxicity is most apparent when Cu, which has beneficial effects, is at its lowest levels. At higher levels of Cu (e.g., 50th or 90th percentiles in BKMR models or third tertile in multivariable linear regression models), the adverse associations of Mn and Pb are less pronounced. These findings underline the importance of considering metals in the context of a mixture.
Interactions between Mn and Cu as well as Pb and Cu are plausible given the potential for overlapping toxicological mechanisms. For example, Mn and Cu are important for brain growth and development but in excess they can both cause cellular oxidative stress and degeneration in the brain as well as induce parkinsonianlike symptoms (Ala et al. 2007;Burton and Guilarte 2009;D'Ambrosi and Rossi 2015). In addition, Pb, Mn, and Cu can disrupt calcium-dependent pathways and cellular redox processes (Burton and Guilarte 2009;Lidsky and Schneider 2003;Scheiber et al. 2014;Tjalkens et al. 2006). However, the epidemiologic literature on interactions between Cu and either Mn or Pb is sparse. One recent epidemiologic study of children in Mexico City used BKMR and reported a similar interaction between Cu and Pb and neurobehavior, where higher concentrations of Pb reduced the beneficial effect of Cu on infant neurodevelopment . Despite several differences in study design, including exposure timing (prenatal versus adolescent), age at outcome assessment (6-24 months versus adolescence), and sources of contaminants (unknown versus ferroalloy industry), this finding of a Cu-Pb interaction using BKMR is consistent with ours.
Some evidence exists for interaction between Cu and Mn in experimental studies. In rodents and nonhuman primates, Mn exposure has altered endogenous Cu homeostasis, whereby Mn exposure either reduced (Fu et al. 2015;Robison et al. 2013) or enhanced Cu brain tissue concentrations (Guilarte et al. 2006(Guilarte et al. , 2008Zheng et al. 2009). Furthermore, Mn and Cu compounds both inhibited NMDA-mediated receptors individually (Guilarte and Chen 2007). From these studies, there appears to be potential for these metals to influence one another's concentrations in the body as well as their potential effects on the brain.
This study used BKMR with hierarchical variable selection in order to guide the selection of exposure biomarkers. In analyses using this approach, saliva Mn was most strongly associated with VIQ score. Saliva Mn was imputed for a large portion of our data (44%); yet a sensitivity analysis using participants with complete data for saliva (n = 330), identified saliva Mn as the biomarker most strongly associated with VIQ (see Figure S6). In addition, we sought to examine the influence of extreme values of saliva Mn in our full sample; however, no outliers were identified by the generalized extreme Studentized deviate (ESD) many-outlier procedure using imputed data sets (Rosner 1983). Although research using saliva as a biomarker for Mn exposure is sparse, recent studies describing analytical methodologies for measuring saliva Mn in humans have reported associations between environmental Mn exposure and saliva Mn levels (Butler et al. 2019;Lucas et al. 2015) and have quantified saliva Mn in occupationally exposed workers (Gil et al. 2011;Wang et al. 2008). Among Mn-exposed welders, studies investigating saliva Mn have reported significant correlations with serum Mn (r = 0:57), airborne Mn (r = 0:65), and years of work experience (r = 0:40) (Fan et al. 2017;Wang et al. 2008;Zhou et al. 2010). In environmental studies of children, saliva Mn was significantly associated with air and soil Mn in areas near ferroalloy activity in the same cohort as the present study (Butler et al. 2019), but not with Mn in drinking water in a Canadian cohort Ntihabose et al. 2018). This latter discrepancy may be attributed to analytical differences in saliva sample processing between studies. These findings warrant further study of saliva as a biomarker of Mn exposure.
There are limitations to this study. As with any crosssectional analysis, we are unable to establish temporality between exposure and outcome, thus precluding our ability to evaluate causality. The temporality of metals and neurological disorders is a common concern given that some neurological disorders, such as ASD, could influence both IQ and metals concentrations. For example, metal concentrations in hair among patients with ASD could be influenced by altered metal metabolism, diet, behavior, or medication use from ASD rather than being a causal contributor to the etiology of ASD (Kalkbrenner et al. 2014). However, participants in this study were generally healthy adolescents with no diagnosed neurological disorder. Selection bias could be possible if participation in the study was based on a factor related both to exposures and outcome. For example, it is possible that subjects who perceived they were at greater risk of exposure to environmental metals or at increased risk of health effects from environmental exposures were more likely to participate in the study, although we think this unlikely. In addition, a diagnosed neurobehavioral disorder was an exclusion criterion in the study. We cannot rule out unmeasured confounding by zinc or other coexposures, which may be associated with levels of Mn, Pb, or Cu, as well as with neurodevelopment (Adamo and Oteiza 2010;Mahaffey 1990). In addition, both prenatal and postnatal metals exposure could be associated with neurodevelopment and correlated with adolescent exposure and thereby have acted as potential confounders in our study. In a subsample of our cohort (n = 195), Mn levels were measured in deciduous teeth as a proxy for prenatal and early postnatal metals exposure. We found that tooth Mn levels reflecting the prenatal, postnatal, and childhood exposure windows were not correlated with adolescent hair Mn levels (r = 0:01, p = 0:86; r = 0:08, p = 0:29; and r = 0:10, p = 0:17, respectively, for the three exposure windows). Although we lacked concurrent samples of teeth and hair, the weak correlations estimated between tooth and hair Mn reduced our concern that Mn associations can be explained entirely by prenatal and/or postnatal metals levels.
Although the use of hair as a biomarker for metal exposure has been scrutinized owing to the potential for exogenous contamination (Skröder et al. 2017), we used a conservative, validated method to extensively wash hair of exogenous metals (Eastman et al. 2013;Jursa et al. 2018). This may partly explain why hair Mn concentrations measured in our study (mean = 128 ng=g; range = 11-1,130 ng=g were lower than the range of concentrations reported in other Mn studies of communities near ferromanganese industry, including among Ohio children (mean = 360 ng=g, range 17-15,970) (Haynes et al. 2018) and Brazilian children (mean = 5,830 ng=g, range = 100-86,680) (Menezes-Filho et al. 2011), and U.S. children living near a Superfund site in Oklahoma (mean 471 ng=g, range 89-2,145) (Wright et al. 2006). We expect hair metals levels in our data to represent primarily endogenous levels. Hair Cu and Cr levels in our study (Cu: median = 9:6 lg=g, range: 1.7-191; Cr: median = 0:05 lg=g, range = 0:006-1:81 lg=g) were in the range of global reference values in hair for children 3-15 years of age (range: 7:2-82:7 lg=g; 0:001-4:56 lg=g) (Mikulewicz et al. 2013).
There are numerous strengths to our study. Our use of BKMR allowed us to account for nonlinearity and the potential for highdimensional interactions between metals in association with adolescent cognition. Further, our capabilities to explore the interactive, nonlinear, and joint associations of the metal mixture on IQ extended beyond information gleaned from multiple linear regression. The use of this method is particularly relevant for metals given that a) many metals are essential but neurotoxic in excess and in turn may have nonlinear dose-response curves; b) metals commonly co-occur in the environment (i.e., have the potential to be correlated); and c) the presence of some metals can influence another's metabolism, distribution, and, possibly, mechanisms of toxicity. Although BKMR requires a larger sample size to sufficiently explore high-dimensional interactions, we employed this method using a large data set with a wealth of biomarkers. This was particularly imperative given the lack of consensus on Mn, Cu, and Cr biomarkers of environmental exposure. Use of the hierarchical variable selection procedure in secondary analyses allowed us to identify which biomarkers for Mn, Cu, and Cr were most associated with VIQ and informs biomarker prioritization for future research. This study evaluated the health effects of multiple biomarkers simultaneously in the context of a mixture and to utilize a data-driven approach to inform the selection of biomarkers. Our findings suggest that saliva Mn may be highly informative and particularly important to collect in settings where inhalation exposures are relevant.

Conclusions
Our study investigated the association between a metal mixture of Mn, Pb, Cu, and Cr and cognition in adolescence, a potentially sensitive window for neurodevelopment. Our findings suggest that joint exposure to these metals during adolescence has negative associations with concurrent cognition and that saliva Mn is strongly related to decreased VIQ. We also observed an inverted U-shaped association between ln-Cu and IQ scores, where both low and high concentrations of Cu measured in hair were associated with lower IQ scores. Future research evaluating the joint, interactive, and individual associations of Mn, Pb, and Cu on neurodevelopmental toxicity would further inform this field. Our findings are relevant for improving the health of communities, particularly those where inhalation of metals is an important exposure pathway due to historic or current industrial activity.