# Meta-analysis of microarray data of rainbow trout fry gonad differentiation modulated by ethynylestradiol

Sophie Depiereux, Florence Le Gac, Bertrand de Meulder, Michael Pierre, Raphaël Helaers, Yann Guiguen, Patrick Kestemont, Eric Depiereux

► To cite this version:

Sophie Depiereux, Florence Le Gac, Bertrand de Meulder, Michael Pierre, Raphaël Helaers, et al. Meta-analysis of microarray data of rainbow trout fry gonad differentiation modulated by ethynylestradiol. PLoS ONE, 2015, 10 (9), pp.1-27. 10.1371/journal.pone.0135799. hal-01216754v3

HAL Id: hal-01216754

https://hal.science/hal-01216754v3

Submitted on 4 Nov 2015

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

RESEARCH ARTICLE

# Meta-Analysis of Microarray Data of Rainbow Trout Fry Gonad Differentiation Modulated by Ethynylestradiol

Sophie Depiereux¹*, Florence Le Gac³, Bertrand De Meulder², Michael Pierre², Raphaël Helaers², Yann Guiguen³, Patrick Kestemont¹, Eric Depiereux²

1 University of Namur (UNamur), Research Unit in Environmental and Evolutionary Biology (URBENARILIS), 61 rue de Bruxelles, B-5000, Namur, Belgium, 2 University of Namur (UNamur), Research Unit in Molecular Biology (URBM-NARILIS), 61 rue de Bruxelles, B-5000, Namur, Belgium, 3 Institut National de la Recherche Agronomique, INRA-LPGP, UPR1037, Campus de Beaulieu, 35042, Rennes, France

* sophie.depiereux@unamur.be

## Abstract

Sex differentiation in fish is a highly labile process easily reversed by the use of exogenous hormonal treatment and has led to environmental concerns since low doses of estrogenic molecules can adversely impact fish reproduction. The goal of this study was to identify pathways altered by treatment with ethynylestradiol (EE2) in developing fish and to find new target genes to be tested further for their possible role in male-to-female sex transdifferentiation. To this end, we have successfully adapted a previously developed bioinformatics workflow to a meta-analysis of two datasets studying sex reversal following exposure to EE2 in juvenile rainbow trout. The meta-analysis consisted of retrieving the intersection of the top gene lists generated for both datasets, performed at different levels of stringency. The intersecting gene lists, enriched in true positive differentially expressed genes (DEGs), were subjected to over-representation analysis (ORA) which allowed identifying several statistically significant enriched pathways altered by EE2 treatment and several new candidate pathways, such as progesterone-mediated oocyte maturation and PPAR signalling. Moreover, several relevant key genes potentially implicated in the early transdifferentiation process were selected. Altogether, the results show that EE2 has a great effect on gene expression in juvenile rainbow trout. The feminization process seems to result from the altered transcription of genes implicated in normal female gonad differentiation, resulting in expression similar to that observed in normal females (i.e. the repression of key testicular markers *cyp17a1*, *cyp11b*, *tbx1*), as well as from other genes (including transcription factors) that respond specifically to the EE2 treatment. The results also showed that the bioinformatics workflow can be applied to different types of microarray platforms and could be generalized to (eco)toxicogenomics studies for environmental risk assessment purposes.

<figure><img src="image_3.png" /></figure>

### OPEN ACCESS

**Citation:** Depiereux S, Le Gac F, De Meulder B, Pierre M, Helaers R, Guiguen Y, et al. (2015) Meta-Analysis of Microarray Data of Rainbow Trout Fry Gonad Differentiation Modulated by Ethynylestradiol. PLoS ONE 10(9): e0135799. doi:10.1371/journal.pone.0135799

**Editor:** Josep V. Planas, Universitat de Barcelona, SPAIN

**Received:** January 11, 2015

**Accepted:** July 27, 2015

**Published:** September 17, 2015

**Copyright:** © 2015 Depiereux et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** The raw data and a normalised expression file are available through the accession number GSE58519 at the Gene Expression Omnibus public data repository.

**Funding:** SD was funded by a grant from the 'Fond de la recherche scientifique', (FNRS): http://www.fnrs.be/. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. This work is supported by the FNRS (SD), FRIA (MP) and FRS-FNRS Télévie (BDM).
# Introduction

## 1.1. Transdifferentiation in fish exposed to exogenous hormones

In gonochoristic species, gonads (i.e. testis or ovary) develop from a germ and somatic cell primordium [1]. Sex determination can be genetically controlled (for instance in XX/XY systems, the expression of Y-linked genes leads to testicular development), and the subsequent differentiation processes are supported by molecular and endocrine factors, such as sex steroid hormones [2]. While these processes are strictly genetically controlled in mammals, where the expression of a main male gene, SRY (Sex-determining Region of Y chromosome), irrefutably leads to testis development, gonadal sex differentiation in fish is more labile, which means that their phenotype may be different from their sex genotype. Based on preliminary studies of steroid treatment influence on fish sex phenotype conducted by Yamamoto in the late 1960s, which highlighted this plasticity, exogenous steroid treatments were used in aquaculture to control the gender of fish stocks [3]. However, although the technique is now widely used, the precise mechanisms underlying this transdifferentiation process are still mostly unknown. In the current environmental context, pollution of rivers with endocrine disrupting compounds (EDCs; molecules that mimic natural endogenous hormones or interfere with the endogenous synthesis of hormones) has given rise to concerns about the sensitivity of fish to exogenous hormones, since very low doses (within ng to µg/L range) of these compounds are able to disrupt gonad morphology, leading to higher rates of intersexuality (gonads displaying male and female features simultaneously) and abnormalities. Several studies have shown that these disruptions could affect the fertility of fish and decrease population fitness [4,5]. Considering that EDCs do not necessarily mimic the endogenous physiological processes [6–8], a clear understanding of these mechanisms is crucial.

## 1.2. Microarray approaches

Modern genomic technologies, such as DNA microarray assays, are powerful tools to investigate initial responses to exogenous agents [9], and a considerable amount of microarray data has been and continues to be generated, for an increasing range of species [10]. The rapid expansion of public genomic databases and the development of bioinformatics tools have considerably enhanced our ability to infer biological meaning from large gene lists. This technique has proven to be efficient for the study of predicted biological pathways and regulatory networks, with wide ranging applications from human clinical biomarker discovery to drug development [11]. More recently, ecotoxicological studies aimed at determining the molecular mode of action (MOA) of environmental pollutants and identifying their biomarkers as indicators of exposure have been published [11,12]. However, certain practices commonly used for human and rodent studies are not fully adapted to the use of microarrays in aquatic ecotoxicology studies. Indeed, the inference of biological meaning from large gene lists relies on the quality of the genomic resources for the species investigated, and especially the annotation of the oligonucleotides spotted and the body of pre-existing genomic knowledge (ontologies, pathways, ...). Though the genomes of a growing number of fish species have been sequenced, the available specific genomic resources are limited, hampering data-mining efforts.

Inherent issues related to microarray data analysis must also be taken into account. Even at the lowest significance threshold, thousands of statistical tests inevitably contaminate the differentially expressed gene (DEG) lists with a high number of false positives and false negatives [13]. Many efforts were made in the 2000s to develop new statistical methods [14] improving the gene variance estimation [15–18]. However, the problem still remains fundamentally unsolved, the combination of existing methods generates millions of possibilities [19] and the
appropriate statistical analysis strategy is not handled by all experimental laboratories. Therefore, gene top lists obtained after statistical analysis depend on the methodology used and can dramatically differ between research teams.

## 1.3. Re-analysis and Meta-analysis

In view of the considerations listed above, most microarray datasets remain underexploited. However, the fast evolution of (i) available data in public genomic databases, (ii) statistical analysis performances and (iii) data-mining bioinformatics tools has made it possible to gain new information through the re-analysis of previously published datasets. Moreover, the meta-analysis of several datasets is a way to yet further enhance analytical performance [19–21]. In a previous work, we successfully adapted a microarray analysis workflow validated on human cancer data [19,22][23] to a toxicogenomic study using rainbow trout as model organism [24]. We measured gene expression in juvenile male fish exposed chronically to several doses of the potent xenoestrogen ethynylestradiol (EE2). This allowed highlighting several pathways and ontologies related to EE2 exposure and to propose several potential specific and/or sensitive biomarkers of the intersex stage in juvenile rainbow trout. Here, we propose to apply this methodology to the re-analysis of a previously published dataset studying the effect of EE2 on a time-course experiment using juvenile male rainbow trout. We postulate that the meta-analysis of both datasets will allow highlighting the pathways and key genes potentially implicated in early male-to-female gonadal transdifferentiation.

## Methods

Two datasets were used in the meta-analysis. The microarray analysis workflow is explained below and summarized in Table 1.

### 2.1. Dataset 1 [24]

In this experiment, all-male rainbow trout (*Oncorhynchus mykiss*) fry (INRA experimental fishfarm, Monts d'Arrée-PEIMA, Sizun, France) were exposed from the onset of first feeding [Day 0 = D0 at 60 days post-fertilization (dpf)] to 136 dpf to 5 nominal concentrations of 17α-ethynylestradiol (purity ≥ 98%, Sigma-Aldrich, Germany): 0 (solvent control), 0.01 µg/L, 0.1 µg/L, 1 µg/L and 10 µg/L, with 3 tanks per condition. The actual EE2 concentrations were measured in each tank at 6 time points using the Quantitative Ethynylestradiol Enzyme Immunoassay (EIA) Kit (Marloie, Belgium) according to the manufacturer's instructions. Mean concentrations of EE2 ± SD were 0.08 ± 0.06 µg/l; 1.62 ± 1.74 µg/l and 9.88 ± 5.06 µg/l. The 0.01 µg/l EE2 concentration was under the detection limit (set at 0.02 µg/l). At the end of the exposure, the gonads were collected from all fish, immediately frozen in liquid nitrogen and stored at -80°C until RNA extraction. The gonads from ten fish were pooled to reach enough material for further analyses. Gene expression measurement was performed on 6 replicates per condition (2 pools of gonads per tank) on a one-color 8x60K oligonucleotide array (Agilent Technologies) designed by the INRA-LPGP microarray platform (Beaulieu Campus, Rennes, France). Expression profiles obtained from the DNA microarray dataset were validated with the measurement of several genes by real-time reverse transcription-polymerase chain reaction (RT-PCR) as described previously [25].

### 2.2. Dataset 2 [26]

In a previous experiment (Guiguen et al., unpublished data), transcript expression was recorded over a time-course analysis on gonads sampled following feminization of genetically
<table><caption>Table 1. Summary of the meta-analysis workflow.</caption><thead><tr><th></th><th>Dataset 1</th><th>Dataset 2</th></tr></thead><tbody><tr><td>Species</td><td>Oncorhynchus mykiss—Genetically all-male population</td><td>Oncorhynchus mykiss—Genetically all-male and all-female populations</td></tr><tr><td>Treatment conditions</td><td>Balneation—EE2-5 concentrations : 0 (control); 0.01 ; 0.1 ; 1 and 10 µg/L</td><td>Feeding—EE2-1 concentration : 20 mg/kg food</td></tr><tr><td>Developmental stage</td><td>Day 0 (D0) at 60 dpf (first feeding) up to 136 dpf</td><td>D0 at 55 dpf (first feeding) up to 111 dpf</td></tr><tr><td>Sampling</td><td>One sampling at 136 dpf</td><td>7 samplings at different developmental stages: D0, D0+7 days (D7), D16, D30, D63, D91, and D111</td></tr><tr><td>Sample size</td><td>Pools of gonads from ten fish per sample. 6 replicates per condition</td><td>Pools from 20 to 100 fish per sample. 2 replicates per condition</td></tr><tr><td>Platform used</td><td>8x60K oligonucleotide array (Agilent technology)</td><td>9,023 trout cDNA platform</td></tr><tr><td>Statistical analysis</td><td>One-way ANOVA, with EE2 concentration as criterion</td><td>Two-way ANOVA with time-course (D0 to D111) and sex as criteria</td></tr><tr><td>Scheffe's contrasts</td><td>4 physiologically relevant contrasts defined as control vs. different EE2 concentrations: (1) Control vs. 0.01 µg/L—CT1; (2) Control vs. 0.1 µg/L—CT2; (3) Control vs. 1 µg/L—CT3; (4) Control vs. 10 µg/L—CT4</td><td>19 physiologically relevant contrasts defined as Male-EE2 (= male test = MT) versus male-control (MC) at the same sampling time: (1) MT-MC D7, (2) MT-MC D16, (3) MT-MC D30, (4) MT-MC D63, (5) MT-MC D91, (6) MT-MC D111; Male-EE2 (= male test = MT) versus female-control (FC) at the same sampling time: (7) MT-FC D7, (8) MT-FC D16, (9) MT-FC D30, (10) MT-FC D63, (11) MT-FC D91, (12) MT-FC D111; Male-control (MC) versus female-control (FC) at the same sampling times: (13) MC-FC D0, (14) MC-FC D7, (15) MC-FC D16, (16) MC-FC D30, (17) MC-FC D63, (18) MC-FC D91, (19) MC-FC D111.</td></tr><tr><td>Intersection I<sup>1</sup></td><td>DEGs from the four contrasts in dataset 1: CT1; CT2; CT3; CT4</td><td>DEGs from contrasts including EE2-treated fish versus male and female control fish: MT-MC D7 to D111; MT-FC D7 to D111.</td></tr><tr><td>Intersection II<sup>2</sup></td><td>DEGs from the four contrasts in dataset 1: CT1; CT2; CT3; CT4</td><td>DEGs from contrasts including EE2-treated fish versus male control fish: MT-MC D7 to D111.</td></tr><tr><td>Intersection III<sup>3</sup></td><td>CT1; CT2. These two contrasts were previously correlated with gonad intersex feature<sup>4</sup>.</td><td>MT-MC D7; MT-MC D16; MT-MC D30. These contrasts cover early developmental stages corresponding to differentiating gonads, with first morphological differences appearing between male and female gonads.</td></tr><tr><td>Final intersection<sup>5</sup></td><td>DEG selected in intersection III : (i) common to pathway/ Ontologies selected with intersections I and II (ii) selected by vizualization techniques (Volcanoplots and kinetic profiles)</td><td>DEG selected in intersection III : (i) common to pathway/ Ontologies selected with intersections I and II (ii) selected by vizualization techniques (Volcanoplots and kinetic profiles)</td></tr></tbody></table>

¹ Larger gene list retrieving genes common to both datasets. This list was used for pathway analysis. We postulated that it represents mechanistic signatures of EE2 treatment on rainbow trout juveniles.

² Shorter gene list retrieving genes common to both datasets. This list was used for ontology analysis. We postulated that it represents mechanistic signatures of EE2 treatment on rainbow trout juveniles.

³ Shortest gene list retrieving genes common to both datasets. We postulated that it focuses on genes potentially involved in early responses to EE2 treatment. This third group was selected based on union intersections, which represents DEGs common to at least one condition in dataset 1 and one condition in dataset 2 (for example DEG in CT1 and MT-MC D7).

⁴ [25]

⁵ Limited number of genes considered true positives after all filtering steps of the analysis workflow. We postulated that these genes are potentially implicated in male-to-female transdifferentiation processes. Gene expression profiles of this final gene list were compared with gene expression between normal (control) male and female juvenile rainbow trout.

doi:10.1371/journal.pone.0135799.t001

all-males population (M-control) with estrogens (M-EE2) or masculinization of genetically all-females population (F-control) with either androgen (F-11β) or an anti-aromatase (F-ATD). The present analysis is limited to estrogen-treated fish. Feminizing treatment of all-male populations (PEIMA, Sizun, France) with 17α-ethynylestradiol (EE2, 20 mg/kg of food, Sigma,
St Louis, MO, USA) was applied for 2 months starting from the first feeding. In each group, 20 to 100 gonads were sampled and pooled in duplicate corresponding to the various stages of development: onset of the free swimming period after complete yolk resorption [Day 0 = D0 at=55 dpf], D0+7 days (D7), D16, D30, D63, D91, and D111. Gonads were immediately frozen in liquid nitrogen and stored at -80°C until RNA extraction. Gene expression analyses were carried out using home-made one-color Nylon DNA microarrays obtained as previously described [6] and containing 9,023 trout cDNA clones (TROUT-AGENAE). Expression profiles obtained from the DNA microarray dataset were validated with the measurement of several genes by real-time reverse transcription-polymerase chain reaction (RT-PCR).

## 2.3. Microarray analysis workflow

The initial microarray analysis workflow [24] was adapted to handle another type of microarray platform (dataset 2) and to fit with the meta-analysis scheme and work hypothesis. Briefly, this workflow includes (i) the analysis of pre-treated data by (ii) ANOVA analysis followed by Scheffe's post hoc pairwise comparisons. Thereafter, (iii) biologically relevant sets of contrasts and gene list intersections between the two datasets were selected at different levels of significance. The remaining gene lists, supposed to be enriched in true positive DEGs, were (iv) submitted to over-representation analysis methods (ORA). (v) A set of most relevant DEG was selected using volcanoplots, membership of enriched pathways and/or GOterms and kinetic profiles in both analyses. (vi) Those genes were compared with normal expressions profiles between control rainbow trout male and female developing fry.

## 2.4. Statistical analysis

Data analyses were performed using the R statistical software version 2.15.3 available on the R-Project website (http://cran.r-project.org) and a set of packages available in the Bioconductor repository (http://www.bioconductor.org). This methodology used for the first dataset [24] was also applied to re-analyse the dataset 2. Therefore, brief descriptions and specificities related to this dataset analysis are provided below. See supporting information for detailed scripts.

2.4.1. **Pre-processing and normalisation procedures.** The dataset was configured in numeric matrix. The expression values were first submitted to a quantile-quantile normalization using the normalize.quantile function from the preprocessCore package (S1 Dataset).

2.4.2. ANOVA. A two-way ANOVA was performed, with the time-course (D0 to D 111) and sex as criteria, using the lmFit function of package Limma, the only statistical method able to perform ANOVA with more than one factor level (here the 7 developmental stages tested) on microarray data. Its predictive power has been well ranged in our previous benchmark of method performances [27].

2.4.3. **Post-hoc comparisons.** A post-hoc evaluation was performed using Scheffe's method with the makeContrasts() and eBayes() functions provided in the Limma package (88 contrasts comparing all experimental conditions). The results were summarized using decideT-ests() and summary() (number of significant genes positively or negatively regulated for each evaluated contrast). The Benjamini-Hochberg procedure was used to adjust the p-values (correction for multiple testing) [28]. Considering the very large number of genes per contrast group, we later focused the analysis of dataset 1 on four groups that were the most relevant with regards to the biological context, comparing reference samples with samples treated with EE2 at different concentrations. The contrast labels were defined as control vs. different EE2 concentrations—[EE2]: (i) CT1 = Control vs. [0.01 µg/L]; (ii) CT2 = Control vs. [0.1 µg/L]; (iii) CT3 = Control vs. [1 µg/L] and (iv) CT4 = Control vs. [10 µg/L].
For the 2nd dataset, a differential gene expression analysis was performed using the two-way ANOVA, considering the developmental stages (D0 to D111) and sex (male and female) as criteria. Eighty-eight contrasts were performed comparing all the experimental conditions. We postulated that gene expression comparisons make sense between fish at the same developmental stage. To this end, we defined 19 physiologically relevant contrasts as follows: (i) we selected the 7 male-EE2 (= male test = MT) versus male-control (MC) contrasts at the same sampling time (MT-MC D7, MT-MC D16, MT-MC D30, MT-MC D63, MT-MC D91, MT-MC D111); (ii) the 7 male-EE2 versus female-control (FC) contrasts at the same sampling times (MT-FC D7, MT-FC D16, MT-FC D30, MT-FC D63, MT-FC D91, MT-FC D111), and (iii) the 8 male-control (MC) versus female-control (FC) contrasts at the same sampling times (MC-FC D0, MC-FC D7, MC-FC D16, MC-FC D30, MC-FC D63, MC-FC D91, MC-FC D111). For each contrast, R software topTable() function was used to extract top-ranked genes and relevant statistics (logFC, Average Expression, t statistic, p-value, adjusted p-value, B statistic).

**2.4.4. Meta-analysis.** In theory, meta-analysis of microarray data consists of analyzing several related datasets at once [22]. The statistical analysis is then made on merged datasets. In our analysis, this scheme was impossible as the platforms used were heterogeneous. Here, the meta-analysis consists of retrieving DEGs common to both analyses. The use of intersections between top gene lists has proven to be effective in filtering relevant genes [19,22,24]. Three intersections (Intersection I, II and III) were defined for further analyses (Table 1).

Intersection I has been proposed to retrieve genes common between the two DEGs lists from the two datasets. We postulated that this represents mechanistic signatures of EE2 treatment on rainbow trout male juveniles. This group therefore contained the intersection between all DEGs in the four contrasts in dataset 1 (CT1 to CT4) and all DEGs in the twelve contrasts including EE2-treated fish in dataset 2 (MT-MC D7 to D111, and MT-FC D7 to D111). This group contained 1,535 DEGs and was further used in the pathway overrepresentation analysis (see result section 3.2.1.).

To further focus on genes differentially expressed in male treated fish, a smaller group of DEGs (Intersection II) was obtained at the intersection between all DEGs in the four contrasts in dataset 1 (CT1 to CT4) and all DEGs in the six contrasts including EE2-treated fish versus male control fish in dataset 2 (MTMC D7 to D111). This group contained 438 DEGs and was further used in the ontology overrepresentation analysis (see result section 3.2.2.).

To further investigate the trans-differentiation process, we aimed to focus on genes potentially involved in the early response to EE2 treatment (Intersection III). The statistical and bioinformatics approaches applied here were aimed at retrieving true positives among contaminated large gene lists. For this, we focused on the DEGs related to the lower doses used in dataset 1 (corresponding to contrasts CT1 and CT2). We have reported previously that these conditions are correlated with gonad intersex feature development in rainbow trout male juveniles exposed to EE2 [24]. In dataset 2, the interesting contrasts were those at the early developmental stages, 7, 16 and 30 days after treatment. In developing rainbow trout, this period corresponds to differentiating gonads, with the first morphological differences appearing between male and female gonads. Day 12 corresponds to the occurrence of oocyte meiosis and D27 to the beginning of ovarian lamellar structure development in females. Gametogenesis starts from 60 dpf [29]. The third group was composed of the union intersection between contrasts CT1 and CT2 in dataset 1 and MT-MC D7, MT-MC D16 and MT-MC D30 in dataset 2. The union intersection included DEGs common to at least one condition in dataset 1 and one condition in dataset 2 (for example DEG in CT1 and MT-MC D7). This retrieved DEGs that could be false negatives in other contrasts considered. To further filter these DEGs, (i) we confronted this list to the pathways and ontologies retrieved in previous steps. Thereafter, (ii) we
used visualization techniques. First (a), we submitted them to a Volcanoplot analysis, then (b) we looked at their kinetic profiles in both analyses. Finally, the intersection of gene lists retrieved by each approach was selected. (iii) The final gene list was then compared with DEGs between normal male and female differentiation.

**2.4.5. Annotations.** Considering that the recently available genome of rainbow trout has not yet a complete and accurate functional annotation, no consensus exist for the annotation of the Agilent 60K array and for the AGENAE_TROUTGeneric2_9216 platforms. Bioinformatics tools have different requirements concerning the gene IDs to enter in the program, but commonly they need homogenous identifiers and no redundancy. In addition, the annotation of oligonucleotides improves daily and can change substantially within several months. Therefore, it is important to use the most accurate annotation at the time of analysis. To be able to use DAVID tool (see below), we used ensembl gene ID for Danio rerio (ENSDARG), the only fish species with pathways available on Kegg [30]. For EASE, as no entry exists for ensembl gene ID, Gene Symbol identifiers were used, because the most accurate. We have updated the latest version of the AGENAE program (January 2013) (Analysis of Breeding Animals' Genome; http://www.agenae.fr/programme_agenae/) for these two annotations, using UniProt (http://www.uniprot.org) or e!Ensembl (http://www.ensembl.org) databases.

**2.4.6. Visualizations.** Several bioinformatics tools were used to gain insights of the dataset, at different levels of organization (gene expression, pathways, ontologies) depending on the analysis stringency of the group of DEG tested. The “Intersection I” group of DEG (see results) was submitted to the web tool DAVID (Database for annotation, Visualization and Integrated Discovery; http://david.abcc.ncifcrf.gov/home.jsp) [31,32], to retrieve pathway maps from KEGG following the method applied by Pierre et al, 2010 [19], and with the genes on the 60K rainbow trout Agilent array set up as the ‘background population’. The “Intersection II” group was submitted to the EASE software (version 2.0)(Expression Analysis Systematic Explorer) available on the DAVID website [33]. EASE is a desktop version of DAVID allowing more freedom toward the identifiers entering the program, which suits our heterogeneous Swissprot annotation (several species orthologs). Gene symbol was selected in the Input Genes setting. Then we choose the “Find over represented gene categories” function. In the “Population file” window we enter the entire array annotation without redundancy. Then we set the “Run Basic Analysis” function. Finally, Volcanoplots were made to visualize the DEG of interests using the plot function in the R software. A volcanoplot is a graph sorting the results by p-value (Y-axis),-log fold change (X-axis) and over/under-expression. In such graphs, the most interesting genes are usually located in the upper left and right corner of the plot, depicting genes with low p-values and high fold changes.

# Results

## 3.1. ANOVA and Scheffe's contrasts

Fig 1 summarizes the microarray meta-analysis workflow and shows an overview of the associated results. The tables of contrasts include redundant genes.

The results for the 10 contrasts, presented by number of differentially expressed genes without redundancy but before annotation, are shown in Fig 1 (left-hand table). The cut-off for DEG inference was set to 0.05 for the adjusted p-values (Benjamini-Hochberg). The four contrasts previously defined in section 2.4.3 for the 1st dataset represent 12,742 non-redundant DEGs (sum of the orange cells in the left-hand table in Fig 1, with subtraction of the redundant genes, namely genes present in several contrasts). For the 2nd dataset, the contrasts represent 737 (in orange), 1,878 (in green) and 2,305 non-redundant DEGs (in blue) respectively (Fig 1).
<figure><img src="image_4.png" /><figcaption>**Fig 1.** Flowchart of the workflow of the microarray meta-analysis methodology and associated results. This flowchart summarizes the bioinformatics workflow used for the meta-analysis of two datasets studying sex reversal following exposure to EE2 in juvenile rainbow trout. Both datasets (dataset 1 and dataset 2, pre-treated data) were subjected to ANOVA analyses followed by Scheffe's post hoc pairwise comparisons. Thereafter, biologically relevant sets of contrasts (coloured and in bold in the flowchart) were selected (differentially expressed genes (DEGs) between control and EE2-treated fish in dataset 1; DEGs between control and EE2-treated fish at the same developmental stage in dataset 2). The meta-analysis consisted of retrieving the intersection of the top DEGs lists generated for both datasets, performed at different levels of stringency. This allowed creating different groups of DEGs (Intersections I, II and III) for ORA analyses downstream, given that pathway or Ontologies research have different constraint concerning the group of genes to enter in the program. It is worth noting that each group of DEGs has different biological meaning described further in the text in the results section. Intersection I (low stringency) contained the intersection between all DEGs in the four contrasts in dataset 1 (CT1 to CT4) and all DEGs in the twelve contrasts including EE2-treated fish in dataset 2 (MT-MC D7 to D111, and MT-FC D7 to D111). This group contained 1,535 DEGs (1,171 DEGs having a functional annotation with Danio rerio Ensembl! identifiers-ENSARG) and was further used in the pathway overrepresentation analysis. A smaller group of DEGs (Intersection II) was obtained by retrieving the intersection between all DEGs in the four contrasts in dataset 1 (CT1 to CT4) and all DEGs in the six contrasts including EE2-treated fish versus male control fish in dataset 2 (MTMC D7 to D111). This group contained 438 DEGs (245 DEGs having a functional annotation with Swissprot identifiers) and was further used in the ontology overrepresentation analysis. The third group (Intersection III—highest stringency— 103 DEGs) was composed of the union intersection between contrasts CT1 and CT2 in dataset 1 and MT-MC D7, MT-MC D16 and MT-MC D30 in dataset 2. From this last list of DEGs, a set of most relevant DEGs was selected using volcanoplots, membership of enriched pathways and/or GOterms and kinetic profiles in both analyses. The 20 remaining genes, considered as potentially implicated in the transdifferentiation processes, were selected for further analyses.</figcaption></figure>

doi:10.1371/journal.pone.0135799.g001

## 3.2. Meta-analysis on selected of groups of DEGs (Intersection I, II and III)

### 3.2.1. Intersection I: pathway over-representation analysis.
The list of Ensembl gene IDs for *Danio rerio* corresponding to the “Intersection I” DEG group was used for pathway
<table><caption>**Table 2. List of significantly enriched pathways.** The list of significantly enriched pathways (P-Value/EASE score < 0.05) was obtained by using the DAVID web tool for the Intersection I group in our meta-analysis (1,535 gene IDs restricted to a set of 1,171 well annotated differentially expressed genes recognized by ENSDARG identifiers in DAVID). This group represents the genes differentially expressed following steroid modulation of immature gonads in rainbow trout by EE2, under several types of treatment. The Count column gives the number of genes of our genelist belonging to the pathway. The EASE--Score column represents the P-Value calculated by DAVID. The Benjamini column represents the Benjamini-corrected P-values (correction for multiple testing).</caption><thead><tr><td>Term</td><td>Count</td><td>EASE-score</td><td>Benjamini</td></tr></thead><tbody><tr><td>Cell cycle</td><td>34</td><td>1.1<sup>E-6</sup></td><td>1.3<sup>E-4</sup></td></tr><tr><td>DNA replication</td><td>15</td><td>2.9<sup>E-5</sup></td><td>1.8<sup>E-3</sup></td></tr><tr><td>Oxidative phosphorylation</td><td>29</td><td>6<sup>E-5</sup></td><td>2.5<sup>E-3</sup></td></tr><tr><td>Pyrimidine metabolism</td><td>20</td><td>3.7<sup>E-3</sup></td><td>1.1<sup>E-1</sup></td></tr><tr><td>Nucleotide excision repair</td><td>12</td><td>8.2<sup>E-3</sup></td><td>1.8<sup>E-1</sup></td></tr><tr><td>Mismatch repair</td><td>8</td><td>1.1<sup>E-2</sup></td><td>2<sup>E-1</sup></td></tr><tr><td>Fatty acid elongation in mitochondria</td><td>5</td><td>1.5<sup>E-2</sup></td><td>2.4<sup>E-1</sup></td></tr><tr><td>Ribosome</td><td>18</td><td>1.6<sup>E-2</sup></td><td>2.2<sup>E-1</sup></td></tr><tr><td>PPAR signaling pathway</td><td>12</td><td>3.4<sup>E-2</sup></td><td>3.8<sup>E-1</sup></td></tr><tr><td>RNA degradation</td><td>12</td><td>4.4<sup>E-2</sup></td><td>4.2<sup>E-1</sup></td></tr><tr><td>Progesterone-mediated oocyte maturation</td><td>16</td><td>4.7<sup>E-2</sup></td><td>4.2<sup>E-1</sup></td></tr><tr><td>Valine, leucine and isoleucine degradation</td><td>10</td><td>4.8<sup>E-2</sup></td><td>3.9<sup>E-1</sup></td></tr></tbody></table>

doi:10.1371/journal.pone.0135799.t002

over-representation analysis. Of the 1,535 DEGs in intersection I, 1,171 were accurately annotated with an ENSDARG ID and recognized by the DAVID web tool. The background was made up of 16,977 ensembl IDs (corresponding to the total number of identifiers on the Agilent 60K platform). A total of 79 pathways containing at least 2 genes were identified from the remaining 1,171 genes (S1 Results). Among these pathways, 12 were detected with a significant enrichment in the list of DEGs submitted, with an EASE score below the threshold of 0.05 (Table 2).

Only three pathways (Cell cycle, DNA replication and oxidative phosphorylation) could be detected when the threshold was set to adjusted p-values, after the highly conservative Benjamini multiple testing corrections (see Table 2). Most of these enriched pathways are involved in Genetic Information Processing, including replication and repair (DNA replication, Nucleotide excision repair, Mismatch repair) but also translation (Ribosome) and degradation (RNA degradation). Several pathways belonged to Metabolism, with energy (Oxidative phosphorylation), nucleotide (pyrimidine metabolism), lipid (fatty acid elongation) and amino acid (valine, leucine and isoleucine degradation) groups represented. Of particular interest were those pathways related to cellular processes: the “Cell cycle” was the most significantly enriched set of genes, with the highest number of highlighted genes (34), but the “Endocrine system” was also identified with PPAR signalling pathway and Progesterone mediated oocyte maturation. Gene enrichment is illustrated with red stars corresponding to the DEGs retrieved in our analysis (Fig 2).

**3.2.2. Intersection II: ontology over-representation analysis.** The list of Ensembl gene IDs for *Danio rerio* corresponding to the “Intersection II” DEG group was used for ontology over-representation analysis. Of the 438 DEGs of this intersection, 245 could be accurately annotated. Interestingly, 58 terms were significantly enriched (EASE score less than 0.05). Only GO Terms with highly significant EASE scores (p-values less than 0.01) were considered. Table 3 provides a summary of the results obtained. An exhaustive presentation of the relationships between all genes and each category is provided in S2 Results. Genes that belonged to the “Intersection II” set could be summarized in four broad terms with high significance: (i) the most significant terms belonged to “Reproduction” (more precisely, sexual reproduction and
<figure><img src="image_5.png" /><figcaption>**Fig 2. Cell cycle.** In our study, the Cell Cycle is an example of a significantly enriched pathway resulting from the DAVID analysis. Red stars indicate genes belonging to the list of DEG submitted to the program. The background was made up of 16,977 Ensembl IDs (corresponding to the total number of identifiers on the Agilent 60K platform).</figcaption></figure>

doi:10.1371/journal.pone.0135799.g002

gametogenesis). Interestingly, several genes were related to male gamete maturation and spermatogenesis; (ii) several enriched terms were related to the cellular process (cell cycle, cell proliferation, ...) and more precisely to mitosis (mitotic cell cycle, M phase, S phase, ...); (iii) several terms belonged to the metabolic process (oxidoreductase activity, DNA metabolism, ...); (iv) DNA replication was overrepresented.

### 3.2.3. Intersection III: membership of pathway/ontologies, volcanoplots and kinetics.

Intersection III group contained 103 DEGs. They were compared with genes belonging to enriched Pathways and/or GO Terms highlighted earlier. Among the 103 DEGs, 44 belong to at least one enriched pathway and/or GO Term. Several of these listed genes appeared to be ubiquitous. For example, *RFC3* can be found in the DNA replication pathway as well as GO Terms related to the Cell cycle, Cell proliferation, and the Chromosome. *Ropn1* appeared in Spermatogenesis, Cell proliferation and the mitochondrial membrane.

Two visualization techniques were then used to select genes of interest. First we selected the most significant genes in terms of Fold change and P-Value by Volcanoplot visualization over
<table><caption>**Table 3. List of significantly enriched GO terms.** This table lists the significant GO terms (P-Value/EASE score < 0.01) retrieved by the EASE software for the Intersection II group in this analysis (245 well annotated DEGs recognized by Gene Symbol identifiers in EASE among 438 DEGs IDs). This group represents the DEGs following a steroid modulation of immature gonads in rainbow trout by EE2, under several types of treatment. The EASE-Score column represents the P-Value calculated by DAVID.</caption><thead><tr><td>GO ID</td><td>Term</td><td>EASE score (P-Value)</td></tr></thead><tbody><tr><td>GO:0019953</td><td>Sexual reproduction</td><td>3.2<sup>E-07</sup></td></tr><tr><td>GO:0000003</td><td>Reproduction</td><td>3.2<sup>E-07</sup></td></tr><tr><td>GO:0007276</td><td>Gametogenesis</td><td>1.7<sup>E-06</sup></td></tr><tr><td>GO:0008283</td><td>Cell proliferation</td><td>5.9<sup>E-05</sup></td></tr><tr><td>GO:0007049</td><td>Cell cycle</td><td>7.2<sup>E-05</sup></td></tr><tr><td>GO:0000278</td><td>Mitotic cell cycle</td><td>1.4<sup>E-04</sup></td></tr><tr><td>GO:0000074</td><td>Regulation of cell cycle</td><td>1.0<sup>E-03</sup></td></tr><tr><td>GO:0006261</td><td>DNA dependent DNA replication</td><td>1.0<sup>E-03</sup></td></tr><tr><td>GO:0005743</td><td>Mitochondrial inner membrane</td><td>1.0<sup>E-03</sup></td></tr><tr><td>GO:0009117</td><td>Nucleotide metabolism</td><td>1.3<sup>E-03</sup></td></tr><tr><td>GO:0019866</td><td>Inner membrane</td><td>2.0<sup>E-03</sup></td></tr><tr><td>GO:0005737</td><td>Cytoplasm</td><td>2.3<sup>E-03</sup></td></tr><tr><td>GO:0000279</td><td>M phase</td><td>3.1<sup>E-03</sup></td></tr><tr><td>GO:0005694</td><td>Chromosome</td><td>3.2<sup>E-03</sup></td></tr><tr><td>GO:0016491</td><td>Oxidoreductase activity</td><td>3.5<sup>E-03</sup></td></tr><tr><td>GO:0005740</td><td>Mitochondrial membrane</td><td>3.9<sup>E-03</sup></td></tr><tr><td>GO:0006220</td><td>Pyrimidine nucleotide metabolism</td><td>4.5<sup>E-03</sup></td></tr><tr><td>GO:0006260</td><td>DNA replication</td><td>5.1<sup>E-03</sup></td></tr><tr><td>GO:0006259</td><td>DNA metabolism</td><td>5.4<sup>E-03</sup></td></tr><tr><td>GO:0000084</td><td>S phase of mitotic cell cycle</td><td>5.8<sup>E-03</sup></td></tr><tr><td>GO:0048232</td><td>Male gamete generation</td><td>8.4<sup>E-03</sup></td></tr><tr><td>GO:0007283</td><td>Spermatogenesis</td><td>8.4<sup>E-03</sup></td></tr><tr><td>GO:0019201</td><td>Nucleobase, nucleoside, nucleotide kinase activity</td><td>9.1<sup>E-03</sup></td></tr></tbody></table>

doi:10.1371/journal.pone.0135799.t003

all the contrasts considered. Second, the kinetic profiles of remaining genes were visualized in both datasets.

Volcanoplots are presented in Fig 3. These graphs show the comparison between the five contrasts considered for the 103 DEGs of the Intersection III group. DEGs were selected from these graphs as follows (i) Location in the right or left corners of the graphs (ii) Membership of at least two contrasts between dataset 1 and dataset 2, and (iii) Consistency of fold change type (positive or negative) between both involved contrasts. In total, 60 DEGs were selected by this approach. Among this list, 23 DEGs belonged to enriched Pathways/GO Terms. For the sake of clarity, only cyp11b is represented on graphs, which was interestingly found at a very similar position in all the contrasts considered (i.e. left corner of the graph).

The next step was to look at the kinetics profiles of the 60 DEGs selected by the Volcanoplot analysis from the two whole datasets considered (CT1 to CT4 in dataset 1 and MTMC D7 to D111 in dataset 2). The kinetics profiles for the 20 genes in dataset 1 and dataset 2 are given in Fig 4A and 4B and Fig 5A–5T, respectively. Selection was based on deleting genes whose maximum expression was under 1.5 in overexpressed genes or over -1.5 in underexpressed genes. This retrieved 42 DEGs in dataset 1 and 22 DEGs in dataset 2. Through the identification of DEGs common to these two lists, we found 20 DEGs of particular interest (Table 4). Among these genes, 7 belonged to significantly enriched Pathways/GO Terms. It is worth noting that, among the 20 DEGs selected, there were several exceptions. TEKT1 and Kindiss have opposing fold changes, as they are overexpressed in dataset 1 and underexpressed in dataset 2. Moreover,
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Fig 3. A-E: Volcano plots.** These graphs represent the 103 DEGs belonging to the 'Intersection 3' group of DEGs selected through the intersection between low doses (contrasts CT1 and CT2) in the Dataset 1 and early stages (contrasts MT-MC D7, D16 and D30) in dataset 2. For sake of clarity, red points</figcaption></figure>
represent DEGs. As an example, CYP11B is represented by its label. The green bars represent the log2(+/-2) fold change and the blue bar represents a p-value threshold of 0.05. (A) Selected DEGs in contrast MTMC D7 in dataset 2. (B) Selected DEGs in contrast MTMC D16 in dataset 2. (C) Selected DEGs in contrast MTMC D30 in dataset 2. (D) Selected DEGs in contrast CT1 in dataset 1. (E) Selected DEGs in contrast CT2 in dataset 1.

doi:10.1371/journal.pone.0135799.g003

Aox1 in dataset 1 and Crsp2 in dataset 2 have a maximum log2 (Fold Change) expression over -1.5. However these genes were selected in Volcanoplots and displayed interesting expression patterns, with a peak at low concentrations in dataset 1 and/or at early time points in dataset 2. They were then kept as interesting genes.

In order to focus on relevant genes potentially related to early transdifferentiation processes, the final step was to compare the final list of 20 DEGs to differential gene expression existing between female and male control fish. For this purpose, we compared the list of 20 DEGs with all the DEGs in the contrasts in dataset 2 comparing DEGs between control females versus control males in the developmental stages considered (MC-FC D7, D16 and D30, Table of contrasts in Fig 1 in blue). Ten out of the 20 genes were also differentially expressed between control females and control males (Fig 5A-5T). Interestingly, the 10 genes also differentially expressed in females (Fig 5A-5J) are mostly underexpressed genes and very similar expression profiles between EE2-treated fish and control females are observed at all developmental stages considered. The 10 remaining genes (Fig 5K-5T), considered specific to EE2 treatment, were mostly overexpressed genes, with peak expression at D30. They shared very similar expression patterns.

# Discussion

## 4.1. Differentially expressed genes and pathways

Both dataset analyses generated a high number of differentially expressed genes (DEG) following EE2 treatment, with nearly half of the array content displaying differential expression patterns in dataset 1 (12,742 non-redundant genes), and one third in dataset 2 (2,974 non-redundant genes). Since EE2 treatment is able to induce complete sex reversal of fish gonads, it is not surprising that this impact was detected at the transcriptional level. A similar magnitude of gene expression disruption was also observed in female fish masculinized by androgen 11β-hydroxyandrostenedione [6]. This highlights the major role of steroid hormones in triggering the gonad transcriptome, and their potency. Interestingly, 1,535 DEGs were found to be common to both analyses. These genes are therefore of particular interest considering their implication in the response of juvenile rainbow trout to EE2 treatment. The bioinformatics workflow applied in this study was previously shown to be effective in retrieving pathways in order to focus on potential novel genes associated with cancer metastasis and hypoxia in humans [19,22,23], and to identify potential biomarkers of intersex in rainbow trout [24].

Intersections between groups of DEGs highlight genes implicated in the response of juvenile rainbow trout to exogenous EE2 treatment. To focus on relevant genes related to the male-to-female transdifferentiation process, we further filter on phenotypic anchoring, which consists of encompassing different levels of biological organization [11]. We choose the CT1 and CT2 contrasts as the most relevant in dataset 1 as they represent mostly fishes in the intersex stage [25]. In dataset 2, the stages covering early differentiation, before gametogenesis, (D0 to D30) were chosen as the most relevant in the context studied [29]. Since sexual dimorphic gene expression patterns have been shown to appear from 35 dpf in rainbow trout [34,35], the earliest impact of EE2 on rainbow trout differentiation processes has to be considered.

GO Terms and biological pathways helps interpreting gene expression data in terms of the affected biological processes, pathways and phenotypic changes and on the mode of action and
A **Expression profiles in dataset 1**

<figure><img src="image_11.png" /></figure>

B **Expression profiles in dataset 1**

<figure><img src="image_12.png" /></figure>
**Fig 4. A-B: Selected genes kinetics profiles over the several concentrations tested in dataset 1.** These graphs show the profile of the 20 DEGs retrieved by the analysis over the 4 concentrations tested in dataset 1. X-axis: EE2 concentrations (0 = control). Y-axis: fold change of the different EE2 concentrations vs CTL. **A:** expressions of the 10 genes also differentially expressed between female and male control fish. **B:** expressions of the 10 genes specific to EE2 treatment.

doi:10.1371/journal.pone.0135799.g004

toxicity of a compound [36] including exposure to EDCs [6,37,38]. However especially in fish, namely here the rainbow trout, ORA approaches lacks of resources, genome annotations being scarce and/or heterogeneous. We thus chose to work by gene homology with *Danio rerio*, the only fish for which accurate pathways maps were available on Kegg [39]. A unified annotation of the platforms were made by retrieving *Danio rerio* ensembl gene IDs (ENSDARG) [24]. Considering that many genes and cellular functions are well conserved across taxa [40], the results here can be extrapolated as the general mode of action of EE2 on the juvenile male fish testis transcriptome. In Table 2, 12 pathways and 23 GO Terms appeared as significantly enriched in our results, 3 being significant after the Benjamini correction, shown to be too conservative for the EASE score in previous studies [19,22–24].

## 4.2. Pathways and ontologies associated with EE2 exposure in juvenile fish

We have previously shown that homology with *Danio rerio* retrieved several pathways associated with testis-ova development in juvenile rainbow trout exposed to low doses of EE2 [24]. The meta-analysis also identified numerous pathways as significantly enriched with differentially expressed genes. The advantage of the meta-analysis that it looks at DEGs resulting from EE2 exposure in juvenile fish in several conditions (i.e. several low and high EE2 doses; short and long exposure time), thus enhancing their relevancy as a fingerprint of EE2 exposure.

The general pattern of pathways and GO Terms obtained with the meta-analysis was similar to that obtained in our previous study [24], with the Cell Cycle and DNA replication being the most significantly enriched pathways, and several others related to Genetic information processing and Metabolism (Nucleotide excision repair, Mismatch repair, RNA degradation, Pyrimidine metabolism, Valine, leucine and isoleucine degradation). This was probably due to the high number of genes common to both analyses. This confirms the pathways proposed by Moggs [36] and Fertuck et al. [41] for E2 and EE2 in their study based on the rodent uterotrophic assay. Two new interesting pathways involved in the Endocrine system, namely the “progesterone-mediated oocyte maturation” and “PPAR signaling” pathways, emerged from this analysis. To our knowledge, this is the first time they have been reported after exposure to EE2 in juvenile fish. Half of the DEGs enriched in the “progesterone-oocyte maturation” pathway are common to the “Cell cycle” and correspond to genes expressed in the germline and involved in cell division (several members of the cyclin family: *ccna1*, *ccnb1*, *ccnb2*, *cdk2*; *mad2l1*, *anacp5*, *cdc16*) [42]. *ccna1* is a key gene involved in the control of meiosis [42]. Several kinases shared with the oocyte meiosis pathways are also represented (*rps6ka1*, *mos*). They are implicated in the phosphorylation of different substrates in the mitogen-activated protein kinase (MAPK) signaling pathway involved in cell growth and differentiation [43]. Other ubiquitous genes coding for chaperones (*hsp90Ab1*, *hsp90a.2.*) or G proteins (*gnai2*, *gnai3*) are represented. They are involved in numerous cellular processes. In view of these results, it is difficult to assess whether the overrepresentation of this pathway is specific to the progesterone-mediated response or is due to the differential expression of ubiquitous genes implicated in cellular processes, and especially cell division. However, during normal male testis development, meiosis occurred later than in females. These early meiotic signatures in males exposed to EE2 might therefore originate from pre-vitellogenic oocyte development. It will be then interesting
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /></figure>
**Fig 5. A-T: Selected genes kinetics profiles over the several sampling time in dataset 2.** These graphs show the profile of the 20 DEGs retrieved by the meta-analysis in dataset 2. Moreover, a comparison between gene expression between male exposed to EE2 versus male control and the expression in a control female versus male control is made. Y-axis: fold change of gene expression between males exposed to 20mg/kg of food of ethynylestradiol (dotted lines) or control females (solid lines) at each developmental stage versus male control at the same developmental stage (contrasts MT-MC D7 to D111 and FC-MC D0 to D111, Fig 1). X-axis: sampling time. Day 0 = D0 at 55 dpf, Day 7 = D0+7 days, D16, D30, D63, D91, and D111.

doi:10.1371/journal.pone.0135799.g005

to study this pathway in greater detail with reference to rainbow trout physiology, once further information becomes available.

The overrepresentation of the peroxisome proliferator-activated receptor (PPAR) signaling pathway is very interesting in the studied context. Peroxisome proliferator-activated receptors (PPARs) are ligand-inducible transcription factors belonging to the nuclear hormone (steroid/ thyroid/retinoid) receptor (NHR) superfamily [44]. PPARs regulate gene expression by binding with RXR (Retinoid X Receptor) as a heterodimeric partner to specific DNA sequence elements termed PPREs (Peroxisome Proliferator Response Elements). Target genes are mostly implicated in lipid transport and metabolism [45] and, in recent years, PPARs have emerged as critical regulators of lipid homeostasis in mammals [46]. Moreover, PPARs could play a major role in testicular somatic cells [45]. To date, three PPAR subtypes have been described in mammals and other species: PPARα, PPARβ and PPARγ [47]. They display different tissue distribution and function [46]. In *Xenopus laevis*, oocytes express high levels of PPARβ and smaller amounts of PPARα [44]. Endogenous ligands are mostly fatty acids and their derivatives, but a wide range of various natural and synthetic ligands of PPARs have been described [46]. Several estrogenic compounds (17b-estradiol, NP, genistein) have been shown to interact with PPARs [47–50]. PPARs are also upregulated in testis of rainbow trout treated with T or 11-KT during spermatogenesis [45]. Moreover, both laboratory and field studies have shown that a range of environmental pollutants (including PCBs, alkylphenols and estrogens) provoke peroxisome proliferation in different fish or bivalve mollusc species via the activation of PPARs [51]. Interestingly, the PPAR/RXR heterodimer-responsive element has been characterized in the promoter regions of the zebrafish and *Fundulus heteroclitus* *cyp19a1* genes [49,52]. This indicates that ligands of PPAR may affect the regulation of aromatase genes and therefore have an impact on the reproductive and developmental processes regulated by E2. It has been shown that several PPAR ligands have an inhibitory effect on aromatase expression in normal human breast adipose tissue [53] and ovarian granulosa cancer cells [54]. All these results argue in favor of a possible alteration of the PPAR signaling pathway by ethynylestradiol, as is suggested by our meta-analysis. Further studies will be necessary to determine whether the alteration observed is directly mediated by EE2, and related to a disruption in lipid metabolism regulation and homeostasis, or is involved in early sex differentiation.

Overall, the GO term enrichment was highly correlated with the pathway analysis, with a majority of terms related to genes involved in cell division. Terms related to sexual reproduction (with an enrichment of genes related to male gamete maturation and spermatogenesis) were more represented and appeared in the top list in terms of significance. The term "Oxido-reductase activity", which covers genes coding for key enzymes involved in steroidogenesis (*cyp17a1*, *cyp19a1*, *hsd11b2*), was absent in the pathway analysis.

### 4.3. Relevant genes potentially implicated in the early male-to-female transdifferentiation process

The final step of our analysis retrieved 20 statistically and physiologically relevant genes potentially implicated in the male-to-female transdifferentiation process in trout. Indeed, those genes are expressed shortly after EE2 treatment, at early developmental stages (7–16 and 30
<table><caption>**Table 4.** Description of the 20 genes selected as physiologically relevant in the gonad transdifferentiation process. Fold change (FC) and pvalues (Pval) are given for the five contrasts considered (CT1 and CT2 in dataset 1; MT-MC D7, D16 and D30 in dataset 2), upon which the genes were selected. Figures in bold represents significant values. The membership of genes to significantly enriched pathway and/or GO Terms is indicated.</caption><thead><tr><th>Gene symbol</th><th>Description</th><th>FC CT1</th><th>Pval CT1</th><th>FC CT2</th><th>Pval CT2</th><th>FC MTMC D7</th><th>Pval MTMC D7</th><th>FC MTMC D16</th><th>Pval MTMC D16</th><th>FC MTMC D30</th><th>Pval MTMC D30</th><th>Pathway</th><th>GO Term</th></tr></thead><tbody><tr><td>CYP11B</td><td>Cytochrome P450 11B</td><td>-3.8</td><td>9.8E-04</td><td>-200.3</td><td>3.3E-15</td><td>-1.93</td><td>1.2E-03</td><td>-2.0</td><td>3.7E-04</td><td>-6.2</td><td>6.0E-11</td><td>-</td><td>-</td></tr><tr><td>Cyp17a1</td><td>Steroid 17-alpha-hydroxylase/17,20 lyase</td><td>-5.2</td><td>2.7E-03</td><td>-10.8</td><td>3.7E-06</td><td>-1.34</td><td>5.5E-01</td><td>-1.5</td><td>9.1E-02</td><td>-3.0</td><td>5.7E-05</td><td>-</td><td>Oxidoreductase activity</td></tr><tr><td>Csrp2</td><td>Cysteine and glycine-rich protein 2</td><td>-1.1</td><td>1.5E-02</td><td>-6.9</td><td>6.7E-04</td><td>-1.39</td><td>5.6E-01</td><td>-2.7</td><td>2.1E-03</td><td>-1.8</td><td>9.7E-02</td><td>-</td><td>Cell proliferation</td></tr><tr><td>Lamb1</td><td>Laminin subunit beta-1</td><td>-1.2</td><td>6.5E-01</td><td>-5.8</td><td>8.9E-09</td><td>-1.89</td><td>1.4E-02</td><td>-2.9</td><td>2.4E-05</td><td>-4.4</td><td>6.3E-08</td><td>-</td><td>-</td></tr><tr><td>Tbx1</td><td>T-box transcription factor 1</td><td>1.1</td><td>9.1E-01</td><td>-5.2</td><td>3.7E-07</td><td>-1.41</td><td>7.1E-01</td><td>-2.2</td><td>7.1E-02</td><td>-3.8</td><td>1.1E-03</td><td>-</td><td>Cytoplasm</td></tr><tr><td>Aox1</td><td>Aldehyde oxidase</td><td>-1.4</td><td>2.8E-01</td><td>-1.6</td><td>4.8E-02</td><td>-3.03</td><td>3.6E-03</td><td>-7.6</td><td>7.2E-07</td><td>-3.0</td><td>1.1E-03</td><td>-</td><td>Oxidoreductase activity</td></tr><tr><td>FOXF1</td><td>Forkhead box protein F1</td><td>1.3</td><td>5.4E-01</td><td>2.4</td><td>2.0E-03</td><td>1.27</td><td>6.6E-01</td><td>2.2</td><td>3.8E-03</td><td>3.0</td><td>3.6E-05</td><td>-</td><td>-</td></tr><tr><td>PLTP</td><td>Phospholipid transfer protein</td><td>1.2</td><td>7.9E-01</td><td>6.0</td><td>2.8E-05</td><td>1.14</td><td>7.7E-01</td><td>1.5</td><td>4.6E-02</td><td>2.4</td><td>6.6E-06</td><td>PPAR signaling</td><td>-</td></tr><tr><td>Vtg1</td><td>Vitellogenin</td><td>3.2</td><td>2.6E-01</td><td>93.4</td><td>7.1E-07</td><td>2.32</td><td>8.3E-01</td><td>2.3</td><td>5.7E-01</td><td>27.2</td><td>3.4E-02</td><td>-</td><td>-</td></tr><tr><td>Lim2</td><td>Lens fiber membrane intrinsic protein</td><td>-5.3</td><td>7.6E-06</td><td>-211.8</td><td>2.5E-17</td><td>-1.20</td><td>8.0E-01</td><td>-1.7</td><td>7.1E-02</td><td>-3.4</td><td>7.1E-04</td><td>-</td><td>-</td></tr><tr><td>P4hb</td><td>Protein disulfide-isomerase</td><td>-1.5</td><td>5.1E-01</td><td>-2.9</td><td>8.8E-03</td><td>-1.88</td><td>5.6E-01</td><td>-5.7</td><td>3.8E-03</td><td>-2.0</td><td>5.3E-01</td><td>-</td><td>-</td></tr><tr><td>ARL6IP1</td><td>ADP-ribosylation factor-like protein 6-interacting protein 1</td><td>1.6</td><td>1.7E-02</td><td>1.7</td><td>2.2E-03</td><td>1.23</td><td>8.4E-01</td><td>1.2</td><td>5.7E-01</td><td>4.1</td><td>1.5E-04</td><td>-</td><td>-</td></tr><tr><td>gsta</td><td>Glutathione S-transferase A</td><td>1.4</td><td>2.6E-01</td><td>2.7</td><td>6.2E-05</td><td>1.12</td><td>9.5E-01</td><td>4.7</td><td>2.2E-02</td><td>1.0</td><td>1.0E+00</td><td>-</td><td>-</td></tr><tr><td>HSD17b6</td><td>17-beta-hydroxysteroid dehydrogenase type 6</td><td>1.4</td><td>2.6E-01</td><td>5.1</td><td>1.3E-08</td><td>1.20</td><td>9.2E-01</td><td>-1.4</td><td>5.0E-01</td><td>3.4</td><td>1.7E-02</td><td>-</td><td>-</td></tr><tr><td>alphaTub84D</td><td>Tubulin alpha chain</td><td>1.6</td><td>2.7E-02</td><td>5.1</td><td>1.7E-08</td><td>1.15</td><td>9.5E-01</td><td>-1.2</td><td>7.5E-01</td><td>107.4</td><td>4.4E-10</td><td>-</td><td>-</td></tr><tr><td>rsph1</td><td>Radial spoke head 1 homolog</td><td>1.5</td><td>2.1E-01</td><td>6.2</td><td>1.4E-08</td><td>-1.01</td><td>9.9E-01</td><td>1.2</td><td>1.8E-01</td><td>2.5</td><td>1.4E-06</td><td>-</td><td>-</td></tr><tr><td>kidins220</td><td>Kinase D-interacting substrate of 220 kDa</td><td>2.1</td><td>1.5E-03</td><td>8.9</td><td>1.9E-12</td><td>1.23</td><td>8.0E-01</td><td>-1.7</td><td>1.0E-01</td><td>-3.0</td><td>5.7E-03</td><td>-</td><td>-</td></tr><tr><td>rsph9</td><td>Radial spoke head protein 9</td><td>2.2</td><td>2.2E-02</td><td>10.2</td><td>1.7E-09</td><td>1.01</td><td>1.0E+00</td><td>-1.0</td><td>8.8E-01</td><td>2.1</td><td>1.8E-05</td><td>-</td><td>Cytoplasm</td></tr><tr><td>Tekt1</td><td>Tektin-1</td><td>2.0</td><td>2.4E-02</td><td>12.1</td><td>5.8E-11</td><td>-1.64</td><td>3.5E-01</td><td>-2.2</td><td>2.5E-02</td><td>-2.9</td><td>1.0E-03</td><td>-</td><td>-</td></tr><tr><td>RSPH3</td><td>Radial spoke head protein 3</td><td>4.7</td><td>3.1E-04</td><td>35.4</td><td>1.6E-11</td><td>-1.09</td><td>9.1E-01</td><td>1.5</td><td>8.5E-02</td><td>3.6</td><td>1.5E-06</td><td>-</td><td>Cytoplasm</td></tr></tbody></table>

doi:10.1371/journal.pone.0135799.t004
days after the treatment, with the fish being at 55 dpf on day 0), and in intersex fish following chronic exposure to low doses of EE2 (namely 0.01 and 0.1 µg/L). Among those genes, ten (kidins220, tekt1, gsta, hsd17b6, aox1, lamb1, cyp17a1, tbx1, cyp11b and lim2) were found to be also differentially expressed during normal female development compared to control males, at early developmental stages (D0 to D30 in dataset 2), with a very similar pattern (Fig 5A-5J). Most of these genes were underexpressed, except gsta and hsd17b6. Only the expression patterns of kidins220 and tekt1 were non-similar, as they were overexpressed in dataset 1 and underexpressed in dataset 2. The other ten genes displayed expression patterns specific to the EE2 treatment (vtg, rsph1-3-9, pltp, foxf1, αtub84d, p4hb, csrp2, arlip6) (Fig 5K-5T). They all showed overexpression compared to control males (except p4hb and csrp2), with a peak at D30 in dataset 2. Female expression was similar to that of control males (except for p4hb). This indicates that alterations in the expression of genes naturally involved in the early female differentiation process may be involved in the early EE2-induced feminization.

**4.3.1. Genes involved in natural early female differentiation.** Three of these genes encode key steroidogenic enzymes involved in 11-oxygenated androgen (including 11-KT) biosynthesis in the testis [55]. The 17-α-hydroxylase/17,20 lyase (encoded by *cyp17a1*) converts progesterone into 17α-hydroxyprogesterone and androstenedione through its hydroxylase and lyase activity, respectively. 17β-hydroxylase (encoded by *hsd17b*) then converts androstenedione into testosterone, which in turn is transformed into 11β-hydroxytestosterone, a precursor of 11-KT, by 11β-hydroxylase (encoded by *cyp11b*) [56]. 11-KT is a key androgen involved in testicular differentiation and spermatogenesis in most fishes [57]. In rainbow trout, *cyp11b* and *cyp17a1* are marker genes of early testicular development, and their expression is upregulated in males gonads with a clear sex dimorphism starting at 35 dpf for *cyp11b*, and later (70 dpf) for *cyp17a1* [29,34]. This overexpression in the testis continues throughout male gametogenesis [58,59]. Conversely, *hsd17b* shows marked overexpression during the early oogenesis period (from d60 to d110) and has been shown to be overexpressed during oocyte development [29,57,58]. The alteration of steroidogenesis following exposure to exogenous hormonal treatment or environmental contaminants has been extensively reported in several fish species [8,38,60-63]. In our meta-analysis, *cyp17a* and *cyp11b* were found to be underexpressed under the treatment. Both *in vivo* and *in vitro* studies have shown rapid, strong underexpression of *cyp11b* and *cyp17a* mRNA following E2 or EE2 treatment during and after sex differentiation [8,60,61,63,64], which supports a direct effect of E2 on the testis. This leads to a decrease in circulating androgen levels, as has been measured in several experiments. *Cyp17a1* and *hsd17b* are underexpressed following a short period of exposure to E2 and EE2 in male fathead minnows, whereas their expression is overexpressed in females [60]. Inversely, *hsd17b* was overexpressed in our results, in keeping with its natural overexpression observed in females (Fig 5H). Several types of 17βHSD enzymes exist and their relative roles have not been well defined [65]. Considering that testosterone is the precursor of E2, this enzyme is also implicated in regulating estrogen availability. It is then unclear why this enzyme is overexpressed. However, decreased levels of 11KT were recorded in exposed fish from dataset 1 (article 1). Altogether, the presence of these three genes among our candidates confirms their important role as key regulators of the feminization process, probably through a decrease in androgen levels. More importantly, feminization seems to rely essentially, though not exclusively, on the quick and strong decrease in *cyp11b* expression, whose expression profiles and significance are consistent in all the contrasts considered in our meta-analysis (Fig 3). Moreover, *cyp11b* has already been identified as a key actor triggering the male-to-female transdifferentiation process [8,64].

Another early testicular marker of testis development (*tbx1*) was found among the relevant genes, and was strongly underexpressed in both analyses (dose-dependently in dataset 1 -Fig 4A; and with a peak at D30 in dataset 2- Fig 5D), which is consistent with previous studies
showing strong downregulation of tbx1 by EE2 [66]. tbx1 is a member of the T-box transcription factor gene family involved in embryogenesis and organogenesis, whose expression is restricted to the somatic cells in the testis and appears to be sexually dimorphic very early in differentiating gonads (from 35 dpf) [66,67]. tbx1 is also expressed in the adult testis, and its expression decreases sharply throughout the reproductive cycle [42,68]. Our results support its role in the control and regulation of testicular differentiation, and suggest the importance of its downregulation during the feminization process.

lamb1 and tekt1 encode proteins belonging to the laminin and tektin protein families which are structural macromolecules found in the basal lamina and cytoskeleton, respectively. The basement membrane related gene (lama5) has been reported to be involved in gonadal development of both sexes in rat and fish [29]. Genes encoding tektins are predominantly expressed in the testis within the germline and during late spermatogenesis (stage IIId and V) and are involved in microtubule cytoskeleton organization, ciliary or flagellar mobility [42]. In mice, this gene seems to play a role in spermatogenesis and tekt1 mRNA are measured in spermatocytes and spermatids within seminiferous tubules [69]. Several other genes sharing similar expression profiles with normal females, have never, to our knowledge, been reported as being directly involved in sex differentiation or feminization. Aox1 belongs to a highly conserved family of proteins with oxidoreductase activity. The physiological functions and substrates of AOXs are still unknown, although they are recognized as major drug metabolizing enzymes [70]. The expression profiles of the Lim2 gene are of particular interest, as it was strongly underexpressed in both datasets, with increased repression between D7 and D111. It was reported in the enriched “Organogenesis” GO Term in our previous study [24]. However, this gene encodes the intrinsic Lens fiber membrane protein, with unknown function. Kidins220 (220 kDa kinase D-interacting substrate) is a substrate of protein kinase D (PKD). This gene encode for a membrane protein mostly expressed in brain and neuroendocrine cells, with no detectable mRNA levels in the rat testis [71]. Finally, glutathione S-transferase A (encoded by gsta) belongs to a family of isoenzymes recognized as classical markers of toxicity, due to their well described role in the detoxification of endogenous and exogenous substances [72]. Numerous studies have reported its overexpression (mostly in the liver) following exposure to xenobiotics [73,74]. Further investigation will be necessary to determine why all of these gene expressions are altered during natural early female development, why they respond to EE2 treatment, and their potential role in early sex differentiation.

### 4.3.2. Genes specific to EE2 exposure.
Among those genes, the well-known biomarker of estrogen exposure VTG was selected, which confirms the validity of our analysis to retrieve important genes responding to EE2 treatment. The gene atub84d encoding a cytoskeleton component was found to be overexpressed during spermatogenesis [42] and the gene product of ARL6IP1 was reported to be involved in protein transport, membrane trafficking or cell signaling during hematopoietic maturation in vitro [75,76] and in neural crest development in zebrafish embryos in vivo [77]. Its in vivo roles during embryonic development are still largely unknown. Three genes encoding members of the radial spoke head protein family (rsph 1, rsph 2 and rsph 9) were selected in our meta-analysis. They were found to be upregulated by the EE2 treatment in both datasets, with a peak of expression between D16 and D91 in dataset 2. These proteins are involved in various functions such as flagella motility, sensory functions and development. rsph3 has been reported among genes preferentially expressed in meiotic and postmeiotic germ cells in the testis of rainbow trout and were upregulated following T supplementation in vivo [78]. p4hb displays a peculiar pattern of expression, with clear underexpression between D7 and D30 closely related to the female curve. However, numerous studies reported its overexpression in the liver of several fish species, including rainbow trout, following estrogen treatment [79–81], and even proposed it as an effective biomarker of xenoestrogen
exposure. A specific profile of expression related to early development in the testis has thus been identified here and must be confirmed. *pltp* encodes for a protein involved in lipid transport in the PPAR signaling pathway. Considering the possible impact of the alteration of this pathway on expression levels of aromatase, the key enzyme that triggers the female developmental processes (see pathway discussion section above), the implication of this pathway in sex (trans)differentiation would constitute a very interesting topic for further study. In view of these results, these gene products do not seem to be directly involved in differentiation processes, or are involved to a still unknown extent. Their expression profiles might be further investigated as biomarkers of estrogen exposure during early differentiation in rainbow trout.

Two other genes in this top list, namely *foxf1* and *csrp2*, are transcription factors whose altered expression can lead to numerous changes in downstream gene transcript expression. *foxf1* belongs to the family of forkhead transcription factors, known as key players in development and metabolism (*foxl2* is a important gene triggering female differentiation and development in fish) [82]. *foxf1* expression during early development in rainbow trout is located mainly in the digestive tract and gills, and is associated with the development of gill filaments and lamellae [83]. The presence of ERE sequences in the promoter region of several fox genes (including *foxf1*) suggests that their transcriptional regulation is under ER control [84]. This is supported, at least for *foxf1*, by our results showing overexpression of this gene in both datasets. Cysteine and glycine-rich protein 2 (*csrp2*) is involved in cell differentiation. Its expression is naturally dimorphic during early differentiation, and it is underexpressed in females. Following EE2 treatment, it appears as strongly repressed in dataset 1 and to a lesser extent in dataset 2. This repression of *csrp2* expression was also reported after exposure to sewage in Chinook salmon (*Oncorhynchus tshawytscha*) [85].

**Conclusions and Perspectives**

In conclusion, our bioinformatics workflow meta-analysis successfully retrieved pathways related to estrogen exposure in fish and suggested novel ones. This led to the formulation of new hypotheses to be tested in order to elucidate the mechanisms underlying the modes of actions of estrogen, particularly in juvenile fish. Comparative homology analysis with *Danio rerio* thus appeared to be a relevant and powerful approach and suggests that its use could be extended to other fish species for environmental risk assessment purposes. Moreover, our approach identified several genes potentially implicated in the early trans-differentiation process. The expression of half of these genes was closely related to genes differentially expressed during natural female differentiation, including several markers repressed during testis development (*cyp11b*, *cyp17a1*, *tbx1*). Of particular interest is the quick and strong underexpression of *cyp11b*, which encodes the enzyme responsible for androgen production, which appears to be a key step in the feminization process. Other genes (including VTG) were found to respond specifically to the estrogen treatment. This underscores that EE2 acts by different mechanisms to trigger feminization induction in males, some closely related to normal mechanisms in developing females, other triggering the expression of genes not usually expressed in a natural feminization process. Although the expression of several of these genes has been verified on *in vitro* models (i.e. *cyp11b* in rainbow trout [63], *cyp17a1* in zebrafish [86], VTG in numerous species, including rainbow trout [87], suggesting a direct response to EE2), these results must still be validated to assess whether the altered expression of all of these genes is the result of a direct steroid action of the estrogen on the testis, or reflects an indirect retro control on the hypothalamus–pituitary axis, or a conjunction of both.
**Supporting Information**

S1 Dataset. Pre-processed data of dataset 2.

(XLSX)

S1 Results. DAVID software output from Intersection I DEG group analysis. Exhaustive list of pathways obtained with Intersection I DEG group from the meta-analysis between dataset 1 and dataset 2.

(DOCX)

S2 Results. EASE software output from Intersection II DEG group analysis. Exhaustive (p<0.01) list of ontologies obtained with Intersection II DEG group from the meta-analysis between dataset 1 and dataset 2. Genes belonging to each GO term are represented.

(XLSX)

**Acknowledgments**

The authors are grateful to J. Montfort and A. Le Cam from the INRA-LPGP/BIOSIT (Rennes, France) functional genomics platform for their help in the microarray experiment procedures. Thanks are also extended to the entire URBM bioinformatics team for their support in several of the statistical and bioinformatics steps of this study. Finally, we thank François Moreews from the Institut de Recherche en Informatique et Systèmes Aléatoires (IRISA) for his valuable help regarding the rainbow trout platform annotation. This work is supported by the FNRS (S. Depiereux), FRIA (M. Pierre) and FRS-FNRS Télévie (B. De Meulder).

**Author Contributions**

Conceived and designed the experiments: SD FLG YG PK ED. Performed the experiments: SD BDM MP RH. Analyzed the data: SD BDM RH. Contributed reagents/materials/analysis tools: SD FLG YG PK ED. Wrote the paper: SD.

**References**

1. Brennan J, Capel B. One tissue, two fates: molecular genetic events that underlie testis versus ovary development. Nat Rev Genet [Internet]. 2004 Jul [cited 2012 Apr 12]; 5(7):509–21. Available: http://www.ncbi.nlm.nih.gov/pubmed/15211353 PMID: 15211353

2. Devlin RH, Nagahama Y. Sex determination and sex differentiation in fish: an overview of genetic, physiological, and environmental influences. Aquaculture. 2002; 208:191–364.

3. Yamamoto T. Sex differentiation. In: Hoar W, Randall D, editors. Fish Physiology. 1969. p. 117–75.

4. Jobling S, Nolan M, Tyler CR, Brighty G, Sumpter JP. Widespread sexual disruption in wild fish. Environ Sci Technol. 1998; 32(17):2498–506.

5. Kidd KA, Blanchfield PJ, Mills KH, Palace VP, Evans RE, Lazorchak JM, et al. Collapse of a fish population after exposure to a synthetic estrogen. Proc Natl Acad Sci U S A [Internet]. 2007 May 22; 104 (21):8897–901. Available: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=1874224&tool=pmcentrez&rendertype=abstract PMID: 17517636

6. Baron D, Montfort J, Houlgatte R, Fostier A, Guiguen Y. Androgen-induced masculinization in rainbow trout results in a marked dysregulation of early gonadal gene expression profiles. BMC Genomics [Internet]. 2007 Jan [cited 2012 May 2]; 8:357. Available: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=2099445&tool=pmcentrez&rendertype=abstract PMID: 17916255

7. Vizziano D, Baron D, Randuineau G, Mahè S, Cauty C, Guiguen Y. Rainbow trout gonadal masculinization induced by inhibition of estrogen synthesis is more physiological than masculinization induced by androgen supplementation. Biol Reprod [Internet]. 2008 May [cited 2012 Mar 10]; 78(5):939–46. Available: http://www.ncbi.nlm.nih.gov/pubmed/18199883 doi: 10.1095/biolreprod.107.065961 PMID: 18199883

8. Vizziano-Cantonnet D, Baron D, Mahè S, Cauty C, Fostier A, Guiguen Y. Estrogen treatment up-regulates female genes but does not suppress all early testicular markers during rainbow trout male-to-
female gonadal transdifferentiation. J Mol Endocrinol [Internet]. 2008 Nov [cited 2012 Mar 10]; 41 (5):277-88. Available: http://www.ncbi.nlm.nih.gov/pubmed/18719050 doi: 10.1677/JME-08-0039 PMID: 18719050

9. Fodor S, Read J, Pirrung M, Stryler L, Lu A, Solas D. Light-directed, spatially addressable parallel chemical synthesis. Science (80-). 1991; 251(4995):767-73. PMID: 1990438

10. Nguyen C, Gidrol X. DNA Microarrays. Boisseau P, Houdy P, Lahmani M, editors. Nanoscience [Internet]. Berlin, Heidelberg: Springer Berlin Heidelberg; 2010 [cited 2012 Mar 29];911-35. Available: http://www.springerlink.com/index/10.1007/978-3-540-88633-4

11. Denslow ND, Garcia-Reyero N, Barber DS. Fish 'n' chips: the use of microarrays for aquatic toxicology. Mol Biosyst [Internet]. 2007 Mar [cited 2012 Mar 30]; 3(3):172-7. Available: http://www.pubmedcentral. nih.gov/articlerender.fcgi?artid=2259425&tool = pmcentrez&rendertype = abstract PMID: 17308663

12. Wang R-L, Biales A, Bencic D, Lattier D, Kostich M, Villeneuve D, et al. DNA microarray application in ecotoxicology: experimental design, microarray scanning, and factors affecting transcriptional profiles in a small fish species. Environ Toxicol Chem [Internet]. 2008 Mar; 27(3):652-63. Available: http:// www.ncbi.nlm.nih.gov/pubmed/17990945 PMID: 17990945

13. Berger F, De Meulder B, Gaigneaux A, Depiereux S, Bareke E, Pierre M, et al. Functional analysis: evaluation of response intensities—tailoring ANOVA for lists of expression subsets. BMC Bioinformatics [Internet]. 2010 Jan; 11:510. Available: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=2964684&tool = pmcentrez&rendertype = abstract doi: 10.1186/1471-2105-11-510 PMID: 20942918

14. Dudoit S, Yang YH, Callow MJ, Speed TP. Statistical methods for identifying differentially expressed genes in replicated cdna microarray experiments. Stat Sin. 2002; 12:111-39.

15. Opgen-rhein R, Strimmer K. Accurate Ranking of Differentially Expressed Genes by a Distribution-Free Shrinkage Approach. Stat Appl Genet Mol Biol. 2007; 6(1).

16. Smyth G. Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Genet Mol Biol. 2004; 3(1): Article 3.

17. Baldi P, Long AD. A Bayesian framework for the analysis of microarray expression data: regularized t-test and statistical inferences of gene changes. Bioinformatics. 2001; 17(6):509-19. PMID: 11395427

18. Berger F, De Hertogh B, Pierre M, Gaigneaux A, Depiereux E. The 'Window t test': a simple and powerful approach to detect differentially expressed genes in microarray datasets. Cent Eur J Biol. 2008; 3:327-44.

19. Pierre M, DeHertogh B, Gaigneaux A, DeMeulder B, Berger F, Bareke E, et al. Meta-analysis of archived DNA microarrays identifies genes regulated by hypoxia and involved in a metastatic phenotype in cancer cells. BMC Cancer [Internet]. 2010 Jan; 10:176. Available: http://www.pubmedcentral. nih.gov/articlerender.fcgi?artid=2880990&tool = pmcentrez&rendertype = abstract doi: 10.1186/1471- 2407-10-176 PMID: 20433688

20. Gur-dedeoglu B, Konu O, Kir S, Ozturk AR, Bozkurt B, Ergul G, et al. A resempling-based meta-analysis for detection of de differential gene expression in breast cancer. BMC Cancer. 2008; 16(8):396.

21. Rolland AD, Lardenois A, Goupil A-S, Lareyre J-J, Chalmel F, Houlgatte R, et al. A transcriptome-wide profiling identifies androgen-responsive genes in rainbow trout testis: relevance to male gonad development and spermatogenesis onset. 2010;

22. Pierre M, DeHertogh B, DeMeulder B, Bareke E, Depiereux S. Enhanced Meta-analysis Highlights Genes Involved in Metastasis from Several Microarray Datasets. J Proteomics Bioinform [Internet]. 2011 [cited 2013 Apr 21]; 04(02):36-43. Available: http://www.omicsonline.com/0974-276X/JPB-04-036.digital/JPB-04-036.html

23. De Meulder B, Berger F, Bareke E, Depiereux S, Michiels C, Depiereux E. Meta-Analysis and Gene Set Analysis of Archived Microarrays Suggest Implication of the Spliceosome in Metastatic and Hypoxic Phenotypes. Culhane AC, editor. PLoS One [Internet]. 2014 Jan 31 [cited 2014 Feb 4]; 9(1): e86699. Available: http://dx.plos.org/10.1371/journal.pone.0086699 doi: 10.1371/journal.pone.0086699 PMID: 24497970

24. Depiereux S, De Meulder B, Bareke E, Berger F, Le Gac F, Depiereux E, et al. Adaptation of a Bioinformatics Microarray Analysis Workflow for a Toxicogenomic Study in Rainbow Trout. PLoS One [Internet]. 2015; 10(7):e0128598. Available: http://dx.plos.org/10.1371/journal.pone.0128598 doi: 10.1371/journal.pone.0128598 PMID: 26186543

25. Depiereux S, Liagre M, Danis L, De Meulder B, Segner H, Depiereux E, et al. Intersex occurrence in rainbow trout (Oncorhynchus mykiss) male fry chronically exposed to ethynylestradiol. PLoS One. 2014; 9(7):e98531. doi: 10.1371/journal.pone.0098531 PMID: 25033040

26. Guiguen Y, Montfort J, Le Cam A, Yano A, Nicol B, Hugot K, et al. Microarray-based analysis of rainbow trout steroid-modulated gonadal differentiation.
27. De Hertogh B, De Meulder B, Berger F, Pierre M, Bareke E, Gaigneaux A, et al. A benchmark for statistical microarray data analysis that preserves actual biological and technical variance. BMC Bioinformatics [Internet]. 2010 Jan; 11:17. Available: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=2831002&tool=pmcentrez&rendertype=abstract doi: 10.1186/1471-2105-11-17 PMID: 20064233

28. Benjamini Y, Hochberg Y. Controlling the False Discovery Rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B. 1995; 57(1):289-300.

29. Baron D, Houlgatte R, Fostier A, Guiguen Y. Large-scale temporal gene expression profiling during gonadal differentiation and early gametogenesis in rainbow trout. Biol Reprod [Internet]. 2005 Nov [cited 2012 Mar 10]; 73(5):959-66. Available: http://www.ncbi.nlm.nih.gov/pubmed/16014816 PMID: 16014816

30. Flicek P, Amode MR, Barrell D, Al E. Ensembl 2012. Nucleic Acids Res. 2012; 40.

31. Dennis G, Sherman BT, Hosack D a, Yang J, Gao W, Lane HC, et al. DAVID: Database for Annotation, Visualization, and Integrated Discovery. Genome Biol [Internet]. 2003 Jan; 4(5):P3. Available: http://www.ncbi.nlm.nih.gov/pubmed/12734009 PMID: 12734009

32. Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc [Internet]. Nature Publishing Group; 2009 Jan [cited 2012 Jul 13]; 4(1):44-57. Available: http://dx.doi.org/10.1038/nprot.2008.211 doi: 10.1038/nprot.2008.211 PMID: 19131956

33. Hosack D a, Dennis G, Sherman BT, Lane HC, Lempicki R a. Identifying biological themes within lists of genes with EASE. Genome Biol [Internet]. 2003 Jan; 4(10):R70. Available: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=328459&tool=pmcentrez&rendertype=abstract PMID: 14519205

34. Vizziano D, Randuineau G, Baron D, Cauty C, Guiguen Y. Characterization of early molecular sex differentiation in rainbow trout, Oncorhynchus mykiss. Dev Dyn [Internet]. 2007 Aug [cited 2012 May 2]; 236(8):2198-206. Available: http://www.ncbi.nlm.nih.gov/pubmed/17584856 PMID: 17584856

35. Cavileer T, Hunter S, Okutsu T, Yoshizaki G, Nagler JJ. Identification of novel genes associated with molecular sex differentiation in the embryonic gonads of rainbow trout (Oncorhynchus mykiss). Sex Dev [Internet]. 2009 Jan [cited 2012 May 2]; 3(4):214-24. Available: http://www.ncbi.nlm.nih.gov/pubmed/19752601 doi: 10.1159/000228722 PMID: 19752601

36. Moggs JG. Molecular responses to xenoestrogens: mechanistic insights from toxicogenomics. Toxicology [Internet]. 2005 Oct 1 [cited 2013 Sep 23]; 213(3):177-93. Available: http://www.ncbi.nlm.nih.gov/pubmed/15996808 PMID: 15996808

37. Koskinen H, Pehkonen P, Vehniäinen E, Krasnov A, Rexroad C, Afanasyev S, et al. Response of rainbow trout transcriptome to model chemical contaminants. Biochem Biophys Res Commun [Internet]. 2004 Jul 30 [cited 2013 May 22]; 320(3):745-53. Available: http://www.ncbi.nlm.nih.gov/pubmed/15240111 PMID: 15240111

38. Wang R-L, Bencic D, Villeneuve DL, Ankley GT, Lazorchak J, Edwards S. A transcriptomics-based biological framework for studying mechanisms of endocrine disruption in small fish species. Aquat Toxicol [Internet]. Elsevier B.V.; 2010 Jul 1 [cited 2013 Jun 19]; 98(3):230-44. Available: http://www.ncbi.nlm.nih.gov/pubmed/20359757 doi: 10.1016/j.aquatox.2010.02.021 PMID: 20359757

39. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res [Internet]. 2000 Jan 1; 28(1):27-30. Available: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=102409&tool=pmcentrez&rendertype=abstract PMID: 10592173

40. Miller KM, Maclean N. Teleost microarrays: development in a broad phylogenetic range reflecting diverse applications. J Fish Biol [Internet]. 2008 Jun [cited 2012 Mar 6]; 72(9):2039-50. Available: http://doi.wiley.com/10.1111/j.1095-8649.2008.01913.x

41. Fertuck KC, Eckel JE, Gennings C, Zacharewski TR. Identification of temporal patterns of gene expression in the uteri of immature, ovariectomized mice following exposure to ethynylestradiol. Physiol Genomics [Internet]. 2003 Oct 17 [cited 2013 Sep 23]; 15(2):127-41. Available: http://www.ncbi.nlm.nih.gov/pubmed/12915738 PMID: 12915738

42. Rolland AD, Lareyre J-J, Goupil A-S, Montfort J, Ricordel M-J, Esquerré D, et al. Expression profiling of rainbow trout testis development identifies evolutionary conserved genes involved in spermatogenesis. BMC Genomics [Internet]. 2009 Jan [cited 2012 Apr 22]; 10:546. Available: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=2786911&tool=pmcentrez&rendertype=abstract doi: 10.1186/1471-2164-10-546 PMID: 19925684

43. Lalancette C, Platts AE, Lu Y, Lu S, Krawetz S a. Computational identification of transcription frameworks of early committed spermatogenic cells. Mol Genet Genomics [Internet]. 2008 Oct [cited 2013 Aug 31]; 280(3):263-74. Available: http://www.ncbi.nlm.nih.gov/pubmed/18615256 doi: 10.1007/s00438-008-0361-2 PMID: 18615256
44. Lemberger T, Desvergne B, Wahli W. Peroxisome proliferator-activated receptors: a nuclear receptor signaling pathway in lipid physiology. Annu Rev Cell Dev Biol [Internet]. 1996 Jan; 12:335-63. Available: http://www.ncbi.nlm.nih.gov/pubmed/8970730 PMID: 8970730

45. Gac F Le, Lareyre JJ, Montfort J, Esquerre D, Houlgatte R. Transcriptional analysis of spermatogenesis regulation by sex steroids in trout. Cybium. 2008; 32(2):119-21.

46. Leaver MJ, Boukouvala E, Antonopoulou E, Diez A, Favre-Krey L, Ezaz MT, et al. Three peroxisome proliferator-activated receptor isotypes from each of two species of marine fish. Endocrinology [Internet]. 2005 Jul [cited 2013 Sep 12]; 146(7):3150-62. Available: http://www.ncbi.nlm.nih.gov/pubmed/15790725 PMID: 15790725

47. Ibabe a, Herrero A, Cajaraville MP. Modulation of peroxisome proliferator-activated receptors (PPARs) by PPAR(alpha)- and PPAR(gamma)-specific ligands and by 17beta-estradiol in isolated zebrafish hepatocytes. Toxicol In Vitro [Internet]. 2005 Sep [cited 2013 Sep 24]; 19(6):725-35. Available: http://www.ncbi.nlm.nih.gov/pubmed/15964169 PMID: 15964169

48. Ruggeri B, Ubaldi M, Lourdusamy A, Soverchia L, Ciccocioppo R, Hardiman G, et al. Variation of the genetic expression pattern after exposure to estradiol-17beta and 4-nonylphenol in male zebrafish (Danio rerio). Gen Comp Endocrinol [Internet]. 2008 Aug [cited 2012 Mar 18]; 158(1):138-44. Available: http://www.ncbi.nlm.nih.gov/pubmed/18602103 doi: 10.1016/j.ygcen.2008.05.012 PMID: 18602103

49. Dong W, Wang L, Thornton C, Scheffler BE, Willett KL. Benzo(a)pyrene decreases brain and ovarian aromatase mRNA expression in Fundulus heteroclitus. Aquat Toxicol [Internet]. 2008 Jul 30 [cited 2013 Feb 22]; 88(4):289-300. Available: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=2530897&tool=pmcentrez&rendertype=abstract doi: 10.1016/j.aquatox.2008.05.006 PMID: 18571745

50. Campbell B, Dickey JT, Swanson P. Endocrine changes during onset of puberty in male spring Chinook salmon, Oncorhynchus tshawytscha. Biol Reprod [Internet]. 2003 Dec [cited 2013 Sep 29]; 69 (6):2109-17. Available: http://www.ncbi.nlm.nih.gov/pubmed/12930714 PMID: 12930714

51. Cajaraville MP, Cancio I, Ibabe A, Orbea A. Peroxisome proliferation as a biomarker in environmental pollution assessment. Microsc Res Tech. 2003; 61(2):191-202. PMID: 12740826

52. Kazeto Y, Ijiri S, Place a R, Zohar Y, Trant JM. The 5'-flanking regions of CYP19A1 and CYP19A2 in zebrafish. Biochem Biophys Res Commun [Internet]. 2001 Nov 2 [cited 2013 Sep 24]; 288(3):503-8. Available: http://www.ncbi.nlm.nih.gov/pubmed/11676471 PMID: 11676471

53. Rubin GL, Duong JH, Clyne CD, Speed CJ, Murata Y, Gong C, et al. Ligands for the Peroxisomal Proliferator-Activated Receptor gamma and the Retinoid X Receptor InhibitAromatase cytochrome P450 (CYP19) expression Mediated by Promoter II in Human Breast Adipose. Endocrinology. 2002; 143 (8):2863-71. PMID: 12130549

54. Mu YM, Yanase T, Nishi Y, Waseda N, Oda T, Tanaka a, et al. Insulin sensitizer, troglitazone, directly inhibits aromatase activity in human ovarian granulosa cells. Biochem Biophys Res Commun [Internet]. 2000 May 19 [cited 2013 Sep 29]; 271(3):710-3. Available: http://www.ncbi.nlm.nih.gov/pubmed/10814527 PMID: 10814527

55. Kusakabe M, Nakamura I, Evans J, Swanson P, Young G. Changes in mRNAs encoding steroidogenic acute regulatory protein, steroidogenic enzymes and receptors for gonadotropins during spermatogenesis in rainbow trout testes. J Endocrinol [Internet]. 2006 Jun [cited 2013 Sep 17]; 189(3):541-54. Available: http://www.ncbi.nlm.nih.gov/pubmed/16731785 PMID: 16731785

56. Villeneuve DL, Blake LS, Brodin JD, Greene KJ, Knoebel I, Miracle AL, et al. Transcription of key genes regulating gonadal steroidogenesis in control and ketoconazole- or vinclozolin-exposed fathead minnows. Toxicol Sci [Internet]. 2007 Aug [cited 2012 Mar 10]; 98(2):395-407. Available: http://www.ncbi.nlm.nih.gov/pubmed/17517826 PMID: 17517826

57. Pandian TJ. Endocrine sex differentiation in fish. CRC Press. Pandian TJ, editor. Boca raton London new York; 2013.

58. Baron D, Houlgatte R, Fostier A, Guiguen Y. Expression profiling of candidate genes during ovary-to-testis trans-differentiation in rainbow trout masculinized by androgens. Gen Comp Endocrinol [Internet]. 2008 Apr 1 [cited 2012 Mar 10]; 156(2):369-78. Available: http://www.ncbi.nlm.nih.gov/pubmed/18299129 doi: 10.1016/j.ygcen.2008.01.016 PMID: 18299129

59. Liu S, Govoroun M, D'Cotta H, Ricordel MJ, Lareyre JJ, McMeel OM, et al. Expression of cytochrome P450(11beta) (11beta-hydroxylase) gene during gonadal sex differentiation and spermatogenesis in rainbow trout, Oncorhynchus mykiss. J Steroid Biochem Mol Biol [Internet]. 2000 Dec 31; 75(4-5):291-8. Available: http://www.ncbi.nlm.nih.gov/pubmed/11282285 PMID: 11282285

60. Filby AL, Thorpe KL, Maack G, Tyler CR. Gene expression profiles revealing the mechanisms of anti-androgen- and estrogen-induced feminization in fish. Aquat Toxicol [Internet]. 2007 Feb 28 [cited 2012 Mar 10]; 81(2):219-31. Available: http://www.ncbi.nlm.nih.gov/pubmed/17222921 PMID: 17222921
61. Filby AL, Santos EM, Thorpe KL, Maack G, Tyler CR. Gene expression profiling for understanding chemical causation of biological effects for complex mixtures: a case study on estrogens. Environ Sci Technol [Internet]. 2007 Dec 1; 41(23):8187-94. Available: http://www.ncbi.nlm.nih.gov/pubmed/18186357 PMID: 18186357

62. Hinfry N, Baudiffier D, Leal MC, Porcher J-M, Aït-Aïssa S, Le Gac F, et al. Characterization of testicular expression of P450 17α-hydroxylase, 17,20-lyase in zebrafish and its perturbation by the pharmaceutical fungicide clotrimazole. *Gen Comp Endocrinol* [Internet]. Elsevier Inc.; 2011 Dec 1 [cited 2013 Jan 25]; 174(3):309-17. Available: http://www.ncbi.nlm.nih.gov/pubmed/21963961 doi: 10.1016/j.ygcen.2011.09.008 PMID: 21963961

63. Baron D, Fostier A, Breton B, Guiguen Y. Androgen and estrogen treatments alter steady state messengers RNA (mRNA) levels of testicular steroidogenic enzymes in the rainbow trout, Oncorhynchus mykiss. Mol Reprod Dev [Internet]. 2005 Aug [cited 2012 Oct 17]; 71(4):471-9. Available: http://www.ncbi.nlm.nih.gov/pubmed/15858796 PMID: 15858796

64. Govoroun M, McMeel OM, Mecherouki H, Smith TJ, Guiguen Y. 17-Estradiol Treatment Decreases Steroidogenic Enzyme Messenger Ribonucleic Acid Levels in the Rainbow Trout Testis. Endocrinology. 2001; 142(5):1841-8. PMID: 11316749

65. Adamski J, Jakob FJ. A guide to 17beta-hydroxysteroid dehydrogenases. Mol Cell Endocrinol [Internet]. 2001 Jan 22; 171(1-2):1-4. Available: http://www.ncbi.nlm.nih.gov/pubmed/11165003

66. Yano A, Nicol B, Guerin A, Guiguen Y. The duplicated rainbow trout (Oncorhynchus mykiss) T-box transcription factors 1, tbx1a and tbx1b, are up-regulated during testicular development. Mol Reprod Dev [Internet]. 2011 Mar [cited 2012 Mar 1]; 78(3):172-80. Available: http://www.ncbi.nlm.nih.gov/pubmed/21308851 doi: 10.1002/mrd.21279 PMID: 21308851

67. Yano A, Guyomard R, Nicol B, Jouanno E, Quillet E, Klopp C, et al. An immune-related gene evolved into the master sex-determining gene in rainbow trout, oncorhynchus mykiss. Curr Biol [Internet]. 2012 Jun 20 [cited 2012 Jul 23];null(null). Available: http://dx.doi.org/10.1016/j.cub.2012.05.045

68. Schulz RW, de França LR, Lareyre J-J, Le Gac F, LeGac F, Chiarini-Garcia H, et al. Spermatogenesis in fish. *Gen Comp Endocrinol* [Internet]. 2010 Feb 1 [cited 2012 Mar 10]; 165(3):390-411. Available: http://www.ncbi.nlm.nih.gov/pubmed/19348807 doi: 10.1016/j.ygcen.2009.02.013 PMID: 19348807

69. Larsson M, Norrander J, Gräslund S, Brundell E, Linck R, Sta S. The spatial and temporal expression of Tekt1, a mouse tektin C homologue, during spermatogenesis suggest that it is involved in the development of the sperm tail basal body and axoneme. Eur J Endocrinol. 2000; 79:718-25.

70. Kurosaki M, Bolis M, Fratelli M, Barzago MM, Pattini L, Perretta G, et al. Structure and evolution of vertebrate aldehyde oxidases: from gene duplication to gene suppression. Cell Mol Life Sci [Internet]. 2013 May [cited 2013 Sep 10]; 70(10):1807-30. Available: http://www.ncbi.nlm.nih.gov/pubmed/23263164 doi: 10.1007/s00018-012-1229-5 PMID: 23263164

71. Iglesias T, Cabrera-Poch N, Mitchell MP, Naven TJ, Rozengurt E, Schiavo G. Identification and cloning of Kidins220, a novel neuronal substrate of protein kinase D. J Biol Chem [Internet]. 2000 Dec 22 [cited 2013 Sep 14]; 275(51):40048-56. Available: http://www.ncbi.nlm.nih.gov/pubmed/10998417 PMID: 10998417

72. Salinas A, Wong M. Glutathione S-transferases—a review. Curr Med Chem. 1999; 6(4):279-309. PMID: 10101214

73. Mayon N, Bertrand A, Leroy D, Malbrouck C, Mandiki SNM, Silvestre F, et al. Multiscale approach of fish responses to different types of environmental contaminations: a case study. Sci Total Environ [Internet]. 2006 Aug 31 [cited 2012 Apr 4]; 367(2-3):715-31. Available: http://www.ncbi.nlm.nih.gov/pubmed/16740295 PMID: 16740295

74. Moens LN, Smolders R, van der Ven K, van Remortel P, Del-Favero J, De Coen WM. Effluent impact assessment using microarray-based analysis in common carp: a systems toxicology approach. Chemosphere [Internet]. 2007 May [cited 2012 Apr 20]; 67(11):2293-304. Available: http://www.ncbi.nlm.nih.gov/pubmed/17267021 PMID: 17267021

75. Pettersson M, Bessonova M, Gu HF, Groop LC, Jönsson JI. Characterization, chromosomal localization, and expression during hematopoietic differentiation of the gene encoding Arl6ip, ADP-ribosylation-like factor-6 interacting protein (ARL6). *Genomics* [Internet]. 2000 Sep 15 [cited 2013 Sep 12]; 68(3):351-4. Available: http://www.ncbi.nlm.nih.gov/pubmed/10995579 PMID: 10995579

76. Ingley E, Williams JH, Walker CE, Tsai S, Colley S, Sayer MS, et al. A novel ADP-ribosylation like factor (ARL-6), interacts with the protein-conducting channel SEC61beta subunit. FEBS Lett [Internet]. 1999 Oct 1; 459(1):69-74. Available: http://www.ncbi.nlm.nih.gov/pubmed/10508919 PMID: 10508919

77. Yang TC, Tu CT, Tsai HJ. ADP-ribosylation-like factor-6 interacting protein is required for neural crest development in zebrafish embryos. International Journal of Developmental Neuroscience. 2008. p. 867-92.
78. Rolland AD, Lardenois A, Goupil A-S, Lareyre J-J, Houlgatte R, Chalmel F, et al. Profiling of androgen response in rainbow trout pubertal testis: relevance to male gonad development and spermatogenesis. PLoS One [Internet]. 2013 Jan [cited 2013 Sep 9]; 8(1):e53302. Available: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=3536770&tool=pmcentrez&rendertype=abstract doi: 10.1371/journal.pone.0053302 PMID: 23301058

79. Tilton SC, Givan S a, Pereira CB, Bailey GS, Williams DE. Toxicogenomic profiling of the hepatic tumor promoters indole-3-carbinol, 17beta-estradiol and beta-naphthoflavone in rainbow trout. Toxicol Sci [Internet]. 2006 Mar [cited 2012 Apr 20]; 90(1):61–72. Available: http://www.ncbi.nlm.nih.gov/pubmed/16192472 PMID: 16192472

80. Katsiadaki I, Williams TD, Ball JS, Bean TP, Sanders MB, Wu H, et al. Hepatic transcriptomic and meta-bolomic responses in the Stickleback (Gasterosteus aculeatus) exposed to ethinyl-estradiol. Aquat Toxicol [Internet]. Elsevier B.V.; 2010 May 5 [cited 2013 Feb 21]; 97(3):174–87. Available: http://www.ncbi.nlm.nih.gov/pubmed/19665239 doi: 10.1016/j.aquatox.2009.07.005 PMID: 19665239

81. Larkin P, Knoebel I, Denslow ND. Differential gene expression analysis in fish exposed to endocrine disrup-
ting compounds. Comp Biochem Physiol Part B Biochem Mol Biol [Internet]. 2003 Oct [cited 2012
Mar 10]; 136(2):149–61. Available: http://linkinghub.elsevier.com/retrieve/pii/S1096495903002288

82. Carlsson P, Mahlapuu M. Forkhead Transcription Factors: Key Players in Development and Metabo-lism. Dev Biol [Internet]. 2002 Oct [cited 2013 Aug 27]; 250(1):1–23. Available: http://linkinghub.elsevier.com/retrieve/pii/S0012160602907803 PMID: 12297093

83. Hidaka Y, Tanaka S, Suzuki M. Sequence and expression of the rainbow trout winged helix/forkhead transcription factor FoxF1. Gene Expr Patterns [Internet]. 2004 Nov [cited 2013 Sep 10]; 5(1):29–34. Available: http://www.ncbi.nlm.nih.gov/pubmed/15533815 PMID: 15533815

84. Lin C-Y, Vega VB, Thomsen JS, Zhang T, Kong SL, Xie M, et al. Whole-genome cartography of estro-
gen receptor alpha binding sites. PLoS Genet [Internet]. 2007 Jun [cited 2012 Apr 28]; 3(6):e87. Avail-
able: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=1885282&tool =
pmcentrez&rendertype = abstract PMID: 17542648

85. Osachoff HL, van Aggelen GC, Mommsen TP, Kennedy CJ. Concentration-response relationships and temporal patterns in hepatic gene expression of Chinook salmon (Oncorhynchus tshawytscha) exposed to sewage. Comp Biochem Physiol Part D Genomics Proteomics [Internet]. Elsevier Inc.; 2013 Mar [cited 2013 Sep 9]; 8(1):32–44. Available: http://www.ncbi.nlm.nih.gov/pubmed/23246653 doi: 10.1016/j.cbd.2012.10.002 PMID: 23246653

86. De Waal PP, Leal MC, García-López A, Liarte S, de Jonge H, Hinfray N, et al. Oestrogen-induced androgen insufficiency results in a reduction of proliferation and differentiation of spermatogonia in the zebrafish testis. J Endocrinol [Internet]. 2009 Aug [cited 2013 Jan 25]; 202(2):287–97. Available: http://www.ncbi.nlm.nih.gov/pubmed/19420009 doi: 10.1677/JOE-09-0050 PMID: 19420009

87. Van den Belt K, Verheyen R, Witters H. Comparison of vitellogenin responses in zebrafish and rainbow trout following exposure to environmental estrogens. Ecotoxicol Environ Saf [Internet]. 2003 Oct [cited 2013 Feb 19]; 56(2):271–81. Available: http://linkinghub.elsevier.com/retrieve/pii/S0147651303000046 PMID: 12927559