Article pubs.acs.org/crt
Benchmark Dose Modeling Estimates of the Concentrations of Inorganic Arsenic That Induce Changes to the Neonatal Transcriptome, Proteome, and Epigenome in a Pregnancy Cohort Julia E. Rager,*,† Scott S. Auerbach,‡ Grace A. Chappell,† Elizabeth Martin,§ Chad M. Thompson,∥ and Rebecca C. Fry§ †
ToxStrategies, Inc., Austin, Texas 78759, United States National Toxicology Program, National Institutes of Health, Research Triangle Park, North Carolina 27709, United States § Department of Environmental Sciences and Engineering, The University of North Carolina at Chapel Hill, Chapel Hill, North Carolina 27516, United States ∥ ToxStrategies, Inc., Houston, Texas 77494, United States ‡
S Supporting Information *
ABSTRACT: Prenatal inorganic arsenic (iAs) exposure influences the expression of critical genes and proteins associated with adverse outcomes in newborns, in part through epigenetic mediators. The doses at which these genomic and epigenomic changes occur have yet to be evaluated in the context of dose−response modeling. The goal of the present study was to estimate iAs doses that correspond to changes in transcriptomic, proteomic, epigenomic, and integrated multi-omic signatures in human cord blood through benchmark dose (BMD) modeling. Genome-wide DNA methylation, microRNA expression, mRNA expression, and protein expression levels in cord blood were modeled against total urinary arsenic (U-tAs) levels from pregnant women exposed to varying levels of iAs. Dose−response relationships were modeled in BMDExpress, and BMDs representing 10% response levels were estimated. Overall, DNA methylation changes were estimated to occur at lower exposure concentrations in comparison to other molecular endpoints. Multi-omic module eigengenes were derived through weighted gene co-expression network analysis, representing co-modulated signatures across transcriptomic, proteomic, and epigenomic profiles. One module eigengene was associated with decreased gestational age occurring alongside increased iAs exposure. Genes/proteins within this module eigengene showed enrichment for organismal development, including potassium voltage-gated channel subfamily Q member 1 (KCNQ1), an imprinted gene showing differential methylation and expression in response to iAs. Modeling of this prioritized multi-omic module eigengene resulted in a BMD(BMDL) of 58(45) μg/L U-tAs, which was estimated to correspond to drinking water arsenic concentrations of 51(40) μg/L. Results are in line with epidemiological evidence supporting effects of prenatal iAs occurring at levels 1/3 the lowest dose tested in order to avoid artificial minimization of BMDs. Other parameters included setting the maximum iterations to 250 and the confidence level to 0.95. Probes were removed if they showed either BMD > the highest dose, BMD/BMDL (BMD lower confidence limit) ratios > 20, and/or BMD < the lowest dose in order to avoid model extrapolation. The model curves were required to have goodness-of-fit p > 0.1 (likelihood ratio test) to be selected as adequately describing dose−response trends in methylation/expression.12,16 Molecules that met these statistical criteria were identified as showing dose-dependent changes associated with U-tAs, in addition to their previous statistical association to U-tAs via multivariate linear regression modeling. A Spearman rank correlation test was used to analyze the potential relationship between the number of dosedependent molecular endpoints and median BMD estimates. Statistical Evaluation of Molecules in Association with Birth Outcomes. All molecules associated with U-tAs level were also evaluated for statistical association with birth outcomes using multivariable regression models similar to those described above. The evaluated birth outcomes included newborn weight, gestational age, weight/gestational age ratio, height, head circumference, and placental weight. These models used the specific -omic profile under evaluation as the dependent variable with each birth outcome of interest (run separately) as the independent variable, while adjusting for mother’s age (continuous variable), mother’s smoking status (categorical variable), alcohol consumption (binary variable), and UtAs (continuous variable). Statistical analyses were carried out using R (v3.2.4) using the fit linear models (lm) and adjust p-values for multiple comparisons (p.adjust) functions. Molecules with a false discovery rate less than 5% (q-value 0.80. This resulted in the construct of a co-expression network that was weighted, emphasizing high correlations and de-emphasizing low correlations between molecules. Modules were then identified as groups of densely interconnected molecules in the weighted network analysis with high topological overlap, measured based on an average linkage hierarchical clustering with a dynamic tree-cutting algorithm. Simply put, the resulting modules represented clusters of highly interconnected molecules with high positive or negative correlations. With these defined modules, module eigengenes were then calculated as the first principal component of a given module.27 In the current analysis, all molecules that were identified as significantly associated with U-tAs that also were annotated to known genes/proteins (total number of molecules = 3762) were analyzed collectively using the WGCNA R package.30 Data from 31 subjects were used here, as these represented the subset of subjects with all four -omics endpoints evaluated. WGCNA resulted in the identification of module eigengenes, which were ultimately utilized for the purposes of BMD modeling. Because module eigengenes are defined as the first principal component of a given module, in the current analysis, module eigengenes represented the collective methylation and/or expression levels of hundreds-to-thousands of molecules that are comodulated in response to iAs. Resulting module eigengenes were analyzed through BMDExpress using the aforementioned methods and correlated against arsenic measures and birth outcome data using Pearson correlation tests (WGCNA R package). Network Analysis of Molecules within Prioritized Module Eigengene. Network analysis was carried out to understand systemslevel responses that occur within a co-modulated set of multi-omic molecules identified through WGCNA. For this analysis, iAsassociated genes and proteins from a prioritized module eigengene were overlaid onto a global molecular interaction network based on published molecular interactions, and networks were algorithmically constructed based on connectivity, as enabled through Ingenuity Pathway Analysis (Ingenuity Systems). Sets of genes/proteins corresponding to biological functions and disease signatures within the constructed networks were then identified using a right-tailed Fisher’s exact test and corrected for multiple testing using the Benjamini−Hochberg procedure, as performed previously.17,19
Table 1. Characteristics of the Study Cohort na (%) or mean, median [range] exposure measures DW iAs (μg As/L) UtAs, SG adj (μg/L) subject characteristics maternal age at delivery (years) alcohol: no consumption alcohol: some consumption smoking status: nonsmokers smoking status: current smokers newborn characteristics newborn sex: male newborn sex: female gestational age (weeks) birth weight (g) APGAR score placental weight (g) length (cm) head circumference (cm)
49.9, 30.5 [ 0.9), suggesting that the general trends observed for BMD values were indicative of biological response sensitivity to dose rather than merely the number of times an endpoint was assessed. BMD Modeling of Multi-omic Parameters through WGCNA. Five module eigengenes were identified through WGCNA (Figure 2A), representing the collective methylation and/or expression levels of hundreds-to-thousands of molecules that were co-modulated in response to iAs. Four of these module eigengenes showed dose−response curves that met our statistical criteria (Supporting Information, Table S1). One 1915
DOI: 10.1021/acs.chemrestox.7b00221 Chem. Res. Toxicol. 2017, 30, 1911−1920
Article
Chemical Research in Toxicology
Figure 2. Module eigengenes representing co-modulated molecules associated with U-tAs across multi-omic signatures. (A) Correlation matrix showing that module eigengenes (MEs) were highly correlated with several measures of arsenic exposure, including U-tAs (left half). Module eigengene correlations to birth outcomes are also shown (right half). (B) BMD curve fit for the prioritized module eigengene (“MEgrey”) displaying the eigengene data points grouped by quintiles (red squares) with the best fitting curve plotted as a blue line. BMD estimates are represented by the black vertical lines closest to the right, and BMDL estimates are represented by the black vertical lines closest to the left.
this module contained KCNQ1 methylation and mRNA expression co-modulations, as well as other example genes of interest identified through network analysis as involved in organismal growth signaling that were co-modulated in response to iAs (Figure 3A,B). Network Analysis of Multi-omic Responses in Prioritized Module Eigengene. Network analysis of the molecules within the prioritized module eigengene (MEgrey) generated multiple significant networks containing known interactions between proteins encoded by differentially methylated/expressed genes, miRNAs, and other molecules (Supporting Information, Table S3). One of the most highly enriched diseases/functions associated with these network-level
module eigengene (“MEgrey”) was significantly correlated with gestational age (p = 0.005), whereas the other module eigengenes (“MEblue”, “MEyellow”, “MEbrown”, and “MEturquoise”) were not as significantly associated with birth outcomes (p > 0.01) (Figure 2A). Importantly, gestational age is a birth outcome measure that was previously identified to be significantly associated with prenatal iAs exposure in the BEAR cohort.18 The MEgrey module eigengene of interest was, thus, prioritized for further analysis and was estimated to have a BMD(BMDL) of 58(45) μg/L U-tAs (Figure 2B). The comodulation of 984 molecules was identified in this module, including 514 CpG methylation sites, 12 miRNAs, 344 mRNAs, and 114 proteins (Supporting Information, Table S2). Notably, 1916
DOI: 10.1021/acs.chemrestox.7b00221 Chem. Res. Toxicol. 2017, 30, 1911−1920
Article
Chemical Research in Toxicology
Figure 4. Network showing organismal development signaling associated with molecules within the prioritized multi-omic module eigengene (“MEgrey”). Molecules with iAs-associated changes are colored, and molecules with associated signaling are white.
profiles, median BMD(BMDL) estimates for CpG site methylation, miRNA expression, mRNA expression, and protein expression changes were estimated to occur at drinking water arsenic levels of 31(16), 38(18), 56(28), and 40(14) μg As/L water, respectively. The prioritized multi-omic module eigengene was estimated to have a BMD(BMDL) of 58(45) μg/L U-tAs, which was estimated to correspond to DW-iAs concentrations of 51(40) μg/L.
Figure 3. Example molecules of interest showing concerted responses to prenatal iAs exposure, grouped within the prioritized module eigengene (“MEgrey”). Example molecules are shown that had methylation and/or expression levels that were (A) positively or (B) negatively correlated to MEgrey, a collective measure of methylation and/or expression levels of 984 co-regulated molecules. Methylation and expression levels are Z-score normalized.
■
DISCUSSION This study evaluated dose−response relationships for transcriptomic, proteomic, and epigenomic responses to prenatal iAs exposure in newborn cord blood from a pregnancy cohort in Gómez Palacio, Mexico, with the goal of estimating the exposure doses at which important molecular changes occur. A multi-omic parameter (i.e., module eigengene) was also modeled that included measures from >900 highly correlated genes/proteins. Analyses resulted in the estimation of BMD(BMDL) values for individual iAs-responsive molecules as well as trends in overall -omic changes. A variety of BMD estimates were provided to serve as a starting point for estimating doses at which prenatal iAs exposure may influence important neonatal transcriptomic, proteomic, and epigenomic signatures. The BMD modeling results showed that, together, CpG methylation changes were estimated to occur at a median BMD(BMDL) of 33(15) μg As/L urine; miRNA expression at 42(18) μg As/L urine; mRNA expression at 64(30) μg As/L urine; and protein expression at 45(12) μg As/L urine. The multi-omic module eigengene that was prioritized for its correlation to gestational age was estimated to have a BMD(BMDL) of 58(45) μg/L U-tAs. These BMD(BMDL) values were estimated to relate to DW-iAs at the following levels: CpG methylation changes at 31(16) μg As/L; miRNA expression changes at 38(18) μg As/L; mRNA expression changes at 56(28) μg As/L; protein expression changes at 40(14) μg As/L; and the multi-omic module eigengene at 51(40) μg/L water. To compare these BMD(BMDL) estimates to existing recommendations, the current WHO recommended limit and U.S. EPA limit for arsenic in drinking water is 10 μg/
responses was organismal development (Benjamini−Hochberg p < 1.6 × 10−10) (Supporting Information, Table S4). A highly significant network related to organismal development was among those that were constructed based on known molecular interactions (Figure 4) and contained 15 molecules related to this functional category. Of these organismal developmentrelated molecules, KCNQ1 had CpG island methylation levels also showing a significant association with several birth outcomes, including gestational age (Supporting Information, Table S1). Other molecules related to organismal development and associated with iAs exposure included neuronal differentiation 1 (NEUROD1), a protein with increased expression; histone deacetylase 5 (HDAC5), a gene with decreased cytosine methylation; and HDAC4, a gene with both increased cytosine methylation and decreased mRNA expression. HDAC4 and HDAC5 methylation was also associated with gestational age (Supporting Information, Table S1). Estimating Drinking Water Levels Associated with BMD Results. The BMD analysis was carried out by modeling molecular endpoints against U-tAs, as U-tAs is recognized as the most reliable indicator of human arsenic exposure occurring within the last several days.31 The amount of iAs in drinking water that correlates with the U-tAs BMD estimates for the prioritized gene of interest was estimated using a linear regression model relating measured U-tAs versus DW-iAs levels. Evaluating overall exposure-associated changes in -omic 1917
DOI: 10.1021/acs.chemrestox.7b00221 Chem. Res. Toxicol. 2017, 30, 1911−1920
Article
Chemical Research in Toxicology L,32,33 and the current U.S. EPA oral reference dose (RfD) is 0.3 μg As/kg·day, based on a no observed adverse effect level (NOAEL) of 9 μg As/L water for hyperpigmentation, keratosis, and possible vascular complications (last revised in 1991).34 These exposure limits were largely based on health effects observed in adults, and the U.S. EPA arsenic regulations are in the process of being updated,35 making the results of the current analysis on the effects of arsenic exposure in newborns of timely importance. BMD modeling findings are in line with previous epidemiological evidence supporting adverse health effects associated with prenatal exposure to arsenic starting at exposures as low as 50 μg iAs/L drinking water. For example, an investigation of a cohort in Bangladesh that included >29 000 pregnancies found that significantly increased risks of fetal loss and infant death occurred at concentrations >50 μg iAs/L water,36 and it reported a dose-dependent relationship between iAs exposure and infant death.36 Another study evaluated 1578 mother−infant pairs from Bangladesh and related birth outcomes to U-tAs.37 At a low level of iAs exposure (
