Androgen-induced masculinization in rainbow trout results in a marked dysregulation of early gonadal gene expression profiles

Background Fish gonadal sex differentiation is affected by sex steroids treatments providing an efficient strategy to control the sexual phenotype of fish for aquaculture purposes. However, the biological effects of such treatments are poorly understood. The aim of this study was to identify the main effects of an androgen masculinizing treatment (11β-hydroxyandrostenedione, 11βOHΔ4, 10 mg/kg of food for 3 months) on gonadal gene expression profiles of an all-female genetic population of trout. To characterize the most important molecular features of this process, we used a large scale gene expression profiling approach using rainbow trout DNA microarrays combined with a detailed gene ontology (GO) analysis. Results 2,474 genes were characterized as up-regulated or down-regulated in trout female gonads masculinized by androgen in comparison with control male or female gonads from untreated all-male and all-female genetic populations. These genes were classified in 13 k-means clusters of temporally correlated expression profiles. Gene ontology (GO) data mining revealed that androgen treatment triggers a marked down-regulation of genes potentially involved in early oogenesis processes (GO 'mitotic cell cycle', 'nucleolus'), an up-regulation of the translation machinery (GO 'ribosome') along with a down-regulation of proteolysis (GO 'proteolysis', 'peptidase' and 'metallopeptidase activity'). Genes considered as muscle fibres markers (GO 'muscle contraction') and genes annotated as structural constituents of the extracellular matrix (GO 'extracellular matrix') or related to meiosis (GO 'chromosome' and 'meiosis') were found significantly enriched in the two clusters of genes specifically up-regulated in androgen-treated female gonads. GO annotations 'Sex differentiation' and 'steroid biosynthesis' were enriched in a cluster of genes with high expression levels only in control males. Interestingly none of these genes were stimulated by the masculinizing androgen treatment. Conclusion This study provides evidence that androgen masculinization results in a marked dysregulation of early gene expression profiles when compared to natural testicular or ovarian differentiation. Based on these results we suggest that, in our experimental conditions, androgen masculinization proceeds mainly through an early inhibition of female development.


Background
The embryonic gonad has the potential to develop into a fully functional organ able to produce the gametes necessary for sexual reproduction. Sex differentiation is a crucial step in this developmental process and is considered as the differentiation from a bipotential gonadal primordium towards a testis or an ovary. In teleostean fish, sex differentiation can be controlled by in vivo treatments with sex steroids (reviewed in [1]) as in reptiles and amphibians and to some extent in birds (reviewed in [2][3][4]). In fish, these steroid treatments are often able to induce fully functional sex-inversed phenotypes and these treatments have been widely used to produce all-male or all-female populations of fish for aquaculture purposes [5]. Many studies have been focused on the role of these hormones during gonadal sex differentiation highlighting for instance the crucial role of estrogens in ovarian differentiation [1]. However, most of the studies performed thus far were focused on a very small number of well characterized genes, proteins or hormones and mostly on natural gonadal differentiation.
Rainbow trout, Oncorhynchus mykiss, has a male heterogametic XY genetic system and we experimentally produced XX and YY males allowing the production of genetically all-male and all-female populations [6]. These all-male or all-female populations provide a unique opportunity to work on numerous animals for which the normal gonadal development as testis or ovary is known a priori. Using the extensive collection of expressed sequenced tags (ESTs) obtained through sequencing projects in trout as a resource [7,8], we designed and built a DNA microarray in order to characterize, on a genome-wide scale, the mechanisms by which 11β-hydroxyandrostenedione (11βOH∆4), a natural androgen in fish [9,10] is able to masculinize the embryonic ovary.
Using this genome-wide approach we characterized 2,474 genes (2372 microarray and 102 real-time RT-PCR gene expression profiles) with a clear differential temporal expression profile in females masculinized by androgen. We classified these genes in 13 different clusters of correlated temporal expression profiles, and searched within these clusters for significant enrichment in Gene Ontology (GO) terms. This strategy allowed us to define a few very clear biological trends potentially explaining how androgen induces masculinization of female fish. Our results clearly demonstrate that masculinization with androgen proceeds through a marked dysregulation of gene expression profiles, including a quick down-regulation of the ovarian pathway. Surprisingly, most of the genes over-expressed during natural testicular differentiation were not restored by the androgen-induced masculinization suggesting that the inhibition of female gonadal development is the main required step sufficient for building a testis.

Results
The complete dataset is available through the National Center for Biotechnology Information (NCBI), in the Gene Expression Omnibus database [11] under the GSE7018 accession number. After statistical filtering, 2,474 expression profiles (2372 microarray and 102 realtime RT-PCR gene expression profiles, data available as supplemental material in Additional file 1) were identified as being characteristic for either natural differentiation (ovarian or testicular differentiation) or androgeninduced masculinization (trans-differentiating gonads). Among these 2,474 expression profiles, 73% (1,805) were associated with genes with significant homologies with well characterized proteins in Swissprot or Prodom databases (the complete list of clones and their annotations is available as supplemental material in Additional file 2).

Biological sample clustering and histology
This analysis was carried out on fish sampled at several stages of development from the onset of the free swimming period (Day 0 = D0), when fish first started to be fed with the androgen treatment until 110 days after the beginning of the treatment (D110). Unsupervised hierarchical clustering of samples ( Fig 1A) reveals 3 main groups of correlated samples according to their global gene expression profiles i.e., late ovarian samples (D60 to D110, correlation coefficient R = 0.78), middle and late gonad samples of androgen-treated fish gonads (D27 to D110, R = 0.37) and middle and late testicular samples (D27 to D110, R = 0.26). All early samples (D0 to D12) cluster together with a weak correlation (R = 0.07), indicating that differences in the expression profiles of these early samples are rather large. Histological analysis of gonads at D12 (Fig. 1B, panels a, c, e) reveals characteristic features of differentiating gonads in control fish with the first appearance of ovarian meiosis and lamellar structures in females, or spermatogonia cysts in males. At this timepoint (D12), gonads of androgen-treated females appear as a thin structure characterized by scattered germ cells in a predominant stroma of conjunctive tissue with fibroblast like cells. After 90 days of treatment (Fig. 1B, panels b, d, f) the gonads of androgen-treated females display a classical testicular organization with cysts of germ cells engaged at various stages of meioses whereas the control male gonads show the same organization but with only gonial mitosis. At D90, the control female gonads contain previtellogenic oocytes surrounded by flattened granulosa cells and oogonia within clearly organized ovarian lamellae.

Global analysis of gene expression profiles
The 2,474 expression profiles were analyzed using a kmeans clustering (with k = 13) in order to individualize clusters of genes with similar expression profiles (Fig. 2). These expression profiles and the 13 k-means clusters are available online as a browseable file [12]. Among these clusters, clusters 1 to 4 are characterized by a specific high expression levels of a very large number of genes (N = 1,204) specific to the late ovarian samples (D60 to D110). These clusters could have been merged as their expression profiles seem very similar. In cluster 5 (N = 184) these high expression levels in the female group after D60 are also present in the androgen-treated group after D27. Cluster 6 (N = 89) contains genes with increasing expres-sion profiles starting from D16 in androgen-treated females and from D27 in control males. Clusters 7 (N = 132) and 8 (N = 77) are characterized by an early (starting from D27, cluster 7) or late (starting from D90, cluster 8) increase in gene expression specific to the androgentreated females. Genes in cluster 9 (N = 214) display a late down-regulation starting from D60 in control female gonads. Cluster 10 (N = 133) is the only cluster that does not show any difference between males, females and androgen-treated females and contains genes with continuously decreasing expression levels from D0 to D110. Cluster 11 (N = 172) contains genes down-regulated both in females (from D60) and androgen-treated females (from D27). Cluster 12 (N = 214) is characterized by gene Classification of gonad samples and histological analysis of some characteristic gonadal stages  expression levels specifically down-regulated (from D27) in androgen-treated females. However, this down-regulation is not maintained after the end of the treatment (i.e. after D90). Cluster 13 is the smallest cluster in terms of the number of genes (N = 55) and it displays high expression levels specifically in males throughout all the sampling times (from D0 to D110). Expression levels in androgen-treated females slightly increase after completion of the treatment (> D90).

Annotation of gene clusters using Gene Ontology (GO)
We then searched for Gene ontology (GO) terms significantly enriched in these groups of correlated expression profiles compared to the overall GO terms found. Among the 1,805 genes with an annotation we selected 1,276 unique genes annotated with an official name, and 1,133 were found to be associated with at least one GO category.
The top 5 significantly enriched GO terms for each of the clusters are described in Table 1 and 2 and the complete analysis is available online (as supplemental material in Additional file 3).
As gene expression in clusters 1 to 5 displayed similar expression profiles in the control female group, these clusters have been analysed both all together and separately ( Table 1). If considered as a homogenous cluster it displays considerable enrichments in GO terms likely to be important characteristics of these late ovarian stages. Among these GO terms the most representatives (see Table 1

A B
5730, see figure 3 for a detailed composition of this GO term). GO term 'mitosis' (ID 7067) is highly over-represented in cluster 3 (see Table 1 and figure 4 for a detailed composition of this GO term) with a 2.1 fold relative enrichment (p value of 4.10-3). Among clusters 1 to 5, only genes from cluster 5 are also highly expressed in the androgen-treated group. The main GO theme of cluster 5 is related to translation, with over-representation of GO 'ribosome' (ID 5840, see figure 5 for a detailed composition of this GO term), 'ribonucleoprotein complex' (ID 30529), and 'ribosome biogenesis and assembly' (ID 42254). Within the GO term 'ribosome', a large proportion of genes (7 out of 13) up-regulated in cluster 5 are also annotated as 'cytosolic large ribosomal subunit' (ID 5842, e.g. rpl5, rpl7, rpl13a, rpl17, rpl19 and rpl21, see Fig. 5). Cluster 6 contains genes that are up-regulated both in control males and androgen-treated females. It contains genes ( Table 2) involved in 'cell surface' (ID 9986), 'T cell activation' (ID 42110) and the 'immune response' (ID 6955). Clusters 7 and 8 are of particular interest as they display an up-regulation only in androgen-treated females. Cluster 7 is characterized by the GO terms 'extracellular matrix' (ID 5578), 'actin binding' and 'muscle contraction' (ID 3779 and ID 6936) including genes like myosins (e.g. myl6 and myl11), tropomyosins (e.g. tpm1, tpm3 and tpm4) and some muscle markers (e.g. calponin2, cnn2 and transgelin, tagln). The GO term 'extracellular matrix' (Fig. 6) contains 30 different genes, and 10 are specifically up-regulated in cluster 7 following androgen treatment in females (4.6 fold enrichment with a p value < 4.10-4). Many of these genes (e.g. col1a1, col1a2, col6a2, and mfap2), are 'structural constituents of the extracellular matrix' (ID 5201), or are involved in 'cell adhesion' (ID 7155) (e.g. col6a2, sparc and postn). The expression profiles of some of these genes related to the extracellular matrix are shown in figure 7 along with the tgfb1 expression profile that also shows a slight up-regulation from D60 to D110 following androgen treatment These GO terms were selected as the 5 most significant GO terms enriched in each K-means cluster (KM) 1 to 5 with a P-value < 0.05 and a relative enrichment RE > 2 (except for the analysis of the fusion of clusters 1 to 5). n/N: number of genes assigned with a specific GO term in the cluster (n) with regards to the number of all genes assigned with the same GO term in the 2,474 analyzed genes (N). Interesting GO terms in regard to gonad differentiation and development are underlined. Genes with GO terms in bold type are detailed in figures 3 to 5. GO categories containing fewer than 5 genes were excluded.
in females. Cluster 8 contains genes related to 'chromosome' (ID 5694), 'condensed chromosome' (ID 793) but also, albeit only a few (2 out of 6), to the 'progression through the first phase of meiosis' ('meiosis I', ID 7126, i.e., rad1 and the synaptonemal complex central element protein 2, syce2). Cluster 9 is characterized by enrichment in the GO terms 'organ morphogenesis' (ID 9887), 'development' (ID 48513) and 'structural constituent of cytoskeleton' (ID 5200). Cluster 10 contains genes involved in 'urogenital system development' (ID 1655) and 'lipid binding' (ID 8289). Clusters 11 and 12 contain genes that are also assigned with the GO term 'extracellular matrix' but more specifically in relation with degradation of this extracellular matrix i.e., GO term 'metallopeptidase' (5.4 and 3.8 fold enrichment respectively) with at least 3 matrix metallopeptidases (e.g. mmp13, mmp14, mmp19) that are directly involved in 'collagen catabolism' (ID 30574). Cluster 13 contains genes related to 'sex differentiation' (ID 7548) (e.g. amh, sox9, dmrt1, gata4, lhx9), and 'steroid biosynthesis' (ID 6694) (e.g. cyp17a1, cyp11b2, star and nr5a2). According to these annotations this cluster is likely to contain other These GO terms were selected as the 5 most significant GO terms enriched in each K-means cluster (KM) 6 to 13 with a P-value < 0.05 and a relative enrichment RE > 2. n/N: number of genes assigned with a specific GO term in the cluster (n) with regards to the number of all genes assigned with the same GO term in the 2,474 analyzed genes (N). Interesting GO terms in regard to gonad differentiation and development are underlined. Genes with GO terms in bold type are detailed in figure 6. GO categories containing fewer than 5 genes were excluded.
important genes also involved in testicular differentiation or in steroidogenesis regulation. Several expression profiles for such potential genes are shown in figure 8.

Validation and enrichment of DNA microarray data by real-time RT-PCR
Expression profiles of 102 genes, involved in early gonad development, were measured by real-time RT-PCR. These genes all belonged to the 13 distinct k-means clusters (see Additional file 1). Among these 102 genes, 84 were only measured by real-time RT-PCR and these gene expression profiles were thus added to the microarray dataset. The remaining 18 genes were common between the real-time RT-PCR dataset and the microarray dataset (see Table 3) and were used to validate our microarray dataset. Among the 18 common genes belonging to 11 out of the 13 distinct k-means clusters (e.g. 1, 2, 3, 4, 5, 6, 8, 9, 10, 11, and 13), 15 expression profiles were found to have a significant correlation between the two techniques (see Table 3 and Fig. 9). Only tfa (transferrin), timp2 (tissue inhibitor of metalloproteinase 2) and bzrp (benzodiazepine receptor, peripheral) expression profiles did not significantly correlate. As a result, among the 13 k-means groups, at least 10 groups (e.g. 1, 2, 3, 4, 5, 6, 8, 9, 11, and 13) contained genes with significant correlation between real-time RT-PCR and DNA microarray measurements.

Discussion
Our global approach based on gene expression profiling clearly reveals that, in our experimental conditions (11βOH∆4, 10 mg/kg of food for 3 months), the androgen masculinization does not induce a natural physiological response since the transcriptome of testicular transdifferentiating gonads is quite different from the one observed during natural testicular differentiation. These differences might be due to the non physiological dosage of androgen used in our experiment. A similar study using a lower dosage may help to clarify this issue, but this study Expression profiles of genes belonging to the Gene Ontology term "nucleolus" (ID 5730) Figure 3 Expression profiles of genes belonging to the Gene Ontology term "nucleolus" (ID 5730). Clones with similar gene names that belong to the same contig (ensemble of clones with overlapping sequences) but displayed very different expression profiles were excluded from the table (when the expression profiles were not too different the gene symbol are marked with an asterisk). Clones with similar gene names and belonging to different contigs were annotated as gene symbol_1 and gene symbol_2 as they could be considered as potentially duplicated genes or differential splicing forms. KM: K-means cluster number. N: number of genes belonging to the GO term in the specified cluster. Gene Ontology "nucleolus" (ID 5730) Gene Ontology "nucleolus" (ID 5730)

-1 +1
was first designed with an androgen dosage that is commonly used in rainbow trout aquaculture conditions. Whether the observed gene expression dysregulations are the reflection of a direct action of the androgens on the gonad, an indirect retro-control on the hypothalamuspituitary axis, or a conjunction of both, remains to be elucidated. However, the synthesis of Gonadotropin Releasing Hormone (GnRH) and of pituitary hormones is established very early during rainbow trout ontogenesis [13] with at least the synthesis of Follicle Stimulating Hormone FSH [14]. Thus indirect feedback effects cannot be totally excluded.
Due to the lack of specific Gene Ontology (GO) annotation for rainbow trout we linked the best blast hits of each clone sequence with a cross-species GO annotation. This strategy relies on the accuracy of the blast homology search and also on the resulting accuracy of the GO anno-tations with regards to their use in a fish species. However, even if this could lead to potential errors on a gene per gene scale, the global analysis and stringent statistical screen that we carried out enabled us to unambiguously assign most clusters with a clear biological theme. Among these GO categories, some were considered as biologically informative -i.e., not too general like for instance GO "physiological process" -and robust as they contain a sufficient number of different genes to support a potential biological meaning. We focused our analysis on these biologically informative GO categories.
With regards to the effects on the gonad, our analysis first reveals that female development is highly affected by the androgen treatment, with a down-regulation of most of the genes involved in early oogenesis stages. However within this analysis we did not characterize any cluster of early female-specific up-regulated genes potentially Expression profiles of genes belonging to the Gene Ontology term "mitosis" (ID 7067) Figure 4 Expression profiles of genes belonging to the Gene Ontology term "mitosis" (ID 7067). Clones with similar gene names that belong to the same contig (ensemble of clones with overlapping sequences) but displayed very different expression profiles were excluded from the table (when the expression profiles were not too different the gene symbol are marked with an asterisk). Clones with similar gene names and belonging to different contigs were annotated as gene symbol_1 and gene symbol_2 as they could be considered as potentially duplicated genes or differential splicing forms. KM: K-means cluster number. N: number of genes belonging to the GO term in the specified cluster.

-1 +1
involved in ovarian differentiation. Expression profiles of some female-specific candidate genes [15] (e.g. foxl2a and foxl2b, cyp19a1, fst, inha) were introduced in our analysis and all these genes were strongly and quickly inhibited by the masculinizing androgen treatment. But when pooled with our DNA microarrays dataset they did not form a tight cluster. This is probably because no additional similar expression profile was found within the DNA microar-ray dataset. This very small number of early femalespecific genes is in agreement with the small number of candidate genes known to be involved in the ovarian differentiation pathway [16] in comparison with the relatively high number of genes that are known to characterize testicular differentiation [17]. In agreement with this view, our analysis clearly characterizes a cluster displaying testicular-specific gene expression profiles, Expression profiles of genes belonging to the Gene Ontology term "ribosome" (ID 5840) Figure 5 Expression profiles of genes belonging to the Gene Ontology term "ribosome" (ID 5840). Clones with similar gene names that belong to the same contig (ensemble of clones with overlapping sequences) but displayed very different expression profiles were excluded from the table (when the expression profiles were not too different the gene symbol are marked with an asterisk). Clones with similar gene names and belonging to different contigs were annotated as gene symbol_1 and gene symbol_2 as they could be considered as potentially duplicated genes or differential splicing forms. KM: K-means cluster number. N: number of genes belonging to the GO term in the specified cluster.

KM N Clone name Gene Symbol Gene Name
Gene Ontology "ribosome" (ID 5840) Gene Ontology "ribosome" (ID 5840) -1 +1 containing both genes known to be involved in testicular differentiation (e.g. amh, sox9, dmrt1, gata4, lhx9) [18] and some potential new players revealed by our analysis. Interestingly the expression levels of all these genes are not restored by the androgen masculinizing treatment, and this could indicate that they are probably not necessary for early testicular differentiation in rainbow trout.
Among the gene clusters specifically up-regulated in females following masculinization with androgens, extracellular matrix, muscle markers/cytoskeleton and meiosis were characterized as the 3 main gene annotations. Simultaneous up-regulation of extracellular matrix protein genes expression and down-regulation of matrix protein-ase genes was detected in gonads of androgen-treated females. At the same time the histological analysis of these gonads showed that they contain a predominant stroma of conjunctive tissue with fibroblast like cells. Matrix protein synthesis and the concomitant decrease in matrix proteinase activity have been well described as a characteristic fibrotic response of an excessive Transforming Growth Factor-beta (TGFβ) production [18,19]. Of special interest in that context is the up-regulation of transforming growth factor-β1 (tgfb1) in gonads of androgen-treated animals. In rat, TGFβ induced morphological changes in Leydig cells, accompanied by an increased secretion of fibronectin, laminin and collagen IV [20]. In fibroblasts treated with TGFβ1 a similar over-expression of genes Expression profiles of genes belonging to the Gene Ontology term "extracellular matrix" (ID 5578) Figure 6 Expression profiles of genes belonging to the Gene Ontology term "extracellular matrix" (ID 5578). Clones with similar gene names that belong to the same contig (ensemble of clones with overlapping sequences) but displayed very different expression profiles were excluded from the table (when the expression profiles were not too different the gene symbol are marked with an asterisk). Clones with similar gene names and belonging to different contigs were annotated as gene symbol_1 and gene symbol_2 as they could be considered as potentially duplicated genes or differential splicing forms. The 15 different contigs with homologies to zona pellucida protein homologs (zp2, zp3 and zp4 genes that are highly represented in clusters 1 to 5 have been removed from that table). KM: K-means cluster number. N: number of genes belonging to the GO term in the specified cluster.  Gene Name Gene Symbol Gene Ontology " Gene Ontology "extracellular extracellular matrix" (ID 5578) matrix" (ID 5578)

-1 +1
Expression profiles of some representative genes within the GO category: extracellular matrix (ID 5578) Figure 7 Expression profiles of some representative genes within the GO category: extracellular matrix (ID 5578).
Expression profile values were extracted from the DNA microarray dataset and normalized to the highest signal value observed among all samples for each gene. This highest signal value was arbitrarily set at 100 and the resulting values are designated the relative signal intensity for the studied gene at the indicated time points. Expression profiles of transforming growth factor-beta1, tgfb1 (AJ007836) were obtained by real-time RT-PCR. tgm2: transglutaminase 2, C polypeptide. col1a1: collagen, type I, alpha 1. col1a2: collagen, type I, alpha 2. fbln1: fibulin 1. sparc: secreted acidic cysteine rich glycoprotein. mmp1a: matrix metalloproteinase 1a. mmp13: matrix metalloproteinase 13.

Males Females Females Androgen
Androgen treated treated females females associated with matrix formation has been detected including many different matrix protein genes, like SPARC (Secreted Protein, Acidic and Rich in Cysteine), MGP (matrix Gla protein), and TGFβ1 itself [21], that we also detected as up-regulated in gonads following androgen treatment. It could then be hypothesized that this late androgen up-regulation of tgfb1 in trout gonads triggers a fibrotic response. Surprisingly, these effects are detected transiently and rather late after the application of the androgen treatment (but concomitantly with tgfb1 up-regulation). Whether this reflects a total dysregulation or an exacerbation of a testicular-specific event remains to be analyzed. However, extracellular matrix deposition is known as a major event for the testicular organization. For instance, LAMA5 (Laminin α5) has been characterized as a structural protein involved in the formation of the basement membrane of the testicular cords [22] and this pro-tein was found to be anti-correlated with Anti-Müllerian Hormone (AMH) [23]. In trout gonads, amh expression is not restored to male levels in androgen-treated females. This may produce a disrupted expression of some structural proteins, like lama5. In the same manner, sparc is highly up-regulated in androgen-treated females. In mouse Sparc gene expression has been identified in pre-Sertoli cells at the time of sex differentiation [24] and this protein has also been postulated to play a crucial role in both Leydig and Sertoli cells differentiation by affecting their morphology [25]. Structural proteins including matrix proteins are then of major importance for a complete and functional testicular differentiation and their up-regulation in trout following an androgen treatment inducing testicular transdifferentiation may be the consequence of a dysregulation of some major regulators of their synthesis like amh or tgfb1.
Expression profiles of some representative genes from cluster 13 Figure 8 Expression profiles of some representative genes from cluster 13. Expression profile values were extracted from the DNA microarray dataset and normalized to the highest signal value observed among all samples for each gene. This highest signal value was arbitrarily set at 100 and the resulting values are designated the relative signal intensity for the studied gene at the indicated time points. angiopoietin-like 7 (angptl7, tcad0004.b.24 and tcac0004.n.11), gonadal soma-derived growth factor (gsdf, tcac0002.i.22), family with sequence similarity 49, member b (fam49b, tcac0003.e.03), claudin 12 (cldn12, tcad0009.h.07). Each gene is depicted by its gene symbol and gene name (according to the zebrafish nomenclature). Symbols for mouse and human homolog genes are also given with their corresponding GenBank accession number. K-means cluster number (KM) is given for each gene with the correlation coefficient r to compare correlation between expression profiles measured by real-time RT-PCR (GenBank # in bold) and microarrays (clone numbers). Significant correlations were determined using the Pearson's correlation coefficient test (**: p < 0.01; ***: p < 0.05; NS: Not Significant).

gsdf
We also detected a high number of genes associated with cytoskeletal reorganization and muscle development that were up-regulated by the treatment. Some of them (e.g. cnn1, myh11, myl6, tagln) are even considered as characteristic smooth muscle markers. . It is therefore suggested that the masculinizing androgen treatment may induce the differentiation and subsequently a disturbed androgen-dependent proliferation of these peritubular myoid cells. These cells are also probably involved in the important extracellular matrix synthesis that occurs concomitantly with this differentiation.
In our experiment, the androgen treatment also induced a precocious spermatogenesis as revealed both by the histological analysis and by the increased expression levels of some genes involved in testicular meiosis. In fish, androgens and particularly, 11-oxygenated androgens, are strongly involved in spermatogenesis regulation [32] and they have been shown to directly induce spermatogenesis in vitro in some species [33]. Similarly, in mammals, three independent studies using Sertoli cell-specific AR-knockout mice (mice knockout for the androgen receptor, AR) demonstrated that the action of androgen is an absolute requirement for the completion of spermatogenesis, particularly in the process of meiosis [34][35][36].
Scatterplots of DNA microarray and real-time RT-PCR measurements for inha (inhibin alpha), amh (anti-müllerian hormone), star (Steroidogenic acute regulatory protein) and tf (transferrin) Figure 9 Scatterplots of DNA microarray and real-time RT-PCR measurements for inha (inhibin alpha), amh (antimüllerian hormone), star (Steroidogenic acute regulatory protein) and tf (transferrin). Each value is represented as the percentage of the highest value in each experiment (DNA microarray or real-time RT-PCR) in all experimental groups (red for females; blue for males; black for androgen-treated females). For each gene, the correlation coefficient R and colorized vectors of microarray and real-time RT-PCR values are given.

Conclusion
This study gives a first comprehensive survey of gene expression during androgen-induced masculinization in female rainbow trout. Our data provide supportive evidences that this treatment results in a marked dysregulation of gene expression levels when compared to natural testicular or ovarian differentiation. In our experimental condition the androgen treatment induces the complete down-regulation of female specific genes, but not the complete restoration of the male-specific gene expression patterns. Instead, some disturbed responses were characterized by an exacerbation of extracellular matrix synthesis and muscle type cell differentiation and proliferation (myoid cells) followed by a precocious meiosis of germ cells. All together, we suggest that androgen masculinization acts mainly through an early inhibition of female development rather than through a direct induction of testicular differentiation.

Animals and samplings
Research involving animal experimentation has been approved by the authors' institution (authorization no. . It conforms to principles for the use and care of laboratory animals and is in compliance with French and European regulations on animal welfare (European Convention for the Protection of Vertebrate Animals Used for Experimental and Other Scientific Purposes, ETS no. 123, January 1991). All-male and all-female rainbow trout populations were obtained at the INRA experimental fish farm (Sizun, France) as previously described [37]. Treatment with androgens (female treated Group 'F11β') was carried out at the onset of the first feeding [Day 0 = D0 at 55 days post-fertilization (55 dpf)], on an all-female population. The androgen, 11β-hydroxyandrostenedione (11βOH∆4, Sigma, St. Louis, MO, USA), was administered by adding it to the food (10 mg/kg food) during 3 months starting from the first feeding and this treatment has been shown to produce 100% sex-inversions [10]. In each group, 20 to 100 gonads were sampled and pooled in duplicates corresponding to the various stages of development: onset of the free swimming period after complete yolk resumption (Day 0 = D0), D0+7 days (D7), occurrence of oocyte meiosis (D12), beginning of ovarian lamellar structures development (D27), occurrence of previtellogenic oocytes (D60), D90 and D110. They were immediately frozen in liquid nitrogen and stored at -80°C until RNA extraction. Additional gonads were sampled at the same time-points for histological analysis, which was performed as previously described [38]. . Spotted DNA was then denatured and UV -cross-linked onto nylon filters. All DNA microarrays used in this study were made at the same time and under the same conditions. These trout microarrays contained 9,216 DNA spots representing 9,120 trout cDNA clones and a set of 96 controls. Among these cDNA clones, 7,584 were issued from a pooled-tissues library and 1,536 from a testis library [7]. Negative controls consisted of 80 spots of an Arabidopsis thaliana cytochrome c554 clone which is devoid of similarity with trout DNA sequences, 8 spots of poly(dA)80 and 8 spots of PCR reaction without template.

DNA microarray hybridizations
Microarrays were hybridized with two types of 33Plabeled probes. The first one was an oligonucleotide with a sequence common to all spotted PCR-products (vector hybridization) in order to determine the amount of target DNA accessible to hybridization in each spot. After stripping, a second hybridization was performed with complex probes made from 1 µg of retrotranscribed total RNA [40][41][42]. Protocols for probes preparation, hybridizations and washes are available online [43]. After stringent washes, arrays were exposed to phosphor-imaging plates and scanned with a FUJI BAS 5000 at 25 µm resolution. Hybridization signals were quantified using ArrayGauge software (Fuji Ltd, Tokyo, Japan).

Real-time RT-PCR
In order to validate and enrich the DNA microarray dataset, expression of 102 genes involved in early gonad development [15] was measured by real-time reverse transcription-polymerase chain reaction (RT-PCR). For cDNA synthesis, 1 µg of total RNA was denatured in the presence of random hexamers (0.5 µg) for 5 min at 70°C, and then chilled on ice. Reverse transcription (RT) was performed at 37°C for 1 h using M-MLV reverse transcriptase (Promega, Madison, WI, USA) as described by the manufacturer. Real-time PCR was carried out as previously described [15] using the iCycler iQTM (Bio-Rad, Hercules, CA, USA) and the SYBER Green PCR master Mix (Eurogentec, Seraing, Belgium). For each target gene, all the samples were analyzed on the same plate in the same PCR assay. PCR data were processed as previously described, each transcript level being normalized by division with the expression values of the constitutive elongation factor 1α (ef1a), which was used as an internal standard [15]. Data were then included in the microarray data matrix for clustering analysis (see next paragraph).

Data analysis
First, non-linear effects such as background, print-tip effects or saturation were corrected by LOWESS [44], using a channel by channel procedure [45]. Each array was individually normalized to the median profile of all arrays. We used the print-tip LOWESS version implemented in the statistical software package R [46]. Data were further corrected for the amount of spotted cDNA. This step is necessary as it has been shown that the signal intensity is proportional to the amount of probe on the surface of the array [40,47]. This effect can be observed both for glass and Nylon surfaces. This effect is corrected by the use of a reference in dual channel arrays, and by an independent measurement of the spotted amount of DNA probe in single channel arrays. On Nylon membranes, this effect is linear and can be corrected by dividing the signal by the amount of probe [40]. Briefly, sample signal intensity of each spot ("S") was divided ("S/V") by the corresponding signal intensity of the same spot obtained with the vector hybridization ("V"). To minimize experimental differences between different complex probe hybridizations, 'S/V' values from each hybridization were divided by the corresponding median value of 'S/V' (quantile normalization).
A triple filtering procedure was then applied to the microarray dataset. The first consisted of filtering background signals due to low amount of spotted DNA. When a "V" spot signal was too weak (vector signal < 3× vector local background), the data of the corresponding cDNA clone was discarded (missing data). The second filtering procedure was applied to eliminate non informative genes that were not measured (sample signal < 3× sample local back-ground) in more than 20% of the samples. Finally, genes exhibiting little variation (coefficient of variation < 0.1) across all arrays were excluded from the analysis [48,49]. After these three filtering steps, 2,372 genes were retained for further analysis.
All data (2372 microarray and 102 real-time RT-PCR gene expression profiles) were then log2-transformed and were analyzed by unsupervised and supervised clustering methods. Hierarchical clustering (Cluster program [50]) investigated the relationships between the genes and between the samples by using centroid linkage clustering with Pearson's uncentered correlation as similarity metric on data that were median-centered on genes. Gene clusters were distinguished using the non-hierarchical unsupervised learning k-means algorithm implemented in the Cluster program [50]. It was run on log2-transformed and gene median-centered data with a maximum cycles parameter of 100. The optimal minimal 'k' number of clusters, corresponding to the stability of the k-means clustering, was empirically set at 13. Indeed, with smaller k numbers, some clusters merged together whereas with greater k numbers, the size of some clusters decreased (less than 50 genes to truly empty clusters). Results (colorized matrix) of hierarchical and k-means clustering analyses were visualized using the Java TreeView program [51]. Functional annotation of genes was performed using Gene Ontology [52] and the GoMiner program [53]. Significance of over-or under-representation was calculated using Fisher's exact test at 0.05% risk.