# The transcriptomic signature of different sexes in two protogynous hermaphrodites: Insights into the molecular network underlying sex phenotype in fish

A. Tsakogiannis, T. Manousaki, Jacques Lagnel, A. Sterioti, M. Pavlidis, N.

Papandroulakis, C. C. Mylonas, C. S. Tsigenopoulos

► To cite this version:

A. Tsakogiannis, T. Manousaki, Jacques Lagnel, A. Sterioti, M. Pavlidis, et al.. The transcriptomic signature of different sexes in two protogynous hermaphrodites: Insights into the molecular network underlying sex phenotype in fish. Scientific Reports, 2018, 8 (1), pp.1-17. 10.1038/s41598-018-21992-9. hal-02620154

HAL Id: hal-02620154

https://hal.inrae.fr/hal-02620154

Submitted on 25 May 2020

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
SCIENTIFIC REPORTS

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

OPEN

# The transcriptomic signature of different sexes in two protogynous hermaphrodites: Insights into the molecular network underlying sex phenotype in fish

Received: 17 August 2017

Accepted: 14 February 2018

Published online: 23 February 2018

A. Tsakogiannis<fnref n="1" />,², T. Manousaki<fnref n="1" />, J. Lagnel<fnref n="1" />, A. Sterioti<fnref n="1" />, M. Pavlidis², N. Papandroulakis<fnref n="1" />, C. C. Mylonas<fnref n="1" /> & C. S. Tsigenopoulos<fnref n="1" />

Sex differentiation is a puzzling problem in fish due to the variety of reproductive systems and the flexibility of their sex determination mechanisms. The Sparidae, a teleost family, reflects this remarkable diversity of sexual mechanisms found in fish. Our aim was to capture the transcriptomic signature of different sexes in two protogynous hermaphrodite sparids, the common pandora *Pagellus erythrinus* and the red porgy *Pagrus pagrus* in order to shed light on the molecular network contributing to either the female or the male phenotype in these organisms. Through RNA sequencing, we investigated sex-specific differences in gene expression in both species' brains and gonads. The analysis revealed common male and female specific genes/pathways between these protogynous fish. Whereas limited sex differences found in the brain indicate a sexually plastic tissue, in contrast, the great amount of sex-biased genes observed in gonads reflects the functional divergence of the transformed tissue to either its male or female character. A common "crew" of well-known molecular players is acting to preserve either sex identity of the gonad in these fish. Lastly, this study lays the ground for a deeper understanding of the complex process of sex differentiation in two species with an evolutionary significant reproductive system.

The differences between sexes have always fascinated, yet troubled philosophers and scientists<fnref n="1" />. Unlike the commonality of sexual reproduction, the ways that species express their sexuality vary considerably in the animal kingdom. Among the animals that reproduce sexually, teleost fishes show a very wide repertoire of reproductive modes, including all known reproductive styles found in vertebrates, from gonochorism (i.e., after sex is determined, it remains unchanged for the entire life cycle) to unisexuality (all-female species) and hermaphroditism (sequential, serial and simultaneous, including outcrossing and selfing species)²,³. These diverse sexual phenotypes in fish are regulated by a variety of sex determination (SD) mechanisms, along a continuum of environmental and genetic factors³,⁴.

The knowledge concerning sex differences at the molecular level results mainly from studies in humans and model vertebrates such as mouse, chicken and African clawed frog⁵⁻⁸. In fish, only few genes have been reported so far to play the role of master regulators in the SD systems of an investigated species. These include *DmY/* *Dmrt1bY* in *Oryzias latipes* and in *Oryzias curvinotus*⁹, *Gsdf* in *Oryzias luzonensis*¹⁰, *Amhy* in *Odontesthes hatcheri*¹¹, *Amhr2* in *Takifugu rubripes*, *T. pardalis* and *T. poecilonotus*¹² and *SdY* in the salmonid family¹³. Recently, a distant cis-regulatory element of *Sox3* necessary for male determination in *Oryzias dancena* (XX/XY SD system) has also been identified¹⁴ and *Dmrt1* has been suggested as the sex-determining gene (SDG) in *Cynoglossus semi-laevis* (ZZ/ZW SD system)¹⁵. To date, studies regarding the sex-specific differences in gene expression have been conducted mainly in SD systems of model fish species that are well characterized at the genomic level, with distinguishable heteromorphic sex chromosomes, exhibiting Genetic Sex Determination (GSD) and gonochorism¹⁶⁻¹⁸.

<footnote n="1">¹Institute of Marine Biology, Biotechnology and Aquaculture (IMBBC), Hellenic Centre for Marine Research (H.C.M.R.), Heraklion, Greece. ²Department of Biology, University of Crete, Heraklion, Greece. Correspondence and requests for materials should be addressed to C.S.T. (email: tsigeno@hcmr.gr)</footnote>
The flexibility governing reproductive mechanisms in fish, with the presence of different systems among even phylogenetically close species makes the picture more complicated. However, sex-linked markers and/or sex-determining regions have already been reported in several species, including some that are of interest to the aquaculture industry, such as the catfish *Clarias gariepinus*¹⁹, the Nile tilapia *Oreochromis niloticus*²⁰, the rainbow trout *Oncorhynchus mykiss*²¹, the turbot *Scophthalmus maximus*²², the half-smooth tongue sole *Cynoglossus semi-laevis*²³, some salmonids²⁴⁻²⁶, Gasterosteidae²⁷, Eigenmannia²⁸ and the European sea bass *Dicentrarchus labrax*²⁹.

Among teleosts, the Sparidae family is considered to be one of the most diversified families regarding reproductive systems, exhibiting nearly all sexual modes known in fish³⁰. In addition to gonochoristic species, various forms of hermaphroditism have been recorded³¹. Therefore, the Sparidae family is a unique model for comparative studies to understand the molecular mechanisms underlying different sexual motifs. Until recently, genetic studies on sex determination/differentiation and sex change in sparids were conducted mainly on the protandrous black porgy *Acanthopagrus schlegelii*³²,³³ and some QTL studies have been carried out in gilthead seabream *Sparus aurata*³⁴. The first work on sex-specific gene expression in a spard was on a rudimentary hermaphrodite sharpnout seabream *Diplodus puntazzo*³⁵. Despite the fact that hermaphroditism is widespread in many fish lineages³⁶, studies on sex-related transcriptomic differences have been conducted in only a few hermaphrodites: the protandrous Asian seabass *Lates calcarifer*³⁷, the protogynous bluehead wrasse *Thalassoma bifasciatum*³⁸ and, very recently the protogynous clownfish *Amphiprion bicinctus*³⁹.

The common pandora *Pagellus erythrinus* and the red porgy *Pagrus pagrus* are two other sparids exhibiting protogynous hermaphroditism. They are both benthopelagic fish, distributed through the Eastern Atlantic Ocean along the west coast of Europe and Africa (from Scandinavia to Angola), the entire Mediterranean Sea and the Black Sea. The red porgy's distribution extends also to the western Atlantic, from the northern USA to Argentina. Typically, both fishes mature first as females, changing to males after two years of age or when attaining a body length of 17–18 cm for common pandora and after three years of age or about 24 cm of body length for red porgy⁴⁰,⁴¹. However, sexual patterns appear to be far more complex. Sex change and the sexual structure of populations of sequential hermaphrodite fishes may be genetically controlled or triggered by environmental/external factors, such as behavioral or/and demographic changes in their social systems⁴², or even a combination of both⁴³. According to Cataudella et al.⁴⁴, most of the Mediterranean sparids have a haploid chromosome number of 24 and recently a linkage map, consisting of 24 LGs, was built for the common pandora⁴⁵. To our knowledge, there is no evidence concerning the presence of heteromorphic sex chromosomes in these spard species or whether their sex determination system is genetic, environmental or an amalgam of both. Because of the intriguing reproductive mode of common pandora and the red porgy, they are important to be studied from an evolutionary point of view, as protogyny has been reported as the most common model of hermaphroditism in teleosts⁴⁶⁻⁴⁸. Moreover, their potential for large-scale aquaculture covers the need for species diversification in the sector⁴⁰.

This study presents the transcriptomic signature of different sexes in these two spard fishes, which share the same protogynous reproductive mode. Using RNA sequencing⁴⁹, we analyzed and compared differences between male and female transcriptomes of common pandora and red porgy. In particular, we investigated gene expression differences between male and female gonads and brains, the two most important tissues in sexual development and reproduction⁵⁰. We obtained a global view of sex-biased expression in these tissues and through comparative analysis we also revealed common male and female specific genes/pathways that are implicated in sex differentiation and sex maintenance. We also carefully investigated the presence of genes reported to be involved in sex determination/differentiation mechanisms in other vertebrates and fish and compared their expression patterns in the two species under study. Our results contribute to the understanding of the complex processes of sex differentiation in these protogynous hermaphrodites.

## Results

**Pre-processing treatment of reads and assembly construction.** Illumina HiSeq 2500 sequencing of red porgy and common pandora samples yielded 685,707,304 (342,853,652 read pairs) and 644,189,592 (322,094,796 read pairs) reads, respectively. The raw sequence data have been submitted to the National Centre for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database (BioProject ID: PRJNA395994).

For red porgy, pre-processing of reads resulted in 247,712,854 high-quality reads used for the initial transcriptome assembly, which consisted of 295,367 putative transcripts (N50: 1,961 nucleotides, mean length: 928 nucleotides) and showed significant similarity to 23,122 out of 26,754 European sea bass genes. After transcript quality assessment, our initial dataset was limited to 98,012 contigs (N50: 2,978; mean length: 1,587 nucleotides) corresponding to 67,132 genes and was significantly similar to 20,130 European sea bass genes. Filtering resulted in elimination of 197,355 possibly spurious transcripts and this restricted dataset constituted the final transcriptome with only slightly reduced numbers of hits to the unique genes of European sea bass.

For common pandora, filtering of reads resulted in 232,768,947 high-quality reads that were used to construct an initial assembly of 396,201 putative transcripts (N50: 1,529 nucleotides, mean length: 824 nucleotides). Additional editing excluded 254,892 possibly spurious transcripts, thus limiting our dataset to 141,309 contigs (N50: 2,450; mean length: 1,254 nucleotides), corresponding to 98,211 genes. On this occasion, the unfiltered transcript dataset showed significant similarity to 23,631 European sea bass genes out of 26,754, while the filtered assembly showed significant similarity to 21,028 genes. Following the same reasoning, the restricted dataset constituted the final transcriptome (Fig. 1).

To evaluate the accuracy of the assembled sequences (transcripts), all the usable sequencing reads were aligned onto the transcripts using Bowtie²⁵. A high percentage of red porgy and common pandora reads (84% and 81%, respectively) were successfully back-mapped on each assembly. Similarities between assemblies and UniProtKB/SwissProt database disclosed 21,217 and 21,251 contigs for red porgy and common pandora, respectively, with coverage higher than 80%. In addition, BUSCO run (Supplementary Table ST1) indicated that the assemblies
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The pipeline followed to build assemblies and differential expression profiles in brains and gonads of two protogynous sparids. Flow chart of the basic steps implemented from raw reads for the selection of final assembled loci and differentially expressed genes in brains and gonads of both species.</figcaption></figure>

were nearly complete with 86% and 82% complete matches of vertebrate orthologs, found in red porgy and common pandora assemblies, respectively.

**Transcriptome annotation.** The blastx search in red porgy and common pandora assembled transcriptomes ($e$-value cut-off $<10^{-5}$) returned 43,123 and 50,106 contigs, with at least one blast hit against Uniprot⁵² and 38,354 and 44,392 contigs with a protein domain match against InterProScan⁵³ search, respectively (Supplementary Table ST2). To identify the possible functions, Gene Ontology (GO) assignments were used to classify the sequences. Total numbers of 38,265 red porgy contigs and 51,026 common pandora contigs with at least one GO term were found (Blast2Go⁵⁴ annotated) (Supplementary Table ST2). The sequences from both species were categorized to thirteen molecular function (MF), 22 biological process (BP) and thirteen cellular component (CC) gene categories in GO level 2 (general function categories) (Supplementary Figures SF1-4). To assess the functional diversity of the assembled transcriptome, GO annotations of red porgy were compared with those of common pandora transcriptome, reflecting a similar functional distribution on GO categories and similar sequence diversity of the two transcriptomes. Profiling of both species (whole transcriptome) was generally very much alike, as shown in the comparative view produced by WEGO⁵⁵ (Fig. 2).

The potential enzymes were characterized based on the chemical reaction they catalyze, using the prediction of Enzyme Commission (EC) numbers for each sequence. For red porgy, among the 38,265 blast2go-annotated contigs, 9,572 contigs were classified to 1,194 different EC numbers, while among the 51,026 common pandora blast2go-annotated contigs, 12,821 were classified to 1,243 different EC numbers (Supplementary Table ST2).

**Global expression pattern of brains and gonads.** The principal component analysis (PCA) plots obtained for the global gene expression pattern observed in brain and gonad tissues of the nine red porgy and ten common pandora individuals showed the formation of three separate clusters for both species: the brains, the male gonads (testes) and the female gonads (ovaries) (Fig. 3). While testes and ovaries showed a clear separation with a remarkable variation within the two groups, male or female brains from both species were clustered together with very similar expression patterns.

Tissue-specific gene ontology profiling for both fish revealed differences in GO terms representation between brains and gonads, in all three categories in GO level 2. When comparing brains and gonads, except from altered percentages of the same main ontology categories (shared between brains and gonads), there were terms such as synapse (CC), cell junction (CC), behavior (BP), rhythmic process (BP), cell aggregation (BP), cell killing (BP), chemoattractant activity (MF) and chemorepellent activity (MF) that were present only in the brains of the two sparids, but were totally absent or present in low numbers in their gonads. This would indicate brain specificity of the underlying genes (Supplementary Figures SF1-4).

For a better understanding of the biological processes and molecular functions underlying the nature of a male or female phenotype, GO term enrichment analysis was carried out in order to detect differentially-expressed (DE) genes over-represented in the different sexes of both species. Here, the term 'enrichment' stands for the comparison between the DE genes of the female (test-set) against the male (reference-set) either for brains or gonads. The same picture of sample clustering in both species was also shown through heatmap plots (Fig. 4). Comparisons between males and females were conducted separately for the brain and gonad samples to produce PCA plots and heatmaps (Supplementary Figures SF5 and SF6). This closer view again showed that sex-dependent expression differences are greater by far in gonads than in brains. In the brain PCA plots of both species, males
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Comparative view of GO annotation. Comparison of Gene Ontology. (GO) categories (assigned by WEGO) between red porgy and common pandora. GO terms for the genes of red porgy (dark purple) and common pandora (pink) were categorized into cellular components, molecular functions and biological processes.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Principal component analysis of global expression profiles in two protogynous sparids. Principal component analysis of overall expression profiles of the gonad and brain samples of male and female individuals (a) Red porgy (PC1 = 49% variance and PC2 = 27% variance), (b) Common pandora (PC1 = 44% variance and PC2 = 22% variance). Blue = males, Red = females and Triangles = Gonads, Circles = Brains). All genes with non-zero total read count in the count matrix used for the PCA plot (a filtered matrix).</figcaption></figure>

and females clustered together, with red porgy brain samples clustered more tightly than those of the common pandora that show a looser clustering.

**Comparing gene expression in male and female brains.** Out of 67,132 red porgy and 98,211 common pandora total loci, 43,330 and 61,048 genes had at least 15 reads in all compared samples (the sum of counts in
<figure><img src="image_6.png" /><figcaption>b. Common pandora</figcaption></figure>

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

**Figure 4.** Heatmap plots for global expression profiles in two protogynous sparids. Heatmaps of the variance stabilized, transformed count data for male and female brains and gonads. (a) Red porgy, (b) Common pandora samples clustered according to their global gene expression patterns. Only the top 500 variable genes were used to produce the heatmaps for visibility reasons.

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

**Figure 5.** Comparative view of the number of genes differently expressed in brains and gonads of males and females of two protogynous sparids. Venn-like diagram showing the genes differentially expressed in brains and gonads of male and female red porgy and common pandora. Light-Blue: male-biased genes. Red: female-biased genes. Dark-Purple (intersection): un-biased genes.

each row of the count matrix >15) and were used for differential expression analysis. In general, brain DE genes showed low counts and low fold-change differences (Supplementary Figure SF7). For red porgy brain samples, 553 genes differed significantly with 365 of them having a male-biased (over-expression) and 188 of them having female-biased expression. Common pandora brain samples showed higher numbers of DE genes and an opposite direction of expression than red porgy brain samples. In total, 4,910 genes showed differential expression between the sexes, of which 1,733 were up-regulated in the male brains and 3,177 of them were up-regulated in the female brains (Fig. 5).

A GO analysis was conducted for all three general function categories (BP, MF, CC). Focusing on biological processes, brain enrichment analysis resulted in 241 GOs over-represented in female (light red-colored in Supplementary Excel Table SET1) and 265 in male common pandora (light blue-colored in Supplementary Excel Table SET1). For red porgy, 58 GOs were over-represented in females (light red-colored in Supplementary Excel Table SET2) and 113 in males (light blue-colored in Supplementary Excel Table SET2).

**Comparing gene expression in male and female gonads.** In gonad samples, we applied the same filter and kept only loci with a sum of expected counts in each row of the count matrix >15. Specifically, 43,692
and 52,126 genes fulfilled this criterion and were examined for differential expression in red porgy and common pandora gonads, respectively. The analysis revealed striking differentiation between testes and ovaries. In total, 18,724 red porgy gonadal genes showed a significant differential expression, of which 15,989 were up-regulated in the testes and 9,951 of them were up-regulated in the ovaries. In the case of common pandora, 25,928 gonadal genes had significantly differential expression, 16,148 them with up-regulation in their testes and 9,780 them with up-regulation in their ovaries (Fig. 5).

Focusing again on biological processes, Enrichment Analysis of DE genes in common pandora gonad tissues resulted in 236 GO terms significantly enriched in testes (light blue-colored in Supplementary Table SET3) while 517 GO terms were found in ovaries (light red-colored in Supplementary Table SET3); likewise, in red porgy 284 GO terms were significantly enriched in testes (light blue-colored in Supplementary Table SET4) and ovaries enriched with 370 GO terms (light red-colored in Supplementary Table SET4).

**Sex-specific candidate genes in gonads.** To have a closer view in the molecular playground causing sex dimorphism in these sparids, we checked the presence and the expression patterns in their gonads, of genes that have been reported as having an important role or participating in sex determination/differentiation or in sex maintenance in other fish or vertebrates. We selected 114 candidate genes in total. Out of this list, 84 genes were present in the transcriptomes of both species, while *Srd5a2* was present only in red porgy transcriptome and *Dmrtpa* was found only in common pandora. Table 1 includes all the genes with significantly differential expression in at least one of the studied species (excluding the not significantly different -NSD- characterized- genes in both species) and compares their expression pattern between them.

**Targeted search of candidate genes for sex determination/differentiation.** To identify genes that play a role in reproduction and more specifically in sex determination and differentiation of fish, we searched the annotated brain and gonad datasets for these GO categories and checked if they were differentially expressed. In particular, DE genes related to broad terms of biological processes, such as *reproduction* and *reproductive process* were searched in gonads and brains of both sparids. A large number of genes related to *reproduction* and a few to *reproductive process* were found in the gonads of both species. In the common pandora, 492 female-biased and 717 male-biased genes were found under the *reproduction* GO term, while thirteen female-biased and eleven male-biased genes were found under the *reproductive process* GO term. In the red porgy, 467 female-biased and 430 male-biased genes belonging to *reproduction* and nine female-biased and twelve male-biased genes belonging to *reproductive process* were detected. The same search in brains retrieved 246 *reproduction* related genes expressed in the common pandora, of which 153 were male-biased and 93 female-biased. Only four male-biased *reproductive process* related genes were found to be significantly expressed in common pandora brains. In the red porgy, 22 male-biased and just three female-biased *reproduction* related genes were found and no *reproductive process* related genes were among the DE genes. Further data mining, focusing on more specific GO terms, revealed nine loci in common pandora and five loci in red porgy containing the GO term *sex determination* and displayed significant differential expression in male and female gonads. A search for loci implicated in *sex differentiation* with differential expression in male and female gonads retrieved fourteen loci for common pandora, while 22 were found for red porgy. In brains, only one *sex determination* related and two *sex differentiation* related loci were found to be differentially expressed in the common pandora, while no differentially expressed loci for these GO terms were found in red porgy (Supplementary Excel Tables SET5–7). Finally, several genes expressed in male brains characterized with the GO term *social behavior* enriched in the red porgy and *inter-male aggressive behavior* and *male courtship behavior* in common pandora (Supplementary Excel Table SET8).

## Discussion

The present study investigated sex-specific gene-expression patterns in gonads and brains of two protogynous hermaphrodites belonging to the Sparidae family, the common pandora and the red porgy. This is the first comparative transcriptomic analysis of two species displaying the same reproductive mode, focusing on genes that are differentially expressed between sexes to pinpoint a common network that forms and/or maintains the male and female phenotype. In both species, male and female brains showed far fewer differences in expression compared to the gonads. Male and female brains were clustered together, forming a dense cluster. Quite the opposite was observed for the gonads, which differ greatly between the two sexes. Testes and ovaries were well separated and there was great variance between the formed groups (Figs 3, 4 and Supplementary Figure SF5).

**Minor sex differences indicated a more homogeneous and sexually plastic brain.** Our study revealed subtle differences in expression between the sexes in the brain of the two protogynous species. This seems to be a common pattern in tissue-based studies, where sex-biased expression tended to be highest in the gonads compared to other tissues⁵⁶. A low degree of differential expression in the brain has also been reported in other fish studies³⁵,³⁹,⁵⁷. Nevertheless, in the current study we examined the brain as a whole, whereas it has been recognized that sex-specific differences in fish are noticeable only when comparing certain regions of the brain⁵⁸. Despite the overall limited differences observed for brains compared to the gonads, common pandora male and female brains demonstrated a greater divergence than those of the red porgy, with many more genes differentially expressed between the sexes. Moreover, the direction of expression is female-biased in the common pandora, in contrast with the male-biased pattern observed here in the red porgy, but also in the sparid sharpsnout seabream³⁵ and other fish examined³⁸,⁵⁷. However, these results were in line with the expression pattern observed in protogynous clownfish brain³⁹.

A closer look into the specific genes that were up-regulated in male or female brains revealed numerous candidates with a well-known role in sex differentiation in other species. Zona pellucida (ZP) genes that showed a
<table><thead><tr><th rowspan="2">Gene name</th><th rowspan="2">Gene description</th><th colspan="2">Gene ID</th><th colspan="4">Expression</th></tr><tr><th>Common pandora</th><th>Red porgy</th><th colspan="2">Common pandora</th><th colspan="2">Red porgy</th></tr><tr><th></th><th></th><th></th><th></th><th>sex</th><th>log2FC</th><th>sex</th><th>log2FC</th></tr></thead><tbody><tr><td>Cyp19a1</td><td>Aromatase a (gonad isoform)</td><td>TR72019|c0_g1</td><td>TR50491|c0_g1</td><td>F</td><td>-8.12</td><td>F</td><td>-5.33</td></tr><tr><td>Gdf9</td><td>Growth and differentiation factor 9</td><td>TR70663|c0_g1</td><td>TR52139|c0_g1</td><td>F</td><td>-7.36</td><td>F</td><td>-8.35</td></tr><tr><td>Sox3</td><td>SRY-related HMG box 3</td><td>TR69218|c2_g2</td><td>TR64334|c1_g1</td><td>F</td><td>-6.91</td><td>F</td><td>-7.45</td></tr><tr><td>LHb</td><td>Luteinizing hormone, beta polypeptide</td><td>TR66625|c0_g1</td><td>TR40794|c1_g1</td><td>F</td><td>-6.54</td><td>F</td><td>-9.67</td></tr><tr><td>Foxl2</td><td>Forkhead box L2</td><td>TR59502|c0_g1</td><td>TR43506|c0_g1</td><td>F</td><td>-6.18</td><td>F</td><td>-7.59</td></tr><tr><td>Figla</td><td>Factor in the germline alpha</td><td>TR54099|c0_g1</td><td>TR49469|c0_g1</td><td>F</td><td>-5.77</td><td>F</td><td>-9.09</td></tr><tr><td>Wnt4a</td><td>Wingless-type MMTV integration site family, member 4a</td><td>TR67106|c0_g2</td><td>TR40636|c0_g1</td><td>F</td><td>-5.61</td><td>F</td><td>-7.02</td></tr><tr><td>Cyp26a1</td><td>Cytochrome P450, family 26, subfamily a, polypeptide 1</td><td>TR72448|c0_g1</td><td>TR49995|c0_g1</td><td>F</td><td>-4.60</td><td>F</td><td>-6.73</td></tr><tr><td>Dmrt2/Dmrt2a</td><td>Doublesex- and mab-3-related transcription factor 2(a)</td><td>TR78139|c1_g1</td><td>TR59371|c2_g1</td><td>F</td><td>-4.52</td><td>F</td><td>-7.92</td></tr><tr><td>Lhr</td><td>Luteinizing hormone receptor</td><td>TR89660|c1_g2</td><td>TR61712|c0_g1</td><td>F</td><td>-4.40</td><td>F</td><td>-3.84</td></tr><tr><td>Stra6</td><td>Stimulated by retinoic acid gene 6</td><td>TR86270|c0_g1</td><td>TR59700|c3_g1</td><td>F</td><td>-3.99</td><td>F</td><td>-5.45</td></tr><tr><td>Arb</td><td>Androgen receptor beta</td><td>TR85026|c0_g1</td><td>TR60273|c1_g1</td><td>F</td><td>-3.77</td><td>F</td><td>-4.16</td></tr><tr><td>Sox10</td><td>SRY-related HMG box 10</td><td>TR79411|c5_g2</td><td>TR64695|c1_g1</td><td>F</td><td>-3.48</td><td>NSD</td><td>(-)</td></tr><tr><td>WT1b</td><td>Wilms tumor protein 1b</td><td>TR89074|c1_g1</td><td>TR57899|c0_g1</td><td>F</td><td>-3.40</td><td>F</td><td>-2.96</td></tr><tr><td>Dax1</td><td>Dosage-sensitive sex reversal, adrenal hypoplasia critical region, on chromosome X, gene 1</td><td>TR19608|c0_g1</td><td>TR64306|c1_g1</td><td>F</td><td>-2.93</td><td>F</td><td>-0.74</td></tr><tr><td>Gata-4</td><td>Gata-binding protein 4</td><td>TR74516|c2_g1</td><td>TR63369|c1_g1</td><td>F</td><td>-2.92</td><td>NSD</td><td>(-)</td></tr><tr><td>Sf1/FTZ-F11</td><td>Steroidogenic factor-1/fushi tarazu factor-1</td><td>TR80163|c1_g2</td><td>TR56092|c1_g1</td><td>F</td><td>-2.88</td><td>F</td><td>-2.29</td></tr><tr><td>Star-like</td><td>Steroidogenic acute regulatory protein</td><td>TR52301|c0_g1</td><td>TR53314|c0_g1</td><td>F</td><td>-2.80</td><td>NSD</td><td>(-)</td></tr><tr><td>WT1</td><td>Wilms tumor protein 1a</td><td>TR89074|c0_g1</td><td>TR67279|c4_g5</td><td>F</td><td>-2.79</td><td>F</td><td>-2.41</td></tr><tr><td>Pdgfba</td><td>Platelet-derived growth factor beta a</td><td>TR86429|c1_g1</td><td>TR58625|c0_g1</td><td>F</td><td>-2.62</td><td>F</td><td>-1.97</td></tr><tr><td>Fst/Fsta</td><td>Follistatin</td><td>TR66770|c0_g3</td><td>TR48089|c0_g1</td><td>F</td><td>-2.37</td><td>F</td><td>-2.67</td></tr><tr><td>Pdgifab</td><td>Platelet-derived growth factor alpha b</td><td>TR75590|c0_g1</td><td>TR64531|c0_g1</td><td>F</td><td>-2.20</td><td>F</td><td>-1.74</td></tr><tr><td>Pdgrfb2</td><td>Platelet-derived growth factor receptor, beta 2</td><td>TR83411|c3_g1</td><td>TR55622|c2_g1</td><td>F</td><td>-2.18</td><td>F</td><td>-3.33</td></tr><tr><td>Hsd17b1</td><td>17β-Hydroxysteroid dehydrogenase type 1</td><td>TR82119|c0_g3</td><td>TR59823|c0_g1</td><td>F</td><td>-2.11</td><td>F</td><td>-6.02</td></tr><tr><td>Hdac2</td><td>Histone deacetylase 2</td><td>TR66800|c1_g1</td><td>TR56322|c1_g2</td><td>F</td><td>-2.11</td><td>F</td><td>-3.04</td></tr><tr><td>Dhh</td><td>Desert hedgehog</td><td>TR84448|c0_g1</td><td>TR59771|c0_g1</td><td>F</td><td>-2.09</td><td>NSD</td><td>(-)</td></tr><tr><td>Esr2b/Esrrb2</td><td>Estrogen receptor 2b</td><td>TR77491|c0_g1</td><td>TR52662|c0_g1</td><td>F</td><td>-1.97</td><td>NSD</td><td>(-)</td></tr><tr><td>Esrrb</td><td>Estrogen related receptor beta</td><td>TR89871|c1_g3</td><td>TR66216|c1_g1</td><td>F</td><td>-1.91</td><td>F</td><td>-4.07</td></tr><tr><td>Ara/Ar</td><td>Androgen receptor/Androgen receptor alpha</td><td>TR83707|c1_g2</td><td>TR60985|c0_g1</td><td>F</td><td>-1.89</td><td>F</td><td>-2.70</td></tr><tr><td>Hdac7</td><td>Histone deacetylase 7</td><td>TR84823|c3_g1</td><td>TR57004|c0_g1</td><td>F</td><td>-1.77</td><td>F</td><td>-1.16</td></tr><tr><td>Srd5a3</td><td>5α-reductase 3</td><td>TR83935|c0_g1</td><td>TR68483|c0_g2</td><td>F</td><td>-1.53</td><td>NSD</td><td>(-)</td></tr><tr><td>Hdac10</td><td>Histone deacetylase 10</td><td>TR88916|c0_g6</td><td>TR67490|c0_g1</td><td>F</td><td>-1.46</td><td>F</td><td>-1.42</td></tr><tr><td>Dnmt3abc1</td><td>DNA methyltransferase 3ab_contig 1</td><td>TR78244|c1_g1</td><td>TR53015|c0_g1</td><td>F</td><td>-1.29</td><td>F</td><td>-1.42</td></tr><tr><td>Cttnb1</td><td>Catenin (cadherin-associated protein), beta 1</td><td>TR64263|c0_g2</td><td>TR46279|c0_g1</td><td>F</td><td>-0.80</td><td>F</td><td>-1.22</td></tr><tr><td>Pdgifra</td><td>Platelet-derived growth factor receptor, alpha</td><td>TR88623|c1_g1</td><td>TR66258|c1_g1</td><td>F</td><td>-0.76</td><td>NSD</td><td>(-)</td></tr><tr><td>Dnmt1</td><td>DNA methyltransferase 1</td><td>TR83768|c0_g1</td><td>TR62712|c0_g1</td><td>F</td><td>-0.72</td><td>F</td><td>-0.85</td></tr><tr><td>DDX11</td><td>DEAD/H box helicase 11</td><td>TR90413|c2_g1</td><td>TR66969|c0_g1</td><td>M</td><td>0.86</td><td>M</td><td>0.80</td></tr><tr><td>Raraa</td><td>Retinoid acid receptor alpha a</td><td>TR82123|c4_g2</td><td>TR36645|c1_g1</td><td>M</td><td>0.95</td><td>M</td><td>1.61</td></tr><tr><td>Tdrd1</td><td>Tudor domain containing 1</td><td>TR91962|c2_g5</td><td>TR68477|c1_g4</td><td>M</td><td>0.96</td><td>NSD</td><td>(-)</td></tr><tr><td>Nr3c1</td><td>Glucocorticoid receptor</td><td>TR84453|c3_g1</td><td>TR63170|c1_g1</td><td>M</td><td>1.13</td><td>M</td><td>1.60</td></tr><tr><td>Tdrd7</td><td>Tudor domain containing 7</td><td>TR88488|c1_g2</td><td>TR63911|c0_g1</td><td>M</td><td>1.25</td><td>M</td><td>1.81</td></tr><tr><td>Esrra1</td><td>Estrogen-related receptor alpha</td><td>TR86271|c1_g1</td><td>TR60672|c1_g2</td><td>M</td><td>1.55</td><td>M</td><td>1.55</td></tr><tr><td>Srd5a1</td><td>Steroid-5-alpha-reductase, alpha polypeptide 1</td><td>TR91999|c2_g1</td><td>TR64324|c0_g1</td><td>M</td><td>1.58</td><td>NSD</td><td>(-)</td></tr><tr><td>Sox9</td><td>SRY-related HMG box 9</td><td>TR69398|c2_g1</td><td>TR59011|c5_g2</td><td>M</td><td>1.91</td><td>M</td><td>1.25</td></tr><tr><td>Dnmt3aa</td><td>DNA methyltransferase 3aa</td><td>TR73384|c7_g1</td><td>TR64267|c0_g1</td><td>M</td><td>2.09</td><td>M</td><td>4.11</td></tr><tr><td>Ep300a</td><td>Histone acetyltransferase—E1A binding protein 300a</td><td>TR84051|c2_g6</td><td>TR60830|c1_g1</td><td>M</td><td>2.77</td><td>M</td><td>2.49</td></tr><tr><td>Amh</td><td>Anti-Müllerian hormone or Müllerian-inhibiting substance</td><td>TR91164|c5_g4</td><td>TR66874|c1_g1</td><td>M</td><td>2.89</td><td>NSD</td><td>(-)</td></tr><tr><td>Hdac8</td><td>Histone deacetylase 8</td><td>TR86781|c2_g5</td><td>TR62362|c1_g1</td><td>M</td><td>3.04</td><td>M</td><td>2.73</td></tr><tr><td>Fstl5</td><td>Follistatin-like 5</td><td>TR77519|c0_g5</td><td>TR60004|c0_g1</td><td>M</td><td>3.22</td><td>M</td><td>3.15</td></tr><tr><td>Rarb</td><td>Retinoid acid receptor beta</td><td>TR81734|c1_g3</td><td>TR65946|c1_g2</td><td>M</td><td>3.58</td><td>M</td><td>2.99</td></tr><tr><td>Dmrt3/Dmrt3a</td><td>Double sex and mab-3 related transcription factor 3(a)</td><td>TR67762|c1_g1</td><td>(-)</td><td>M</td><td>3.80</td><td>(-)</td><td></td></tr><tr><td>Cyp11c1/Cyp11b</td><td>Cytochrome P450, family 11, subfamily C, polypeptide 1/11beta-hydroxylase</td><td>TR72664|c0_g1</td><td>TR47689|c0_g1</td><td>M</td><td>5.34</td><td>NSD</td><td>(-)</td></tr><tr><td>Dmrt1</td><td>Doublesex- and mab-3-related transcription factor 1</td><td>TR70670|c0_g1</td><td>TR67635|c0_g1</td><td>M</td><td>6.74</td><td>M</td><td>6.21</td></tr></tbody><tfoot><tr><td colspan="8">Continued</td></tr></tfoot></table>
<table><caption>**Table 1.** Expression pattern of gonadal genes implicated in sex determination/differentiation of common pandora and red porgy. **M**: male-biased, **F**: female-biased, **NSD**: not significantly different.</caption><thead><tr><th rowspan="2">Gene name</th><th rowspan="2">Gene description</th><th colspan="2">Gene ID</th><th colspan="4">Expression</th></tr><tr><th>Common pandora</th><th>Red porgy</th><th colspan="2">Common pandora</th><th colspan="2">Red porgy</th></tr><tr><th></th><th></th><th></th><th></th><th>sex</th><th>log2FC</th><th>sex</th><th>log2FC</th></tr></thead><tbody><tr><td>Hsd11b3</td><td>11β-Hydroxysteroid dehydrogenase type 3</td><td>TR71766|c4_g1</td><td>TR60165|c5_g1</td><td>NSD</td><td>(-)</td><td>F</td><td>-6.62</td></tr><tr><td>Fgf20b</td><td>Fibroblast growth factor 20-like/b</td><td>TR73946|c0_g1</td><td>TR56801|c1_g1</td><td>NSD</td><td>(-)</td><td>F</td><td>-3.46</td></tr><tr><td>Gsdf</td><td>Gonadal soma derived factor</td><td>TR72782|c0_g1</td><td>TR61049|c0_g1</td><td>NSD</td><td>(-)</td><td>F</td><td>-3.23</td></tr><tr><td>Rarab</td><td>Retinoid acid receptor alpha b</td><td>TR84655|c1_g1</td><td>TR63371|c2_g1</td><td>NSD</td><td>(-)</td><td>F</td><td>-2.38</td></tr><tr><td>Cyp11a2</td><td>Cytochrome P450, family 11, subfamily A, polypeptide 2/11alpha-hydroxylase</td><td>TR88875|c1_g5</td><td>TR61267|c0_g1</td><td>NSD</td><td>(-)</td><td>F</td><td>-2.06</td></tr><tr><td>Hdac11</td><td>Histone deacetylase 11</td><td>TR84940|c1_g2</td><td>TR66619|c2_g2</td><td>NSD</td><td>(-)</td><td>F</td><td>-1.44</td></tr><tr><td>KAT2b</td><td>Histone acetyltransferase—K(lysine) acetyltransferase 2b</td><td>TR81329|c0_g2</td><td>TR59334|c0_g1</td><td>NSD</td><td>(-)</td><td>F</td><td>-1.43</td></tr><tr><td>Nr3c2</td><td>Mineralocorticoid receptor</td><td>TR84031|c1_g1</td><td>TR58538|c2_g1</td><td>NSD</td><td>(-)</td><td>F</td><td>-1.06</td></tr><tr><td>Amhr2</td><td>Anti-Müllerian hormone receptor 2</td><td>TR90446|c1_g4</td><td>TR66380|c1_g1</td><td>NSD</td><td>(-)</td><td>F</td><td>-0.84</td></tr><tr><td>Fgf9</td><td>Fibroblast growth factor 9</td><td>TR91211|c1_g5</td><td>TR66253|c2_g4</td><td>NSD</td><td>(-)</td><td>M</td><td>1.86</td></tr><tr><td>Esr1/Er</td><td>Estrogen receptor</td><td>TR84515|c0_g1</td><td>TR68522|c4_g2</td><td>NSD</td><td>(-)</td><td>M</td><td>1.92</td></tr><tr><td>Pdgfa1</td><td>Platelet-derived growth factor alpha a</td><td>TR87951|c1_g6</td><td>TR66514|c0_g1</td><td>NSD</td><td>(-)</td><td>M</td><td>2.65</td></tr><tr><td>Srd5a2</td><td>5α-reductase 2</td><td>(-)</td><td>TR59720|c7_g2</td><td>(-)</td><td>(-)</td><td>M</td><td>3.17</td></tr></tbody></table>

sex-biased expression in both brains and gonads were over-expressed in female brains. Zona pellucida proteins are glycoproteins of the fish chorion that are encoded by multiple gene families. In humans and mice, but also in some teleosts⁵⁹⁻⁶¹, they are expressed in the ovary. In addition, several factors of the *Sox* gene family, as well as several members of the Forkhead box transcription factor (*Fox*) family were detected, which might play important roles in the brain of these protogynous fish. The comparative analysis revealed commonly over-expressed genes in the male brains of both species, pinpointing candidates that might make a significant contribution to the expression of the male character of the brain in protogynous sparisids. Among these, the neurobeachin (*Nbea*) is required in synaptic function, regulating neurotransmitter receptor trafficking to synapses⁶². A male over-expression shared between species also displayed the *PRKCE* (Protein kinase C epsilon type) and *EphA4* (*Ephrin type-A receptor* 4). The first is implicated in platelet activation and positive regulation of the mitogen-activated protein kinase (MAPK) cascade, while the second is a receptor tyrosine kinase implicated in the ephrin receptor signaling pathway and positive regulation of Jun kinase activity. The role of such kinases in mammalian sex determination has been reported⁶³, while mice studies underline the role of the phosphorylated form of MAPK (pMAPK) as a marker of neuronal activity in a reproductive context⁶⁴.

In the common pandora and red porgy brains we also found the early growth response factor 1 (*Egr1*), the early growth response factor 4 (*Egr4*), as well as *Fos* (Proto-oncogene c-Fos) and all had a male-biased expression profile. The *Fos* gene, which encodes a nuclear phosphoprotein, reflects immediate neural activity, whereas *Egr1* is believed to regulate later-acting genes, such as gonadotropin-releasing hormone 1 (GnRH1), a primary signal essential for reproduction⁶⁵. The *Egr1* and *c-Fos*, both characterized as immediate-early genes, are related to changes in brain activity and genomic responses to social cues, behavior and mating information⁶⁶.

Overall, in the brain of both protogynous species there were some sex-specific expression differences, although it is not known whether these differences mediate sexually dimorphic events. The more “homogeneous” transcriptome profile of the brain, initially perceived by the lack of clustering in male and female brain, became more apparent in GO term enrichment analysis, where GO terms related to female development and differentiation were over-represented in male brains and vice-versa. Sex differences in the brain, even in the human brain, are still poorly defined. As such, the issue of “brain sex differentiation” remains a controversial subject of research, especially in fish. Furthermore, the well-established phenomenon in mammals of fetal or neonatal brain exposure to sex hormones producing irreversible differences of brain structure and function, appears not to be the case in fish. Instead, fish brains are not irreversibly sexualized, conserving properties of the embryonic mammalian brain, such as high neurogenic activity, throughout life⁶⁷. The notion that fish exhibit much more plasticity with regard to brain sex differentiation has been strengthened by the findings of the present study.

**Major sex differences in gonadal gene expression patterns.** The gonads is the tissue that differs the most between sexes⁶⁸. Despite the lability of gonadal tissue in sex-changers and the fact that testes and ovaries share a common precursor, they clearly represent two distinct types of tissue upon the establishment of sex. The functional divergence of the two gonadal types is reflected in their transcriptomic profiles, in terms of the number of genes differentially expressed, as well as the expression magnitude (i.e. fold-change differences). The observation of almost double the number of up-regulated genes in males compared to females indicates a male-biased expression tendency. More and more transcriptomic studies dealing with sex-specific expression differences in teleosts have pointed out male-biased over-expression tendencies¹⁶,³⁵,³⁸,³⁹,⁵⁷,⁶⁹,⁷⁰. In zebrafish, a study exclusively addressed the “male sex drive” hypothesis, originally proposed by Singh and Kulathinal⁷¹, suggesting that male-biased evolutionary pressures might have resulted in a “transcriptome masculinization” of the species. The present study, with the two protogynous hermaphrodites, along with the previous one on rudimentary sharpsnout seabream³⁵, agree with a skewed ratio of male up-regulated versus female up-regulated transcripts in gonads, in
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Comparative view of known sex-related genes, commonly (intersection) and uniquely expressed in common pandora (left) and red porgy (right). Above and below the central line are the male-specific (Blue) and the female-specific (Red) genes. The Y-axis serves as an indication of expression magnitude (in terms of Log2Fold-Change differences). The commonly expressed genes (intersection) are placed according to their Log2Fold-Change values in common pandora and had an analogous expression magnitude in red porgy. The genes uniquely expressed in each species are placed according to their expression values in that particular species.</figcaption></figure>

favor of the first. In cases such as the protogynous hermaphrodites studied here, with no signs of heterogametic chromosomes, such a tendency may explain a great deal about the formation and maintenance of the male and female phenotypes.

**Expression profile of important molecular players in protogynous hermaphrodite gonads.**

Studies on vertebrates and other fish species have brought to light several candidate genes with proven or suspected relationship in sex determination and differentiation. As there is scanty literature-based information regarding the sex-related genetic cascade in these protogynous species, we compiled a list of well-studied sex-biased genes and searched for their expression patterns in the gonadal transcriptomes of both sparids. Despite the plethora of different sexual systems in fish, there is some evidence of conservation of the genetic cascade, with its "basic" components (molecular players) found repeatedly, not only among fishes but also generally in vertebrates<sup>72,73</sup>. However, the specific role and place for each component of this at least partially conserved "crew" is not prominent and might change in different species and/or sexual systems<sup>57,74,75</sup>.

Based on the literature, we profiled the expression of known candidate (Fig. 6) molecular-players/genes establishing the common female (*Cyp19a1*, *Sox3*, *Foxl2*, *Figla*, *Wnt4a*, *Gdf9*, *Fst*, *Lhb*, *Ara*, *Arb*, etc.) and male identity (*Dmrt1*, *Sox9*, *Dnmt3aa*, *Rarb*, *Raraa*, *Hdac8*, etc.) of the gonads of these protogynous sparids. Additionally, we focused on those contributing in species-specific manner either to female (*Gata-4*, *Sox10*, *Amhr2*, *Gsdf*, *Srd5a3* etc.) or to male (*Amh*, *Dmrt3a*, *Srd5a1*, *Srd5a2*, etc.) characters. Starting with the formation and maintenance of the female phenotype, some genes over-expressed, showed high Log₂ fold changes in the ovaries of both sparids. The *Figla* gene that appears to have a decisive role in the female phenotype of both species is such an example. This transcription factor up-regulates ZP genes, which code for products used to build the vitelline envelope surrounding eggs in teleosts (zona matrix)⁷⁶⁻⁷⁸. Moreover, it has been suggested that the *Figla* factor displays an anti-parallel direction of expression to *Dmrt1* during the transition from male to female in a protandrous⁷⁹ and vice versa in a protogynous teleosts⁸⁰. The same expression pattern was showed for *Foxl2*, a prominent marker of ovarian development across vertebrates⁸¹. The *Foxl2* has been reported as a female-biased gene in several gonochoristic⁸²⁻⁸⁵ and, recently, in some hermaphroditic species³⁸,³⁹.

The gene of ovarian aromatase is probably the one with the most protagonistic role in gonadal fate. Together with estrogens it has been hypothesized to control not only the ovarian, but also the testicular differentiation. When it is over-expressed it appears to trigger and maintain ovarian differentiation and once its expression reduces it appears to allow for testicular differentiation to take place⁸⁶. Here, a strong female-biased expression of gonadal aromatase gene (*Cyp19a1a*) was observed in both species. In fish, *Cyp19a1a* has an essential role in ovarian differentiation. Its expression was not only reported in gonochorists⁷⁵,⁸⁷⁻⁸⁹, but also in other hermaphrodite fish³⁹,⁹⁰⁻⁹⁴ which repeatedly highlights its involvement with the gonadal sex differentiation and sex change. It has been also suggested that the expression of this conserved and well-characterized protein is controlled not only by a "classical" transcriptional activator such as *Foxl2*, but also by epigenetic modifications⁹⁵. Gonadal aromatase gene and its regulation through methylation modifications has also been involved in thermo-sensitive
sex differentiation systems in fish⁹²,⁹⁶. Thus, changes in methylation levels of the *Cyp19a1a* promoter might be an additional control factor in sex-changing species⁹⁵. The precise role of aromatase in the sex determination/differentiation systems of these species awaits further and deeper studies. The *Cyp19a1a* promoter contains binding sites for *Wilms tumor 1 protein* (*Wt1*)⁸⁶. In mammals, there is only one *Wt1* gene critical for development of the urogenital system, whilst in teleosts two *Wt1* (*Wt1a*, *Wt1b*) genes have been identified and both are over-expressed in ovaries⁹⁷.

The Luteinizing hormone (Lh) and its receptor (Lhr) showed a strong female over-expression in both species examined here. According to studies in medaka, both genes are involved in fish ovulation⁹⁸. A common female-biased expression also displayed the fushi tarazu factor-1 (*Ftz-F1*) or steroidogenic factor-1 (*Sf-1*). Interestingly, the opposite has been described in two other hermaphrodites³⁸,³⁹ and in some well-known gonochorists⁹⁹,¹⁰⁰, where *Sf-1* gene has been reported as a male-biased gene. It has been hypothesized that *Sf-1/Ftz-F1* along with some additional partners, such as *Foxl2* and *cAMP*, co-regulate the transcription of aromatase, as the later contains binding sites for *Sf-1* genes⁸⁶. Two members of the Wnt/b-catenin pathway, *Wnt4* (wingless-type MMTV integration site family, member 4) and *Ctnnb1* (Catenin, beta 1) and *Fst* (follistatin) well recognized markers of ovarian differentiation linked to this pathway, all seem to contribute to the female identity of the gonads in both spарids studied here.

Three factors of the retinoic acid (RA) signaling pathway were found in the transcriptome of both species and of those three factors *Cyp26a1* over-expressed in their ovaries of both spарids. The role of this pathway in sex determination/differentiation and particularly in ovarian differentiation¹⁰¹⁻¹⁰³ is further supported by respective findings among ovary-enriched GO terms. Recently, a RA homeostatic mechanism stimulated by retinoic acid gene 8 (*Stra8*) via an independent pathway has been shown to implicate that *Cyp26a1* and *Aldh1a2* mediate the meiotic entry of germ cells in Nile tilapia¹⁰². While *Stra8* was absent from the transcriptomes of the studied spарids, *Stra6* was present and over-expressed in the ovaries of both species. A search for genes encoding RA receptors (*Raraa*, *Rarab* and *Rarb*) revealed that all three exhibited sex-biased expression profiles. The first two were commonly male-biased, while *Rarb* was over-expressed in red porgy ovaries.

Four members of the TGF-b superfamily (*Gdf9*, *Gsdf*, *Amh*, *Amhr2*) were found in both species' transcriptomes and displayed sex-biased expression. In particular, growth differentiation factor-9 (*Gdf9*), which is required for ovarian folliculogenesis¹⁰⁴ was up-regulated in the ovaries of both species. Our findings are in line with those of another protogynous fish, the bluehead wrasse³⁸. The gonadal soma derived factor (*Gsdf*), a gene probably unique to teleosts and recently identified as the male-specific master sex-determining gene in *Oryzias luzonensis*¹⁰⁵ was uniquely over-expressed in the ovaries of red porgy, possibly indicating a different role for the *Gsdf* gene in this protogynous species. A well-known member of the TGF-b family, the anti-Mullerian hormone (*Amh*) seems to contribute to the common pandora male phenotype, while it showed no sex-biased differences in red porgy. Interestingly, its receptor, *Amhr2*, is female-biased in the later species. The *Amh* and *Amhr2* are known master SD genes in the Patagonian pejerrey *Odontesthes hatcheri*¹¹ and in the tiger pufferfish *Takifugu rubripes*¹², respectively. The *Amh* role in teleost sexual differentiation, which lack Müllerian ducts, is not clear, though it has been proposed that proliferation and differentiation of spermatogonia are the possible common functions of this gene among vertebrates¹⁰⁵,¹⁰⁶. The same gene has also been found among the male-specific genes in some gonochoristic⁵⁷,¹⁰⁷ and hermaphroditic³⁵,³⁹ fish and it is involved in the sex change of the protandrous black porgy³⁰.

Two major components of the *Sox* genes¹⁰⁸ (*Sox3* and *Sox9*) were found among the commonly expressed gonadal genes between the two protogynous species. The present study demonstrated that *Sox3* has a high female-specific expression in both common pandora and red porgy. In mammals, *Sox3* is implicated in the induction of testis formation¹⁰⁹ and spermatogenesis¹¹⁰. However, its role in fish sex determination and differentiation is still ambiguous as in some studies it is linked to the male phenotype¹⁴,¹¹¹, while in others to the female³⁸,⁵⁷,¹¹². The other *Sox* gene, *Sox9*, is a male-biased gene in both species. The direction of expression of *Sox9* agrees with the male-specific role of this gene in mammals¹¹³. Nevertheless, *Sox9* has been reported to have a male-biased expression in rudimentary sharpsnout seabream³⁵, the protogynous bluehead wrasse³⁸ and the protandrous clownfish³⁹. On the contrary, in the medaka *Sox9* is expressed in ovaries¹¹⁴. Interestingly, in zebrafish *Sox9* has two homologues, one expressed in oocytes and the other in testicular Sertoli cells¹¹⁵. Taken together, the results suggest that *Sox3* and *Sox9* are important in the common pandora and red porgy female and male phenotype formation respectively and suggest both genes as handy molecular players that could contribute to either sex side, in teleosts.

Despite major differences in the genetic control of sexual development in different organisms, *Dmrt* genes are actively involved in sex-specific differentiation in vertebrates, invertebrates and all teleosts that have been studied so far¹¹⁶⁻¹¹⁸. A family member, *Dmrt2/Dmrt2a* which is linked to neurogenesis¹¹⁹ and left-right patterning¹²⁰ during embryonic development, displayed a female-biased expression shared between species. This agrees with the patterns observed in the Atlantic cod¹²¹ and four cichlids species⁵⁷. On the other hand, the most studied member of this gene family, *Dmrt1*, was highly up-regulated in the testes of both protogynous species studied. Hence, the results of this study indicate that *Dmrt1* together with *Sox9* could play a fundamental role in expressing and preserving the male character in gonad. The role of *Dmrt1* in sex determination and spermatogenesis in vertebrates and in testicular development and maintenance in teleosts is well established¹¹⁷. In line with the current study results, characteristic male-biased expression of *Dmrt1* has been reported in numerous gonochoristic⁹⁹,¹²²⁻¹²⁴ and other hermaphrodite teleosts³³,³⁵,³⁸,³⁹,⁹¹. *Dmrt3/Dmrt3a*, another male-specific gene in general¹²⁵, appears to contribute to male phenotype of common pandora, as it was strongly expressed in their testes, but it was absent from the red porgy transcriptome.

Another interesting finding from our study was the expression profiles of androgen and estrogen receptor genes (*Ara*, *Arb*, *Er/Esr1* and *Esr2b*), as well as estrogen-related receptors (*Esrra* and *Esrrb*) in the gonadal transcriptomes of both species. Surprisingly, both androgen receptors were up-regulated in the ovaries of both species
and *Esrra* displayed a common male up-regulation. Although it is still unclear how the endocrine mechanism works to control sex change, steroid hormones are believed to regulate the process¹²⁶,¹²⁷.

Finally, genes of epigenetic maintainers, such as the DNA-modifying enzymes DNA methyltransferases (Dnmt3aa, Dnmt1, Dnmt3ab) and histone tail-modifying and de-modifying enzymes, such as histone acetyl transferases (Ep300a, Kat2b) or histone deacetylases (Hdac8, Hdac2, Hdac10, Hdac11) were found in transcriptome of both protogynous sparids and showed a sex-biased expression. Further, enrichment analysis revealed genes involved in epigenetic related mechanisms, supporting their implication in sex determination/differentiation and apparently sex maintenance. The relationship of the epigenetic mechanisms to sex determination and differentiation in vertebrates has been clearly stated and these mechanisms are considered to be the link between environmentally controlled gene expression and sex reversal⁹⁵.

## Conclusions

In this study, we sequenced the total mRNA from the gonads and brains of both sexes in two protogynous hermaphrodite sparids, the common pandora and the red porgy. This approach permitted us to obtain a global view of sex-specific expression patterns, as well as it allowed us to compare sex-biased gene expression patterns within teleosts having the same reproductive mode. Focusing on the pathways and genes implicated in sex determination/differentiation, we tried to unveil the molecular pathway by which a protogynous hermaphrodite fish develops a masculine or a feminine character. We found minor sex-related expression differences indicating a more homogeneous and sexually plastic brain, whereas there was a plethora of sex biased gene expression in the gonads. This indicates functional divergence of the transformed tissue either to male or female character, even though both types of gonadal tissue originated from the same precursor. We observed the implicated pathways behind sex-biased expression and saw the recruitment of known sex-related genes either to male or female type of gonads in these fish. Taken together, most of the studied genes form part of the cascade of sex determination, differentiation and reproduction across teleosts. In this study, we focused on genes that are active when sex is maintained (sex-maintainers). Comparing related species with the same reproductive style, we saw different combinations of genes with conserved sex-linked roles and some handy molecular players, in a partially conserved network formulating the male and female phenotype. Finally, this analysis lays the ground for understanding the complex process of sex differentiation in these protogynous sparids and makes future comparative analyses of fish with alternative reproductive modes feasible. Future work in gene expression data during sex-reversal will uncover which genes are triggered to block female path and start to transform the gonadal tissue in the opposite direction.

## Methods

**Sample collection.** Animal care was carried out according to the “Guidelines for the treatment of animals in behavioural research and teaching”¹²⁸. The fish used in the study were maintained in registered and approved facilities for the maintenance and carrying of animal experiments and were reared and sampled according to the guidelines of the Directive 2010/63/EU for the protection of animals used for experimental and other scientific purposes (Official Journal L276/33)¹²⁹ and the guidelines of Metcalfe and Craig (2011)¹³⁰. In addition, the experimental sampling protocol was approved by the IMBBC's aquaculture department committee and all methods were carried out in accordance with relevant guidelines and regulations approved by the Ministry of Rural Development and Food and the Regional Directorate of Veterinary Medicine for certified experimental installations (EL 91-BIO-04) and experimental animal breeding (AQUALABS, EL 91-BIO-03). Members of the laboratory include certified technicians by the Federation for Laboratory Animal Science Associations (FELASA).

The fish used for the experiments were either wild fish caught alive, or broodstock fish of wild origin in the Institute of Marine Biology, Biotechnology and Aquaculture (HCMR) (Supplementary Excel Tables SET9). In both cases, mature fish were euthanized using the commercially available slaughter method for marine fish (i.e. immersion in ice-slurry) and the brain and gonads were dissected immediately. Only one of the two gonadal lobes was kept and stored until further processing. The fish were examined macroscopically for sexual maturation, based on the presence of releasable sperm for males or the presence of vitellogenic oocytes for females and/or visual inspection of the gonad morphology. For red porgy, nine brain and nine gonad tissues were sampled from five female and four male individuals, while for common pandora ten brain and ten gonad tissues were sampled from five females and five males. In addition, for each species a pool of several hundred larvae at the “hatching” stage (1 days post hatch) were sampled to increase the possibility of capturing the majority of the species transcripts. All tissues and larvae were collected in a sterile and RNase-free way and stored immediately in RNALater (Applied Biosystems, Foster City, CA, USA). Tissues for RNA extraction were stored at 4°C overnight according to the manufacturer's recommendations and then transferred to -80°C until further processing.

**RNA extraction and sequencing.** To reduce biases regarding different expression in different parts of the same organ, we homogenized the whole brain, one lobe of testis and one lobe of ovaries. All tissue samples, including larvae, were ground with liquid nitrogen using pestle and mortar, homogenized in TRIzol® reagent (Invitrogen) and total RNA was extracted from the TRIzol® homogenate according to the manufacturer's instructions. The quantity of the isolated RNA was measured spectrophotometrically with NanoDrop® ND-1000 (Thermo Scientific), while its quality was tested on an agarose gel (electrophoresis in 1.5% w/v) and further on an Agilent Technologies 2100 Bioanalyzer (Agilent Technologies). The majority of the samples had an RNA Integrity Number (RIN) value higher than 8. This was not the case for one testes from common pandora and it was excluded from the further analysis. In contrast to their high quality agarose image, all ovaries showed lower RIN values. The low RIN values obtained from ovary samples could be due to the abundance of low molecular weight 5S RNA expression. In a previous study of sharpsnout sea bream³⁵ where the RNA extractions were
carried out with a commercial kit and the female gonads were in an early maturing stage, the 5S rRNA was in such abundance that it masked the 18S and 28S rRNA peaks thus hampering the correct estimation and calculation of RNA integrity numbers (RIN). This low molecular weight band of 120 bp, belonging to the smallest rRNA molecule, was also detected in the ovaries and intersex gonads of thicklip grey mullet *Chelon labrosus*, where the role of 5S rRNA in fish sex differentiation and its potentiality as a sex marker in fish gonads has been thoroughly described¹³¹. This type of RNA profile was also observed in female gonads of bluehead wrasse³⁸ and very recently a study described the use of ribosome RNA profile not only as a sex identifier but also as an indicator of ovarian developmental stage¹³². Finally together with two larval samples, 39 samples were used for mRNA paired-end library construction with the Illumina TruSeq™ RNA Sample Preparation Kits v2 following the manufacturer's protocol (Poly-A mRNA isolation with oligo-dT beads, mRNA fragmentation, followed by transcription into first-strand cDNA using reverse transcriptase and random hexamer primers) and sequenced as 100 bp paired reads in three lanes (1,5 lanes per species) of a HiSeq. 2500™ following the protocols of Illumina Inc. (San Diego, CA).

**Read pre-processing.** All computations were performed at the high-performance computing bioinformatics platform of HCMR. Read quality was assessed with FastQC¹³³ and subjected to quality control following a pipeline using Scythe¹³⁴, Sickle¹³⁵, Trimmomatic¹³⁶ and Prinseq (prinseq-lite-0.20.3)¹³⁷. We first used Scythe - A Bayesian adapter trimmer (version 0.994 BETA)¹³⁴, which uses a Naive Bayesian approach to classify contaminant substrings in sequence reads; it considers quality information which can make it robust in picking out 3'-end adapters (which often include poor-quality bases) and we therefore ran it before any quality-based trimming (with the prior contamination rate set in 0.1 '-p 0.1'). Low-quality (parameters 'pe -g -t sanger -q 20 -l 45') nucleotide reads trimming was implemented with Sickle¹³⁵, a sliding-window, adaptive, quality-based trimming tool for FastQ files. Using Trimmomatic¹³⁶, some extra filtering steps and also both 5' and 3' adaptor removal were performed (parameters 'PE -phred33 ILLUMINACLIP:adapter_file.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:25 MINLEN:45 CROP:99'). Finally, low-complexity sequences (threshold entropy value of 30) and poly A/T 5' tail (minimum of 5 A/T) trimming was performed using PrinSeq¹³⁷. Lastly, the reads of each pair of all the libraries were combined using custom Perl scripts.

**Transcriptome assembly construction.** The *de novo* assemblies for red porgy and common pandora were performed using Trinity package¹³⁸ (version 2.0.6) with default parameters (default kmer 25, min length 200 nucleotides). To increase transcriptome coverage, *de novo* assemblies were carried out on a pool of the pre-processed reads from all samples. To assess the assembled transcripts and exclude the spurious ones, we pooled all filtered reads and mapped them to the selected assemblies (Pagrus_Trinity_assembly.fasta and Pagellus_Trinity_assembly.fasta) with Bowtie²⁵¹ within RSEM (version 1.2.19)¹³⁹, using the script available in trinity utilities *align_and_estimate_abundance.pl*. Putative transcripts with less than 1% of a locus reads mapped to that particular isoform (IsoPct < 1) were excluded (as proposed in)¹⁴⁰. The same was done for those with Fragments Per Kilobase of transcript per Million mapped read (FPKM) values less than 0.3. The choice of FPKM threshold was based on BLASTn similarity searches (*e*-value threshold 10⁻¹⁰) against European sea bass cDNA sequences (http://seabass.mpipz.mpg.de). We first contacted a blast search of the unfiltered assemblies against the European sea bass cDNA sequences and then repeated the blast and tested three differently filtered assemblies (a) FPKM_cutoff = 0.3, isopct_cutoff = 1.0, (b) FPKM_cutoff = 0.5, isopct_cutoff = 1.0 and (c) FPKM_cutoff = 1, isopct_cutoff = 1.0. To assess the assembled transcripts quality, we followed three methods: (a) pooling all the reads and mapping them to the assembly using Bowtie²⁵¹ and Samtools¹⁴¹ to estimate the percentage of mapped reads, (b) examining the similarity between the Trinity assembly and swissprot database using blastx search and (c) quantifying the completeness of the transcriptome assembly and determining how many of the 3,023 genes for vertebrates (OrthoDBv8 was used from BUSCO) were present in our two species transcriptomes, using BUSCO (Benchmarking Universal Single-Copy Orthologs) software¹⁴² that searches for highly conserved, near-universal, single copy orthologs. All BLAST¹⁴³ searches were performed in parallel using NOBLAST program¹⁴⁴.

**Transcriptomes functional annotation.** For the annotation of the assembled transcripts, we conducted a BLASTx similarity search against UniProtKB/SwissProt database⁵² (*e*-value threshold 10⁻⁵) keeping the top 20 hits. BLASTx was done in parallel using NOBLAST¹⁴⁴. For gene ontology (GO) mapping, Blast2GO⁵⁴ was used locally to recover all the GO terms associated to the hits obtained by the blast search. After the mapping step, results were subjected to GO annotation, a process of selecting GO terms and assigning them to the query sequences for describing biological processes, molecular functions and cellular components. Gene ontology graphs were generated in "level 2" for GO plot visualization and pie charts for the whole transcriptome profiling. Further, brains and gonads of each species were analyzed separately, using in-house R scripts. The Blast2GO output file was input into the BGI WEGO program⁵⁵ and GO annotations were plotted. The sequences were additionally annotated using InterPro¹⁴⁵; a local InterProScan (version 5, InterPro release 48.0)⁵³ was run in parallel (splitting the query set) on the longest open reading frame (ORF) of the contigs using a custom Perl script based on the EMBOSS program "getorf"¹⁴⁶. GO terms corresponding to these InterPro domains, were merged with the already existing GO terms derived from the blastx against UniProt. The GO terms were modulated using the annotation augmentation tool ANNEX¹⁴⁷ followed by GOSlim¹⁴⁸. Enzyme classification (EC) codes were obtained through the direct mapping of GO terms to the corresponding enzyme codes. Sequences having EC numbers were further characterized by Kyoto Encyclopedia of Genes and Genomes (KEGG)¹⁴⁹ metabolic pathway annotations using custom Perl/R scripts.
**Differential expression analysis.** Differential expression analysis was performed at the "gene" (component/unigene) level by pairwise comparisons between male and female for brain and gonad samples. The *de novo* transcriptome assembly of each species served as a reference for read mapping. For the quantification, the cleaned reads of gonads and brains (not larvae) of each sample of a species were aligned to the species' transcriptome assembly with Bowtie2<sup>51</sup> and abundance was estimated with RSEM, as implemented in the trinity script *align_and_estimate_abundance.pl*. The estimated expected counts for each sample, when the sum of counts in each row >15, were extracted<sup>150</sup>. Loci with less than fifteen reads mapped on them, were excluded prior to differential expression analysis to improve the statistical power. Count data were imported in R (version 3.2.0) and the analysis of differential expression conducted in DESeq. 2 (version: 1.8.1)<sup>151</sup>, an R bioconductor package<sup>152</sup>. Samples of both species were grouped according to sex and expression was compared separately for brains and gonads, following the developers' manual (false discovery rate, FDR, threshold of 0.05). Principal component analysis (PCA)<sup>153</sup> and the heatmap.2 function in the gplots package<sup>154</sup> were used to visualize global similarities and differences among either the brain or the gonadal samples. The global expression analysis conducted with sex and tissue being the (*design/variable/intgroup*) factors influencing the counts in a multifactor analysis. Each species global (brains and gonads) PCA plot was produced including all the genes in the filtered (sum of counts in each row >15) count matrix. Genes with an adjusted *p* value less than 0.05 were considered as significantly differentially expressed.

**Functional GO and pathway enrichment analysis.** Functional-enrichment analysis was performed to identify GO terms and metabolic pathways significantly enriched in differentially expressed genes in brains and gonads, in the two sparids. In each species, male and female up-regulated gonad and brain gene sets were tested for enriched GO terms, in a comparison of female versus male (female/test-set vs male/reference-set). Particularly for gonads, we restricted the gene-set to only those with $|log₂FC| > 2$. We also conducted enrichment analysis for common (between species) male-biased and female-biased gene set using both 'red porgy' and 'common pandora' total datasets as reference set of the annotated transcripts. All enrichment analyses were implemented in BLAST2GO with default settings enabling for a two-tailed test and reducing the output table to the most specific terms with FDR = 0.05.

**Targeted search for sex-related candidate genes.** To identify transcripts related with the GO terms "reproduction" (GO:0000003), "reproductive process" (GO:0022414), "sex determination" (GO:0007530) and "sex differentiation" (GO:0007548), both annotated transcriptomes were searched for the presence of these GO terms and the retrieved genes were checked for whether they exhibit a male or female-biased expression in brains and gonads. Genes known to be involved in sex determination/differentiation in vertebrates and other fish taxa were downloaded for tilapia, medaka and spotted gar *Lepisosteus oculatus* from Ensembl (release 84)<sup>155</sup>, resulting in a set of 114 candidate genes. To find out if orthologs of these genes were present in our species transcriptomes, we used a reciprocal BLASTn best-hit approach<sup>156</sup>. The potentially orthologous genes were then tested based on the differences in their expression patterns between testes and ovaries of both species by using the output files of DeSeq. 2. Those that were differentially expressed in both species or were species-specific were identified.

**Availability of supporting data.** The datasets generated and/or analysed during the current study are available from the corresponding author upon request. Raw sequence reads can be found in the SRA database under BioProject ID: PRJNA395994.

**References**

1. Mittwoch, U. Sex determination in mythology and history. *Arq. Bras. Endocrinol. Metabol.* **49**, 7-13 (2005).

2. Heule, C., Salzburger, W. & Böhne, A. Genetics of Sexual Development: An Evolutionary Playground for Fish. *Genetics* **196**, 579-591 (2014).

3. Desjardins, J. K. & Fernald, R. D. Fish sex: why so diverse? *Curr. Opin. Neurobiol.* **19**, 648-653 (2009).

4. Tobias, U. & Heikki Helanterä. From the Origin of Sex-Determining Factors to the Evolution of Sex-Determining Systems. *Q. Rev. Biol.* **86**, 163-180 (2011).

5. Yoshimoto, S. et al. A W-linked DM-domain gene, DM-W, participates in primary ovary development in Xenopus laevis. *Proc. Natl. Acad. Sci. USA* **105**, 2469-2474 (2008).

6. Sinclair, A. H. et al. A gene from the human sex-determining region encodes a protein with homology to a conserved DNA-binding motif. *Nature* **346**, 240-244 (1990).

7. Koopman, P., Gubbay, J., Vivian, N., Goodfellow, P. & Lovell-Badge, R. Male development of chromosomally female mice transgenic for Sry. *Nature* **351**, 117-121 (1991).

8. Smith, C. A. et al. The avian Z-linked gene DMRT1 is required for male sex determination in the chicken. *Nature* **461**, 267-271 (2009).

9. Matsuda, M. et al. DMY is a Y-specific DM-domain gene required for male development in the medaka fish. *Nature* **417** (2002).

10. Myosho, T. et al. Tracing the Emergence of a Novel Sex-Determining Gene in Medaka, Oryzias luzonensis. *Genetics* **191**, 163-170 (2012).

11. Hattori, R. S. et al. A Y-linked anti-Müllerian hormone duplication takes over a critical role in sex determination. *Proc. Natl. Acad. Sci. USA* **109**, 2955-2959 (2012).

12. Kamiya, T. et al. A Trans-Species Missense SNP in Amhr2 Is Associated with Sex Determination in the Tiger Pufferfish, Takifugu rubripes (Fugu). *PLoS Genet* **8**, e1002798 (2012).

13. Yano, A. et al. An Immune-Related Gene Evolved into the Master Sex-Determining Gene in Rainbow Trout, Oncorhynchus mykiss. *Curr. Biol.* **22**, 1423-1428 (2012).

14. Takehana, Y. et al. Co-option of Sox3 as the male-determining factor on the Y chromosome in the fish Oryzias dancena. *Nat. Commun.* **5**, 4157 (2014).

15. Chen, S. et al. Whole-genome sequence of a flatfish provides insights into ZW sex chromosome evolution and adaptation to a benthic lifestyle. *Nat Genet* **46**, 253-260 (2014).

16. Tao, W. et al. Characterization of gonadal transcriptomes from Nile tilapia (Oreochromis niloticus) reveals differentially expressed genes. *PLoS One* **8**, e63604 (2013).
17. Nicol, B. & Guiguen, Y. Expression Profiling of Wnt Signaling Genes during Gonadal Differentiation and Gametogenesis in Rainbow Trout. *Sex. Dev.* **5**, 318–329 (2011).

18. Kobayashi, T., Kajiura-Kobayashi, H., Guan, G. & Nagahama, Y. Sexual dimorphic expression of DMRT1 and Sox9a during gonadal differentiation and hormone-induced sex reversal in the teleost fish Nile tilapia (Oreochromis niloticus). *Dev. Dyn.* **237**, 297–306 (2008).

19. Kovács, B., Egedi, S., Bártfai, R. & Orbán, L. Male-specific DNA markers from African catfish (Clarias gariepinus). *Genetica* **110**, 267–276 (2000).

20. Ezaz, M. T. et al. Isolation and Physical Mapping of Sex-Linked AFLP Markers in Nile Tilapia (Oreochromis niloticus L.). *Mar. Biotechnol.* **6**, 435–445 (2004).

21. Felip, A., Young, W. P., Wheeler, P. A. & Thorgaard, G. H. An AFLP-based approach for the identification of sex-linked markers in rainbow trout (Oncorhynchus mykiss). *Aquaculture* **247**, 35–43 (2005).

22. Martínez, P. et al. Identification of the Major Sex-Determining Region of Turbot (Scophthalmus maximus). *Genetics* **183**, 1443–1452 (2009).

23. Shao, C.-W. et al. Construction of Two BAC Libraries from Half-Smooth Tongue Sole Cynoglossus semilaevis and Identification of Clones Containing Candidate Sex-DeterminationGenes. *Mar. Biotechnol.* **12**, 558–568 (2010).

24. Phillips, R. B., Park, L. K. & Naish, K. A. Assignment of Chinook Salmon (Oncorhynchus tshawytscha) Linkage Groups to Specific Chromosomes Reveals a Karyotype with Multiple Rearrangements of the Chromosome Arms of Rainbow Trout (Oncorhynchus mykiss). *G3 Genes* [Genomes][Genetics **3**, 2289–2295 (2013).

25. Iturra, P., Lam, N., dela Fuente, M., Vergara, N. & Medrano, J. F. Characterization of sex chromosomes in rainbow trout and coho salmon using fluorescence *in situ* hybridization (FISH). *Genetica* **111**, 125–131 (2001).

26. Larson, W. A., McKinney, G. J., Seeb, J. E. & Seeb, L. W. Identification and Characterization of Sex-Associated Loci in Sockeye Salmon Using Genotyping-by-Sequencing and Comparison with a Sex-Determining Assay Based on the sdY Gene. *J. Hered.* **107**, 559–566 (2016).

27. Ross, J. A., Urton, J. R., Boland, J., Shapiro, M. D. & Peichel, C. L. Turnover of Sex Chromosomes in the Stickleback Fishes (Gasterosteidae). *PLoS Genet.* **5**, e1000391 (2009).

28. Henning, F., Moysés, C. B., Calcagnotto, D., Meyer, A. & de Almeida-Toledo, L. F. Independent fusions and recent origins of sex chromosomes in the evolution and diversification of glass knife fishes (Eigenmannia). *Heredity* (Edinb). **106**, 391–400 (2011).

29. Palaiokostas, C. et al. A new SNP-based vision of the genetics of sex determination in European sea bass (Dicentrarchus labrax). *Genet. Sel. Evol.* **47**, 68 (2015).

30. Mylonas, C. C., Zohar, Y., Pankhurst, N. & Kagawa, H. Reproduction and broodstock management. in *Sparidae: Biology and Aquaculture of Gilthead Sea Bream and Other Species* (eds Pavlidis, M. A. & Mylonas, C. C.) 95–131 (Blackwell Publishing Ltd, 2011).

31. Buxton, C. D. & Garratt, P. A. Alternative reproductive styles in seabreams (Pisces: Sparidae). *Environ. Biol. Fishes* **28**, 113–124 (1990).

32. Wu, G.-C. & Chang, C.-F. The switch of secondary sex determination in protandrous black porgy, Acanthopagrus schlegeli. *Fish Physiol. Biochem.* **39**, 33–38 (2013).

33. He, C. L. et al. Differential Dmrt1 transcripts in gonads of the protandrous black porgy, Acanthopagrus schlegeli. *Cytogenet. Genome Res.* **101**, 309–313 (2003).

34. Loukovitis, D. et al. Quantitative trait loci for body growth and sex determination in the hermaphrodite teleost fish Sparus aurata L. *Anim Genet* **43**, 753–759 (2012).

35. Manousaki, T. et al. The sex-specific transcriptome of the hermaphrodite sparid sharpsnout seabream (Diplodus puntazzo). *BMC Genomics* **15**, 655 (2014).

36. Erisman, B. E., Petersen, C. W., Hastings, P. A. & Warner, R. R. Phylogenetic Perspectives on the Evolution of Functional Hermaphroditism in Teleost Fishes. *Integr. Comp. Biol.* https://doi.org/10.1093/icb/ict077 (2013).

37. Ravi, P., Jiang, J., Liew, W. C. & Orbán, L. Small-scale transcriptomics reveals differences among gonadal stages in Asian seabass (Lates calcarifer). *Reprod. Biol. Endocrinol.* **12**, 1–14 (2014).

38. Liu, H. et al. Large-scale transcriptome sequencing reveals novel expression patterns for key sex-related genes in a sex-changing fish. *Biol. Sex Differ.* **6**, 26 (2015).

39. Casas, L. et al. Sex Change in Clownfish: Molecular Insights from Transcriptome Analysis. *Sci. Rep.* **6**, 35461 (2016).

40. Basurco, B., Lovatelli, A. & Garcia, B. Current status of Sparidae aquaculture. in *Sparidae: Biology and Aquaculture of Gilthead Sea Bream and Other Species* (eds Pavlidis, M. A. & Mylonas, C. C.) 1–50 (Blackwell Publishing Ltd).

41. Bauchot, M. L. & Hureau, J. C. *Sparidae. Fishes of the North-eastern Atlantic and the Mediterranean* **2** (1986).

42. Shapiro, D. Y. Plasticity of gonadal development and protandry in fishes. *J. Exp. Zool.* **261**, 194–203 (1992).

43. Warner, R. R., Robertson, D. R. & Leigh, E. G. Sex change and sexual selection. *Science* (80-). **190**, 633 (1975).

44. Cataudella, S., Riz, P. P. & Sola, L. A chromosome study of eight Mediterranean species of Sparidae (pisces, perciformes). *Genetica* **54**, 155–159 (1980).

45. Manousaki, T. et al. Exploring a Nonmodel Teleost Genome Through RAD Sequencing—Linkage Mapping in Common Pandora, Pagellus erythrinus and Comparative Genomic Analysis. *G3 Genes* [Genomes][Genetics **6**, 509–519 (2016).

46. Avise, J. C. & Mank, J. E. Evolutionary Perspectives on Hermaphroditism in Fishes. *Sex. Dev.* **3**, 152–163 (2009).

47. Godwin, J. Neuroendocrinology of sexual plasticity in teleost fishes. *Front. Neuroendocrinol.* **31**, 203–216 (2010).

48. Larson, E. T. Neuroendocrine regulation in sex-changing fishes. in *Hormones and Reproduction of Vertebrates: Fishes* (eds Norris, D. O. & Lopez, K. H.) **1**, 149–168 (Elsevier, 2011).

49. Wang, Z., Gerstein, M. & Snyder, M. RNA-Seq: a revolutionary tool for transcriptomics. *Nat Rev Genet* **10** (2009).

50. Wilson, C. A. & Davies, D. C. The control of sexual differentiation of the reproductive system and brain. *Reproduction* **133**, 331–359 (2007).

51. Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. *Nat Methods* **9** (2012).

52. Consortium, T. U. UniProt: a hub for protein information. *Nucleic Acids Res.* **43**, D204–D212 (2015).

53. Quevillon, E. et al. InterProScan: protein domains identifier. *Nucleic Acids Res.* **33**, W116–W120 (2005).

54. Conesa, A. et al. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. *Bioinformatics* **21**, 3674–3676 (2005).

55. Ye, J. et al. WEGO: a web tool for plotting GO annotations. *Nucleic Acids Res.* **34**, W293–W297 (2006).

56. Grath, S. & Parsch, J. Sex-Biased Gene Expression. *Annu. Rev. Genet.* https://doi.org/10.1146/annurev-genet-120215-035429 (2015).

57. Böhne, A., Sengstag, T. & Salzburger, W. Comparative Transcriptomics in East African Cichlids Reveals Sex- and Species-Specific Expression and New Candidates for Sex Differentiation in Fishes. *Genome Biol Evol* **6**, 2567–2585 (2014).

58. O'Connell, L. A., Ding, J. H. & Hofmann, H. A. Sex differences and similarities in the neuroendocrine regulation of social behavior in an African cichlid fish. *Horm. Behav.* **64**, 468–476 (2013).

59. Kanamori, A., Naruse, K., Mitani, H., Shima, A. & Hori, H. Genomic organization of ZP domain containing egg envelope genes in medaka (Oryzias latipes). *Gene* **305**, 35–45 (2003).

60. Modig, C. et al. Molecular Characterization and Expression Pattern of Zona Pellucida Proteins in Gilthead Seabream (Sparus aurata). *Biol. Reprod.* **75**, 717–725 (2006).
61. Hyllner, S. J., Westerlund, L., Olsson, P. E. & Schopen, A. Cloning of rainbow trout egg envelope proteins: members of a unique group of structural proteins. *Biol. Reprod.* **64**, 805–811 (2001).

62. Nair, R. et al. Neurobeachin regulates neurotransmitter receptor trafficking to synapses. *J. Cell Biol.* **200**, 61 (2013).

63. Windley, S. P. & Wilhelm, D. Signaling Pathways Involved in Mammalian Sex Determination and Gonad Development. *Sex. Dev.* **9**, 297–315 (2015).

64. Taziaux, M., Keller, M., Balthazart, J. & Bakker, J. Rapid activation of phosphorylated MAPK after sexual stimulation in male mice. *Neuroreport* **22**, 294–298 (2011).

65. Burmeister, S. S. & Fernald, R. D. Evolutionary conservation of the egr-1 immediate-early gene response in a teleost. *J. Comp. Neurol.* **481**, 220–232 (2005).

66. Desjardins, J. K., Klausner, J. Q. & Fernald, R. D. Female genomic response to mate information. *Proc. Natl. Acad. Sci. USA* **107**, 21176–21180 (2010).

67. Le Page, Y. et al. Aromatase, brain sexualization and plasticity: the fish paradigm. *Eur. J. Neurosci.* **32**, 2105–2115 (2010).

68. Mank, J. E., Hultin-Rosenberg, L., Webster, M. T. & Ellegren, H. The unique genomic properties of sex-biased genes: Insights from avian microarray data. *BMC Genomics* **9**, 148 (2008).

69. Xu, D. et al. The testis and ovary transcriptomes of the rock bream (Oplegnathus fasciatus): A bony fish with a unique neo Y chromosome. *Genomics Data* **7**, 210–213 (2016).

70. Machado, M. P., Matos, I., Grosso, A. R., Schartl, M. & Coelho, M. M. Non-canonical expression patterns and evolutionary rates of sex-biased genes in a seasonal fish. *Mol. Reprod. Dev.* n/a-n/a https://doi.org/10.1002/mrd.22752 (2016).

71. Singh, R. S. & Kulathinal, R. J. Male sex drive and the masculinization of the genome. *BioEssays* **27**, 518–525 (2005).

72. Volff, J. N., Nanda, I., Schmid, M. & Schartl, M. Governing Sex Determination in Fish: Regulatory Putsches and Ephemeral Dictators. *Sex. Dev.* **1**, 85–99 (2007).

73. Smith, C. & Wootton, R. J. The remarkable reproductive diversity of teleost fishes. *Fish Fish.* **17**, 1208–1215 (2016).

74. Herpin, A. et al. Divergent Expression Regulation of Gonad Development Genes in Medaka Shows Incomplete Conservation of the Downstream Regulatory Network of Vertebrate Sex Determination. *Mol. Biol. Evol.* **30**, 2328–2346 (2013).

75. Böhne, A., Heule, C., Boileau, N. & Salzburger, W. Expression and Sequence Evolution of Aromatase cyp19a1 and Other Sexual Development Genes in East African Cichlid Fishes. *Mol. Biol. Evol.* **30**, 2268–2285 (2013).

76. Bleil, J. D. & Wassarman, P. M. Structure and function of the zona pellucida: Identification and characterization of the proteins of the mouse oocyte's zona pellucida. *Dev. Biol.* **76**, 185–202 (1980).

77. Shimizu, S., Tsuji, M. & Dean, J. *In vitro* biosynthesis of three sulfated glycoproteins of murine zonae pellucidae by oocytes grown in follicle culture. *J. Biol. Chem.* **258**, 5858–5863 (1983).

78. Liang, L., Soyal, S. M. & Dean, J. FIGalpha, a germ cell specific transcription factor involved in the coordinate expression of the zona pellucida genes. *Development* **124**, 4939 (1997).

79. Wu, G.-C. et al. Testicular dmrt1 Is Involved in the Sexual Fate of the Ovotestis in the Protandrous Black Porgy. *Biol. Reprod.* **86**, 1–11 (2012). 41.

80. Miyake, Y., Sakai, Y. & Kuniyoshi, H. Molecular Cloning and Expression Profile of Sex-Specific Genes, Figla and Dmrt1, in the Protognous Hermaphroditic Fish, Halichoeres Poecilopterus. *Zoolog. Sci.* **29**, 690–710 (2012).

81. Benayoun, B. A., Dipietromaria, A., Bazin, C. & Veitia, R. A. FOXL2: At the Crossroads of Female Sex Determination and Ovarian Function. in *Forkhead Transcription Factors: Vital Elements in Biology and Medicine* (ed. Maiese, K.) 207–226, https://doi.org/10.1007/978-1-4419-1599-3_16 (Springer New York, 2010).

82. Nakamoto, M. et al. Gonadal sex differentiation and expression of Sox9a2, Dmrt1 and Foxl2 in Oryzias luzonensis. *genesis* **47**, 289–299 (2009).

83. Wang, D. S., Kobayashi, T., Zhou, L. Y. & Nagahama, Y. Molecular cloning and gene expression of Foxl2 in the Nile tilapia, Oreochromis niloticus. *Biochem. Biophys. Res. Commun.* **320**, 83–89 (2004).

84. Nakamoto, M., Matsuda, M., Wang, D.-S., Nagahama, Y. & Shibata, N. Molecular cloning and analysis of gonadal expression of Foxl2 in the medaka, Oryzias latipes. *Biochem. Biophys. Res. Commun.* **344**, 353–361 (2006).

85. Liu, Z. et al. Molecular cloning of doublesex and mab-3-related transcription factor 1, forkhead transcription factor gene 2 and two types of cytochrome P450 aromatase in Southern catfish and their possible roles in sex differentiation. *J. Endocrinol.* **194**, 223–241 (2007).

86. Guiguen, Y., Fostier, A., Piferrer, F. & Chang, C.-F. Ovarian aromatase and estrogens: A pivotal role for gonadal sex differentiation and sex change in fish. *Gen. Comp. Endocrinol.* **165**, 352–366 (2010).

87. Kitano, T., Takamune, K., Kobayashi, T., Nagahama, Y. & Abe, S. I. Suppression of P450 aromatase gene expression in sex-reversed males produced by rearing genetically female larvae at a high water temperature during a period of sex differentiation in the Japanese flounder (Paralichthys olivaceus). *J. Mol. Endocrinol.* **23**, 167–176 (1999).

88. Uchida, D., Yamashita, M., Kitano, T. & Iguchi, T. An aromatase inhibitor or high water temperature induce oocyte apoptosis and depletion of P450 aromatase activity in the gonads of genetic female zebrafish during sex-reversal. *Comp. Biochem. Physiol. Part A Mol. Integr. Physiol.* **137**, 11–20 (2004).

89. Karube, M. et al. Characterization and expression profile of the ovarian cytochrome P-450 aromatase (cyp19A1) gene during thermolabile sex determination in pejerrey, Odontesthes bonariensis. *J. Exp. Zool. Part A Ecol. Genet. Physiol.* **307A**, 625–636 (2007).

90. Wu, G.-C. et al. Sex differentiation and sex change in the protandrous black porgy, Acanthopagrus schlegeli. *Gen. Comp. Endocrinol.* **167**, 417–421 (2010).

91. Zhou, L. & Gui, J. F. Molecular mechanisms underlying sex change in hermaphroditic groupers. *Fish Physiol Biochem* **36**, 181–193 (2010).

92. Baroiller, J. F. & D'Cotta, H. The Reversible Sex of Gonochoristic Fish: Insights and Consequences. *Sex. Dev.* **10**, 242–266 (2016).

93. Shao, C. et al. Epigenetic modification and inheritance in sexual reversal of fish. *Genome Res.* https://doi.org/10.1101/gr.162172.113 (2014).

94. Kobayashi, Y., Horiguchi, R., Miura, S. & Nakamura, M. Sex- and tissue-specific expression of P450 aromatase (cyp19a1a) in the yellowtail clownfish, Amphiprion clarkii. *Comp. Biochem. Physiol. Part A Mol. Integr. Physiol.* **155**, 237–244 (2010).

95. Piferrer, F. Epigenetics of sex determination and gonadogenesis. *Dev Dyn* **242**, 360–370 (2013).

96. Ospina-Álvarez, N. & Piferrer, F. Temperature-Dependent Sex Determination in Fish Revisited: Prevalence, a Single Sex Ratio Response Pattern and Possible Effects of Climate Change. *PLoS One* **3**, e2837 (2008).

97. Sandra, G.-E. & Norma, M.-M. Sexual determination and differentiation in teleost fish. *Rev. Fish Biol. Fish.* **20**, 101–121 (2009).

98. Ogiwara, K., Fujimori, C., Rajapakse, S. & Takahashi, T. Characterization of Luteinizing Hormone and Luteinizing Hormone Receptor and Their Indispensable Role in the Ovulatory Process of the Medaka. *PLoS One* **8**, e54482 (2013).

99. Ijiri, S. et al. Sexual Dimorphic Expression of Genes in Gonads During Early Differentiation of a Teleost Fish, the Nile Tilapia Oreochromis niloticus. *Biol. Reprod.* **78**, 333–341 (2008).

100. Baron, D., Houlgatte, R., Fostier, A. & Guiguen, Y. Expression profiling of candidate genes during ovary-to-testis transdifferentiation in rainbow trout masculinized by androgens. *Gen. Comp. Endocrinol.* **156**, 369–378 (2008).

101. Rodríguez-Mari, A. et al. Retinoic acid metabolic genes, meiosis and gonadal sex differentiation in zebrafish. *PLoS One* **8** (2013).

102. Feng, R. et al. Retinoic acid homeostasis through aldh1a2 and cyp26a1 mediates meiotic entry in Nile tilapia (Oreochromis niloticus). *Sci. Rep.* **5**, 10131 (2015).
103. Yamaguchi, T. & Kitano, T. High temperature induces cyp26b1 mRNA expression and delays meiotic initiation of germ cells by increasing cortisol levels during gonadal sex differentiation in Japanese flounder. *Biochem. Biophys. Res. Commun.* **419**, 287–292 (2012).

104. Dong, J. *et al*. Growth differentiation factor-9 is required during early ovarian folliculogenesis. *Nature* **383**, 531–535 (1996).

105. Yoshinaga, N. *et al*. Sexually dimorphic expression of a teleost homologue of Müllerian inhibiting substance during gonadal sex differentiation in Japanese flounder, *Paralichthys olivaceus*. *Biochem. Biophys. Res. Commun.* **322**, 508–513 (2004).

106. von Hofsten, J. & Olsson, P. E. Zebrafish sex determination and differentiation: involvement of FTZ-F1 genes. *Reprod Biol Endocrinol* **3**, 63 (2005).

107. Forconi, M. *et al*. Characterization of Sex Determination and Sex Differentiation Genes in Latimeria. *PLoS One* **8**, e56006 (2013).

108. Pevny, L. H. & Lovell-Badge, R. Sox genes find their feet. *Curr. Opin. Genet.* **Dev. 7**, 338–344 (1997).

109. Sutton, E. *et al*. Identification of SOX3 as an XX male sex reversal gene in mice and humans. *J. Clin. Invest.* **121**, 328–341 (2011).

110. Raverot, G., Weiss, J., Park, S. Y., Hurley, L. & Jameson, J. L. Sox3 expression in undifferentiated spermatogonia is required for the progression of spermatogenesis. *Dev. Biol.* **283**, 215–225 (2005).

111. Shin, H. S., An, K. W., Park, M. S., Jeong, M. H. & Choi, C. Y. Quantitative mRNA expression of sox3 and DMRT1 during sex reversal and expression profiles after GnRHα administration in black porgy, *Acanthopagrus schlegeli*. *Comp. Biochem. Physiol. Part B Biochem. Mol. Biol.* **154**, 150–156 (2009).

112. Yao, B., Zhou, L., Wang, Y., Xia, W. & Gui, J.-F. Differential expression and dynamic changes of SOX3 during gametogenesis and sex reversal in protogynous hermaphroditic fish. *J. Exp. Zool. Part A Ecol. Genet. Physiol.* **307A**, 207–219 (2007).

113. Jiang, T., Hou, C.-C., She, Z.-Y. & Yang, W.-X. The SOX gene family: function and regulation in testis determination and male fertility maintenance. *Mol. Biol. Rep.* **40**, 2187–2194 (2013).

114. Yokoi, H. *et al*. sox9 in a teleost fish, medaka (Oryzias latipes): Evidence for diversified function of Sox9 in gonad differentiation. *Mol. Reprod. Dev.* **63**, 5–16 (2002).

115. Chiang, E. F. L. *et al*. Two Sox9 Genes on Duplicated Zebrafish Chromosomes: Expression of Similar Transcription Activators in Distinct Sites. *Dev. Biol.* **231**, 149–163 (2001).

116. Raymond, C. S. *et al*. Evidence for evolutionary conservation of sex-determining genes. *Nature* **391**, 691–695 (1998).

117. Herpin, A. & Schartl, M. DMRT1 genes at the crossroads: a widespread and central class of sexual development factors in fish. *FEBS J.* **278**, 1010–1019 (2011).

118. Kopp, A. DMRT genes in the development and evolution of sexual dimorphism. *Trends Genet.* **28**, 175–184 (2012).

119. Yoshizawa, A. *et al*. Zebrafish Dmrta2 regulates neurogenesis in the telencephalon. *Genes to Cells* **16**, 1097–1109 (2011).

120. Lourenço, R., Lopes, S. S. & Saúde, L. Left-Right Function of dmrt2 Genes Is Not Conserved between Zebrafish and Mouse. *PLoS One* **5**, e14438 (2011).

121. Johnsen, H. & Andersen, Ø. Sex dimorphic expression of five dmrt genes identified in the Atlantic cod genome. The fish-specific *dmrt2b* diverged from *dmrt2a* before the fish whole-genome duplication. *Gene* **505**, 221–232 (2012).

122. Marchand, O. *et al*. DMRT1 expression during gonadal differentiation and spermatogenesis in the rainbow trout, Oncorhynchus mykiss. *Biochim. Biophys. Acta - Gene Struct. Expr.* **1493**, 180–187 (2000).

123. Berbejillo, J. *et al*. Expression and phylogeny of candidate genes for sex differentiation in a primitive fish species, the Siberian sturgeon, *Acipenser baerii*. *Mol. Reprod. Dev.* **79**, 504–516 (2012).

124. Fernandino, J. I. *et al*. Dimorphic Expression of dmrt1and cyp19a1 (Ovarian Aromatase) during Early Gonadal Development in Pejerrey, Odontesthes bonariensis. *Sex. Dev.* **2**, 316–324 (2008).

125. Picard, M. A.-L. *et al*. The roles of Dmrt (Double sex/Male-abnormal-3 Related Transcription factor) genes in sex determination and differentiation mechanisms: Ubiquity and diversity across the animal kingdom. *C. R. Biol.* **338**, 451–462 (2015).

126. Higa, M., Ogasawara, K., Sakaguchi, A., Nagahama, Y. & Nakamura, M. Role of steriod hormones in sex change of protogynous wrasse. *Fish Physiol. Biochem.* **28**, 149–150 (2003).

127. Wu, G.-C., Du, J.-L., Lee, Y.-H., Lee, M.-F. & Chang, C.-F. Current Status of Genetic and Endocrine Factors in the Sex Change of Protandrous Black Porgy, *Acanthopagrus schlegeli* (Teleostean). *Ann. N. Y. Acad. Sci.* **1040**, 206–214 (2005).

128. Guidelines for the treatment of animals in behavioural research and teaching. *Anim. Behav.* **53**, 229–234 (1997).

129. EU, Directive 2010/63/EU of the European parliament and the council of 22 September 2010 on the protection of animals used for scientific purposes. Official Journal of the European Union L 276/33, Animal protection (2010).

130. Metcalfe, J. D. & Craig, J. F. Ethical justification for the use and treatment of fishes in research: an update. *J. Fish Biol.* **78**, 393–394 (2011).

131. Diaz de Cerio, O., Rojo-Bartolomé, I., Bizarro, C., Ortiz-Zarragoitia, M. & Cancio, I. 5S rRNA and Accompanying Proteins in Gonads: Powerful Markers to Identify Sex and Reproductive Endocrine Disruption in Fish. *Environ. Sci. Technol.* **46**, 7763–7771 (2012).

132. Shen, Z.-G., Yao, H., Guo, L., Li, X.-X. & Wang, H.-P. Ribosome RNA Profiling to Quantify Ovarian Development and Identify Sex in Fish. *Sci. Rep.* **7**, 4196 (2017).

133. FastQC. In. http://www.bioinformatics.babraham.ac.uk/projects/fastqc.

134. Scythe. In. https://github.com/vsbuffalo/scythe.

135. Joshi, N. A. & Fass, J. N. sickle. In. https://github.com/najoshi/sickle.

136. Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. *Bioinformatics* **30** (2014).

137. Schmieder, R. & Edwards, R. Quality control and preprocessing of metagenomic datasets. *Bioinformatics*. https://doi.org/10.1093/bioinformatics/btr026 (2011).

138. Grabherr, M. G. *et al*. Trinity: reconstructing a full-length transcriptome without a genome from RNA-Seq data. *Nat. Biotechnol.* **29**, 644–652 (2011).

139. Li, B. & Dewey, C. N. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. *BMC Bioinformatics* **12** (2011).

140. Haas, B. J. *et al*. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. *Nat. Protoc.* **8**, 1494–1512 (2013).

141. Li, H. *et al*. The Sequence Alignment/Map format and SAMtools. *Bioinformatics* **25** (2009).

142. Simão, F. A., Waterhouse, R. M., Ioannidis, P., Kriventseva, E. V & Zdobnov, E. M. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. *Bioinformatics* https://doi.org/10.1093/bioinformatics/btv351 (2015).

143. Altschul, S. F., Gish, W., Miller, W., Myers, E. W. & Lipman, D. J. Basic local alignment search tool. *J. Mol. Biol.* **215**, 403–410 (1990).

144. Lagnel, J., Tsigenopoulos, C. S. & Iliopoulos, I. NOBLAST and JAMBLAST: New Options for BLAST and a Java Application Manager for BLAST results. *Bioinformatics* **25**, 824–826 (2009).

145. Hunter, S. *et al*. InterPro: the integrative protein signature database. *Nucleic Acids Res.* **37**, D211–D215 (2009).

146. Rice, P., Longden, I. & Bleasby, A. EMBOSS: The European Molecular Biology Open Software Suite. *Trends Genet.* **16**, 276–277 (2016).

147. Myhre, S., Tveit, H., Mollestad, T. & Lægreid, A. Additional Gene Ontology structure for improved biological reasoning. *Bioinformatics* **22**, 2020–2027 (2006).

148. GO Slim. In. http://geneontology.org/page/go-slim-and-subset-guide.

149. Kanehisa, M. & Goto, S. KEGG: Kyoto Encyclopedia of Genes and Genomes. *Nucleic Acids Res.* **28**, 27–30 (2000).
150. Bourgon, R., Gentleman, R. & Huber, W. Independent filtering increases detection power for high-throughput experiments. *Proc. Natl. Acad. Sci. USA* **107**, 9546–9551 (2010).

151. Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq. *2. Genome Biol* **15** (2014).

152. Gentleman, R. C. et al. Bioconductor: open software development for computational biology and bioinformatics. *Genome Biol* **5** (2004).

153. Jolliffe, I. T. *Principal component analysis*. (Springer, 2002).

154. Warnes, G. R. et al. gplots: Various R programming tools for plotting data. *R Packag. version* **2** (2009).

155. Yates, A. et al. Ensembl 2016. *Nucleic Acids Res.* https://doi.org/10.1093/nar/gkv1157 (2015).

156. Overbeek, R., Fonstein, M., D'Souza, M., Pusch, G. D. & Maltsev, N. The use of gene clusters to infer functional coupling. *Proc. Natl. Acad. Sci. USA* **96**, 2896–2901 (1999).

## Acknowledgements

Financial support for this study has been provided by the Ministry of Education and Religious Affairs, under the Call “ARISTEIA I” of the National Strategic Reference Framework 2007–2013 (SPARCOMP, #36), co-funded by the EU and the Hellenic Republic through the European Social Fund. The authors wish to thank Mrs. Margaret Eleftheriou for proof-reading and English editing of the manuscript. We would also like to thank Haris Apostolidis, Constantinos Vatikiotis and Irini Sigelaki for their valuable help in sampling.

## Author Contributions

A.T., T.M., C.C.M., M.P. and C.S.T. designed experiments. A.T., C.S.T., N.P., A.S., C.C.M. conducted sampling. A.T. performed laboratory experiments, analyzed data with J.L. and drafted the manuscript. C.S.T. conceived and managed the project. All authors participated in data interpretation, read and approved the final version of the manuscript.

## Additional Information

Supplementary information accompanies this paper at https://doi.org/10.1038/s41598-018-21992-9.

**Competing Interests:** The authors declare no competing interests.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2018