Human development is a trajectory reflecting dynamic processes mediated by constant interplay between genetic background and environmental exposures. Toxic chemicals can alter this trajectory to produce maladaptive cognitive and behavioral phenotypes. Although only a small subset of children exposed to neurotoxicants present a clinically defined developmental disorder, [1], many more exposed individuals probably exhibit subclinical effects, such as reduced IQ. Variability in outcome with the same levels of observed exposures may be explained, in part, by genetic variability that can alter susceptibility to a neurotoxicant. While cognitive development has genetic contributions [2], multiple genome-wide association studies (GWAS) have failed to identify common variants that alone explain the variance in cognition [3–5]. Interactions between genetic and environmental factors may account for a large proportion of the variance and may provide useful information on toxicant pathogenesis [6]. In a cohort of Mexican mother-infant pairs, apolipoprotein E genotypes are associated with 24-month Mental Development Index (MDI) of the Bayley Scale, which might modify lead neurotoxicity [7].

The weight of evidence supports an association between early life exposure to lead (Pb) and impaired cognitive function in children [8–10]. The central nervous system (CNS) is highly vulnerable to chemical injury during development because critical processes—including neuronal growth, synaptic network formation, neuronal migration, and receptor development—peak during this time [11, 12]. Previous studies have shown that Pb crosses the placenta and accumulates in fetal tissues [13, 14]. While the health effects of early Pb exposure are well described, it remains unclear whether certain individuals are genetically predisposed to adverse outcomes following Pb exposure, as few studies have examined gene–environment (GxE) interactive effects on neurodevelopment. This knowledge is central to understanding the mechanisms of action and to creating effective interventions for prevention and treatment of toxic exposures.

The objective of our study was to examine GxE interactions on neurodevelopmental outcomes in two populations exposed to Pb. We examined prenatal exposure because many genes are expressed only during that specific developmental stage, and also because prenatal exposure has consistently been associated with adverse neurodevelopmental outcomes [15, 16]. Thus, genetic variants that produce GxE interactions may only do so when the exposure timing corresponds to a critical developmental window during which the relevant gene is highly expressed. We conducted a genome-wide GxE interaction study (GWIS) in two birth cohorts of Mexican and Bangladeshi children to identify genetic loci that modify the effect of prenatal Pb exposure on neurodevelopmental scores at approximately 2 years of age. We also performed integrative analyses using in vitro transcriptomic data of Pb-induced changes in human neural stem cells (hNSCs).

We present the first multi-ethnic genome-wide study to discover gene–Pb exposure interactions associated with neurodevelopmental outcomes. Although there are distinct ethnic and cultural backgrounds, socioeconomic conditions, and Pb exposure profiles in the two cohorts examined, we identified one locus on chromosome 8 reaching genome-wide significance for joint effects of two independent SNPs (close to linkage equilibrium) on mental composite scores. Within each cohort, the two SNPs had similar main effects and GxE interactive effects, and coefficients for Pb were also similar in magnitude. Results of joint effects of the two SNPs were consistent with effects of individual SNPs and had markedly increased statistical significance. SNPs rs10503970 and rs9642758 are separated by about 855 kb, comprising the entire genomic region of UNC5D, from ~227 kb upstream to ~67 kb downstream, respectively. Moreover, both SNPs were in either strong or moderate linkage disequilibrium (LD) with certain SNPs within UNC5D which might have potential impacts via LD-linked SNPs (Additional file 3: Figure S2). However, there were no known eQTL or functions related to these LD-linked SNPs. Further, this locus contains few genes besides UNC5D, with the closest one (KCNU1) located over 900 kb downstream of SNP rs9642758. Thus, UNC5D appears to be the candidate gene of this locus for GxE interaction on mental development. Similarly, in a study of prostate cancer, despite a weak association between prostate cancer risk and a missense change in HSD3B1 genes, the joint effects of two polymorphisms in the HSD3B gene family demonstrate stronger evidence of association [37].

UNC5D encodes a member of the human dependence receptor UNC5 family for netrin. Netrins control guidance of CNS commissural axons and peripheral motor axons via association with either DCC (DCC netrin 1 receptor) or some UNC5 receptors, leading to axon attraction or repulsion, respectively [38]. UNC5D is specifically expressed in layer 4 of primary sensory areas of the developing neocortex, and thus may be involved in axon guidance by mediating axon repulsion of neuronal growth cones during neurodevelopment [32]. Temporal trajectory profiling from the Human Brain Transcriptome project suggests that expression of UNC5D peaks during mid to late pregnancy, gradually decreases, and is then maintained consistently at a relatively low level after birth (Additional file 8: Figure S5). UNC5D also acts as a dependence receptor required for apoptosis induction when not associated with netrin [39].

A chromosomal breakpoint of apparently balanced chromosome rearrangements was identified in intron 5 of UNC5D in a family with neurodevelopmental disorders [40]. The proband was noted with a developmental delay at 18 months of age affecting walking and speech, as well as with comprehension and behavioral difficulties at 9 years of age. A younger sister with the same translocation had language delay as well as low-to-average fine and gross motor skills at 2 years of age.

In our GxE interaction models, the main effects for Pb levels in prenatal blood were associated with lower mental composite scores in Mexico cohort, but not in the Bangladesh cohort. Pb is a well-known neurotoxicant; nonetheless, some studies have not observed significant adverse associations [41]. A prospective study revealed that fetal lead exposure was associated with the Bayley Mental Developmental Index (MDI) at 3 and 6 months, but no statistically significant relationships were observed with Bayley MDI scores at 3 years of age [42]. There may be several plausible explanations for a lack of a main Pb effect in the 2–4 year old Bangladeshi children, including an exposure range that is too narrow, or co-exposures to other neurotoxicants that may overwhelm a possible lead effect. However, the direction of residual and/or unmeasured confounding by co-exposures could be difficult to predict given the complex correlations between neurotoxicants in different areas of Bangladesh.

We found that minor alleles of SNPs rs10503970 and rs9642758 were significantly associated with lower neurodevelopment scores. However, the coefficient of GxE interaction was positive, which is in the opposite direction of SNPs and Pb exposure (Table 4). We expected to observe a synergy between risk factors resulting in even lower neurodevelopment scores, but our findings contradict this expectation. A direct interpretation of these findings is that, in contrast to the adverse effects of Pb exposure in the presence of wildtype genotypes, Pb exposure counteracts the detrimental effects of minor alleles on mental and motor composite scores at these sites. A possible explanation is that minor alleles may directly or indirectly change the spatial–temporal trajectories of UNC5D expression in the brain, resulting in deleterious effects on neurodevelopment. Pb neurotoxicity can arise by targeting multiple sites at multiple levels—molecular, cellular, inter-cellular, and tissue—which may also perturb spatial–temporal trajectories of UNC5D and offset detrimental effects of minor alleles.

In addition, there was a significant interaction between Pb levels and two SNPs on mental development scores in two genetically different populations and with different distributions of genotypes, strongly suggesting that UNC5D is an important candidate within the biological pathways of Pb neurotoxicity. However, our genome-wide study could not provide further evidence, so future functional research should focus on investigating and comparing the effects of these SNPs on biological functions of UNC5D, especially on spatial–temporal trajectories of brain tissue. Additionally, future research should also examine functional changes of UNC5D induced by Pb exposure.

By cross-evaluation of top gene lists between genome-wide GxE analyses and in vitro global expression profiling induced by Pb treatment, we identified one common gene, SLC1A5, that ranked high (15th) in in vitro transcriptomic assays but ranked relatively low (368th) in GxE analyses for mental composite score (intronic SNP rs1644343). SLC1A5, a member of the solute carrier 1 (SLC1) family, encodes a plasma membrane transporter for the neurotransmitter glutamate [43]. Major functions of glutamate transporters in the brain, expressed in both neurons and astrocytes, are removal of released glutamate from the synaptic cleft and initiation of the recycling cascade for restoring released glutamate in synaptic vesicles. Thus, glutamate transporter malfunction has been implicated in various nervous system diseases [44–46].

Pb is a ubiquitous environmental heavy metal pollutant whose toxic effects are principally manifested in the CNS of immature organisms [47]. Even at low doses, Pb is a potent neurodevelopmental toxicant associated with cognitive deficits and related behavioral disturbances [48–52]. The molecular basis of observed changes include selective blockade of glutamatergic synapses with a subsequent decrease in long-term potentiation, which underlies learning and memory processes [53, 54]. A possible mechanism is that all glutamate transporters possess specific cysteine residues with redox-active-SH groups, which regulate glutamate transport via the s-glutathionylation process, and these cysteine residues are vulnerable to biological and chemical oxidants, such as heavy metals. Reduced transporter function leads to neurotoxic levels of extracellular glutamate, resulting in neuronal damage via oxidative and excitotoxic mechanisms, which often operate together [55]. Beyond effects at the protein level, our SLC1A5 results provide evidence that Pb toxicity also has effects at transcriptomic and genomic levels on glutamate transporters.

This study incorporated data from cohorts with different experimental designs and different biomarkers of prenatal Pb exposure. It is therefore not surprising that few common genes (only SLC1A5) could be identified from top-ranked GWAS gene lists and from transcriptomic studies (in vitro RNAseq), because there are subsets of genes from different cohorts that might directly or indirectly interact with each other at the network level under a general biological scheme of Pb neurotoxicity. Using stringent criteria of direct gene–gene interactions between studies with exclusion of interactions within studies, a large 52-gene network of direct gene–gene interactions revealed from integrative analysis supported our hypothesis (Fig. 2). This network had five hub genes (≥5 direct interactions), including two genes from GWIS (CHUK and TWIST1) and three genes from in vitro assays (PRKCD, HSPA5, and XBP1). There were direct interactions among hub genes, further dividing them into two subgroups, XBP1–TWIST1 (21 genes) and PRKCD–CHUK–HSPA5 (31 genes), connected by an interaction between two GWIS-identified genes, PRKCQ and PRKG1.

We found that the network contained three members of the protein kinase C family—PRKCD (PKCδ), PRKCH (PKCη), and PRKCQ (PKCθ)—that are known signal transducers in the CNS, with roles regulating synaptic development and signaling involved in learning and memory [54]. PKC is well established as one of the major molecular targets of Pb toxicity, and the effect of Pb on PKCs depends on Pb concentrations and periods of exposure in different regions of the brain. As a biological and physiological analogue to Ca²⁺, Pb can directly activate Ca²⁺-dependent PKCs (α, β1, β2, and γ) at ~1000 times less than the concentration of required calcium [56]. At higher concentrations estimated to exist in the brain tissue of environmentally exposed humans (μM range), Pb inhibits kinase activity of all PKCs in in vitro assays [57, 58]. Interestingly, all three PKCs identified from the network belonged to a PKC subset called novel PKCs (δ, ε, η, and θ), which are Ca²⁺-independent. Among them, PRKCD was a hub gene with the highest density of connections in the network (15 direct interactions with GWIS-identified genes), and PRKCD expression was significantly increased by Pb treatment (fold change = 1.4; p = 4.0 × 10⁻⁴; FDR = 0.055). PKCδ is actively involved in regulation of cell apoptosis, especially through a mitochondrial dysfunction pathway, in response to a large and diverse array of apoptotic stimuli [59]. Oxidative stress, often induced by metal ion exposure, results in overexpression and activation of PKCδ, and activated PKCδ translocates to mitochondria and alters mitochondrial functions, such as decreased membrane potential and release of cytochrome C [60].

Besides PRKCD, the other two hub genes identified from in vitro assays, HSPA5 and XBP1, were key genes in ER stress pathways, a major cellular target of oxidative stress manifested at the maternal–fetal interface from early pregnancy onwards [61]. Growing evidence indicates that some of the toxic effects of Pb are oxidative stress via both ER stress and mitochondrial dysfunction [62]. In vitro studies reveal that Pb exposure upregulates gene and protein expression of HSPA5 (GRP78), a molecular chaperone that binds transiently to proteins traversing through the ER and facilitates their folding, assembly, and transport [63]. Consistently, we observed increased expression of HSPA5 in in vitro assays (fold change = 1.2; p = 1.62 × 10⁻⁶; FDR = 7.18 × 10⁻⁴). Moreover, Pb directly binds to HSPA5, and Pb-bound HSPA5 stimulates the unfolded protein response of ER stress, which can also be activated by XBP1 [63]. Further, among the remaining 16 network genes identified from in vitro assays, 10 genes were involved in oxidative stress (e.g., CALR, CLU, and HSP90B1, which encode chaperones involved in ER stress) as revealed by literature review and functional annotations at UCSC Genome Bioinformatics (http://genome.ucsc.edu/).

Although evidence suggests that two hub genes identified by GWIS participate in modulation of oxidative stress via the NF-κB pathway [64], both CHUK and TWIST1 play important roles in early neurogenesis by governing cell fate decisions of neural differentiation at the neural crest stage [65, 66]. In contrast to the majority of genes from in vitro assays involved in processes relevant to oxidative stress, many GWIS-identified genes within the network had functions related to neurodevelopment, such as axonogenesis, neurogenesis, and synaptic transmission (Additional file 10: Table S3). Specifically, SLC1A2, also a member of the SLC1 family, encodes another plasma membrane transporter for the neurotransmitter glutamate, and GIRK2 encodes an N-methyl-d-aspartate type of excitatory amino acid receptor. All of these genes are involved in glutamatergic synapses, which are direct targets for Pb’s effects in the brain [54]. Taken together, through a common theme of Pb exposure, enrichment of oxidative stress genes (mainly from in vitro assays) and neurodevelopmental genes (mainly from GWIS) within the large network of integrative analysis suggests that certain genetic polymorphisms within/near genes relevant to neurophysiological functions might modify Pb exposure effects via oxidative stress on neurodevelopmental outcomes.

The major limitation of this study was that Pb exposure was only assessed from blood collected at birth in both cohorts. Several GxE interactions, such as UNC5D and SCL1A5, also play important roles in neurodevelopment during early- or mid-stage pregnancy. On the other hand, neurodevelopment was assessed at 2 years of age, and postnatal Pb exposure, which also is a time critical for neurodevelopment, was not measured [67]. Lack of complete information on timing and duration of Pb exposure might further reduce the power of detecting GxE interactions within two different ethnic groups with relatively small sample sizes. Moreover, study cohorts had significantly different profiles of metal exposure. The Bangladesh population is well known for high-level arsenic exposure in addition to Pb exposure [68]. In a Mexican cohort similar to the present study, we previously found that high manganese co-exposure increased Pb toxicity among young children [49]. Despite the importance of examining joint exposures to toxicants, most previous toxicological studies have focused on a single agent and either do not measure or do not adjust for potential confounders or modifying effects of other chemicals. For the purpose of meta-analysis across cohorts, we also only focused on Pb exposure because it was commonly measured in both cohorts.

In conclusion, GWIS revealed several candidates containing genetic polymorphisms that modified the response to prenatal Pb exposure on neurodevelopmental outcomes. Although their biological functions were generally related to neurodevelopment, these candidates seemed to play different roles at different developmental stages. Follow-up integrative analysis suggested that some of the candidates might have functions via Pb-induced oxidative stress. Our findings highlight the importance of exposure timing that corresponds to a critical developmental window in GxE interaction research. Future research should emphasize comprehensive measurement of exposure at various developmental stages.