Arsenic- and cadmium-induced toxicogenomic response in mouse embryos undergoing neurulation

Abstract

Arsenic (As) and cadmium (Cd) are well-characterized teratogens in animal models inducing embryotoxicity and neural tube defects (NTDs) when exposed during neurulation. Toxicological research is needed to resolve the specific biological processes and associated molecular pathways underlying metal-induced toxicity during this timeframe in gestational development. In this study, we investigated the dose-dependent effects of As and Cd on gene expression in C57BL/6J mouse embryos exposed in utero during neurulation (GD8) to identify significantly altered genes and corresponding biological processes associated with embryotoxicity. We quantitatively examined the toxicogenomic dose-response relationship at the gene level. Our results suggest that As and Cd induce dose-dependent gene expression alterations representing shared (cell cycle, response to UV, glutathione metabolism, RNA processing) and unique (alcohol/sugar metabolism) biological processes, which serve as robust indicators of metal-induced developmental toxicity and indicate underlying embryotoxic effects. Our observations also correlate well with previously identified impacts of As and Cd on specific genes associated with metal-induced toxicity (Cdkn1a, Mt1). In summary, we have identified in a quantitative manner As and Cd induced dose-dependent effects on gene expression in mouse embryos during a peak window of sensitivity to embryotoxicity and NTDs in the sensitive C57BL/6J strain.

Figure 1. Analysis of Genes Altered by Arsenic and Cadmium

This flow diagram represents the statistical approach used to explore dose- and time-dependent effects of As and Cd exposure on gene expression and corresponding GO biological processes. Linear models assessing for dose (As or Cd) and time effects were used to identify genes that were significantly altered by metal exposure (ANOVA, B_Dose, p<0.0001). Enrichment analyses of GO biological processes were conducted for all genes significantly differentially expressed due to As or Cd exposure to identify biological processes altered by metals. Quantitative analysis of enriched GO biological processes was employed using GO-Quant, enabling dose- and time-dependent calculations in genes within GO subsets. In addition, K-means cluster gene analysis was used to assess unique dose-response relationships due to As (or Cd) exposure and corresponding enrichment of GO biological processes. We assessed the 116 genes altered by both As and Cd for unique dose-response relationships using K-means cluster gene analysis, and we determined the quantitative relationship within these clusters across dose and time (average FC). Enrichment of transcription factor binding sites and associated GO biological processes were also evaluated.

Figure 2. Dose Dependent Metal-Induced Gene Expression Alterations

(A) The venn diagram illustrates 1960 and 775 genes identified to be significantly altered by As and Cd, respectively (F-test, p<0.0001), and overlapping genes intersecting between these two populations (116 total). Using principal components analysis (PCA), we determined similarities in gene expression alterations between exposure groups across all genes altered by As and/or Cd (B). Within genes identified to be significantly altered by As (C) and Cd (D) (F-test, p<0.0001), post-hoc analyses (t-test) were completed to identify the total number of differentially expressed genes with each dose group in comparison with their respective concurrent control (8h or 12h). The percentage of genes differentially expressed (t-test, p<0.05) are displayed in comparison to the total number of genes identified to be significantly altered across all dose groups (black bar). Additionally, the proportion of up (red) and down (green)-regulated genes are indicated by color.

Figure 3. Enrichment of GO Biological Processes Altered by As and Cd

Using a cross scatter-plot, we illustrate overlapping and non-overlapping enriched GO biological processes altered by As and Cd. Genes significantly impacted by As and Cd (p<0.0001), were separately analyzed for enrichment of GO biological processes (MAPPFinder). Assessing over 4500 GO biological processes, we identified 159 (As) and 38 (Cd) GO biological processes meeting our criteria of enrichment, (p<0.02), Z score >2, and ≥3 genes altered). Using the product of the −log (p-value) and Z score, we obtained a measure of ranking the significance of GO categories to allow direct comparisons between As and Cd. The dashed line represents the minimum criteria to be identified as significant (~2.6). Symbols are differentiated by GO hierarchy. Overlapping altered GO biological processes are represented as having enrichments scores >2.6 for both As and Cd (Quadrant I). Enriched categories unique to As or Cd are identified in Quadrant II and IV, respectively.

Figure 4. Quantitative Analysis of Dose- and Time-Dependent Response in Enriched GO Biological Processes Impacted by As and/or Cd

Selected enriched GO-based biological processes/pathways were assessed across dose and time using GO-Quant. We quantitatively describe genes across dose and time within each GO category by showing the absolute average FC of impact for As (A) or Cd (B). The absolute average FC associated with all genes impacted by As or Cd (p<0.0001) is also included to make comparisons with all genes impacted. Overlapping enriched GO biological processes are labeled by the same symbol and color between figures A and B.

Figure 5. Enriched GO Biological Processes Altered by As Dependent on Dose-Response Relationship

K-means cluster analysis was conducted on all 1960 genes found to be significantly altered by As (p<0.0001) (TIGR MEV). Two relationships were determined, Cluster I and II. Separate GO analyses were conducted for each cluster. Significant GO gene categories (Biological Process Level 4) determined using both MappFinder (M) and/or DAVID (D) are shown. Enriched categories were determined by the following criteria: p-value (<0.02), Z score >2 and a minimum of 3 genes for MappFinder and p<0.05, >2 genes altered for DAVID. The total number of genes altered within each GO category as determined by MappFinder or DAVID is shown in parentheses. The absolute average ratio of gene expression response for genes within each enriched GO category was calculated using GO-Quant across dose and time and depicted by color using TIGR MEV (Figure 4C). Color coding was used to differentiate between GO biological process families within the GO hierarchy and to make comparisons between As and Cd (Figure 4 and 5).

Figure 6. Enriched GO Biological Processes Altered by Cd Dependent on Dose-Response Relationship

K-means cluster analysis was conducted on all 775 genes identified to be significantly altered with Cd. In general, four relationships were determined, clusters I-IV. Separate GO analyses were conducted for each cluster. Significant GO gene categories (Biological Process Level 4) determined using both MappFinder (M) and/or DAVID (D) are shown. Enriched categories were determined by the following criteria: p-value (<0.02), Z score >2 and a minimum of 3 genes for MappFinder and p<0.05, >2 genes altered for DAVID as indicated by letters M or D. The number of genes altered within each GO category as determined by MappFinder or DAVID are shown in parentheses. Color coding was used to differentiate between GO biological process families within the GO hierarchy and to make comparisons between As and Cd (Figure 4 and 5).

Figure 7. Dose Dependent Gene Expression Alterations Altered by both As and Cd

Genes identified to be altered by As and Cd (p<0.0001) were analyzed using K-means clustering analysis across dose and time. Three dose-response relationships were identified (I-III) (A). Quantitative analysis of each cluster was conducted by calculating the average FC associated with all genes within each cluster for As (B) or Cd (C) across dose and time. Transcription factor binding motifs (TFBS) were identified using oPPOSUM (D). TFBS (Z score >2.5, ≥10 genes represented) were selected as enriched (*). Approximately 70% of commonly altered genes were linked with the TFBS database. In addition, TFBS analysis was conducted within each cluster (I-III) to identify if overrepresented TFBS were enriched within each cluster. Associations with selected GO biological processes were determined using DAVID. All enriched TFBS and GO associations are displayed for each gene directly to the right of its name.