# The Gene Toolkit Implicated in Functional Sex in Sparidae Hermaphrodites: Inferences From Comparative Transcriptomics

Alexandros Tsakogiannis, Tereza Manousaki, Jacques Lagnel, Nikolaos Papanikolaou, Nikos Papandroulakis, Constantinos C. Mylonas, Costas S. Tsigenopoulos

► To cite this version:

Alexandros Tsakogiannis, Tereza Manousaki, Jacques Lagnel, Nikolaos Papanikolaou, Nikos Papandroulakis, et al.. The Gene Toolkit Implicated in Functional Sex in Sparidae Hermaphrodites: Inferences From Comparative Transcriptomics. Frontiers in Genetics, 2019, 9, 10.3389/fgene.2018.00749. [hal-02620563](https://hal.inrae.fr/hal-02620563)

HAL Id: hal-02620563

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

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
<figure><img src="image_2.png" /></figure>

# The Gene Toolkit Implicated in Functional Sex in Sparidae Hermaphrodites: Inferences From Comparative Transcriptomics

Alexandros Tsakogiannis<fnref n="1" />,<fnref n="2" />, Tereza Manousaki<fnref n="1" />, Jacques Lagnel<fnref n="1" />, Nikolaos Papanikolaou³, Nikos Papandroulakis<fnref n="1" />, Constantinos C. Mylonas<fnref n="1" /> and Costas S. Tsigenopoulos<fnref n="1" />*

Sex-biased gene expression is the mode through which sex dimorphism arises from a nearly identical genome, especially in organisms without genetic sex determination. Teleost fishes show great variations in the way the sex phenotype forms. Among them, Sparidae, that might be considered as a model family displays a remarkable diversity of reproductive modes. In this study, we sequenced and analyzed the sex-biased transcriptome in gonads and brain (the tissues with the most profound role in sexual development and reproduction) of two sparids with different reproductive modes: the gonochoristic common dentex, *Dentex dentex*, and the protandrous hermaphrodite gilthead seabream, *Sparus aurata*. Through comparative analysis with other protogynous and rudimentary protandrous spard transcriptomes already available, we put forward common male and female-specific genes and pathways that are probably implicated in sex-maintenance in this fish family. Our results contribute to the understanding of the complex processes behind the establishment of the functional sex, especially in hermaphrodite species and set the groundwork for future experiments by providing a gene toolkit that can improve efforts to control phenotypic sex in finfish in the ever-increasingly important field of aquaculture.

**Keywords:** sex-biased gene expression, hermaphroditism, RNA-seq, comparative transcriptomics, gonads, brain

## OPEN ACCESS

**Edited by:**
Nguyen Hong Nguyen,
University of the Sunshine Coast,
Australia

**Reviewed by:**
Shikai Liu,
Ocean University of China, China
Jun Hong Xia,
Sun Yat-sen University, China

**\*Correspondence:**
Costas S. Tsigenopoulos
tsigeno@hcmr.gr

**Specialty section:**
This article was submitted to
Livestock Genomics,
a section of the journal
Frontiers in Genetics

**Received:** 31 July 2018
**Accepted:** 31 December 2018

**Published:** 18 January 2019

**Citation:**
Tsakogiannis A, Manousaki T, Lagnel J, Papanikolaou N, Papandroulakis N, Mylonas CC and Tsigenopoulos CS (2019) The Gene Toolkit Implicated in Functional Sex in Sparidae Hermaphrodites: Inferences From Comparative Transcriptomics. *Front. Genet.* 9:749.
doi: 10.3389/fgene.2018.00749

## INTRODUCTION

Sexual dimorphism is widespread in nature and describes the complete set of differences between sexes (Jazin and Cahill, 2010; Angelopoulou et al., 2012; Gamble and Zarkower, 2012). Despite differences in morphology and behavior, females and males are genetically nearly identical. However, many genes are more actively transcribed in one sex than the other, and sex-biased expression is considered to be the way by which sexual dimorphisms can arise from the genome (Ranz et al., 2003; Ellegren and Parsch, 2007; Mank et al., 2008; Parsch and Ellegren, 2013). The establishment of sexual phenotype is accomplished through the fine-tuned action and collaboration of numerous genes whose expression is prompted or repressed at preset crucial time periods during development (Angelopoulou et al., 2012).

<footnote n="1">¹ Hellenic Centre for Marine Research, Institute of Marine Biology, Biotechnology and Aquaculture, Heraklion, Greece,</footnote>

<footnote n="2">² Department of Biology, University of Crete, Heraklion, Greece, ³ Department of Medicine, University of Crete, Heraklion, Greece</footnote>
Teleost fishes show great variations in the way the sex phenotype forms. Sexual differentiation in this group of vertebrates is evolutionarily flexible, as different sexual patterns exist, ranging from gonochorism to hermaphroditism. Studies in fish have already redefined our understanding of the complexity and plasticity of the sex determination (SD) process, yet detailed genetic information on SD genes is currently available for only a limited number of species (e.g., *DmY/Dmrt1bY* in *Oryzias latipes* and in *O. curvinotus* (Matsuda et al., 2002), *Gsdf* in *O. luzonensis* (Myosho et al., 2012), *Amhy* in *Odontesthes hatcheri* (Hattori et al., 2012), *Amhr2* in *Takifugu rubripes*, *T. pardalis*, and *T. poecilonotus* (Kamiya et al., 2012), and *SdY* in *Oncorhynchus mykiss* (Yano et al., 2012). Furthermore, 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* (Kovács et al., 2000), the Nile tilapia *Oreochromis niloticus* (Ezaz et al., 2004), the rainbow trout *Oncorhynchus mykiss* (Felip et al., 2005), the turbot *Scophthalmus maximus* (Martínez et al., 2009), the half-smooth tongue sole *Cynoglossus semilaevis* (Shao et al., 2010), some salmonids (Iturra et al., 2001; Phillips et al., 2013; Larson et al., 2016), Gasterosteidae (Ross et al., 2009), species of the genus Eigenmanna (Henning et al., 2011), and the European sea bass *Dicentrarchus labrax* (Palaiokostas et al., 2015).

Sex control constitutes one important and favorably targeted area of research in aquaculture, as it has been one of the most important of the facilitators for large-scale industrial production (Budd et al., 2015). A large body of research has targeted sexual development in commercially important fish species in order to understand and control sex and be able to regulate sexual differentiation and maturation (Budd et al., 2015). Among aquacultured species, the Sparidae present a morphologically and ecologically highly diverse family of Perciform fishes. An exceptional characteristic of the species in this family is that they exhibit nearly every reproductive motif known in fishes, rendering them ideal for comparative studies. In addition to gonochoristic species (e.g., *Dentex dentex*), in which individuals retain the same sex through their lifetime, various forms of hermaphroditism have been recorded (Buxton and Garratt, 1990). Some species change from functional females to functional males (protogyny) (e.g., *Pagrus pagrus* and *Pagellus erythrinus*), while others function first as males and then change to females (protandry) (e.g., *Sparus aurata*). Furthermore, some species exhibit rudimentary hermaphroditism (e.g., *Diplodus puntazzo*), so-called “late gonochorism” (Buxton and Garratt, 1990), as individuals possess both male and female gonadal tissue during sex differentiation and before puberty, and one of the two gonadal types develops fully only after sexual maturation, with no evidence of further sex reversal.

With the advent of next-generation sequencing technologies, such as RNA-Seq (Wang et al., 2009), whole-transcriptome approaches have been used to study the molecular mechanisms of sexual differentiation in fish (Zhang et al., 2011; Forconi et al., 2013; Sun et al., 2013; Tao et al., 2013; Böhne et al., 2014; Liu et al., 2015; Casas et al., 2016). Regarding sparids, most studies on sex determination/differentiation and sex change have been conducted in protandrous species, such as the black

porgy *Acanthopagrus schlegelii* (He et al., 2003; Wu and Chang, 2013) and gilthead seabream *Sparus aurata* (Loukovitis et al., 2012). The main focus of these studies is to track differences in gene expression between males and females and discuss the expression patterns of important genes in the general context of sex determination/differentiation and sex-change. Recently we sequenced and analyzed the brain and gonad transcriptome of the rudimentary hermaphrodite sharpsnout seabream *Diplodus puntazzo* (Manousaki et al., 2014) and of two protogynous sparids, the red porgy *Pagrus pagrus* and the common pandora *Pagellus erythrinus* (Tsakogiannis et al., 2018). In this study, we sequenced and analyzed the sex-biased transcriptome of two sparids with different reproductive modes: the gonochoristic common dentex *Dentex dentex* and the protandrous gilthead seabream *Sparus aurata*. Then, we explored the sex-specific expression patterns in all five transcriptomes from these five sparid species exhibiting four different reproductive modes. Using RNA sequencing, we analyzed and compared differences between male and female transcriptomes in gonads and brain, the tissues with the most profound role in sexual development and reproduction (Vilain and McCabe, 1998; Wilson and Davies, 2007) to obtain a global view of sex-biased expression in these tissues, and through comparative analysis to reveal common male and female-specific genes and pathways that are possibly involved in the development of the functional sex phenotype. Our results contribute to the understanding of the complex processes of sex differentiation, especially on hermaphrodite species.

# MATERIALS AND METHODS

## Ethics Statement

Animal welfare was achieved according to the “Guidelines for the treatment of animals in behavioral research and teaching” (Anonymous, 1997) (see also Tsakogiannis et al., 2018). All fish utilized in the study were kept in registered and authorized facilities to maintain and perform animal experiments; rearing and sampling followed the guidelines of the Directive 2010/63/EU for the protection of animals used for experimental and other scientific purposes (Official Journal L276/33) (EU, 2010). In addition, experimental sampling protocols were approved by the IMBBC’s aquaculture department committee and methods were in accordance with relevant guidelines and regulations approved by the Hellenic 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). Laboratory personnel include accredited technicians by the Federation for Laboratory Animal Science Associations (FELASA).

## Sample Collection

Sample collection information for the rudimentary sharpsnout seabream, and the protogynous common pandora and red porgy has been described in detail in Manousaki et al. (2014) and Tsakogiannis et al. (2018), respectively. The same procedure was followed for sample collection of the common dentex and
the gilthead seabream. The fish used for the experiments were broodstock fish of wild origin in the Institute of Marine Biology, Biotechnology and Aquaculture (IMBBC, HCMR). The 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 and kept in RNA-later (Qiagen) following provider's guidelines. The fish were examined macroscopically for sexual maturation, based on the presence of releasable sperm for males or the presence of vitellogenin oocytes for females and/or visual inspection of the gonad morphology. Specifically, for common dentex, brains and gonads from five mature females (average weight = 1.04 kg, age = 40 months old) and six mature males (average weight = 1.76 kg, age = 46 months old) were sampled (six biological replicates for the male gonads, five biological replicates for the male brains and five biological replicates for the females per tissue), while for gilthead seabream the same tissues were sampled from four mature males (average weight = 0.68 kg, age = 28 months old) and four mature females (average weight = 2.2 kg, age = 52 months old) (four biological replicates for each sex per tissue). Especially for common dentex, a pool of few hundred larvae at the "hatching" stage (1 day post-hatch, dph) was also sampled to enrich the species first transcriptome assembly (Table 1). In total, common dentex and gilthead searbeam sequences were obtained from 38 samples which were used for RNA extraction, library preparation, and sequencing.

## RNA Extraction and Sequencing

Each one of the samples was 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. RNA quantity was measured spectrophotometrically with NanoDrop® ND-1000 (Thermo Scientific), while its quality was tested first on agarose gel (electrophoresis in 1.5% w/v) and then on Agilent Technologies 2100 Bioanalyzer (Agilent Technologies). Finally, each one of the 38 RNA samples of common dentex and gilthead seabream was used for mRNA paired-end library construction with the Illumina TruSeq™ RNA Sample Preparation Kits v2 following manufacturer's guidelines. For common dentex, 22 cDNA libraries were sequenced in one and a half lane of an Illumina HiSeq2500™ instrument and run with a 100 bp paired-end manner. All 16 gilthead seabream libraries were sequenced on an Illumina HiSeq4000™ instrument with the 150 bp paired-end method in one lane following the protocols of Illumina Inc. (San Diego, CA). The raw sequences obtained for each sample are presented in Table 1.

## Bioinformatic Analysis

Raw reads were pre-processed following the pipeline published in Ilias et al. (2015). Initially, we used Scythe [https://github.com/vsbuffalo/scythe] to identify adapter contamination setting prior contamination rate set in 0.1 "-p 0.1," Sickle to trim low quality reads (parameters "pe -g -t sanger -q 20 -l 45"), Trimmomatic (Bolger et al., 2014) (parameters "PE -phred33 ILLUMINACLIP:adapter_file.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:25 MINLEN:45 CROP:99")

<table><caption>**TABLE 1** | Samples and raw sequences of common dentex and gilthead seabream.</caption><thead><tr><th></th><th>Biological replicates</th><th>Samples</th><th>Sex</th><th>Tissue</th><th>Raw reads</th></tr></thead><tbody><tr><td rowspan="14">Common dentex</td><td rowspan="4">5 female brains</td><td>DdFB1</td><td>Female</td><td>Brain</td><td>19,384,462</td></tr><tr><td>DdFB2</td><td>Female</td><td>Brain</td><td>16,808,874</td></tr><tr><td>DdFB3</td><td>Female</td><td>Brain</td><td>25,074,840</td></tr><tr><td>DdFB4</td><td>Female</td><td>Brain</td><td>29,206,520</td></tr><tr><td rowspan="4">5 female gonads</td><td>DdFB5</td><td>Female</td><td>Brain</td><td>23,474,036</td></tr><tr><td>DdFG1</td><td>Female</td><td>Gonad</td><td>27,705,642</td></tr><tr><td>DdFG2</td><td>Female</td><td>Gonad</td><td>24,365,178</td></tr><tr><td>DdFG3</td><td>Female</td><td>Gonad</td><td>22,076,722</td></tr><tr><td rowspan="5">5 male brains</td><td>DdFG4</td><td>Female</td><td>Gonad</td><td>21,219,416</td></tr><tr><td>DdFG5</td><td>Female</td><td>Gonad</td><td>23,457,354</td></tr><tr><td>DdMB1</td><td>Male</td><td>Brain</td><td>24,070,998</td></tr><tr><td>DdMB2</td><td>Male</td><td>Brain</td><td>23,322,846</td></tr><tr><td>DdMB4</td><td>Male</td><td>Brain</td><td>20,056,858</td></tr><tr><td rowspan="5">6 male gonads</td><td>DdMB5</td><td>Male</td><td>Brain</td><td>15,196,762</td></tr><tr><td>DdMB6</td><td>Male</td><td>Brain</td><td>13,886,812</td></tr><tr><td>DdMG1</td><td>Male</td><td>Gonad</td><td>25,897,816</td></tr><tr><td>DdMG2</td><td>Male</td><td>Gonad</td><td>34,311,430</td></tr><tr><td>DdMG3</td><td>Male</td><td>Gonad</td><td>31,033,676</td></tr><tr><td rowspan="6">Gilthead seabream</td><td rowspan="2">1 larvae pool</td><td>DdHatch</td><td>NA</td><td>Whole body</td><td>20,765,736</td></tr><tr><td>22</td><td>-</td><td>-</td><td>512,822,392</td></tr><tr><td rowspan="2">4 female brains</td><td>SaFB1</td><td>Female</td><td>Brain</td><td>17,280,880</td></tr><tr><td>SaFB2</td><td>Female</td><td>Brain</td><td>17,739,060</td></tr><tr><td rowspan="5">4 female gonads</td><td>SaFB3</td><td>Female</td><td>Brain</td><td>20,802,438</td></tr><tr><td>SaFB4</td><td>Female</td><td>Brain</td><td>17,691,814</td></tr><tr><td>SaFG1</td><td>Female</td><td>Gonad</td><td>34,194,856</td></tr><tr><td>SaFG2</td><td>Female</td><td>Gonad</td><td>33,241,228</td></tr><tr><td>SaFG3</td><td>Female</td><td>Gonad</td><td>34,271,716</td></tr><tr><td rowspan="5">4 male brains</td><td>SaFG4</td><td>Female</td><td>Gonad</td><td>32,293,436</td></tr><tr><td>SaMB1</td><td>Male</td><td>Brain</td><td>20,938,852</td></tr><tr><td>SaMB2</td><td>Male</td><td>Brain</td><td>28,548,174</td></tr><tr><td>SaMB3</td><td>Male</td><td>Brain</td><td>25,808,438</td></tr><tr><td>SaMB4</td><td>Male</td><td>Brain</td><td>28,432,110</td></tr><tr><td rowspan="4">4 male gonads</td><td>SaMG1</td><td>Male</td><td>Gonad</td><td>33,073,506</td></tr><tr><td>SaMG2</td><td>Male</td><td>Gonad</td><td>31,482,028</td></tr><tr><td>SaMG3</td><td>Male</td><td>Gonad</td><td>32,654,702</td></tr><tr><td>SaMG4</td><td>Male</td><td>Gonad</td><td>31,861,858</td></tr><tr><td>Total</td><td></td><td>16</td><td></td><td></td><td>440,315,096</td></tr></tbody></table>

and Prinseq (Schmieder and Edwards, 2011) to remove poly A/T tails and low complexity reads.

The two new assemblies were produced using Trinity v2.0.6 (Grabherr et al., 2011) with default parameters (default kmer
25, min length 200 nucleotides). Downstream assembly filtering was conducted using the same criteria as in Manousaki et al. (2014) (IsoPct < 1 and FPKM < 0.3). To assess the quality of all five assemblies we used three criteria: (1) assessing the percentage of mapped reads of each species back to the assembly using Bowtie2 (Langmead and Salzberg, 2012) and Samtools (Li et al., 2009); (2) examining the similarity between the Trinity assembly and Swissprot database (UniProt Consortium, 2015) using blastx (Altschul et al., 1990) searches within the parallel NOBLAST program (Lagnel et al., 2009); and (3) quantifying the completeness of each assembly using BUSCO (search for “vertebrates” in OrthoDBv8) (Simão et al., 2015). The newly obtained assemblies were then functionally characterized as described in Tsakogiannis et al. (2018).

Differential expression (DE) analysis was performed at gene level by pairwise comparisons between male and female for brain and gonad samples in each species separately. The pre-processed reads of each species were aligned to the respective assembly with Bowtie2 (Langmead and Salzberg, 2012) and abundance was estimated with RSEM (Li and Dewey, 2011), 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 (Bourgon et al., 2010). To accomplish the same consensus in differential expression (DE) analysis among all species, we re-analyzed sharpsnout seabream to bring into line with the other four species. Count data were imported in R (version 3.2.0) and the analysis of differential expression conducted in DESeq2 (version: 1.8.1) (Love et al., 2014), an R bioconductor package (Gentleman et al., 2004). Each species brain and gonadal samples were grouped according to sex and expression was compared (male vs. female), following the developers' manual (false discovery rate, FDR, threshold of 0.05). Genes with an adjusted *p*-value <0.05 were considered significantly differentially expressed. Principal component analysis (PCA) (Jolliffe, 2002) and the heatmap.2 function in the gplots package (Warnes et al., 2009) were used to visualize global similarities and differences among the gonadal samples. Each species PCA plot was produced including all DE genes (Differentially Expressed genes) in the filtered (sum of counts in each row >15) count matrix.

Functional-enrichment analysis was performed to identify gene ontology (GO) terms and metabolic pathways significantly enriched in DE genes in gonads. Here, the term ‘enrichment’ stands for the comparison between the DE genes of the female (test-set) against the male (reference-set) and the gene-set was restricted to only those with |log₂FC| > 2. All enrichment analyses were implemented in BLAST2GO (Conesa et al., 2005) with default settings enabling a two-tailed test (FDR = 0.05).

To track down the expression patterns in known sex-related candidate genes in vertebrates and other fish taxa, we downloaded the respective genes of tilapia, medaka and spotted gar from Ensembl (release 84) (Yates et al., 2015) resulting in a set of 114 candidate genes (see also Tsakogiannis et al., 2018). To see if orthologs of these genes were present in our species transcriptomes, we implemented a reciprocal BLASTn best-hit approach (Overbeek et al., 1999). Potentially orthologous genes differentially

expressed in all species or were species-specific were identified.

The final step included a comparative analysis for male and female-biased orthologous gonadal genes. Following an in-house built pipeline, using perl scripts and the Vennt tool (http://drpowell.github.io/vennt/) we combined the expression data of *D. dentex*, *S. aurata*, *P. pagrus*, *P. erythrinus*, and *D. puntazzo* gonadal orthologous genes.

We performed a BLASTn between each pair of the five transcriptomes (e-value threshold 10⁻¹⁰), keeping the best reciprocal hits (bidirectional best blast hits), an approach followed also in Tsakogiannis et al. (2018). We then selected all transcripts that have orthologs in all five species. Next, we combined the expression information with the existing orthologous genes. Following that, we visualized and analyzed the data using the freely available Vennt tool (version v0.8.4, GPL v.3, http://drpowell.github.io/vennt), which generates dynamic Venn diagrams for differential gene expression. Based on the Vennt results, we concluded that an FDR value lower than 0.05 produces robust results. Finally, we plotted a Venn diagram for all five species using the R package “venn” (version 1.5, GPL v.3, https://cran.r-project.org/web/packages/venn/index.html). A more detailed description of parameters used in bioinformatic analysis can be found in **Supplementary Excel Table SET1**.

# RESULTS

## Pre-processing Treatment of Reads and Assembly Construction

Illumina sequencing of common dentex and gilthead seabream yielded 512,882,392 and 440,315,096 reads, respectively. The raw sequence data has been submitted to the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database (common dentex: BioProject ID PRJNA481721 and gilthead seabream: BioProject ID PRJNA415944).

Reads' pre-processing led to 272,911,232 and 368,540,942 high-quality paired reads for common dentex and gilthead seabream, respectively. High-quality reads were used to construct the initial transcriptome assemblies, consisting of 247,693 transcripts for common dentex (N50: 1,891 nucleotides, mean length: 870 nucleotides) and 292,100 for gilthead seabream (N50: 1,939 nucleotides, mean length: 872 nucleotides).

Following transcript quality control, the initial assemblies were limited to 129,012 contigs (N50: 2,592; mean length: 1,259 nucleotides) corresponding to 98,747 genes for common dentex and 119,677 contigs (N50: 2,806; mean length: 1,321 nucleotides) corresponding to 92,160 genes for gilthead seabream. The restricted final datasets combined with the respective datasets of the previously published sparids are shown in a comparative representation in **Figure 1**.

To evaluate the accuracy of the assembled transcripts, clean reads were aligned against the transcripts, reaching relatively high percentages in both common dentex (85%) and gilthead seabream (80%). Then, a sequence similarity search between each species and UniProtKB/SwissProt database disclosed 19,429
<figure><img src="image_3.png" /></figure>

**FIGURE 1** | Pipeline for building assemblies and differential expression profiles in brains and gonads of the two species sequenced here in addition to sharpsnout seabream (Manousaki et al., 2014), common pandora and red porgy (Tsakogiannis et al., 2018). Flow chart of the basic steps implemented from raw reads to the selection of final assembled loci and DE genes in brains and gonads of all species.

and 16,690 contigs with hits for common dentex and gilthead seabream, respectively.

Finally, BUSCO run (Supplementary Table ST1) indicated that both assemblies were almost complete with 85% complete matches of vertebrate orthologs, found in both common dentex and gilthead seabream assemblies.

## Transcriptome Annotation

A blastx search of both assembled transcripts against UNIPROT (UniProt Consortium, 2015; *e*-value cut-off <10⁻⁵) returned 26,509 and 24,326 genes with at least one hit for common dentex and gilthead seabream, respectively, and 24,248 and 23,961 genes with a protein domain match against InterProScan (Quevillon et al., 2005), respectively.

To identify possible functions, Gene Ontology (GO) assignments were used to classify the sequences. A total of 21,705 and 23,169 common dentex and gilthead seabream genes, respectively, were attributed with at least one GO term. The potential enzymes were characterized based on the chemical reaction they catalyze using the prediction of Enzyme Commission (EC) numbers for each sequence. In the case of common dentex 4,616 genes were classified to 1,037 different

EC numbers, while in the case of gilthead seabream, 6,050 genes were classified to 1,215 different EC numbers (Table 2).

## Global Expression Pattern of Brains and Gonads

The global gene expression pattern observed for brain and gonad tissues for the common dentex and gilthead seabream is summarized by implementing a principal component analysis (PCA) (Supplementary Figure SF1). The plots show the formation of three separate clusters for each species, representing the brains, the male gonads and the female gonads. While male and female gonads show a clear separation with a remarkable variation within the two groups, male or female brains are clustered together exhibiting very similar expression patterns. The same pattern of samples clustering is also shown through heatmap plots (Supplementary Figure SF2).

## Sex-Biased Gene Expression in Male and Female Brains

The analysis revealed minor differences in expression between sexes in the brain. For the purpose of the comparison, we kept only loci with a sum of expected counts in each row of the count
matrix >15 leading to 35,440 and 45,793 brain genes for common dentex and gilthead seabream, respectively. In general, brain DE genes showed low counts and low fold-change differences. Common dentex showed no differential expression in the brain. In gilthead seabream, 54 genes were differentially expressed with 16 reported as male-biased and 38 as female-biased (Figure 2).

## Sex-Biased Gene Expression in Male and Female Gonads

The analysis revealed great differences in expression between male gonads and female gonads. For the purpose of the

comparison, we kept only loci with a sum of expected counts in each row of the count matrix >15 keeping 39,587 and 41,368 for common dentex and gilthead seabream, respectively. More particularly, in the case of common dentex 22,762 genes were significantly differentially expressed with 15,264 over-expressed in males and 7,498 in females, while in the case of gilthead seabream 29,281 genes were differentially expressed with 18,751 over-expressed in males and 10,530 over-expressed in females (Figure 3).

For a better understanding of the biological processes underlying the nature of the functional male

<table><caption>TABLE 2 | Annotation summary of common dentex and gilthead seabream transcriptomes.</caption><thead><tr><th rowspan="3"></th><th colspan="8">Annotation summary</th></tr><tr><th colspan="4">Common dentex</th><th colspan="4">Gilthead seabream</th></tr><tr><th>contigs</th><th>% contigs</th><th>genes</th><th>% genes</th><th>contigs</th><th>% contigs</th><th>genes</th><th>% genes</th></tr></thead><tbody><tr><td>BLASTx against Uniprot</td><td>129,012</td><td></td><td>98,747</td><td></td><td>119,677</td><td></td><td>92,160</td><td></td></tr><tr><td>InterPro</td><td>46,471</td><td>36</td><td>26,509</td><td>27</td><td>39,049</td><td>33</td><td>24,326</td><td>26</td></tr><tr><td>With IPR</td><td>42,450</td><td>33</td><td>24,248</td><td>25</td><td>38,263</td><td>32</td><td>23,961</td><td>26</td></tr><tr><td>With ≥ 1GO</td><td>38,378</td><td>30</td><td>21,986</td><td>22</td><td>34,251</td><td>29</td><td>21,462</td><td>23</td></tr><tr><td>Blast2Go annotated</td><td>30,720</td><td>24</td><td>17,717</td><td>18</td><td>27,437</td><td>23</td><td>17,102</td><td>19</td></tr><tr><td>EC</td><td>34,595</td><td>27</td><td>21,705</td><td>22</td><td>37,535</td><td>31</td><td>23,169</td><td>25</td></tr><tr><td></td><td>7,136</td><td>6</td><td>4,611</td><td>5</td><td>9,877</td><td>8</td><td>6,050</td><td>7</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>FIGURE 2 | (A) Comparative view of the number of DE genes in brain between males and females in the two species. Venn-like diagram showing the DE genes in the two sexes in the brain of (A. 1.) common dentex and (A. 2.) gilthead seabream. Light-Blue: male-biased genes. Red: female-biased genes. Dark-Purple (intersection): un-biased genes. (B) Principal component analysis (PCA) plots in (B. 1.) common dentex and (B. 2.) gilthead seabream and (C) Heatmaps of the variance stabilized, transformed count data plots, for the DE genes in (C. 1.) common dentex and (C. 2.) gilthead seabream (M, Male; F, Female).</figcaption></figure>
Gonads

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

**FIGURE 3** | (A) Comparative view of the number of DE genes between males and females in the two species gonads. Venn-like diagram showing the DE genes in the two sexes in the gonads in (A. 1.) common dentex and (A. 2.) gilthead seabream. Light-Blue: male-biased genes. Red: female-biased genes. Dark-Purple (intersection): un-biased genes. (B) Principal component analysis (PCA) plots in (B. 1.) common dentex and (B. 2.) gilthead seabream and (C) Heatmaps of the variance stabilized, transformed count data plots, for the DE genes in (C. 1.) common dentex and (C. 2.) gilthead seabream (M, Male; F, Female).

or female gonadal phenotype, GO term enrichment analysis was implemented. Enrichment Analysis of DE genes in common dentex revealed 230 GO terms significantly enriched in testis (light blue colored in **Supplementary Excel Table SET2**) and 316 GO in ovary (light red colored in **Supplementary Excel Table SET2**). In gilthead seabream, 366 GO terms were significantly enriched in testes (light blue colored in **Supplementary Excel Table SET3**) and 298 GO in ovaries (light red colored in **Supplementary Excel Table SET3**).

Expression Profile of
Reproduction-Related and Key- Candidate
Genes in the Gonads

In order to identify genes that play an important role in reproduction we searched the annotated brain and gonad datasets for relevant GO categories and checked if they were differentially expressed. More particularly, DE genes related with broad biological processes terms, such as “reproduction” and “reproductive process” were searched in gonads of both sparids. A large number of “reproduction”-related and few “reproductive process”-related genes were found in the gonads of both species. In common dentex,

133 female-biased and 77 male-biased genes belonging to the GO term “reproduction,” and four female-biased and 14 male-biased genes belonging to “reproductive process” were detected. In gilthead seabream, 191 female-biased and 166 male-biased genes belonging to the GO term “reproduction,” and 19 female-biased and 20 male-biased genes belonging to “reproductive process” were detected (Supplementary Excel Tables SET4–12).

In an effort to identify the molecular machinery employed by each species to express the functional male or the female phenotype of the gonad, we aimed at analyzing and comparing the expression patterns of orthologous genes known for their involvement in the fish and the vertebrate sex determination/differentiation cascade. We selected 114 candidate genes in total. Combining data from the two newly sequenced species and the three sparids analyzed in Manousaki et al. (2014) and Tsakogiannis et al. (2018), we obtained a comparative view of expression of those genes in all five species (Table 3). Out of this list, 56 genes were present in the transcriptomes of all five species. Table 3 includes the 69 mostly discussed genes in other fish or vertebrates that were differentially expressed in at least one of the five species, and compares their expression pattern between them.
<table><caption>TABLE 3 | Expression pattern of key- gonadal genes implicated in sex determination/differentiation of the five spaid species.</caption><thead><tr><th rowspan="2">Gene name</th><th rowspan="2">Gene description</th><th colspan="5">Expression</th><th colspan="5">Gene ID</th></tr><tr><th>P. pagrus</th><th>P. erythrinus</th><th>D. dentex</th><th>D. puntazzo</th><th>S. aurata</th><th>P. pagrus</th><th>P. erythrinus</th><th>D. dentex</th><th>D. puntazzo</th><th>S. aurata</th></tr></thead><tbody><tr><td>Aldh1a2</td><td>Aldehyde dehydrogenase 1 family, member A2</td><td>NSD</td><td>NSD</td><td>M</td><td>M</td><td>M</td><td>TR49490|c0_g1</td><td>TR70446|c0_g1</td><td>c52437_g1</td><td>comp94934_c0</td><td>TRINITY_DN57044_c0_g1</td></tr><tr><td>Amh</td><td>Anti-Müllerian hormone or Müllerian-inhibiting substance</td><td>NSD</td><td>M</td><td>(-)</td><td>M</td><td>NSD</td><td>TR66874|c1_g1</td><td>TR91164|c5_g4</td><td>NA</td><td>comp107049_c0</td><td>TRINITY_DN166272_c0_g1</td></tr><tr><td>Amhr2</td><td>Anti-Müllerian hormone receptor 2</td><td>F</td><td>NSD</td><td>M</td><td>M</td><td>M</td><td>TR66380|c1_g1</td><td>TR90446|c1_g4</td><td>c70823_g1</td><td>comp111277_c1</td><td>TRINITY_DN58797_c0_g1</td></tr><tr><td>Ar</td><td>Androgen receptor / Androgen receptor alpha</td><td>F</td><td>NSD</td><td>M</td><td>M</td><td>F</td><td>TR60985|c0_g1</td><td>TR83707|c1_g2</td><td>c72407_g2</td><td>comp108409_c1</td><td>TRINITY_DN41747_c0_g2</td></tr><tr><td>Arb</td><td>Androgen receptor beta</td><td>F</td><td>F</td><td>NSD</td><td>NSD</td><td>NSD</td><td>TR60273|c1_g1</td><td>TR85026|c0_g1</td><td>c73222_g3</td><td>comp122243_c0</td><td>TRINITY_DN60348_c4_g1</td></tr><tr><td>Ctnnb1</td><td>Catenin (cadherin-associated protein), beta 1</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>TR46279|c0_g1</td><td>TR64263|c0_g2</td><td>c56056_g2</td><td>comp112822_c5</td><td>TRINITY_DN28847_c0_g1</td></tr><tr><td>Cyp11a2</td><td>Cytochrome P450, family 11, subfamily A, polypeptide 2 / 11alpha-hydroxylase</td><td>F</td><td>NSD</td><td>NSD</td><td>F</td><td>M</td><td>TR61267|c0_g1</td><td>TR88875|c1_g5</td><td>c71096_g1</td><td>comp111704_c1</td><td>TRINITY_DN61832_c0_g2</td></tr><tr><td>Cyp11b</td><td>Cytochrome P450, family 11, subfamily C, polypeptide 1 / 11beta-hydroxylase</td><td>NSD</td><td>M</td><td>M</td><td>M</td><td>M</td><td>TR47689|c0_g1</td><td>TR72664|c0_g1</td><td>c64871_g2</td><td>comp108630_c0</td><td>TRINITY_DN6958_c0_g1</td></tr><tr><td>Cyp19a1</td><td>Aromatase a (gonad isoform)</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>TR50491|c0_g1</td><td>TR72019|c0_g1</td><td>c56116_g1</td><td>comp117077_c0</td><td>TRINITY_DN50138_c0_g2</td></tr><tr><td>Cyp26a1</td><td>Cytochrome P450, family 26, subfamily a, polypeptide 1</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>TR49995|c0_g1</td><td>TR72448|c0_g1</td><td>c62822_g1</td><td>comp104920_c0</td><td>TRINITY_DN35779_c0_g1</td></tr><tr><td>Dax1</td><td>Dosage-sensitive sex reversal, adrenal hypoplasia critical region, on chromosome X, gene 1</td><td>F</td><td>F</td><td>NSD</td><td>M</td><td>F</td><td>TR64306|c1_g1</td><td>TR19608|c0_g1</td><td>c55395_g1</td><td>comp87968_c1</td><td>TRINITY_DN25170_c0_g1</td></tr><tr><td>Ddx11</td><td>DEAD/H box helicase 11</td><td>M</td><td>M</td><td>NSD</td><td>NSD</td><td>NSD</td><td>TR66969|c0_g1</td><td>TR90413|c2_g1</td><td>c75063_g2</td><td>comp114371_c9</td><td>TRINITY_DN59935_c0_g1</td></tr><tr><td>Ddx4</td><td>Dead Box 4 Vasa Homolog</td><td>NSD</td><td>NSD</td><td>M</td><td>NSD</td><td>M</td><td>TR52887|c0_g1</td><td>TR69725|c0_g1</td><td>c59427_g1</td><td>comp116347_c5</td><td>TRINITY_DN64824_c0_g2</td></tr></tbody><tfoot><tr><td colspan="12" style="text-align: right;">(Continued)</td></tr></tfoot></table>
<table><caption>TABLE 3 | Continued</caption><thead><tr><th rowspan="2">Gene name</th><th rowspan="2">Gene description</th><th colspan="5">Expression</th><th colspan="5">Gene ID</th></tr><tr><th>P. pagrus</th><th>P. erythrinus</th><th>D. dentex</th><th>D. puntazzo</th><th>S. aurata</th><th>P. pagrus</th><th>P. erythrinus</th><th>D. dentex</th><th>D. puntazzo</th><th>S. aurata</th></tr></thead><tbody><tr><td>Dhh</td><td>Desert hedgehog</td><td>NSD</td><td>F</td><td>F</td><td>F</td><td>F</td><td>TR59771|c0_g1</td><td>TR84448|c0_g1</td><td>c62990_g1</td><td>comp116302_c0</td><td>TRINITY_DN25379_c0_g1</td></tr><tr><td>Dmrt1</td><td>Doublesex- and mab-3-related transcription factor 1</td><td>M</td><td>M</td><td>M</td><td>M</td><td>M</td><td>TR67635|c0_g1</td><td>TR70670|c0_g1</td><td>c71436_g2</td><td>comp110196_c1</td><td>TRINITY_DN51794_c0_g2</td></tr><tr><td>Dmrt2a</td><td>Doublesex- and mab-3-related transcription factor 2(a)</td><td>F</td><td>F</td><td>F</td><td>F</td><td>(-)</td><td>TR59371|c2_g1</td><td>TR78139|c1_g1</td><td>c64698_g2</td><td>comp112381_c0</td><td>NA</td></tr><tr><td>Dmrt3a</td><td>Double sex and mab-3 related transcription factor 3(a)</td><td>(-)</td><td>M</td><td>M</td><td>M</td><td>(-)</td><td>NA</td><td>TR67762|c1_g1</td><td>c65860_g1</td><td>comp88038_c1</td><td>NA</td></tr><tr><td>Dmrt1</td><td>DNA methyltransferase 1</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>TR62712|c0_g1</td><td>TR83768|c0_g1</td><td>c61611_g1</td><td>comp115925_c0</td><td>TRINITY_DN51749_c0_g1</td></tr><tr><td>Dmrt3aa</td><td>DNA methyltransferase 3aa</td><td>M</td><td>M</td><td>M</td><td>M</td><td>M</td><td>TR64267|c0_g1</td><td>TR73384|c7_g1</td><td>c70794_g1</td><td>comp120818_c3</td><td>TRINITY_DN58191_c0_g1</td></tr><tr><td>Dmrt3ab</td><td>DNA methyltransferase 3ab</td><td>F</td><td>F</td><td>M</td><td>M</td><td>NSD</td><td>TR53015|c0_g1</td><td>TR78244|c1_g1</td><td>c71723_g5</td><td>comp120102_c11</td><td>TRINITY_DN57426_c3_g1</td></tr><tr><td>Ep300a</td><td>Histone acetyltransferase—E1A binding protein 300a</td><td>M</td><td>M</td><td>NSD</td><td>NSD</td><td>M</td><td>TR60830|c1_g1</td><td>TR84051|c2_g6</td><td>c63663_g10</td><td>comp114282_c1</td><td>TRINITY_DN61531_c11_g1</td></tr><tr><td>Ep300b</td><td>Histone acetyltransferase—E1A binding protein 300b</td><td>NSD</td><td>NSD</td><td>(-)</td><td>F</td><td>F</td><td>TR62718|c3_g1</td><td>TR88233|c0_g14</td><td>c63663_g1</td><td>NA</td><td>TRINITY_DN61531_c7_g1</td></tr><tr><td>Er</td><td>Estrogen receptor</td><td>M</td><td>NSD</td><td>M</td><td>NSD</td><td>M</td><td>TR68522|c4_g2</td><td>TR84515|c0_g1</td><td>c71439_g2</td><td>comp117911_c10</td><td>TRINITY_DN43934_c0_g1</td></tr><tr><td>Esr2a</td><td>Type II estrogen receptor</td><td>NSD</td><td>NSD</td><td>M</td><td>NSD</td><td>F</td><td>TR65583|c1_g1</td><td>TR85716|c0_g1</td><td>c72247_g1</td><td>comp112065_c1</td><td>TRINITY_DN59537_c6_g2</td></tr><tr><td>Esr2b</td><td>Estrogen receptor 2b</td><td>(-)</td><td>(-)</td><td>M</td><td>(-)</td><td>(-)</td><td>NA</td><td>NA</td><td>NA</td><td>comp108758_c2</td><td>NA</td></tr><tr><td>Esra</td><td>Estrogen-related receptor alpha</td><td>M</td><td>M</td><td>M</td><td>M</td><td>F</td><td>TR60672|c1_g2</td><td>TR86271|c1_g1</td><td>c68941_g2</td><td>comp112646_c1</td><td>TRINITY_DN60587_c0_g1</td></tr><tr><td>Esrb</td><td>Estrogen-related receptor beta</td><td>F</td><td>F</td><td>F</td><td>NSD</td><td>F</td><td>TR66216|c1_g1</td><td>TR89871|c1_g3</td><td>c71889_g1</td><td>comp117190_c7</td><td>TRINITY_DN55451_c1_g1</td></tr></tbody></table>

(Continued)
<table><caption>TABLE 3 | Continued</caption><thead><tr><th rowspan="2">Gene name</th><th rowspan="2">Gene description</th><th colspan="5">Expression</th><th colspan="5">Gene ID</th></tr><tr><th>P. pagrus</th><th>P. erythrinus</th><th>D. dentex</th><th>D. puntazzo</th><th>S. aurata</th><th>P. pagrus</th><th>P. erythrinus</th><th>D. dentex</th><th>D. puntazzo</th><th>S. aurata</th></tr></thead><tbody><tr><td>Esrb2</td><td>Estrogen receptor 2b</td><td>NSD</td><td>F</td><td>NSD</td><td>NSD</td><td>M</td><td>TR52662|c0_g1</td><td>TR77491|c0_g1</td><td>c61415_g2</td><td>comp108758_c0</td><td>TRINITY_DN53665_c0_g1</td></tr><tr><td>Fg/20a</td><td>Fibroblast growth factor 20a</td><td>NSD</td><td>NSD</td><td>NSD</td><td>M</td><td>(-)</td><td>TR56801|c3_g1</td><td>TR72320|c1_g1</td><td>c37432_g2</td><td>comp106110_c1</td><td>NA</td></tr><tr><td>Fg/20b</td><td>Fibroblast growth factor 20b</td><td>F</td><td>NSD</td><td>M</td><td>M</td><td>F</td><td>TR56801|c1_g1</td><td>TR73946|c0_g1</td><td>c62141_g1</td><td>comp85925_c1</td><td>TRINITY_DN108566_c0_g1</td></tr><tr><td>Fg/9</td><td>Fibroblast growth factor 9</td><td>M</td><td>NSD</td><td>NSD</td><td>(-)</td><td>(-)</td><td>TR66253|c2_g4</td><td>TR91211|c1_g5</td><td>c72708_g3</td><td>NA</td><td>NA</td></tr><tr><td>Figla</td><td>Factor in the germline alpha</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>TR49469|c0_g1</td><td>TR54099|c0_g1</td><td>c52499_g1</td><td>comp119797_c11</td><td>TRINITY_DN14620_c1_g1</td></tr><tr><td>Foxl2</td><td>Forkhead box L2</td><td>F</td><td>F</td><td>F</td><td>NSD</td><td>F</td><td>TR43506|c0_g1</td><td>TR59502|c0_g1</td><td>c59874_g1</td><td>comp95964_c0</td><td>TRINITY_DN52797_c0_g1</td></tr><tr><td>Fsta</td><td>Follistatin</td><td>F</td><td>F</td><td>M</td><td>M</td><td>(-)</td><td>TR48089|c0_g1</td><td>TR66770|c0_g3</td><td>c56990_g1</td><td>comp112642_c0</td><td>NA</td></tr><tr><td>Fstl3</td><td>Follistatin-like 3</td><td>NSD</td><td>NSD</td><td>M</td><td>M</td><td>M</td><td>TR64396|c0_g3</td><td>TR80964|c2_g1</td><td>c67341_g1</td><td>comp103373_c0</td><td>TRINITY_DN44642_c0_g1</td></tr><tr><td>Fstl4</td><td>Follistatin-like 4</td><td>NSD</td><td>NSD</td><td>F</td><td>M</td><td>NSD</td><td>TR66843|c3_g1</td><td>TR80451|c1_g1</td><td>c40110_g1</td><td>comp118767_c9</td><td>TRINITY_DN59298_c0_g1</td></tr><tr><td>Fstl5</td><td>Follistatin-like 5</td><td>M</td><td>M</td><td>M</td><td>NSD</td><td>NSD</td><td>TR60004|c0_g1</td><td>TR77519|c0_g5</td><td>c72310_g1</td><td>NA</td><td>TRINITY_DN64086_c1_g2</td></tr><tr><td>FTZ-F1</td><td>Steroidogenic factor-1/fushi tarazu factor-1</td><td>F</td><td>NSD</td><td>M</td><td>M</td><td>M</td><td>TR56092|c1_g1</td><td>TR80163|c1_g2</td><td>c67456_g3</td><td>comp120045_c3</td><td>TRINITY_DN50710_c0_g2</td></tr><tr><td>Gata-4</td><td>Gata-binding protein 4</td><td>NSD</td><td>F</td><td>NSD</td><td>M</td><td>M</td><td>TR63369|c1_g1</td><td>TR74516|c2_g1</td><td>c68365_g1</td><td>comp114090_c1</td><td>TRINITY_DN49474_c0_g2</td></tr><tr><td>Gdf9</td><td>Growth and differentiation factor 9</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>TR52139|c0_g1</td><td>TR70663|c0_g1</td><td>c59003_g1</td><td>comp69471_c0</td><td>TRINITY_DN8454_c0_g1</td></tr><tr><td>Gsdf</td><td>Gonadal soma derived factor</td><td>F</td><td>NSD</td><td>M</td><td>(-)</td><td>M</td><td>TR61049|c0_g1</td><td>TR72782|c0_g1</td><td>c67953_g1</td><td>NA</td><td>TRINITY_DN56147_c0_g1</td></tr><tr><td>Hdac10</td><td>Histone deacetylase 10</td><td>F</td><td>NSD</td><td>F</td><td>F</td><td>NSD</td><td>TR67490|c0_g1</td><td>TR88916|c0_g6</td><td>c74106_g2</td><td>comp104913_c0</td><td>TRINITY_DN55862_c0_g1</td></tr><tr><td>Hdac11</td><td>Histone deacetylase 11</td><td>F</td><td>NSD</td><td>M</td><td>M</td><td>F</td><td>TR66619|c2_g2</td><td>TR84940|c1_g2</td><td>c74380_g7</td><td>comp120490_c0</td><td>TRINITY_DN63334_c0_g1</td></tr><tr><td>Hdac2</td><td>Histone deacetylase 2</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>TR56322|c1_g2</td><td>TR66800|c1_g1</td><td>c70637_g5</td><td>comp104608_c0</td><td>TRINITY_DN47239_c0_g1</td></tr><tr><td>Hdac7</td><td>Histone deacetylase 7</td><td>F</td><td>F</td><td>NSD</td><td>M</td><td>F</td><td>TR57004|c0_g1</td><td>TR84823|c3_g1</td><td>c67700_g1</td><td>comp118197_c1</td><td>TRINITY_DN56857_c0_g1</td></tr></tbody></table>

(Continued)
<table><caption>TABLE 3 | Continued</caption><thead><tr><th rowspan="2">Gene name</th><th rowspan="2">Gene description</th><th colspan="5">Expression</th><th colspan="5">Gene ID</th></tr><tr><th>P. pagrus</th><th>P. erythrinus</th><th>D. dentex</th><th>D. puntazzo</th><th>S. aurata</th><th>P. pagrus</th><th>P. erythrinus</th><th>D. dentex</th><th>D. puntazzo</th><th>S. aurata</th></tr></thead><tbody><tr><td>Hdac8</td><td>Histone deacetylase 8</td><td>M</td><td>M</td><td>M</td><td>M</td><td>M</td><td>TR62362|c1_g1</td><td>TR8678|c2_g5</td><td>c67655_g1</td><td>comp111646_c0</td><td>TRINITY_DN59663_c0_g1</td></tr><tr><td>Lhb</td><td>Luteinizing hormone, beta polypeptide</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>TR40794|c1_g1</td><td>TR66625|c0_g1</td><td>c21134_g1</td><td>comp104837_c0</td><td>TRINITY_DN49327_c0_g1</td></tr><tr><td>Lhr</td><td>Luteinizing hormone receptor</td><td>F</td><td>NSD</td><td>M</td><td>M</td><td>M</td><td>TR61712|c0_g1</td><td>TR89660|c1_g2</td><td>c69306_g1</td><td>comp101199_c1</td><td>TRINITY_DN55453_c0_g1</td></tr><tr><td>Nr3c1</td><td>Glucocorticoid receptor</td><td>M</td><td>M</td><td>M</td><td>M</td><td>NSD</td><td>TR63170|c1_g1</td><td>TR84453|c3_g1</td><td>c69943_g1</td><td>comp110357_c2</td><td>TRINITY_DN30611_c0_g3</td></tr><tr><td>Nr3c2</td><td>Mineralocorticoid receptor</td><td>F</td><td>NSD</td><td>M</td><td>M</td><td>F</td><td>TR58538|c2_g1</td><td>TR84031|c1_g1</td><td>c70382_g1</td><td>comp109470_c2</td><td>TRINITY_DN45316_c0_g1</td></tr><tr><td>Pdgtfa</td><td>Platelet-derived growth factor alpha</td><td>M</td><td>NSD</td><td>NSD</td><td>M</td><td>M</td><td>TR66514|c0_g1</td><td>TR87951|c1_g6</td><td>c73083_g1</td><td>comp113241_c2</td><td>TRINITY_DN57763_c0_g1</td></tr><tr><td>Pdgtfb</td><td>Platelet-derived growth factor beta</td><td>F</td><td>F</td><td>NSD</td><td>M</td><td>(-)</td><td>TR64531|c0_g1</td><td>TR75590|c0_g1</td><td>c66673_g2</td><td>comp107178_c1</td><td>NA</td></tr><tr><td>Pdgtfa</td><td>Platelet-derived growth factor beta a</td><td>F</td><td>NSD</td><td>M</td><td>M</td><td>NSD</td><td>TR58625|c0_g1</td><td>TR86429|c1_g1</td><td>c70504_g1</td><td>comp117159_c4</td><td>TRINITY_DN61779_c0_g1</td></tr><tr><td>Pdgtfb</td><td>Platelet-derived growth factor beta b</td><td>NSD</td><td>NSD</td><td>M</td><td>M</td><td>M</td><td>TR65607|c0_g1</td><td>TR87425|c0_g1</td><td>c67490_g1</td><td>comp113875_c0</td><td>TRINITY_DN40994_c0_g2</td></tr><tr><td>Pdgtfr</td><td>Platelet-derived growth factor receptor, alpha</td><td>NSD</td><td>NSD</td><td>M</td><td>M</td><td>M</td><td>TR66258|c1_g1</td><td>TR88623|c1_g1</td><td>c71284_g1</td><td>comp116299_c0</td><td>TRINITY_DN59458_c0_g1</td></tr><tr><td>Pdgtfr1</td><td>Platelet-derived growth factor receptor, beta 1</td><td>NSD</td><td>NSD</td><td>M</td><td>M</td><td>M</td><td>TR64788|c0_g1</td><td>TR80831|c3_g2</td><td>c71970_g3</td><td>comp104719_c1</td><td>TRINITY_DN52984_c0_g1</td></tr><tr><td>Pdgtfr2</td><td>Platelet-derived growth factor receptor, beta 2</td><td>F</td><td>F</td><td>F</td><td>NSD</td><td>NSD</td><td>TR55622|c2_g1</td><td>TR83411|c3_g1</td><td>c73150_g1</td><td>comp113376_c0</td><td>TRINITY_DN57277_c0_g1</td></tr><tr><td>Raraa</td><td>Retinoid acid receptor alpha a</td><td>M</td><td>M</td><td>M</td><td>M</td><td>M</td><td>TR36645|c1_g1</td><td>TR82123|c4_g2</td><td>c64412_g3</td><td>comp113587_c1</td><td>TRINITY_DN36265_c0_g1</td></tr><tr><td>Rarab</td><td>Retinoid acid receptor alpha b</td><td>F</td><td>NSD</td><td>M</td><td>NSD</td><td>NSD</td><td>TR63371|c2_g1</td><td>TR84655|c1_g1</td><td>c68281_g8</td><td>comp120127_c2</td><td>TRINITY_DN45351_c0_g2</td></tr><tr><td>Rarb</td><td>Retinoid acid receptor beta</td><td>M</td><td>M</td><td>M</td><td>M</td><td>M</td><td>TR65946|c1_g2</td><td>TR81734|c1_g3</td><td>c72322_g2</td><td>comp119441_c0</td><td>TRINITY_DN63608_c1_g1</td></tr></tbody></table>

(Continued)
<table><caption>TABLE 3 | Continued</caption><thead><tr><th rowspan="2">Gene name</th><th rowspan="2">Gene description</th><th colspan="5">Expression</th><th colspan="5">Gene ID</th></tr><tr><th>P. pagrus</th><th>P. erythrinus</th><th>D. dentex</th><th>D. puntazzo</th><th>S. aurata</th><th>P. pagrus</th><th>P. erythrinus</th><th>D. dentex</th><th>D. puntazzo</th><th>S. aurata</th></tr></thead><tbody><tr><td>Rspo1</td><td>R-spondin-1</td><td>NSD</td><td>NSD</td><td>M</td><td>(-)</td><td>(-)</td><td>TR44782|c0_g2</td><td>TR43218|c0_g2</td><td>c50451_g1</td><td>NA</td><td>NA</td></tr><tr><td>Sox3</td><td>SRY-related HMG box 3</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>TR64334|c1_g1</td><td>TR69218|c2_g2</td><td>c60997_g2</td><td>comp87875_c1</td><td>TRINITY_DN61608_c6_g1</td></tr><tr><td>Star-like</td><td>Steroidogenic acute regulatory protein</td><td>NSD</td><td>F</td><td>M</td><td>M</td><td>NSD</td><td>TR53314|c0_g1</td><td>TR52301|c0_g1</td><td>c62942_g2</td><td>comp76502_c0</td><td>TRINITY_DN51715_c0_g2</td></tr><tr><td>Stra6</td><td>Stimulated by retinoic acid gene 6</td><td>F</td><td>F</td><td>F</td><td>F</td><td>F</td><td>TR59700|c3_g1</td><td>TR86270|c0_g1</td><td>c68491_g1</td><td>comp105129_c1</td><td>TRINITY_DN43771_c0_g1</td></tr><tr><td>Tdrd1</td><td>Tudor domain containing 1</td><td>NSD</td><td>M</td><td>M</td><td>M</td><td>M</td><td>TR68477|c1_g4</td><td>TR91962|c2_g5</td><td>c70984_g2</td><td>comp113472_c0</td><td>TRINITY_DN59161_c0_g1</td></tr><tr><td>Tdrd7</td><td>Tudor domain containing 7</td><td>M</td><td>M</td><td>M</td><td>M</td><td>M</td><td>TR63911|c0_g1</td><td>TR88488|c1_g2</td><td>c73678_g1</td><td>comp115842_c3</td><td>TRINITY_DN61324_c1_g1</td></tr><tr><td>Wnt4a</td><td>Wingless-type MMTV integration site family, member 4a</td><td>F</td><td>F</td><td>F</td><td>(-)</td><td>(-)</td><td>TR40636|c0_g1</td><td>TR67106|c0_g2</td><td>c15210_g2</td><td>NA</td><td>NA</td></tr><tr><td>WT7</td><td>Wilms tumor protein 1a</td><td>F</td><td>F</td><td>NSD</td><td>M</td><td>M</td><td>TR67279|c4_g5</td><td>TR89074|c0_g1</td><td>c64908_g1</td><td>comp117315_c1</td><td>TRINITY_DN59519_c2_g1</td></tr><tr><td>WT7b</td><td>Wilms tumor protein 1b</td><td>F</td><td>F</td><td>NSD</td><td>(-)</td><td>M</td><td>TR57899|c0_g1</td><td>TR89074|c1_g1</td><td>c67340_g1</td><td>NA</td><td>TRINITY_DN57126_c1_g1</td></tr></tbody></table>

M, male-biased; F, female-biased; NSD, not significantly different.
## Comparative Analysis of Sex-Biased Gene Expression in the Gonads of All Five Sparid Species

Apart from the 114 genes that were identified based on literature and public datasets, we looked for sex-biased differences using a comparative whole-transcriptome approach, based on the orthology assignment of the annotated transcriptomes of the gonads of all the species (Figure 4) (homologous genes between the five studied species in Supplementary Excel Table SET13). We compared the expression data of the annotated gonadal orthologous genes for the five species. Those multi-group comparisons helped us to identify common (in two, three, four, and all five species) and unique (species-specific) male- and female-biased genes in gonads. To gain greater insight into the genes that consistently play a role in gonadal identity we focused on genes that were differentially expressed in all five species (common/shared genes). The analysis revealed 1,418 annotated genes expressed in the testes (Figure 4A, intersection of all five species; Supplementary Excel Table SET14) and 1,665 annotated genes expressed in the ovaries (Figure 4B, intersection of all five species; Supplementary Excel Table SET15) of all five studied sparids. Moreover, this analysis allowed us to make comparisons of the testis- or ovary- specific genes between different reproductive modes. As such, one could focus on the protogynous-specific genes, with a shared over-expression in only the two protogynous species vs. all the other species. The analysis revealed 349 protogynous-specific orthologous genes overexpressed in testes (Supplementary Excel Table SET16) and 578 protogynous-specific orthologs overexpressed in ovaries (Supplementary Excel Table SET17). Considering that the sharpsnout seabream is a rudimentary or late gonochoristic species, the comparative analysis depicted 1,098 gonochoristic-specific male-biased orhologs (Supplementary Excel Table SET18), and 238 gonochoristic-specific female-biased orthologs (Supplementary Excel Table SET19). In an effort to gain a more profound insight into the protandrous-specific orthologs, the analysis produced 773 overexpressed genes in the testes (Supplementary Excel Table SET20) and 1,116 overexpressed genes in the ovaries of only the protandrous gilthead seabream (Supplementary Excel Table SET21). The maximum number of male-biased shared orthologs (1,098 genes) was found in the gonochoristic common dentex and rudimentary hermaphrodite sharpsnout seabream, when the maximum number of female-biased shared orthologs (578 genes) was found in the protogynous common pandora and red porgy.

## DISCUSSION

The present study reports new transcriptome assemblies for two species of the Sparidae family, and compares the sex-specific gene-expression patterns in gonads and brains in four hermaphrodites and one gonochorist belonging to this family. This is the first comparative transcriptomic analysis of five phylogenetically close fish species displaying alternative reproductive modes. Focusing on genes that are expressed

differentially between sexes, we pinpoint a common network that forms and/or maintains the functional male and female phenotype.

## Slight Differences Between Sexes in Brain. A “Sexually” Un-biased Tissue?

The analysis of the two newly sequence sparid species in this study confirms the insignificant sex-biased expression differences in the brain, reported previously in Manousaki et al. (2014) and Tsakogiannis et al. (2018). Generally, the DE genes in the brain are limited both in number and in expression magnitude. A limited number of DE genes between sexes were found in gilthead seabream brain, similar to those found in the rudimentary hermaphrodite sharpsnout seabream, whereas a total lack of sex-biased expression was observed in the gonochoristic common dentex. Conversely, the greatest differences in expression between males and females were observed in the brain of the two protogynous species. Regarding these protogynous species, it has been observed that the red porgy brain samples are more tightly clustered than those of common pandora that demonstrate a more relaxed clustering (Tsakogiannis et al., 2018). Moreover, the direction of expression was female-biased in the protogynous common pandora, the rudimentary hermaphrodite sharpsnout seabream and the protandrous gilthead seabream. These results are in line with the expression pattern observed in the protandrous clownfish Amphiprion bicinctus brain (Casas et al., 2016). In contrast, a male-biased pattern was observed here in the red porgy, as seen also in the protogynous bluehead wrasse Thalassoma bifasciatum (Liu et al., 2015) and four cichlid species (Böhne et al., 2014).

An interesting pattern arises when comparing the protogynous (highest sex-biased expression in brain), the protandrous and rudimentary hermaphrodite species (very low sex-biased expression in brain), and the gonochoristic species (no sex-biased expression in brain). However, the unbiased brain of the common dentex is not a common pattern in gonochoristic species, as this has not been previously reported in four gonochoristic East African cichlid species showing a low, but not zero sex-biased expression in their brain (Böhne et al., 2014). Generally, we found no consensus pattern among fish from the present and previous studies, indicating species-specific profiles in the brain. This has been previously observed and is linked to the great variation of sex-dependent behavioral and physiological traits in fishes (Maehiro et al., 2014).

A closer look at 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. Some noteworthy genes that are putatively associated with male-like behavior were the vasotocin-neurophysin VT 1-like and isotocin-neurophysin IT 1-like, found among male-biased genes in the brain of gilthead seabream. We further searched in the transcriptomes of common dentex, gilthead seabream and sharpsnout seabream for the early growth response factor 1 (Egr1), the early growth response factor 4 (Egr4), as well as the proto-oncogene c-Fos (Fos), that were observed in the brain of the two protogynous species, the common pandora
A

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

B

<figure><img src="image_7.png" /><figcaption>**FIGURE 4** | Multigroup comparisons of overexpressed genes in the gonads of all five sparids. The graphic shows the numbers of shared genes between species within the two groups of DE genes (A) testis overexpressed or male-biased and (B) ovary overexpressed or female-biased.</figcaption></figure>

and the red porgy (Tsakogiannis et al., 2018). Interestingly, in gilthead seabream the *Fos* gene was female-biased and thus it might be indicative of the different behavior-related reproductive strategies between protandrous and protogynous species. These genes were also present in the brain transcriptomes of the rest of the species, but were unbiased. The *Egr1* and *c-Fos* genes, both characterized as immediate-early genes, are related to changes in brain activity and genomic responses to social cues, behavior and mating information (Desjardins et al., 2010).

## Common Patterns of Gonadal Gene Expression in All Five Sparids

The highest degree of sex-biased gene expression is observed in the gonads (Mank et al., 2008). Despite the flexibility of gonadal tissue in sex-changers, the degree of difference between testes and ovaries is the same between the hermaphrodites and the gonochorist. Testes and ovaries are derived from the same sexually undifferentiated gonadal primordium, but once the sex has been established, even in hermaphrodites they represent two morphologically and functionally different organs. The transcriptomic profiles of the functional male and the female gonad differ in terms of the number of DE genes, as well as the expression magnitude (i.e., fold-change differences).

The male-biased expression tendency in gonads, applies in all studied sparids and it seems to be a common pattern in fish (Tao et al., 2013; Böhne et al., 2014; Manousaki et al., 2014; Liu et al., 2015; Casas et al., 2016; Machado et al., 2016; Xu et al., 2016; Tsakogiannis et al., 2018). Important signaling pathways such as WNT (wingless-type MMTV integration site family), NOTCH, TGFβ (transforming growth factor beta), nuclear receptors, for example the retinoic acid receptor (RAR), and RTKs (receptor tyrosine kinases) such as FGFR (fibroblast growth factor receptor), EGFR (epidermal growth factor receptor), and IR (insulin receptor) were all found to be among the significantly overrepresented terms of common dentex and gilthead seabream

as well as in the rest of the species (Manousaki et al., 2014; Tsakogiannis et al., 2018).

The lists of enriched GOs were searched for biological processes related to the context of gonadal differentiation. Searching for terms related to hormone regulation we discovered several relevant terms over-represented in common dentex and gilthead seabream gonads such as: “positive regulation of growth hormone receptor signaling pathway,” “positive regulation of hormone metabolic process,” “growth hormone secretion” in testes of gilthead seabream and “positive regulation of peptide hormone secretion” in the testes of common dentex. Furthermore, androgen and estrogen related terms such as “positive regulation of androgen secretion,” “androgen receptor signaling pathway,” “androgen metabolic process,” “estrogen biosynthetic process” and “response to estrogen” were found to be enriched in the gonads of both common dentex and gilthead seabream. Among the female-specific GOs in both common dentex and gilthead seabream were the “positive regulation of ovarian follicle development,” while testis enriched terms in gilthead seabream were “sperm motility,” “spermatid development.” Interestingly, in common dentex female gonads the “hermaphrodite genitalia development” term was found.

Findings such as “response to retinoic acid,” “regulation of retinoic acid biosynthetic process,” “positive regulation of retinoic acid receptor signaling pathway” or “retinoic acid receptor signaling pathway” among both species gonads’ enriched terms, indicate the role of the retinoic acid (RA) metabolism in gonadal development.

Furthermore, “negative regulation of canonical Wnt signaling pathway” and “regulation of Wnt signaling pathway, calcium modulating pathway” terms were enriched in ovaries of common dentex and gilthead seabream, respectively. Additionally, the term “cellular response to transforming growth factor beta stimulus” and “positive regulation of transforming growth factor beta receptor signaling pathway” was enriched in male gonads
of common dentex and gilthead seabream, respectively. Known and described mostly in mammals (Windley and Wilhelm, 2015), these pathways apparently play an essential role in the context of sex differentiation and gonadal differentiation and development in the studied teleost species.

Lastly, in both common dentex and gilthead seabream, many epigenetic pathways are over-represented and seem to co-regulate reproductive processes. The enriched GO terms cover two major epigenetic mechanisms for gene expression regulation such as DNA methylation and histone modifications (acetylation, deacetylation).

## The Gonadal Transcriptomic Profile of Different Sex- Modes

A transcriptome-wide comparative analysis provided the opportunity to examine sex-biased expression profiles in five species displaying alternative reproductive modes and identify several interesting genes with sex mode-specific expression patterns (Figure 5). Focusing on the gonochoristic-specific (common dentex and sharsnout seabream) male-biased genes, several growth factor genes like *Fgf20*, *Fgfr2*, and *Fgfr1l* were found. Testis development requires the repression of *Wnt4* by *Fgf* signaling (Jameson et al., 2012). As we already know from mammals, *Fgf9* has an important function in male development, creating a positive feedback cycle with *Sox9*, inhibiting the WNT4 pathway in testis. The *Fgf9* gene has not yet been detected in teleosts and seems to be replaced by *Fgf20*, one of the genes involved in male sexual development also in *Latimeria* (Forconi et al., 2013). Among the gonochoristic-specific male-biased genes was also the *Foxc1a* that is required for primordial germ cell migration and antral follicle development in mice (Mattiske et al., 2006).

The analysis also revealed genes encoding for transcription factors such as *Tgfb1*, *Tgfb3*, *Tcf7*, and *Mafb*, where the former is considered to be expressed by the Sertoli cells and may be important for spermatogenesis (Memon et al., 2008), while, regarding the latter, studies in mice show that it is indispensable for the fetal testis morphogenesis and the maintenance of spermatogenesis (Shawki et al., 2018). The female-biased gonochoristic-specific genes included, among others, the TATA binding proteins *Taf8* and *Taf15*, the cytochoromic enzymes *Cyp20a* and *Spd1l* that is also a female-biased gene in the guppy, *Poecilia reticulata* (Sharma et al., 2014). Additionally, *Rec8* was found among the female-biased genes. This gene is present in all sequenced genomes of teleosts, including tilapia, but it is yet to be determined whether RA regulates germ cells meiotic initiation via *Rec8* in teleosts (Feng et al., 2015).

The protandrous-specific (gilthead seabream) male-biased genes included several transcription factors like *Supt20h*, *Maff*, *Brf2*, *Taf4*, *Taf7*. *Taf7* is one of genes shown to contribute to male fertility in mouse genetic studies and is expressed in spermatogenic cells (Cheng et al., 2007). A protandrous-specific gene expressed in testis were also the *Wnt8*, the testis expressed sequence 2, *Tex2*, and the *Cyp17a*. Recently, studies on mice detected *Cyp17* in germ cells: spermatogonia, pachytene spermatocytes, spermatids and sperm, and the loss of this enzyme

is associated with sperm abnormalities and infertility. This suggests that, in addition to steroidogenesis, *Cyp17* may have a role in sperm structure and function (Hinfray et al., 2013).

On the other hand, the protandrous-specific female-biased genes included the estrogen receptors *Esrrg* and *Esr2*, the fibroblast growth factors *Fgf12* and *Frs2*, the forkhead box proteins *Foxp2* and *Foxn4* as well as the transcription factors *Etv6*, *Sp8b*, and *Sox8*. Interestingly, *Sox8*, a male determining factor, expressed uniquely in the ovaries of the protandrous gilthead seabream.

Among the protogynous-specific (common pandora and red porgy) male-biased genes were the testis expressed sequence 10, *Tex10*, the transcription factor *Ebf3*, the sperm associated protein *Spaca4* and the transforming growth factor b association protein *Tgfbrap1*. The protogynous-specific female-biased genes included the platelet-derived growth factor *Pdfgfb*, the platelet activating factor receptor *Ptafr* and the platelet receptor *C10orf54*. Of interest is that several transcription growth factors like *Tcf7*, *Hes5*, *Mafb* as well as *Sox7* were all among the female-biased genes in the two protogynous species while they were all male-biased in the gonochoristic common dentex and the late—gonochoristic species, sharpsnout seabream.

## Genes With a Conserved Role to Identify the Male and Female Character of the Gonad in Sparidae

The analysis provided us with more than a thousand distinct annotated orthologous genes in all five studied species, which might be regarded as potential candidates of gonadal sex identity in this group of teleosts and with a conserved role in identifying the male and female functional character of the gonad in Sparidae.

To gain more insight into the genes that have a conserved role in preserving the gonadal identity of the established sex at this group of teleosts we focused on a particular group, the genes that were DE in all five species (common/shared genes). Among the 1,418 annotated genes expressed in the testes, several genes coding for transcription factors like *Jund*, *E2f1*, *E2f2*, *Ebf3*, *E4f1*, and *Rfx3* and transcription activators like *Smarca4*, were found. There were also some testis associated proteins *Spag6*, *Spag8*, *Spag17* and spermatogenesis associated proteins like *Spata4* and *Spata22*, *Spef2* (the sperm flagellar protein 2), as well as some apoptotic factors like *Tp53bp2*, *Ppp1r13b*, and the *Aifm3* (apoptosis inducing factor). Among the 1,665 annotated genes expressed in the ovaries, there were the transcription factors *Btf3*, *Btf3l4*, *E2f4*, *Hes7*, *Brf1*, and *Ovol1* that are characterized by a specific and sustained high expression level in females compared to males, a finding that is also observed in the gonads of rainbow trout (Vizziano et al., 2007), including also the well-known *Sox* factors *Sox3* and *Sox11*. Additionally, the zona pellucida genes *Zp3* and *Zp4* were commonly over-expressed in ovaries of all five species.

Studies on other fish and vertebrates have identified several candidate genes with proven or presumed relationship in sex determination and differentiation. We compiled a list with some of the most-studied sex-biased genes and compared their
<figure><img src="image_8.png" /><figcaption>**FIGURE 5** | Comparative view of the DE genes overexpressed in testes (blue) and ovaries (pink) specific only to gonochoristic species (Common dentex and Sharpsnout seabream), to protogynous (Common pandora and Red porgy), to protandrous (Gilthead seabream) and shared among all five species.</figcaption></figure>

expression patterns in the gonads of the studied sparids. The majority of the studied genes were present in the transcriptome of all five Sparidae species. Interestingly, there were some that showed a consensus bias in expression toward one sex. This might prove not only their association with the sex differentiation mechanism in fish but also their conserved roles in promoting either the female or the male character of the gonad in these species or, reductively, in this group of teleosts. The expression profiling brought to light well-known genes whose expression pattern jointly identifies the female (*Cyp19a1*, *Sox3*, *Figla*, *Gdf9*, *Cyp26a*, *Ctnnb1*, *Dnmt1*, *Stra6*) and male phenotype (*Dmrt1*, *Sox9*, *Dnmt3aa*, *Rarb*, *Raraa*, *Hdac8*, *Tdrd7*) of the gonad between these sparids.

Starting with the female phenotype of the gonad, that has long been assumed as the default state at least in mammals (Herpin and Schartl, 2011b), gonadal aromatase (*Cyp19a1a*) gene is a protagonist in the gene network controlling the gonadal fate. Here, a strong female-biased expression of gonadal aromatase gene (*Cyp19a1a*) was observed in all studied sparids. In fish, *Cyp19a1a* has a fundamental role in ovarian differentiation, both in gonochorists (Kitano et al., 1999; Uchida et al., 2004; Karube et al., 2007; Böhne et al., 2013), but also in hermaphrodites

where its connection with sex-reversal is underlined (Kobayashi et al., 2010; Wu et al., 2010; Zhou and Gui, 2010; Shao et al., 2014; Baroiller and D'Cotta, 2016; Casas et al., 2016). The *Cyp19a1a* together with estrogens seem to control both ovarian and testicular differentiation in an adverse way. Over-expression of *Cyp19a1a* prompts and maintains female development, while its reduced expression allows for male development to occur (Guiguen et al., 2010). The *Cyp19a1a* that seems to act as a convergent read-out signifying commitment to the female pathway and responding to temperature in the European sea bass, represents a potential molecular link between the environment and the sex determination network (Navarro-Martín et al., 2011; Capel, 2017).

A factor promoting the formation and maintenance of the female phenotype, as it was commonly and strongly expressed in female gonads of all species of the present study, is that of *Figla* (*factor in germ-line alpha*). The *Figla* gene encodes a transcription factor that can coordinate the expression of the oocyte-specific zona pellucida genes and folliculogenesis (Liang et al., 1997). Besides its well-accepted role in promoting ovarian development, *Figla* has also the ability to repress spermatogenesis and, thus acts as an antagonistic factor of male genetic hierarchies (Qiu
et al., 2015). In hermaphrodites, its anti-parallel direction of expression to *Dmrt1* during the transition from male to female in a protandrous (Wu et al., 2012) and *vice versa* in a protogynous teleost (Miyake et al., 2012) has been suggested.

Likewise, the Luteinizing hormone (*Lh*) showed a strong female over-expression in all species examined here. Our findings, combined with studies in medaka, that have showed that *Lh* gene is involved in fish ovulation (Ogiwara et al., 2013), indicate that this gene has a conserved role to foster ovulation in wide range of fish even in phylogenetically remote species.

Another female-specific transcription factor is the growth differentiation factor-9 (*Gdf9*). *Gdf9* displayed sex-biased expression and was up-regulated in the ovaries of all studied species. This member of (TGFβ) superfamily, is an oocyte-specific factor in different vertebrates including fish (Elvin et al., 1999; Johnson et al., 2005; Liu and Ge, 2007). In mice, it is required in ovarian folicullogenesis (Dong et al., 1996), while in zebrafish it might also act as female-state protector by suppressing apoptosis; a mechanism responsible for degeneration of the ovarian tissue to pave the way to testicular development (Chen et al., 2017).

The common female “armory” is also assisted by a basic component of the *Wnt/b-catenin* pathway, the catenin beta 1 (*Ctnnb1*) gene. The present study confirms that the *Wnt/b-catenin* pathway has an active role in gonadal differentiation in all studied species (enrichment analysis). *Ctnnb1* showed a consistent pattern with a clear female-biased expression in these sparids thus proposing a conserved and acute role of *Ctnnb1* in preserving gonad’s female identity in this fish lineage.

Similarly, the role of retinoic acid (RA) signaling pathway in gonadal differentiation has been well-stated (Yamaguchi and Kitano, 2012; Rodríguez-Marí et al., 2013; Feng et al., 2015). The implication of this pathway in fish sex determination remains largely unknown, and data from different species show that the enzymes responsible for RA catabolism may vary (Feng et al., 2015). In the current study, two pathway members, *Cyp26a* and *Stra6* showed a female-biased expression in all studied sparids. On the other hand, a search for genes encoding RA receptors revealed two retinoic acid receptors (*Raraa* and *Rarb*) that were male-biased in all studied sparids.

Evidence suggests that both Wnt/b-catenin and retinoic acid pathways are implicated in the sex differentiation process in the studied species. In the present study it is shown that these pathways are not exclusively female-specific but rather “modules” of each of them could be used to assist the development of both sexes. Indeed, the female specificity of the Wnt/b-catenin pathway in fish, as is also known in mammals (Chassot et al., 2014), is not a standard (Nicol et al., 2012).

Two recognized members of the *Sox* gene family, *Sox3* and *Sox9*, were found among the commonly expressed gonadal genes in all the studied species. The present study demonstrates that *Sox3* had a high female-specific expression in the gonadal transcriptomes of all five sparids. *Sox* genes are highly conserved across evolution and encode a group of proteins that are well-known transcriptional regulators, involved in the decision of cell fates during development and in the control of diverse developmental processes (Pevny and Lovell-Badge, 1997). The

*Sox3* ability to act as a trigger and its suitability for generating new sex-determining genes (Marshall Graves and Peichel, 2010) renders it a unique candidate gene that could be used as a master regulator but also as a more downstream factor in the SD cascade. The exact function of *Sox3* in gonadal differentiation in this group of teleosts remains to be elucidated, whereas in the newly analyzed species, common dentex and gilthead seabream, it is verified that it clearly supports female development.

The other *Sox* gene, *Sox9*, was found among the commonly expressed male-biased gene in all five fish studied here. Information from two other hermaphroditic fish, the protogynous bluehead wrasse (Liu et al., 2015) and the protandrous clownfish (Casas et al., 2016), support *Sox9* male-specific expression. In contrast, the medaka *Sox9* is expressed in ovarian oocytes (Yokoi et al., 2002) whereas in zebrafish this gene has two homologs, one expressed in oocytes of the ovary and the other in testicular Sertoli cells (Chiang et al., 2001). Here, the medaka *Sox9* ortholog, while found in the gonadal transcriptomes of all sparids, did not show any sex-biased expression. Interestingly, the tilapia *Oreochromis niloticus* and the *Tetraodon nigroviridis* orthologs of *Sox9a* were unanimously male-biased. This gene in mammalian systems is involved in testis differentiation, as well as in male fertility maintenance (Jiang et al., 2013). The redundancy that exists among *Sox* transcription factors with different *Sox* members activating similar target genes (Windley and Wilhelm, 2015) and also the expansion of this gene family in teleost fishes after the teleost-specific WGD (Voldoire et al., 2017), render them convenient molecular tools. Thus, it is possible that slightly different genes or gene copies, following a lineage- or even species-specific functional divergence, could play the role that *Sox9* plays in mammals. Taken together, the results suggest that *Sox3* and *Sox9* have a conserved role in the formation of the functional female and the male phenotypes of the gonad, respectively, in Sparidae, but also suggest that both genes are “convenient” molecular players that could contribute to either sex side in teleosts.

If the sex determination/differentiation process could be viewed as a struggle for prevalence between the male and female hierarchies (Herpin and Schartl, 2011b), then *Dmrt1* is an evolutionary conserved “molecular soldier” that stands in the battlefront, fighting for the male side in both vertebrates and invertebrates (Raymond et al., 2000; Marshall Graves and Peichel, 2010; Herpin and Schartl, 2011a; Matson and Zarkower, 2012). In the present study, *Dmrt1* was a highly over-expressed gene in testes of all studied sparids. In accordance with our results, its characteristic male-biased expression has been reported in numerous other gonochorists (Marchand et al., 2000; Fernandino et al., 2008; Ijiri et al., 2008; Berbejillo et al., 2012) as well as hermaphrodite teleosts (He et al., 2003; Zhou and Gui, 2010; Liu et al., 2015; Casas et al., 2016). Thus, the *Dmrt1* role in the regulation of male development seems to be conserved not only in gonochorism but also in hermaphroditism. Characterized by a conserved DNA-binding motif [Doublesex- and Mab-3-related (DM) domain] and acting as a transcription factor, *Dmrt1* over-expression not only endorses male development by activating male-promoting genes, but also maintains male fates by repressing multiple female-promoting genes (Matson
et al., 2011). An important function of *Dmrt1* is to actively
and continuously suppress female gonad fates, thus maintaining
the male fate (Herpin and Schartl, 2011b; Matson et al., 2011).
The study of Shao et al. (2014) (Shao et al., 2014) implicates
*Dmrt1* with methylation-related epigenetic modifications and
suggests that is a critical gene that responds to environmental
change and triggers the sex reversal cascade in tongue sole.
Such action designates *Dmrt1* as a key player in “flexible”
reproductive mechanisms like those in the studied spard fishes
for which sex is naturally reversible. Apparently, *Dmrt1* together
with *Sox9* seems to signify the male character of the gonad in
this fish lineage. It is now evident that the environment can
affect the chromatin structure and DNA accessibility without
causing actual nucleotide changes (Jaenisch and Bird, 2003).
Genes encoding for epigenetic factors were abundant in the
transcriptomes of all sparids. 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 (*Hdac2*, *Hdac8*, *Hdac10*, *Hdac11*) were
found in the gonads and were sex-biased. We are still far from
understanding exactly how the epigenetic marks work, but since
expression of most of the afore-mentioned genes is regulated
from such factors, there is a strong indication that they play
an integral role in shaping and/or maintaining the gonadal
phenotype.

CONCLUSIONS

In this study, we reported the transcriptome of gonads and
brains of both sexes in two fishes of the Sparidae family, and
the results were studied jointly with published transcriptomes
of three other spard species showing different reproduction
modes (gonochoristic and protogynous hermaphrodites). This
approach permitted us to obtain a global view of the
transcriptome and to compare sex-specific expression patterns,
as well as sex-biased gene expression patterns within fishes
having alternative reproductive modes. The minor sex-related
expression differences in the brain, when studied as a whole,
might indicate a sexually un-biased tissue. On the contrary, in
the gonads there was a plethora of differentially expressed genes
between the sexes. This indicates functional divergence of the
ovary and testis, even though both gonads originated from the
same precursor tissue and despite the plasticity of the gonadal
tissue in sex-changing species. Focusing on pathways and genes,
we attempted to unveil the molecular means through which
these fish species develop into males or females, and observed
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 differentiation

and reproduction across teleosts. The knowledge obtained in
this study and tools developed during the process have laid the
groundwork for future experiments that can improve the sex
control of this species and help the in-depth understanding of
the complex process of sex differentiation in the more flexible
multi-component systems such as these studied here.

AVAILABILITY OF SUPPORTING DATA

The datasets generated and/or analyzed during the current study
are available from the corresponding author on request. The
raw sequence data has been submitted to the National Center
for Biotechnology Information (NCBI) Sequence Read Archive
(SRA) database (common dentex: BioProject ID PRJNA481721,
and gilthead seabream: BioProject ID PRJNA415944). The
transcriptomic sequences can be retrieved under the following
accession numbers: SRR7536480- SRR7536499, SRR7536500-
SRR7536502 for common dentex and SRR6223527-
SRR6223532, SRR6223535-SRR6223542 for gilthead
seabream.

AUTHOR CONTRIBUTIONS

AT, TM, CM, and CT designed experiments. AT, CT, NPs, and
CM conducted sampling. AT performed laboratory experiments.
AT, TM, NPu, and JL analyzed data. AT wrote the manuscript.
CT conceived and managed the project. All authors read and
approved the final version of the manuscript.

ACKNOWLEDGMENTS

We would like to thank Prof. Michalis Pavlidis for his very helpful
support and advice. We 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.
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. Last, the
authors express gratitude to the reviewers for their thorough
reading of our manuscript and their intuitive comments and
suggestions.

SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found
online at: https://www.frontiersin.org/articles/10.3389/fgene.
2018.00749/full#supplementary-material

REFERENCES

Altschul, S. F., Gish, W., Miller, W., Myers, E. W., and Lipman, D. J. (1990). Basic local alignment search tool. *J. Mol. Biol.* 215, 403–410. doi: 10.1016/S0022-2836(05)80360-2

Angelopoulou, R., Lavranos, G., and Manolakou, P. (2012). Sex determination strategies in 2012: towards a common regulatory model? *Reprod. Biol. Endocrinol.* 10:13. doi: 10.1186/1477-7827-10-13

Anonymous (1997). Guidelines for the treatment of animals in behavioural research and teaching. *Anim. Behav.* 53, 229–234. doi: 10.1006/anbe.1996.0293
Baroiller, J. F., and D'Cotta, H. (2016). The reversible sex of gonochoristic fish: insights and consequences. *Sex. Dev.* 10, 242–266. doi: 10.1159/000452362

Berbejillo, J., Martinez-Bengochea, A., Bedo, G., Brunet, F., Volff, J.-N., and Vizziano-Cantonnet, D. (2012). Expression and phylogeny of candidate genes for sex differentiation in a primitive fish species, the Siberian sturgeon, *Acipenser baeri*. *Mol. Reprod. Dev.* 79, 504–516. doi: 10.1002/mrd.22053

Böhne, A., Heule, C., Boileau, N., and Salzburger, W. (2013). Expression and sequence evolution of aromatase cyp19a1 and other sexual development genes in east african cichlid fishes. *Mol. Biol. Evol.* 30, 2268–2285. doi: 10.1093/molbev/mst124

Böhne, A., Sengstag, T., and Salzburger, W. (2014). 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. doi: 10.1093/gbe/evu200

Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. *Bioinformatics* 30, 2114–2120. doi: 10.1093/bioinformatics/btu170

Bourgon, R., Gentleman, R., and Huber, W. (2010). Independent filtering increases detection power for high-throughput experiments. *Proc. Natl. Acad. Sci. U.S.A.* 107, 9546–9551. doi: 10.1073/pnas.0914005107

Budd, A., Banh, Q., Domingos, J., and Jerry, D. (2015). Sex control in fish: approaches, challenges and opportunities for aquaculture. *J. Mar. Sci. Eng.* 3, 329–355. doi: 10.3390/jmse3020329

Buxton, C. D., and Garratt, P. A. (1990). Alternative reproductive styles in seabreams (Pisces: Sparidae). *Environ. Biol. Fishes* 28, 113–124. doi: 10.1007/BF00751031

Capel, B. (2017). Vertebrate sex determination: evolutionary plasticity of a fundamental switch. *Nat. Rev. Genet.* 18, 675–689. doi: 10.1038/nrg.2017.60

Casas, L., Saborido-Rey, F., Ryu, T., Michell, C., Ravasi, T., and Irigoien, X. (2016). Sex change in clownfish: molecular insights from transcriptome analysis. *Sci. Rep.* 6:35461. doi: 10.1038/srep35461

Chassot, A.-A., Gillot, I., and Chaboissier, M.-C. (2014). R-spondin1, WNT4, and the CTNNB1 signaling pathway: strict control over ovarian differentiation. *Reproduction* 148, R97–R110. doi: 10.1530/REP-14-0177

Chen, W., Liu, L., and Ge, W. (2017). Expression analysis of growth differentiation factor 9 (Gdf9/gdf9), anti-müllerian hormone (Amh/amh) and aromatase (Cyp19a1a/cyp19a1a) during gonadal differentiation of the zebrafish, *Danio rerio*. *Biol. Reprod.* 96, 401–413. doi: 10.1095/biolreprod.116.144964

Cheng, Y., Buffone, M. G., Kouadio, M., Goodheart, M., Page, D. C., Gerton, G. L., et al. (2007). Abnormal sperm in mice lacking the Taf7l gene. *Mol. Cell. Biol.* 27, 2582–2589. doi: 10.1128/MCB.01722-06

Chiang, E. F. L., Pai, C. I., Wyatt, M., Yan, Y. L., Postlethwait, J., and Chung, B. (2001). Two Sox9 genes on duplicated zebrafish chromosomes: expression of similar transcription activators in distinct sites. *Dev. Biol.* 231, 149–163. doi: 10.1006/dbio.2000.0129

Conesa, A., Götz, S., García-Gómez, J. M., Terol, J., Talón, M., and Robles, M. (2005). Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. *Bioinformatics* 21, 3674–3676. doi: 10.1093/bioinformatics/bti610

Desjardins, J. K., Klausner, J. Q., and Fernald, R. D. (2010). Female genomic response to mate information. *Proc. Natl. Acad. Sci. U.S.A.* 107, 21176–21180. doi: 10.1073/pnas.1010442107

Dong, J., Albertini, D. F., Nishimori, K., Kumar, T. R., Lu, N., and Matzuk, M. M. (1996). Growth differentiation factor-9 is required during early ovarian folliculogenesis. *Nature* 383, 531–535. doi: 10.1038/383531a0

Ellegren, H., and Parsch, J. (2007). The evolution of sex-biased genes and sex-biased gene expression. *Nat. Rev. Genet.* 8, 689–698. doi: 10.1038/nr.g2167

Elvin, J. A., Clark, A. T., Wang, P., Wolfman, N. M., and Matzuk, M. M. (1999). Paracrine actions of growth differentiation factor-9 in the mammalian ovary. *Mol. Endocrinol.* 13, 1035–1048. doi: 10.1210/mend.13.6.0310

EU (2010). *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.

Ezaz, M. T., Harvey, S. C., Boonphakdee, C., Teale, A. J., McAndrew, B. J., and Penman, D. J. (2004). Isolation and physical mapping of sex-linked AFLP

markers in Nile tilapia (*Oreochromis niloticus* L.). *Mar. Biotechnol.* 6, 435–445. doi: 10.1007/s10126-004-3004-6

Felip, A., Young, W. P., Wheeler, P. A., and Thorgaard, G. H. (2005). An AFLP-based approach for the identification of sex-linked markers in rainbow trout (*Oncorhynchus mykiss*). *Aquaculture* 247, 35–43. doi: 10.1016/j.aquaculture.2005.02.014

Feng, R., Fang, L., Cheng, Y., He, X., Jiang, W., Dong, R., et al. (2015). Retinoic acid homeostasis through aldh1a2 and cyp26a1 mediates meiotic entry in Nile tilapia (*Oreochromis niloticus*). *Sci. Rep.* 5:10131. doi: 10.1038/srep10131

Fernandino, J. I., Hattori, R. S., Shinoda, T., Kimura, H., Strobl-Mazzulla, P. H., Strüssmann, C. A., et al. (2008). Dimorphic expression of dmrt1and cyp19a1 (Ovarian Aromatase) during early gonadal development in pejerrey, *Odontesthes bonariensis*. *Sex. Dev.* 2, 316–324. doi: 10.1159/000195681

Forconi, M., Canapa, A., Barucca, M., Biscotti, M. A., Capriglione, T., Buonocore, F., et al. (2013). Characterization of sex determination and sex differentiation genes in latimeria. *PLoS ONE* 8:e56006. doi: 10.1371/journal.pone.0065006

Gamble, T., and Zarkower, D. (2012). Sex determination. *Curr. Biol.* 22, R257–R262. doi: 10.1016/j.cub.2012.02.054

Gentleman, R. C., Carey, V. J., Bates, D. M., Bolstad, B., Dettling, M., Dudoit, S., et al. (2004). Bioconductor: open software development for computational biology and bioinformatics. *Genome Biol.* 5:R80. doi: 10.1186/gb-2004-5-10-r80

Grabherr, M. G., Haas, B. J., Yassour, M., Levin, J. Z., Thompson, D. A., Amit, I., et al. (2011). Trinity: reconstructing a full-length transcriptome without a genome from RNA-Seq data. *Nat. Biotechnol.* 29, 644–652. doi: 10.1038/nbt.1883

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

Hattori, R. S., Murai, Y., Oura, M., Masuda, S., Majhi, S. K., Sakamoto, T., et al. (2012). A Y-linked anti-Mullerian hormone duplication takes over a critical role in sex determination. *Proc Natl Acad Sci U.S.A.* 109, 2955–2959. doi: 10.1073/pnas.1018392109

He, C. L., Du, J. L., Wu, G. C., Lee, Y. H., Sun, L. T., and Chang, C. F. (2003). Differential Dmrt1 transcripts in gonads of the protandrous black porgy, *Acanthopagrus schlegeli*. *Cytogenet. Genome Res.* 101, 309–313. doi: 10.1159/000074354

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

Herpin, A., and Schartl, M. (2011a). Dmrt1 genes at the crossroads: a widespread and central class of sexual development factors in fish. *FEBS J.* 278, 1010–1019. doi: 10.1111/j.1742-4658.2011.08030.x

Herpin, A., and Schartl, M. (2011b). Sex determination: switch and suppress. *Curr. Biol.* 21, R656–R659. doi: 10.1016/j.cub.2011.07.026

Hinfry, N., Nóbrega, R. H., Caulier, M., Baudiffier, D., Maillot-Maréchal, E., Chadili, E., et al. (2013). Cyp17a1 and Cyp19a1 in the zebrafish testis are differentially affected by oestradiol. *J. Endocrinol.* 216, 375–388. doi: 10.1530/JOE-12-0509

Ijiri, S., Kaneko, H., Kobayashi, T., Wang, D.-S., Sakai, F., Paul-Prasanth, B., et al. (2008). Sexual dimorphic expression of genes in gonads during early differentiation of a teleost fish, the Nile tilapia *Oreochromis niloticus*. *Biol. Reprod.* 78, 333–341. doi: 10.1095/biolreprod.107.064246

Ilias, A., Lagnel, J., Kapantaidaki, D. E., Roditakis, E., Tsigenopoulos, C. S., Vontas, J., et al. (2015). Transcription analysis of neonicotinoid resistance in Mediterranean (MED) populations of *B. tabaci* reveal novel cytochrome P450s, but no nAChR mutations associated with the phenotype. *BMC Genomics* 16:939. doi: 10.1186/s12864-015-2161-5

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

Jaenisch, R., and Bird, A. (2003). Epigenetic regulation of gene expression: how the genome integrates intrinsic and environmental signals. *Nat. Genet.* 33, 245–254. doi: 10.1038/ng1089
Jameson, S. A., Lin, Y.-T., and Capel, B. (2012). Testis development requires the repression of Wnt4 by Fgf signaling. *Dev. Biol.* 370, 24–32. doi: 10.1016/j.ydbio.2012.06.009

Jazin, E., and Cahill, L. (2010). Sex differences in molecular neuroscience: from fruit flies to humans. *Nat. Rev. Neurosci.* 11, 9–17. doi: 10.1038/nrn2754

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

Johnson, P. A., Dickens, M. J., Kent, T. R., and Giles, J. R. (2005). Expression and function of growth differentiation factor-9 in an oviparous species, *Gallus domesticus*. *Biol. Reprod.* 72, 1095–1100. doi: 10.1095/biolreprod.104.036822

Jolliffe, I. T. (2002). *Principal Component Analysis*, 2nd Edn. Berlin: Springer.

Kamiya, T., Kai, W., Tasumi, S., Oka, A., Matsunaga, T., Mizuno, N., et al. (2012). A trans-species missense SNP in Amhr2 is associated with sex determination in the tiger pufferfish, *Takifugu rubripes* (Fugu). *PLoS Genet.* 8:e1002798. doi: 10.1371/journal.pgen.1002798

Karube, M., Fernandino, J. I., Strobl-Mazzulla, P., Strüssmann, C. A., Yoshizaki, G., Somoza, G. M., et al. (2007). 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.* 307, 625–636. doi: 10.1002/jez.416

Kitano, T., Takamune, K., Kobayashi, T., Nagahama, Y., and Abe, S. I. (1999). 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. doi: 10.1677/jme.0.0230167

Kobayashi, Y., Horiguchi, R., Miura, S., and Nakamura, M. (2010). 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. doi: 10.1016/j.cbpa.2009.11.004

Kovács, B., Egedi, S., Bártfai, R., and Orbán, L. (2000). Male-specific DNA markers from African catfish (*Clarias gariepinus*). *Genetica* 110, 267–276. doi: 10.1023/A:1012739318941

Lagnel, J., Tsigenopoulos, C. S., and Iliopoulos, I. (2009). NOBLAST and JAMBLAST: new options for BLAST and a Java Application Manager for BLAST results. *Bioinformatics* 25, 824–826. doi: 10.1093/bioinformatics/btp067

Langmead, B., and Salzberg, S. L. (2012). Fast gapped-read alignment with Bowtie 2. *Nat Methods* 9, 357–359. doi: 10.1038/nmeth.1923

Larson, W. A., McKinney, G. J., Seeb, J. E., and Seeb, L. W. (2016). 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. doi: 10.1093/jhered/esw043

Li, B., and Dewey, C. N. (2011). RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. *BMC Bioinformatics* 12:323. doi: 10.1186/1471-2105-12-323

Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., et al. (2009). The sequence alignment/map format and SAMtools. *Bioinformatics* 25, 2078–2079. doi: 10.1093/bioinformatics/btp352

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

Liu, H., Lamm, M. S., Rutherford, K., Black, M. A., Godwin, J. R., and Gemmell, N. J. (2015). Large-scale transcriptome sequencing reveals novel expression patterns for key sex-related genes in a sex-changing fish. *Biol. Sex Differ.* 6:26. doi: 10.1186/s13293-015-0044-8

Liu, L., and Ge, W. (2007). Growth differentiation factor 9 and its spatiotemporal expression and regulation in the zebrafish ovary. *Biol. Reprod.* 76, 294–302. doi: 10.1095/biolreprod.106.054668

Loukovitis, D., Sarropoulou, E., Batargias, C., Apostolidis, A. P., Kotoulas, G., Tsigenopoulos, C. S., et al. (2012). Quantitative trait loci for body growth and sex determination in the hermaphrodite teleost fish *Sparus aurata* L. *Anim. Genet.* 43, 753–759. doi: 10.1111/j.1365-2052.2012.02346.x

Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. *Genome Biol.* 15:550. doi: 10.1186/s13059-014-0550-8

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

Maehiro, S., Takeuchi, A., Yamashita, J., Hiraki, T., Kawabata, Y., Nakasone, K., et al. (2014). Sexually dimorphic expression of the sex chromosome-linked genes cntfa and pdlim3a in the medaka brain. *Biochem. Biophys. Res. Commun.* 445, 113–119. doi: 10.1016/j.bbrc.2014.01.131

Mank, J. E., Hultin-Rosenberg, L., Webster, M. T., and Ellegren, H. (2008). The unique genomic properties of sex-biased genes: insights from avian microarray data. *BMC Genomics* 9:148. doi: 10.1186/1471-2164-9-148

Manousaki, T., Tsakogiannis, A., Lagnel, J., Sarropoulou, E., Xiang, J., Papandroulakis, N., et al. (2014). The sex-specific transcriptome of the hermaphrodite spard sharpsnout seabream (*Diplodus puntazzo*). *BMC Genomics* 15:655. doi: 10.1186/1471-2164-15-655

Marchand, O., Govoroun, M., D'Cotta, H., McMeel, O., Lareyre, J.-J., Bernot, A., et al. (2000). DMRT1 expression during gonadal differentiation and spermatogenesis in the rainbow trout, *Oncorhynchus mykiss*. *Biochim. Biophys. Acta* 1493, 180–187. doi: 10.1016/S0167-4781(00)00186-X.

Marshall Graves, J. A., and Peichel, C. L. (2010). Are homologies in vertebrate sex determination due to shared ancestry or to limited options? *Genome Biol.* 11:205. doi: 10.1186/gb-2010-11-4-205

Martínez, P., Bouza, C., Hermida, M., Fernández, J., Toro, M. A., Vera, M., et al. (2009). Identification of the major sex-determining region of turbot (*Scophthalmus maximus*). *Genetics* 183, 1443–1452. doi: 10.1534/genetics.109.107979

Matson, C. K., Murphy, M. W., Sarver, A. L., Griswold, M. D., Bardwell, V. J., and Zarkower, D. (2011). DMRT1 prevents female reprogramming in the postnatal mammalian testis. *Nature* 476, 101–104. doi: 10.1038/nature10239

Matson, C. K., and Zarkower, D. (2012). Sex and the singular DM domain: insights into sexual regulation, evolution and plasticity. *Nat. Rev. Genet.* 13, 163–174. doi: 10.1038/nrg3161

Matsuda, M., Nagahama, Y., Shinomiya, A., Sato, T., Matsuda, C., Kobayashi, T., et al. (2002). DMY is a Y-specific DM-domain gene required for male development in the medaka fish. *Nature* 417, 559–563. doi: 10.1038/nature751

Mattiske, D., Kume, T., and Hogan, B. L. M. (2006). The mouse forkhead gene Foxc1 is required for primordial germ cell migration and antral follicle development. *Dev. Biol.* 290, 447–458. doi: 10.1016/j.ydbio.2005.12.007

Memon, M. A., Anway, M. D., Covert, T. R., Uzumcu, M., and Skinner, M. K. (2008). Transforming growth factor beta (TGFβ1, TGFβ2 and TGFβ3) null-mutant phenotypes in embryonic gonadal development. *Mol. Cell. Endocrinol.* 294, 70–80. doi: 10.1016/j.mce.2008.08.017

Miyake, Y., Sakai, Y., and Kuniyoshi, H. (2012). Molecular cloning and expression profile of sex-specific genes, figla and Dmrt1, in the protogynous hermaphroditic fish, *Halichoeres Poecilopterus*. *Zool. Sci.* 29, 690–710. doi: 10.2108/zsj.29.690

Myosho, T., Otake, H., Masuyama, H., Matsuda, M., Kuroki, Y., Fujiyama, A., et al. (2012). Tracing the emergence of a novel sex-determining gene in Medaka, *Oryzias luzonensis*. *Genetics* 191, 163–170. doi: 10.1534/genetics.111.137497

Navarro-Martín, L., Viñas, J., Ribas, L., Díaz, N., Gutiérrez, A., Di Croce, L., et al. (2011). DNA methylation of the gonadal aromatase (cyp19a) promoter is involved in temperature-dependent sex ratio shifts in the European sea bass. *PLOS Genet.* 7:e1002447. doi: 10.1371/journal.pgen.1002447

Nicol, B., Guerin, A., Fostier, A., and Guiguen, Y. (2012). Ovary-predominant wnt4 expression during gonadal differentiation is not conserved in the rainbow trout (*Oncorhynchus mykiss*). *Mol. Reprod. Dev.* 79, 51–63. doi: 10.1002/mrd.21404

Ogiwara, K., Fujimori, C., Rajapakse, S., and Takahashi, T. (2013). Characterization of luteinizing hormone and luteinizing hormone receptor and their indispensable role in the ovulatory process of the Medaka. *PLoS ONE* 8:e54482. doi: 10.1371/journal.pone.0054482

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

Palaiokostas, C., Bekaert, M., Taggart, J., Gharbi, K., McAndrew, B., Chatain, B., et al. (2015). A new SNP-based vision of the genetics of sex determination in European sea bass (*Dicentrarchus labrax*). *Genet. Sel. Evol.* 47:68. doi: 10.1186/s12711-015-0148-y

Parsch, J., and Ellegren, H. (2013). The evolutionary causes and consequences of sex-biased gene expression. *Nat. Rev. Genet.* 14, 83–87. doi: 10.1038/nrg3376

Pevny, L. H., and Lovell-Badge, R. (1997). Sox genes find their feet. *Curr. Opin. Genet.* 7, 338–344. doi: 10.1016/S0959-437X(97)80147-5
Phillips, R. B., Park, L. K., and Naish, K. A. (2013). 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 3, 2289–2295. doi: 10.1534/g3.113.008078

Qiu, Y., Sun, S., Charkraborty, T., Wu, L., Sun, L., Wei, J., et al. (2015). Figla favors ovarian differentiation by antagonizing spermatogenesis in a teleosts, Nile tilapia (*Oreochromis niloticus*). *PLoS ONE* 10:e0123900. doi: 10.1371/journal.pone.0123900

Quevillon, E., Silventoinen, V., Pillai, S., Harte, N., Mulder, N., Apweiler, R., et al. (2005). InterProScan: protein domains identifier. *Nucleic Acids Res.* 33, W116–W120. doi: 10.1093/nar/gki442

Ranz, J. M., Castillo-Davis, C. I., Meiklejohn, C. D., and Hartl, D. L. (2003). Sex-dependent gene expression and evolution of the Drosophila transcriptome. *Science* 300, 1742–1745. doi: 10.1126/science.1085881

Raymond, C. S., Murphy, M. W., O'Sullivan, M. G., Bardwell, V. J., and Zarkower, D. (2000). Dmrt1, a gene related to worm and fly sexual regulators, is required for mammalian testis differentiation. *Genes Dev.* 14, 2587–2595. doi: 10.1101/gad.834100

Rodriguez-Mari, A., Cañestro, C., BreMiller, R. A., Catchen, J. M., Yan, Y. L., Postlethwait, J. H., et al. (2013). Retinoic acid metabolic genes, meiosis, and gonadal sex differentiation in zebrafish. *PLoS ONE* 8:e73951. doi: 10.1371/journal.pone.0073951

Ross, J. A., Urton, J. R., Boland, J., Shapiro, M. D., and Peichel, C. L. (2009). Turnover of sex chromosomes in the stickleback fishes (Gasterosteidae). *PLoS Genet.* 5:e1000391. doi: 10.1371/journal.pgen.1000391

Schmieder, R., and Edwards, R. (2011). Quality control and preprocessing of metagenomic datasets. *Bioinformatics* 27, 863–864. doi: 10.1093/bioinformatics/btr026

Shao, C., Li, Q., Chen, S., Zhang, P., Lian, J., Hu, Q., et al. (2014). Epigenetic modification and inheritance in sexual reversal of fish. *Genome Res.* 24, 604–615. doi: 10.1101/gr.162172.113

Shao, C.-W., Chen, S.-L., Scheuring, C. F., Xu, J.-Y., Sha, Z.-X., Dong, X.-L., et al. (2010). Construction of Two BAC libraries from half-smooth tongue sole cynoglossus semilaevis and identification of clones containing candidate sex-determination genes. *Mar. Biotechnol.* 12, 558–568. doi: 10.1007/s10126-009-9242-x

Sharma, E., Künstner, A., Fraser, B. A., Zipprich, G., Kottler, V. A., Henz, S. R., et al. (2014). Transcriptome assemblies for studying sex-biased gene expression in the guppy, *Poecilia reticulata*. *BMC Genomics* 15:400. doi: 10.1186/1471-2164-15-400

Shawki, H. H., Oishi, H., Usui, T., Kitadate, Y., Basha, W. A., Abdellatif, A. M., et al. (2018). MAFB is dispensable for the fetal testis morphogenesis and the maintenance of spermatogenesis in adult mice. *PLoS ONE* 13:e0190800. doi: 10.1371/journal.pone.0190800

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

Sun, F., Liu, S., Gao, X., Jiang, Y., Perera, D., Wang, X., et al. (2013). Male-biased genes in catfish as revealed by RNA-Seq analysis of the testis transcriptome. *PLoS ONE* 8:e68452. doi: 10.1371/journal.pone.0068452

Tao, W., Yuan, J., Zhou, L., Sun, L., Sun, Y., Yang, S., et al. (2013). Characterization of gonadal transcriptomes from Nile tilapia (*Oreochromis niloticus*) reveals differentially expressed genes. *PLoS ONE* 8:e63604. doi: 10.1371/journal.pone.0063604

Tsakogiannis, A., Manousaki, T., Lagnel, J., Sterioti, A., Pavlidis, M., Papandroulakis, N., et al. (2018). The transcriptomic signature of different sexes in two protogynous hermaphrodites: insights into the molecular network underlying sex phenotype in fish. *Sci. Rep.* 8:3564. doi: 10.1038/s41598-018-21992-9

Uchida, D., Yamashita, M., Kitano, T., and Iguchi, T. (2004). 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. doi: 10.1016/S1095-6433(03)00178-8

UniProt Consortium (2015). UniProt: a hub for protein information. *Nucleic Acids Res.* 43, D204–D212. doi: 10.1093/nar/gku989

Vilain, E., and McCabe, E. R. B. (1998). Mammalian sex determination: from gonads to brain. *Mol. Genet. Metab.* 65, 74–84. doi: 10.1006/mgme.1998.2749

Vizziano, D., Randuineau, G., Baron, D., Cauty, C., and Guiguen, Y. (2007). Characterization of early molecular sex differentiation in rainbow trout, *Oncorhynchus mykiss*. *Dev. Dyn.* 236, 2198–2206. doi: 10.1002/dvdy.21212

Voldoire, E., Brunet, F., Naville, M., Volff, J.-N., and Galiana, D. (2017). Expansion by whole genome duplication and evolution of the sox gene family in teleost fish. *PLoS ONE* 12:e0180936. doi: 10.1371/journal.pone.0180936

Wang, Z., Gerstein, M., and Snyder, M. (2009). RNA-Seq: a revolutionary tool for transcriptomics. *Nat. Rev. Genet.* 10, 57–63. doi: 10.1038/nrg2484

Warnes, G. R., Bolker, B., Bonebakker, L., Gentleman, R., Huber, W., Liaw, A., et al. (2009). *gplots: Various R Programming Tools for Plotting Data*. R Package version 2. Available online at: https://CRAN.R-project.org/package=gplots

Wilson, C. A., and Davies, D. C. (2007). The control of sexual differentiation of the reproductive system and brain. *Reproduction* 133, 331–359. doi: 10.1530/rep-06-0078

Windley, S. P., and Wilhelm, D. (2015). Signaling pathways involved in mammalian sex determination and gonad development. *Sex. Dev.* 9, 297–315. doi: 10.1159/000444065

Wu, G.-C., and Chang, C.-F. (2013). The switch of secondary sex determination in protandrous black porgy, *Acanthopagrus schlegeli*. *Fish Physiol. Biochem.* 39, 33–38. doi: 10.1007/s10695-012-9618-0

Wu, G.-C., Tomy, S., Lee, M.-F., Lee, Y.-H., Yueh, W.-S., Lin, C.-J., et al. (2010). Sex differentiation and sex change in the protandrous black porgy, *Acanthopagrus schlegeli*. *Gen. Comp. Endocrinol.* 167, 417–421. doi: 10.1016/jygcn.2009.11.003

Wu, G. C., Chiu, P. C., Lin, C. J., Lyu, Y. S., Lan, D. S., and Chang, C. F. (2012). Testicular dmrt1 is involved in the sexual fate of the ovotestis in the protandrous black porgy. *Biol. Reprod.* 86:41. doi: 10.1095/biolreprod.111.095695

Xu, D., Shen, K.-N., Fan, Z., Huang, W., You, F., Lou, B., et al. (2016). 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. doi: 10.1016/j.gdata.2016.01.007

Yamaguchi, T., and Kitano, T. (2012). 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. doi: 10.1016/j.bbrc.2012.02.012

Yano, A., Guyomard, R., Nicol, B., Jouanno, E., Quillet, E., Klopp, C., et al. (2012). An immune-related gene evolved into the master sex-determining gene in rainbow trout, *Oncorhynchus mykiss*. *Curr. Biol.* 22, 1423–1428. doi: 10.1016/j.cub.2012.05.045

Yates, A., Akanni, W., Amode, M. R., Barrell, D., Billis, K., Carvalho-Silva, D., et al. (2015). Ensembl 2016. *Nucleic Acids Res.* 44, D710–D716. doi: 10.1093/nar/gkv1157

Yokoi, H., Kobayashi, T., Tanaka, M., Nagahama, Y., Wakamatsu, Y., Takeda, H., et al. (2002). sox9 in a teleost fish, medaka (Oryzias latipes): evidence for diversified function of Sox9 in gonad differentiation. *Mol. Reprod. Dev.* 63, 5–16. doi: 10.1002/mrd.10169

Zhang, Z., Wang, Y., Wang, S., Liu, J., Warren, W., Mitreva, M., et al. (2011). Transcriptome analysis of female and male *Xiphophorus maculatus* Jp 163 A. *PLoS ONE* 6:e18379. doi: 10.1371/journal.pone.0018379

Zhou, L., and Gui, J. F. (2010). Molecular mechanisms underlying sex change in hermaphroditic groupers. *Fish Physiol. Biochem.* 36, 181–193. doi: 10.1007/s10695-008-9219-0

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Tsakogiannis, Manousaki, Lagnel, Papanikolaou, Papandroulakis, Mylonas and Tsigenopoulos. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.