Identification of estrogen-responsive genes is an essential step toward understanding mechanisms of estrogen action during mammary gland development. To identify these genes, 16 prepubertal heifers were used in a 2 × 2 factorial experiment, with ovarian status (intact or ovariectomized) as the first factor and estrogen treatment as the second (control or estradiol). Heifers were ovariectomized at ∼4.5 mo of age, and estrogen treatments were initiated 1 mo later. After 3 days of treatment, gene expression was analyzed in the parenchyma and fat pad of the bovine mammary gland using a high-density oligonucleotide microarray. Oligonucelotide probes represented 40,808 tentative consensus sequences from TIGR Bos taurus Gene Index and 4,575 singleton expressed sequence tags derived from libraries of pooled mammary gland and gut tissues. Microarray data were analyzed by use of the SAS mixed procedure, with an experiment-wide permutation-based significance level of P < 0.1. Considerable differences in basal gene expression were noted between mammary parenchyma and fat pad. A total of 124 estrogen-responsive genes were identified, with most responding only in the parenchyma or the fat pad. The majority of genes identified were not previously reported to be estrogen responsive. These undoubtedly include genes that are regulated indirectly but also include known estrogen-targeted genes and novel genes with potential estrogen-responsive elements in their promoter regions. The distinctive expression patterns regulated by estrogen in parenchyma and fat pad shed light on the need for both tissues to obtain normal mammary development.
- gene expression
- prepubertal heifer
- tissue type
the mammary gland is an organ that primarily develops postnatally under hormonal control. During the period from birth to puberty, the mammary parenchyma undergoes rapid, allometric, hyperplastic growth (8, 36), characterized by the proliferation and extension of mammary ducts into the surrounding mammary fat pad. In the absence of estrogen (4, 24, 40) or mammary estrogen receptors (5, 12), this early ductal growth does not occur. Proliferation of epithelial cells appears to be mediated by paracrine factors derived from the epithelial or stromal compartments of the mammary gland (1, 5, 8, 26). To gain a fuller understanding of estrogen regulation of mammary growth, it is necessary to evaluate estrogenic effects on gene expression in both the parenchyma and the fat pad. In vivo, these effects of estrogen may be direct or may be secondary to effects on estrogen-target genes within the mammary gland or target genes within peripheral tissues.
A biological response to estrogenic stimulation is initiated by interaction of estrogen with cognate receptors. Two estrogen receptor (ER) isoforms, ERα and ERβ, mediate a broad range of physiological and pathological effects (28, 33). In the classical pathway, estrogen-bound ER binds to estrogen-responsive elements (EREs) in the promoter region and nucleates the assembly of transcriptional complexes. ER can also be recruited to non-ERE sites via co-activators or co-repressors. Additionally, estrogen receptors in the plasma membrane have been found to regulate important biological functions, signaling through G protein activation and cross-activation of other membrane receptors such as the insulin-like growth factor I receptor (22, 33). The functional consequences of these various actions are transcriptional activity and expression of target genes. This in turn may alter the expression of additional genes, whose activity is modulated by the function of targeted genes. Genes that are regulated directly or indirectly by estrogen are referred to as estrogen-responsive genes.
The roles of estrogen in many biological processes, such as sexual development, mammary development, reproductive function and behavior, likely result from the cumulative effects of all estrogen-responsive genes (ERGs) and their interactions. Therefore, cataloguing ERGs in responsive tissues represents a first step in understanding mechanisms of estrogen action. This would allow the organization of ERGs into functional categories, shedding light on the actions of estrogen in complex organs such as the developing mammary gland. In humans and rodents, many attempts have been made to identify ERGs using microarray, serial analysis of gene expression (SAGE), and in silico prediction (11, 14, 18, 37). These genes are catalogued in two public databases. Over 1,500 ERGs, identified in human and rodent tissues, are catalogued in the Dragon Estrogen Responsive Genes Database (ERGDB) (38) (http://defiant.i2r.a-star.edu.sg/projects/Ergdb-v2/index.htm), whereas a more restricted listing of estrogen-targeted genes is catalogued in the ERTargetDB (19) (http://bioinformatics.med.ohio-state.edu/ERTargetDB/). A systematic search for ERGs in domestic animals has not been reported.
In this study, we performed a global evaluation of ERGs in bovine mammary gland, using a high-density oligonucleotide microarray. Our intent was to provide information needed to understand regulation of mammary growth and potential interactions between mammary parenchyma and fat pad. We examined the effect of estrogen status on gene expression profiles in the parenchyma and fat pad portions of mammary glands in both intact and ovariectomized prepubertal heifers. We were especially interested in differential responses of the two tissues to estrogen. The data obtained from this study will help to prioritize genes to pursue in the study of reciprocal signaling between parenchyma and stroma during ruminant mammary gland development. Comparison of the list of ERGs in ruminants with those in humans and rodents may also provide insight into the molecular basis for species differences in mammary gland development and parenchymal architecture.
MATERIALS AND METHODS
Animals and experimental design.
Sixteen 3-mo-old Holstein heifers (average wt, 106 kg) were purchased and then maintained at the Cornell Teaching and Research Center. The Cornell University Animal Care and Use Committee approved all procedures used in this study. Heifers were blocked by body weight and used in a 2 × 2 factorial design with ovarian status (intact or bilateral ovariectomy) as the first factor and estrogen treatment as the second (estrogen or control).
Bilateral ovariectomy was performed, under general anesthesia, when heifers reached ∼150 kg (4.5 mo old). Estrogen was supplied as 17β-estradiol dissolved in corn oil (Sigma, St. Louis, MO). Treatment was initiated 30 days after ovariectomy, using daily subcutaneous injection of 17β-estradiol (0.1 mg/kg body wt) for 3 consecutive days. Similarly, control animals received equal amounts of corn oil. Injections were administered at 24-h intervals, and heifers were killed ∼6 h after the final injection, or 54 h after initiation of estrogen treatment. This time was selected because previous studies demonstrated a significant proliferative response to estrogen within 2–3 days (8, 42), and our intent was to evaluate changes in gene expression coinciding with increased proliferation of mammary epithelial cells. Heifers were killed at the Cornell University abattoir by stunning with a captive bolt and exsanguination. At slaughter, udders were removed and weighed, and samples were obtained from midregions of the parenchyma and fat pad tissues of the left mammary glands. Tissues were snap frozen in liquid N2 until extracted for total RNA.
Isolation of total RNA and generation of biotin-labeled cRNA.
Total RNA was purified using an RNeasy Mini kit (Qiagen, Valencia, CA). Residual genomic DNA in the total RNA samples was removed by incubation with 4–10 units DNase I per 100 μg total RNA (Ambion, Austin, TX) at 37°C for 30 min. RNA integrity was verified using the Bioanalyzer 1000 (Agilent Technologies, Palo Alto, CA), and the concentration was determined using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Rockland, DE).
Biotin-labeled cRNA was generated in a two-step procedure. Initially, first-strand cDNAs were synthesized from 4.0 μg of total RNA using 1.0 unit of SuperScript II RT (Invitrogen, Carlsbad, CA) in the presence of 100 pmol T7 promoter Oligo dT primer. After second-strand synthesis, the DNA was purified with a DNA Clean and Concentrator-5 kit (Zymo Research, Orange, CA). The purified double-stranded cDNA was transcribed into cRNA using the MEGAscript in vitro transcription kit (Ambion). The in vitro transcription reaction was carried out in a total volume of 23.0 μl consisting of 3.0 μl of double-stranded cDNA, 2.3 μl of 10× Ambion reaction buffer, 2.3 μl of 10× Ambion T7 enzyme mix, and 15.4 μl of NTP labeling mix (7.5 mM ATP, 7.5 mM GTP, 5.625 mM UTP, 5.625 mM CTP, 1.875 mM biotin-16-UTP, and 1.875 mM biotin-11-CTP). The in vitro transcription reaction was incubated at 37°C for ∼16 h in a thermocycler. The cRNA was purified with an RNeasy Mini kit (Qiagen). Generally, 40–60 μg of cRNA were synthesized using 4.0 μg of input total RNA. The size range of the cRNA, expected to be between 300 and 3,000 bp with the maximum intensity centered at least 1,000 bp, was verified by Agilent Bioanalyzer analysis. The biotinylated cRNA was fragmented to 50–200 bp.
Oligonucleotide microarray, hybridization, image acquisition, and data analysis.
A total of 45,383 unique bovine sequences, including 40,808 tentative consensus sequences (TCs) from The Institute for Genomic Research (TIGR) Bos taurus Gene Index (BtGI; http://www.tigr.org), and 4,575 singleton expressed sequence tags were used to construct high-density bovine oligonucleotide DNA microarrays. Oligonucleotides (oligos) were synthesized in situ using photo deprotection chemistry with the Maskless Array Synthesizer system (NimbleGen, Madison, WI). Each oligo synthesized represented a feature on the hybridization surface of the microarray, with four features for every oligo. The feature size for the microarrays was ∼16 × 16 μm, and there were ∼345,000 features within a 17 × 13 mm array area. Overall, the features were populated with two different 60-mer oligos designed for each TC and one oligo for each singleton (86,191 distinct oligos in total). There were eight features for each TC (2 oligos with 4 features per oligo) and four features for each singleton (1 oligo with 4 features per oligo).
The microarrays were prehybridized with 1× MES hybridization buffer (100 mM 2-morpholinoethanesulfonic acid, 1.0 M Na+, 20 mM EDTA, 0.01% Tween 20), 40 μg of herring sperm DNA, and 200 μg of acetylated bovine serum albumin (BSA) at 45°C for 15 min followed by hybridization with 10 μg of denatured and fragmented cRNA per microarray at 45°C for 16–20 h with constant rotation. After hybridization, the microarrays were immediately washed extensively with nonstringent buffer (6× SSPE, 0.01% Tween 20) at room temperature followed by a stringent wash (100 mM MES salt and free acid solution, 0.1 M Na+, 0.01% Tween 20) at 45°C. After the final rinse with the nonstringent buffer, the microarrays were stained with 1× stain buffer (100 mM MES, 1 M Na+, 0.05% Tween 20, 50 mg/ml BSA, and 1 mg/ml Cy3-streptavidin) at room temperature for 25 min. The stain buffer was removed, and the microarrays were rinsed once more with nonstringent buffer. The microarrays were immediately dried under a stream of argon gas and scanned using an Axon GenePix 4000B scanner (Molecular Devices, Union City, CA) at 5-μm resolution.
Real-time RT-PCR analysis was carried out with the iQ SYBR Green Supermix kit (Bio-Rad, Hercules, CA) using 200 nM each amplification primer (see Table 2) and the first-strand cDNA (100 ng of the input total RNA equivalents) in a 25-μl reaction volume. The amplification was carried out on the iCycler iQ Real Time PCR Detection system (Bio-Rad) with the following profile: 95°C for 60 s and 40 cycles of 94°C for 15 s, 60°C for 30 s, and 72°C for 30 s. Melting curve analysis was performed for each primer pair. Reaction products were sequenced using a CEQ8000 automated sequencer and DTCS Quickstart Chemistry (Beckman Coulter, Fullerton, CA) to verify amplification of the appropriate transcripts. Expression levels of β-actin were used as endogenous controls within each sample. Relative gene expression was calculated by use of the 2−ΔΔCT method (where CT is threshold cycle) (23).
In addition to using real-time RT-PCR for validation of microarray results, we assayed expression levels of proliferating cell nuclear antigen (pcna) and ATP-binding cassette family A, member 3 (abc3), using the conditions described above, except that hybridization temperatures were 54°C (pcna) or 57°C (abc3) and sample concentrations were determined from a standard curve of known target cDNA copy numbers (1 × 102 to 1 × 107 molecules), which was analyzed simultaneously with samples. Standards for each gene were prepared from PCR amplicons purified using the QIAquick purification kit (Qiagen). Product concentrations were determined using the 2100 Bioanalyzer (Agilent Technologies). Primers for pcna amplification were 5′-TCGTCTCAGGCGTTCATAGTC-3′ (sense) and 5′-AACATGGTGGCGGAGTCG-3′ (antisense), yielding a product of 118 bp. Primers for abc3 amplification were 5′-GCCACCTTCCTCGTTGTC-3′ (sense) and 5′-AAGTTGCTCACTGCCATCC-3′ (antisense), yielding a product of 125 bp. For pcna, PCR efficiency was 97%, and the linear correlation coefficient was 0.997. For abc3, PCR efficiency was 92%, and the linear correlation coefficient was 0.993. Reaction products were sequenced to confirm specificity of amplification.
Microarray data were extracted from the raw images with the use of NimbleScan software (NimbleGen, Madison, WI). Preliminary data assessment was conducted using SAS (34) procedures. A normal distribution was obtained by converting the data to a log2 scale (16). Evaluation of data by probe set was performed, and hybridizations for probes having a coefficient of variation on the log scale <0.1 were excluded from analyses. A total of 2,630 probes were excluded from 32 hybridizations of 86,191 unique probe sets, or <0.1% of 2,758,112 probe-hybridization combinations were removed from the analysis.
Two mixed linear models were used in the analysis, based on the methods of Wolfinger et al. (41). The first model was fit to account for systematic differences that are fit across all DNA probes. Specifically, (1) where yijkl is the base-2 logarithm of the observed intensity for probe i, tissue j, cow k, and slide l. Here, μ represents the overall mean across all probes and slides, pi is the fixed probe effect, tj is the fixed issue (fat pad or parenchyma) effect, ck is the random cow effect, sl is the random slide effect, and εijkl is the random residual. Due to computational constraints, an analysis using SAS was not possible. Variance components were estimated, and solutions were obtained using MTDFREML (6). Fraction of variance accounted for by the cow effects was ∼72%, for slide effects was <1%, and for random residual variance was 18%.
Next, we let r represent the realized residuals from the model (Eq. 1), calculated by subtracting the fitted values for each effect in the model from the observed value. Additionally, residual effects were adjusted to have a mean of zero for each probe. Then, the second-level mixed linear model was fit independently for each probe using SAS Proc Mixed. Specifically, for probe i (2) where μi is the overall mean; cij is the random cow effect; (t·s·e)iklm is the fixed interaction of tissue k (fat pad or parenchyma), surgery l (intact or ovariectomized), and estrogen treatment m (estrogen or control); and γijklmn is the random residual. Notice that all effects were fit independently for each probe. Contrast statements corresponding to the following effects were calculated for every probe as linear combinations of the solutions for the three-way interaction (tissue, surgery, and estrogen main effects) and two-way interactions of tissue and estrogen, surgery and estrogen, and surgery and tissue. Additional contrasts were calculated for every probe to estimate estrogen effects within the following: fat pad, parenchyma, intact, and ovariectomized effects.
Experiment-wise P values were determined using permutation. The permutation algorithm used was based on the proposed method of Xie et al. (43). Essentially, the least significant fraction of one-half of the data, based on the minimum nominal P values for each probe, was identified. These probes were considered to be equally expressed across the contrasts evaluated. The test statistics generated by permuting these data were presumed to represent data from the null hypothesis (i.e., no treatment effects). The assignment of data from a slide to 1 of the 32 treatments (i.e., the original definition of a slide) was made at random across all of the equally expressed probes for a permutation. The values from such permutations were combined to represent one sample from the appropriate distribution, so the sample sizes represented the same number of probes. The P values corresponding to maximum significance observed by contrast were stored. This process was repeated 10,000 times, and adjusted P values were determined as the percentile of the nominal P value in this list of 10,000 values for each contrast separately. Hierarchical clustering of genes that were differentially expressed was performed using Pearson correlation in the program GeneSpring version 7.3 (Agilent Technologies).
Analysis of treatment effects on transcript abundance of pcna and abc3, generated by real-time RT-PCR, was by two-way ANOVA. The impact of estrogen treatment was tested as a main effect. An effect of ovariectomy was tested by pairwise comparison of intact (oil) and ovariectomized (oil) treatments using Scheffé's test. A one-sided test was used to test the hypothesis that estrogen increased and ovariectomy decreased abundance of pcna and abc3 transcripts, reflecting an anticipated impact on epithelial cell proliferation.
Treatment-induced effects on proliferation of mammary epithelial cells.
Fifty-four-hour treatment of intact and ovariectomized heifers with exogenous estrogen increased the quantity of pcna and abc3 mRNA (P < 0.05, Fig. 1). However, there was no effect of ovariectomy on abundance of these transcripts. Because these genes are related to cell proliferation, the data strongly support an estrogenic induction of mammary cell proliferation within the 54-h period of estrogen treatment.
Global analysis of tissue-specific genes and ERGs in the mammary parenchyma and fat pad of both intact and ovariectomized heifers.
A SAS Mixed Model procedure with permutation testing was used to identify significantly altered genes. Three main effects (estrogen, tissue, surgery) and eight interactions (estrogen by tissue, estrogen by intact, estrogen by ovariectomy, estrogen by parenchyma, estrogen by fat pad, tissue by intact, tissue by ovariectomy, and estrogen by tissue by surgery) were analyzed. The data were analyzed based on measurements for hybridization to individual oligos. A heatmap of gene expression is provided in Fig. 2. Expression of targets for 1,360 oligos was influenced by tissue or by estrogen treatment (P < 0.1, Supplemental Data File 1; the online version of this article contains the supplemental data). The majority of these transcripts displayed tissue specificity, being differentially expressed in the mammary parenchyma or fat pad (Fig. 2A). Additionally, targets for 180 oligos were found to be responsive to estrogen treatment (Fig. 2B). Most estrogen-responsive targets clustered into groups that were characterized by increased transcript expression in one or both tissues (clusters 1–4), but a small number of transcripts were downregulated in response to estrogen treatment (cluster 5).
From the 180 oligos whose targets were responsive to estrogen treatment, 124 ERGs were identified with significance at a permutation-corrected P < 0.1 with regard to estrogen-related interactions (Table 1). A gene was declared differentially expressed when one of the oligos representing the gene was found to be significant for the main effect or interaction. GenBank search (basic local alignment search tool; BLAST) was conducted for individual oligos to remove redundancy after they were identified as significant for a particular contrast. When a gene was represented by multiple sequences and/or multiple oligos, fold changes and P values were averaged for all oligos or sequences representing those genes that met the significance threshold (P < 0.1). The oligos and tentative consensus sequences that contributed to the expression data summarized in Table 1 are provided in Supplemental Data File 2.
As summarized in Table 1, the main effect of estrogen alone (est*) produced the greatest number of genes with significantly altered levels of expression, with 76 genes that met the threshold for statistical significance. This was followed by the effect of estrogen within the parenchyma (par*est) with 71 genes, and estrogen within the fat pad (fp*est) with 38 genes. Of the 109 genes significantly influenced by estrogen within the parenchyma and fat pad, only greb1 (TC280992, a homolog to human GREB1) responded significantly in parenchyma and fat pad. However, an additional 16 genes were commonly influenced in both tissue types, as evidenced by a significant estrogen main effect (Table 1, bottom 16 rows). The genes that were influenced by estrogen administration were categorized and depicted (see ⇓Figs. 4 and 5). It is clear that the parenchyma and fat pad responded differently to estrogen stimulation. Although the broad categorization of responsive genes was similar between the two tissues (see Fig. 4), there were appreciable differences in the affected genes and physiological response of the two tissues (see Fig. 5).
All 124 ERGs were searched against two public databases that catalogue ERGs, ERGDB and ERTargetDB. Thirty of these genes were previously reported to be estrogen responsive in other species, including well-documented ERGs such as stanniocalcin 1, α-1-anti-proteinase, progesterone receptor, nucleobindin 2, insulin-like growth factor I (igf1), and tissue factor pathway inhibitor (tfpi). However, the remaining 95 genes were not previously recognized to be estrogen responsive or involved in estrogen action. These novel estrogen target genes included 2,3,7,8-tetrachlorodibenzo-p-dioxin-induced poly-ADP-ribose polymerase (tiparp), bone morphogenetic protein (BMP)-binding endothelial regulator (bmper), and betaine-homocysteine methyltransferase 2 (btmt2). Although the vast majority of these 124 ERGs appeared to be upregulated by estrogen, a few genes such as rgs2, bhmt2, and galnt6 were downregulated by estrogen.
In contrast to the identification of 124 genes whose expression was influenced by treatment with exogenous estrogen, no genes in the mammary parenchyma or fat pad were found to be significantly influenced by the reduction in endogenous estrogens induced in this study by bilateral ovariectomy.
Fourteen genes that represented different levels of transcript abundance and different functional classes were selected for quantitative real-time RT-PCR confirmation (Table 2). Correlation analysis showed very good agreement between real-time RT-PCR and microarray analyses. Microarray results for all 14 genes tested were confirmed by real-time RT-PCR with regard to direction and significance of change. Linear regression analysis (n = 192) demonstrated a strong positive relationship between the two technological platforms with R2 = 0.75 (Fig. 3).
Estrogen exerts a broad range of physiological and pathological effects in different tissues mainly via interactions with its receptors ERα and ERβ and their binding to nuclear response elements, although transcriptional regulation in the absence of DNA binding and nongenomic effects via membrane receptors have been demonstrated (22, 33). As a mitogen, estrogen is involved in the regulation of cell proliferation and differentiation in many tissues. In the current study, we confirmed that estrogen stimulated proliferation of mammary epithelial cells in both intact and ovariectomized heifers. Estrogen increased abundance of transcripts for pcna and abc3 in mammary parenchyma. PCNA belongs to a class of proteins that play a key role in ensuring processivity of DNA polymerases (25). While pcna is not limited to actively dividing cells, cells in S-phase express the protein at rates ∼10-fold that of growth-arrested cells (7). As in many tissues, sorting mammary cells based on their ability to efflux Hoechst 33342 has resulted in identification of a side population of cells with enhanced expression of ATP-binding cassette protein G2 (ABCG2) that is markedly enriched for somatic stem cells (2). We evaluated expression of abc3, another ATP-binding cassette transporter (35), which appears to be more closely coupled to mammary proliferation than does abcg2 (A. V. Capuco, unpublished data) and may serve as a putative marker for bovine mammary stem/progenitor cells. On the basis of pcna and abc3 expression data from the current study, it can be concluded that estrogen administration increased proliferation of cells in the mammary parenchyma. This conclusion is further supported by increased bromodeoxyuridine (BrdU) labeling of mammary epithelial cells in response to exogenous estrogen (27a). The utility of assessing pcna or abc3 transcript abundance as indexes of mammary cell proliferation is confirmed by Pearson correlation coefficients with BrdU labeling of 0.94 and 0.98 for pcna and abc3, respectively. Although we were unable to demonstrate an effect of ovariectomy on pcna and abc3 transcript abundance, ovariectomy tended to decrease BrdU labeling of mammary epithelial cells, and systemic effects of treatments were further confirmed by decreased uterine weight in ovariectomized animals and increased uterine weight in estrogen-treated heifers (27a).
A liberal level of significance (P ≤ 0.1) was applied to our microarray analysis to protect against type II statistical error. Nonetheless, we were unable to detect an effect of ovariectomy on transcript profiles for the mammary parenchyma or fat pad. Because ovariectomy has been documented to inhibit prepubertal mammary development, the lack of an observed effect on mRNA expression in the mammary gland may be the result of several factors. First, an effect on the transcriptome may not have been evident due to the age of heifers when ovariectomized or the timing of expression profiling after ovariectomy. Berry et al. (4) reported that mammary glands of heifers ovariectomized before ∼6 wk of age were more severely affected by ovariectomy than animals ovariectomized between 8 and 12 wk of age. Because heifers in the current study were ovariectomized at ∼18 wk of age, the impact of ovariectomy on the mammary gland may have been lessened. As the first study to evaluate the impact of ovariectomy on transcript profiling in the bovine mammary gland, we were unable to assess whether the time after ovariectomy was appropriate to maximize potential effects on gene expression in the bovine gland. In the study by Berry et al. (4), pronounced effects of ovariectomy were evident when mammary growth was evaluated 3–5 mo after ovariectomy. In the current experiment, we attempted to evaluate changes in mammary gene expression at a time that was closer to ovariectomy, yet sufficiently delayed to permit a decline in circulating estrogens. Indeed, plasma estradiol-17β was below the level of detection, by radioimmunoassay, 1 wk after ovariectomy (27a). If ovariectomy elicited an effect on gene expression within the time course of our experiment, our data suggest that either ovariectomy alters mammary growth at the level of protein expression rather than mRNA expression or that the microarray technology employed was not capable of detecting a biologically significant change in mRNA expression. The former alternative is feasible, and we are currently evaluating the impact of treatment on micro-RNAs that can provide a mechanism for altering gene expression at the protein level (44). Our inability to detect changes in the transcriptome by microarray analysis may be because the technology is not sufficiently sensitive, or because the effect of treatment may have been targeted to a limited cell population, such as somatic stem cells or progenitor cells. Within the context of this study, we found that real-time RT-PCR was more sensitive than microarray analysis. For example, although estrogen did not influence levels of pcna and abc3 mRNA when assessed by microarray analysis, this effect was significant when assessed by real-time RT-PCR. Additional study is necessary to evaluate the impact of ovariectomy on gene expression and mammary gland development.
In bovine mammary gland, although ERα is expressed in significant quantity, ERβ mRNA is present at a very low concentration and no ERβ protein is detectable, indicating a relative lack of a role for ERβ in mammary gland development (10). Evidence indicates that, in response to estrogen, most genes regulated by ERα are distinct from those regulated by ERβ (20). Because ERβ seemingly is not expressed in bovine mammary glands, it is likely that the 124 ERGs identified in the current study represent genes that are directly regulated by ERα or by genes that are affected downstream. Potential downstream regulation includes genes that are downstream within targeted cells, as well as genes subject to estrogen-responsive, juxtacrine, paracrine, or systemic effects.
IGF-I has long been recognized to be the product of an ERG (30). IGF-I of mammary fat pad origin was proposed to be a paracrine mediator of estrogen action, and locally produced IGF-I is important in promoting ductal growth (32, 39). Our finding that igf1 was significantly upregulated by estrogen (est* and fp*est) provided additional evidence for a biological role in mammary gland development. It was documented that exogenous IGF-I can stimulate the phosphorylation of ER and induce the expression of ERGs (21). We are interested in further exploring the molecular basis of cross-talk between IGF-I and ER signaling pathways.
In this study, we identified 71 ERGs in parenchyma (par*est) and 38 ERGs in fat pad (fp*est). The genes influenced by estrogen in fat pad include cathespin L (ctsl), calcipressin 1 (dscr 1, a potent inhibitor of calcineurin), chondrolectin (chodl), collagens such as col6a3 and cola1l, monoamine oxidase a and b (maoa, maob), igf1, stem cell growth factor precursor (clec11a), tfpi, thrombospondin 2 (thbs2), and a few metabolic enzymes. The genes influenced by estrogen in parenchyma include centromere protein F (cenpf, a kinetochore protein closely associated with increased cell proliferation), kinases such as map3k8 and phosphatases, suprabasin and dermokine (2 physically linked genes associated with differentiating keratinocytes in stratified epithelia), rasd1, stc1, slpi, and keratin 1. The distinct expression profiles influenced by estrogen in parenchyma and fat pad were readily apparent by the limited number of genes that were commonly regulated in both tissues. Among the 17 commonly regulated genes were 2 genes related to cell signaling (greb1, cdca7l), 2 genes related to proliferative growth (bcl2l1, hmgb3), 3 genes related to cell metabolism (a4galt, atp1b3, sc4mol), 3 genes related to cell morphology (ckap4, sdc2, dynlt1), and 2 genes in the miscellaneous category that are chaperones or involved in protein folding (hspa5, pdia5). The importance of tissue or cell type in analysis of gene expression studies is emphasized not only by the differences in response to estrogen treatment but by the considerable differences in basal gene expression between the two tissues.
Pertinent to estrogenic stimulation of mammary growth was increased expression of greb1, an early ERG of unknown function and proposed mediator of the growth response to estrogen (15, 31), and bcl2l1, an anti-apoptotic regulator, in both tissues. Increased expression of genes that are potential stem cell markers or regulators (hmgb3, abc3, clec11a) in one or both of the tissues is suggestive of a targeted effect of estrogens on mammary stem cells or progenitor cells.
It is established that estrogens play an important role in promoting ductal growth during the prepubertal period and that the terminal ductular units in bovine mammary parenchyma are the primary sites of estrogen action (8). It is not surprising to see that genes associated with cell cycle controls were regulated by estrogen mainly in parenchyma (Figs. 4 and 5). The genes regulated by estrogen in fat pad, although they may not be directly involved in cell proliferation, may be important in providing signals for estrogen action in parenchyma via growth factors such as IGF-I and by initiating remodeling of the extracellular matrix (ECM). For example, thbs2 is an adhesive glycoprotein that inhibits angiogenesis and induces apoptosis by mediating cell-cell and cell-matrix interactions or by interacting with matrix metalloproteinases (MMPs), some growth factors, and matrix components such as fibronectin, fibrinogen and collagen. In murine skin tumors, thbs2 was highly expressed in fibroblasts and other stroma cells but was undetectable in epithelial cells (17). Tfpi is a serine protease inhibitor that has modest inhibitory effects on MMPs. In this study, the observation that estrogen treatment significantly upregulated expression of thbs2, tfpi, ctsl, collagens, and serpine2 (nexin or plasminogen activator inhibitor type 1) in fat pad provides compelling evidence that direct or indirect estrogenic effects on the fat pad are targeted toward ECM deposition and remodeling. Because ductal growth occurs within a sheath of loose connective tissue that penetrates between the surrounding adipocytes, seemingly isolating epithelial cells from the adipocytes (8), these effects on the ECM should lay the foundation for ductal growth (see Fig. 5). The molecular mechanisms by which estrogen induces tissue- or cell-specific gene expression are unclear. However, cell type-specific distribution of ER co-regulators (9, 13, 27) and receptor subtypes could account for some of the observed differential regulation, and response to indirect mediators, such as IGF-I, may account for others.
It is also interesting to note that mechanisms of estrogen action in bovine mammary gland may involve gene family members different from those earmarked as ERGs in human and mouse. For example, keratin 19 was proven to be an estrogen target gene in humans and rodents, but keratin 1, rather than keratin 19, was responsive to estrogen in bovine mammary gland. Similarly, dusp5 (not dusp3), ctsl (not ctsd), gria3 (not gria2), cenpf (not cenpa), galnt6 (not galnt4), hmg3 (not hmg2), and pdzk3 (not pdzk1) in bovine mammary gland were responsive to estrogen treatment.
The majority of ERGs identified in the present study were not previously identified as estrogen responsive in other systems. For example, BMP-binding endothelial cell precursor-derived regulator (bmper, a novel BMP-binding protein) was originally cloned in mice (29) and was shown to directly interact with BMP2, BMP4, and BMP6 and inhibit BMP-dependent development events such as smad activation and endothelial cell differentiation. BMP6 is closely correlated with tumor differentiation and skeletal metastasis. Previous studies demonstrated that estradiol could selectively increase the expression of BMP6. Activation through binding of ERα to the half-site ERE response element in the BMP6 promoter has been demonstrated (45). However, bmper has not been recognized as an ERG. To identify potential EREs in this and other novel ERGs identified in our study, we conducted in silico mapping and an ERE search for selected novel ERGs. The most 5′-ends (200-bp sequences) of these expressed sequence tags (ESTs) were mapped by BLAST on the latest bovine genome assembly (Btau_2.0) to retrieve 5- to 30-kb genomic sequences upstream of these ESTs. These genomic fragments were further analyzed using GenScan (http://genes.mit.edu/GENSCAN.html) and Dragon ERE Finder (3). Table 3 lists some of these search results. An unknown EST (TC269874) was shown to be upregulated by estrogen treatment, particularly in fat pad (fp*est). The EST itself is a noncoding sequence and may represent the 3′-untranslated region of a novel gene. GenScan identified at least four predicted peptide sequences in the 30-kb genomic fragment upstream from the most 5′-sequence of the TC, and multiple putative EREs existed in this region (Table 3). This gene warrants further study.
The microarray data are indicative of tissue specificity with regard to basal gene expression and supportive of tissue-specific changes in gene expression that provide for increased mammary epithelial cell proliferation in response to estrogen treatment (Fig. 5). Altered expressions of genes that are indicative of enhanced cell proliferation include increased cenpf, phlda2, and pbk and decreased qscn6. Putative genetic markers for stem/progenitor cells (hmgb3, abc3) in parenchyma suggest an impact of estrogen on mammary stem cells. Ductal elongation occurs within sheaths of connective tissue (Fig. 5), and our microarray data support increased ECM turnover in the parenchyma (e.g., increased hpse, serpina1, capn8) and increased deposition of ECM in the fat pad (e.g., increased col6a3, serpine2, bmper). ECM and connective tissue growth in the fat pad may be necessary to permit the penetration of mammary ducts into the surrounding mammary fat pad. Effects of estrogen were clearly tissue dependent, as a limited number of genes were commonly affected in both tissues. Expression of greb1 was commonly expressed in mammary parenchyma and fat pad. Although the function of greb1 is unknown, it is an early responder to estrogen, exists as multiple isoforms, and is hypothesized to play a role in the proliferative response to estrogen and cancer progression (15, 31). Finally, the role of the mammary fat pad as a source of paracrine stimulators of epithelial cell proliferation is suggested by increased expression of the stem cell growth factor precursor, clec11a, and igf1 in fat pad of estrogen-treated heifers. These growth factors may serve as estromedins, whereby the mammary fat pad may influence mammary epithelial development.
Future work could focus on the use of laser microdissection to compare the estrogen-induced expression profiles of individual cell types. This will enable further elucidation of cross-talk between different cell types or regions within a given tissue architecture. Additional work is necessary to characterize the mechanism by which ovariectomy impacts mammary gland development. In this regard, we are characterizing the impact of ovariectomy on micro-RNA expression within the bovine mammary gland.
Financial support for expression studies was provided by the United States Department of Agriculture-Agricultural Research Service (CRIS no. 1265-3100-086-00D). This research was also funded in part by a grant from the Cornell Center for Advanced Technology in Biotechnology, which is supported by New York State Science and Technology Foundation and by industrial partners. In addition, we thank Mike Fowler, Land O'Lakes Animal Milk Products Company, for support of this experiment.
Mention of trade names or commercial products in this publication is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the United States Department of Agriculture.
We thank Joy Castano and Larry Wood for excellent technical assistance and Dr. Luiz Coutinho for assistance with cluster analysis. We also acknowledge Dr. Bruce J. Aronow for helpful suggestions with data processing and presentation.
Address for reprint requests and other correspondence: A. V. Capuco, USDA-ARS, Bovine Functional Genomics Lab, Bldg. 200, Rm. 14, BARC-East, Beltsville, MD 20705 (e-mail:).
Article published online before print. See web site for date of publication (http://physiolgenomics.physiology.org).
- Copyright © 2006 the American Physiological Society