Estrogen-Dependent Gene Expression in the Mouse Ovary

Estrogen (E) plays a pivotal role in regulating the female reproductive system, particularly the ovary. However, the number and type of ovarian genes influenced by estrogen remain to be fully elucidated. In this study, we have utilized wild-type (WT) and aromatase knockout (ArKO; estrogen free) mouse ovaries as an in vivo model to profile estrogen dependent genes. RNA from each individual ovary (n = 3) was analyzed by a microarray-based screen using Illumina Sentrix Mouse WG-6 BeadChip (45,281 transcripts). Comparative analysis (GeneSpring) showed differential expression profiles of 450 genes influenced by E, with 291 genes up-regulated and 159 down-regulated by 2-fold or greater in the ArKO ovary compared to WT. Genes previously reported to be E regulated in ArKO ovaries were confirmed, in addition to novel genes not previously reported to be expressed or regulated by E in the ovary. Of genes involved in 5 diverse functional processes (hormonal processes, reproduction, sex differentiation and determination, apoptosis and cellular processes) 78 had estrogen-responsive elements (ERE). These analyses define the transcriptome regulated by E in the mouse ovary. Further analysis and investigation will increase our knowledge pertaining to how E influences follicular development and other ovarian functions.


Introduction
The importance of estrogen (E) in female reproductive endocrinology and in ovarian function has been well documented [1]. Estrogen signalling is primarily transduced by estrogen receptors (ER) a and b [2]. ER are members of a conserved superfamily of ligand activated transcription factors. The effects of E on ER are exerted through a complex array of convergent and divergent signaling pathways that mediate genomic events involved in regulation of mitogenesis, differentiation and apoptosis [3,4]. The interaction of E and ER with specific DNA sequences called estrogen responsive elements (EREs), constitutes a primary genomic signaling pathway [3,4]. The ERE-bound ER recruits an ensemble of co-factors responsible for the alteration of local chromatin structure and interaction with the basal transcription machinery.
Much effort has been invested in developing techniques to identify genes of interest, such as, Northern blot, semiquantitative RT-PCR and serial analysis of gene expression (SAGE). Several candidate genes, such as cyclin D2 [5], growth and differentiation factor 9 (GDF-9) [5], forkhead box transcription factor (Foxo1) [6], inhibina and inhibinbB [7] were shown to be regulated by E using these approaches. However, many downstream gene targets of E were unlikely to have been discovered using this approach. The aromatase knockout (ArKO) female mouse is deficient in aromatase activity postnatally and therefore is an excellent model to allow us to define E-dependent ovarian genes in adulthood. Using RT-PCR we identified genes usually associated with male reproduction such as Sox9, DAX1, liver homolog-1 (Lrh-1) and Mullerian-inhibiting substance [8] as being significantly increased [8] in the ArKO ovary. Steroidogenic enzymes such as 17a-OHase, 17ß-Hsd1, and 17b-Hsd3 mRNA's were also increased and the patterns of expression of the steroidogenic enzymes responsible for androgen biosynthesis in the ArKO ovaries correlated with increased serum testosterone levels in ArKO females [8].
The development of microarray technology now enables the simultaneous measurement of thousands of gene transcripts in a biological sample [9]. Therefore, in order to identify E-dependent genes in the ovary, this study utilized a microarray approach profiling wildtype (WT) and ArKO ovaries. The study aimed to 1) confirm the preliminary observations on estrogen dependent genes, 2) provide novel information about genes that can be used to unravel the mechanism of E in maintaining the female gonad, and 3) provide new insights into the regulation of ovarian follicular development.

Animals
Wild-type (WT) and ArKO mice on a J129/C57B6 background were maintained under specific pathogen-free (SPF) conditions, on a 12L:12D regimen and fed ad libitum a soy free mouse chow (Glen Forrest Stockfeeders, Western Australia). All animal procedures were approved by a Monash University Animal Ethics Committee (Project number: MMCB2002/37) and were carried out in accordance with the Australian Code of Practice for the Care and Use of Animals for Scientific Purposes.

Experimental design
The experimental approach is summarized in Figure 1. WT and ArKO mice, 16 weeks of age (n = 3/grp) were killed by CO 2 asphyxiation. The abdomen was opened and both ovaries from each animal were collected and snap frozen in liquid nitrogen. RNA was extracted from individual ovaries using a phenolchloroform-based method (Ultraspec; Fisher Biotech, Subiaco, Western Australia, Australia). Quantification of RNA concentration and purity was measured using the NanoDrop spectrophotometer (Thermo Scientific). The quality of the mouse RNA was ascertained using the Agilent Bioanalyser 2100 using the NanoChip protocol. A total of 500ng RNA was amplified and labeled using the Illumina TotalPrep RNA Amplification kit (Ambion) as per the manufacturer's instructions. A total of 1.5ug of labelled cRNA was then prepared for hybridisation to the Sentrix Mouse-6 Expression Beadchip (v1.1) by preparing a probe cocktail (cRNA @ 0.05ug/ul) that includes GEX-HYB Hybridisation Buffer (supplied with the beadchip). A total hybridisation volume of 30ul is prepared for each sample and 30ul loaded into a single array on the Sentrix Mouse-6 Expression Beadchip (v1.1). The Sentrix Mouse-6 Expression Beadchip (v1.1) allows for six samples and targets 48,500 unique well-documented RefSeq transcripts. A total of 6 different labelled samples can be loaded into 6 individual arrays per beadchip. The chip is hybridised at 58uC for 16 hours in an oven with a rocking platform. After hybridisation, the chip is washed using protocols outlined in the Illumina manual. Upon completion of the washing, the chips are then coupled with Cy3 and scanned in the Illumina BeadArray Reader. The scanner operating software, BeadStudio, converts the signal on the array into a text file for analysis. All data is MIAME compliant and the raw data has been deposited in a MIAME compliant database ArrayExpress (ArrayExpress accession: E-MEXP-2824). The Illumina microarray was performed by the Australian Genome Research Facility (AGRF) Melbourne, Australia.

Statistical analysis for microarray data
Prior to differential gene expression analysis, quality control diagnostic measurements were run on the raw data generated from the Illumina chips. 6 Sentrix Mouse-6 Expression Beadchip (v1.1) samples were normalised using lumi package (part of bioconductor) in R statistical software and were subsequently imported into GeneSpring GX software for further analysis. The significant E-dependent differentially expressed genes (DEG) were obtained from comparing KO samples to WT samples and were defined by a fold difference of 2.0 and a P-value cut-off of 0.05 (derived from ANOVA). Genes that met these parameters Figure 1. Flowchart of microarray data ranking and analysis. A) 450 E-dependent DEG list were identified by microarray in ovaries of ArKO vs WT as having 62-fold change expression with p-value,0.05. B) The E-dependent DEG list was annotated using Gene Ontology and analysed for molecular function, cellular component and biological processeses. C) E-dependent DEG from the list were ranked based upon fold change of up/ down regulation (Fold change); significance of the change (P-value); 5 major GeneGo biological processes genes: Hormonal processes, Reproductive processes, Sex determination and differentiation, Apoptosis and, Cellular processes. D) Genes that possess EREs. Genes identified from the list using these methods were compiled to form an E-dependent DEG with ERE Shortlist which was used for downstream analysis of gene networks and pathways affected by E using GeneGo pathway analysis. The E-dependent DEG list can be found in Table S1  Processes categories. A) Pie chart shows the distribution of the 388 E-dependent DEG that were matched to a Molecular Function using GO [6]. B) Pie chart shows the distribution of the 319 E-dependent DEG that were matched to a Cellular Component using GO [6]. C) Pie chart shows the distribution of the 302 E-dependent DEG that were matched to a Biological Processes using GO [6]. doi:10.1371/journal.pone.0014672.g002 Bioinformatic analysis and downstream annotation E-dependent DEG were ranked according to the methods described in the results section. Specific methods for Gene Ontology, GeneGo pathways analysis and Dragon estrogen response element (ERE) Finder version 6 are detailed below.
1. Gene Ontology (GO) [6] Analysis. E-dependent DEG lists were annotated according to the GO database (http://www. geneontology.org) using ontology categories for Molecular Function, Cellular Component and Biological Processes. The representation of E-dependent DEG within each ontology category was measured and statistical significance of the overlap was done using a hypergeometric p-value without multiple testing corrections, to determine the likelihood of coincidental overlap with ontology groups. A p-value lower than 0.05 indicates over representation of genes from the E-dependent DEG list within that particular category, suggestive of a functional effect. Annotation of the E-dependent DEG lists using GO was conducted by the AGRF Melbourne, Australia.
2. Dragon ERE Locator version 6. Nucleotide search from National Center for Biotechnology Information (http://www. ncbi.nlm.nih.gov/nucleotide/) was used to retrieve the sequence information of the differentially expressed genes within the selected Biological Processes categories in the FASTA format. These sequences were then examined for the presence of ERE binding sites using the latest Dragon ERE Finder version 6.0 (http://apps.sanbi.ac.za/ere/index.php) [10] based upon information present in the transcriptional regulation, from patterns to profiles (TRANSFAC) database [11]. The detection algorithm was tested on several large datasets and achieved a sensitivity of 83% [10]. An ERE is defined as a site which contains the 17 basepairs consensus sequence with its flanking 20 nucleotides on either side against the promoter sequences of Eukaryotic Promoter Database (EPD) [10,12,13].
3. GeneGo pathway analysis. The E-dependent DEG list containing the 450 unique E-dependent DEG, complete with Illumina transcript identifiers, were uploaded from a Microsoft Excel spreadsheet onto Metacore 5.0 software (GeneGo pathways analysis) (http://www.genego.com). GeneGo recognizes the Illumina identifiers and maps the E-dependent DEG to the MetaCore TM data analysis suite, generating maps to describe common pathways or molecular connections between Edependent DEG on the list. Graphical representations of the molecular relationships between genes were generated using the GeneGo pathway analysis, based upon processes showing significant (P,0.05) association.

E-dependent Differentially expressed genes (DEG) in ArKO compared to WTs
Microarray analysis of whole genome expression in ArKO ovary compared to WT ovary identified 450 differentially expressed transcripts with a 62-fold expression difference that was significant (p,0.05 according to t-tests, n = 3) ( Table 1). This list is referred to as the E-dependent DEG list and can be found in full in the supplementary data. The 450 E-dependent DEGs represent ,0.94% of the total number of transcripts present on the Illumina BeadChip. Of these 450 E-dependent DEG, 291 were up regulated (within a range of 2 to 47-fold change), and 159 DEG were down regulated (range between 22 to 215.7-fold change) in the absence of E.

GO annotation of the E-dependent DEG list
We hypothesised that some of the E-dependent DEG would be involved in downstream molecular pathways important for maintaining the female gonad. In order to identify which types of global cellular processes or specific molecular functions were responsive to E, the E-dependent DEG list was annotated using the GO [6] database. Each of the 450 genes was assigned to Molecular Function, Cellular Component and Biological Process categories as designated by the GO database. GO categories significantly over-represented in the E-dependent DEG list are those that match a greater number of the Edependent DEG than would be expected by chance (based upon comparison to the total number of genes in the genome assigned to the category). Figure 2A summarises the annotation of the E-dependent DEG list with a Molecular Function using GO. 388 of the 450 genes were matched to categories and the distribution of these 388 Edependent DEG into GO Molecular Function categories. The five GO categories with the greatest proportion of E-dependent DEG were Binding, Catalytic activity, Signal Transducer activity, Transporter activity and Enzyme regulator activity. Figure 2B summarises the annotation of the E-dependent DEG list with a Cellular Component using GO. 319 of the 450 genes were matched to categories and the distribution of these 319 Edependent DEG into GO Cellular Component categories. The five GO categories with the greatest proportion of E-dependent DEG were Cell, Extracellular region, Organelle, Protein complex and Cellular component unknown. Figure 2C summarises the annotation of the E-dependent DEG list to GO Biological Processes. 302 of the 450 genes were matched to categories and the distribution of these 302 E-dependent DEG into GO Biological Process categories. The GO categories with the greatest proportion of E-dependent DEGs were Cellular process, Physiological process, Development, Regulation of Biological Process and Response to Stimulus.
Ranking of genes on the E-dependent DEG list to create an E-dependent DEG with ERE Shortlist As described in Table 1, 450 E-dependent DEG were identified in ArKO ovaries that showed expression 62-fold compared to WT, with a significance of p,0.05 (E-dependent DEG list, Table S1). In order to narrow the E-dependent DEG list to those most likely to represent direct biological E target genes, four approaches were used to create a ranked shortlist of E-dependent DEG (Figure 1). The approaches used were, 1) Fold Change: the E-dependent DEG list was sorted according to fold change in expression, to identify the genes with the greatest up regulation and down regulation; 2) Pvalue: the E-dependent DEG list was sorted according to the statistical significance of the expression difference; 3) GeneGo analysis: 5 major GeneGo biological processes genes: Hormonal processes, Reproductive processes, Sex determination and differentiation, Apoptosis and, Cellular processes and 4) Genes containing EREs: Dragon ERE Locator version 6 was used to identify genes that contain ERE promoters from the shortlisted genes from step 1, 2 and 3. Using these criteria the list was reduced fold change in expression. The twenty genes found to be most up regulated and the twenty found to be most down regulated in the DEG List based upon fold change values are summarised in Table 2 and Table 3, respectively. The fold change values range from 47 to 6.3 for the most up regulated genes and from 215.7 to 24.8 for the most down regulated genes.
2. Genes from the E-dependent DEG list showing the most significant changes in gene expression. The twenty genes found to be most significantly up regulated compared to control and the twenty found to be most significantly down regulated compared to control, based upon the lowest p-values, are summarised in Table 4 and Table 5, respectively. The p-values range from 0.000001 to 0.000214 for the most up regulated genes and from 0.000004 to 0.000203 for the most down regulated genes.
3. GeneGo. Having narrowed the E target genes identified using the Illumina microarray down to a workable shortlist of 78 E-dependent DEG with ERE (Table S2); the molecular interactions between shortlisted genes were determined. Of the total set of 78 genes, ,37 (,47%) were not previously documented as ovarian genes, as determined by an extensive literature review and cross-checking with online databases (PubMed and Ovarian Kaleidoscope Database). Therefore, the identification of known and unknown genes could reveal how E affects regulation in folliculogenesis. We used GeneGo Pathways Analysis software to analyse the E-dependent DEG Shortlist and create 'networks' best connecting the 78 unique genes. Table 6 lists the top five most relevant networks identified by the GeneGo software. All the molecules listed within each network are connected based upon their known relationships or functions within the GeneGo Pathways Knowledge Base. The gene content of the E-dependent DEG Shortlist is used as the input list for generation of biological networks using the Analyze Networks (AN) algorithm. This is a variant of the shortest paths algorithm with main parameters of 1) relative enrichment with the E-dependent DEG Shortlist, and 2) relative saturation of networks with canonical pathways (g-Score). These networks are built upon the E-dependent DEG Shortlist and unique for the Edependent DEG Shortlist (Detailed network object legend, Figure  S1). For instance, Network A (Figure 3) with the highest g-Score (43.64) ( Table 6) contained 11 genes. The top processes associated with this network included 'organ development', 'anatomical structure morphogenesis' and 'organ morphogenesis'. Both Network B and C (Figure 4 and 5) were noted as having genes that are related to 'response to stimulus' and 'inflammatory responses'. Similarly, Network D ( Figure 6) further emphasised processes associated with organ development, system development and anatomical structure development. Network E (Figure 7) was associated with cellular processes, such as, regulation of multicellular organismal process, regulation of cell migration and regulation of cellular component movement. All Networks: A to E (Figure 3, 4, 5, 6, 7) complemented processes that were identified in the GO analysis.

Discussion
In this study, a genome wide transcriptional profiling microarray analysis was used for the first time to specify E-dependent genes or associated targets of E action in the mouse ovary. First, our studies led us to identify a set of 450 E-dependent ovarian genes whose biological roles may have significant relevance to functional and structural development of the ovary in mice. Second, we identified a set of 78 of these genes were most likely to have direct biological E targets and out of these, approximately 37 genes were not previously described in the ovary. Female ArKO mice are infertile due to failure to ovulate [14]. They have elevated levels of circulating gonadotrophins and testosterone, severely underdeveloped uteri [14,15,16] and increased adiposity [17,18].  The molecules that make up the top 5 Networks (A-E) identified using the E-dependent DEG with ERE Shortlist and ranked by GeneGo according to the g-Score are shown. The g-Score modifies the z-Score based on the number of Canonical Pathways used to build the network. If a network has a high g-Score, it is saturated with expressed genes (from z-Score) and it contains many Canonical Pathways. Sorting the The ovaries possess hemorrhagic cysts and apparent sex-reversal of the ovarian somatic cells; that is granulosa cells become Sertolilike cells [8,20,21]. The ArKO ovaries also displayed an infiltration of macrophages and contain substantial collagen deposition, a characteristic of tissues that are fibrotic [14]. It was logical to ask the question: were there changes in genes associated with these processes?

Masculinization
One of the major genes associated with masculinization, Sox9 [19] was confirmed to be expressed in the ArKO ovary from our microarray analysis results. Early in gonad development soon after the expression of SRY begins, Sox9 expression is strongly upregulated in Sertoli cells; Sox9 expression is maintained in the testis, whereas it is downregulated in the ovary [21]. The presence of the Sertoli cell marker Sox9 in the ArKO ovary confirmed the presence of testicular cell-like cells [19,20]. 17a-Hydroxylase (17a-OHase) and 17b-Hydroxysteroid Dehydrogenase type-3 (17b-Hsd3, the Leydig cell-specific) are two important steroidogenic enzymes in the production of testosterone. ArKO mice also showed increased levels of 17a-OHase and 17b-Hsd3 at 10 weeks of age, when Sertoli-and Leydig-like cells are beginning to populate the ovary [19] which were correlated with high serum testosterone levels [8]. Based on the known mechanisms maintaining Sertoli cell development and function in the normal (Wt testis), we identified two very intriguing novel targets of E that were significantly increased in the ArKO ovary, Claudin-11 (Cldn11) and platelet-derived growth factor receptor a (Pdgfra). Cldn11 is a protein component in tight junctions between Sertoli cells which are important for the maintenance of the blood-testis barrier [22,23] and Pdgfra was identified previously as a key player downstream of Sry in testis organogenesis and Leydig cell differentiation [24]. In the ArKO ovary, granulosa cells become Sertoli-like cells [8,19,20] and interstitial cells were identified ultrastructurally as Leydig cells [19]. However, to date, the exact role of E in the somatic cell trans/de-differentiation is unclear, although the fact that granulosa and interstitial cells express ER [25] points to a direct action on the gonadal somatic cells to maintain a ovarian phenotype. It is possible that the oocyte, which also expresses ER [26], plays a part in maintaining the ovarian phenotype with many recorded cases of premature oocyte loss leading to transdifferentiation of granulosa cells into Sertoli-like cells [8,19]. We therefore, hypothesized that E actively down regulates genes involved in masculinization to suppress the testicular phenotype, and that this action could be either direct on the somatic cells or via the oocyte or both.

Folliculogenesis
From the morphological and stereological assessments made in previous studies [14,19], follicular growth beyond the preantral stage is unsuccessful in the absence of E, with most succumbing to atresia, developing into hemorrhagic cysts or taking on a testicular phenotype (as discussed above). We hypothesized that this phenomenon occurs due to the loss of oocytes in the ArKO ovary and this was supported by the array of genes that were identified in this study. For example, Paired basic amino acid cleaving system 4 isoform 1 (Pace4) or Proprotein convertase subtilisin/kexin type 6 (Pcsk6) is a proprotein convertase required for producing active forms of important TGFb ligands in the ovary [27]. The levels of this enzyme were increased in the ArKO ovary. Steady state mRNA levels of Pcsk6 in granulosa cells are normally required to be suppressed by the oocyte during preantral to antral follicle transition [28]. Thus, Pcsk6 represents another granulosa cell transcript regulated by oocytes as well as E, but exhibiting a novel expression pattern during follicular development. Cytochrome P450 lanosterol 14 alpha-demethylase (Cyp51), which was decreased in the ArKO ovary, has long been known for reactions involved in drug metabolism and synthesis of cholesterol, steroids and other lipids. Recent studies on Cyp51 show that it promotes the initiation of meiosis in early folliculogenesis [29,30]. The decrease in expression of Cyp51 in the ArKO ovary implies that E has a role in controlling early follicular development. This observation is also supported in primate studies by Zachos et al. [31], in which fetal baboons deprived of E showed a 50% reduction in the number of primordial follicles which was restored by E 2 administration. Similar data for the ArKO mouse ovary has not been reported, although it has been shown that E influences the breakdown of oocyte nests prior to the formation of primordial follicles [32,33].

Novel Network Processes
A significant over representation of genes involved in organ development and morphogenesis, inflammatory responses and cellular processes were identified in both the Gene Ontology annotation and GeneGo Pathway analyses of differentially expressed genes in ArKO ovary. 1. Hedgehog signalling. One of the affected pathways in the ArkO ovary is the Hedgehog-Patched (HH) signalling pathway, which is involved in ovarian follicle development. Expression of HH ligands, Sonic hedgehog (Shh), Indian hedgehog (Ihh) and Desert hedgehog (Dhh) are developmentally regulated in the adult ovary [36]. Other components of the HH pathway such as Hedgehog receptors (PTCH1 and PTCH2) and signal transducer Smoothened (SMO) are expressed by granulosa cells of the follicle [34]. Hedghog signalling in the mammalian ovary plays a role in the communication between granulosa cells and the developing theca cells [35]. A recent study by Ren and colleagues [36], showed that Amhr2 cre/+ SmoM2 mutant mice that constitutively expresses the dominant active form of SMO in the ovary, were anovulatory because of a dramatic reduction in the muscle cells of the theca layer surrounding the developing follicles. These muscle cells within the theca layer play a role in the release of the oocyte at the time of ovulation. One of the smooth muscle markers, actin, gamma 2, smooth muscle, enteric (Actg2) was found to be significantly decreased in the Amhr2 cre/+ SmoM2 mutant mice [36] and this was also the case in the ArKO ovary (Table 4), suggesting E might play a role in regulating the smooth muscle cells within the ovary during ovulation. In our study, Ihh mRNA expression was found to be decreased in the ArKO mouse ovary (Table 4). The reduced Ihh production in the ArKO ovary implies an inhibition of HH signalling and a hormonal regulation of the Hedgehog system ( Figure 3).
2. Inflammatory response and fibrosis. It was not surprising to see the inflammatory responses network filling two of the top 5 network processes (Figure 4 and 5) in the ArKO ovary. Previous studies have shown an increased infiltration of macrophages and deposition of collagen in the ArKO ovary, suggesting that in the absence of E, the ArKO ovary begins to degenerate [19]. Mast cells, macrophages and granulocytes usually accumulate in rat ovaries during the preovulatory period [37]. Their presence in ArKO ovaries indicates that the chemokine signals for ovulation including an increase in the abundance of the immune mediators have been activated, but the follicles have been unable to respond. The increased deposition of collagen in the absence of E suggests that matrix metalloproteinases (MMPs), including the collagenases, have not been activated. This was reflected in the network processes ( Figure 6), as MMP-2, MMP-9, tissue inhibitors of MMP 1 (TIMP1) and TIMP2 were either directly or indirectly affected in the absence of E. Detailed studies on the hormonal regulation, site of synthesis, and absolute requirement on individual MMPs and TIMPs, for successful ovulation remain to be elucidated. The improved ovarian morphology in ArKO mice treated with E [8], correlated with a decrease in 3. WNT signalling. WNT signaling influences cellular differentiation and proliferation in the embryonic and adult ovaries. In this study, we identified a significant decrease in Secreted frizzled-related protein 4 (SFRP4) mRNA expression in the ArKO mouse ovary. SFRP4 is a member of the Frizzled like cystein rich domain family and a modulator of non canonical WNT signalling [38]. It has been shown that SFRP4 transcripts colocalize in the ovary to sites of Fz-1, Wnt4, and Fz-4 expression, suggesting that SFRP4 may regulate signaling through these factors in the ovary [38,39,40]. Both SFRP4 and Fz-1 are expressed in granulosa cells of large follicles suggesting a potential regulation by SFRP-4 of Fz-1 signals that may impact granulosa cell differentiation or expression of genes involved in follicle rupture and ovulation [38,40]. Wnt4 is a member of the WNT family of intracellular growth and differentiation factors, which regulate several key developmental steps. Deficiency of Wnt4 leads to partial female to male sex reversal and a marked reduction in oocyte number [41]. In addition, Fz-4 deficient mice showed defective corpora lutea [42]. It is therefore, tempting to speculate that E is regulating SFRP4 to maintain the ovarian phenotype, The ArKO model will be a useful tool to further study the effect of loss of estrogen on WNT signaling. Other genes of interest This study identified genes with known or suspected roles in ovarian functions as stated above. Interestingly, there were 37 genes from the E-dependent DEG Shortlist that possessed EREs that were completely novel to ovarian function. The roles of these genes in ovarian function will require further research and verification. There were also genes (DAX1 [8], MIS [8], ERa, [8], ERb [8], GDF9 [43], and Cyclin D2 [5]) that were expected to be regulated (up or down) in the absence of E as a result of previous studies, which were not detected in the present study. It is likely that some genes will be missed, because they may be affected by E at earlier time-points or their expression may be below the level of detection in this study.

Conclusion
Our success in identifying E-dependent ovarian genes validates the usefulness of whole-genome arrays. We identified more than double the number of E-dependent ovarian/oocyte factors previously reported. The E-dependent ovarian genes identified in this study provide attractive candidates for studying female infertility [44]. In addition, analyses can be easily expanded by inclusion of additional expression profiles, such as, a set of E 2 replacement ArKO ovarian genes or refined by identifying which ER subtype is responsible for specific E-dependent ovarian genes.
Such study will be required in order to define more precisely mechanisms of ovarian function, as well as provides new targets for the development of contraceptives and highlights genes previously uncharacterised in the ovary, that may be involved in the pathophysiology of female reproductive disorders. Furthermore, these E-dependent expression datasets will enable further gene discovery focusing on microRNAs and non-coding RNAs.