De Novo Assembly, Characterization and Comparative Transcriptome Analysis of the Gonads of Jade Perch (Scortum barcoo)

Simple Summary Jade perch (Scortum barcoo) has been affirmed as an economically important species for aquaculture. However, research on its sex and reproductive regulation is still lacking. Herein, the first gonad transcriptomes of juvenile jade perch were identified through Illumina Novaseq technology. According to the comparative transcriptome data of ovaries and testis, a large number of differential expression genes (DEGs) mainly involved in the steroidogenesis pathway, gonad development and differentiation, gametogenesis and gamete maturation were identified. As predicted, the expression patterns of these DEGs in jade perch were similar to well-known sex-related genes in other species. These results will provide valuable information for further studies on sex determination and gonadal development in teleost fish. Abstract Due to the high meat yield and rich nutritional content, jade perch (Scortum barcoo) has become an important commercial aquaculture species in China. Jade perch has a slow growth rate, taking 3–4 years to reach sexual maturity, and has almost no difference in body size between males and females. However, the study of its gonad development and reproduction regulation is still blank, which limited the yield increase. Herein, the gonad transcriptomes of juvenile males and females of S. barcoo were identified for the first time. A total of 107,060 unigenes were successfully annotated. By comparing male and female gonad transcriptomes, a total of 23,849 differentially expressed genes (DEGs) were identified, of which 9517 were downregulated, and 14,332 were upregulated in the testis. In addition, a large number of DEGs involved in sex differentiation, gonadal development and differentiation and gametogenesis were identified, and the differential expression patterns of some genes were further verified using real-time fluorescence quantitative PCR. The results of this study will provide a valuable resource for further studies on sex determination and gonadal development of S. barcoo.


Introduction
Jade perch (Scortum barcoo), also known as barcoo grunter, belongs to the genus Scortum, the family Theraponidae and the order Perciformis [1]. It is indigenous to Australia and normally found in the Lake Eyre Basin and Barcoo River of Australia [2]. Jade perch has a high meat yield and rich nutritional value. Currently, the main sources of omega-3 are deep-sea fish and algae, which are difficult to obtain. Thus, jade perch is believed to be an alternative source of omega-3 due to its high level of omega-3 [3]. Moreover, jade perch grow rapidly and can reach commercial specification in 6-10 month under artificial aquaculture conditions [4]. Jade perch adapt well to high-density aquaculture and are suitable for cage culture and factory recirculating aquaculture [5]. Since being introduced into China in 2001, it has become an extremely important and valuable fish in China [6]. Under natural conditions, jade perch can reach sexual maturity at 3-4 years old. In the non-breeding season, the individual differences between mature male and female fish are small. During the breeding season, the belly of the jade perch accumulates a lot of fat, inducing abdominal enlargement, so it is also difficult to distinguish its sex from the external morphology. The lack of sexual dimorphism in body size and the inability to discriminate between males and females based on a single trait of jade perch greatly limit the selection for reproductive production, sex control and genetic improvement [7,8], which limits the yield increase of jade perch. Thus, a better understanding of the sex-related genes of this species should be developed.
Stephan first proposed the universal phenomenon of sex transition is a common phenomenon in fishes [9]. Later, researchers found sex transition in Monopterus javanensis and Sparidae fishes [10,11], which started the research on sex determination and sex transition in fishes. The sex of fish is regulated by both genetic and environmental factors [12]. Unlike environmental factors, genetic factors determine the initial sex of fish [13]. After sex determination, the development process of sex-related traits also involves more complex gene interactions in which sex-differentiation genes play a key role [14]. At present, many sexdetermining genes and sex-determination pathways have been identified in fish, including dmy, sdY, amhr2, amhy, gsdf, dmrt1, the Wnt signaling pathway, the transforming growth factor-β (TGF-β) gene family and the estrogen signaling pathway [15,16]. In the process of aquaculture, it is of great significance to study the mechanism of sex control in fish. First, it is beneficial to achieve the artificial control of fish sex at the biological molecular level in production, which is crucial for reproductive management and breeding practices. Secondly, it is of great significance to clarify the theories of sex determination and differentiation in fish. Up to now, previous studies mainly focused on breeding, morphology, rearing conditions, feeding efficiencies and immune responses of jade perch [4,5,7,[17][18][19][20], whereas less attention was paid to sex determination and differentiation. Therefore, it is necessary to explore functional genes related to sex determination and differentiation in jade perch in order to better understand the mechanisms of sex regulation and effectively guide jade perch breeding for sex control.
In this study, the first comparative transcriptome analysis of jade perch gonads was generated using Illumina-based transcriptome sequencing. The transcriptome comparative analysis results revealed a great number of sex-related genes, which were further identified and discussed following quantitative real-time PCR (qRT-PCR). Our study aims to fill the blanks in jade perch gonadal transcriptome data, identify candidate sex-related genes and provide a theoretical basis for the subsequent exploration of jade perch sexdetermining genes.

Sample Collection
In our study, 6 seven-month-old jade perch individuals (3 females and 3 males), the gender of which was identified by the morphological observation of the gonads, were obtained from Shaoguan, Guangdong province. The fishes were anesthetized with 300 mg/L MS222 (Sigma, Saint Louis, MO, USA). Gonads from each individual were excised within 60 s of euthanasia. The gonads to be used for RNA extraction were stored in liquid nitrogen, and those intended for gonadal histology analysis were fixed in Bouin's solution for 24 h at room temperature and then transferred to 75% ethanol. All animal experiments were conducted in accordance with the guidelines and approval of the respective Animal Research and Ethics Committees of Sun Yat-Sen University.

RNA Extraction and Library Construction
The total RNA was extracted using an RNA-isolating Total RNA Extraction Reagent (Vazyme, Nanjing, China), following the manufacturer's instructions. The quantity and quality of the RNA were examined using Nanodrop 2000 (Thermo Scientific, Wilmington, DE, USA) and 1.0% agarose gel electrophoresis, respectively. Then, the integrity and quantity of the RNA were checked with a precise concentration using an Agilent 4200 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). RNA with qualified quality was entered into the library construction process. The jade perch transcriptome libraries were constructed using the NEBNext ® Ultra™ RNA Library Prep Kit for Illumina ® (NEB, Ipswich, MA, USA) kit. The mRNA of jade perch was enriched using magnetic beads with Oligo (dT). Under a high-temperature environment and the presence of metal ions, the RNA was fragmented, and the first cDNA strand was synthesized by random hexamers. Then the second cDNA strand was synthesized by adding an enzyme, a buffer solution and dNTPs (dATP, dTTP, dGTP and dCTP). The synthesized double-stranded cDNA was purified using magnetic beads, followed by end repair and the addition of A. The ligation of sequencing adapters and fragment size sorting was performed using magnetic sorting beads. The sorted fragments were enriched by PCR, and the final PCR products were purified to obtain the final library.

Library Sequencing, De Novo Assembly and Annotation
The qualified libraries were sequenced using the Illumina Novaseq6000 (Illumina, Inc., San Diego, CA, USA) high-throughput sequencing platform with a sequencing strategy of PE150 (Pair-End 150), and the amount of sequencing data for each gonadal sample was not less than 6 Gb.
The results of Illumina high-throughput sequencing were converted to raw reads after base calling using CASAVA (v1.8.2) software. In order to filter the splice sequences and low-quality sequences in the raw data and ensure the quality of analysis data, Fastp v0.20.1 software [31] was used to filter the raw reads, and high-quality clean reads were obtained. The parameters for filtering data in Fastp v0.20.1 software were set to "-l 150 -q 20", and the rest were default parameters. All subsequent analyses were based on clean reads. Then, the de novo assembly was carried out by Trinity (v2.4.0) software [32]. The integrity of sequence assembly was assessed by BUSCO v5.2.2 [33].
Three forward and 3 reverse reading frames were used to predict the coding region of the unigene, and then 6 possible coding protein sequences were predicted. The obtained coding protein sequences were compared with non-redundant protein sequences (Nr) (http s://www.ncbi.nlm.nih.gov/, accessed on 20 July 2022) and the UniProt protein database (https://www.uniprot.org/, accessed on 25 July 2022). Finally, the encoding method with the best matching (the maximum score) was selected as the encoding method for this gene.
Based on homology searches, unigenes were annotated against major public databases, including the Clusters of Eukaryotic Orthologous Group database (KOG; http://www. ncbi.nlm.nih.gov/KOG, accessed on 28 June 2022), the Nr database, the UniProt protein database and the Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.geno me.jp/kegg, accessed on 18 July 2022) using BLASTP or Diamond. Blast2GO was used to further annotate the Nr results [34]. The Gene Ontology (GO; http://geneontology.org, accessed on 20 July 2022) analysis was performed according to the function of the protein sequence. The GO enrichment bar chart could visually reflect the information of genes in 3 aspects: biological process, cellular component and molecular function.

Identification of DEGs and Enrichment Analysis
Using HISAT2 v2.1.0 [35], the clean reads of each gonadal sample are first aligned to the non-redundant assembled transcripts. According to the alignment results, StringTie v1.3.3b [36] was used to further calculate the transcripts per million (TPM), fragments per kilobase of exon model per million mapped fragments (FPKM), coverage and other gene-expression-related values of each transcript. Python prepDE.py was used to convert the output results of StringTie into the "edgeR" package (V3.6) recognition format [37,38], and edgeR was used for gene differential expression analysis. "p value < 0.05, |log2FC|>2" was used for screening condition to screen significantly DEGs. The number of DEGs was counted by statistical significance, and then against the background of DEGs, the clusterProfiler program of the R package [39] was applied based on Fisher's exact test and Benjamini correction. The pathways with p-values < 0.05 were significantly enriched in DEGs compared with the whole genome.

Validation of DEGs Using Quantitative Real-Time PCR (qRT-PCR)
To validate the reliability of DEGs from RNA-seq, a total of 18 putatively sex-related DEGs were selected for quantitative real-time PCR (qRT-PCR) analysis. The primers of the DEGs and the reference gene β-actin were designed using Primer 5.0 software ( Table 1). The cDNA templates were synthesized using the HiScript III RT SuperMix for qPCR (+gDNA wiper; Vazyme, Nanjing, China). The qRT-PCR was conducted using the SYBR Green qPCR Mix (GDSBIO, Guangzhou, China) and performed on a Roche LightCycler 480 real-time PCR system (Roche, Basel, Switzerland) as per the following process: pre-denaturation at 95 • C for 3 min followed by 45 cycles at 95 • C for 10 s, 58 • C for 20 s and 72 • C for 15 s and a final extension performed at 72 • C for 5 min, ending with a dissociation curve process. Each gene was tested for 3 independent biological replicates and 3 technical replicates, and the expression of each gene was normalized using the β-actin by 2 −∆∆CT method [40].

Gonadal Histology
Firstly, the gonads were successively dehydrated in 85%, 95% and 100% ethanol. After dehydration, the gonads were placed in xylene to clear and embedded in paraffin for 3 h. The embedded tissue blocks were serially sectioned at a thickness of 5-7 µm. After drying at 42 • C overnight, the slices were stained with hematoxylin and eosin. The clear slides were quickly sealed with neutral resin and air-dried overnight at room temperature. Stained gonads sections were observed and photographed using a light microscope (Nikon Eclipse Ni-U, Shinagawa-ku, Japan).

Overview of Illumina Sequencing and Assembly Results
A total of six cDNA libraries, including three testes and three ovaries, were constructed by RNA-seq. After quality control and filtering, a total of 41.18 Gb of clean reads were generated using the Illumina HiSeq platform, with an average of 6.86 Gb and a range of 5.78-7.92 Gb per sample. The mean value of Q20 was 97.71%, and that of Q30 was 93.46% (Table 2). After being assembled de novo, a sum of 107,060 unigenes was assembled, and the total length of these unigenes was 107,070,530 bp. The sequence length distribution of the unigenes ranged from 185 to 18,090 bp, with an average length and N50 of 1000 bp and 2336 bp, respectively (Table 3). Based on the predicted length statistics of unigenes, a total of 81,664 (76.28%) unigenes were 1-1000 bp in length, and 3202 (2.99%) of them were more than 5000 bp in length ( Figure 1).

Unigenes Annotation
To further study the specific molecular processes in jade perch, all of the predictable unigenes were annotated in different databases by significant similarity to protein sequences. Then, a total of 28,514 unigenes were annotated successfully. Among all the annotated unigenes, 28,165 (98.78%) unigenes were matched with the Uniport database, which had the highest percentage of annotated unigenes. However, this was only 15,592 (54.68%) unigenes compared to the KOG database, which may be related to the genomic information available limitation for jade perch (Table 3). Species distribution statistics were performed according to the most similar genes compared to the Nr database. The species with the largest number of genes indicated that the species contained the largest number of genes in the Nr database and was relatively similar to that of jade perch. The results showed that Siniperca chuatsi (4985, 17.87%), Chelmon rostratus (1836, 6.58%), Morone saxatilis (1454, 5.21%), Larimichthys crocea (1372, 4.92%) and Liparis tanakae (1024, 3.67%; Figure 2).
In addition, the GO, KOG and KEGG databases were used for functional prediction and classification of all unigenes. In the GO database, a total of 71,984 (67.24%) unigenes were annotated and classified into three functional categories. Among these three functional groups, the terms "cellular process" (36.15%), "cellular anatomical entity" (88.10%) and "binding" (45.87%) were dominant in the biological process, cellular

Unigenes Annotation
To further study the specific molecular processes in jade perch, all of the predictable unigenes were annotated in different databases by significant similarity to protein sequences. Then, a total of 28,514 unigenes were annotated successfully. Among all the annotated unigenes, 28,165 (98.78%) unigenes were matched with the Uniport database, which had the highest percentage of annotated unigenes. However, this was only 15,592 (54.68%) unigenes compared to the KOG database, which may be related to the genomic information available limitation for jade perch (Table 3). Species distribution statistics were performed according to the most similar genes compared to the Nr database. The species with the largest number of genes indicated that the species contained the largest number of genes in the Nr database and was relatively similar to that of jade perch. The results showed that Siniperca chuatsi (4985, 17.87%), Chelmon rostratus (1836, 6.58%), Morone saxatilis (1454, 5.21%), Larimichthys crocea (1372, 4.92%) and Liparis tanakae (1024, 3.67%; Figure 2).
In addition, the GO, KOG and KEGG databases were used for functional prediction and classification of all unigenes. In the GO database, a total of 71,984 (67.24%) unigenes were annotated and classified into three functional categories. Among these three functional groups, the terms "cellular process" (36.15%), "cellular anatomical entity" (88.10%) and "binding" (45.87%) were dominant in the biological process, cellular component and molecular function aspects, respectively ( Figure 3A). Based on the KEGG database, a sum of 18,855 unigenes was clustered into five different functional categories. The top three distributions were "signal transduction" (2755 unigenes), "global and overview maps" (1742 unigenes) and "endocrine system" (1398 unigenes; Figure 3B). After KOG annotation, a total of 13,725 unigenes were grouped into 25 classifications, and the largest distribution was "general function prediction only" (2835 unigenes), followed by "signal transduction mechanisms" (2591 unigenes). component and molecular function aspects, respectively ( Figure 3A). Based on the KEGG database, a sum of 18,855 unigenes was clustered into five different functional categories. The top three distributions were "signal transduction" (2755 unigenes), "global and overview maps" (1742 unigenes) and "endocrine system" (1398 unigenes; Figure 3B). After KOG annotation, a total of 13,725 unigenes were grouped into 25 classifications, and the largest distribution was "general function prediction only" (2835 unigenes), followed by "signal transduction mechanisms" (2591 unigenes).

Differential Genes Expression Analysis
The TPM values were used to normalize the levels of gene expression. A total of 23,849 DEGs were identified between the testes and ovaries, with 14,332 (60.09%) DEGs significantly high expressed in testes and 9517 (39.91%) high expressed in ovaries ( Table  4). The DEGs information is shown in the volcano plot ( Figure 4). The enrichment analysis of these DEGs was further performed using GO and KEGG pathway analyses. The results

Differential Genes Expression Analysis
The TPM values were used to normalize the levels of gene expression. A total of 23,849 DEGs were identified between the testes and ovaries, with 14,332 (60.09%) DEGs significantly high expressed in testes and 9517 (39.91%) high expressed in ovaries ( Table 4). The DEGs information is shown in the volcano plot ( Figure 4). The enrichment analysis of these DEGs was further performed using GO and KEGG pathway analyses. The results of the KEGG enrichment analysis revealed that the neuroactive ligand-receptor interaction was the most representative, followed by the calcium signaling pathway and the cytokinecytokine receptor interaction ( Figure 5).    . KEGG enrichment analysis of the top 20 pathways. X-axis: the ratio of the number of differential genes annotated in the KEGG pathways to the whole number of differential genes; Yaxis: the name of the enriched KEGG pathways. The sizes of the dots represent the count of the differential genes, and the color of the dots represents the significant enrichment from red to blue.  Figure 5. KEGG enrichment analysis of the top 20 pathways. X-axis: the ratio of the number of differential genes annotated in the KEGG pathways to the whole number of differential genes; Y-axis: the name of the enriched KEGG pathways. The sizes of the dots represent the count of the differential genes, and the color of the dots represents the significant enrichment from red to blue.

Validation of Transcriptomic Results by qRT-PCR
A sum of 10 testes-biased and eight ovaries-biased genes was selected to validate the transcriptomic expression profiles by qRT-PCR analysis. As expected, the expression patterns of qRT-PCR were consistent with the RNA-seq results ( Figure 6). DEGs such as hsd11b1, gsdf, sf1, dmrtb1, dmrt3, nanos2, amh, sox9, dscam and cyp11a1 were ovariesupregulated ( Figure 6A), whereas unigenes such as figla, hsd17b1, dmrt2b, bmp15, nanos3, foxl2, foxh1 and gmnn were testes-upregulated ( Figure 6B). To validate the reliability and accuracy of DEG transcriptomic expression levels, the consistent tendencies of expression levels between the RNA-Seq and qRT-PCR data were evaluated using a correlation analysis ( Figure 6C).

Gonads Histology Analysis
The testes and ovaries of 7 month old jade perch were used for the histological analysis. The testes of jade perch were full of primary spermatogonia, which did not show proliferation or differentiation, and no spermatocyte was found ( Figure 7A). Meanwhile, primary growth oocyte appeared in the ovary ( Figure 7B).

Gonads Histology Analysis
The testes and ovaries of 7 month old jade perch were used for the histological analysis. The testes of jade perch were full of primary spermatogonia, which did not show proliferation or differentiation, and no spermatocyte was found ( Figure 7A). Meanwhile, primary growth oocyte appeared in the ovary ( Figure 7B).

Discussion
It is believed that sex determination is the most paramount for successful reproduction in animals. During the process of sex determination, a set of sexual functional genes is involved in many complex biological processes to support the gonads' differentiation into ovaries or testes. To figure out the complicated regulatory mechanism of sex determination, transcriptome sequencing has been demonstrated in more and more species and is considered a highly effective method. Jade perch, which has a high meat yield and rich nutritional value, is a significant commercial aquaculture species. Up to now, the gene expression profile and the molecular mechanisms during sex determination and gonad differentiation have not been revealed in jade perch. In our study, transcriptome analysis and comparative analysis have been conducted between 7 month old male and female juveniles and the expression patterns of DEGs in the differentiating gonads were established for jade perch.
As a key enzyme in the steroid hormone synthesis pathway, StAR is responsible for transporting cholesterol, the initial substrate for steroid synthesis, from the cytoplasm to the mitochondria [45]. Changes in the expression of StAR are able to modulate various signals of steroidogenesis, enabling a rapid response by steroidogenic tissues. In fish, StAR is thought to play a role in testicular development, spermatogenesis and spermatogenesis by regulating androgen production. The expression of the cyp11b2 gene, serum testosterone and 11-KT levels in StAR2 knockout male tilapia were significantly lower than those in control fish [46]. Mutant StAR1 in zebrafish significantly increases testosterone and estrogen levels and inhibits oocyte maturation [47]. In jade perch, the expression of the StAR1 and StAR2 genes was higher in the testes, which agreed with the expression profile of spotted scat and bastard halibut (Paralichthys olivaceus) [24,48].

Discussion
It is believed that sex determination is the most paramount for successful reproduction in animals. During the process of sex determination, a set of sexual functional genes is involved in many complex biological processes to support the gonads' differentiation into ovaries or testes. To figure out the complicated regulatory mechanism of sex determination, transcriptome sequencing has been demonstrated in more and more species and is considered a highly effective method. Jade perch, which has a high meat yield and rich nutritional value, is a significant commercial aquaculture species. Up to now, the gene expression profile and the molecular mechanisms during sex determination and gonad differentiation have not been revealed in jade perch. In our study, transcriptome analysis and comparative analysis have been conducted between 7 month old male and female juveniles and the expression patterns of DEGs in the differentiating gonads were established for jade perch.

DEGs Involved in Steroidogenesis Pathway
In fish, sex steroid hormones, especially 17β-estradiol (E 2 ) and 11-ketotestosterone , are involved in sex determination, sex differentiation and sex maintenance [41,42]. Due to the different expression profiles of genes that encode steroid-metabolizing enzymes, the serum levels of E 2 and 11-KT have sexual dimorphism between females and males [24,43]. At present, cyp19a1a, cyp11b2, cyp11a2, hsd17b and hsd11b1 are considered essential steroidogenesis genes [44].
As a key enzyme in the steroid hormone synthesis pathway, StAR is responsible for transporting cholesterol, the initial substrate for steroid synthesis, from the cytoplasm to the mitochondria [45]. Changes in the expression of StAR are able to modulate various signals of steroidogenesis, enabling a rapid response by steroidogenic tissues. In fish, StAR is thought to play a role in testicular development, spermatogenesis and spermatogenesis by regulating androgen production. The expression of the cyp11b2 gene, serum testosterone and 11-KT levels in StAR2 knockout male tilapia were significantly lower than those in control fish [46]. Mutant StAR1 in zebrafish significantly increases testosterone and estrogen levels and inhibits oocyte maturation [47]. In jade perch, the expression of the StAR1 and StAR2 genes was higher in the testes, which agreed with the expression profile of spotted scat and bastard halibut (Paralichthys olivaceus) [24,48].
In male fish, androstenedione is catalyzed by 17β-HSD (hsd17b3), P450c11 (cyp11b2) and 11β-HSD (hsd11b2) to form 11-KT, while in female fish, androstenedione is converted to E 2 by P450c19 (cyp19a1a/b) [41,49,50]. As a crucial gene encoding key rate-limiting enzymes that catalyze the synthesis of fish-specific 11-KT, cyp11b2 plays an important role in the development of the testes. Cyp11b2 was found to be highly expressed during spermatogenesis in catfish (Clarias gariepinus), rainbow trout (Oncorhynchus mykiss), Japanese eel (Anguilla japonica), tilapia, Atlantic salmon (Salmo salar) and medaka (Oryzias latipes), but hardly detectable at all stages of ovarian development, suggesting that cyp11b2 plays an important role in spermatogenesis [51][52][53][54][55][56][57]. Moreover, in the development of fish ovaries, the cyp19a1a gene, which can directly catalyze the switching of androstenedione to E 2 , plays an essential role in the process of primitive gonads differentiating into ovaries [58,59]. After the expression of the mutant cyp19a1a gene in the female medaka, the gonads still develop into ovaries in the early stage, but a partial female-to-male sex reversal occurs in the later stage [60]. In zebrafish and Nile tilapia, mutant cyp19a1a results in female-to-male sex reversal [61][62][63][64]. Similarly, the expression of cyp11b2 in males was higher than that in females in jade perch, whereas the cyp19a1a gene was highly expressed in the ovaries. Therefore, the results of these genes indicated that these DEGs were crucial in the steroid hormone synthesis pathway and gonad development in fish.

DEGs Involved in Gonad Differentiation and Development
Sex determination is influenced by genetics, the environment or both. The development of testes and ovaries is crucial for gametogenesis and reproduction. Therefore, understanding the molecular mechanisms of sex determination is an important direction in reproductive biology. Sex determination and gonadal differentiation in fish are diverse. In recent years, nearly 20 different genes have been confirmed to be involved in teleost sex determination and gonadal differentiation, including the Dmrt gene family, the Sox gene family and the transforming growth factor-β (TGF-β) gene family [15,16].
The Dmrt1 gene is specifically expressed in the male gonads of semi-smooth tongue sole at the critical stage of sex determination, and its expression level is significantly upregulated in phenotypic males after female sex reversal [65]. Dmrt1 deficient Thai betta fish (Betta splendens) undergo sex reversal from male to female [66]. Steroidogenic factor-1 (sf1) is a key gene in vertebrate sex determination. It is expressed in the adrenal gland, gonad, spleen, ventromedial hypothalamus and pituitary gonadotrophins and regulates a variety of steroidogenic enzymes, including aromatase [67]. Protandrous black porgy fish (Acanthopagrus schlegeli) has a high expression of sf1 in the testis and a low expression in the ovary before sexual maturity, and it has an antagonistic effect on the development of oocytes [68]. Transcription factor binding sites for dmrt1 and sf1 were revealed using a computer-simulated promoter analysis 5 upstream of the spotted scat gsdf gene. Results showed that dmrt1 activated gsdf expression in spotted scat in a dose-dependent manner in the presence of sf1 [69]. In jade perch, the expression of dmrt1, dmrt3 and sf1 were higher in testes, which suggests a key role in the testes' development.
In recent years, members of the TGF-β signaling pathway have also been confirmed by more and more researchers to be sex-determining genes in vertebrates, especially in fish [70]. Mullerian ducts are absent in fish, but the anti-Mullerian hormone (amh) gene is widely present in fish and participates in fish germ cell proliferation. In zebrafish, mutations in the amh gene have been found to increase the proportion of females [71]. In addition, the amh copy gene amhy has been confirmed to be the sex-determining gene in Patagonian pejerrey (Odontesthes hatcheri), Nile tilapia and Northern pike (Esox lucius) by gene editing knockout and transgenic methods [72][73][74]. Myosho et al. found a specific copy of gsdf on the Y chromosome in Oryzias luzonensis, and the up-regulation of this gene in the reproductive primordium of the early embryo would lead to the differentiation of gonads into males [75].
Through genome sequencing analysis, it was found in sablefish (Anoplopoma fimbria) that there were specific sequences of different lengths in the upstream region of the gsdf gene on the X and Y chromosomes [76]. Similarly, A specific copy of the gsdf gene was found on the Y chromosome of Atlantic halibut (Hippoglossus hippoglossus) [77], suggesting that gsdf and its adjacent upstream regulatory elements are responsible for sex determination. Growth and differentiation factor 9 (gdf9) and bone morphogenetic protein 15 (bmp15) are involved in ovarian development in bony fish by regulating the expression of gonadotropin and its receptors. In European sea bass (Dicentrarchus labrax) and Japanese flounder, changes in the expression levels of gdf9 and bmp15 affect the expression of follicle-stimulating hormone receptors (fshr) and luteinuclein receptors (lhr) [78,79]. Overexpression of bmp15 in cultured zebrafish eggs could inhibit gonadotropin-induced ovarian maturation [80]. The expression pattern of these genes in jade perch was in accord with the roles of sex determination, i.e., the mRNA levels of gsdf and amh were higher in males, while the bmp15 and gdf9 genes had higher expression levels in females.
The nanos2 gene is specifically expressed in male germ cells [81]. The Nanos2 gene is specifically expressed in tilapia testes but not expressed in eggs and ovaries [82]. The nanos3 gene regulates the migration of primordial germ cells and plays an important role in the production and maintenance of oocytes and the development and differentiation of germ cells [83]. Previous studies have confirmed that the nanos3 gene is highly expressed in the ovary of zebrafish (Danio rerio) and large yellow croaker (Larimichthys crocea) and involved in the development and maintenance of germline stem cells (GSCs), whereas it is expressed to a lesser extent in the testes [84,85]. In jade perch, nanos2 expression was higher in males, and nanos3 expression was higher in females, which suggests its important role in germ cell development and differentiation.

DEGs Involved in Gametogenesis and Gamete Maturation
In the praxis of fish culture, the ultimate goal of breeding farmed species is to promote the development of testes and ovaries and then to obtain healthy mature gametes for the generation of offspring. In recent years, genes including zona pellucida nuclear receptor subfamily 0 group B member 1 (nr0b1), forkhead box (Fox) and sperm-binding proteins (zps) have been identified as associated genes of oogenesis, oocyte maturation, spermatogenesis and sperm motility [21,29].
Dax1 (DSS-AHC critical region on the X chromosome 1, gene 1), which belongs to the nuclear receptor family, is a special orphan protein encoded by the nr0b1 gene [86,87]. Dax1 has been identified as an important factor in the regulation of embryonic cell pluripotency [88]. The expression pattern of dax1 varies among different fish species. It is found that dax1 begins expression in the bipotential gonads of zebrafish 13 days after hatching and has significant expression in the ovary after the morphological differentiation of the gonads. Moreover, in dax1 −/− zebrafish, oocyte proliferation is reduced, resulting in arrested ovarian development and sex reversal to testes [89]. In medaka, dax1 is highly expressed in the ovary of adult fish but is hardly expressed in the testes [90]. On the contrary, in tilapia, the dax1 gene is continuously expressed from the critical period of sex determination and differentiation (5 days post-hatching) into adult fish, and there is no significant difference in the expression level between female and male gonads [91]. Shi et al. proposed that dax1 is involved in ovarian development in spotted scat by regulating Nr5A1-mediated cyp19a1a expression [92]. According to the transcriptome data in this study, the nr0b1 gene was more highly expressed in jade perch testes, which showed a different pattern.
The Fox family is a class of transcriptional regulatory nuclear proteins, and members of the Fox family are involved in the regulation of multiple biological processes. Researchers indicated that fox genes might play important roles in sex determination and gonadal development in teleosts [93]. In zebrafish, foxh1 expression is higher in females than in males, and zebrafish eggs with foxh1 gene deletion develop abnormally [94,95]. However, in tilapia, foxh1 is expressed in the cytoplasm of oocytes in the ovary [93], and female fishes with foxh1 knockout are arrested in follicular development, which eventually leads to female infertility [96]. In teleost fish, zona pellucid (ZP) protein provides mechanical protection to the oocyte during fertilization, mainly structurally [97,98]. Zp3, a member of the zp superfamily, is a major female-specific gene that forms the extracellular matrix around the oocyte and mediates sperm binding. The Zp3 gene is highly expressed in the ovaries of Spot-fin porcupine fish and Spinibarbus hollandi, and the researchers hypothesized that zp3 genes may be involved in ovarian folliculogenesis [21,29]. Herein, the foxh1 and zp3 genes had higher expression levels in female jade perch, suggesting the important roles of these genes in oogenesis.

Conclusions
This work is the first research on the gonad transcriptome of juvenile jade perch. Herein, we assembled 107,060 unigenes based on gonad transcriptome data. A high number of DEGs, which are mainly involved in the steroidogenesis pathway, gonad development and differentiation, gametogenesis and gamete maturation, were identified by comparative transcriptome analysis. For the proven DEGs associated with gonadal development and reproduction, we found similar expression profiles in jade perch, suggesting that they play a conserved role in gonadal development and gametogenesis. In addition, we also identified and analyzed the expression profiles of other genes that may be involved in gonad development and gametogenesis. The results of this study will provide valuable information for further studies on sex determination and gonadal development in teleost fish.