Independent and combined associations of multiple-heavy-metal exposure with lung function: a population-based study in US children

Previous research has found relationships between some single metals and lung function parameters. However, the role of simultaneous multi-metal exposure is poorly understood. The crucial period throughout childhood, when people are most susceptible to environmental dangers, has also been largely ignored. The study aimed to evaluate the joint and individual associations of 12 selected urinary metals with pediatric lung function measures using multi-pollutant approaches. A total of 1227 children aged 6–17  years from the National Health and Nutrition Examination Survey database of the 2007–2012 cycles were used. The metal exposure indicators were 12 urine metals adjusted for urine creatinine, including arsenic (As), barium (Ba), cadmium (Cd), cesium (Cs), cobalt (Co), mercury (Hg), molybdenum (Mo), lead (Pb), antimony (Sb), thallium (Tl), tungsten (Tu), and uranium (Ur). The outcomes of interest were lung function indices, including the 1st second of a forceful exhalation (FEV1), forced vital capacity (FVC), forced expiratory flow between 25 and 7% of vital capacity (FEF25–75%), and peak expiratory flow (PEF). Multivariate linear regression, quantile g-computation (QG-C), and Bayesian kernel machine regression models (BKMR) Yiting Chen and Anda Zhao should be regarded as joint First Authors. Supplementary Information The online version contains supplementary material available at https:// doi. org/ 10. 1007/ s1065302301565-0. Y. Chen · R. Li · W. Kang · S. Tong · S. Li (*) School of Public Health, Shanghai Jiao Tong University School of Medicine, 227 South Chongqing Road, Huangpu District, Shanghai, China e-mail: lsh9907@163.com Y. Chen e-mail: Emmachen9@163.com A. Zhao Shanghai Ninth People’s Hospital, Shanghai Jiao Tong University School of Medicine, Shanghai, China J. Wu · Y. Yin · S. Tong Shanghai Children’s Medical Center, Shanghai Jiao Tong University School of Medicine, Shanghai, China S. Tong School of Public Health, Institute of Environment and Population Health, Anhui Medical University, Hefei, China S. Tong School of Public Health and Social Work, Queensland University of Technology, Brisbane, Australia S. Li MOE-Shanghai Key Laboratory of Children’s Environmental Health, Shanghai Jiao Tong University School of Medicine, Shanghai, China J. Chen (*) College of Public Health, Shanghai University of Medicine & Health Sciences, 279 Zhouzhu Highway, Pudong New Area, Shanghai 201318, China e-mail: chenjy22@sumhs.edu.cn 5214 Environ Geochem Health (2023) 45:5213–5230 1 3 Vol:. (1234567890) were adopted. A significantly negative overall effect of metal mixtures on FEV1 (β = − 161.70, 95% CI − 218.12, − 105.27; p < 0.001), FVC (β = − 182.69, 95% CI − 246.33, − 119.06; p < 0.001), FEF25–75% (β = − 178.86 (95% CI − 274.47, − 83.26; p < 0.001), and PEF (β = − 424.17, 95% CI − 556.55, − 291.80; p < 0.001) was observed. Pb had the largest negative contribution to the negative associations, with posterior inclusion probabilities (PIPs) of 1 for FEV1, FVC, and FEF25–75%, and 0.9966 for PEF. And Pb’s relationship with lung function metrics showed to be nonlinear, with an approximate “L” shape. Potential interactions between Pb and Cd in lung function decline were observed. Ba was positively associated with lung function metrics. Metal mixtures were negatively associated with pediatric lung function. Pb might be a crucial element. Our findings highlight the need for prioritizing children’s environmental health to protect them from later respiratory disorders and to guide future research into the toxic mechanisms of metal-mediated lung function injury in the pediatric population.


Introduction
Pulmonary function is an early predictor for the short-or long-term morbidity and mortality of a wide variety of diseases, including not only potential respiratory health status (Bui et al., 2018), but a vast range of non-respiratory diseases, i.e., cardiovascular diseases, metabolic disorders, and premature death (Agustí et al., 2017;Cheng et al., 2021;Choi et al., 2020). Since the crucial period for lung development is more likely to occur during childhood, identifying modifiable risk factors for poor lung function in youngsters is therefore essential for an individual's later quality of life. Environmental factors are qualifiable and have greater preventive potential; timely and effective awareness and interventions during childhood may have substantial benefits for long-term health.
Heavy metals are a group of ubiquitous contaminants that are present in air, water, food, and soil around the globe. They typically enter the human body through ingestion, dermal contact, or inhalation (Sedman et al., 1994). Mounting evidence has suggested metals are harmful to adults' lung function, particularly for occupational workers . Despite evidence that childhood periods are a critical time window for lung development and likely a vulnerable period for metal exposure, there is limited study on the childhood population, with only seven relevant studies to our knowledge (Feng et al., 2022;Leung et al., 2013;Madrigal et al., 2018Madrigal et al., , 2021Pan et al., 2020;Zeng et al., 2017;Zheng et al., 2013). And six out of seven of the previous researches looked at no more than five blood/urine metal constituents simultaneously (Feng et al., 2022;Leung et al., 2013;Madrigal et al., 2018;Pan et al., 2020;Zeng et al., 2017;Zheng et al., 2013). Only in a recently published study among asthmatic children, 17 toenail metals were analyzed, in which different blood metals were found to be negatively associated with varying pulmonary function variables in singlepollutant models (Madrigal et al., 2021).
Overall, the most studied metals were lead (Pb), mercury (Hg), and cadmium (Cd), whereas no consensus has yet been reached (Feng et al., 2022;Leung et al., 2013;Madrigal et al., 2018Madrigal et al., , 2021Pan et al., 2020;Zeng et al., 2017). Specifically, Hg did not appear to be associated with spirometric indicators (Feng et al., 2022;Madrigal et al., 2021;Pan et al., 2020); Cd exposure was found to be significantly associated with declines in the 1st second of a forceful exhalation (FEV 1 ) (Feng et al., 2022;Madrigal et al., 2021); four other studies suggested null associations (Leung et al., 2013;Madrigal et al., 2018;Pan et al., 2020;Zeng et al., 2017); and blood Pb exposure was shown to be associated with FEV 1 and FVC impairment in US children and adolescents (Feng et al., 2022) and Chinese preschoolers (Leung et al., 2013); however, a study based on 6-17-yearold children (n = 1,234) suggested that urinary Pb was only associated with lower forced expiratory flow between 25 and 75% of vital capacity (FEF 25-75% ) (Madrigal et al., 2018); the absence of association between blood Pb and spirometric metrics was indicated (Madrigal et al., 2018(Madrigal et al., , 2021Zeng et al., 2017). The discrepant results may be due to a singular metal's effect being the primary focus, rather than considering metal interactions and confounding effects of other metals. Metal interaction was only noted in two cross-sectional investigations of children, where interactions between blood Pb and Cd on pulmonary function were diverse, suggesting a positive and absent connection (Pan et al., 2020;Zeng et al., 2017). The real-life exposure pattern is more likely to be mixture exposure, and metals with a high frequency of co-exposure can have synergistic, additive, antagonistic, or potentiating effects (Bhagat et al., 2021). Therefore, uncertain interactions, highly correlated metal exposures, and potentially nonlinear or nonadditive exposure-outcome relationships can all be worrisome. Conventional multivariable parametric regression approaches may not be sufficient in these situations.
In this study, by adopting advanced statistical techniques, herein the quantile g-computation model (QG-C) (Keil et al., 2020) and Bayesian kernel machine regression (BKMR) (Bobb et al., 2015), the individual and combined effects of urinary metals on lung health outcomes would be examined based on a nationally representative sample of children and adolescents who participated in the National Health and Nutrition Examination Survey (NHANES) in the USA from survey cycles 2007-2012. In addition, twelve metals were chosen for our analysis based on their prevalence in the environment, toxicity, potential for biomagnification, and accessibility, including arsenic (As), barium (Ba), Cd, cesium (Cs), cobalt (Co), Hg, molybdenum (Mo), Pb, antimony (Sb), thallium (Tl), tungsten (Tu), and uranium (Ur) (Kumar et al., 2019;Madrigal et al., 2018Madrigal et al., , 2021. More specifically, we wanted to see: (1) whether the joint urine mixture of metals is related to adverse pulmonary function consequences; (2) the exposure-response associations between metal exposures and each lung function parameter; and (3) which individual metal exposure has a greater effect when it occurs as a component of mixtures, or how the parts of a mixture interact. Our findings are expected to inspire fresh perspectives on the connection between lung function and exposure to metal mixtures, and highlight the value of minimizing adverse early life environmental exposures that may have an impact on the long-term health of the susceptible population.

Study population
The NHANES was conducted by the Centers for Disease Control and Prevention's National Center for Health Statistics (NCHS) using a complex sampling frame to achieve a sample representative of the US population. The NHANES study protocols were approved by the NCHS Institutional Review Board, and participants provided written informed consent and/or assent (Statistics, 2017). A total of 30,442 participants were evaluated during the survey cycle from 2007 to 2012, and the data were publicly available. There were 6791 participants, aged 6-17 years. And of these, 1522 participants had complete data on urinary metals and qualified lung function parameters. We also excluded participants with missing data on covariates, including family income-to-poverty ratio (PIR), serum cotinine, body mass index (BMI), and wheezing status (n = 314). Finally, 1227 participants aged 6-17 years who had complete data on covariates were selected for inclusion (Fig. 1).

Outcome measurements
Spirometry participants had to be 6 years old or older. Participants with the following conditions will be excluded: chest pain at the time of the examination; a physical condition with forceful expiration; supplemental oxygen use; recent surgery on the eyes, chest, or abdomen; a recent heart attack, stroke, exposure to tuberculosis, or bloody cough; as well as a history of a detached retina, a collapsed lung, or an aneurysm (NIfHSCN., 2011). According to the recommendations of the American Thoracic Society (ATS) (Miller et al., 2005b), spirometry was carried out while standing using a standardized protocol. Those who were unable to stand up straight were permitted to finish the test while they were seated. Participants aged 6-10 years were asked to exhale for a minimum of 3 s because of agedependent discrepancies in elastic recoil, while those aged 11 years and older were asked to exhale for a minimum of 6 s. Participants repeated the test up to eight times, or until they achieved a reasonable and reproducible spirogram. Each participant's top three spirometry readings were recorded and evaluated according to ATS standards (Miller et al., 2005b). We limited our analysis to only including spirometry data from participants with FEV 1 and FVC values rated A or B. Four pulmonary function parameters-FVC, FEV 1 , FEF 25-75% , and peak expiratory flow (PEF)were the focused outcomes of our data analysis.

Heavy metal measurements
Urinary metal measurement, which is noninvasive, sensitive, and rapid, is an alternate option to  (O'Brien et al., 2016). Urinary heavy metals were therefore used to assess the effect of multi-metal exposure on lung function parameters. Urinary heavy metals were tested in one-third of those aged 6 and up. Urine samples were prepared, saved, and delivered to the laboratory for examination. The following 12 elements were detected in urine by inductively coupled plasma-mass spectrometry (ICP-MS), a multi-element analytical method: beryllium (Be), Co, Mo, Cd, Sb, Cs, Ba, Tu, platinum (Pt), Tl, Pb, and Ur. Utilizing inductively coupled plasma dynamic reaction cell-mass spectrometry (ICP-DRC-MS), the levels of Hg and As in the urine were measured. Because 80% of the data were below the limit of detection (LOD), beryllium and platinum were excluded. As a result, 12 different types of heavy metals-As, Ba, Cd, Cs, Co, Hg, Mo, Pb, Sb, Tl, Tu, and Ur-were examined in our analysis. LOD values were divided by the square root of 2 to impute values below the LOD. To account for the subjects' urine dilution, urinary metal concentrations were corrected for urinary creatinine levels.

Covariates
We collected sociodemographic and basic health information, such as age, sex, race, BMI, PIR, and current wheezing in past 12 months, using standardized questionnaires administered in person. PIR is a ratio of family income to poverty threshold, based on the ratio of the family household income to the poverty level established by the US Department of Health and Human Services, with higher PIR values indicating higher household income (in category: < 1, 1-3, > 3). Current wheezing was a binary variable (yes/no) answered to the question "Had wheezing or whistling in chest in the past 12 months?". BMI was calculated using ageand gender-specific growth charts, and it was then classified as underweight (BMI < 5th percentile), normal weight (BMI ≥ 5th to < 85th percentile), overweight (BMI ≥ 85th to < 95th percentile), and obese (BMI ≥ 95th percentile) (Kuczmarski et al., 2002). LOD for serum cotinine was 0.015 ng/mL, and cotinine levels were split into two groups: below LOD (< 0.015 ng/mL), and above LOD (≥ 0.015 ng/mL) (Hu et al., 2021).

Statistical analysis
Descriptive statistics were used to analyze the sample's demographic and clinical characteristics. The frequency and proportion expressions, the mean and standard deviation (SD), were used for categorical and continuous variables, respectively. Urinary metal levels did not present normality as per the Shapiro-Wilk test (p < 0.05) (Table S1), so they were log-transformed to fit a normal distribution for further statistical analysis. We used the Spearman rank correlation coefficient to evaluate the correlation between metals. Three statistical methods, including multivariate linear regression, QG-C, and BKMR, were used to assess the individual and combined effects of multiple-metal exposures on lung function parameters.
First off, linear regression was implemented, and relationships between metal concentrations as a continuous variable and pulmonary function in children were examined. Initial analyses looked at the univariate associations between each urine metal and each endpoint of interest. Then, covariates including age, sex, race, BMI, PIR, and current wheezing in the past 12 months, BMI, serum cotinine levels were added to each model. Finally, covariate-adjusted multiple-metal models were built to evaluate the relationship between each metal and the pulmonary function parameters. The variance inflation factors (VIFs) in the covariate-adjusted multiple-metal models were less than 3 (Table S2), suggesting there is no evidence for significant collinearity.
Second, the combined effects of metals on children's lung function measurements were evaluated, considering all the covariates. Thus, a chemical mixture can have a large overall collective impact even if each component has a minor effect. QG-C combines weighted quantile sum (WQS) regression and g-computation, which can result in a more precise approximation of the mixture's total effect (Keil et al., 2020). It gives estimates for how a one-quantile increase in all exposures in the mixture will have a simultaneous effect on the outcomes of interest. Weights in the QG-C model can be specified as negative or positive, which represent the proportion of a specific exposure's negative or positive partial effect. In our investigation, the quantile was determined to have a metal content rise of one quartile.
Finally, BKMR was used to evaluate the overall relationship between metal mixtures and lung function metrics as well as potential interactions and nonlinear dose responses for each metal. All the analyses were controlled for the covariates mentioned above. The joint effect is calculated by evaluating the expected change in lung function parameters caused by concurrent changes in all the mixture's components from their median levels. When exposure levels to all the other metals were kept at the median, 25th, or 75th percentile, potential interactions between metals and metal-specific exposure-response curves were depicted. To fit the BKMR model, 10,000 iterations of the Markov chain Monte Carlo (MCMC) method were applied. The posterior inclusion probabilities (PIPs) were calculated to determine if metal was significant, and the threshold of 0.5 was used. And here is a description of the BKMR model. Y i is defined as the outcome of interests; z i = (z i1 ,…, z iM ) T is a vector of m exposures, h(.) symbolizes the hypothetical exposure-response function, x i T is a vector comprising a collection of confounders, β represents the effect of the covariates, and the residuals ε i ~ N (0, σ 2 ) are assumed to be independent and identically (iid) normally distributed with a common variance (Bobb et al., 2018).
To achieve representative and unbiased results when analyzing data from a complex survey, sample weights are often recommended. A sample weight was given to each respondent that took non-response into consideration as well as the likelihood of selection into the subsample. If the variables used to compute the sampling weights are further adjusted in linear regression models, it sometimes influences the precision of the results and introduces overadjustment bias (Gelman, 2007). And since complex survey sampling weights cannot be accounted for by either the QG-C or the BKMR models, we presented our main results without incorporating sampling To examine the robustness of our results, a sensitivity analysis was conducted by rerunning the multiple linear regression models after taking the sample weights into account. The subsample weight was used to compute 6-year sample weights in accordance with NCHS recommendations (Johnson et al., 2013). Additionally, children who were currently wheezing were excluded, and the analyses were rerun to verify the sensitivity. Furthermore, because a child's lung function rises linearly with age until the boy's and girl's teenage growth spurts, which occur at about 10 and 12 years, respectively (Hu et al., 2021;Wang et al., 1993). Sensitivity analysis was performed to demonstrate the relationships of metals with lung function parameters in different age subgroups (6-10 and 11-17).
The statistical significance level was set at p < 0.05 for all statistical analyses, which were performed using R software (version 4.4.3). The QG-C and BKMR models were established using the R packages qgcomp and bkmr, respectively. Table 1 shows the demographic features of the 1227 population. The participants' average age was 11.74 ± 3.37, and there were slightly fewer females than males (49.7% vs. 50.3%). The majority of participants were non-Hispanic whites (30.6%), had a poverty level above 1 (68.6%), were without wheezing symptoms (87.8%), had a cotinine level above the LOD (73.6%), and had a normal BMI (59.0%).

Urinary metal concentrations distribution and correlation
The detectable rates of most urinary metals were > 85%, except for Sb and Cd, with a detectable rate of 78.200% and 65.300%, respectively. The geometric mean and distribution of the urine metal exposures are shown in Table S3. The minimum and maximum geometric means of the metal exposures of Ur and Mo were 0.008 μg/g creatinine and 63.077 μg/g creatinine, respectively. The Spearman correlation analysis showed that 12 urinary metals were significantly correlated among themselves (p < 0.01) (Fig. S1), and the correlation coefficients ranged between 0.03 and 0.70.
The association between metal exposure and lung function parameters Table 2 presents the associations between 12 logtransformed urinary metals and lung function parameters using a multivariate linear regression model. All metals showed a significant inverse association with pulmonary functions, but in single-metal models except for Ba, significant relationships were only revealed with FEF 25-75% (β = − 79.978, 95% CI − 158.783, − 1.17, p = 0.047) and PEF (β = − 215.441, 95% CI − 352.871, − 78.011, p = 0.002). After the adjustment of all covariates, most significant cases remained while the associations of Ba and Ur with pulmonary function measurements were not observed. The linear regression models including full metals were created, considering all possible covariates in order to assess the simultaneous effects of the multi-metal on pulmonary function. Interestingly, in multi-metal regression models, significant negative associations were observed between Pb and Tu with four pulmonary function parameters; Cd was only negatively associated with FEV 1 (β = − 52.493, 95% CI − 100.552, − 4.435, p = 0.032) and FVC (β = − 70.495, 95% CI − 124.517, − 16.472, p = 0.011), while Ba was positively correlated with FEV 1 (β = 57.462, 95% CI 20.058, 94.867, p = 0.003) and FVC (β = 67.895, 95% CI 25.847, 109.942, p = 0.002). We found similar relationships between individual metal exposure and lung function metrics when we considered the NHANES' complex survey design and sampling weights (Table S4), excluded those children with current wheezing (Table S5), and examined different age subgroups (Table S6).

QG-C
The QG-C model was also utilized to determine the total mixture effects of 12 metals. Figure S2 (Table 3). Figure 2 shows the weights representing the proportion of each metal's positive or negative partial effect in the QG-C model. Pb took the biggest negative weights in the relationship with all four lung function metrics, followed by Tu and Cd who accounted for the main

BKMR
We examined the nonlinear dose responses, joint effects, and possible interactions of metals on lung function indices using the BKMR model. Figure 3 estimates the joint effect of the metal mixture by evaluating the effects of simultaneously increasing quantiles of all metals on lung function metrics relative to when all metals were fixed at their median. The overall adverse effects of metal mixtures on lung function parameters were observed. After correcting for the confounders, Fig. 4 shows the exposure-response curves for each individual metal when all other metals are set at their median concentrations. Most of the exposures demonstrated a nearly linear association. Pb was found to have significantly inverse exposure-response correlations with all four pulmonary function assessments, with an approximate "L"-shaped association with FEV 1 , FVC, and PEF. Nonlinear associations were also indicated in the association of Cd with the endpoints of interest. And Ba showed positive relationships with pulmonary function, which was in line with the results of the CI, confidence interval; As, arsenic; Ba, barium; Cd, cadmium; Co, cobalt; Cs, cesium; Hg, mercury; Mo, molybdenum; Pb, lead; Sb, antimony; Tl, thallium; Tu, tungsten Ur, uranium; FEV 1 , forced expiratory volume in one second; FVC, forced vital capacity; FEF 25-75% : forced expiratory flow between 25 and 75% of vital capacity; PEF, peak expiratory flow a Model adjusted for NHANES cycles, age, gender, race, PIR, wheeze, cotinine, BMI b Model adjusted for all the covariates and analytic metals   Table S7, where the PIP of Pb was the highest for all the 4 lung function parameters. Figure S3 depicts the independent influence of 12 metals on lung function after the other metals in the combination were set to the 25th, 50th, and 75th percentiles. Similar to the results of multivariable linear regression and the QG-C model, Pb had the largest contribution to the negative relationships, while Ba was the main contributor to the positive associations. Figures S4-7 show the probable pairwise interactions, displaying the exposure-response curve for each single metal at specific percentiles (the 25th, 50th, and 75th percentiles) of a second metal. There was no indication of chemical interactions in the parallel exposure-response relationships for most metals. The negative effects of Cd on lung function appeared to be reduced at higher levels of Pb.

Discussion
According to our knowledge, this is the first study to use advanced multi-pollutant models to assess the individual and joint associations of 12 urinary metals with four major pulmonary function parameters (FEV 1 , FVC, FEF 25-75% , and PEF) in the pediatric population (6-17 years old). Overall, we observed a trend of declining pulmonary function with increasing concentrations of the metal mixture. And Pb had the highest weight for the negative associations and revealed an almost "L"-shaped association with FEV 1 , FVC, and PEF. Additionally, there were suggestive antagonistic interactions between Pb and Cd. And Ba was positively associated with lung function metrics. Broadly, we found that Pb was the most highly weighted of the 12 metals studied, and it was inversely associated with all four lung function measures. Our findings are congruent with those of earlier cross-sectional researches in healthy Chinese youngsters and American children and teens (Feng et al., 2022;Pan et al., 2020). The association between blood Pb, Hg, and Cd levels and pulmonary parameters was investigated in 221 healthy children, and only Pb levels were significantly associated with FEV 1 and FVC impairment (Pan et al., 2020). And blood Pb was also found to be the most damaging to FEV 1 and FVC decline when the combined effects of polycyclic aromatic hydrocarbons (PAHs) and blood metals (Pb, Cd, and Hg) on health pulmonary function indices in 6-19-year-old children from the NHANES survey were analyzed (Feng et al., 2022). And unlike most previous studies that merely focused on large airway function (FEV 1 and FVC), our study examined the individual and joint effects of 12 urinary metals, adding evidence on small airway function metrics (FEF 25-75% and PEF). PEF is a flowlimited and effort-independent metric that describes lung and airway mechanics (Tantucci et al., 2002); and FEF 25-75% represented small airway function and was demonstrated to be more sensitive than FEV 1 % in reflecting airway hyperresponsiveness, lower airway inflammation, and allergy sensitization (Ciprandi & Cirillo, 2019;Qin et al., 2021).
Interestingly, discrepancies were found when comparing our study to a prior one that included 1234 children aged 6-17 and used data from the NHANES 2011-2012 survey cycle. With multivariate linear regression, the negative effects of urine Pb were only identified in FEF 25-75% (Madrigal et al., 2018). Divergent methodologies and sample sizes may be responsible for the inconsistencies. A study with 206 preschoolers revealed no association between Pb and lung function levels (Zeng et al., 2017). And exposure scenarios could also induce divergence; for instance, associations between 17 toenail metal concentrations and pulmonary function parameters among children with asthma (n = 39) in Chicago were assessed, and no association of toenail Pb with pulmonary function was found (Madrigal et al., 2021). Aside from the limited sample size, the reliability of many toenail Fig. 3 Overall risk (95% CI) of chemical exposures on a FEV 1 , b FVC, c FEF 25-75% , and d PEF, when comparing all the urinary heavy metals at different percentiles with their median level. The results were estimated by BKMR models adjusted for all covariates. CI, confidence interval; FEV 1 , forced expiratory volume in one second; FVC, forced vital capacity; FEF 25-75% : forced expiratory flow between 25 and 75% of vital capacity; PEF, peak expiratory flow; BKMR, Bayesian kernel machine regression metals has not yet been fully established (White et al., 2018), and an earlier study suggested no correlation between levels of Pb in toenails and urine (Kuiper et al., 2014). The disparity of biomarkers in urine and toenails could reflect the varying exposure timings and accumulations, as well as lead exposure from various sources. The evidence for this is, however, limited, and further research will be needed to clarify this.
Additionally, we discovered that Cd was negatively associated with FEV 1 and FVC and made a significant contribution to the combination effect, while the null association was seen with FEF 25-75% and PEF when considering the effects of several other metals. Cd is a known lung toxicant that can accumulate in tissues and organs via blood circulation, causing respiratory diseases when inhaled, ingested, or absorbed via the skin. Several rat models also provided supportive evidence that exposure to Cd could alter lung structure and function (Lai & Diamond, 1992;Surolia et al., 2015). Epidemiology researches have frequently indicated that Cd exposure is linked with impaired pulmonary function in occupational workers (Moitra et al., 2013(Moitra et al., , 2015, general adults (Sobel et al., 2022;Yoon et al., 2015), and smokers (Mannino et al., 2004;Rokadia & Agarwal, 2013), but there have been few studies on children. Toenail Cd was found to have a significant negative association with FEV 1 in Chicago school-aged asthmatic children (n = 39), but not with FVC (Madrigal et al., 2021). Interestingly, two similar studies with NHANES data, one in children aged 6-17, demonstrated no association of both blood (n = 1234) and urinary (n = 408) Cd with spirometry metrics without taking other metals' effects into consideration (Madrigal et al., 2018), while another study involving 6-19 years children (n = 1734) suggested an inverse relationship with FEV 1 (β = − 89.8, 95% CI − 155.3, − 24.3, p = 0.007), Fig. 4 Univariate exposure-response functions (shaded area indicated 95%CI) between heavy metals and a FEV 1 , b FVC, c FEF 25-75% , and d PEF, while fixing other urinary heavy metal exposures at their median level. The results were estimated by BKMR models adjusted for all covariates. CI, confidence interval; As, arsenic; Ba, barium; Cd, cadmium; Co, cobalt; Cs, cesium; Hg, mercury; Mo, molybdenum; Pb, lead; Sb, antimony; Tl, thallium; Tu, tungsten Ur, uranium; FEV 1 , forced expiratory volume in one second; FVC, forced vital capacity; FEF 25-75% : forced expiratory flow between 25 and 75% of vital capacity; PEF, peak expiratory flow; BKMR, Bayesian kernel machine regression but not with FVC (β = − 71.2, 95% CI − 145.9, 3.5, p = 0.062) when considering the influence of polycyclic aromatic hydrocarbons (PAHs) and other metals, therein Pb and Hg, in multivariate regressions (Feng et al., 2022). Therefore, the methodology differences may explain some of our discrepant results. Additionally, three cross-sectional studies based on the Chinese population reported no association of Cd with FEV 1 and FVC (Leung et al., 2013;Pan et al., 2020;Zeng et al., 2017). Of those, one was among children with a mean age of 12.5 years without a diagnosis of heavy metal exposure (n = 221) (Pan et al., 2020), and one was among preschoolers living in areas with electronic waste recycling facilities (Zeng et al., 2017). And another study with preschool children (n = 893) further indicated no association of Cd with FEF 25-75% and PEF (Leung et al., 2013). Differences in sample size, age, environment, and statistical models may all account for the discrepant results between previous studies and ours.
It is noteworthy that, in our study, Pb and Cd had antagonistic interactions. In contrast to our findings, a newly published cross-sectional study revealed a slight synergistic effect of Pb and Cd on children's lung function (Pan et al., 2020), and another reported an absence interaction (Zeng et al., 2017). The antagonized toxic effect of Pb and Cd co-exposure revealed in the current study could be due to possible competition between Pb and Cd (Pandya et al., 2010). More antagonistic interactions were found in which the presence of Cd mediated Pb uptake and the development of Cd had a significant impact on Pb in mouse kidneys (Smith et al., 2012). Moreover, the urinary half-lives of the metals studied varied from 1 to 4 days for Pb (Nordberg et al., 2014), and up to around 14 years for Cd (Suwazono et al., 2009). Thus, Cd and levels imply chronic exposure, whereas Pb may indicate acute exposures, i.e., within the last few weeks. Our findings may suggest that acute metal exposure may be a stronger predictor of pulmonary function decline. Still, few other studies have addressed the precise nature of the interacting contributions, and additional studies using different samples or other methods are needed to generate empirical evidence.
Tu was first shown to be strongly linked with decreased lung function parameters among children. Thus, the effect of Tu has historically been assessed in factory workers through welding processes, which involve the potential hazards of inhalation exposures that may lead to acute or chronic respiratory diseases (Meo & Al-Khlaiwi, 2003). To our best knowledge, only one recent study has taken Tu in a general population of American Indian adults; in both models of single metals and metal mixtures, associations for Tu were mostly null, with no consistent pattern of relationship for Tu across the different spirometry endpoints (Sobel et al., 2022). Currently, there is no consensus regarding the association of Tu with lung function parameters, and since this is the first study to be conducted in a general population of children, more epidemiological and toxicological evidence is warranted.
Unexpectedly, Ba was found to be positively associated with lung function metrics. The respiratory health impact of Ba was limited and discordant. Similar to our findings, a case-control study in Chinese adults (n = 1102) suggested that urinary barium was likely to have a protective effect on asthma (OR, 0.44, 95% CI 0.27, 0.71) (Huang et al., 2016). The null finding was reported by a study among college students (n = 21) regarding Ba as one of the chemical constituents of particulate matter in relation to pulmonary function . A possible explanation for this difference may be explained by their limited sample size. Two animal models have shown the negative effect of Ba on bronchopulmonary reactivity (Hicks et al., 1986) and inflammatory cell infiltration in the lung (Ueha et al., 2020). Ba does not lead to bioaccumulation in human tissues, but if insoluble compounds are inhaled, their levels can build up in the lungs to cause a benign condition known as "baritosis" (Peana et al., 2021). However, there is most likely a good safety margin between even the highest natural barium intake and the threshold for toxic effects.
In addition, in BKMR models, with the increase in the urinary metal mixture concentration, FEF 25-75% initially decreased and then slowly recovered. One explanation is that, at higher concentrations, the heavy metal mixture may stimulate the body's production of antioxidants, which might shield the lung tissue from the oxidative harm the heavy metals cause (Song et al., 2019;Zhang et al., 2023). Another possibility is that the heavy metal mixture at higher concentrations induces an adaptive response in the lungs, resulting in the activation of cellular pathways involved in repairing damaged lung tissue or the induction of growth factors that promote lung tissue regeneration (Jonckheere et al., 2022;Mitra et al., 2022;Miyata & van Eeden, 2011). However, the observed increase in FEF 25-75% at higher concentrations of heavy metals may not necessarily be beneficial, as it could also be a compensatory mechanism to counteract the negative effects of heavy metal exposure on lung function. Ur and Mo were slightly positively correlated with lung function parameters, which may also be related to the regulatory role of Ur and Mo in the immune system (Caicedo et al., 2010;Miller et al., 2005a;Thompson González et al., 2022). Further research is needed to better understand the mechanisms underlying the observed association between heavy metal exposure and lung function.
Although the specific underlying molecular mechanisms are yet unknown, results from epidemiological and experimental studies suggest that the relationship between multiple metals and lung function is primarily due to the formation of reactive oxygen species (ROS) and oxidative stress brought on by metal exposures. For example, the longitudinal study of 1243 workers with a 4-year follow-up has similar results to the current study; in addition to the effects of Pb on the reduction of FEV 1 , they further suggested that high Pb exposure was associated with an increase in 8-iso-PGF2α, supporting the effect of Pb exposure on elevated lipid peroxidation and redox imbalance, which could eventually lead to morphological and functional changes in lung epithelial cells (Wei et al., 2020). Oxidative stress is the well-documented mechanism mediating the association between Pb and Cd exposure and lung function impairment (Atapaththu et al., 2016;Moitra et al., 2015;Wei et al., 2020). And heavy metals were suggested to increase 8-iso-PGF2α and malonaldehyde (MDA) generation (Ashrap et al., 2021;Jomova & Valko, 2011), and these two biomarkers were commonly used as oxidative stress indicators signaling lipid peroxidation. And sustained exposure to ROS may cause DNA damage in epithelial cells, initiating apoptotic pathways via the overexpression of p53 and TGF, as well as the secretion of numerous factors from fibroblasts (Plataki et al., 2005). Increased apoptosis of pulmonary epithelial cells may result in a loss of cell turnover equilibrium, and fibroblastinduced factors, leading to aberrant reepithelialization (Kuwano, 2007). The above series of mechanistic pathways would result in decreased lung compliance and restricted lung function (Gutteridge & Halliwell, 2000). However, there is still much to learn about the complete picture of putative mechanistic pathways for metal mixtures and pediatric lung function.
Our study's findings need to be interpreted in the context of certain limitations. First, after excluding the population that does not fulfill the research's inclusion requirements, the population's representativeness and the extensibility of the results may be influenced. Second, substituting values below the lower LOD with LOD divided by the square root of 2 is likely to introduce bias (Nie et al., 2010). Third, we made multiple comparisons in both the linear regression models and the BKMR models, which increased the possibility of false positive results. Fourth, because there was only one measurement of the metal concentration per individual at a single point in time, urinary metals could only reflect short-term exposure to metals from all sources, which could raise some potential concerns about exposure misclassification. Multiple samples taken from people at various times can enhance the accuracy of exposure assessment. Fifth, urinary Pb was adopted instead of traditionally used blood Pb levels to indicate lead exposure. Urinary Pb levels and blood Pb levels have previously been proved to be strongly correlated (Li et al., 2018;Sallsten et al., 2022), and urinary Pb may also be used to measure environmental exposure to lead (Higashikawa et al., 2000). Last, although several significant confounders were controlled for, the confounding effects of other factors must not be neglected. Finally, because our study was cross-sectional, it was not possible to draw conclusions about causality. Accordingly, large-scale experimental or prospective design investigations should be conducted to confirm the causal links between metal mixtures and lung function.

Conclusion
Generally, our study observed a negative association between metal mixtures and lung function in children. Among the metal mixtures, Pb contributed primarily to the decrement of all the lung function metrics, and the interaction between Pb and Cd was observed.
The mutual verifications by using the three different models provided public health implications for safeguarding children from exposure to heavy metals as well as for preventing the development of both respiratory and non-respiratory morbidities, which would be of inspiring significance in guiding future research on the toxic mechanisms of metal-mediated lung function injury in the pediatric population.