Histological and transcriptomic effects of 17α-methyltestosterone on zebrafish gonad development

Sex hormones play important roles in teleost ovarian and testicular development. In zebrafish, ovarian differentiation appears to be dictated by an oocyte-derived signal via Cyp19a1a aromatase-mediated estrogen production. Androgens and aromatase inhibitors can induce female-to-male sex reversal, however, the mechanisms underlying gonadal masculinisation are poorly understood. We used histological analyses together with RNA sequencing to characterise zebrafish gonadal transcriptomes and investigate the effects of 17α-methyltestosterone on gonadal differentiation. At a morphological level, 17α-methyltestosterone (MT) masculinised gonads and accelerated spermatogenesis, and these changes were paralleled in masculinisation and de-feminisation of gonadal transcriptomes. MT treatment upregulated expression of genes involved in male sex determination and differentiation (amh, dmrt1, gsdf and wt1a) and those involved in 11-oxygenated androgen production (cyp11c1 and hsd11b2). It also repressed expression of ovarian development and folliculogenesis genes (bmp15, gdf9, figla, zp2.1 and zp3b). Furthermore, MT treatment altered epigenetic modification of histones in zebrafish gonads. Contrary to expectations, higher levels of cyp19a1a or foxl2 expression in control ovaries compared to MT-treated testes and control testes were not statistically significant during early gonad development (40 dpf). Our study suggests that both androgen production and aromatase inhibition are important for androgen-induced gonadal masculinisation and natural testicular differentiation in zebrafish.

Little is known about the direct transcriptional effects of MT during gonad masculinisation. In this study, we treated juvenile Tg(vas:egfp) transgenic zebrafish with 100 ng/L MT from 20 days post fertilisation (dpf ) to 40 dpf and 60 dpf. These time points were selected to reflect the onset (20 dpf ) and completion (40 dpf ) of juvenile ovary-to-testis transformation as well as the onset of female puberty (60 dpf ) in zebrafish [68,69,75]. RNA sequencing (RNA-Seq) was used to profile global gene expression patterns in MT-treated gonads. The gonadal transcriptomes of MT-treated zebrafish were compared with those of control zebrafish testes and ovaries to provide insights into the molecular basis for MT-induced gonadal masculinisation in zebrafish.

Ethics statement
This study was approved by the University of Otago Animal Ethics Committee (AEC No. 101/09). All experiments were performed in accordance with the Good Practice Guide for the use of animals in research, testing and teaching.

Zebrafish husbandry
Zebrafish were maintained according to Westerfield [76]. We used larval and juvenile transgenic zebrafish expressing an enhanced green fluorescent protein (EGFP) under the control of the vasa promoter, Tg(vas:egfp), derived from the domesticated Tübingen/AB strains [73]. EGFP expression in Tg(vas:egfp) transgenic zebrafish enables visualisation and isolation of gonads before the gonads can be unequivocally distinguished from other tissues [77]. Higher levels of fluorescence are detected in ovaries than testes [69,73], which can be used to distinguish the phenotypic sex of the fish.
17α-methyltestosterone treatment 17α-methyltestosterone (MT) was purchased from Sigma Aldrich (Sigma-Aldrich Sweden AB, Stockholm, Sweden) and dissolved in 100% ethanol to prepare stock solutions of 50 mg/L. At 18 to 19 dpf, juvenile zebrafish were transferred from 4 L tanks into petri dishes where they were size-sorted by visual inspection. Zebrafish that were unusually small (<5.5 mm) or large (>8.5 mm) were removed. The remaining juvenile zebrafish (between 5.5 to 8.5 mm in length) were transferred into 4 L tanks containing 3.5 L of system water at a density of 20 individuals per tank.
For the MT exposure, the juvenile zebrafish were exposed to system water containing either 100 ng/L of MT dissolved in 0.0001% ethanol (100 ng/L MT) or 0.0001% ethanol alone (solvent control). Each exposure group consisted of 3 biological replicates comprising 20 juvenile zebrafish each. The exposure was performed continuously under a semi-static system for 20 or 40 days extending from 20 dpf to 40 dpf or 60 dpf. Water with equivalent concentrations of MT and ethanol was used to replace half of the water in each beaker every second day for the MT treatment and solvent control groups respectively.

Determination of gonadal morphology and sex ratios
At the termination of the treatment at 40 dpf and 60 dpf, the juvenile zebrafish were sacrificed via snap chilling in ice water. EGFP expression of each zebrafish was observed under a Leica M205 FA fluorescence dissecting stereo microscope (Leica Microsystems, Bannockburn, Illinois, USA) to determine the gonadal sex of each fish and sex ratios for each exposure group. Between 18 to 24 fish per exposure group (comprising 6 to 8 fish from each of the 3 replicates per time point) were selected for histological analysis to confirm gonadal sex as determined by EGFP expression and to determine the developmental stages of each gonad. Whole zebrafish were fixed overnight in 10% neutral buffered formalin, dehydrated and embedded in paraffin. Haematoxylin and eosin-stained sections were examined under light microscopy to determine gonadal sex. Ovaries were classified as Stage I or Stage II [75], based on the absence or presence of Stage II oocytes containing cortical alveoli, respectively [78]. Testes were categorised into transition gonads, immature testes or mature testes. Gonads with basophilic apoptotic bodies and infiltration of stromal tissue were classified as transition gonads [68]. Testes were classified as immature or mature depending on the absence or presence of spermatozoa [79]. Spermatids differentiate into spermatozoa, through a series of morphological changes inclusive of nuclear compaction and flagellum formation, which are released into the seminiferous tubule lumen at the end of spermiogenesis [80].

RNA isolation
The juvenile zebrafish gonads were dissected and placed into aliquots of RNAlater (Invitrogen, Life Technologies GmbH / ThermoFisher Scientific, Darmstadt, Germany) or Lysis buffer (RLT buffer from an RNeasy Mini kit, Qiagen GmbH, Hilden, Germany or RA1 lysis buffer from a Nucleospin RNA II kit, Macherey-Nagel GmbH & Co. KG, Düren, Germany) and stored at −80°C for subsequent RNA isolation. The gonads from ten individuals per exposure group were pooled together for RNA extraction. Gonads were homogenized by passing through a 20 gauge needle or using 2 cycles of the TissueLyser II system (Qiagen GmbH, Hilden, Germany) set at 20 Hz for 1 min per cycle.
RNA was isolated using an RNeasy Mini kit (Qiagen GmbH, Hilden, Germany) or a Nucleospin RNA II kit (Macherey-Nagel GmbH & Co. KG, Düren, Germany) and eluted in 53ul of RNase-free water.
Three biological replicates of gonad pools consisting of 10 individuals each were generated for each exposure group (MT exposure and solvent control), sex (male or female) and time point (40 dpf or 60 dpf) (Fig. 1). The RNA samples were quantified and integrity was assessed using the Agilent 2100 Bioanalyser (Agilent Technologies, Palo Alto, California, USA). RNA samples with RNA Integrity Numbers (RIN) ≥7.0 (majority with RIN ≥ 8.0 to 9.0) were used for RNA sequencing (RNA-Seq).

Read annotation, mapping, assembly and quantification
The quality of the raw reads was assessed using the FastQC software [81] and the sequencing QC report tool Fig. 1 General overview schematic of the zebrafish gonad RNA-Seq experiment. RNA was isolated from pools of gonads dissected from 10 individuals and subjected to RNA-Seq analyses. The gonads were pooled according to sex (male, female and masculinised male), developmental period tested (40 dpf and 60 dpf) and treatment group (solvent control and 17α-methyltestosterone treatment). Gonad RNA pools were analysed in triplicate in CLC Genomics Workbench 6.0.2 software (Qiagen Bioinformatics GmbH, Hilden, Germany). CLC Genomics Workbench was used to trim low quality sequences (phred-based error probability threshold of 0.05) and to remove TruSeq adapter sequences from the raw reads. Trimmed reads from each pool were mapped and aligned to the Zebrafish Zv9 reference sequence (ftp://ftp.ensembl.org/pub/current_fasta/danio_rerio/dna/) from the ENSEMBL database. Annotation of the Zebrafish ZV9 reference genome was performed using the respective GTF file (ftp://ftp.ensembl.org/pub/current_gtf/danio_rerio/) using the CLC Genomics Workbench 6.0.2 ' Annotate with GFF/GTF' plug-in. In order to correct for differences in transcript length and library size, transcript abundance was normalised using the reads per kilobase per million mapped reads (RPKM) method [82]. Transcripts with RPKM value greater than or equal to 1.0 (RPKM ≥1.0) were regarded as expressed. Transcripts with RPKM value greater than or equal to 5.0 (RPKM ≥5.0) were regarded as reliably expressed.

Identification of differentially expressed transcripts
Differentially expressed transcripts were identified in the zebrafish gonad transcriptomes across different conditions (sex, age, MT treatment and gonad phenotype) using the proportions-based beta-binomial Baggerley's test [83]. The RPKM data was transformed by adding a constant (1.0) and then normalised using a quantile normalisation. Transcripts responsive to MT were identified via pairwise comparisons of MT-treated gonads with age-matched control ovaries and testes. Transcripts were regarded as differentially expressed if they complied with (1) the normalised fold change threshold of greater than or equal to two for upregulated transcripts and smaller than or equal to minus two for downregulated transcripts (≥ 2.0-fold) and (2) a false discovery rate (FDR) corrected p-value of ≤0.05 (p ≤ 0.05) [84].
To investigate the overall gene expression patterns across the treatment conditions and phenotypic sexes, we performed hierarchical clustering using Pearson correlation distance with complete linkage, and Principal Component Analyses (PCA), as implemented and visualised in CLC Genomics Workbench. The sets of differentially expressed transcripts that overlapped or differed between treatments were visualised using Venn diagrams plotted using the Venny software (http://bioinfogp.cnb. csic.es/tools/venny/).

Gene ontology enrichment analysis of differentially expressed genes
Gene Ontology Consortium (GO) [85] functional annotation terms and categories were assigned to the differentially expressed genes with the ' Add Annotations' tool of CLC Genomics Workbench. The GO terms, IDs and annotations were downloaded from the Gene Ontology database in the 20 May 2013 release of the ZFIN Zebrafish GO gene association file (http://www.geneontology.org/ GO.downloads.annotations.shtml).

Functional enrichment and network pathway analysis
Biological functions and metabolic pathways significantly overrepresented among the differentially expressed genes were identified using Metacore (GeneGo, Thomson Reuters, Carlsbad, California, USA). This approach utilizes Fisher's exact test with an FDR correction for multiple testing. Differentially expressed genes were defined as genes where (1) the absolute value for normalised fold change was greater than or equal to two and (2) the FDR adjusted p-value was less than 0.05 [84].
Validation of differentially expressed genes using quantitative RT-PCR Quantitative real-time PCR (qRT-PCR) was conducted on 12 differentially expressed transcripts identified from the RNA-seq data. TaqMan Gene Expression Assays (Applied Biosystems, Thermo Fisher Scientific, Waltham, Massachusetts, USA) specific to our genes of interest were used with a Stratagene Mx3000P (Agilent Technologies, Santa Clara, California, USA) Real Time-PCR thermal cycler. Each qPCR reaction was performed in triplicate. β-actin 1, ribosomal protein L13 alpha (rpl13α) and eukaryotic translation elongation factor 1 alpha 1, like 1 (eef1a1l1) were tested for their utility as the reference genes for this study. Previous studies have shown that the expression levels of these genes remain fairly constant across different developmental periods, treatment conditions, sexes and tissue types in zebrafish [86,87]. In this study, minor differences in β-actin 1 expression were found between different pooled ovary samples however the expression levels of rpl13α and eef1a1l1 were consistent across different tissue types (adult testes and ovaries) as previously reported [88] (data not shown). eef1a1l1 gave the highest average expression stability values with the geNorm algorithm [89].The comparative CT method (ΔΔCt) was used to determine relative gene expression compared to the reference gene eef1a1l1. The female ovary groups (40CO and 60CO) were used as the calibrator for calculation of relative expression. Relative expression was expressed as fold change (Fold Change = 2 − ΔΔCt). Complete details of the genes selected, TaqMan Gene Expression Assays and qPCR cycling conditions are provided in Additional file 1.

Sex ratios and gonad histology of control zebrafish
Sex ratios were determined via visualisation of EGFP fluorescence intensity of Tg(vas:egfp) zebrafish under a fluorescence microscope. Gonads exhibiting low and high EGFP fluorescence were classified as testes and ovaries respectively [69,73]. Male and female zebrafish were present in the control at 40 dpf and 60 dpf. The control sex ratios were 32% males at 40 dpf and 45% males at 60 dpf (Additional file 2). The mortality rates were very low (≤10%) for the control with no significant difference (p-value > 0.05) observed between replicates (Additional file 2).
To determine the effects of MT on zebrafish gonad development, we treated juvenile zebrafish with 100 ng/L MT from 20 dpf to 40 dpf and 60 dpf. All MT-treated zebrafish had testes except for one individual at 40 dpf (40 dpf: 59/60 individuals and 60 dpf: 98/98 individuals) (Additional file 2) suggesting that MT treatment had stimulated gonadal masculinisation.
The EGFP fluorescence signal was considerably fainter in 40 dpf testes than for 40 dpf ovaries. A few gonads (4/4 males) appeared to be undergoing gonadal transformation into testes at 40 dpf and were classified as presumptive testes (Fig. 2a). Several residual bodies indicative of apoptotic oocytes were present. The remaining oocytes in these gonads had lost their tight contact and showed signs of degeneration including shrinkage, irregularity of shape and basophilic cytoplasm. Stromal tissue had infiltrated the inter-oocyte spaces. In a few individuals, the germ cells had re-arranged to form tubule-like structures.
Regardless of the stage of ovarian development, intense EGFP fluorescence signals were detected in ovaries of 40 dpf and 60 dpf zebrafish. At 40 dpf, all control females analysed using histological examination possessed ovaries at early stages of development (14 /14 individuals). These immature ovaries were predominantly filled with tightly packed Stage IB oocytes although pre-meiotic Stage IA oogonia were located at the caudal edges of the ovaries (Fig. 2b). Most MT-treated zebrafish had mature testes at 40 dpf (17/22 individuals) (Fig. 2e). One was undergoing juvenile ovary-to-testis transformation (Fig. 2c). The remaining three individuals possessed immature testes (Fig. 2d). In contrast, control males had testes which had just completed juvenile ovaryto-testis transformation (4/4 individuals) (Fig. 2a).
Similar to 40 dpf control ovaries, the 40 dpf MTtreated ovary contained Stage IB oocytes (Fig. 2f ). Early range-finding experiments testing different periods of MT exposure suggested that occasionally a few individuals exposed to MT from 20 dpf to 40 dpf or 50 dpf were refractory to the masculinising effects of MT treatment (data not shown).
Testes at 60 dpf barely exhibited any EGFP fluorescence. There were three types of 60 dpf control testes: presumptive (2/10 males), immature (4/10 males) and mature testes (4/10 males) (Fig. 3ac). 60 dpf presumptive testes resembled 40 dpf presumptive testes which were populated with a few degenerative oocytes and large volumes of stromal tissue (Fig. 3a). Most 60 dpf control testes (8/10 individuals) were at a more advanced stage of differentiation than those at 40 dpf (Fig. 3b and c). Well-differentiated spermatogonial cysts containing spermatogonia, spermatocytes and spermatids were observed in 60 dpf control immature and mature testes ( Fig. 3b and c). The lumens of mature 60 dpf control testes were filled with mature spermatozoa (Fig. 3c).
Considerable individual differences in the extent of ovarian development were seen in 60 dpf control females. A few ovaries (6/13 females) consisted primarily of Stage IB oocytes (Fig. 3d). Others (7/13 females) contained oocytes at more advanced stages of development such as Stage II cortical alveolar oocytes ( Fig. 3e and f ).

Global sex differences in gonadal gene expression patterns
To elucidate the genes involved in zebrafish sexual differentiation, we compared the transcriptomes of developing zebrafish testes and ovaries using RNA Sequencing. Gonads were classified as testes or ovaries on the basis of EGFP expression in the juvenile Tg(vas:egfp) zebrafish as described above. Greater numbers of gene transcripts were present in the zebrafish testis than zebrafish ovary at both 40 dpf and 60 dpf (Additional file 3). About 12,000 transcripts were consistently expressed (RPKM ≥ 5) in zebrafish testes (Table 2). In contrast,~8000 transcripts were expressed in zebrafish ovaries at 40 dpf and 60 dpf (Table 2).
Developing zebrafish testes and ovaries showed striking differences in the numbers of differentially expressed transcripts at 40 dpf and 60 dpf. At 40 dpf, over 5000 transcripts were differentially expressed (≥ 2-fold, FDRadjusted p-value ≤0.05) between testes and ovaries ( Fig. 4a  and b). Of these, half were preferentially expressed in testes and half were expressed more highly in ovaries. Significantly greater numbers of sexually dimorphic transcripts were identified at 60 dpf than at 40 dpf. At 60 dpf, the number of sex differentially expressed transcripts had risen to almost 8500 ( Fig. 4c and d). Similar to what was observed at 40 dpf, the proportion of testis-enriched transcripts and ovary-enriched transcripts were approximately the same at 60 dpf.

Characterisation of the testicular differentiation programme
Transcripts with pro-male roles in fishes (dmrt1, amh, gsdf and sox9a) and genes encoding steroidogenic enzymes for 11-oxygenated androgen production (cyp11c1, cyp17a1, hsd11b2, nr5a1a and star) and androgen receptivity (ar) were more highly expressed in control testes than control ovaries during gonadal transformation and testicular differentiation (Table 3, Additional file 4). We also found upregulation of genes involved in gammaaminobutyric acid (gabrr2b and glsb) and dopamine (drd2l) signalling, immune response (il1b, mif and irf9) and encoding histone variants (h1f0, h1fx, h2afy2, histh1l and hist2h2l) in control testes (CT) compared to control There were developmental stage-specific differences in expression levels of testis-biased transcripts (Table 3).
Gata4, pdgfra and tp53 expression in testes were significantly higher in testes than ovaries at 40 dpf (early testicular development) but not 60 dpf (Table 3, Additional file 4). The magnitude of testis-biased expression of several pro-male genes (dmrt1, gsdf, wt1a, nr0b1 and amh) was significantly higher at 40 dpf than 60 dpf (Table 3, Additional file 4). Genes encoding sperm tail proteins septin3, septin8b and odf3b had been upregulated in testes compared to ovaries at 60 dpf but not 40 dpf (Table 3, Additional file 4, Additional file 6 ).

Characterisation of the ovarian differentiation program
Transcripts implicated in teleost fish female sex determination, ovarian development, oogenesis and folliculogenesis (foxl2, cyp19a1a, lhx8a, figla, bmp15, gdf9, zp21. and zp3b) were more highly expressed in control ovaries than control testes (Table 4, Additional file 4). The gene names and IDs can be found in Additional file 5.
We observed stage-dependent differences in ovarian expression patterns. Genes involved in fish folliculogenesis (bmp15, figla, gdf9, lhx8a, zp2.1, zp3b and vldr) had been more highly expressed in ovaries than testes from 40 dpf (Table 4). However, the expression levels of several pro-female transcripts (foxl2, cyp19a1a, esr2a and vtg1) had only been significantly upregulated in ovaries compared to testes at 60 dpf (Table 4). We also found higher expression of several genes involved in histone modification and DNA methylation (ehmt1a, ehmt2 and dnmt1) and encoding histone variants (h2afx and h1m) in ovaries compared to testes ( Table 4).
The expression patterns of a few transcripts with established roles in female sex determination and sex differentiation (sox9b, ctnnb1, ctnnbip1, wnt4a and rspo1) were not identified as sexually dimorphic in this study (Additional file 4).

Methyltestosterone induces masculinisation and de-feminisation of the gonadal transcriptome
Substantially fewer transcripts were differentially expressed between MT-treated testes and control testes in comparison with control ovaries, particularly at 60 dpf (Table 5). Although more than five thousand transcripts were differentially expressed in MT-treated gonads compared to control ovaries at 40 dpf (Table 5), only about a thousand differentially expressed transcripts separated the transcriptomes of MT-treated gonads from control testes transcriptomes ( Table 5). The transcripts differentially expressed in 40 dpf MT-treated testes compared to control ovaries overlapped considerably with those differentially expressed between control testes and control ovaries ( Fig. 4a and b).
At 60 dpf, the number of transcripts differentially expressed in MT-treated testes relative to control ovaries rose to seven and a half thousand (Table 5). At the same time, the number of transcripts differentially expressed between the MT-treated gonads and control testes fell to 20 (Table 5).
Marked changes indicating masculinisation and defeminisation of the gonad transcriptome were observed in testes treated with MT to 40 dpf and 60 dpf. Among the two and a half thousand transcripts with upregulated expression in 40 dpf MT-treated testes relative to 40 dpf control ovaries (Table 5), 61% were more highly expressed in control testes than control ovaries i.e., testis-enriched (Fig. 4a). Many of these transcripts have established roles in vertebrate testicular development and spermatogenesis ( Table 6, Additional file 7).
Of nearly three thousand transcripts showing downregulated expression in 40 dpf MT-treated testes compared to 40 dpf control ovaries (Table 4), 74% were more highly expressed in control ovaries than control testes i.e., ovaryenriched (Fig. 4b). Genes previously implicated in fish folliculogenesis were included among the transcripts more highly expressed in ovaries compared to MT-treated testes ( Table 7, Additional file 7). The proportions of transcripts with testis-enriched and ovary-enriched expression in 60 dpf MT-treated testes rose to 92% of the three and a half thousand upregulated transcripts ( Fig. 4c and Table 5) and 87% of the four thousand downregulated transcripts ( Fig.  4d and Table 5), respectively. This confirms that the correlation between the expression profiles of MT-treated gonads and control testes increased with developmental stage.

Methyltestosterone induces acceleration of spermatogenesis
The overlap between the testis-biased transcripts and transcripts showing upregulated expression in MTtreated testes relative to ovaries at 40 dpf was higher for 60 dpf control testes (96%) than 40 dpf control testes (61%) (Fig. 5a). A similar pattern was observed for the transcripts which exhibited downregulated expression in 40 dpf MT-treated testes. A greater proportion of transcripts exhibiting downregulated expression in 40 dpf MT-treated testes was shared with 60 dpf control testes (89%) than 40 dpf control testes (74%) (Fig. 5b). A comparison between the transcriptomes of 40 dpf MT-treated testes and 60 dpf control testes revealed a difference of just 136 differentially expressed transcripts (Additional file 8).
The pattern continued at 60 dpf where 92% of transcripts exhibiting upregulated expression in 60 dpf MT-treated testes overlapped with 60 dpf testis-biased transcripts compared to 55% for 40 dpf testis-biased transcripts (Fig. 5c). 87% of transcripts with downregulated expression in 60 dpf MT-treated testes compared to 60 dpf control ovaries was shared with 60 dpf control testis as compared to 56% for 40 dpf control testis (Fig. 5d).
The minor differences in the genetic cascades of 40 dpf MT-treated testes and 60 dpf control testes concurs with the gonad histological data which suggested that 17α-methyltestosterone had accelerated testicular development ( Fig. 2 and Fig. 3).

Methyltestosterone switches gonad transcription to the male program
The pathways enriched with genes exhibiting sexually dimorphic expression in zebrafish testes and ovaries were markedly different (Additional file 9). Similarly, the pathways overrepresented in MT-treated testes were considerably different from those in ovaries (Additional file 9). There was, however, substantial overlap between the pathways overrepresented in testes and MT-masculinized gonads, particularly at 60 dpf (Additional file 9). 'Cell cycle' , 'reproduction -progesterone-mediated oocyte maturation' and 'transcription' pathways were highly enriched in zebrafish ovaries during early ovarian development (40 dpf ) (Additional file 9). Zebrafish testes showed an enrichment of genes involved in 'cytoskeleton remodelling' , and 'immune response' , 'TGFβ signalling' and 'G-protein signalling' pathways at 40 dpf (Additional file 9).
While 40 dpf zebrafish ovaries were highly enriched in genes for 'cell cycle' , 'development' and 'transcription' pathways, MT-treated testes were enriched in genes involved in 'immune response' , 'cytoskeleton remodelling' and different 'cell cycle' pathways (Additional file 9).
Mitotic 'cell cycle' and nucleotide metabolism pathways were overrepresented in 40 dpf MT-treated testes (Additional file 9). In contrast, there was significant overrepresentation of pathways involved in 'apoptosis' , Brca1 mediated regulation of transcription, the synthesis of cortisone, cortisol and androgen signalling in 40 dpf control testes (Additional file 9).
To validate the RNA-Seq results, qPCR was conducted on 12 differentially expressed genes from four pairwise comparisons based on sex and treatment: 40CT vs. 40CO, 40MT vs. 40CO, 60CT vs. 60CO and 60MT vs. 60CO (Additional file 10). The direction of fold-change for the two platforms agreed for most of the genes (10/ 12 genes tested) (Additional file 10) confirming the accuracy of the RNA-Seq data.

Discussion
In this study, we investigated the effects of MT treatment on zebrafish gonadal transcriptomes during gonadal differentiation. We found that MT treatment masculinises gonadal transcriptomes and is accompanied by morphological changes consistent with male gonad development. MT activated pro-male gene expression (dmrt1, amh and gsdf) and induced steroidogenic enzymes required for production of 11-oxygenated androgens (cyp11c1 and hsd11b2). MT treatment also repressed gonadal expression of pro-female genes cyp19a1a, foxl2, bmp15 and gdf9, particularly during folliculogenesis at 60 dpf. Taken together, our data suggest that MT-induced masculinisation involves both activation of androgen receptor-mediated pathways and inhibition of aromatase.

MT activates pro-male gene expression
We found strong upregulation of dmrt1 expression in both control and MT-treated testes in this study, consistent with its essential role in zebrafish testis development [90]. In tilapia, upregulation of dmrt1 expression is associated with androgen-induced gonadal masculinisation [91]. Similarly, dmrt1 expression increased following MT treatment of XX medaka [92], and upon temperaturedependent induction of gonadal masculinisation in European seabass [93], tilapia [94], medaka [95] and pejerrey [96]. We suggest that the activation of dmrt1 may be one of the key steps underpinning hormonal and environmental gonadal masculinisation in fish.
MT treatment increased amh expression in zebrafish gonads. Amh operates downstream of dmrt1 in zebrafish [90] and is highly expressed in zebrafish transforming testes where it may inhibit cyp19a1a expression [13,97]. Defective amh/amhr signalling in medaka causes excessive germ cell proliferation and male-to-female sex reversal in half of XY hotei mutants [98,99]. The data are collectively consistent with the idea that amh signalling operates downstream of MT to repress female patterns of germ cell proliferation in teleosts.
MT treatment upregulated gsdf gonadal expression in zebrafish, consistent with a recent study in medaka demonstrating that MT induced gsdf expression in gonads [92]. In medaka, gsdf is involved in early testicular differentiation and is a key regulator of male-specific gene expression [100][101][102][103]. Gsdf expression has not been characterised in zebrafish. Our results suggest gsdf may be important for zebrafish testicular differentiation.
Interestingly, we found higher expression of gabrr2b and drd2l in control and MT-treated testes than control ovaries. Gamma-aminobutyric acid (GABA) and dopamine are important regulatory neurotransmitters in the hypothalamus-pituitary-gonadal axis. In teleost fish, both GABA and dopamine modulate production of Table 5 Identification of genes differentially expressed (≥ 2-fold, FDR-adjusted p-value ≤0.05) in the gonads in response to methyltestosterone treatment at 40 dpf and 60 dpf using Baggerley's test [83]. CO  hypothalamic gonadotrophin releasing hormone and pituitary gonadotrophins [104][105][106][107][108][109][110]. GABA A receptor has been implicated in Leydig cell proliferation in mice [111,112]. GABA may regulate androgen production in rodent testes [113][114][115]. Dopamine  type 2 receptor has been reported in mammalian testes and male germ cells where it regulates sperm capacitation and motility [116,117]. Although dopamine antagonists are used with gonadotrophin analogues (eg. Ovaprim) for induced spawning in aquaculture [118,119], the role of dopamine signalling in fish spermatogenesis has not been studied. Our expression data suggests that GABAergic and dopamine signalling may be important for teleost testicular development.
Genes encoding the cytokines mif and il1b had been upregulated in control and MT-testes compared to control ovaries. GO pathways involved in immune response pathways had been significantly upregulated in control testes and MT-treated testes compared to control ovaries. MIF-JAB1 and interleukin 1 signalling pathways were significantly over-represented in testes. MIF and IL-1β are cytokines involved in innate immunity and inflammation. Both act as paracrine factors responsible for regulation of testosterone production by Leydig cells in rat testes [120][121][122][123][124]. MIF is produced by Leydig cells and Sertoli cells where it is involved in the cross-talk between Leydig cells and testicular seminiferous tubule somatic cells, spermagonial cell migration and may be involved in regulation of sex hormone production and spermiogenesis [125][126][127][128]. IL-1β is expressed in mice testicular germ and somatic cells [129]. Interleukin 1 is Fig. 5 Comparison of differentially expressed genes in MT-treated zebrafish testes with testis-enriched genes of 40 dpf and 60 dpf control testes. n = 3 pools of 10 individuals each with the exception of 40MT which comprised of 2 pools. a Venn diagram of numbers of genes which showed upregulated expression in 40 dpf MT-treated testes (40MT) overlapped with numbers of genes upregulated in 40 dpf testes (40CT) and 60 dpf testes (60CT). b Venn diagram of numbers of genes which showed downregulated expression in 40 dpf MT-treated testes (40MT) overlapped with numbers of genes downregulated in 40 dpf testes (40CT) and 60 dpf testes (60CT). c Venn diagram of numbers of genes which showed upregulated expression in 60 dpf MT-treated testes (60MT) overlapped with numbers of genes upregulated in 40 dpf testes (40CT) and 60 dpf testes (60CT). d Venn diagram of numbers of genes which showed downregulated expression in 60 dpf MT-treated testes (60MT) overlapped with numbers of genes downregulated in 40 dpf testes (40CT) and 60 dpf testes (60CT) important for Sertoli cell proliferation [130,131]. Our data suggests that mif and il1b may mediate androgeninduced sex reversal in zebrafish.
Unexpectedly, we found that a few irf genes which included irf9 were more highly expressed in control and MT-treated testes than control ovaries. A truncated form of irf9 was identified as the master sex determining gene in rainbow trout [132]. Our expression data supports Yano et al. 2014's hypothesis that interferon signalling may be involved in teleost testicular development and spermatogenesis. It may be worthwhile to study possible roles of interferon in teleost testes.

MT represses pro-female gene expression
MT treatment caused strong inhibition of bmp15 expression in zebrafish gonads. Bmp15 is necessary for activation of cyp19a1a expression and therefore estrogen production in granulosa cells surrounding Stage II oocytes [142]. Loss of cyp19a1a expression in juvenile zebrafish ovaries leads to germ cell apoptosis and female-to-male sex reversal [70]. Our results suggest that MT treatment may suppress ovarian development by preventing bmp15-mediated cyp19a1a expression in granulosa cells. MT treatment also strongly inhibited genes involved in Wnt signalling pathways. Wnt signalling is important for ovarian differentiation in zebrafish [143], medaka [144,145] and rainbow trout [146] so its repression is consistent with masculinisation.
We found that MT downregulated expression of genes involved in chromatin histone modification and DNA methylation pathways in epigenetics (dnmt1, hdac11, ehmt1a and ehmt2) compared to control ovaries. Interestingly, hdac11 and ehmt2 expression had been upregulated with heat-induced masculinisation in European seabass [93]. The expression levels of the histone variants h2afx and h1m had been lower in MT-treated testes than ovaries, similar to reported female-biased expression in Olive flounder [147].

How do sex steroid hormones cause gonadal masculinisation?
Consistent with observations during MT-induced gonadal masculinisation in tilapia [40], we found cyp11a1, hsd3b1 and cyp19a1a transcripts were reduced when juvenile zebrafish were treated with MT. We also found that MT treatment actively promoted cyp11c1 and hsd11b2 expression in zebrafish gonads. Cyp11c1 works together with hsd11b2 to convert androgens into 11-oxygenated androgens. This contrasts with studies in rainbow trout which found downregulated cyp11c1 expression during androgen-induced gonadal masculinisation using 11βhydroxyandrostenedione [41,45]. Despite both being 11-KT precursors, MT and 11β-hydroxyandrostenedione produce different intermediate products with different androgenic potencies. It is possible that MT induces gonadal masculinisation by different mechanisms to those employed by 11β-hydroxyandrostenedione, perhaps via these differing intermediate products.
Our results suggest that androgen production may also be important for natural juvenile ovary-to-testis gonadal transformation in zebrafish. Cyp11c1 and hsd11b2 expression was higher in transforming testes than juvenile ovaries during juvenile ovary-to-testis transformation at 40 dpf. This concurs with the higher levels of cyp11c1 expression previously reported in juvenile ovotestes compared to juvenile ovaries amid zebrafish juvenile ovary-to-testis transformation [13,143]. This indicates that the steroidogenic enzymes required for androgen production are active during early stages of juvenile ovary-to-testis transformation. We found that pathways involved in cortisol biosynthesis from cholesterol and cortisone biosynthesis and metabolism and androgen receptor nuclear signalling were higher in transforming testes than MT-treated testes at 40 dpf. Cortisol and androgen signalling pathways were shown to be linked in high temperature-induced masculinisation of pejerrey via hsd11b2 [148]. This suggests that hsd11b2-mediated androgen production may be important for testicular differentiation in zebrafish.
Taken together, our data show that increased 11oxygenated androgen production and decreased aromatase expression may be equally important for MT-induced gonadal masculinisation in zebrafish.

MT alters germ cell proliferation rates in juvenile ovaries
The juvenile ovary-to-testis transformation process in zebrafish is characterised by oocyte apoptosis [71]. In zebrafish, genes in pro-apoptotic pathways (tp53) have been implicated in testicular differentiation [149], while anti-apoptotic pathways (NF-κB) promote ovarian differentiation [150]. Partial or complete germ cell depletion in gonads invariably leads to gonadal masculinisation [88,149,151,152].
Gonadal hormones can regulate germ cell apoptosis in vertebrates [153,154] including teleost fish [5]. In zebrafish, 17α-ethinylestradiol and fadrozole can alter germ cell apoptosis and proliferation during gonad differentiation [54]. Several genes that regulate cell apoptosis and proliferation are also involved in gonadal sexual differentiation [155]. A threshold of germ cells or meiotic oocytes is required for ovarian development in zebrafish [152,156] such that the phenotypic sex in zebrafish depends critically on the number of germ cells in gonads during early development [88,152,157]. How germ cell numbers maintain female identity is not fully understood.
We observed that MT treatment downregulates genes involved in DNA replication in early S phase, Estrogen Receptor 1 (Esr1) regulation of G1/S phase, mitogenic action of Estradiol /Esr1 and ligand-dependent activation of the Estrogen receptor 1/ Specificity Protein (Esr1/Sp) pathway in 40 dpf juvenile ovaries. We suggest that androgens may induce female-to-male sex reversal via modulation of germ cell proliferation rates. MT-induced repression of estrogen-activated mitotic pathways in juvenile ovaries leads us to hypothesise that MT may trigger female-to-male sex reversal via inhibition of estrogen-activated proliferation of oocytes in zebrafish ovaries. This would reduce oocyte numbers below the threshold, leading to activation of male-specific gene expression in gonadal somatic cells.
In this study, we used transgenic zebrafish derived from domesticated zebrafish strains. Domesticated zebrafish strains use polygenic sex determination [58,59] and are sensitive to environmental and hormonal factors [32,53,61,63,158]. The sensitivity of domesticated zebrafish strains to hormonal effects on gonad development render them good models for determining the mode of action of androgens on induction of testicular differentiation. Evidence comparing Atlantic Silverside utilising genetic sex determination (GSD) and temperature sex determination (TSD) suggest that expression of cyp19a1a, a key gene in ovarian development, was earlier, stronger, more sexually dimorphic and less temperature sensitive in GSD strains than TSD strains [159]. It is possible that the thresholds of environmental and hormonal sensitivity for key sex determining loci may vary between wild (ZZ/ZW sex determination) and domesticated (polygenic sex determination) zebrafish strains. Different mechanisms may be activated for androgen-induced sex reversal of wild zebrafish populations. Future studies using wild strains would facilitate elucidation of the pathways involved. Due to the polygenic mode of sex determination and lack of sex-linked markers in domesticated zebrafish strains, we were unable to distinguish between male and female zebrafish prior to morphological gonadal differentiation. Therefore, the MT treatment was conducted on zebrafish with unknown proportions of male and female sexual genotypes. This resulted in production of populations comprising males and neo-males. The development and use of all-male and allfemale zebrafish lines [160] in future studies would help further clarify whether MT treatment has different effects on genetic males and females.