2001 Progress Report: Molecular Characterization of a Biological Threshold in Developmental ToxicityEPA Grant Number: R827445
Title: Molecular Characterization of a Biological Threshold in Developmental Toxicity
Investigators: Knudsen, Thomas B. , Charlap, Jeffrey H. , Craig, Robert C.
Current Investigators: Knudsen, Thomas B.
Institution: Jefferson Medical College
EPA Project Officer: Hahn, Intaek
Project Period: October 1, 1999 through September 30, 2002 (Extended to September 30, 2003)
Project Period Covered by this Report: October 1, 2000 through September 30, 2001
Project Amount: $207,170
RFA: Children's Vulnerability to Toxic Substances in the Environment (1999) RFA Text | Recipients Lists
Research Category: Children's Health , Health Effects , Human Health , Health
Objective:The goal of this research project is to characterize, in molecular terms, the biological basis of a threshold in developmental toxicity and place this into context of a quantitative dose-response model for risk assessment. Vulnerability to xenobiotic exposure is partly determined by the pattern of genes expressed in target cells. By profiling embryonic gene expression with cDNA microarrays, critical changes at the low end of the dose response curve for developmental toxicity will be visualized. The research prototype is the purine analogue 2CdA (2-chloro-2'-deoxyadenosine), an ocular teratogen that increases the risk for microphthalmia in early mouse embryos in a manner predetermined by the p53 tumor suppressor genotype and p53 protein induction at 4.5 hours after exposure. The risk for 2CdA-induced microphthalmia is substantially lowered when pregnant dams are co-treated with PK11195, a potent antagonist of the mitochondrial peripheral-type benzodiazepine receptor (Charlap, et al., 2002). Using functional genomics and computational biology, this study focuses on enumerating and classifying responses of the early mouse embryo to teratogen exposure with respect to key benchmarks of dose and time.
Progress Summary:Studies during the first year of this grant focused on the dose-dependent effects of 2CdA exposure on day 8 mouse embryos. Benchmark doses (BMD) were calculated from microphthalmia rates fit to a probit model using BMDS software (www.epa.gov/ncea/bmds): 1.5 mg/kg is the highest dose with no observable adverse effect (NOAEL); 2.5 mg/kg is the estimated dose for an extra 5 percent risk (BMD5) of microphthalmia; and 5.0 mg/kg induces microphthalmia in 24.4 percent of CD-1 mouse fetuses. Test doses 1.25- and 2.5 mg/kg flank the threshold for developmental toxicity, which is modeled as the BMDL (replaces NOAEL) = 2.0 mg/kg. During the second year of this grant, microarray data were obtained for the headfold of day 8 mouse embryos collected at three time intervals after low-dose 2CdA exposure. Treatment of pregnant CD-1 mice with 2.5 mg/kg 2CdA on day 8 of gestation provided exposed embryos. It should be emphasized that this dose, unlike 5.0 mg/kg, causes little or no apoptosis in the early embryo. Dams exposed to 2.5 mg/kg 2CdA and 4.0 mg/kg PK11195 provided embryos for intervention samples. Data were collected for two sets of experiments, one comparing 2CdA-exposed samples (test) with controls (reference) and the other comparing 2CdA + PK11195 co-treated samples (test) with 2CdA alone (reference). When PK11195 treatment was delayed beyond the time of 2CdA exposure, 4.5 hours emerged as the longest delay that would still rescue fetuses from an increased risk for microphthalmia. Therefore, RNA was isolated from the headfold (4-6 somite pair stage) at 3.0 hours, 4.5 hours, and 6.0 hours post-exposure. Control RNA samples were collected from untreated embryos at corresponding times.
Methods: Ten to 20 headfolds were pooled in each microarray sample. Quality assurance of RNA was determined by absorbance A260/280 ratio of >1.7 and 1 percent agarose gel electrophoresis to rule out degradation of the isolated RNA and/or DNA contamination, and RT-PCR amplification of positive control genes (beta-actin, 16S ribosomal RNA). For each measurement, RNA (1 µg) was converted to cDNA with incorporation of biotin-11-dCTP or fluorescein-11-dCTP to generate 800-1,200 base labeled nucleotide fragments from the 3'-end of the transcript. Approximately 15-20 ng labeled target RNA was competitively hybridized (56-59oC) to glass microarray probes printed with sequence-verified cDNAs (PerkinElmer Life Sciences) representing more than 10 different tissue sources (80 percent brain-derived, more than 40 percent full-length cDNA sequence). Histochemical detection used peroxidase-conjugated anti-fluorescein and cyanine-3-tyramide (FL-Cy3) reaction followed by peroxidase-conjugated streptavidin and cyanine-5-tyramide (biotin-Cy5) reaction. Each array was repeated with reversal of the labeling scheme (e.g., FL-test / biotin-reference v. biotin-test / FL-reference) and independent sample replication. The log-intensity ratio of test/reference genes was normalized with locally weighted regression, standardized using absolute deviation of data points from their mean, and equalized between paired swaps for each gene in the matrix. Data reduction and visualization used GeneSpring software (v 4.1.5, Silicon Genetics, Redwood City, CA). Workflow was designed to select outlier genes that passed several quality control criteria, including an absolute signal threshold for minimal detection and concordance for the test/reference ratio between replica samples.
Figure 1. Data display shows the effect of normalization on the spread of 22 microarrays (log2 scale). Box and whiskers plot shows median (line) 25th-75th percentiles (box) and 5th-95th percentiles (whiskers) with outliers as individual points. Adjacent pairs are replica swaps. Left panel plots log2 ratio before normalization. Right panel plots these data after normalization. The alignment of median (zero) and scaling of data brought 90 percent of data points within 0.25. Thus, in each case the outliers are defined as those genes (~240) falling below or above the 5th-95th percentiles, respectively at a level of 0.5-fold.
Results: The experimental design tested embryonic exposure scenarios that have weak but discrete adverse effect on fetal morphology. Consistent with this design, the normalized microarray data set revealed weak effects of toxicant exposure on the profile of embryonic gene expression. Normalization of the microarray data imposed two major constraints on the individual samples (Figure 1). It centers the normal expression centered to a mean value of 0.0, and it scales the variance to an absolute deviation of 0.3. Assuming the differentially expressed genes on a genome-wide level are relatively few in number, this normalization method should be independent of outliers. Consistent with this assumption, outlier genes explained only part of the structure in the expression data. Some outlier elements (genes) did not adequately distinguish parameter-dependent variance from parameter-independent variance; furthermore, some non-outlier genes fit the data structure very well. A systematic approach was developed to enrich the data set for a core of genes having parameter-dependent variance. The approach was to conditionally filter outlier genes to anchor them to the basis vector for each experimental parameter and to decompose the resulting covariance matrix into principal components for each parameter. Hierarchical clustering of the data was applied to distinguish parameter-dependent expression variance from parameter-independent expression variance. Appositional expansion of this gene list used a minimum correlation coefficient of 0.95 to incorporate genes not included in the outlier domain. The resulting gene list numbered 331 comprised of 81 genes from the dose response, 118 from the time course, and 161 from PK11195 intervention. Graphical representation of the principal components of response for this core response signature is shown in Figure 2.
Figure 2. Core response signature of 331 genes influenced by low dose 2CdA exposure in day 8 mouse embryos. Left panel: 2CdA dose response analyzed in day 8 mouse embryos at 3.0 hours post-exposure. Two principal components of dose response appeared between 0.625 (nonteratogenic) and 5.0 mg/kg (teratogenic) exposures. Middle panel: 2CdA time course analyzed in the headfold after 2.5 mg/kg exposure. Three principal components were observed of response between 3.0 hours and 6.0 hours post-exposure. Right panel: 2CdA time course analyzed in the headfold between 3.0 hours and 6.0 hours post-exposure for embryos rescued with PK11195.
Genes associated with the dose vector showed two principal components of response that were diametrically opposed to one another. The two behaviors diverged with increasing dosage. Perhaps importantly, the point of divergence was approximate to the BMDL for microphthalmia at 2.0 mg/kg. This implies a novel cellular regulation invoked as the toxicant exposure crosses the exposure threshold for developmental toxicity. Genes associated with the time vector showed three principal components of response. Two behaviors were diametrically opposed to each other, and a third was intermediate between them. Temporal changes following 2CdA exposure should be interpreted within the context of p53 protein induction. This is a critical step in the pathogenesis of microphthalmia with 2CdA. Because p53 begins to accumulate in embryos between 3.0 hours and 4.5 hours after 2CdA exposure, all three temporal behaviors are invoked before p53 protein induction. Three principal components were observed of response for genes associated with the time vector during PK11195 rescue. Two behaviors were diametrically opposed to each other, and a third was intermediate. Changes invoked with 2CdA + PK11195 co-treatment, unlike the straight 2CdA time course, were more or less confined to 3.0 hours and 4.5 hours post-exposure times. This evidence is among the first to indicate that regulation of the embryonic transcriptome is directly relevant to disease susceptibility during the critical period following toxicant exposure.
Because changes in gene expression reflect the primary vulnerability of
embryonic systems to toxic insult as well as cellular adjustments that may or
may not contribute to the disease process, the 2CdA core response was narrowed
down to genes that comprehensively describe the covariance between experimental
parameters. The solution to this problem derived from basic elements of
symmetry. A Pearson covariance matrix was generated from the expression data of
core response genes and to reveal orthogonal symmetry. Decomposition of this
matrix was undertaken by serially deleting genes disturbing regularity (Figure
3). The orthogonal blocks of genes that correlate (blue) or anticorrelate
(yellow) with one another represent principal components of response. Assuming
genes arrive at their position due to parameter-dependent expression variance,
the regularity between parameters can be regarded as a metric for coherent
expression patterns and variation across conditions.
Figure 3. Covariance matrix of 162 genes in the 2CdA core response. The Pearson correlation between any two genes is represented by a colored square at their point of intersection. Strong correlations (+1.0, syn-) are shown by yellow color intensity and anticorrelations (-1.0, con-) as blue. These may be contrasted with the diagonal element, which represents the singular value with perfect correlation (same genes). Letters correspond to the data set used for comparison: (A) all three parameters of dose, time, and intervention; (B) 2CdA time course; (C) dose response; and (D) PK11195 intervention.
Several general characteristics are worth noting in the covariance display. One pertains to the regularity contributed by gene expression data derived from each experimental parameter. In Figure 3A, the covariance matrix is based on expression profiles of the core response for all three experimental parameters studied (dose, time, intervention). The pattern closely matches Figure 3B, which shows covariance as a function of the 2.5 mg/kg 2CdA time course. Patterns generated for covariance as a function of the two experimental parameters were incomplete. Data derived for the 2CdA dose response (Figure 3C) or PK11195 co-treatment time plot (Figure 3D) each matched to distinct domains of the overall structure, but not the structure in its entirety. It is postulated that the temporal response of the embryonic transcriptome best describes the principal components of the neodisease state. If this is true, then the temporal response signature of the precursor target cell populations to a low-dose exposure may be a better class predictor of developmental toxicity than the dose response signature at one particular time point. This has clear implications for risk assessment of developing systems, which are dynamic in nature, and hence for children's vulnerability to toxicant exposures. A second generalization that can be made from results of the covariance display pertains to the inter-relationships of gene clusters. The microarray data reduction scheme clustered genes into orthogonal sectors having positive (blue) and negative (yellow) correlation. Sector boundaries did not always carry between parameters. For example, anticorrelations in the dose sector carried as a function of time but not intervention.
Figure 4. Hierarchical clustering of the 2CdA core response. Dendritogram derived from replica microarray samples on day 8 of gestation (DOSE is mg/kg referenced to corresponding control sample; TIME is hour post-exposure referenced to control sample; INTRV designates 2CdA co-treatment with PK11195 referenced to corresponding sample exposed to 2.5 mg/kg 2CdA alone). Genes were color-coded for up regulation (red) or down regulation (green). Color intensity gives the magnitude relative to normal expression (black).
Hierarchical clustering revealed complex relativity within the genes used to
construct the covariance matrix across conditions. One issue worth noting is the
data display comparing earlier samples, which used RNA isolated from the whole
day 8 mouse embryo, versus the more recent samples that were enriched for
headfold and that used a second generation labeling and detection system (Figure
4). The lower cellular complexity and improved reagents increased the
sensitivity of the assay system as might well have been expected. Coherent
changes in expression ratios became evident in the embryonic headfold at a
low-dose level (2.5 mg/kg). Analysis of the PK11195 co-treated samples revealed
changes of a complex, indicating that this agent was not merely weakening the
responsive to 2CdA. These changes more likely reflect the novel downstream
action of PK11195. Consistent with this interpretation, PK11195 co-treatment
significantly rescued fetuses from microphthalmia up to 4.5 hours after 5.0
mg/kg 2CdA exposure (Figure 5). It is concluded that responses of the embryonic
transcriptome provide useful information toward understanding chemical mode of
action and perhaps the mechanisms underlying vulnerability of specific
developing systems to toxicant exposures.
Figure 5. The incidence of malformed fetuses during time delay of 4.0 mg/kg PK11195 given after 5.0 mg/kg 2CdA. Control (CON) and 5.0 mg/kg 2CdA alone (2CdA) groups are given for reference. PK11195 had no adverse effect alone. Distributions of malformation incidences are shown as the mean (line), 25th-75th percentiles (box), and 5th-95th percentiles (whiskers) for litters (n = 6 per group). ANOVA revealed significant differences (P = 0.020); post hoc analysis with Dunnett's Multiple Comparison Test localized group differences to 2CdA alone and 6.0 hours co- treatment groups relative tothe control (P < 0.05, asterisk).
Normal embryonic development is dependent on multiple genetic signals and responses that control cellular decisions to proliferate, differentiate, migrate, or die. Animal development uses approximately 17 distinct signaling pathways that integrate morphogenetic processes with the more fundamental regulatory processes that govern cellular growth and metabolism. Critical questions pertain to the relativity between gene expression clusters and coherent pathways as well as to how altered expression of coherent pathways can be used to predict chemical mode of action in developmental toxicity. When the genes associated with 2CdA's core response signature were annotated and classified by linkages to cell signaling pathways, there was striking representation of the tyrosine kinase (RTK) signaling pathway and its putative downstream targets. Early embryos rely heavily on RTK signaling to regulate cell growth, shape, and motility. Many birth defects, including microphthalmia, have been linked to altered function of the RTK pathway in humans and experimental animals. These features prompt a more thorough investigation of the RTK pathway in exposed embryos using bioinformatics and computational biology. RTK signals flow from plasma membrane receptor > adaptor protein > guanine nucleotide exchange factor > Ras (small G-protein) > kinase cascade > transcription factors. Signal diversification and integration with parallel signal transduction pathways is accomplished at multiple levels in this pathway. A variety of distinct growth factor receptors, splice variants, docking proteins, and kinases are known that ultimately determine the intracellular flow of signals. For simplicity, the main flow of signals associated with fibroblast growth factor (FGF) was used. FGF receptors FGFR1 and FGFR2 stood out by their differentially sensitive to the exposure conditions studied in this project. RT-PCR analysis confirmed the expression of these receptors and further suggested presence of at least three splice variants of FGFR1 in the mouse embryo. FGFR3 was detected but unaffected by exposure, and FGFR4 was not detected.
Ligand (FGF) binding to RTKs results in dimerization and activation of
protein tyrosine kinase activity. Dimers in turn crossphosphorylate each other
to initiate transient assembly of an intracellular signaling complex that relays
the signal into the cytoplasm. Two of the main intracellular signaling proteins
that recognize the phosphorylated receptors are adaptor protein (Grb2) and the
sevenless guanine nucleotide exchange factor (SOS). Two homologues of these
genes (Grb2, SOS2) were differentially expressed in response to 2CdA exposure.
The adaptor proteins couple RTK to Ras proteins, which are 21-KDa GTP-binding
proteins anchored to the inner face of the plasma membrane. Ras-p21 protein
cycles between inactive (RasGDP) and active (RasGTP) states. The former is
promoted by GTPase activating proteins (GAPs) and the latter by guanine
nucleotide exchange factors (GEFs). Several Ras-related genes were
differentially expressed following 2CdA exposure, including a GEF (KIAA0351) and
regulators of Ras-p21 activation (Krev-1, Rsu-1, rap1GAP). Expression patterns
of these genes is generally consistent with block regulation of elements from
FGFR1 > > Ras in this intracellular signaling cascade. Ras-GTP has several
potential downstream target kinases that couple the flow of signals to genes for
proteins that mediate reorganization of the actin cytoskeleton, intracellular
vesicle transport, and stress responses. The most widely represented downstream
target identified in the response to 2CdA pertains to the Rac-Rho pathway (RhoG,
activated p21cdc42Hs kinase, Rac3, LFP40, p122 RhoGAP, SRGAP3, N-WASP, Arp2,
HSN). Most of these genes showed syn-expression with the FGFR1 >> Ras
segment, a pattern generally consistent with genes promoting reorganization of
the actin cytoskeleton. Initial analysis identified genes involved in the Rab
(rab2, rab5, Rip-1, RanBP5) and Raf (Raf, p38/MAPK12, MAPKAP kinase) pathways
that couple RasGTP with intracellular vesicle transport and the MAP kinase
Figure 6. Impact of PK11195 on the 2CdA molecular benchmark signature corresponding to a 5 percent increased risk for microphthalmia (left panel) and intervention with PK11195 (right panel). Left panel: syn and con clusters of expression ratios (log2 2CdA/control ratio, normal expression = 0) for 96 genes in the day 8 embryonic headfold as a function of time after exposure to 2.5 mg/kg 2CdA. The "green cluster" included the major block of FGFR1 >> Ras >> Rac/Rho pathway genes. Ratios are plotted on log2 scale. Right panel: Same clusters for headfold samples co-treated with PK11195 (log2) expressed as a ratio of PK11195 + 2CdA samples relative to 2CdA alone (e.g., above and beyond the effect of 2CdA alone).
Cluster analysis grouped the RTK elements into distinct syn and con behaviors relative to the time metric following 2.5 mg/kg 2CdA exposure. Using the centroid for each cluster to model similarly expressed genes at a minimum correlation coefficient of 0.95 and across a variety of conditions, a molecular benchmark signature of 96 genes (Figure 6) was captured. The expanded gene set included many transcription factors and post-transcriptional regulatory elements that may link with RTK signals. Several interesting genes had syn-expression with the FGFR >> Ras >> Rac-Rho pathway, including FGFR1, amyloid protein A4 (695 and 751 amino acid splice variants), DiGeorge's critical region 2 (DCGR2), and Treacher-Collins syndrome (TCOF-1). These genes were transiently up regulated at 3.0 hours after 2CdA and further up regulated with PK11195 co-treatment during the post-exposure preceding irrevocable pathogenesis of the eye. Alteration at the critical time (4.5 hours) was not further influenced by PK11195, suggesting that the former changes are adaptive and the latter pathogenic. On the other hand, several interesting genes had con-expression with the FGFR >> Ras >> Rac-Rho pathway, including FGFR2, Pax6/aniridia, and retinitis pigmentosa-3 (RP3), which are all associated with ocular malformations in humans and mice. Down regulation at 3.0 hours after 2CdA continued with PK11195 co-treatment, and this again may reflect critical adaptive signals prior to the time of irrevocable pathogenesis.
Future activities for this project include:
(1) Run additional microarrays for headfold samples derived from embryos exposed to a low dose of 2CdA (2.5 mg/kg) and combined exposure to Ro5 4864, an agonist of the Bzrp. This will allow a comparison of the molecular benchmark signature for the antagonist (PK11195) and agonist (Ro5 4864) responses.
(2) Use bioinformatics and computational biology to continue identifying potential linkages.
(3) Conduct PCR validation of specific biomarkers to confirm theoretical linkages between coherent expression clusters across the different experimental conditions.