|Home | About | Journals | Submit | Contact Us | Français|
The goal of this study was to identify a set of hepatic genes regulated by ligand-dependent activation of the estrogen receptor in juvenile rainbow trout (Oncorhynchus mykiss). A custom rainbow trout oligo DNA microarray, which contains probes targeting approximately 1450 genes relevant to carcinogenesis, toxicology, endocrinology, and stress physiology was utilized to identify transcriptional fingerprints of in vivo dietary exposure to 17β-estradiol (E2), tamoxifen (TAM), estradiol + tamoxifen (E2 + TAM), diethylstilbestrol (DES), dehydroepiandrosterone (DHEA), dihydrotestosterone (DHT), and cortisol (CORT). Estrogen exposure altered the expression of up to 49 genes involved in reproduction, immune response, cell growth, transcriptional regulation, protein synthesis and modification, drug metabolism, redox regulation, and signal transduction. E2, DES, and DHEA regulated 18 genes in common, mostly those associated with vitellogenesis, cell proliferation, and signal transduction. Interestingly, DHEA uniquely regulated several complement component genes of importance to immune response. While the effect of TAM on E2-induced changes in gene expression was mostly antagonistic, TAM alone increased expression of VTG1 and other genes associated with egg development and immune response. Few genes responded to CORT treatment, and DHT significantly altered expression of only one gene targeted by the OSUrbt array. Hierarchical cluster and principal components analyses revealed distinct patterns of gene expression corresponding to estrogens and non-estrogens, though unique patterns could also be detected for individual chemicals. A set of estrogen-responsive genes has been identified that can serve as a biomarker of environmental exposure to xenoestrogens.
Estrogens are essential for the normal growth, differentiation, development, and function of reproductive tissues. Moreover, these steroid hormones are important regulators of various physiological processes in other tissues considered nonclassical estrogen targets including bone, liver, kidney, cardiovascular system, spleen, lung, and brain. Most biological actions of estrogens in humans and other vertebrates are mediated by two distinct, differentially expressed intracellular receptors, estrogen receptor alpha (ERα) and estrogen receptor beta (ERβ), which are considered unusually promiscuous receptors having affinity for numerous environmental contaminants and pharmaceuticals. Human and animal exposures to xenoestrogens have been linked to reproductive tract and development abnormalities and increased risk of breast, ovary, and liver cancer (Damstra et al., 2002; Flötotto et al., 2001; Martin et al., 2007). Simple, standardized in vitro assays are available for preliminary screening of chemicals for possible carcinogenic or endocrine-disrupting properties including various ER gene reporter assays (Harris et al., 1997; Sanseverino et al., 2005). However, with recent advances in genomic technology, sophisticated experiments are now possible in which DNA microarrays are utilized to determine the expression level of hundreds to thousands of genes simultaneously. The potential of “omics” technologies applied to the field of toxicology is considerable and includes drug development, risk assessment of novel drugs, novel applications for existing drugs, discovery of diagnostic markers, and compound screening for potential toxic, carcinogenic, and/or endocrine-disrupting effects (Benninghoff, 2007; Hamadeh et al., 2002; Shirai and Asamoto, 2003). Researchers anticipate that DNA microarrays may be routinely used in the identification of a set of differentially regulated genes, a gene expression fingerprint, that would serve as a biomarker of specific chemical exposure and could be used in predicting adverse health effects based upon the apparent mode of action discerned from the gene expression data (Hamadeh et al., 2002). Several studies have described distinct gene expression patterns resulting from exposure to chemicals from multiple compound classes including heavy metals, polycyclic aromatic hydrocarbons, and xenoestrogens (Bartosiewicz et al., 2001; Hamadeh et al., 2002; Koskinen et al., 2004; Terasaka et al., 2004; Vezina et al., 2004; Wang et al., 2004). The knowledge of the potential mechanism of action of a suspect chemical can be very useful in the evaluation of potential hazards that xenobiotics may pose to human and environmental health.
Rainbow trout have been widely used as a research model in environmental toxicology studies for decades, and this species has been proposed as an appropriate animal model for studies in comparative ecotoxicogenomics (Denslow et al., 2007). Previously, researchers described the design, manufacture, and application of a custom oligonucleotide microarray (OSUrbt array) for trout (Tilton et al., 2005). The OSUrbt version 2 microarray contains 1676 elements representing approximately 1450 genes important for carcinogenesis, environmental toxicology, comparative immunology, stress physiology, and endocrinology. In the present study, the OSUrbt microarray was used to identify a set of hepatic genes regulated by ligand-dependent activation of the estrogen receptor in juvenile rainbow trout. We hypothesized that (1) the transcriptional profile resulting from exposure to model estrogen compounds would be distinct from that of non-estrogens, (2) a subset of genes commonly regulated by all estrogen treatments could be identified, and (3) the selective estrogen receptor modulator (SERM) tamoxifen (TAM) would antagonize estrogen-induced changes in expression of these fingerprint genes.
Mt Shasta strain rainbow trout were hatched and reared at the Sinnhuber Aquatic Research Laboratory at Oregon State University in Corvallis, Oregon. Fish were maintained in flow-through 375-l tanks with activated carbon water filtration, water temperature at 12°C, and a 12:12 h light:dark cycle. Each experimental group consisted of three tanks containing 15 juvenile trout approximately 15 months old with an average weight of 140 g. Fish were fed a maintenance ration (2.8% of body weight) of Oregon Test Diet, a semipurified casein-based diet with menhaden oil as the lipid (Lee et al., 1991), for 2 weeks prior to experimental treatments. All procedures for treatment, handling, maintenance, and euthanasia of trout used in this study were approved by the Oregon State University Institutional Animal Care and Use Committee.
Eight experimental treatments were selected to examine hepatic gene expression in response to estrogens, a SERM, and non-estrogens. Treatments and concentrations (calculated with respect to diet wet weight) are as follows: control (CON), 5 ppm 17β-estradiol (E2), 50 ppm TAM, 5 ppm E2 + 50 ppm TAM (E2 + TAM), 2 ppm diethylstilbestrol (DES), 500 ppm dehydroepiandrosterone (DHEA), 5 ppm dihydrotestosterone (DHT), and 5 ppm cortisol (CORT). Test chemicals were purchased from Sigma-Aldrich (St Louis, MO) and were added directly to the oil portion of the custom diet. Experimental treatments were administered for 14 days, and feeding occurred 5 days per week. The concentrations of E2 and DHEA were selected based upon their ability to significantly induce vitellogenin (Vtg) gene and protein expression (Tilton et al., 2006), and the concentration of DES was chosen based upon reported relative binding affinity for the trout estrogen receptor (Matthews et al., 2000). To examine potential agonistic and antagonistic effects of TAM, a 10-fold higher concentration of that chemical compared to E2 was tested alone and coadministered with 5 ppm E2. Finally, DHT and CORT were administered at concentrations equivalent to that of E2.
On day 15, trout were euthanized with an overdose (250 ppm) of tricane methanesulfonate. Four males were sampled from every experimental tank. Immediately after sacrifice, an approximate 200 mg portion of each liver was quick frozen in TRIzol reagent (Invitrogen, Carlsbad, CA) and stored at − 80°C. Body weight and liver weight were recorded at the time of sacrifice to calculate the hepatosomatic index (HSI = liver weight/body weight × 100), and blood sera E2 levels were measured using a commercial enzyme immunoassay kit (Cayman Chemical, Ann Arobor, MI).
Total hepatic RNA was extracted in TRIzol reagent according to supplier’s instructions. Equal amounts of RNA from each individual liver sample were subpooled according to tank (four males pooled per tank) resulting in three biological replicates for each treatment. A reference RNA pool was made by combining equal amounts of RNA from all individual male CON liver samples (12 liver samples in total). Following cleanup with the RNeasey Mini kit (Qiagen, Valencia, CA), RNA quantity and quality were ascertained using the NanoDrop ND-1000 UV-Vis Spectrophotometer (NanoDrop Technologies, Inc., Wilmington, DE) and the Bioanalyzer 2100 (Agilent, Palo Alto, CA), respectively.
Details on the development, manufacture, and quality control assessment of the OSUrbt microarray have been provided previously (Tilton et al., 2005) and will be only summarized here. Array printing and quality control analysis were conducted at the Center for Genome Research and Biocomputing at Oregon State University. Each element was printed in duplicate onto Corning UltraGap slides (Acton, MA). Additionally, 16 replicate sets (one for each array block) of 10 SpotReport Alien Oligos (Stratagene, La Jolla, CA) were also printed on each array. Inclusion of these control features allows for standardization across arrays and estimation of cDNA labeling efficiency, hybridization sensitivity, and intra- and interarray variability. Buffer-only spots were included as negative controls. Altogether, each array consists of 4096 spots. Array print quality (shape and alignment of features, missing features, etc.) was assessed by scanning for red reflectance using a ScanArray 4000 scanner (PerkinElmer, Boston, MA). One array per printing batch was subjected to quality control hybridization (SpotQC, Integrated DNA Technologies) to confirm absence or presence of DNA in appropriate features. Arrays were stored for less than 6 months prior to hybridization.
For detection of gene expression on the OSUrbt array, the Genisphere 3DNA Array 350 kit (Hatfield, PA) was used in a dye-swap, reference sample design following the supplier’s protocol. During reverse transcription (SuperScript II, Invitrogen) of total RNA (7 μg) using the supplied Genisphere oligo d(T) primers, each biological sample was tagged separately with a capture sequence for one of two fluorescent dendromer reagents, Cy3 or Cy5. Each reverse transcription reaction also included spiked-in mRNA corresponding to the SpotReport Alien Oligo controls. Corresponding reference RNA was reverse transcribed with the capture sequence for the opposite reagent (Cy3 or Cy5). To account for any differences in dye labeling efficiency, technical replicates were performed for two of the three biological replicates so that the capture sequence tags for Cy3 and Cy5 were swapped between the sample and reference RNA (replicate A, reference is labeled Cy3 and sample is labeled Cy5; replicate B, reference is labeled Cy5 and sample is labeled Cy3). Five arrays per experimental treatment were processed, including the CON treatment.
Following reverse transcription, equal amounts of capture sequence–tagged sample and reference cDNA were combined, concentrated using the Millipore Microcon YM-30 (Billerica, MA), and treated with RNase cocktail (Ambion, Austin, TX) to degrade any remaining RNA. Hybridization of cDNA and capture reagents to the OSUrbt arrays was performed as described previously (Tilton et al., 2005) with a few modifications. Sterile 20× standard saline citrate (SSC) was prepared with 3 M sodium chloride and 0.3 M sodium citrate at a pH of 7.0; appropriate dilutions were made with sterile water for hybridization wash buffers. Sodium dodecyl sulfate (SDS) was also added to some hybridization wash buffers at a 0.1% (wt/vol) final concentration. Prior to hybridization, arrays were prewashed with gentle mixing on a rocking platform twice in 0.1% SDS for 5 min at room temperature (RT), twice in 2× SSC, 0.1% SDS for 10 min at 47°C, once in 0.1× SSC at RT for 5 min, and once in ultra-pure water for 5 min. Slides were dried by centrifugation. cDNA was hybridized to the arrays in formamide-based hybridization buffer (25% formamide, 4× SSC, 0.5% SDS and 2× Denhardt’s solution, Genisphere, Inc.) for 16 h at 47°C in a Hybex Microarray Incubation System (SciGene Corporation, Sunnyvale, CA). After hybridization, the arrays were rinsed in 2× SSC, 0.1% SDS for 10 min at 47°C in the Hybex water bath. Then, arrays were washed with gentle mixing twice in 2× SSC, 0.1% SDS for 5 min at 47°C, twice in 1× SSC for 5 min at RT, and twice in 0.1× SSC for 5 min at RT. Slides were dried by centrifugation and prewarmed to 49°C. Next, the 3DNA Cy3/Cy5 capture reagent mixture was hybridized to the arrays for 3 h at 49°C in the dark. After this second hybridization, the arrays were rinsed, washed, and dried again as previously described in SSC buffers containing 100 μM DTT dithiothreitol in the dark. All arrays were treated with DyeSaver2 Anti-fade Coating (Genisphere, Inc.) to preserve fluorescence of the capture reagent.
Within 24 h of hybridization, array images at a resolution of 5 μm were obtained using the ScanArray 4000 (PerkinElmer) with 90% laser power at 543 and 633 nm excitation wavelengths for Cy3 and Cy5, respectively. Photomultiplier tube settings were adjusted using the intensities of the SpotReport Alien Oligo controls to balance global intensities for the Cy3 and Cy5 fluors and to normalize across slides.
Array image files were processed, and spot intensities were quantified using QuantArray software (PerkinElmer); raw median signal and background values were obtained and uploaded to BioArray Software Environment (Saal et al., 2002) for database management and data normalization. For each array feature, raw median signal values were background corrected, LOWESS normalized, and exported to Excel (Microsoft Corp., Redmond, WA). Fold change in gene expression was determined by calculating the ratio of fluorescence intensity for the sample compared to the reference for each array feature. Average fold change values were calculated by determining the geometric mean for array technical replicates (duplicate spots for each element per array), treatment technical replicates (dye swaps), and finally for the biological replicates (n = 3) so that triplicate fold change expression values were available for every array element; these values were log2 transformed and uploaded to TM4 MultiExperiment Viewer (TMeV) (Saeed et al., 2003). Statistical analyses of gene expression were performed in TMeV using the normalized, geometric mean expression values for each biological replicate to compare experimental treatments to the control; a statistically significant change in gene expression was inferred when p < 0.05 (Welch’s t-test, between subjects and assuming unequal variances). Unsupervised, bidirectional hierarchical clustering and principal components analyses were performed in TMeV. Normalized data were also exported to GraphPad Prism 4 (GraphPad Software, Inc., San Diego, CA) for pairwise correlation analysis and for statistical comparisons of quantitative PCR expression data.
Array feature annotation was performed by querying the Dana-Farber Cancer Institute (DFCI) R.trout Gene Index (http://compbio.dfci.-harvard.edu/tgi/) for the closest EST match to the array 70-mer sequence. Matching EST sequences were then BLASTX queried in the National Center for Biotechnology Information genome database. The top hit (lowest E-score) was selected as the matching gene. If an EST had no significant (E-value < 10−6) BLASTX hit, then the most significant BLASTN hit is shown. Additional information on the matching gene and BLASTX hit information are provided in Supplementary Table 3 including the DFCI accession number, the accession number for the corresponding GenBank nucleotide sequence and the species, E-value, % identity, and SwissProt accession number for the most significant BLASTX hit. For the proteins encoded by the putative trout homolog mRNAs, functional information was inferred from annotations in the Gene Ontology, Online Mendelian Inheritance in Man (OMIM), and SwissProt Protein Knowledgebase databases.
To validate changes in gene expression detected on the OSUrbt array, mRNA levels of select genes were evaluated by the quantitative real-time reverse transcriptase (qRT) PCR. Total hepatic RNA was used from the previous isolation as described above and treated with DNase (Invitrogen) following the supplier’s directions. Total RNA (1 μg) was reverse transcribed (Superscript II, Invitrogen) according to the supplier’s protocol with oligo d(T)18 primer and a final reaction volume of 50 μl. PCR was conducted using 1 μl of synthesized cDNA template, 0.3 μM each primer, and 1× SYBR Green PCR master mix (DyNAmo qPCR kit, Finnzymes, Finland) with a final reaction volume of 20 μl. Primer pairs were designed so that the PCR product contained the original 70-mer oligonucleotide sequence (unless noted otherwise), and the primer sequences are shown in Table 1. PCR was performed with a DNA Engine Cycler and Opticon 2 Detector (MJ Research, Waltham, MA) for 35 cycles with denaturation at 94°C for 10 s, annealing at the optimum temperature determined for each primer pair (54°C – 62°C) for 20 s and extension at 72°C for 18 s. For DNA quantification, the C(t) value was set to ensure that values were within the linear range of the fluorescence curve. PCR standards for each gene were prepared by gel purification of PCR products (QIAX II, Qiagen, Valencia, CA), quantified using the PicoGreen dsDNA Quantification Kit (Molecular Probes, Eugene, OR), and serially diluted for final concentrations ranging from 0.001 to 100 ng DNA.
Glyceraldehyde-3-phosphate dehydrogenase (GADPH) and β-actin were selected as genes for normalizing the qRT-PCR data since neither microarray nor qRT-PCR analyses revealed significant treatment-induced changes in mRNA levels of these two housekeeping genes. All expression values obtained by qRT-PCR were normalized by the mean fold change of both genes. These normalization values were calculated from the ratio of pg cDNA for a biological replicate over the average pg cDNA for all three CON biological replicates for both β-actin and GADPH genes as follows: HSKGA = (ActA/Actavg) × (GADPHA/GADPHavg) where HSKGA indicates the mean fold change of the two housekeeping genes, A indicates the biological replicate, and avg indicates the average of three biological replicates. For comparison to microarray expression values, fold-change ratios were calculated for treated samples compared to the CON treatment.
Blood samples from each donor animal were obtained at sacrifice. The plasma fraction of each sample was obtained by centrifugation at 850X g and then snap frozen at − 80°C. Blood plasma Vtg was measured following an ELISA procedure previously described (Donohoe and Curtis, 1996; Shilling et al., 2001) using a rabbit anti-chum salmon Vtg antibody graciously donated by A. Hara at Hokkaido University. Briefly, samples were incubated with the VTG antibody (1:1500) in 96-microwell plates for 24 h at 4°C and then were transferred to preblocked plates (1% bovine serum albumin) coated with purified rainbow trout VTG (25 ng per well) and incubated again for 24 h at 4°C. To develop the assay, the plates were first incubated with biotinylated donkey anti-rabbit IgG (1:1500, GE Healthcare, Buckinghamshire, U.K.) for 2 h at 37°C, then with streptavidin horseradish peroxidase conjugate (1:600, GE Healthcare) for 2 h at 37°C, and finally with 0.01% 3,3′5,5′-tetramethylbenzidine and 0.01% hydrogen peroxide in 0.5 M sodium acetate, pH 6.0, for 10 min at RT. Color development was stopped with 2 M sulfuric acid. Absorbance at 450 nm was measured on a SpectraMax 190 spectrophotometer with SoftMax Pro 4.0 software (Molecular Devices, Sunnyvale, CA). A Vtg standard curve (6.25–3200 ng/ml) was included in every assay plate, and the Vtg concentration of each sample was calculated from the linear portion of this curve (GraphPad Software). The interassay coefficient of variance was 8.4% with an assay detection limit of 6.25 ng/ml. Equal volumes (100 μl) of blood plasma from four males were pooled by tank so that each pool represented a biological replicate, as in the microarray experiment. Each sample was assayed in duplicate with a minimum of three dilutions (1:10 to 1:500 for samples with expected low Vtg levels; 1:5000 to 1:50 000 for samples with expected high Vtg levels) to ensure that absorbance values of unknowns were within the linear part of the standard curve.
The 14-day dietary exposure to the experimental treatments described above did not significantly alter body weight, although hepatosomatic indices were significantly increased in trout fed E2, E2 + TAM, and DHEA compared to control animals (Supplementary Fig. 1); HSI values were not significantly altered in animals treated with DES, TAM, DHT, or CORT. Blood sera levels of E2 in animals fed E2, E2 + TAM, or DHEA were also markedly elevated (Supplementary Fig. 2); DHEA did not significantly cross-react with the E2-specific antibody in the Enzyme immunoassay (data not shown).
All raw gene expression data are available at the Gene Expression Omnibus online microarray data repository (http://www.ncbi.nlm.nih.gov/geo/, accession GSE8226). Calculated fold-change ratios of experimental samples compared to reference after background subtraction and LOWESS normalization are provided in Supplementary Table 1. Quality control analysis of array hybridization using the SpotReport Alien oligos showed that, in general, intra- and interarray variability was low and that hybridization was consistent and reproducible (Supplementary Fig. 3).
Multiple criteria were used to reduce the original raw data sets to subsets of genes considered as significantly regulated by any one of the experimental treatments. Results of the stepwise application of these criteria are shown in Table 2. These criteria included filters for statistical significance (p value < 0.05 by Welch’s t-test with unequal variance), spot consistency (> 1.5-fold change in 9 out of the 10 biological and technical replicate features) and mean fold change (> 2-fold). Many more genes were identified as statistically significant by the Welch’s t-test than were included in the final gene list (Table 2). Inclusion of the 2-fold change criterion eliminated numerous genes that were modestly induced or repressed by the experimental treatments.
Genes that passed all three stringency criteria are listed in three tables categorized by the type of response or lack thereof to E2 exposure. Genes induced or repressed by E2 are listed in Tables Tables33 and and4,4, respectively. Genes that were not altered by E2 exposure, but were induced or repressed by one of the other experimental treatments, are listed in Table 5. Gene descriptions are provided based on sequence homology using the most significant (E < 10−6) BLASTX hit. While most OSUrbt array features represent unique genes, some elements represent distinct oligonucleotide sequences for the same gene (i.e., OmyOSU208 and OmyOSU222 array features both targeted vitellogenin). Thus, the average of the expression values for all features matching the indicated gene is reported (Tables (Tables33–5); mean fold change values for each array element are available in Supplementary Table 3.
Of the approximate 1450 genes represented on the OSUrbt array, fewer than 4% were significantly regulated by any of the experimental treatments. E2 exposure altered the expression of 49 genes, more than any of the other compounds tested. More transcripts were significantly induced by E2 (38 genes) than were repressed (11 genes), and this pattern was similar for most of the experimental treatments (Tables (Tables33 and and4).4). E2, DES, and DHEA commonly altered expression of 18 genes (Fig. 1), though DHEA exposure differentially regulated 18 additional genes that were not similarly affected by E2 or DES (Table 5). TAM induced expression of five transcripts that were also upregulated by E2 (Fig. 1). When TAM was coadministered with E2, fewer total genes were differentially regulated compared to E2 alone. DHT, a nonaromatizable androgen, significantly altered the expression of only one transcript, which was also induced by E2. Exposure to the non-estrogenic corticosteroid hormone CORT had no effect on expression of genes that were regulated by E2, though expression levels of three other mRNAs were uniquely induced (Table 5). Some genes were identified that surpassed both the spot consistency and fold change criteria but did not pass the statistical criterion due to high variability among the biological replicates. These genes have been included in Tables Tables33–5 for informational purposes but were not included in the treatment group tallies.
Genes have been categorized by function based on putative trout homologs using the Gene Ontology, OMIM, and SwissProt Protein Knowledgebase databases. The functional categories listed in Tables Tables33–5 were defined by the authors based on information obtained for each gene from the above sources. The biological processes associated with hepatic genes differentially regulated by dietary exposure to estrogenic and non-estrogenic treatments are indicated within each gene list table. The transcripts differentially regulated by E2 represent numerous biological processes important for signal transduction, transcription, translation, protein modification and protein transport, redox regulation, drug metabolism, stress and immune response, extracellular matrix, and cell cycle, growth, and proliferation. Genes associated with vitellogenesis and immune response were most strongly induced by estrogen exposure (> 30-fold increase), whereas the strongest repressions of gene expression (> 4-fold decrease) were observed for transcripts associated with cell cycle, growth, and proliferation (e.g., cdc2l1), transcription (e.g., dmrt2), and drug metabolism (CYP2K5). Finally, several transcripts encoding proteins of unknown biological function were differentially regulated in E2-exposed animals, including RTN9A-1 and RTN9A-2.
Of the 18 mRNAs commonly regulated by E2, DES, and DHEA, some encoded transcripts for proteins associated with vitellogenesis and egg development, signal transduction and transcription, protein transport, immune response, and cell proliferation. While 22 genes were determined to be uniquely regulated by E2, treatment with DES and DHEA influenced expression of most of these genes in a similar manner though to a lesser degree (e.g., klf9, dmrt5, tcpbp, chmtxn, sdha, and TM4SF). Some unique targets for E2 included a number of transcripts involved in signal transduction, drug metabolism, and protein synthesis, folding, and transport. DHEA exposure exclusively altered expression of 18 genes, including mRNAs encoding proteins associated with ion binding and transport (icta and fth2) and redox regulation (cox6A1). Also, multiple transcripts encoding complement components C3–3, C3–4, and C5 were significantly repressed in DHEA-exposed animals.
Genes induced by TAM were similar to those upregulated by E2, including mRNAs encoding proteins involved in vitellogenesis, transcription, and immune response, although the degree of response was much lower in TAM-exposed animals. TAM uniquely repressed expression of PGDS by approximately 5-fold compared to control. Combined treatment of E2 + TAM resulted in many fewer differentially expressed genes compared to E2 treatment alone; for most genes determined to be significantly induced or repressed by cotreatment of E2 + TAM, the degree of change was reduced compared to E2 exposure alone.
Transcripts induced by CORT exposure include hp1 and hp2, involved in the acute phase immune response, and gth2b, which is an important reproductive hormone.
Several analytical approaches were used to examine the relationships of transcript profiles among the estrogen and non-estrogen treatments. First, values for all 1676 array elements were charted in scatter-plot graphs to compare directly the E2-induced gene expression profile with each of the other experimental treatments (Fig. 2). A correlation coefficient (r value) for each comparison was calculated from the linear regression of the two data sets. Pairwise analyses revealed strong correlations between E2 and DES, DHEA, and E2 + TAM treatments (r = 0.891, 0.803, and 0.844, respectively) and a moderate correlation between E2 and TAM treatments (r = 0.687) (Figs. 2A–D). Alternatively, the E2 expression profile was poorly correlated with profiles for DHT and CORT treatments (r < 0.001 and = 0.085, respectively) (Figs. 2E and 2F).
Hierarchical cluster analyses were used to visualize patterns of gene expression for estrogen and non-estrogen treatments for either all array elements (Fig. 3A) or for a subset of genes that were significantly regulated by any one of the experimental treatments (Fig. 3B). Bidirectional hierarchical clustering, in which expression values were sorted according to similarity of expression among all genes (left tree) and among all samples (top tree), revealed distinct patterns of expression corresponding to two primary nodes in the sample tree. In a subset of genes that were regulated by at least one of the experimental treatments, the left node includes estrogen treatments E2, DES, and DHEA, as well as the E2 + TAM treatment group (Fig. 3B). The right node includes CON, DHT, and CORT samples as well as TAM samples within a separate subnode. Distinct patterns of gene expression are evident for each experimental condition.
Principal components analysis was used to reduce the dimensionality of the gene expression data sets so that general relationships between experimental treatments could be discerned more easily (Fig. 4). This analysis supports the observations from previous correlation analyses indicating that estrogen treatments have gene expression profiles distinct from non-estrogens. CON, CORT, and DHT treatments grouped closely together (Fig. 4, hatched gray circle), whereas E2, DES, and E2 + TAM treatments were grouped together distant from and were the non-estrogen group (Fig. 4, solid gray circle). TAM and DHEA treatments were isolated within the PCA plot, indicating that gene expression signatures for these compounds are distinct from each other and from the estrogen and non-estrogen groups.
Ten genes were selected for validation by qRT-PCR to represent several observed patterns of gene expression: highly induced by estrogens (VTG1, ESR2), moderately induced by estrogens (cstD and ESR1), repressed by estrogens (tcpbp and CYP2K5), repressed by TAM (PGDS), repressed by DHEA (C5), induced by CORT (hp1), and unchanged (CYP1A). In general, patterns of gene expression evaluated by qRT-PCR were highly similar to expression values obtained from the OSUrbt array (Fig. 5). For some genes, qRT-PCR detected a greater magnitude change in gene expression, most notably for VTG1. One notable discrepancy between array and qRT-PCR results was observed for ESR2. PCR analysis showed no significant treatment-induced changes in ESR2 expression, contrary to measurements with the array. Otherwise, no major discrepancies between the qPR-PCR and array data were noted for the remaining genes evaluated by PCR.
Levels of circulating Vtg protein detected in blood plasma were strongly correlated with hepatic Vtg mRNA expression detected by the microarray (Fig. 6). As was observed for Vtg transcript levels, TAM treatment significantly increased plasma Vtg protein, though not to the extent of the other estrogen treatments E2, DES, and DHEA. Moreover, TAM cotreatment did not significantly reduce E2-induced elevation of plasma Vtg protein, similar to results obtained for Vtg transcript levels with the microarray and qRT-PCR. DHT and CORT treatment did not significantly alter blood plasma Vtg protein levels.
In the present study, we identified a transcriptional profile representing in vivo exposure to estrogen in the liver of male juvenile rainbow trout. Distinct patterns of gene expression were evident for estrogen and non-estrogen treatment groups as well as for each specific chemical treatment. Eighteen genes were identified as robust and specific indicators of estrogen exposure. This gene set included known targets for xenoestrogens involved in oocyte development and reproduction (VTG1, ZRP, VEP) and transcription (ESR1, ESR2) as well as some lesser known or novel estrogen-regulated transcripts including genes involved in stress response (VHSV4, VHSV6), signal transduction (ikk1), protein transport (Sec61αB), cell growth and proliferation (cd82, BMP7, BRMS1L, and Cdc2l1), nucleoside metabolism (UrdPase 1), nitrogen metabolism (BTD), and two genes with unknown function (RTN9-A1, RTN9-A2). Given the rigorous criteria used to select genes as differentially regulated by experimental treatment, it is likely that this biomarker set is quite conservative. Though of potential biological interest, moderately responsive genes that did not meet the 2-fold arbitrary limit established for this study would likely not be robust indicators of xenoestrogen exposure.
A number of other studies have employed macro- or microarray platforms to investigate the hepatic transcriptional response to estrogen exposure in several animal models including various fish species such as rainbow trout (Hook et al., 2006; Tilton et al., 2006), European flounder (Williams et al., 2007), common carp (Moens et al., 2006), largemouth bass (Larkin et al., 2002) and sheepshead minnow (Larkin et al., 2003), and rodents (Boverhof et al., 2004; Kato et al., 2004). For studies employing fish animal models, VTG is invariably the most responsive gene to estrogen exposure, confirming the suitability of this gene as a sensitive biomarker of estrogen exposure (Sumpter and Jobling, 1995). Other transcripts common to most estrogen-regulated gene lists include mRNAs encoding other proteins involved in oocyte development, collectively referred to as the zona radiata proteins (also termed choriogenin and vitelline envelope proteins) (Arukwe and Goksøyr, 2003). However, upon careful inspection of these published gene lists, it is clear that a majority of the estrogen-regulated genes are different among these various studies. These differences may be related to the specific chemicals, species, exposure routes, or time points studied. Moreover, the state of genome annotation for each species studied as well as the diversity of array platforms and sample analysis techniques used could be contributing factors to discrepancies observed among the studies. Nonetheless, many of the biological processes involved in the transcriptional response to estrogen are very similar among these experiments. As one example, several studies report estrogen-induced changes in genes involved in lipid metabolism and/or transport including fatty acid–binding protein (this study; Hook et al., 2006; Tilton et al., 2006; Williams et al., 2007) and apolipoprotein B (Boverhof et al., 2004; Hook et al., 2006; Larkin et al., 2002; Moens et al.,2006; Tilton et al., 2006). Although the liver is not considered a classical estrogen-responsive tissue in rodents, exposure of ovariectomized mice to the pharmaceutical xenoestrogen ethynylestradiol elicited changes in transcript levels of numerous genes belonging to functional groups similar to those identified in fish studies, notably cell growth and proliferation, cytoskeleton and extracellular matrix, oxidative metabolism, and lipid metabolism and transport (Boverhof et al., 2004). Although there appears to be some conservation of the biological response to estrogen among these diverse organisms, the applicability of a specific estrogen transcriptional fingerprint across species is uncertain.
This study is the first to utilize a toxicogenomics approach to examine the hepatic transcriptional response to a combined exposure of E2 and TAM. The first clinically relevant SERM, TAM is a multifunctional drug with different actions at various target tissues including agonist action in bone and liver, antagonist action in the breast and brain, and pleiotropic activity in other tissues such as the endometrium (reviewed in Diel, 2002). The selectivity of SERMs may occur via multiple mechanisms, including differential binding of ER subtypes, recruitment of different coactivators or corepressors, and regulation of upstream promoter regions of target genes (Diel, 2002; Dutertre and Smith, 2000; Shang and Brown, 2002). This diversity of molecular mechanism of action can lead to a broad range of transcriptional responses. In the present study, E2-induced changes in expression of the 18 identified estrogen-responsive genes were decreased in animals co-treated with TAM, indicating that the observed transcriptional changes likely involved interaction with the hepatic ER. It is important to note that, in this study with a test concentration of TAM 10-fold greater than that of E2, only a partial antagonism of the E2-induced transcriptional profile was observed. A complete block of the E2-induced response may have been observed if a higher concentration had been used or if TAM had been administered prior to the E2 treatment. Indeed, Marlatt et al. (2006) observed that a 250-fold excess concentration of TAM was required to block completely the Vtg response to E2 in juvenile male trout. In the present study, we also observed that in vivo exposure to TAM alone markedly increased expression of several known estrogen target genes, including VTG1, VEP, and VHSV4. This observation of dual, contrasting transcriptional activities is consistent with results of a previous study in trout in which TAM elevated hepatic expression of VTG and ERα, yet antagonized E2-induced expression of these genes (Vetillard and Bailhache, 2006). However, plasma Vtg was not altered in juvenile trout that were exposed for 20 days to 500 ppm waterborne TAM (Marlatt et al., 2006); this difference could be related to the length or route of exposure. One potential mechanism for dual transcriptional activities of TAM in trout may be via interaction with different ER subtypes, of which there are two highly expressed in trout liver (ERα1 and ERβ2) (Nagler et al., 2007).
Although DHEA has been reported to have weak affinity for the trout ER (Matthews et al., 2000), the estrogenic response to DHEA was most likely due to its metabolism to estradiol, rather than direct interaction of DHEA with the trout hepatic ER. Average blood plasma E2 concentration in trout fed DHEA was more than 200-fold higher than in control animals and was nearly equivalent to plasma E2 levels in trout administered E2 directly. However, the gene expression profile of DHEA exposure was quite distinct from that of E2, indicating that DHEA may have effects on hepatic transcription independent of its E2 metabolite. Eighteen transcripts were identified as uniquely regulated by this compound, including genes involved in protein modification (sae2), cell growth (ccnd1), metabolism (dlat), platelet aggregation (cd9), ion binding and transport (icta and fth2), and others. In particular, DHEA treatment markedly repressed expression of several genes belonging to the complement pathway including C3–3, C3–4, C5, and to a lesser extent c1r/c1s. The complement system plays a critical role in innate immunity, particularly in pathogen defense (Sunyer et al., 2005). Perturbations of the complement system could lead to increased risk of infection and has been associated with some diseases involving the immune system, such as asthma (Wills-Karp and Koehl, 2005), lupus (Bao and Quigg, 2007), and HIV (Datta and Rappaport, 2006).
Although the concentration of DES was selected based upon its reported higher affinity for the trout ER compared to E2 (Matthews et al., 2000), results of this study suggest that the hepatic transcriptional response to DES was generally lower in comparison to E2. Tissue levels of the test chemicals were not measured; thus, it is possible that the reduced response to DES could be related to the internal dose received. Another notable observation of the present study was the general lack of a hepatic transcriptional response to DHT. Blum et al. (2004) observed changes in expression of 11 hepatic genes in largemouth bass exposed to a 10-fold higher concentration of DHT than was administered in the present study. Although the DHT dose chosen for this trout study was selected to match that of E2, this concentration may have been insufficient to induce a robust transcriptional response in our animal model.
In the present study, qRT-PCR was used to quantify mRNA expression of genes selected to represent patterns of estrogen induction, estrogen repression, and lack of response to estrogen. In general, qRT-PCR results matched well with the microarray data. However, for some genes such as VTG1, ctsD, and hp1, the extent of induction or repression was under-estimated as measured by the OSUrbt array. This phenomenon has been observed in a number of microarray studies (Boverhof et al., 2004; Hook et al., 2006) and is likely a consequence of background and saturation signal intensity limitations. Thus, microarray experiments are most useful for semiquantitative assessment of gene expression.
One noticeable inconsistency in the validation of our microarray results by PCR was for the ESR2 gene. The 70-mer oligonucleotide sequence for this array feature aligned very near the 3′ end of the sequence for ESR2 (accession AJ289883). Because a useful primer pair surrounding this aligned sequence could not be designed, a previously published primer pair targeting a different region of this gene was used for validation. Thus, the qRT-PCR assay utilized in this study did not truly mimic the OSUrbt array measurement of ESR2 gene expression. Moreover, it is possible that DNA for similar genes, such as one of the other ER isoforms (Nagler et al., 2007) or a splice variant of the rtERα (Pakdel et al., 2000), significantly cross-hybridized to this array feature. Hook et al. (2006) also had difficulty reproducing microarray expression data for rainbow trout ER by qRT-PCR, which they attributed to high affinity of the array sequence to a different ER isoform. The lack of perfect correlation between microarray and qRT-PCR methodologies could present a problem for experiments designed to identify novel genes in biological processes or for hypothesis generation. However, when employed in risk assessment or environmental monitoring studies, the overall pattern of gene expression corresponding to a particular chemical class is of greater importance than the specific identity of each differentially regulated feature.
The transcriptional profile identified for estrogen exposure in the present study can serve as a biomarker of environmental exposure to xenoestrogens. The identification of 18 genes commonly regulated by estrogens that represent a number of biological functions demonstrates that a toxicogenomics approach to chemical evaluation can provide more information about the potential mechanism of action of an unknown chemical exposure compared to traditional single endpoint biomarkers. While a complete genome array may be preferred in experiments designed for hypothesis generation, this economical, targeted microarray allows for large-scale studies utilizing tens to hundreds of arrays and increases the practical sample size of treatment groups allowing for more statistically robust experiments. In future research, this targeted array and the transcription profile of estrogen exposure obtained in this study will be used to categorize chemicals with unknown mode of action and to address the issue of predictive utility of transcriptional fingerprints in animals with unknown exposures.
The authors wish to acknowledge the assistance of Eric Johnson, Greg Gonnerman, and Cari Buchner at the Sinnhuber Aquatic Research Laboratory for care of the animals used in this study. We also thank Caprice Rosato and the Center for Genome Research and Biocomputing for assistance in the manufacture of the OSUrbt custom microarray. Technical assistance provided by Dr. Susan Tilton, Marilyn Henderson, and other members of the Williams laboratory is also gratefully acknowledged.
FUNDING National Institute of Environmental Health and Safety (ES03850, ES07060, ES00210, and ES013534).
Parts of this manuscript were presented at 14th North America International Society for the Study of Xenobiotics Meeting, Rio Grande, Puerto Rico, 22–26 October 2006.