Associations of Pre- and Postnatal Air Pollution Exposures with Child Behavioral Problems and Cognitive Performance: A U.S. Multi-Cohort Study

As part of the Environmental Influences on Child Health Outcomes (ECHO) study, we aimed to maximize power in the current analysis to estimate effects from air pollution exposures on child behavioral and cognitive outcomes and explore nuanced research questions of effect modification by child sex by pooling all three U.S. pregnancy cohorts that compose the ECHO-PATHWAYS Consortium: the Conditions Affecting Neurocognitive Development and Learning in Early Childhood (CANDLE) study, the Infant Development and Environment Study (TIDES), and the Global Alliance to Prevent Prematurity and Stillbirth (GAPPS). Based in Shelby County, Tennessee, the CANDLE study initially aimed to identify risk factors for child neurodevelopment.⁴⁴ Women were considered eligible if they were 16–40 y old, had medically low-risk singleton pregnancies, and planned to deliver at a participating study hospital. From 2006 to 2011, 1,503 women were recruited in their second trimester from either the general community or affiliated medical group clinics. More than three quarters of the enrolled participants ( $n=\mathrm{1,143}$) completed a clinic visit when the resulting children were age 4–6 y. The TIDES study was designed to examine the impacts of exposure to endocrine disrupting chemicals, notably phthalates, on child health and development.⁴⁵ From 2010 to 2012, recruitment commenced in academic medical centers in four cities: San Francisco, California; Rochester, New York; Minneapolis, Minnesota; and Seattle, Washington. Pregnant women in their first trimester were considered eligible if they were age 18 y or older, were English-speaking, had singleton pregnancies without any serious threat, and planned to deliver at a participating study hospital. There were 749 women who were retained in the study throughout the pregnancy and delivered a live birth, and 551 mother–child dyads completed the 6-y visit. Last, the GAPPS study was established to inform evidence-based treatments and interventions to reduce preterm birth and stillbirth through development of a biorepository (www.gapps.org). Families who had participated in the GAPPS study and met eligibility criteria (consented to contact for future study; had prenatal questionnaire data and biospecimens available; child age was eligible for the 4–6-y visit) were invited to participate in the ECHO-PATHWAYS Consortium. From July 2017 to April 2020, 439 mother–child dyads from Seattle and Yakima, Washington, were enrolled and completed the follow-up visit at child age 4–6 y.

The current study included 1,967 CANDLE, TIDES, and GAPPS children who completed behavioral and cognitive assessments at clinical visits at 4–6 y of age and had valid residential addresses in the pre- and/or postnatal windows reported by parents. Each woman provided informed consent upon enrollment in original cohorts. This analysis uses previously collected data from the three cohorts and was approved by the Human Subjects Division at the University of Washington.

Child behaviors were assessed using the Child Behavior Checklist (CBCL), administered at a visit that occurred at age 4–6 y. The CBCL involves caregiver reporting of a wide range of emotional and behavior problems in children and is broadly used in both research and clinical settings.^46,47 One of two CBCL versions were administered, depending on the child’s age: the CBCL preschool form (ages 1.5–5 y)⁴⁸ or the CBCL school-age form (ages 6–18 y).⁴⁹ All CANDLE participants completed the preschool form, all TIDES participants completed the school-age forms, and GAPPS families completed a mix of both forms. The preschool form includes report of the frequency of 99 child behaviors in the past 2 months, whereas the school-age form includes 112 behaviors in the past 6 months. Caregivers rate these items on a scale of Not True (0), Somewhat or Sometimes True (1), to Very True or Often True (2). Although additional behaviors appropriate for children up to age 18 y are added in the school-age form, given that the children in our study were essentially in the same developmental stage, caregivers reported similar types and counts of behaviors across forms. Therefore, we combined the data collected by two CBCL forms. Because the standardized t-scores estimated from these two CBCL forms differ in whether child sex was adjusted,⁵⁰ we computed the raw score for the Total Problems scale as the primary outcome and further calculated the standardized z-score using all ECHO-PATHWAYS participants with behavioral problems measured by either CBCL form as the reference to verify our findings.

Child cognitive performance was examined at the same visit as the behavior measurement. The IQs of the CANDLE, TIDES, and GAPPS children were assessed using the Stanford-Binet Intelligence Scales, Fifth Edition (SB-5),^51,52 the abbreviated five-subtest version of the Wechsler Intelligence Scale for Children, Fifth Edition (WISC-V),^53,54 and the Wechsler Preschool & Primary Scale of Intelligence, Fourth Edition (WPPSI-IV, age 4:0–7:7 version),^55,56 respectively. The three IQ batteries are examiner-administered, highly reliable, and valid measures of intellectual and cognitive abilities in childhood. All examiners were trained on the administration and scoring by licensed psychologists. They participated in didactic instruction and guided practice, interrater reliability exercises, as well as weekly supervision by psychologists post training. Considering that Full-Scale IQ is less frequently included in large population studies due to its heavy burden on examiners and participants, we aimed to maximize use of these data as others have before us.^57–59 In CANDLE, full-scale IQ was derived from the 10 subtests in the SB-5 addressing five cognitive factors with verbal and nonverbal tests for each factor: knowledge, fluid reasoning, quantitative reasoning, visual-spatial processing, and working memory. In TIDES, Full-Scale IQ was estimated using the Tellegen and Briggs formula,⁶⁰ incorporating five WISC-V domains—verbal comprehension, visual spatial, fluid reasoning, working memory, and processing speed. The calculation of Full-Scale IQ from the WPPSI-IV in GAPPS included scores in the five cognitive domains, similar to the five WISC-V domains in TIDES. Although the specific tests used to capture different domains of cognitive performance may vary across the IQ instruments, full-scale IQ provides a standardized metric of overall performance across all the subtests, with a high correlation across instruments.⁶¹

Detailed residential addresses were collected from CANDLE participants at enrollment and updated at each subsequent point of contact. The availability of address data varied by site in the TIDES study: all participants reported residential addresses at enrollment, at the age 4 y visit, and at the age 6 y visit; participants in Rochester and San Francisco were contacted for one more update between enrollment and the age 4 y visit. GAPPS participants were asked to provide comprehensive address histories at the age 4–6 y visit retrospectively from enrollment to present. Point-based ${\mathrm{NO}}_2$ and ${\mathrm{PM}}_{2.5}$ exposures were estimated from a spatiotemporal model on a 2-wk scale.⁶² This model used monitoring data from regulatory networks, further enhanced with air pollution measurements from intensive research cohort-specific monitors. A geographic information system was used to identify covariates representing land use characteristics that could reflect spatial variability in air pollution distributions, and the dimension-reduced regression covariates were obtained using partial least squares from more than 400 geographic variables. The space–time features of pollution concentrations were decomposed into spatially varying long-term averages, spatially varying seasonal and long-term trends, and spatially correlated but temporally independent residuals, and these components were fitted jointly in a likelihood-based spatiotemporal extension of universal kriging. We estimated biweekly ${\mathrm{NO}}_2$ and ${\mathrm{PM}}_{2.5}$ predictions from region-specific models (three regions for the ${\mathrm{NO}}_2$ models and nine regions for the ${\mathrm{PM}}_{2.5}$ models), and averaged the exposure concentrations over each trimester, the whole pregnancy, and the two postnatal windows from childbirth to 2 y old and from 2 to 4 y old. The 2–4-y ${\mathrm{PM}}_{2.5}$ estimates were missing in 150 participants whose 4-y-old birthday was beyond the prediction window of the spatiotemporal model (30 December 1998 to 4 July 2017). Moving was accounted for by calculating the time-weighted averages of ${\mathrm{NO}}_2$ and ${\mathrm{PM}}_{2.5}$ in the relevant windows based on the reported move date. For families who moved between two points of contact and did not report a move date, we estimated the move date as the midpoint of the two contact dates. We refer readers to a recent paper by Kirwa et al. (2021)⁶³ for a thorough discussion of different air pollution prediction models.

Several indicators for maternal, child, and family characteristics, including multilevel social determinants, were harmonized across cohorts. Maternal characteristics included age at enrollment; region- and inflation-adjusted household annual income⁶⁴; household members (2–3 vs. 4 vs. 5 vs. $\ge$ 6); education level ( $<\text{high school}$ vs. high school/Graduate Equivalency Diploma vs. college/technical school vs. graduate or professional degree); marital status (married/living as married vs. single/living as single); pregnancy smoking (smoker vs. nonsmoker); pregnancy alcohol consumption (ever vs. never); pregnancy vitamin supplement intake (ever vs. never); prepregnancy body mass index (BMI); IQ measured by Wechsler Abbreviated Scale of Intelligence [the composite score of four subtests (Vocabulary, Similarities, Block Design, and Matrix Reasoning) from the first edition⁶⁵ in CANDLE, the composite score of two subtests (Vocabulary and Matrix Reasoning) from the second edition^66,67 in TIDES and GAPPS]; depression, assessed at the visit when child outcome assessments were taken, by either the Center for Epidemiologic Studies Depression Scale (CES-D)⁶⁸ or the Patient-Reported Outcomes Measurement Information System (PROMIS)⁶⁹ and quantified as t-scores; and breastfeeding practice (ever vs. never). Child characteristics included age at cognitive and behavioral assessments; sex (male vs. female); birth order (first born vs. not first born); birth year; and secondhand smoking exposure (anyone living in the child’s home smoked vs. no one smoked). The indices in two of the three domains of the Childhood Opportunity Index (COI) were calculated based on the residential address history in pre- and postnatal windows: A larger index in the social and economic subscale reflected higher neighborhood-level socioeconomic status (SES), and a larger index in the educational subscale indicated better early childhood education quality.⁷⁰ We also included parent-reported child race (White vs. Black, vs. others) as a confounder, given that race is not only a proxy for racial residential segregation, which directly associates with air pollution exposures, but also a strong predictor of socioeconomic position that results in health disparities.⁷¹ American Indian/Alaska Native, Asian, and Native Hawaiian/other Pacific Islanders were grouped together to improve statistical power and enhance the data harmonizability across cohorts.

We conducted descriptive analyses in the entire sample and by cohort to summarize the characteristics of the participants, the air pollution exposures, and the child cognitive and behavioral assessments. Linear regressions with robust standard errors were performed to estimate the associations of individual exposures ( ${\mathrm{PM}}_{2.5}$ or ${\mathrm{NO}}_2$) in each window with child Total Problems score and IQ based on observations with complete data pooled from three cohorts. To enable comparisons with studies with relatively low air pollution levels and variabilities, effect estimates were rescaled to a 2-unit incremental increase, which approximates the interquartile ranges (IQRs) for long-term exposures in the six study sites. Based on substantive knowledge from existing literature, we developed a staged adjustment approach comprising three models, allowing us to explore the influence of increasing levels of adjustment on results. We further created a directed acyclic graph (DAG; Figure S1) to help visualize the relationships. Model 1 (the minimal model) controlled for basic demographics—child sex, child age at outcome assessments, and study site. An indicator of CBCL form was additionally included in the analysis of Total Problems score. We defined Model 2 as the primary model, which was further adjusted for major confounders or precision variables, including child race, maternal education, log-transformed region- and inflation-adjusted household income, household members, an interaction between household members and income to account for nonproportional financial needs of a household grow with additional members,⁷² marital status, maternal age at enrollment, birth order, pregnancy smoking, pregnancy alcohol consumption, maternal depression, maternal IQ, child secondhand smoking exposure, and COI subscales (economics and education). Adjustment for covariates that are only associated with the outcomes (i.e., precision variables) in a linear setting is desirable, because it will improve the precision of the effect estimates by reducing residual variance.⁷³ Model 3, an extended model, included additional adjustments for four covariates that may be proxies for unmeasured confounders⁷⁴: Prepregnancy BMI, prenatal vitamin supplement intake, and breastfeeding may serve as proxies for maternal preventative behaviors, and child year of birth may act as a surrogate for birth cohort effects.

In a secondary analysis, we introduced cross product terms of child sex and individual air pollution exposures in each window in separate primary models (i.e., Model 2) and estimated the interaction p-values as well as sex-specific associations. Several sensitivity analyses were performed to verify the robustness of our main findings. First, we implemented complete case analysis as the primary approach. Although the data for most covariates were nearly complete in the analytic data set, 16.2% of the women did not have IQ measurement, and more than half of the missingness in this covariate came from the GAPPS study. To address the issue of missing data in maternal IQ and other covariates, we began by verifying the assumption that the chance of being a complete case solely depends on the observed covariates, where data were missing at random (MAR).⁷⁵ We constructed the Receiver Operating Characteristic (ROC) curve by computing the metrics of sensitivity and specificity at various threshold settings via 5-fold cross-validation⁷⁶ and further estimated the area under the ROC curve—a measure indicating the probability of models with the observed covariates to correctly distinguish complete cases from noncomplete cases.⁷⁷ The area under the ROC curve indicated that this MAR assumption was plausible, and both complete case analysis and multiple imputation would give unbiased results. We thus employed multiple imputation by chained equations (MICE) in the minimal, primary, and extended models.^75,78 Each missing value was imputed 10 times with 100 iterations between each round of imputation using predictive mean matching.⁷⁹

Second, to investigate the validity of combining data of behavioral outcomes measured by different CBCL forms, we replaced the raw Total Problems score with the standardized $z$-score. Third, considering that children born earlier than 34 wk are at higher risks of substantially lower IQ, attention, and executive function in comparison with those with a gestational week 37 and above,⁸⁰ we restricted the sample to those born at gestational week 34 or later. Fourth, to investigate whether the associations of interest were confounded by the exposures in the other windows, we simultaneously included ${\mathrm{NO}}_2$ or ${\mathrm{PM}}_{2.5}$ estimates in all three trimesters in the primary models and further controlled for postnatal exposures. In addition, we mutually adjusted for exposures over the whole pregnancy and postnatal pollution estimates. Fifth, to further confirm the identified critical exposure windows, we fit fully adjusted constrained distributed lag models (DLM) of biweekly air pollution predictions with varied degrees of freedom (df) from 4 to 9. The prenatal windows were restricted to 38 wk, and those born prior to gestational week 36 were imputed with the measurements in the latest available biweekly window. The prenatal and postnatal biweekly exposures were modeled separately with adjustments for average exposures in the other period(s), considering the different biological mechanisms linking child neurodevelopment with air pollution exposures in utero and after birth. Sixth, to assess the linearity of the relationships between air pollution exposures and outcomes, we conducted generalized additive models with full adjustments in the overall sample, as well as in strata by child sex. Finally, to better understand the impacts of inherent heterogeneity across the cohorts and sites and potential modified confounding by site on the findings, we implemented three additional analyses: a leave-one-cohort-out analysis, a leave-one-site-out analysis, and a set of three models compared to the results of fixed effects models (i.e., the primary analysis), including fixed effects models with site–covariate interactions, mixed-effects models with a fixed intercept and random slopes by site, and mixed-effects models with site–covariate interactions. All analyses were conducted in R (version 3.6.5; R Development Core Team).

The inclusion of participants from enrollment to outcome assessment as well as the analytic sample sizes are illustrated in Figure 1. Among the total analytic sample of 1,967 mother–child dyads, the CANDLE, TIDES, and GAPPS studies contributed 53%, 27%, and 20% of the population, respectively (Table 1). Children were on average 5.2 y of age (SD: 1.0) at outcome measurements, and there was an approximately equal number of boys and girls. Nearly half (48%) were identified as White by their parents and 39% as Black. One-fifth of the children were living with at least one relative who smoked. Mothers had an average age of 28.5 y (SD: 6.0) at enrollment, 56% had a college degree and above, and 70% were married or living as married. The median region- and inflation-adjusted household income was $\$\mathrm{55,800}\text{\hspace{0.17em}USD}$ (IQR: 62,500). Maternal IQ was normally distributed, with a mean of 100 (SD: 17.4). The average t-score of CES-D or PROMIS was 48.5 (SD: 7.4), indicating a typical level of maternal depression on average. Compared across cohorts, the CANDLE cohort comprised a larger proportion of Black participants and low-income families, whereas the TIDES participants were relatively more educated, with a higher household income. The analytic sample was similar to the overall sample of participants at enrollment (Table S1). The distribution of child Total Problems score was slightly right skewed (Figure S2) with a median of 18 (IQR: 22), whereas child IQ was relatively normally distributed with a mean of 102.6 (SD: 15.3). Variations in the distribution of child outcomes by cohort were observed (Table 2).

${\mathrm{NO}}_2$ levels were relatively stable and did not show a clear pattern over time, with a mean ranging from 8.4 (SD: 3.8) to 9.0 (SD: 3.9) ppb in different pre- and postnatal windows (Table 3). The correlation between ${\mathrm{NO}}_2$ in nonoverlapping windows was highest between overall pregnancy and age 0–2 y (Spearman correlation: 0.84) and lowest between the first and third trimester (Spearman correlation: 0.32) (Table S2). Prenatal ${\mathrm{PM}}_{2.5}$ was marginally higher than the exposures in the two postnatal windows. The concentrations were 9.0 (SD: 2.3) ${\mathrm{\mu }\mathrm{g}/\mathrm{m}}^3$, 8.8 (SD: 2.0) ${\mathrm{\mu }\mathrm{g}/\mathrm{m}}^3$, and 8.4 (SD: 1.8) ${\mathrm{\mu }\mathrm{g}/\mathrm{m}}^3$ averaged over the pregnancy, 0–2 y, and 2–4 y, respectively. ${\mathrm{PM}}_{2.5}$ aggregated in shorter windows had a greater variation than the exposures in the longer windows. There was a medium to high correlation of ${\mathrm{PM}}_{2.5}$ across windows (Spearman correlation: 0.50–0.89), mainly driven by between-cohort correlations. The variation in ${\mathrm{PM}}_{2.5}$ was lower than that in ${\mathrm{NO}}_2$, and ${\mathrm{PM}}_{2.5}$ and ${\mathrm{NO}}_2$ in the same period were moderately correlated (Spearman correlation: 0.01–0.47). The ${\mathrm{NO}}_2$ concentrations in Seattle, San Francisco, and Minneapolis were greater than those in Memphis, Rochester, and Yakima (Figure S3). In contrast, Memphis had the highest ${\mathrm{PM}}_{2.5}$ level, followed by San Francisco, Minneapolis, and Rochester. Seattle and Yakima in Washington state had the lowest ${\mathrm{PM}}_{2.5}$ level.

From the interaction models with child sex, we found a stratum-specific association between higher second trimester ${\mathrm{PM}}_{2.5}$ and more behavioral problems in girls [ $\mathrm{\beta }$: 1.50 (95% CI: 0.19, 2.81) in girls, $\mathrm{\beta }$: $-0.35$ (95% CI: $-1.89$, 1.18) in boys, $P_{\text{interaction}}$: 0.026], and a stratum-specific association between higher second trimester ${\mathrm{PM}}_{2.5}$ and a lower IQ in boys [ $\mathrm{\beta }$: $-0.07$ (95% CI: $-0.92$, 0.79) in girls, $\mathrm{\beta }$: $-1.19$ (95% CI: $-2.18$, $-0.2$) in boys, $P_{\text{interaction}}$: 0.040], although the results were null from the primary analysis (Figure 2; Table S3). In both postnatal windows, there was suggestive evidence of stronger estimated effects of ${\mathrm{PM}}_{2.5}$ on Total Problems score in girls than in boys; nevertheless, the insignificant interaction terms indicated that the effect difference between groups might be due to chance [0–2 y: $\mathrm{\beta }$: 3.50 (95% CI: 0.61, 6.39) in girls, $\mathrm{\beta }$: 1.72 (95% CI: $-1.25$, 4.70) in boys, $P_{\text{interaction}}$: 0.101; 2–4 y: $\mathrm{\beta }$: 4.80 (95% CI: 1.25, 8.36) in girls, $\mathrm{\beta }$: 2.74 (95% CI: $-0.74$, 6.23) in boys, $P_{\text{interaction}}$: 0.121]. No sex difference in other associations was detected.

First, to handle the missing data in all covariates in the primary model, we constructed the ROC curve to verify the plausibility of the MAR assumption. An average area under curve of 0.891 (95% CI: 0.868, 0.914) was calculated, suggesting that we had an 89.1% chance of correctly distinguishing a complete case from one with missing data using the model with only covariates (Figure S4). We then implemented MICE to impute the missingness and show the results in Table S4. In comparison with the primary results, the point estimates obtained from MICE for the associations with child Total Problems score were mostly attenuated with a higher precision, whereas there was no clear pattern for the changes in associations with child IQ. The positive associations between ${\mathrm{NO}}_2$ in the first trimester, the second trimester, and the whole pregnancy and Total Problems score achieved statistical significance in both the complete case analysis and MICE [first trimester: $\mathrm{\beta }$: 0.51 (95% CI: 0.03, 0.98); second trimester: $\mathrm{\beta }$: 0.61 (95% CI: 0.11, 1.11); whole pregnancy: $\mathrm{\beta }$: 0.81 (95% CI: 0.17, 1.45)]. The positive association between first trimester ${\mathrm{PM}}_{2.5}$ and IQ also gained a higher precision in MICE [ $\mathrm{\beta }$: 0.70 (95% CI: 0.05, 1.35)]. Second, we estimated a correlation of 0.99 between the raw Total Problems score and the standardized z-score, and the conclusions remained the same after we replaced the standardized z-score as the outcome (Table S5). Third, when restricting the analytic sample to participants born at gestational week 34 or later, the estimated direct effects of air pollution exposures on both outcomes were very similar to the findings from the primary analysis except for a slightly larger association with significance between first trimester ${\mathrm{PM}}_{2.5}$ and IQ [ $\mathrm{\beta }$: 0.88 (95% CI: 0.06, 1.70)] (Table S6). Fourth, after simultaneously adjusting for exposures across windows, we derived significant associations between second trimester ${\mathrm{NO}}_2$ and Total Problems score with smaller effect sizes than those from single exposure models [adjustments for ${\mathrm{NO}}_2$ in other trimesters: $\mathrm{\beta }$: 0.77 (95% CI: 0.01, 1.53); adjustments for ${\mathrm{NO}}_2$ in other trimesters and postnatal windows: $\mathrm{\beta }$: 0.79 (95% CI: 0.04, 1.54)] (Table S7). Stronger positive associations between first trimester ${\mathrm{PM}}_{2.5}$ and IQ were also detected when ${\mathrm{PM}}_{2.5}$ in other windows were included [adjustments for ${\mathrm{PM}}_{2.5}$ in other trimesters: $\mathrm{\beta }$: 0.84 (95% CI: 0.02, 1.66); adjustments for ${\mathrm{PM}}_{2.5}$ in other trimesters and postnatal windows: $\mathrm{\beta }$: 1.29 (95% CI: 0.30, 2.29)]. Other significant associations in the primary analysis were attenuated to null. Constrained DLM results depended on the df and were generally not consistent with our primary findings. The DLM with a df of 8 identified a critical window of gestational week 4–5 where higher ${\mathrm{PM}}_{2.5}$ exposures were related to more behavioral problems (Figure S5), which agreed with the significantly positive association found in the first trimester. DLMs using a df of 4 to 9 also indicated positive associations between ${\mathrm{PM}}_{2.5}$ and IQ in slightly different windows within the range of gestation week 6–13, and a similar relationship has been evident in the analysis with the exposure in the first trimester, although the results are insignificant. Apart from these findings, DLMs either showed associations that contradicted our hypotheses and the primary results or displayed associations that were null in our primary analysis, particularly in DLMs with a high df. In addition, the smooth effect curves generated from GAMs in the overall sample (Figure 3) and in each sex stratum (Figure S6) generally did not indicate significant departures from the conclusions of the primary and secondary analyses, although the association between first trimester ${\mathrm{PM}}_{2.5}$ and Total Problems score appeared to be predominantly driven by the exposure outliers in the high end. Last, leaving out the CANDLE cohort (Memphis) presented the biggest impact on results of the leave-one-cohort-out (Table S8) and leave-one-site-out (Table S9) analyses: All the significant associations detected in the primary analysis became null. However, the exclusion of most other study sites did not cause meaningful changes in the result. Comparing the fully adjusted fixed effects models (i.e., the primary analysis) to fixed effects models with site–covariate interactions as well as mixed-effects models both with and without site–covariate interactions, results were generally consistent. We found attenuation in the association of first trimester ${\mathrm{PM}}_{2.5}$ with Total Problems score and the associations of ${\mathrm{PM}}_{2.5}$ in age 2–4 y with both outcomes in the mixed-effects models, but inclusion of site–covariate interactions corrected these attenuations, raising the possibility of site-specific confounding that is only apparent when one of the sites with a smaller population is upweighted in the mixed-effects model (Table S10).