In Vitro Effects of Lead on Gene Expression in Neural Stem Cells and Associations between Up-regulated Genes and Cognitive Scores in Children

NSCs derived from NIH-approved H9 (WA09) human embryonic stem cells were purchased from Life Technologies and cultured according to the supplier’s protocol. An aqueous solution of 1 mM Pb acetate trihydrate (cat. no. 316512; Sigma Aldrich) stock was used in all experiments. Transfection of siRNAs was performed with DharmaFECT 1 (ThermoFisher) following the manufacturer’s protocol. All siRNAs were obtained from Sigma—nontargeting control (SIC001), si-NRF2-1 (SASI_Hs01_00182393), si-NRF2-2 (SASI_Hs02_00341015), and si-KEAP1 (SASI_Hs01_00080908). All experiments were performed in passage 3 cells.

For the MTT [3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyl tetrazolium bromide] assay, cells were seeded 24 hr prior to exposure at 1 × 10⁴ per well of a 96 well plate. Exposure to 0, 0.5, 1, 2, 5, and 10 μM Pb was performed in 8 replicate wells. The assay was performed according to the MTT manufacturer’s protocol (Sigma Aldrich). Briefly, 0.05 mg of MTT was added to each well for 3 hr. Formazan crystals were solubilized in 10% Triton™ X-100 plus and 0.1 N HCl in anhydrous isopropanol after repeated pipetting. Absorbance was read at 570 nm, and background at 690 nm removed. Mean absorbance, which is correlated with cell number, is reported along with the standard error of the mean of eight replicate samples. For growth assay, human neural stem cells (hNSCs) were seeded in 24 well plates at 5 × 10⁴ per well and treated with the control vehicle phosphate-buffered saline (PBS) or human recombinant SPP1 protein (Eton Bioscience) at 50 or 250 ng/mL. The next day, hNSCs were exposed to 2 μM Pb for 3 days. Cell counting was done by hemocytometer with Trypan blue staining to exclude dead cells. Six replicates were done for each condition.

Poly-adenylated RNA species were isolated from 1 μg of total RNA and converted to a cDNA library for RNA sequencing using the TruSeq® RNA v2 kit (Illumina). Sample preparation involves isolating poly-adenylated RNA, RNA fragmentation, cDNA synthesis, ligation of adapters, PCR amplification using DNA barcodes, and library validation and quantification. Four samples were multiplexed into a single lane of the Illumina HiSeq 2000 for paired-end reads of 100 bp. Sequencing was performed at the Bauer Core Illumina Sequencing Facility (FAS Center for Systems Biology, Cambridge, MA).

Low-quality reads (< 25 phred), adaptors, and poly-A tails were trimmed with cutadapt (Martin 2011). Read pairs with one or more reads shorter than 20 bp were removed. Quality of reads was assessed using FASTQC (Babraham Bioinformatics). Reads were aligned to human genome build 19 using Tophat2 (Trapnell et al. 2009) and compiled into count tables using HTseq-count (International HapMap Consortium 2003). Counts were normalized in edgeR (Robinson and Oshlack 2010). Differential expression was determined by a generalized linear model. Differentially regulated transcripts were identified following a Benjamini–Hochberg multiple testing correction (q < 0.05) that had a greater than ± 0.2 fold change and a minimum count per million mapped of one.

RNA was reverse transcribed using SuperScript™ III reverse transcriptase and oligo-dT (Life Technologies). The resulting cDNA was amplified using 2× SYBR® mix (Qiagen) and 3 mM of each primer in a StepOne Plus Thermocycler (Applied Biosystems) in quantitative reverse transcriptase polymerase chain reaction (qRT-PCR). Melt curves were checked for single-length amplification products. Fold changes were calculated using the 2-ΔΔCt method. GAPDH is the housekeeping gene used for normalization in all qPCR assays. All primers used in this study and their respective sources or design are listed in Table S1.

SPP1 levels were assessed in whole cell extract using standard Western blotting procedures with 1:1,000 Anti-Osteopontin antibody (EPR3688; abcam). Relative protein concentrations were quantified in Image-J (NIH; see https://imagej.nih.gov/ij/). SPP1 levels in cell culture media was assessed using the human osteopontin (OPN) Quantikine® ELISA Kit (DOST00; R&D Systems). Media was sampled 60 hr post-transfection, after 48 hr of contact with the cells.

NSCs were expanded to approximately 8 million cells, of which half were enriched for NRF2 using siRNA knockdown of NRF2’s negative regulator KEAP1 and the other half were transfected with a nontargeting siRNA control. Samples were prepared following the SimpleChIP® Plus Enzymatic Chromatin IP Kit protocol supplied by the manufacturer (Cell Signaling Technology). Briefly, 48 hr post-transfection, NRF2 was cross-linked to DNA using 1.5% formaldehyde. Nuclei were collected and lysed by sonication. Chromatin DNA was digested with micrococcal nuclease for 18 min into small fragments (150–900 bp). Nuclear extracts were incubated overnight with NRF2 antibody (cat. no. 12721; Cell Signaling Technology) and antibody-bound complexes were captured by SureBeads™ protein G magnetic beads (Bio-Rad). Bound DNA was purified and underwent quantitation by PCR using primers for putative SPP1 ARE, NQO1 ARE (Chorley et al. 2012) and RPL30-exon 3 (Cell Signaling Technology).

The cohort of infants analyzed in this study is a subset of the Early Life Exposure in Mexico to ENvironmental Toxicants (ELEMENT) prospective birth cohort, which was designed to assess the roles of environmental and social stressors in birth outcomes as well as infant and child development. The characteristics of the cohort are provided in Table S2. Between 2007 and 2011, mothers were recruited during pregnancy and only one child for each mother was included in the study. Relevant to this study, prenatal Pb exposure was assessed from maternal blood during the second trimester of pregnancy. Infant neurodevelopment was assessed at 24 months of age using a Spanish version of the Bayley Scales of Infant and Toddler Development, Third Edition (Bayley 2005). Three primary outcome indices are derived from the assessment: the Cognitive Development Index (CDI), the Language Development Index (LDI) and the Psychomotor Development Index (PDI). Detailed information on the study design and data collection procedures have been published previously (Ettinger et al. 2009; González-Cossío et al. 1997; Hernandez-Avila et al. 2002). The human subjects committees of the National Institutes of Public Health in Mexico, Harvard T.H. Chan School of Public Health, Icahn School of Medicine at Mt. Sinai, and participating hospitals approved all study materials and procedures. Women and children who were old enough signed informed consent letters before enrollment. Consent obtained at enrollment applies to the research described in this study.

In the second trimester of each expectant mother, maternal venous blood was collected in trace element-free tubes and frozen. Samples were shipped at 4°C to the Trace Metals Laboratory at the Harvard T.H. Chan School of Public Health, Boston, MA. Samples were processed in a dedicated trace metal clean room outfitted with a Class 100 clean hood using glassware cleaned for 24 hr in 10% HNO3 and rinsed several times with 18Ω Milli-Q water. Approximately 1 g of blood from each mother was digested in 2 mL concentrated nitric acid for 24 hr, and subsequently overnight in 30% hydrogen peroxide (1 mL per 1 g of blood). Samples were diluted to 10 mL with deionized water. Acid-digested samples were analyzed for total Pb using dynamic reaction cell inductively coupled plasma mass spectrometry (DRC-ICP-MS, Perkin Elmer). Final values are the average of five replicate measurements for each individual sample.

Infant neurodevelopment was assessed at 24 months of age using a Spanish version of the Bayley Scales of Infant and Toddler Development (Bayley 2005). Three primary outcome indices were derived from the assessment: the CDI, a composite variable of test scores pertaining to cognition; the LDI, a composite variable of test scores pertaining to expressive and receptive language; and the PDI score, a composite variable of test scores pertaining to fine and gross motor skills.

DNA was extracted from umbilical cord blood samples in the ELEMENT cohort using the Puregene DNA isolation kit (Gentra Systems) and stored at –20°C. Samples were genotyped using the high density Illumina 1 Million Duo chip at the Center for Applied Genomics of the Children’s Hospital of Philadelphia.

Following the quality control assessment, genotypes of 16 SNPs within the SPP1 transcribed locus or within the 10 kb flanking regions were available for 462 infants in the ELEMENT cohort. Linear regression analyses for both main effect and interaction were performed using PLINK (version 1.07) (Purcell et al. 2007). For main effect, linear regression analyses of minor allele copy number on the three outcomes (CDI, LDI, and PDI) were adjusted for sex, gestational age, maternal age, marital status, presence of siblings, maternal education (high school vs. no high school), age at Bayley Scale assessment, and genome-wide principal components 1 and 2. For interaction analyses, natural log transformation of second trimester maternal blood Pb and a multiplicative interaction term of minor allele copy number and natural log transformation of second trimester maternal blood Pb were included in regression analyses. The p-value cut off for statistical significance was determined using the method proposed by Li and Ji (2005), which takes into account that each SNP test is not an independent test given the linkage disequilibrium (LD) among neighboring SNPs. An LD map and haplotypes were generated using the genotyped data using LD-Plus (Vanderbilt University).

To better understand the effects of Pb on NSCs, we performed global transcriptional profiling in human NSCs exposed to Pb. We chose to use RNA sequencing (RNA-seq) as we reasoned that the sensitivity of the method might allow for the identification of subtle yet significant changes in gene expression. Because we were particularly interested in the effect of Pb on prenatal neurodevelopment, we used human NSCs that were initially generated from an embryonic stem cell line (H09 line). We exposed cultured human NSCs to Pb or vehicle control (Figure 1A). The Pb concentration used in this study (1 μM, or 20.7 μg/dL) is about four times the current CDC level of concern for blood Pb and is within the range in exposed human populations (Pirkle et al. 1994; Zheng et al. 2008). Exposure of NSCs cells to 1 μM Pb for 24 hr resulted in a slight decrease (5%) in cell number compared to that in the control vehicle-treated cells (Figure 1B). This is consistent with previous studies showing the inhibitory effect of Pb on NSC proliferation (Breier et al. 2008; Huang and Schneider 2004).

Using total RNAs from control and Pb-treated NSCs, we constructed RNA-seq libraries, each with a unique barcode that allows multiplexing. To minimize the variation of sequencing runs, we pooled barcoded RNA-seq libraries for next generation deep sequencing. We obtained an average of ~ 38 million reads per sample and tested for differential expression of GRCh37 Ensembl-annotated genes. Following a stringent Benjamini–Hochberg multiple testing correction (α < 0.05), we identified a total of 19 differentially expressed genes (3 down-regulated and 16 up-regulated) in Pb-treated NSCs, as shown in Table 1 and in the volcano plot in Figure 1C. Full results from the differential expression analysis are included in Excel File Table S1.

Among the most statistically significant up-regulated genes in Pb-treated NSCs are NQO1 and HMOX1, which are well known targets involved in the cellular response to oxidative stress. The cellular oxidative stress response is mediated by the master transcriptional factor NRF2 (Kensler et al. 2007). NRF2 works to activate transcription by binding to the antioxidant responsive elements (AREs) in target genes. Many genes in addition to NQO1 and HMOX1 also contain ARE elements and are targets of NRF2. We thus examined the rest of the Pb-up-regulated genes and found that at least 10 out of 16 genes had been previously identified as direct targets of NRF2 (Table 1). The expression changes of many of these genes induced by Pb were confirmed by qPCR. As shown in Figure 1D, there is a remarkable consistency between the levels of gene expression measured by RNA-seq and qPCR. The induction of NQO1, HMOX1 and the other known NRF2 target genes strongly suggests that Pb elicits oxidative stress and activates NRF2 in NSCs. Using NQO1 expression as a surrogate marker for NRF2 activation, we determined the dose response of NSCs to Pb. As shown in Figure 1E, Pb as low as 0.1 μM induced a significant increase (33%) in NQO1 expression, indicating that NRF2 activation may be particularly sensitive to Pb exposure in NSCs.

Four genes (SPP1, F2RL2, EGF, and SLC7A8) that were not previously known as NRF2 targets are up-regulated by Pb in NSCs (Figure 1C and Table 1). SPP1, also known as OPN, is an extracellular matrix protein that has been shown to be up-regulated in neuro-injury and is implicated as a potent neuroprotectant (Meller et al. 2005; Topkoru et al. 2013). Because of its potential connection to neural function, SPP1 was chosen for the follow-up characterization. The up-regulation of SPP1 by Pb in NSC cells was confirmed by qRT-PCR (Figure 2A). The SPP1 gene has three splice variants, all of which were up-regulated upon Pb exposure (Figure 2B), indicating that the mechanism involved in SPP1 up-regulation by Pb is not splicing variant specific. Comparison of Ct values indicates SPP1-A is the dominant form in NSCs, with mRNA levels ~ 10-fold higher than SPP1-B and ~ 20-fold higher than SPP1-C (data not shown). Dose response showed that Pb induced SPP1 mRNA expression at 0.1 μM and that the effect maximized at around 2 μM (Figure 2C). The extent of SPP1 induction by Pb was higher at 20 hr than at later time points (Figure 2D), suggesting a potential negative feedback regulation of SPP1 mRNA expression. Consistent with mRNA up-regulation, Western blotting showed that total SPP1 protein level is increased in Pb-treated NSCs (Figure 2E). Since SPP1 is a secreted protein, we measured the amount of SPP1 protein in the culture medium of Pb-exposed and unexposed NSCs. As shown in Figure 2F, there was more SPP1 protein in the media of NSCs cell culture of Pb-exposed NSCs than in the control cells. After 60 hr of Pb exposure, the level of SPP1 in culture media of Pb-exposed cells was 1.8-fold higher than in unexposed cells.

Studies have shown that SPP1 is pro-proliferative (Kalluri and Dempsey 2012) and mediates the survival and proliferation of neural stem cells (Rabenstein et al. 2015). SPP1 up-regulation by Pb thus may constitute a mechanism to protect NSCs from Pb toxicity. We have shown that Pb inhibits hNSC proliferation (Figure 1B). We thus determined whether SPP1 attenuates the inhibitory effect of Pb on hNSCs proliferation. As shown in Figure 2G, addition of recombinant human SPP1 protein at 250 ng/mL in the culturing medium significantly increased the growth of hNSCs in the presence of Pb (90% vs. 79%). Even at lower concentration (50 ng/mL), SPP1 still increased hNSC cell growth in Pb-treated hNSCs, though to a lesser extent. These data support a neuroprotective role of SPP1 in reducing Pb toxicity in hNSCs.

We next determined whether SPP1 up-regulation by Pb in NSCs is part of the NRF2-mediated oxidative stress response. We exposed NSCs to the canonical NRF2 activator DL-Sulforaphane. As shown in Figure 3A, 1 μM DL-Sulforaphane significantly induced the expression of NQO1 and SPP1 expression. NRF2 is normally sequestered and degraded in the cytoplasm by its negative regulator KEAP1 (Kensler et al. 2007). Upon oxidative stress, NRF2 dissociates from KEAP1, accumulates and then translocates to the nucleus. Thus, NRF2 can be activated by inactivation of KEAP1. A pooled siRNA-mediated knocked down of KEAP1 by > 70% (Figure 3B) led to an increase in both SPP1 mRNA (9.6-fold 48 hr post siRNA transfection) and in the secreted SPP1 protein (3.9-fold, 60 hr post siRNA transfection) (Figure 3C,D). Together, these results indicate that NRF2 activation leads to increased SPP1 expression.

To test whether NRF2 is required for SPP1 up-regulation, we knocked down NRF2 in NSCs and subsequently subjected the NSCs to Pb treatment. As shown in Figure 3E, two siRNAs both efficiently knocked down NRF2. Pb exposure itself did not affect either the baseline NRF2 expression or the knockdown efficiency (Figure 3E). However, knockdown of NRF2 significantly attenuates Pb-induced SPP1 up-regulation (Figure 3F). The effect is likely specific as two different NRF2 siRNAs produced similar attenuation on SPP1 up-regulation. These data indicate that up-regulation of SPP1 by Pb exposure in NSCs is mediated by NRF2.

NRF2 controls target gene expression by binding to specific DNA sequences known as AREs within the promoters of target genes. Analyses of AREs in the promoters of canonical NRF2 target genes have identified a consensus sequence motif (RTKAYnnnGCR) that is required for NRF2 binding (Erickson et al. 2002). Using an ARE position weight matrix (Wang et al. 2007), we examined the promoter of SPP1 and identified a putative ARE sequence ~ 600 bp up-stream of the transcription start site (Figure 4A). To test whether NRF2 directly interacts with the putative SPP1 ARE, we performed a ChIP assay. We used KEAP1 knockdown to increase the NRF2 signal in NSCs. The NRF2 protein was immunoprecipitated and bound DNA fragments were examined to detect the presence of the putative SPP1 ARE sequence. We compared the signal level in control and KEAP1-knockdown cells. As shown in Figure 4B, the qPCR signal specific to the SPP1 ARE in KEAP1-knockdown NSCs was significantly higher than that in control NSCs. A similar increase was also observed for a canonical ARE located up-stream of NQO1, whereas no significant increase was observed for a non-NRF2 target sequence, RPL30-exon 3. In addition, control ChIP using rabbit IgG showed very little pull down of SPP1 ARE and there was no difference between KEAP1-knockdown and control cells (see Figure S1). Together these data indicate that the SPP1 gene contains a functional ARE and is a direct target of NRF2.

Because NSCs play an important role in early brain development, we directly examined the association of SPP1 genetic variants with neurodevelopmental outcomes affected by Pb exposure in children. We took advantage of existing genotyping data from genome-wide association studies (GWAS) in the ELEMENT cohort (Z.W., unpublished data, 2016), which was designed to assess the roles of environmental and social stressors in birth outcomes. A total of 16 common SPP1 SNPs (minor allele frequency > 5%) were genotyped in the ELEMENT cohort (n = 462). The relative genomic location and LD of the SNPs are shown in Figure 5A. We performed regression analyses to examine the main effect of the SNPs on the three cognitive outcomes: CDI, LDI, and Psychomotor Development Index (PDI), as well as the SNP interaction with Pb exposure (Table 2). From these analyses, we identified the SNP rs12641001 with a statistically significant main effect association with CDI (p = 0.005) (Figure 5B). Rs12641001 does not show a statistically significant interaction with second trimester Pb exposure. According to the model each copy of the minor allele, T, increases CDI by 2.6 points. Upon examining the LD map (Figure 5A), we found the T-allele of rs12641001 tags two haplotypes that spans the SPP1 promoter region and the first part of the transcript. One of the tagging SNPs (rs2728127) is suggestively associated with CDI with a raw p-value of 0.03471. Rs12641001 is suggestively associated with LDI (p = 0.078) but has no evident association with PDI (p = 0.4089). Therefore, rs12641001 has the largest and most significant association with CDI in children (Figure 5B).