De novo gonad transcriptome analysis of the common littoral shrimp Palaemon serratus: novel insights into sex-related genes

The common littoral shrimp Palaemon serratus (Pennant, 1777) is a crustacean decapod with a geographical distribution ranging the Atlantic Ocean (from Scotland and Denmark to Mauritania, including Azores, Madeira and Canary Islands) and the entire Mediterranean Sea and the Black Sea [1]. This species inhabits the intertidal and subtidal soft-sediment of estuaries in the reproductive season, and rocky bottoms covered with seagrass and algae [2, 3]. Palaemon serratus fishing activity is crucial in some European communities, mainly around the British Isles, France and northern Spain [4]. In Galicia (NW Spain) the volume of catches varies from 47.6 to 90.7 tons traded per year, what equals to worth of 2 million euros per year in this region (data obtained from https://www.pescadegalicia.gal/↗ on 26 Sep 2018, Xunta de Galicia). The high commercial value of this species could possibly lead to overfishing [4]. Implementation of proper management measures that ensure a sustainable exploitation will prevent depletion or genetic deterioration of wild fisheries [5, 6]. Aquaculture practices might improve P. serratus production at once reduce the fishing pressure over the wild populations.

In the field of aquaculture, reproductive traits are considered economically significant. Hence, understanding sexual and reproductive development is necessary to obtain successful and sustainable cultures, to increase seed quality or to breed genetically improved lines [7, 8]. For instance, as sex dimorphism in growth is common in crustaceans, monosex aquaculture of commercially relevant species is especially interesting. In monosex populations the yield is increased because energy from reproduction is invested in growth, resulting in larger-size individuals than in sex-mixed cultures [9, 10]. A better knowledge about the genetics of crustacean sexual development facilitates the application of biotechnological strategies, such as sex-change induction, benefiting productivity [7].

Sexual development includes sex determination and sex differentiation processes. In Decapoda, sex is determined by the initiation of a genetic cascade triggered by a master sex-determining gene. Downstream genes in this cascade act as sex-regulator genes, leading to the sex differentiation pathway, which in turn results in a sex-specific phenotype development [11]. Due to the lack of genomic information in crustacean decapods, sex-determining genes have not been identified and even sex-related genes have been rarely reported [9]. Several genes are considered preliminary candidates to be implicated in decapod sex determination but since these genes have been identified through homology screening, the list is heavily biased by genes characterized in model species as Drosophila melanogaster, Caenorhabditis elegans and Mammalia (e.g. Sxl, Tra, Tra-2, Dsx, Fem-1 or Sry; see review in [11]). Among these candidates, it is noteworthy to highlight that Dmrt genes (doublesex and male abnormal-3-related transcription factors) have been noted as the only gene family with a conserved function in sex determination across metazoans [11, 12], so they are particularly intriguing. About sex differentiation, the insulin-like androgenic factor (IAG) is a well characterized hormone with a conserved central role in Malacostraca. Iag expression in the androgenic gland (AG) leads sexual differentiation to maleness by governing the onset of testicular and secondary-sex characteristics development in males. Upstream in the sex differentiation pathway, an array of neuropeptides secreted by the eyestalk regulates the Iag expression [11]. It was shown that AG-implanted females become males, and conversely AG-ablated males turned into females (see review in [13]). As this surgical procedure implies a high mortality rate, it has been recently achieved that Iag-silenced males shifted phenotypically into females in the prawn Macrobrachium rosenbergii [14]. Thus, the biotechnological manipulation of the expression of sex determination or sex differentiation genes has the potential to improve the aquaculture production, such as creating monosex populations among other possibilities. Nevertheless, prior to the implementation of any genetic manipulation technique, it is necessary a depth-understanding about the genetic factors underlying the sexual development in P. serratus.

RNA-Seq greatly enhances the capability for gene discovery in non-model organisms where genomic data is not available. Transcriptome profiling using high-throughput sequencing allows the identification of transcripts involved in biological processes [15]. Comparative transcriptomics and differential expression analyses (DEA) between female and male reproductive tissues enable the detection of transcripts with sex-biased and sex-specific expression. In fact, sex determination and sex differentiation genes have been identified in several commercial decapod species thanks to transcriptomic analyses of certain tissues, including gonads [9, 16–18].

Genetic studies in P. serratus are scarce and mainly focused in its population genetics and cytogenetics. Solely one transcriptomic work is available for this species, providing data to study larval development and metamorphosis [19]. Regarding to sex determination in P. serratus, it is only known that heteromorphic sex chromosomes are absent [20]. No sex determination or sex differentiation genes or pathways have been reported for this shrimp. Accordingly, the aim of the present study was to identify candidate genes to be involved in the sexual development of P. serratus. For this purpose, an ovary and testis transcriptome was assembled and annotated from Illumina high-throughput sequencing reads, and genes with differential expression between ovarian and testicular tissues were studied. To the best of our knowledge, this is the first work that address the transcriptome profile analysis of both male and female P. serratus gonads, providing new data about sex-related genes in this species.

Palaemon serratus is a relevant commercial species in some countries such as UK, Ireland, France, Spain and Portugal [21]. The lack of genomic information about P. serratus hinders the application of potential aquaculture techniques, especially those focused on reproductive traits related to sex dimorphism. Therefore, in the present work we attempted to unravel sex-related genes featured in sex determination, sex differentiation and/or gonadal development using a RNA-Seq approach. This is the first transcriptome analysis of the gonads of a Palaemon species and the first work that provides data about sex-related genes in P. serratus. Within Palaemonidae, sex-related genes have been only studied in two Macrobrachium species [8, 9, 22]. Here, a reference gonadal transcriptome was obtained using ovary and testis reads. Statistics indicated a high quality de novo assembly but only the 29.57% of the transcriptome could be annotated due to the scarcity of crustacean genomic sequences. Precisely, the most expressed up-regulated genes both in ovary and in testis were not annotated. As the non-annotated transcripts correspond to unknown novel transcripts or to unreviewed transcripts, these unannotated differentially expressed sequences deserve attention in future functional analyses about sex-related genes in this shrimp. Additionally, the complete transcriptome annotation is provided, largely increasing the sequence resources available for this species.

Sexual development includes several processes orchestrated by a variety of regulators. Overall, sex determination and sex differentiation are intricated processes not always clearly distinguishable because their signalling cascades can be integrated [11]. Sex determination mechanisms are widely divergent in animals, even between closely related species. This variability is due to the rapid evolution that sex-biased genes experience [11, 23] and is reflected to crustaceans, where is not uncommon for findings in sex-related genes to differ among species. This can be linked to there is no conserved sex determination pathway in decapods and it likely evolved independently several times, making difficult to trace master sex-regulators [24]. Back to P. serratus, heteromorphic sex chromosomes are absent [20] and none sex determination system has been documented. It was suggested that there is a sex chromosome dosage compensation mechanism involving Msl3 gene in the tissues of the palaemonid Macrobrachium nipponense [25] as in Drosophila. Unfortunately, Msl3 was not a DEG between females and males in P. serratus, giving no hint to whether a heterogametic sex exists in this shrimp or not and what is it. Hence, we focused our efforts to study genes described as ‘sex-related’, mainly in crustaceans.

Orthologs of sex determination genes in the model arthropod Drosophila were found in our transcriptome database. Sxl and Tra-2 orthologs were detected without sex differential expression meanwhile Tra, Dsx and Fru were absent. In Drosophila sex is ruled by the genetic pathway Sxl-Tra/Tra-2-Dsx/Fru, being Sxl the master sex determinant gene [12, 26]. It has been proposed that crustaceans may adopt the Drosophila sex determination pathway given the findings reported in some species as Penaeus monodon [27], Macrobrachium nipponense [28], Penaeus chinensis [29], Penaeus vannamei [30] and Eriocheir sinensis [31]. However, our data are in line with those that stand that these genes do not act in decapods as they do in insects (see review in [7]), as suggested for the lobster Sagmariasus verreauxi [32] or for the crab Scylla paramamosain [33].

In the nematode Caenorhabditis elegans, Fem-1 is a component of the sex determination signalling pathway that promotes the male phenotype [34, 35]. There are studies in decapods that pointed out the putative role of Fem-1 genes in male sex determination [18, 36]. Orthologs of the three members of the Fem-1 family were found in P. serratus but none of them was up-regulated in the testis. Fem-1a and Fem-1b were not DEGs between sexes whilst Fem-1b was found to be slightly up-regulated in the ovary. Particularly to Fem-1b, [37] detected a higher expression of this gene in the testes than in the ovaries of the prawn Macrobrachium nipponense, which is the opposite of what we found in P. serratus. Also in M. nipponense, [38] reported a ovary-specific Fem-1 gene that could be involved in sex determination or differentiation and in ovarian maturation in this species. Fem-1b ortholog was up-regulated in the ovary of P. serratus but it was not an ovary-specific gene given that it also exhibits a considerable expression in testis. We conclude that similarly to in Scylla paramamosain [17] and in Penaeus vannamei [39], our results seem to indicate that whether Fem-1 genes are involved in the sexual development in P. serratus has yet to be established given that they are expressed in both gonads.

The male-determining gene in most mammals is the Y chromosome Sry gen [40, 41]. SRY along with SF-1 induces testicular development through the activation of the transcription factor SOX9 [42]. SOX9 up-regulates via FGF9 the expression of the Dmrt1 gene, which is the major male sex differentiation gene, promoting testis development and maintenance [43]. Expression of Sry, Sf-1, Sox9 and Fgf9 was not detected in the gonads of P. serratus. Some members of the Sox family were DEGs, Sox8 was up-regulated in ovaries and Sox5, Sox14 and Sox15 were specifically expressed in testes. The up-regulation of Sox8 in females was unexpected as this gene has been related with testis development [12, 32]. Both Sox5 and Sox14 were previously identified as genes involved in male sex differentiation with expression in testis tissues [15, 44–46] but Sox15 has never been defined as a testis gene before. Nevertheless, the most relevant finding in P. serratus gonadal transcriptome is the testis-specific expression of the Dmrt1 gene.

In some vertebrate species Dmrt1 has been qualified as a sex-determining gene [47, 48]. Moreover, Dmrt1 paralogs are the master sex determinant genes in medaka fish [49], frogs [50] and, recently in Sagmariasus verreauxi [11] as the first invertebrate species in which Dmrt1 determinates the sex. DMRTs are transcription factors characterized by the presence of a DM domain DNA binding motif. The relationship between DM domain genes and sex has been deeply investigated and as it is thought that their ancestral function was likely to determine gonadal sex and they subsequently expanded to control sexual dimorphism in other tissues [43]. Dmrt is the only gene family with a conserved function in sex determination across Animalia [11] and orthologs have been identified with a testis-restricted expression in the transcriptome of a few decapod species as the crabs Eriocheir sinensis [16] and Scylla paramamosain [17]. Keeping this in mind, the testis-specific Dmrt1 ortholog found in P. serratus should be considered the best candidate gene to be involved in the sex determination of this species. Future efforts should be directed to functionally characterize the Dmrt1 gene and to pursue upstream regulators and downstream targets. Another Dmrt gene was found in the transcriptome with an extremely low expression (TPM < 1), the Dmrt11E gene that has been previously detected in some decapods. This gene exhibited a testis-biased expression in Macrobrachium rosenbergii [51] and an androgenic gland-biased expression in Sagmariasus verreauxi [32], and in both cases it was suggested that Dmrt11E is a male differentiation regulator via IAG. As the Dmrt11E expression in P. serratus gonads is very low, another tissue should be the primary site of expression instead of testis, likely the AG. Owing to its proved relationship with the IAG in other species, expression of Dmrt11E should be studied in different organs of P. serratus.

The IAG is the key regulator of male sex differentiation in the members of Malacostraca and its expression takes place exclusively in the androgenic gland (AG) of males. The Iag gene was characterized in several crustaceans, e.g. in the prawn Penaeus monodon [52], in the shrimps Penaeus vannamei [53], Macrobrachium lar, Palaemon paucides and P. pacificus [54], or in the spiny lobsters Sagmariasus verreauxi and Jasus edwardsii [32]. The expression level of the Iag gene in our gonadal transcriptome was very low (TPM = 0.84), likely because in Palaemon species the AG is located along the sperm ducts and not in the testis [54], and they should not be dissected along with the testis. A better knowledge about Iag is crucial because its genetic manipulation-based biotechnology has the potential to dramatically transform the entire aquaculture industry [55]. Monosex culture has the potential to enhance the production because energy for reproduction is allocated to growing, so the individuals reach higher sizes [54]. As female shrimps grow larger and faster, all-female population cultures are preferred for P. serratus. It has been proved in different decapod species that AG removal feminizes males [13], but this surgical procedure frequently entails mortality, so get monosex cultures by genetic manipulation is highly attractive. All-male populations were achieved for Macrobrachium rosenbergii silencing Iag using RNA interference [14]. To obtain all-female populations in P. serratus we suggest exploring the manipulation of the Iag gene to induce female sex-reversal. We provide the first Iag sequence for P. serratus and, even though is a partial sequence, it has the potential to pave the way to further biotechnological approaches that enable the production of female monosex cultures. These aquaculture strategies may enhance P. serratus production and at the same time prevent the genetic deterioration of the wild stocks caused by overfishing.

Genes referred in literature as ‘testis development’ genes were also investigated. KIFC1 is a C-terminal kinesin motor protein that participates in acrosome biogenesis and nuclear reshaping during spermiogenesis in the palaemonids Macrobrachium nipponense [56] and Palaemon modestus [57] among other crustacean species. Kifc1 gene showed a high expression in the testis of M. nipponense and P. modestus but in both species this gene was also being expressed in other tissues, likely taking part in vesicle transportation processes [56], so is not strange that Kifc1 was not a DEG between ovary and testis in P. serratus. Temporal and spatial expression of Kifc1 during spermiogenesis in P. serratus could be address since its protein is vital in the formation of the acroframosome, an exclusive structure of caridean shrimp spermatids. Concerning to DMC1, it plays a major role in meiotic recombination and has been associated to spermatogenesis in crustaceans [27]. Unlike in the crawfish Procambarus clarkii [18] and the crab Portunus trituberculatus [58], Dmc1 was not up-regulated in the testis but in the ovary of Palaemon serratus. It is important to highlight that Dmc1 was a gonad up-regulated gene, most likely because it is expressed in meiotic germ cells [59], but its role in the spermatogenesis of P. serratus should not be directly attributed without further testing. Another gene related with germ cell development is vasa, an ATP-dependent RNA helicase. Since vasa plays a role in both oogenesis and spermatogenesis its expression was exclusively detected in the gonads of the shrimp Penaeus vannamei [60] or of the crab Scylla paramamosain [61]. In this sense it was not surprising that vasa was a gonad up-regulated gene respect to non-gonad tissues but not a DEG between the ovary and the testis of P. serratus.

Regarding to female sex determination, Foxl2 gene encodes a conserved forkhead transcription factor preferentially expressed in the ovary of vertebrates, controlling the ovarian differentiation and maintenance by repression of testis-specific genes [62, 63]. If Dmrt1 is present, Foxl2 expression is repressed, but in the absence of Dmrt1, Foxl2 inhibits the male developmental pathway and promotes the female. The expression of Foxl2 showed a changing pattern among crustacean species, i.e. up-regulated in ovary [18], not DEG between sexes [31] or even up-regulated in testis [64]. Foxl2 was not a DEG between ovary and testis in P. serratus, which supports that the role of Foxl2 in sex determination in invertebrates remains unclear [65, 66]. Dmrt1 also acts repressing the RSPO1-WNT4-β-catenin signalling pathway, another female sex determination cascade that promotes ovary development in vertebrates independently and complementary to the Foxl2-leading pathway [67–69]. Given that Rspo1 was not found in the transcriptome of P. serratus and Wnt4 and β-catenin were detected as not DEGs, the existence and function of this pathway is unknown in this shrimp, as it was already advanced for other decapods with orthologs found [31, 33]. Thereby, vertebrate pathways leadered by Foxl2 and Rspo1 do not seem to determine female development in P. serratus in the light of our data. However, further experiments should confirm these lacks of function in sexual development.

Vitellogenesis, the production and accumulation of yolk, is crucial to oogenesis and ovarian maturation. In oviparous vertebrates vitellogenin synthesis is enhanced by 17β-estradiol E2, with the estrogen receptor (ER) and the HSP90 acting as mediators [70–73]. The elements of the E2-ER-HSP90 pathway were found in the transcriptome of decapod species [7, 64, 74] but only Hsp90 showed a higher expression in ovaries than in testes in the crab Scylla paramamosain [17]. An estrogen-related receptor gene (Err) was found in P. serratus without differential expression and Hsp90 was up-regulated in the testicular tissue, so our results agree that more studies are necessary to clarify if the E2-ER/ERR-HSP90 pathway exists in crustaceans and whether vitellogenesis it is regulated by estrogen-like hormones as it is in vertebrates. Another regulatory pathway that stimulates ovarian development and vitellogenesis in some decapods involves methyl farnesoate (MF), a crustacean juvenile hormone analogue [75, 76]. Farnesoic acid-O-methyl transferase (FAOMeT) encondes the enzyme that catalyzes formation of MF and it was up-regulated in the ovaries of P. serratus respect to testes and non-gonad tissues, indicating the putative role of this hormone in the ovarian maturation of the species. Vitellogenin (Vg) and vitellogenin receptor (VgR), the main vitellogenesis genes, were up-regulated in the ovary of P. serratus as in multiple decapod species (e.g. [24] or [58]). The expression of VgR was higher than the expression of Vg, likely because the hepatopancreas is considered the primary site of VG production instead of the ovary while the VGR allow VG uptake from the hemolymph by oocytes. Vitelline membrane outer layer protein 1-like gene (Vmo1) was strangely up-regulated in males, a finding also reported for the crab Eriocheir sinensis [31].

Other ‘ovary development’ genes were also explored in the gonadal transcriptome of P. serratus. Since prostaglandins (PGs) have been described as factors that promote ovary development in crustaceans [77], PG genes were studied in P. serratus. HPGDS and PGES2 showed an ovary-biased expression, so they might have an implication in female gonad development [78, 79]. Although PTGR1 is involved in the ovary development in the shrimp Penaeus monodon [80], it was not preferentially expressed in the ovaries of Palaemon serratus. Cathepsins have been also related with ovarian development in some crustaceans [22, 81, 82]. Cathepsin D is a needed protein for the formation of the yolk in vertebrates [83] and its gene was the only cathepsin gene up-regulated in the ovary of P. serratus. Other up-regulated genes in the ovary of P. serratus that are required for ovary maturation in other crustaceans were chorion peroxidase [31], profilin [84], Smad4 [85], Gnrhr [86] and Pgmrc1 [87]. Expression of Hsc70 was also up-regulated in the ovary of P. serratus, suggesting a putative role in reproductive events as in Macrobrachium rosengergii, where Hsc70 was enriched in the ovary [88]. Fst is known to be involved in folliculogenesis and ovary development in vertebrates [89], but it was not a DEG between the gonads of P. serratus, similarly to Eriocheir sinensis [31]. The expression of Fst should be examined throughout the stages of the ovarian development to evaluate whether it participates in it or not. Cytochrome P450 aromatase (Cyp19a) is also essential in the female gonadal development in vertebrates, converting androgens in estrogens. Several genes belonging to cytochrome P450 superfamily were detected in the transcriptome of P. serratus, but none of them was aromatase, an absence also found in the palaemonid Macrobrachium rosenbergii [90].

Lastly, two genes encoding two members of the CHH-superfamily were detected in the gonadal transcriptome with a very low and not differential expression between sexes: the crustacean hyperglycemic hormone (Chh) and the molt inhibiting hormone (Mih). CHH neuropeptides are multifunctional hormones with roles in reproduction, regulating AG proliferation or MF production among other activities (see review in [77]). The eyestalk is the preferential site of production of these neurohormones, but it was repeatedly demonstrated that these genes are also expressed in multiple tissues [91–93], including gonads in some species as in P. serratus. Likewise, as recent studies have demonstrated that ecdysteroids regulate vitellogenesis, ovarian maturation and spermatogenesis in decapods (see review in [78]), it is also interesting to highlight the expression of the ecdysone receptor (EcR) gene in both female and male P. serratus gonads.

This study encompasses the first large-scale RNA-Seq and comprehensive transcriptome analysis of Palaemon serratus gonads. More than 442 million of clean reads were obtained, 39,186 transcripts were assembled and annotated and 11,087 out of them were found to be DEGs between ovary and testis. Sex-related genes were identified and their expression between sexes was studied. A wide inventory of sex-related genes is provided and thoroughly discussed in the framework of previous findings in other crustacean species. This is the first time that sex-related genes have been addressed in a Palaemon species, so this transcriptomic analysis will facilitate further experimental research that aimed to delve into the sex determination, sex differentiation and gonadal development mechanisms of P. serratus and close species. The candidate genes to be involved in sexual development might also shed light about the evolution of sex-regulators in crustaceans. Furthermore, we report some particularly interesting genes towards investigating future aquaculture applications for P. serratus.

Specimens of P. serratus used in this study were collected inshore from the Ártabro Gulf (43°22′00″N, 8°28′00′W) in the northwest of Spain using a fish trap. Animals were carried alive to the Aquarium Finisterrae dependencies (A Coruña, Spain) where they were kept at 18 °C in an aerated aquarium while they were sorted into sexes. According to [3], sex was determined by the presence (in males) or absence (in females) of the masculine appendix on the endopodite of the second pleopod. Shrimps (3–4.5 g body weight) were anesthetized on ice for 5 min prior to be sacrificed by dissection, and then gonads from three adult females and three adult males were quickly removed and directly immersed in liquid nitrogen. The development stage of the gonads was fully mature in all individuals. RNA isolation and library construction were carried out at AllGenetics & Biology SL (A Coruña, Spain). Total RNA was extracted from the six samples by grinding them with a mortar and pestle under liquid nitrogen. The resulting powder was used for the extraction using NZYTech’s Total RNA isolation kit (NZYTech). Pure RNA was eluted in a final volume of 30 μL and then quantified and quality-checked in an Agilent 2100 Bioanalyzer (Agilent) using the Agilent RNA 6000 Kit (Agilent). Illumina’s TruSeq Stranded mRNA Library Prep Kit (Illumina) was used to prepare six cDNA libraries following the manufacturer’s instructions. Thus, one library per sample was prepared, i.e. three libraries from ovary tissues and three libraries from testis tissues were prepared, or in other words three biological replicates per sex. The fragment size distribution and concentration of the libraries were checked in the Agilent 2100 Bioanalyzer using the Agilent DNA1000 Kit (Agilent). Libraries were quantified using the Qubit dsDNA BR Assay Kit (ThermoFisher Scientific). All libraries were pooled in equimolar amounts, according to the quantification data. The pool was sequenced in two different lanes of a HiSeq 4000 PE100 platform (Illumina).

Raw reads quality control was performed using FastQC v0.11.5 [94]. Trimmomatic v0.35 [95] was used for raw reads trimming. Primer/adaptor sequences were removed and the first 15 bp of the reads were cut. Trimmed reads shorter than 40 bp were discarded. Both female and male gonads trimmed reads were assembled together to obtain a single transcriptome that included ovary and testis transcripts. Trinity software v2.4.0 [96] was used for de novo assembly of these high-quality short reads using default parameters settings (Kmer = 25). Assembled transcriptome completeness was assessed with BUSCO v3.0.2 [97] using the Arthropoda database as reference. EvidentialGene tr2aacds pipeline [98] was used to reduce the transcriptome redundancy. Then BUSCO was run again to check the duplication level. Gene expression given as Transcript Per Million (TPM) was calculated by mapping back all ovary and testis trimmed reads together on the gonads assembled transcriptome (parameters: minimum allowed length fraction = 0.75, similarity fraction = 0.95 and maximum number of matching contigs = 4) using the RNA-Seq tool of the CLC Genomics Workbench v11.0 (Qiagen). Only transcripts with a TPM ≥ 1 were included in the definitive transcriptome. Ribosomal and mitochondrial contigs were identified by BLASTn [99] against the NCBI non-redundant (nr) database and they were also excluded from the final transcriptome. Functional annotation was carried out using the Trinotate v3.1.1 workflow (https://trinotate.github.io↗). In detail, sequence similarity of the assembled transcripts was evaluated using BLASTx [99] against the UniProtKB/Swiss-Prot database (E-value cut-off of 1e-5). TransDecoder v5.3.0 [100] was used to identify putative protein coding regions, including homology options as retention criteria for the candidates ORFs. Predicted ORFs were identified by BLASTp [99] queries against the UniProtKB/Swiss-Prot database (E-value cut-off of 1e-5). Protein functional domains were identified using HMMER3 [101] against the PFAM domain database. Signal peptide and transmembrane domains prediction was performed with SignalP v4.1 [102] and TMHMM v2.0 [103], respectively. WEGO [104] was used to plot the Gene Ontology (GO) functional classification and distribution of the annotated transcripts.

Trimmed reads of each gonad sample were mapped back separately on this reduced high-quality set of ovary and testis transcripts (parameters: minimum allowed length fraction = 0.75, similarity fraction = 0.95 and maximum number of matching contigs = 4) using the RNA-Seq tool of the CLC Genomics Workbench v11.0 (Qiagen). Thus, expression of each transcript in each sample was calculated as Transcript Per Million (TPM). Subsequently, a DEA between ovary and testis samples was carried out in order to identify DEGs between ovaries and testes using the Differential Expression for RNA-Seq tool of the CLC Genomics Workbench v11.0 (Qiagen). GO enrichment analysis was conducted to reveal GO terms significantly enriched in DEGs using the CLC Genomics Workbench v11.0 (Qiagen). Transcriptomic SRA data (NCBI Sequence Read Archive) from other P. serratus tissues, larvae (SRR4341161–2) and muscle (SRR4341163–4), was used to identify differentially expressed genes (DEGs) between gonad tissues (ovary and testis) and non-gonad tissues (larvae and muscle). Firstly, SRA reads quality was checked with FastQC v0.11.5 [94] and Trimmomatic v0.35 [95] was used to trim the reads as follows: HEADCROP:15 TRAILING:25 MINLEN:40. SRA trimmed reads were mapped on the non-redundant and reduced gonads transcriptome to calculate the gene expression (TPM) of the assembled transcripts in the larvae and muscle samples. Pairwise differential expression analyses (DEAs) were performed between gonad and non-gonad samples (ovaries vs. larvae, ovaries vs muscle, testes vs larvae and testes vs muscle) using the CLC Genomics Workbench v11.0 (Qiagen). Up-regulated expressed genes in gonad tissues were tagged with a ‘G’ (gonad up-regulated) in the transcriptome annotation table.