Study Population
The DCHS is a South African, population-based birth cohort that enrolled pregnant women from two primary health care clinics in peri-urban communities: Mbekweni and TC Newman. T These clinics serve two demographically distinct populations, specifically a majority Black African ancestry community and a majority mixed ancestry community [28]. The DCHS has followed infants from birth until at least 5 years of age [28]. The current study population is composed of 142 mother-child pairs with measures available for cord blood DNA methylation, genotype data, and Bayley Scores of Infant and Toddler Development, 3rd edition (BSID-III) in at least one of the following domains: cognitive development, general adaptive behavior, language, and motor outcomes. Inclusion was also limited to mother-child pairs with measures available for relevant covariates which included principal components (PCs) from genome-wide genotype data, maternal age, maternal smoking, maternal alcohol use, birth weight, sex, and socioeconomic status (SES)(Table 1). Smoking status was determined by maternal urine cotinine levels collected prenatally, while alcohol use was measured via the Alcohol, Smoking and Substance Involvement Screening Test (ASSIST), a tool introduced by the World Health Organization (WHO) and which has shown good validity in LMIC settings [29]. Socioeconomic status was captured as a validated score comprising four socioeconomic indicators: maternal educational attainment, employment status, household income and assets. [29] .
The DCHS staff obtain written consent from mothers on an annual basis and the study was approved by the Ethics Committee of the Faculty of Health Sciences, University of Cape Town, by Stellenbosch University and the Western Cape Provincial Research committee [28].
DNA methylation measurements
As described previously [14], DNA was measured from cord blood collected at time of delivery [30]. DNA methylation measures were obtained with both the Illumina Infinium HumanMethylation450 BeadChips (n=156) and the MethylationEPIC BeadChips (n=160). Pre-processing and statistics were done using R 3.5.1 and raw iDat files were imported into Rstudio where intensity values were converted into beta values. The 450K and EPIC datasets were merged using the minfi R package [31]. Background subtraction, color correction and normalization were performed using the preprocessFunnorm function [32] . Following sample and probe filtering, 273 samples and 409,033 probes remained for downstream analysis. Of these samples, 142 had genotype data, at least one BSID-III score measured at 2 years of age, and data available for all relevant covariates (Table 1). Batch effects were removed using ComBat from the R package sva [33] . Cord blood cell type composition was predicted using the most recent cord blood reference data set [34] .
Neurodevelopment measurements
The BSID-III is a widely used tool for assessing neurodevelopment in children up to 42 months of age. We included scores from across four distinct domains: cognitive development, language skills, motor function, and adaptive behavior [29]. The BSID-III has been validated in LMIC settings and previous research reports its use in the DCHS specifically [29] . As described previously [29], the DCHS assessed neurodevelopment using BSID-III at two years of age. The DCHS BSID-III assessment was conducted by a trained professional and incorporates direct observation of the child as well as caregiver input. Composite scores for cognitive, motor, language, and general adaptive behavior domains were scaled to have a mean of 100 and standard deviation of 15 as per standard use of the tool.
Assessment of Indoor Air Pollution Exposure
As described previously [35–37], PM10 was measured using a personal air sampling pump (AirChek 52; SKC, Eighty Four, PA, USA), connected to a styrene filter cassette (37 mm cassette blank; SKC) with a gravimetrically pre-weighted filter (PVC filter 37 mmx5 mm with support pad; SKC) left in the home for 24 hours during the 2nd trimester of pregnancy [35,36]. Filters were weighed after sampling and the National Institute for Occupational Safety and Health method 0600 was used to calculate an average PM10 concentration over 24 hours [37,38]. These 24-hour average PM10 measurements were used for our analyses.
Statistical Analysis
We used different high-dimensional mediation analysis approaches to investigate whether differential DNAm mediates the association between prenatal indoor air pollution exposure and neurodevelopmental delay. Previous research using this subsample of the DCHS has found a significant total effect of PM10 on neurodevelopment in the cognitive domain, but not in the BSID-III general adaptive behavior, language, or motor function domains [5]. Therefore, we treated the BSID-III cognitive domain as our primary outcome and the other BSID-III neurodevelopment domains were used as secondary outcomes to evaluate consistency of results across domains.
Mediation analyses rely on the following three assumptions: 1) no exposure (PM10) - mediator (DNAm) confounding, 2) no mediator (DNAm) - outcome (BSID-III Score) confounding and 3) no exposure (PM10) -outcome (BSID-III Score) confounding [39]. To fulfill these assumptions to the best of our knowledge, we constructed three directed acyclic graphs (DAGs) to visualize each of these three paths (Figure 1). Confounders were selected based on existing literature [40,41] and a minimal sufficient adjustment set was identified for each pathway. Exposure-mediator models were adjusted for SES score, genetic ancestry, and maternal smoking and mediator-outcome models were adjusted for maternal alcohol use, maternal age, SES score, child sex, genetic ancestry, and maternal smoking. We adjusted for genetic ancestry by including the first five genotype PCs to account for population stratification [14]. Models were also adjusted for the first three cell type principal components (PCs), which explained >90% of cell type heterogeneity [14,42,43]. Birth weight as a proxy for gestational age was recognized as another possible mediator of the exposure-mediator association and is a possible mechanism through which prenatal indoor air pollution exposure could impact DNAm, therefore we did not control for birth weight in our analyses.
To assess the role of DNAm as a potential mediator of the association between prenatal exposure to PM10 and neurodevelopment at age two years, we employed three well-documented methods for high-dimensional mediation analysis: HIMA (high-dimensional mediation analysis), DACT (divide-aggregate composite-null), and gHMA (gene-based high-dimensional mediation analysis) [25–27]. As these methods are novel, there is no consensus as to which is the gold standard [44]. Therefore, we have incorporated an analysis pipeline based on a combination of these methods with traditional causal mediation analysis in order to assess the robustness of results (Figure S1). It should be noted that all models, used both for high-dimensional mediation analysis and causal mediation analysis, were adjusted for the appropriate confounders as defined above.
HIMA, a high-dimensional mediation analysis method introduced by Zhang et al. (2016) [25], employs a dimension reduction technique followed by minimax concave penalty (MCP) - penalized estimation of mediation effects and joint significance testing for mediation effects in order to identify significant mediators. Dimension reduction is performed using the sure independence screening (SIS) method which is built on a correlation learning framework that filters out features that are weakly correlated with the response variable [45]. The HIMA joint significance test rejects the null hypothesis of no mediation only when both the exposure-mediator () and mediator-outcome effects () are significant [25] .
DACT leverages epigenome-wide multiple testing to estimate the proportions of the composite null hypothesis to improve power [26,44]. A preliminary step for the DACT method is to create two linear models pertaining to each CpG site: the exposure-mediator association () and the mediator-outcome association (). While HIMA uses a screening technique to reduce dimensionality, DACT does not involve a screening step by default. We performed a pre-screen of CpG sites based on their association with the exposure (PM10) and the outcome (BSID-III scores of neurodevelopment); only CpG sites with p<0.05 for both associations were included in the downstream analyses. (Figure S1). Given previous findings indicating a negative association between prenatal indoor air pollution exposure and neurodevelopment in the DCHS cohort [5], we chose to additionally filter our sites by only allowing a negative natural indirect effect (NIE) defined by (acting in the same direction as the association between indoor air pollution exposure and neurodevelopment). As we chose to pre-screen our CpG sites, we used the Efron correction feature of the DACT package to estimate the proportions of the composite null (Figure S2) as opposed to Jin and Cai correction which is recommended if performing epigenome-wide mediation effect testing with DACT [26]. Of note, neither of these methods account for correlation between mediators. As such, we calculated Pearson correlation coefficients for all sets of identified CpG sites to better understand if and how these CpG sites were correlated with one another.
The gHMA method was developed by Fang et al. (2020) [27] and focuses on genes as functional units as opposed to individual CpG sites. gHMA is primarily composed of three components: 1) linear mediation analysis, 2) nonlinear mediation analysis and 3) an omnibus test of mediation effects. Significance testing results for both the linear and nonlinear mediation analysis steps are combined using the gHMA omnibus (gHMA-O) test. As the true relationship between mediators and outcomes are often not well understood in practice, gHMA-O transforms and combines p-values from the linear and nonlinear analyses in order to construct the gHMA-O test statistic, which is used to assess mediation effects at the gene level [27] . CpG sites were annotated by closest gene using the Bioconductor package hiAnnotator (https://bioconductor.org/packages/release/bioc/html/hiAnnotator.html) and the Ensembl gene predictions (ensGene, version of Apr-06-2014; http://hgdownload.soe.ucsc.edu/goldenPath/
hg19/database/ensGene.txt.gz) as previously described elsewhere [46].
P-values for CpG sites which passed the screening step and were tested using DACT and HIMA were corrected for multiple testing using the Benjamini-Hochberg false discovery rate correction (BH FDR) [47]. Due to the fact that gHMA assesses differentially methylated gene regions as opposed to individual CpG sites, gHMA p-values were FDR corrected for the total number of gene regions tested, rather than the number of distinct CpG sites. HIMA and DACT CpG sites that remained significant at a false discovery rate of 0.05 were then validated via traditional causal mediation analysis using the function mediate from the R package mediation to obtain estimates of natural indirect effect (NIE), direct effect (DE), total effect (TE), and proportion mediated (PM) [48].
A recent DCHS study identified three CpG sites (cg26971411 [SPTBN4], cg00490349 [intergenic], cg15660740 [intergenic]) associated with neurodevelopment measured by BSID-III [14]. We also examined these CpG sites in our causal mediation analysis step to investigate whether they mediate the association between prenatal PM10 exposure and neurodevelopment.
Maternal HIV has been previously linked to neurodevelopment; however, the association between maternal HIV and DNAm is not well understood [49,50]. Therefore, we also conducted a sensitivity analysis to determine the effect of including maternal HIV as a potential confounder of the mediator-outcome association.