Transcriptomic, proteomic, and metabolomic analyses identify candidate pathways linking maternal cadmium exposure to altered neurodevelopment and behavior

Cadmium (Cd) is a ubiquitous toxic heavy metal of major public concern. Despite inefficient placental transfer, maternal Cd exposure impairs fetal growth and development. Increasing evidence from animal models and humans suggests maternal Cd exposure negatively impacts neurodevelopment; however, the underlying molecular mechanisms are unclear. To address this, we utilized multiple -omics approaches in a mouse model of maternal Cd exposure to identify pathways altered in the developing brain. Offspring maternally exposed to Cd presented with enlarged brains proportional to body weights at birth and altered behavior at adulthood. RNA-seq in newborn brains identified exposure-associated increases in Hox gene and myelin marker expression and suggested perturbed retinoic acid (RA) signaling. Proteomic analysis showed altered levels of proteins involved in cellular energy pathways, hypoxic response, and RA signaling. Consistent with transcriptomic and proteomic analyses, we identified increased levels of retinoids in maternally-exposed newborn brains. Metabolomic analyses identified metabolites with significantly altered abundance, supportive of changes to cellular energy pathways and hypoxia. Finally, maternal Cd exposure reduced mitochondrial DNA levels in newborn brains. The identification of multiple pathways perturbed in the developing brain provides a basis for future studies determining the mechanistic links between maternal Cd exposure and altered neurodevelopment and behavior.


Results
Maternal Cd exposure increases proportional brain weight at birth. The exposure model and hybrid mating strategy utilized in this study are described in detail in previous work 15 . Briefly, F 0 mothers were administered one of three Cd doses chronically through their drinking water before and during pregnancy: 0 ppm (control), 1 ppm (low, environmentally relevant dose), and 50 ppm (high dose consistent with other rodent studies). For ease of annotation, the generation of the mice and the Cd dose to which they were exposed will be herein referred to as: F (generation) (Cd dose) . For example, F 0 1ppm refers to mothers exposed to 1 ppm Cd and F 1 50ppm refers to offspring maternally exposed to 50 ppm Cd. A hybrid mating scheme between C57BL/6J ('B') and CAST/EiJ ('C') inbred mouse strains was utilized to analyze parent-of-origin effects in a separate analysis (Hudson and Cowley, unpublished). Herein, F 1 mice will be referred to by their maternal × paternal cross, i.e. B × C (B mother × C father) or C × B (C mother × B father). F 1 mice were dissected within 24 h of birth and herein will be referred to as 'newborn' .
The only statistically significant differences observed in raw brain weights at birth following maternal Cd exposure were in B × C F 1 1ppm males (increase) and C × B F 1 50ppm males (decrease); the remaining groups showed no significant differences (Fig. 1A). We previously showed that maternal Cd exposure leads to a significant reduction in body weight at birth 15 . When brain weight was normalized to body weight to account for these differences, a significant increase in the proportional weight of the brain was observed at birth in the F 1 50ppm pups of both crosses and both sexes, as well as the F 1 1ppm C × B males (Fig. 1B). Differences in raw or proportional brain weights were not observed at 6 months, with the exception of a significant increase in proportional brain weight that persisted in C × B F 1 50ppm females ( Supplementary Fig. 1). To identify any significant cell damage due to maternal Cd exposure, hematoxylin and eosin staining was performed on serial coronal sections of whole newborn F 1 0ppm and F 1 50ppm B × C male and female brains. No overt differences in cell damage or death were identified in corresponding regions between individuals of different maternal Cd doses (data not shown).
Maternal Cd exposure programs perturbed anxiety-related behavior during adulthood. Due to the significant increase in proportional brain weight at birth and the previously reported association of early life exposure to Cd with altered behavior, Social Interaction (SI) and Open Field (OF) behavioral tests were performed in aged F 1 mice 2-4 weeks prior to sacrifice (6 months) to test whether altered behavior in adulthood can be programmed by maternal Cd exposure alone. The SI test is used to assess sociability in the presence of a novel animal 27 , while the OF test is used to measure activity, anxiety, and willingness to explore 28 .
Maternal Cd exposure had a more profound and diverse effect on sociability behaviors assessed in the SI test compared to the anxiety-related behaviors assessed in the OF test. F 1 50ppm B × C females and both exposure groups of B × C F 1 males investigated the novel animal more frequently than control mice, and F 1 50ppm B × C females and F 1 1ppm B × C males spent a longer duration of time near the novel animal than control mice (Fig. 2, Supplementary Table 1). Maternally-exposed B × C animals appeared to initially investigate the novel animal more quickly than controls, though this difference was not statistically significant (Supplementary Table 1 www.nature.com/scientificreports/  www.nature.com/scientificreports/ contrast to B × C animals, maternally-exposed C × B animals appeared to investigate the novel animal less frequently and took longer to initially investigate than controls, though these data were not statistically significant (Supplementary Table 1). In the OF testing, C × B F 1 1ppm males showed a significant increase in activity (as indicated by an increased number of visits to both the perimeter and center), willingness to explore the center of the box, and took less time to initially explore the center compared to control C × B males (Supplementary Table 1). This was a similar, yet non-significant trend in other maternally-exposed mice (Supplementary Table 1). B × C F 1 1ppm females visited the perimeter of the box significantly more than control B × C females, indicative of more anxious behavior (Supplementary Table 1). These behavioral data indicate that, despite the consistent effect of maternal Cd exposure on newborn brain weight, the effect of exposure on offspring behavior at adulthood may vary depending on genetic background. This might be in part due to some significant differences in baseline behaviors between the four control groups (Supplementary Table 1).
Due to the robust Cd-associated behavioral changes observed in B × C F 1 females, as well as our previously published work that identified Cd-associated transcriptomic changes in the newborn heart-which would therefore permit us to identify conserved gene expression changes in multiple organs due to maternal Cd exposure, these animals were selected for further molecular analyses. Subsequent data presented here were obtained in these animals unless otherwise specified.
Transcriptomic analysis of newborn brains identifies maternal Cd-associated activation of myelin markers. In order to understand how maternal Cd exposure alters pathways that could contribute to enlarged proportional brain weight at birth and perturbed behavior at adulthood, RNA-seq was performed on whole brains of newborn F 1 0ppm and F 1 50ppm B × C females (n = 4 for each dose with each individual birthed in a different litter). After quality control filtering, over 13,000 genes mapped to both parental genomes (14,417 aligned to B genome, 13,277 aligned to C genome, Supplementary Table 2). From these, 15 differentially expressed genes (DEGs) that were significant after multiple testing correction and common to both parental genome alignments were identified (Table 1, Fig. 3A). Select DEGs were validated with qRT-PCR, including on F 1 1ppm B × C females (n = 8 for each dose, Fig. 3B). Twelve of the fifteen DEGs are associated with oligodendrocytes or myelination (Table 1), a process performed by oligodendrocytes after birth in the mouse 29 . Eight of the fifteen DEGs are part of the myelin transcriptome 30 , and seven of these eight DEGs had increased expression due to maternal Cd exposure (Table 1). Term enrichment analysis for the 15 DEGs showed a significant enrichment for terms related to retinoic acid (RA) signaling (Supplementary Table 2, terms with adj p val < 0.05 highlighted in red), which is required prior to the onset of myelination 31 . Based on these observations, we hypothesized that the process of myelination is occurring earlier in maternally-exposed F 1 mice. However, we were unable to detect myelin in the prefrontal cortex of these animals at birth using histological approaches (Black-Gold II stain, data not shown), consistent with previous work showing that myelin is challenging to visualize in rodents during early postnatal life due to its low abundance and despite the presence of molecular markers of myelin 29,32 . Maternal Cd exposure perturbs Hox gene expression in the newborn brain. RNA-seq identified Hoxb8 as the most significantly differentially expressed gene in the brain due to maternal Cd exposure (Table 1). Hoxb8 demonstrated negligible expression in control brains, consistent with a more posterior expression pattern www.nature.com/scientificreports/ during normal development 44 , yet was significantly increased in the brains of mice maternally-exposed to Cd. Further evaluation of the RNA-seq data suggested that other Hox genes may be similarly affected, although their differences in expression did not achieve statistical significance using this approach (Supplementary Table 2). To further assess the extent of Hox gene misregulation by maternal Cd exposure, the expression of all Hox 1-9 genes contained in the A-D clusters was quantified with qRT-PCR in newborn B × C female brains exposed maternally to 0 ppm, 1 ppm, or 50 ppm (n = 8 for each group). Consistent detectable amplification (see Materials and Methods for criteria) in at least one of the three treatment groups was only found in genes in the Hoxa and Hoxb clusters (Fig. 4). Expression data were normalized to the 1 ppm group in each respective gene due to undetectable amplification of several Hox genes in multiple control samples (see Materials and Methods for detailed explanation). 50 ppm maternal Cd exposure significantly increased the expression of Hoxb2 and more posteriorly expressed Hoxa and Hoxb genes (i.e. Hoxa5 through Hoxa6 and Hoxb5 through Hoxb8) when compared to controls, and 1 ppm maternal Cd exposure significantly increased the expression of Hoxa6 and Hoxb2 when compared to controls (Fig. 4).  . qRT-PCR of Hox genes in newborn B × C female brains exposed maternally to Cd. Individuals with undetectable Ct values were assigned expression values of 0%. Expression data shown relative to mean of 1 ppm (normalized to 100%) for each gene due to several genes having undetectable levels in controls. *p < 0.05, **p < 0.01, ***p < 0.001 (one-way ANOVA with post-hoc Dunnett's test comparing 1 ppm and 50 ppm to 0 ppm). www.nature.com/scientificreports/ Proteomic analysis of newborn brains identifies maternal Cd-associated changes in the abundance of proteins related to cellular energy, hypoxia, histone 1 subunits, myelin, and retinoic acid. To determine whether the observed changes in gene expression correlated with altered abundances at the protein level, a discovery-based mass spectrometry proteomics approach was performed on newborn B × C female whole brains. 46 proteins were identified as significantly different in abundance in the F 1 50ppm brains compared to controls after correcting for multiple testing and filtering hits based on criteria described in the Materials and Methods ( Table 2, raw data in Supplementary Table 3). In general, a dose-dependent increase or decrease in protein abundance was seen as a result of maternal Cd exposure; however, this was only significant for 6/46 proteins in the F 1 1ppm brains when compared to controls. Among the 46 differentially abundant proteins, there was a significant enrichment for terms related to cellular energy and metabolism pathways and hypoxia (Supplementary Table 4), consistent with a maternal Cd exposureinduced Fe deficiency that we reported previously 15 . Additionally, three histone 1 subunits were significantly decreased in F 1 50ppm brains compared to controls ( Table 2). Consistent with findings from RNA-seq, 20 of the 46 differentially abundant proteins in the F 1 50ppm brains are known to be part of the myelin proteome 45 (Supplementary Table 4), and all twenty were significantly increased in abundance ( Table 2). The most significant term in the ESCAPE enrichment analysis predicted an increase in Sox2 (Supplementary Table 4), which is essential for the proliferation and differentiation of oligodendrocytes during postnatal myelination 37 . Consistent with the findings of RNA-seq, there was a significant enrichment for terms related to RA (Supplementary Table 4).
Maternal Cd exposure leads to reduced mitochondrial content in the brain at birth. Previous work by another group found severely damaged mitochondria in electron micrographs of brains from rats whose mothers were exposed to Cd during pregnancy and lactation 24 . Given this and our proteomic analysis showing alterations in proteins involved in cellular energy pathways, we quantified mitochondrial DNA (mtDNA) content in newborn B × C female brains. mtDNA content, as quantified through qPCR, was significantly reduced by approximately 25% in newborn B × C female brains as a result of both 1 ppm and 50 ppm maternal Cd exposure (Fig. 5).
Metabolomic analyses of newborn brains in response to maternal Cd exposure. We previously showed maternal Cd exposure leads to a significant iron deficiency in newborn pups 15 . Considering the iron deficiency and the proteomic data containing a significant enrichment for terms related to hypoxia and cellular energy and metabolism pathways, we hypothesized that the iron deficiency was leading to a hypoxic environment, resulting in an increase in anaerobic respiration in the brain to meet glucose demands. We therefore quantified the metabolite lactate, a major byproduct of anaerobic respiration that increases in hypoxic conditions 46 in a small subset of B × C female newborn brains. There were no statistically significant differences in lactate quantified as mM per mg of dry extract weight between the three treatment groups, potentially due to our limited sample size (Supplementary Fig. 2A). We hypothesized that myelination may be occurring earlier in neurodevelopment due to maternal Cd exposure, which may alter the lipid proportional weight in the brain. To account for differences in composition of extract weights sampled for each individual as a proportion of their total brain weight, lactate data were normalized to proportion of brain weight. These data were more suggestive of a dosedependent increase in lactate, though these data were not statistically significant ( Supplementary Fig. 2B).
We then performed a mass spectrometry discovery-based approach for identifying metabolites that were differentially abundant in F 1 1ppm and F 1 50ppm B × C female newborn brains when compared to F 1 0ppm controls. Five metabolites were identified in the F 1 1ppm B × C females, while 14 metabolites were identified in the F 1 50ppm B × C females (Supplementary Table 5). However, after multiple testing correction, only 3 metabolites remained significantly different in the F 1 50ppm B × C females: 2-hydroxyvaleric acid, 2-hydroxycaproic acid, and xanthine ( Fig. 6  Maternal Cd exposure leads to increased levels of retinoic acid in the brain at birth. Given our observations of altered levels of RA-associated transcripts and proteins in brains of maternal Cd-exposed mice, and the known role of RA in behavior 47 , myelination 31 , and regulation of Hox gene expression 44 , we hypothesized that RA levels in the brain were perturbed at birth as a result of maternal Cd exposure. To test this, we quantified three retinoid species in newborn B × C male and female whole brains: all-trans RA, retinol (RA substrate), and retinyl ester (retinol storage). Maternal Cd exposure led to significantly increased levels of RA in the brain of both F 1 1ppm and F 1 50ppm pups (Fig. 7A). The F 1 1ppm pups tended to have higher brain levels of RA than the F 1 50ppm pups, though this was not statistically significant. Retinol levels were significantly increased as a result of 50 ppm maternal Cd exposure in females or when the data from both sexes were pooled (Fig. 7B). Retinyl ester was significantly increased due to 1 and 50 ppm maternal Cd doses in males and when the data for sexes were pooled, but was only significantly increased in F 1 50ppm female brains (Fig. 7C). Taken together, these data suggest increased biosynthesis and/or reduced degradation of RA.

Discussion
Exposure to Cd during early life is associated with impaired neurodevelopment and behavior 18,19,[21][22][23] , but the underlying molecular mechanisms responsible for these changes are poorly understood. In the present study, we utilized a mouse model of maternal Cd exposure and demonstrated that maternal exposure alone is sufficient to program altered behavior among F 1 mice during adulthood. Behavioral endpoints for both males and females showed considerable variability within each group, which in females may partly be explained by not controlling for estrus cycle stage. Nonetheless, maternally-exposed mice overall demonstrated perturbation of anxietyrelated behavior. To inform on the molecular mechanisms underlying these behavioral defects, we identified www.nature.com/scientificreports/ exposure-associated transcriptomic, proteomic, and metabolomic changes in the brain at birth. Together, these data suggest that maternal Cd exposure may cause premature myelination in the brain and inappropriate regulation of Hox gene expression. These data also suggest that Cd perturbs RA signaling in the developing brain, Table 2. Significantly differentially abundant proteins in the F 1 50ppm B × C female newborn brains compared to controls, ranked by adj p val in 50 ppm group. Abundance ratios for these proteins in the F 1 1ppm female newborn brains included. a = found in myelin proteome, b = involved in RA pathways, c = histone 1 subunit.  www.nature.com/scientificreports/ which we subsequently confirmed with targeted analyses. A summary of perturbed pathways identified in this work and our subsequent hypotheses of mechanisms that may link maternal Cd exposure to altered neurodevelopment and behavior are shown in Fig. 8. This work provides a foundation for further studies to understand the mechanisms of maternal Cd exposure. Our mouse model also provides the opportunity to study the adverse effects of Cd on other organ systems. Here we report only impacts on neurodevelopment and behavior, but have previously shown that maternal Cd exposure can affect cardiovascular development and program hypertension in a subset of adult F 1 animals 15 . Two doses of Cd were used in this study: 1 ppm Cd, which reflects a realistic level of routine Cd exposure considering the average levels of Cd naturally occurring in the Earth's crust, and 50 ppm Cd, as this dose results in circulating blood Cd levels consistent with those found in people living in heavily Cd polluted areas 15 . In some of our experiments, the low dose of Cd had more pronounced effects than the high dose; this supports the emerging frequency of non-monotonic dose-responses in toxicological studies 48 . Rather than starting exposure at conception, our mouse model began maternal Cd exposure prior to conception in order to more closely reflect typical patterns of chronic human exposure. Additionally, our mouse model used ingestion of Cd-laced water, rather than injection, to more closely reflect a common route of Cd exposure in humans.
In rodent brains, myelination is an event that primarily occurs during postnatal development 29 . Our RNA-seq and proteomic analyses are suggestive of an earlier onset of myelination due to maternal Cd exposure. Twelve of the fifteen maternal Cd-associated DEGs are involved in myelination or oligodendrocyte differentiation during myelination, and two classical myelin markers (Mbp and Mag) were significantly increased in newborn pup brains in a dose-dependent manner. In addition, 20 proteins that are part of the myelin proteome were found to be significantly increased in abundance due to 50 ppm maternal Cd exposure. Finally, we previously showed that maternal Cd exposure increases circulating fetal copper levels 15 , an element that is required for myelination 49 . To date, the link between Cd and myelin has only been made following postnatal or adult exposure to Cd and suggests a negative effect on myelination 50,51 . However, one study in rats that induced a delay in early postnatal myelination through lysophosphatidylcholine injection into the ventral hippocampus of 10 day old rat pups noted opposite behavioral effects to those seen in this study (i.e. decreased activity and increased anxietyrelated behaviors) 52 . There is emerging evidence that myelin plasticity, remodeling, and altered developmental myelination affects behavior 35,53,54 , thereby supporting a potential mechanistic link through which maternal Cd exposure impacts offspring behavior. Myelin proteins and other molecular markers of myelination can be found in the rodent brain as early as 7 days postnatally 29 , but myelin is difficult to detect using histological approaches, consistent with our inability to visualize myelin in newborn maternally-exposed pup brains with black-gold staining (data not shown).
We also show that Hox gene expression is inappropriately regulated following maternal Cd exposure. The Hox gene family is conserved across metazoan species and is essential during development to confer anterior to posterior identity to tissues. Hox gene expression patterns are uniquely linked to their position along their respective chromosomes: 3' Hox genes are expressed in the anterior end of the embryo, while 5' Hox genes are expressed more posteriorly. Hox expression is regulated in part by a RA gradient in the developing embryo. Mammals have four clusters of Hox genes (A, B, C, D), each cluster located on a different chromosome, and the genes are numbered 1-13 based on their 3' to 5' location within their cluster 44 . Hox 1-4 paralogs are expressed in the mammalian central nervous system, primarily in the developing hindbrain; the remaining Hox 5-13 paralogs are expressed in the spinal cord and other posterior tissues 55 . In this study, we have shown that expression of the Hoxa and Hoxb clusters is perturbed due to maternal Cd exposure, with regards to 2-4 paralogs that www.nature.com/scientificreports/ are normally expressed in the brain as well as more posteriorly expressed Hox genes. Hox1-4 paralog expression in the developing central nervous system is essential for linking neurons involved with motor activity to their behavioral output 56 . Hoxb8, identified through RNA-seq as the most significantly differentially expressed gene in the brain due to maternal Cd exposure, has been shown to affect behavior; mutant mice display excessive and pathological grooming behaviors, and Hoxb8 expression has been found in obsessive-compulsive disorder neural circuits in humans 55 . In vitro, direct exposure to Cd has been shown to induce expression of Hoxb8 57 . Additionally, altered Hox expression may explain how maternal Cd exposure can lead to congenital malformations 58,59 , as mutations in Hox genes can lead to drastic developmental malformations 60 . The brain is the most energy expensive organ in the body. Despite comprising only 2% of the body weight, the brain utilizes 20% of the body's energy while at rest. Proper circulating levels of nutrients such as glucose and oxygen are essential to meet this demand 61 . We have previously shown that maternal Cd exposure leads to systemic alterations of essential trace elements in the offspring at birth, including a severe iron deficiency, which suggests that cellular energy pathways may be negatively impacted as a result 15 . The proteomic analysis and reduced mitochondrial content seen in newborn brains here are consistent with perturbed cellular energy and hypoxic conditions due to maternal Cd exposure. Because of the important role of the brain, this organ is www.nature.com/scientificreports/ typically spared during intrauterine growth restriction. In response to hypoxic or nutrient-poor conditions, the fetus adapts by rerouting its circulation to conserve oxygen and nutrient supply for the brain 62 . This brain-sparing phenomenon may at least partially explain why maternal Cd exposure did not significantly alter the raw weight of brains at birth, despite the mice showing a reduction in birth weight 15 . Using a discovery-based metabolomics approach, we identified three significantly increased metabolites in newborn brains due to 50 ppm maternal Cd exposure. Two of these metabolites, 2-hydroxyvaleric acid and 2-hydroxycaproic acid, are both medium chain alpha-hydroxy fatty acids that are not well studied in the context of mammalian neurodevelopment. Elevated levels of 2-hydroxyvaleric acid are associated with lactic acidosis 63 , a condition characterized by excess lactic acid 64 and consistent with our hypothesis of elevated lactate in the brain due to a hypoxic environment caused by maternal Cd exposure. 2-hydroxyvaleric acid is also used in a staining method for measuring tissue-specific lactate dehydrogenase activity 65 , and the elevation in 2-hydroxyvaleric acid would be consistent with the elevated Ldha identified through our proteomic analysis and further supporting our hypothesis regarding elevated lactate. Elevated 2-hydroxycaproic acid levels are seen in response to iron depravation 66 and as a biomarker for cardiovascular mortality 67 , consistent with our previously published work 15 . 2-hydroxycaproic acid is an inhibitor of aminoacylase 1 67 (ACY1); ACY1 is strongly expressed in the brain and its deficiency, though rare, is characterized by impaired neurodevelopment 68 . The third metabolite identified by metabolomics (xanthine) has been studied in greater detail. Xanthine is a purine derivative involved in a number of enzymatic reactions and is potentially toxic 63 . Xanthine accumulates in mitochondrial dysfunction 69 and is produced in the Mitochondrial DNA Depletion Syndrome disease pathway, which is characterized by decreased levels of mtDNA 70 as was seen in this study as a result of maternal Cd exposure. Molybdenum cofactor deficiency is a severe neurological disorder that results from buildup of xanthine due to a nonfunctional molybdenumdependent cofactor that is required for certain enzymes 63,71 ; we previously showed that maternal Cd exposure leads to a significant decrease in molybdenum in newborns 15 , and molybdenum cofactor production may further be impaired by a maternal Cd-induced iron deficiency 15,72 . The maternal Cd-induced increase in xanthine and the increase in hypoxanthine-guanine phosphoribosyltransferase (Hprt, as identified through our proteomic analysis) may be suggestive of an overall impairment or alteration to purine metabolism, which is emerging as an important metabolic pathway for proper brain function and development 73 . Finally, two enzymes involved in xanthine metabolism (xanthine oxidase and xanthine dehydrogenase) may play a role in RA metabolism and synthesis 74,75 .
Using multiple -omic approaches, we identified RA signaling as a candidate pathway that is perturbed by maternal Cd exposure in the developing brain. To test this, we quantified retinoid species and confirmed www.nature.com/scientificreports/ increased levels in newborn brains as a result of maternal Cd exposure. It is possible that xanthine was increased in maternally-exposed brains due to xanthine oxidase and/or xanthine dehydrogenase being sequestered for RA synthesis, though this will need to verified experimentally. RA is required for the onset of myelination 31 and Hox gene expression is controlled in part by a rostrocaudally-decreasing gradient of RA along the body axis 56 . Abnormal levels of RA during development, either excess or reduced, lead to embryological defects, affect hindbrain development, and alter Hox gene expression 55 . Excess RA during development can lead to profound behavioral abnormalities and mental retardation in the offspring 47 . In vitro and adult exposure to Cd has been shown to induce RA signaling, and Cd teratogenicity is suspected to be at least in part due to increased RA as a result of disruption to RA-metabolizing genes 76 , consistent with the data in this study. The transcriptomic, proteomic, and metabolomic analyses of newborn brains are likely to be affected by the cellular heterogeneity of whole brains. Brain heterogeneity can bias sequencing analyses 77 and may exclude region-or cell-specific changes and limit differential expression analysis to abundant transcripts, consistent with the low number of differentially expressed genes identified in this study. Methods such as surrogate variable analysis (SVA) and remove unwanted variation (RUV) can potentially control for sources of heterogeneity 78,79 . However, one of our key findings is that genes involved in myelination are altered in abundance and given that this is a biological process that occurs throughout the brain, we report the global results in this study without implementing a statistical method to control for heterogeneity. The proteomic analysis can only quantify levels of proteins that are abundant in tissue, and therefore can miss detection of proteins found at lower levels in the cell (e.g., transcription factors 80 ). This may explain the lack of overlap between genes identified as differentially expressed through RNA-seq and the proteins found to be differentially abundant. Despite these caveats, we were able to detect an enrichment for RA-associated terms in both the transcriptomic and proteomic datasets, which was supported by the retinoid species quantification. In addition, 5 proteins with roles in cellular energy metabolism that were found to be significantly increased in abundance here (Ldha, Pfkl, Aldoc, Anxa2, Aldoa) were also found to be significantly increased at the level of gene expression in newborn hearts as a result of maternal Cd exposure 15 , suggesting that maternal Cd exposure-associated alterations to cellular energy metabolism may be apparent in multiple organs.
In summary, we have shown global changes to the transcriptome, proteome, and metabolome in newborn mouse brains as a result of maternal Cd exposure and have generated multiple hypotheses for future mechanistic work to link maternal Cd exposure to altered neurodevelopment (Fig. 8). As a proof of principle that our multiple -omics approach can identify novel candidate pathways perturbed by maternal Cd exposure, we demonstrated We propose that maternal Cd exposure increases retinoic acid (RA) synthesis and/or degradation in the fetal brain. This leads to aberrant Hox expression, altered (early onset of) myelination, and altered behavior. Aberrant Hox expression may explain the previously known role of Cd as a teratogen. Additionally, we propose that maternal Cd exposure damages the mitochondria in the fetal brain and therefore impacts energy production in the brain. Metabolites produced during this process may have feedback on RA production. Finally, we propose that maternal Cd exposure alters levels of essential trace elements, which is sensed by the body as a state of malnutrition and may lead to a hypoxic environment during development due to the significantly reduced levels of iron in circulation. This also alters energy production in the brain and may explain the enlarged proportional brain weights at birth, consistent with the phenomenon known as "brain-sparing". These hypotheses lay the foundation for future mechanistic work to elucidate molecular pathways occurring at each arrow. Data described here or in Hudson et al. 15 that support the hypotheses are italicized in smaller text below the corresponding hypothesis in each box. www.nature.com/scientificreports/ maternal Cd-associated changes to retinoid species in the brain, providing empirical support for the alterations to RA signaling suggested by transcriptomic and proteomic data. Overall, our results support findings of other studies that suggest a neurotoxic effect of early life Cd exposure and provide new insights into understanding the link between maternal Cd exposure and impaired neurodevelopment.

Materials and methods
Animal husbandry, Cd exposure, and tissue collection. Animal work and tissue collection were performed as described previously 15 . Briefly, 5-to 7-week-old C57BL/6J ('B') and CAST/EiJ ('C') females (the F 0 generation) were exposed through their drinking water to 0 ppm, 1 ppm, or 50 ppm Cd in the form of CdCl 2 (Sigma-Aldrich, catalog number 202908) for a period of 5 weeks, then mated with a male of the opposite strain. Cd exposure continued throughout pregnancy and was discontinued at litter birth. A hybrid mating scheme (the F 1 generation) was employed to enable analyses of allele-specific gene expression and DNA methylation for a separate study (Hudson & Cowley, manuscript in preparation). The genotypes of F 1 mice are referred to as 'B × C' (B mother × C father) and 'C × B' (C mother × B father). F 1 animals were sacrificed and dissected within 24 h of birth or at 6 months of age; sample sizes can be found in a previous publication that studied the effect of maternal Cd exposure on cardiovascular features in these mice 15  Behavioral assessments. F 1 mice were subjected to two behavioral tests at approximately 5.5 months (+ /-2 weeks) of age over two consecutive days: open field (OF) and social interaction (SI). Behavioral testing was performed as previously described 81,82 . Females were not controlled for stage of estrous cycle. Boxes were disinfected between the testing of each subject. All testing was performed during the afternoon between 3 pm-6 pm under dimmed lighting in a secluded procedure room, video recorded for 30 min, and scored using TopScan software (Clever Sys Inc.). A randomly selected subset of 2 videos per treatment group and per sex was scored by hand to validate TopScan scoring.
SI was performed in the same box one day after OF testing using an unexposed B mouse of the same sex and of similar age. Pure B strain mice that were the same sex as the hybrid F 1 test animal were used as novel test mice in place of B × C or C × B hybrids for the SI experiments due to the increased aggression shown by hybrid male mice to each other.
Sample sizes used in behavioral experiments, source data, and statistical summaries for all behavioral parameters measured can be found in Supplementary Table 1. The number of litters represented in animals subjected to behavioral testing can be found in Supplementary Table 1 of a previous publication 15 , as the animals described here were subjected to blood pressure testing 1-2 weeks after behavioral testing. A one-way analysis of variance (ANOVA) followed by a post-hoc Tukey test using R software was performed on the data from the OF and SI experiments for the four control groups (B × C females, B × C males, C × B females, C × B males) to identify any significant differences in behaviors of unexposed mice. All other statistical analyses were performed using an ANOVA followed by a post-hoc Dunnett's multiple comparison test using R software, comparing animals of the same sex and genetic background from the 1 ppm and 50 ppm groups to 0 ppm controls.
Histology. Newborn whole brains or heads were immediately immersed in 10% neutral buffered formalin (NBF), fixed for 48 h, transferred to a 30% sucrose solution until sunk, then submitted to the NC State's College of Veterinary Medicine Histology Laboratory. There, the brains were either embedded in paraffin and stained with Hematoxylin & Eosin (H&E), or the heads were frozen and cryosectioned for Black Gold II staining.
Formalin fixed paraffin embedded (FFPE) brains were processed routinely, oriented for 5 µm sectioning in the coronal plane, and stained with H&E stain. H&E staining was performed on every fifth serial coronal section taken from the rostral to caudal end to reduce the number of slides processed. A veterinary board-certified pathologist (ACVP) identified anatomical landmarks throughout the brain on the slides in order for accurate comparisons between the same regions across individuals.
Frozen heads were cryosectioned on the coronal plane at a thickness of 20 µm to obtain sections representing the prefrontal cortex. Further processing and black-gold staining was performed on slides containing comparable regions in the prefrontal cortex across individuals according to manufacturer directions (Millipore Black-Gold II Myelin Staining Kit, catalog number AG105).
Nucleic acid isolation. Frozen whole brains obtained from eight 0 ppm, eight 1 ppm, and eight 50 ppm newborn female B × C pups representing four litters per exposure group were pulverized using a hammer and liquid nitrogen. 7-10 mg of tissue was then used to extract DNA and RNA simultaneously using the AllPrep DNA/RNA/miRNA kit (Qiagen, 80204). RNA was extracted using the recommended protocol for fatty tissues. RNA and DNA were quantified using a Nanodrop 2000. RNA purity and size integrity were determined at the NCSU Genomic Sciences Laboratory (GSL) using an Agilent 2100 Bioanalyzer with an RNA 6000 Nano Chip (Agilent Technologies).

RNA-seq.
Total RNA from the brains of four 0 ppm and four 50 ppm maternally exposed female F 1 mice (each individual representing a different litter) were submitted to the NCSU GSL for indexed library construction and sequencing as described previously 15  www.nature.com/scientificreports/ The quality of raw RNA-seq data was assessed using the FastQC application and the initial 10 bases with poor quality were trimmed. Alignment was performed using STAR short read aligner 83 to the respective mouse strain (C57BL/6J, CAST/EiJ) reference genomes. C57BL/6J and CAST/EiJ RNA-seq data were aligned to mm38 version 87 and CAST/EiJ version 1.86 reference genome downloaded from the Ensembl website, respectively. The number of reads mapped to a genome feature was determined using htseq-count script from HTSeq python package 84 . Independent analyses were conducted for the two strains. The two raw count data matrices were imported to R statistical computing environment for further analysis 85 . Genes that had no count in more than 2/3 of the replicate samples were excluded from the analysis. Data normalization based on dispersion and differential analysis was conducted using the DESeq2 R package 86 . A generalized linear model (~ Treatment) was fitted between the expression count and a treatment group (0 ppm or 50 ppm maternal Cd exposure). Finally, a list of differentially expressed genes was generated between 50 ppm vs. 0 ppm maternal Cd exposure after applying Benjamini-Hochberg multiple testing correction 87 (p adj < 0.05). The entire dataset of all genes with detectable and good quality reads aligned to both genomes can be found in Supplementary Table 2. The heatmap of the 15 significantly differentially expressed genes shown in Fig. 3 was generated through http:// www. heatm apper. ca/ using average linkage clustering method and Pearson distance measurement method.
Data generated by RNA-seq in this study are deposited in Gene Expression Omnibus (GEO), accession number GSE174433.
qRT-PCR and qPCR. 100 ng of total RNA extracted as described above from 24 B × C female newborn F 1 brains was used to synthesize first strand cDNA as described previously 15 ; the 8 samples used in RNA-seq were included in these 24 samples. qRT-PCR was performed on a Real-Time 7300 machine (Applied Biosystems) using SsoAdvanced Universal SYBR Green Supermix (Bio-Rad, 1725271) using the conditions and methods of analysis described previously 15 . The primer sequences used for the Hox genes in Fig. 4 are provided in Supplementary Table 6. Additional Hox genes were tested, but did not have consistent detectable amplification in at least one of the three maternal Cd dose groups. We defined consistent detectable amplification as at least 3 of the 8 pups' cDNA samples per treatment group generating C t values for each technical replicate (i.e. for one pup, generating three C t values for the three technical replicates, and this occurring for at least 3 out of 8 pups in a group).
The ΔΔC t method of qRT-PCR gene expression analysis requires a comparison between a sample's ΔC t value to the average ΔC t value of the control group, where ΔC t = C t (gene of interest) − C t (reference gene). Due to undetectable amplification of several Hox genes in multiple control samples (which resulted in undetermined C t values in up to 6 of the 8 control samples), the average ΔC t value of the 1 ppm group was used as the "control" ΔC t value for calculating the ΔΔC t value. The expression values are therefore given relative to the 1 ppm group (Fig. 4). Samples with undetermined C t values were then given a relative expression value of 0%.
DNA was extracted from 24 B × C female newborn F 1 brains through the Qiagen kit described above. Mitochondrial DNA (mtDNA) was quantified by qPCR as previously described 88,89 .
To provide greater transparency and allow for better reproducibility of results, this study was performed in compliance with MIQE standards 90 (see Supplementary Table 7).
Proteomics. Frozen homogenized whole brains obtained from newborn B × C maternally exposed female pups were submitted to the Molecular Education, Technology, and Research Innovation Center (MET-RIC) at NCSU for proteomic analysis. The brain samples used in the proteomic analysis originated from the four 0 ppm and four 50 ppm females represented in the RNA-seq analysis; additionally, four 1 ppm females representing four different litters were used for proteomic analysis.
Samples were quantitated for protein content using a Pierce bicinchoninic acid (BCA) kit and normalized to 50 μg of protein by diluting the appropriate amount of sample in 50 mM ammonium bicarbonate with 1% (w/w) sodium deoxycholate. Normalized samples were reduced with dithiothreitol (DTT) and alkylated with iodoacetamide (IAM) to break disulfide bonds and prevent reformation. Following this, a filter aided sample preparation (FASP) protocol 91,92 was used to purify the samples which were then treated with trypsin at a 50:1 protein:trypsin concentration. Samples were incubated for 4 h and then aliquoted for LC-MS analysis.
Chromatographic separation was achieved using a Thermo Scientific EASY nLC 1200 system (Bremen, Germany). Pico-frit columns were purchased from New Objective (Woburn, MA) and bomb packed to a length of 20-30 cm with reverse phase ReproSil-Pur 120 C-18-AQ 3 µm particles (Dr. Maisch, Germany). Two microliters of sample was injected and subsequently separated using a gradient of mobile phase A (98% water, 2% acetonitrile, and 0.1% formic acid) and mobile phase B (80% acetonitrile, 20% water 0.1% formic acid). The LC method consisted of a 120 min gradient with a linear ramp from 0% B to 40% B, a 1 min ramp to 100% B which was held for 6 min (123-129 min), followed by equilibration of the column at 0% B (130-140). A flow rate of 300 nL/min was used.
Orbitrap tandem mass spectrometry was performed using a Thermo Scientific Q-Exactive HF (Bremen, Germany) in a top 20 data dependent acquisition (DDA) mode, where the 20 most abundant precursors are selected for fragmentation per full scan. MS1 and MS2 scans were performed at a resolving power of 120,000 and 15,000 at m/z 200, respectively. A dynamic exclusion window of 20 s was used to avoid repeated interrogation of abundant species. Automatic gain control (AGC) was 1e6 and 1e5 for MS1 and MS2 scans, respectively. Samples were run in random order, and a quality control BSA digest was run every fifth injection to ensure proper LC-MS/MS reproducibility.
Resulting raw data were loaded into Proteome Discoverer (version 2.0). Spectrum files were run through the spectrum selector node to appropriately identify peaks and this data was collated using the Sequest HT node and aligned with a FASTA file which contained all proteins indexed by Uniprot and assigned a taxonomy ID of www.nature.com/scientificreports/ 10,090 (mus musculus). Oxidation (dynamic) and carbamidomethylation (static) modifications were considered in Sequest HT with a maximum of 2 potential missed cleavage sites and peptide lengths between 6 and 144 amino acids. The Percolator node was used for false discovery rate (FDR) calculations. Data was then filtered to remove the following: abundance ratio adj. p-value (50 ppm vs 0 ppm) > 0.05, # unique peptides > 1, Abundances (grouped) CV% 50 ppm > 20, Abundances (grouped) CV% 0 ppm > 25, proteins with a combination of low PEP score and CV% > 15, PEP score < 7, and proteins that had a 'Found in Sample' value of any value lower than 'High' for at least 2 individuals.
Proteomics raw data can be found in Supplementary Table 3. Lactate quantification and metabolomics. Four homogenized newborn B × C female brains were tested for each of the three treatment groups, and each sample was tested twice. One 50 ppm sample was destroyed during processing, resulting in a final n = 3 for that group. To each mouse brain sample, 100 µL of acidified mobile phase containing ITSD (50 ug/mL 13C3-sodium lactate, Cambridge Isotope Laboratories) and 500 µL of cold methanol were added prior to homogenization using an Omni tissue homogenizer with disposable tips. After 10 s of homogenization, the samples were allowed to sit at −20 °C for 30 min before centrifugation for 15 min (4 °C, 3000×g) for protein precipitation. The supernatant was removed and transferred to a fresh 1.5 mL Eppendorf tube and taken to dryness on a SpeedVac. The samples were then resuspended in 100 µL mobile phase A and centrifuged for 15 min (4 °C, 3000×g). The supernatant was transferred directly to an LCMS autosampler vial for analysis. Quality control (QC) samples were prepared by pooling and mixing equal volumes of each sample. The QC sample and blank samples were injected at regular intervals throughout the sequence. For quantification of lactate, a 100 ug/mL stock solution of sodium lactate (Cambridge Isotope Laboratories) was prepared using acidified mobile phase containing ITSD (50 ug/mL 13C3-sodium lactate, Cambridge Isotope Laboratories). Seven calibration standards ranging from 3.2 ng/mL to 50 µg/mL were prepared by serially diluting the working standard solution with acidified mobile phase containing ITSD. Chromatographic separation was achieved using a Waters Acquity BEH C18 column (2.1 × 100 mm) holding for 5 min at 100% mobile phase A (5% MeOH in water + 0.15% formic acid) followed by gradient elution (mobile phase B 100% MeCN) to 95% mobile phase B (flow rate 350 µL/min, 2 µL injections in duplicate). High resolution mass spectrometry data was acquired using a Thermo Orbitrap ID-X mass spectrometer in negative mode (spray voltage 3.0 kV, vaporizer temperature 300 °C) with a mass range of m/z 70-800 and a dwell time of 0.6 s based on the chromatographic peak width of 6 s allowing 10 scans across each peak for accuracy in quantification. MS 1 data was collected with a resolution of 120,000 and an AGC target of 2.0e5. MS 2 data was collected for global metabolite profiling (AGC target 1.0e5, resolution 30 K, and stepped HCD collision energy of 20,35,50).

RNA-seq and proteomic enrichment analyses.
Targeted data processing. Peak integration and quantification was performed in TraceFinder 4.1. A standard curve for lactate was constructed using extracted ion chromatogram peak areas from MS1 data and the slope of the curve was calculated in reference to the ISTD using the area of the peak divided by the internal standard peak area with a 1/x weighting. The concentration of lactate in the study samples was calculated in an identical manner relative to the regression line. The lactate calibration curve had an R 2 value of 0.9994 for the linear range of 80 ng/mL to 50 ug/mL.
Untargeted data processing. Preprocessing and annotation in CD 3.0. Raw data files were uploaded into Compound Discoverer 3.0 (Thermo Fisher Scientific) and processed using a workflow to find and identify differences between samples. This workflow performed retention time alignment, unknown compound detection and compound grouping across all samples. Elemental compositions were predicted for all compounds and chemical background was hidden using blank samples. Compound annotations were assigned using ddMS2 data in comparison with the mzCloud database (Thermo Fisher Scientific) as well as molecular formula comparisons with selected ChemSpider databases (KEGG, HMDB, Mass Bank) which were further ranked using the mzLogic algorithm (Thermo Fisher Scientific).
Compound annotations. A total of 1,471 compounds were identified in the data set. Compounds were filtered in order to exclude background resulting in 1,217 compounds of which 162 were annotated in Compound Discoverer. Manual curation of these 162 compounds to remove duplicate and erroneous annotations resulted in 101 compounds.
Statistical analysis in MetaboAnalyst. The filtered data set (99 compounds) with peak areas per file were formatted for further statistical analysis in MetaboAnalyst 95 . Sample normalization, data transformation, and data scaling were set to normalization by sum, log transformation, and pareto scaling, respectively. Results were analyzed using univariate (T-test/ANOVA, Fold Change, Volcano Plot, Correlation Analysis) and multivariate (PCA, PLS-DA) statistics as well as by hierarchical clustering (Heatmap) and supervised feature selection (Random Forest). Normalized data for each sample was exported from MetaboAnalyst into an excel file and the mean value of duplicate values for each sample for the 99 compounds was calculated. Outliers, ANOVAs, and post-hoc Dunnett's test were calculated as described below in the Statistical Analysis section, with a final step Retinoid quantification. Whole brains of newborn B × C maternally-exposed pups were removed within 30 s of pup sacrifice under yellow light, placed in a black 1.7 mL microcentrifuge tube to prevent exposure to white light, then immediately weighed and flash frozen on dry ice. Sample sizes were as follows: sixteen 0 ppm females, five 1 ppm females, eleven 50 ppm females, eleven 0 ppm males, seven 1 ppm males, and eight 50 ppm males. Retinoid concentrations were measured as described previously [96][97][98] . Briefly, an average of 85.2 ± 9.2 mg (mean ± s.d.) (range 67.1-117.1 mg) per sample of brain was homogenized in ground glass homogenizers in 1.0 mL normal saline (0.9% NaCl). Extraction of retinoids was performed under yellow lights using a twostep liquid-liquid extraction that has been described in detail previously using 4,4-dimethyl-RA as an internal standard for RA and retinyl acetate as an internal standard for retinol and total retinyl ester 96,98 . RA levels were determined by liquid chromatography-multistage tandem mass spectrometry (LC-MRM3), which is an LC-MS/ MS method utilizing two distinct fragmentation events for enhanced selectivity 96 . RA was quantified using a Shimadzu Prominence UFLC XR liquid chromatography system (Shimadzu, Columbia, MD) coupled to an AB Sciex 6500 + QTRAP hybrid triple quadrupole mass spectrometer (AB Sciex, Framingham, MA) using atmospheric pressure chemical ionization (APCI) operated in positive ion mode as previously described 96 . Retinol and retinyl ester were quantified via HPLC-UV using an AQUITY H-Class UPLC with a PDA detector (Waters Corporation, Milford, MA) operated in single wavelength monitoring mode at 325 nm according to previously published methodology 97,98 . Retinoids in tissues are expressed as mol/g tissue.

Statistical analysis.
Unless otherwise noted, all statistical analyses were performed using an ANOVA followed by a post-hoc Dunnett's multiple comparison test using R software (version 4.0.2, https:// cran.r-proje ct. org), comparing animals of the same sex and genetic background from the 1 ppm and 50 ppm groups to 0 ppm controls. Outliers, as defined by a point which falls more than 1.5 times the interquartile range above the third quartile or below the first quartile, were omitted from the analysis if found to be present. Data are presented as the mean ± standard error of the mean. *p < 0.05, **p < 0.01, ***p < 0.001. Source data from each figure and data referenced in the main text can be found in Supplementary Table 8.