Journal of Oceanology and Limnology   2021, Vol. 39 issue(4): 1500-1512     PDF       
http://dx.doi.org/10.1007/s00343-020-0251-y
Institute of Oceanology, Chinese Academy of Sciences
0

Article Information

ZHU Qiaoying, HAN Chong, PENG Cheng, ZHOU Xingni, WANG Chongwei, HAN Linqiang, LI Shuisheng, LI Guifeng, LIN Haoran, ZHANG Yong
Identification of potential sex-related genes in Siniperca chuatsi
Journal of Oceanology and Limnology, 39(4): 1500-1512
http://dx.doi.org/10.1007/s00343-020-0251-y

Article History

Received Jul. 17, 2020
accepted in principle Aug. 31, 2020
accepted for publication Oct. 13, 2021
Identification of potential sex-related genes in Siniperca chuatsi
Qiaoying ZHU1,2, Chong HAN1,2, Cheng PENG3, Xingni ZHOU1,2, Chongwei WANG1,2, Linqiang HAN4, Shuisheng LI1,2, Guifeng LI1, Haoran LIN1,2, Yong ZHANG1,2     
1 State Key Laboratory of Biocontrol, Guangdong Provincial Key Laboratory for Aquatic Economic Animals and Southern Marine Science and Engineering Guangdong Laboratory(Zhuhai), School of Life Sciences, Sun Yat-Sen University, Guangzhou 510275, China;
2 Laboratory for Marine Fisheries Science and Food Production Processes, Qingdao National Laboratory for Marine Science and Technology, Qingdao 266373, China;
3 Guangdong Key Laboratory of Animal Conservation and Resource Utilization, Guangdong Public Laboratory of Wild Animal Conservation and Utilization, Guangdong Institute of Applied Biological Resources, Guangzhou 510260, China;
4 Guangdong Liangshi Aquatic Seed Industry Co., Ltd., Foshan 528100, China
Abstract: Mandarin fish (Siniperca chuatsi) is an economically important freshwater fish cultured in China. In this species, females grow faster than males. However, due to the lack of available genomic and transcriptome information, the mechanisms of sex differentiation remain poorly understood. In this study, Illumina high-throughput sequencing technology was used to sequence four cDNA libraries; the tissues examined included hypothalamus, pituitary gland, ovary, and testis. A total of 134 124 high-quality unigenes were obtained on average length of 1 361 bp and N50 of 3 312 bp. A search of all-unigene sequences against NR, NT, SwissProt, KOG, KEGG, GO, and InterPro databases resulted in 59 688 (44.50%), 76 329 (56.91%), 50 432 (37.60%), 45 741 (34.10%), 48 760 (36.35%), 5 241 (3.91%), and 46 099 (34.37%) annotations, respectively. In a comparison of ovarian and testicular libraries, 15 289 ovary-biased genes and 10 035 testis-biased genes were identified, including a series of genes related to sex determination and differentiation, such as cyp19a1a, foxl2, sox9, dmrt1, amh, and others. In addition, 49 495 SSRs and 85 899 SNPs were detected in transcriptome data. Quantitative real-time PCR results of 15 sex-related functional genes indicated that RNA-seq data was reliable. This study will contribute to a better understanding of the molecular mechanisms of sex differentiation and development in Mandarin fish.
Keywords: transcriptome    Mandarin fish    HPG axis    female and male    
1 INTRODUCTION

In almost all vertebrates, sexual reproduction is an important biological process, essential to maintain long-term survival and procreation. In order to ensure the continuation of species, vertebrates have evolved complex patterns of sex determination and differentiation. Fish are the most abundant vertebrates, distributed nearly in all aquatic environments, which contain more than 30 000 species (Froese and Pauly, 2017). Fish have evolved all the strategies for sexual reproduction known in vertebrates, including parthenogenesis (Schartl et al., 1995), monoecism (Warner, 1984) and dioecism. The mechanisms of fish sex determination include two forms: genetic sex determination (GSD) and environmental sex determination (ESD) (Janzen, 1995). However, sex determination and differentiation remain poorly understood as they are extremely plastic and changeable, and even fish with GSD can be easily affected by environmental factors and exogenous hormones (Ospina-Álvarez and Piferrer, 2008). Currently, only a few master sex determination genes have been confirmed in several fish species: dmy in medaka (Oryzias latipes) (Matsuda et al., 2002; Nanda et al., 2002), gsdf in medaka (Oryzias luzonensis) (Myosho et al., 2012), irf9y in rainbow trout (Oncorhynchus mykiss) (Yano et al., 2012), amhr2 in fugu (Takifugu rubripes) (Kamiya et al., 2012), dmrt1 in half smooth tongue sole (Cynogossus semilaevis) (Chen et al., 2014) and amhy in Patagonian silverside (Odontesthes hatcheri) (Hattori et al., 2012), and Nile tilapia (Oreochromis niloticus) (Li et al., 2015).

With the fast development of next-generation sequencing technologies, transcriptome sequencing is increasingly used to acquire gene expression data and to understand the associated mechanisms of regulation. Gonad transcriptome analysis has been used to identify genes related to sex differentiation and determination in many fish such as Nile tilapia (Oreochromis niloticus) (Tao et al., 2013), Japanese flounder (Paralichthys olivaceus) (Fan et al., 2014), Russian sturgeon (Acipenser gueldenstaedtii) (Hagihara et al., 2014), spotted knifejaw (Oplegnathus punctatus) (Du et al., 2017), fugu (Takifugu tubripes) (Wang et al., 2017b) and channel catfish (Ictalurus punctatus) (Sun et al., 2013). These data provided a lot of gonadal transcriptomic information and identified numerous sex-related genes, facilitating further studies on fish sex determination and gonadal differentiation. Moreover, previous studies have also demonstrated that the brain controls reproduction through the hypothalamic-pituitary-gonadal (HPG) axis and it influences gonad development (Weltzien et al., 2004; Sreenivasan et al., 2008). Thus, the hypothalamus and the pituitary gland are considered the tissues of choice to study the regulation mechanisms of sex determination and gonadal differentiation.

Mandarin fish, also known as Chinese perch (Siniperca chuasti), is one of the most important commercial fish, and has been widely cultured in China. It is in great demand from Chinese market because of its high palatability and high content of essential amino acids and unsaturated fatty acids (Chu et al., 2010; Zhang et al., 2011). Mandarin fish is a demersal piscivore that feeds solely on live prey fish (Liang et al., 1998) and, although it grows fast, it is susceptible to diseases. For this reason, previous studies mainly focused on its growth (Wang et al., 2016; Tu et al., 2017), immunity (Wang et al., 2012, 2014, 2017a) and prey preference (He et al., 2013, 2018). In addition, female Mandarin fish grow faster than the males, therefore culturing all-female populations is a very effective approach to boost aquaculture production. However, only a few sexrelated genes have been cloned, such as the aromatase P450 gene (cyp19a) (Zou et al., 2017) and various growth hormone genes (Lu et al., 2008). In this study, the whole HPG axis transcriptome of Mandarin fish was sequenced using an Illumina HiSeq 2000 platform and numerous differentially expressed genes were identified in ovary and testis samples. This study was designed to provide a comprehensive understanding of the reproductive axis in Mandarin fish and to identify sex-related genes participating in sex determination and in gonadal differentiation. These data provide a valuable genomic resource and a large number of molecular markers to be used for further research on sex determination and gonadal differentiation in Mandarin fish.

2 MATERIAL AND METHOD 2.1 Sampling

Samples of the hypothalamus, pituitary gland, testis, and ovary were collected from freshly caught Mandarin fish (weight>500 g), frozen in liquid nitrogen and stored at -80 ℃, to construct four independent libraries. The hypothalamus and pituitary gland samples of six Mandarin fish (3 males and 3 females) were mixed into a single sample. The testis and ovary samples from three different individuals were also mixed into one sample.

2.2 RNA extraction and library construction

Total RNA was extracted from fresh frozen tissue using a TRIzol Reagent kit (Invitrogen, CA, USA) following the extraction protocol. Total RNA concentration, RIN value, 28S/18S and fragment size were measured using an Agilent 2100 Bioanalyzer (Agilent RNA 6000 Nano Kit), while purity was determined by ultraviolet spectrophotometer (NanoDropTM).

After total RNA was extracted, mRNA with polyA tail was enriched using magnetic beads with Oligo (dT), and an appropriate amount of interrupting reagent was added to the mRNA to fragment it at high temperature. One strand of cDNA was synthesized using the interrupted mRNA as template, and then double stranded cDNA was synthesized. After kit purification, recycling, and adhesive end repairing, base "A" was added to the 3ʹ end of cDNA and the joint was connected. Following fragment size selection, PCR amplification was performed. The constructed library was tested by an Agilent 2100 Bioanalyzer and ABI Step One Plus Real-Time PCR System, and by real-time PCR. Finally, the cDNA library was sequenced using Illumina HiSeq™. All procedures were performed in Beijing Genomics Institute, Shenzhen, China.

2.3 De-novo assembly and functional annotation

To ensure the reliability of the results, reads with joint contamination, unknown base N content greater than 5%, and low quality were removed before data analysis, and the filtered reads were named "clean reads" and stored in FASTQ format. Then, a de novo transcriptome assembly was carried out using Trinity software.

Functional annotation was performed through homology searches against the main public databases. All unigenes were compared with the sequences in the NCBI nucleotide (NT, ftp://ftp.ncbi.nlm.nih.gov/blast/db) database using BLASTn, the NCBI nonredundant (NR, ftp://ftp.ncbi.nlm.nih.gov/blast/db) protein database, the Clusters of eukaryotic Orthologous Groups (KOG, http://www.ncbi.nlm.nih.gov/KOG) database, the Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg) and the SwissProt (http://ftp.ebi.ac.uk/pub/databases/swissprot) protein database using BLASTx or Diamond.

Blast2GO was used to identify the NR annotation results, and the number of unigenes related to each gene ontology (GO, http://geneontology.org) was calculated based on biological processes, cell composition and molecular function. InterProScan5 was used to annotate the InterPro (http://www.ebi.ac.uk/interpro) database.

2.4 Identification of differentially expressed genes (DEGs)

Based on the assembly results, clean reads of each sample were compared with each other using Bowtie2 software, and the gene expression level of each sample was calculated using RSEM. The differentially expressed genes (DEGs) between samples were detected using the PossionDis algorithm, which was defined as false discovery rate (FDR)≤0.001, and the gene with expression difference of more than 2 times. DEGs were classified and enriched for the GO function and KEGG pathway classification.

2.5 Data validation by qRT-PCR

The reliability of transcriptome data was verified by qRT-PCR using LightCycler 480 SYBR GreenⅠMaster (Roche) and Light cycler 480 RealTime PCR system (Roche). One denaturation cycle was performed at 95 ℃ for 10 min, and 40 amplification cycles were performed at 95 ℃ for 10 s, 60 ℃ for 20 s and 72 ℃ for 20 s. The relative expressions of 15 DEGs in testis and ovaries, including dmrt1, spata4, fshr, nanos2, sox9, cyp19a1a, foxl2, hsd17b1, hsd17b10, hsd17b12a, hsd17b8, spata2, spata5, star3, and star5, were conformed using beta-actin as the reference gene. All samples were triplicated and double-tailed student t-test was performed in confidence value of 95% (P≤0.05) to determine the significance of gene expression. All primers were designed using Primer Premier 6 (Table 1).

Table 1 Primers of qRT-PCR
2.6 Identification of SSRs and SNPs

Simple sequence repeat (SSR) sequences were identified using MISAv1.0 (http://pgrc.ipkgatersleben.de/misa) in order to search for single nucleotide to hexanucleotide repeats. Parameters were set as follows: single base repeat at least 12 times, double base repeat 6 times, triple base repeat 5 times, quadruple base repeat 4 times, penta base repeat 3 times, hexa base repeat 2 times. When the distance between two microsatellites is less than 100 bp, a composite microsatellite is formed. HISAT (http://ccb.jhu.edu/software/hisat/index.shtml) was used to compare clean reads to unigenes in order to obtain single nucleotide polymorphism (SNP) sequences, and GATK (https://www.broadinstitute.org/gatk) was used to filter out low-quality SNPs with FS filter set at >30.0 and QD filter set at < 2.0.

3 RESULT 3.1 Overview of transcriptome assembly quality

A total database of 27.56 Gb was analyzed using an Illumina HiSeq platform. Following the removal of assembly and redundancy, 134 124 unigenes were obtained, with total length, average length and GC content of 182 89 703 bp, 1 361 bp and 45.75%, respectively. The values of N50, N70, and N90 were 3 312 bp, 1 888 bp, and 451 bp respectively (Table 2). All unigenes were more than 300 bp in length, and 19 453 of them (14.50%) were more than 3 000 bp in length (Fig. 1).

Table 2 Overview of the overall quality of the transcriptome
Fig.1 Length distribution of all-unigene X-axis: length of unigene; Y-axis: corresponding number of unigene.
3.2 Gene function annotation

Through sequence matching against seven public databases, an annotation analysis was conducted on the remaining 134 124 non-redundant unigenes, and 82 290 sequences were annotated. The percentage of annotated unigenes in the NT database (76 329; 56.91%) was the highest, followed by the NR database (59 688; 44.50%), while the sequences annotated in the GO database were the lowest (5 241; 3.91%) (Table 3). A total of 37 028 genes were annotated in five public databases: NR, KOG, KEGG, SwissProt, and Inter Pro (Fig. 2a). The percentage of sequences consistent with all unigene BLASTx hits against other populations in the NR database showed that Larimichthys crocea (42.93%) has the largest amount of homologous sequences to Siniperca chuatsi, followed by Stegastes partitus (18.36%), Oreochromis niloticus (4.82%), and Notothenia coriiceps (3.72%) (Fig. 2b).

Table 3 Annotated gene distribution in seven public databases
Fig.2 Venn diagrams depicting functional annotation and species distribution of homologous genes annotated in the NR database

To predict the functional classification of genes, possible pathways, and gene interactions, all unigenes were annotated in GO, KOG, and KEGG databases.

A total of 18 247 unigenes were assigned to three GO categories: biological processes, cellular component and molecular function. In the category of biological processes, cellular (2 764) and metabolic processes (2 234) were the most represented items. In the cellular component category, the most represented terms were cell (2 100) and cell part (2 080). In addition, in the molecular function category, the most abundant terms related to binding (2 543) and catalytic activity (1 866) (Fig. 3a).

Fig.3 Gene function distribution annotated in GO (a), KOG (b), and KEGG (c) databases

A total of 45 741 (34.10%) KOG-annotated genes were classified into 25 molecular families, with the most abundant distribution observed in "signal transduction mechanisms" (12 936 unigenes), followed by "general function prediction only" (11 266 unigenes); the smallest distribution was "coenzyme transport and metabolism", with only 296 unigenes (Fig. 3b).

A total of 45 979 (27.92%) unigenes annotated in KEGG were divided into six different functional groups, among which the three exhibiting the largest distribution were "signal transduction" (9 131 unigenes), "global and overview maps" (4 820 unigenes) and "cancer: overview" (4 638 unigenes) (Fig. 3c).

3.3 Simple sequence repeats (SSRs) and single nucleotide polymorphisms (SNPs)

A total of 49 495 SSR markers were identified from 31 707 sequences. Most were mono-nucleotide (12 482 unigenes), di-nucleotide (23 987 unigenes) and tri-nucleotide repeats (11 085 unigenes). The AC/ GT motif was most abundant in dinucleotide sequences (17 494 unigenes), followed by the AG/CT motif (4 359 unigenes). The AGG/CCT motif was most abundant in trinucleotide sequences (3 317 unigenes), followed by the AGC/CTG motif (1 893 unigenes). In addition, there were 1 292 quadnucleotide, 372 penta-nucleotide, and 277 hexanucleotide sequences (Fig. 4a).

Fig.4 Distribution of SSRs and SNPs in the transcriptome of Siniperca chuatsi

In addition, a total of 85 899 SNPs were detected in transcriptome data. 24 426, 21 923, 14 414, and 25 136 SNP sites were detected from the four transcriptional libraries derived from hypothalamus, pituitary gland, ovary, and testis samples, respectively. The frequency of A-G transition was higher than C-T, and the frequency of A-T transversion was higher than other transversion frequencies (Fig. 4b; Table 4).

Table 4 Overview of SNP detection results
3.4 Differentially expressed genes in ovary and testis samples

A total of 25 324 DEGs were identified between the ovary and testis samples. Among them, 15 289 (60.37%) were significantly higher in the ovaries (ovary-biased genes) and 10 035 (39.62%) in the testes (testis-biased genes). The remaining 47 035 genes were expressed in both testes and ovaries. Enrichment analysis of the GO signaling pathway showed that the annotation of DEGs to biological processes was predominant, followed by cellular components and molecular functions (Supplementary Fig.S1). In biological processes, 13, 13, and 3 unigenes were annotated to reproduction, reproduction processes, and hormone secretion, respectively (Supplementary Table S1). Enrichment analysis of the KEGG signaling pathway showed that DEGs were mainly distributed in metabolic pathways, endocytosis, and other processes (Supplementary Fig. S2). In addition, different numbers of DEGs were annotated to the following pathways: 150 to GnRH, 259 to Wnt, and 130 to TGF-beta. Among the pathways related to steroids, estrogen signaling, ovarian teroidogenesis, steroid hormone biosynthesis, and steroid biosynthesis were 189, 77, 50, and 27, respectively (Supplementary Table S2).

Based on the comparison of the functional annotation of Mandarin fish transcriptome sequences from various databases with the published data of other species, numerous genes related to the HPG axis and sex differentiation were identified (Supplementary Table S3), such as kisspeptins (kiss), gonadotropin-releasing hormone (gnrh), gonadotropins (gths), forkhead transcription factor2 (foxl2), cytochrome P450 (cyp19a), doublesex and mab-3 related transcription factor 1 (dmrt1), SRY-box transcription factor 9 (sox9), anti-Mullerian hormone (amh).

3.5 Validation of transcriptomic data

Fifteen randomly selected DEGs between ovary and testis, including dmrt1, spata4, fshr, nanos2, sox9, cyp19a1a, foxl2, hsd17b1, hsd17b10, hsd17b12a, hsd17b8, spata2, spata5, star3, and star5, were used to verify the validation of transcriptomic data by qRT-PCR, and the results were consistent with the transcriptome data (Fig. 5).

Fig.5 Validation of gonad transcriptome data by qRT-PCR
4 DISCUSSION

Mandarin fish, commonly known as Chinese perch, displays obvious sexual growth dimorphism that females grow faster than males from juvenile to commercial size. This means that all-female Mandarin fish cultures have greater economic benefits and broader prospects (Wang et al., 2006). However, due to the lack of genomic and transcriptome data and studies on sex-related genes, the molecular mechanisms of sex differentiation in this species are poorly understood. In this study, 134 124 unigenes were obtained by Illumina high-throughput transcriptome sequencing with N50, N70 and N90 values of 3 312, 1 888, and 451, respectively, indicating that the transcriptome data were of high quality and reliable. A large number of sex-related genes in the HPG axis of Mandarin fish were found, providing potential candidate genes for future research on reproduction, development and sexual differentiation.

The HPG axis is the regulatory center of reproductive processes, and its signaling pathway is initiated by gonadotropin-releasing hormone (gnrh). This hormone displays dose-dependent induction by kisspeptin (kiss) (Novaira et al., 2009; Li et al., 2019b), and once it binds to the gnrh receptor (gnrhr), it stimulates the pituitary gland to secrete gonadotropins, including the follicle-stimulating hormone (fsh) and the luteinizing hormone (lh), to promote the production of sex hormones and the development of gametes. Due to the occurrence of fish-specific genome duplication (FSGD) events (Robinson-Rechavi and Laudet, 2001; RobinsonRechavi et al., 2001), most fish species have two kinds of kisspeptin, and this is also true for Mandarin fish (Selvaraj et al., 2010). On the contrary, only kiss2 was found in Solea senegalensis (Mechaly et al., 2011), Tetraodon nigroviridis and Gasterosteus aculeatus, suggesting that kiss2 may play a more important role. Previous studies confirmed that there are at least two or more gnrhs in teleost fish and that gnrh3 is unique to fish (Hildahl et al., 2011). As for gnrh receptors, the existence of two or more gnrhrs has been confirmed in teleost fish, and up to five gnrhrs have been identified in Fugu rubripes (Moncaut et al., 2005) and in Tetraodon nigroviridis (Ikemoto and Park, 2005). In the transcriptome of Mandarin fish, three types of gnrh and three types of gnrhr were found. The selectivity of GnRH receptors to the GnRH ligand varies in different species (Lethimonier et al., 2004), and the specific binding of the GnRH ligand to the receptor in Mandarin fish remains to be explored.

Gths can increase intracellular cAMP and promote the activation of the PKA subunit that regulates the expression of steroid hormone genes (Selstam et al., 1976; Reinhart et al., 1999; Stocco et al., 2005). It was suggested that fsh may play an important role in the cyp19a1a control during ovarian differentiation through the synthesis of cAMP second messenger (Yamaguchi et al., 2007; Guiguen et al., 2010). The major sex hormones active in teleost fish are testosterone (T), 11-ketotestosterone (11-kt) and 17-estradiol (E2). Their synthesis starts with cholesterol. Cholesterol is transferred from the outer mitochondrial membrane to the inner membrane under the action of a steroidogenic acute regulatory protein (star) (Miller, 2007), and then is converted to androstenedione under the catalysis of P450scc, CYP17A1 and 3β-HSD. Subsequently, 17β-HSD3 catalyzes the conversion of androstenedione to T, which is further catalyzed to 11-Ketotestosterone by CYP11 and 11β-HSD2, and CYP19 catalyzes the formation of E2 from androstenedione T (Simpson et al., 1994). Our data showed that in the process of androgen synthesis, most of the genes coding for speed-limit enzymes were highly expressed in males, including cyp11a1, cyp17a1, cyp11b, and hsd11β2. In contrast, during estrogen synthesis, CYP19 aromatase is the key rate-limiting enzyme, encoded by cyp19a1 (including cyp19a1a and cyp19a1b). In particular, cyp19a1a encodes gonadal aromatase, which is highly expressed in females and is mainly involved in estrogen production in Mandarin fish. The function of sex hormones is accomplished by binding to receptors. Several isoforms of er and ar have been found in Mandarin fish, including era, erb1, erb2, ara, and arb, which are similar to isoforms observed in other fish (Ogino et al., 2018).

Since the activity of CYP19 aromatase directly determines the male to female hormone ratio, cyp19a1a is considered the central factor for sex differentiation in fish (Li et al., 2019a). Most of the sex-related genes found in fish are related to cyp19a1a. Foxl2 is the earliest marker of ovarian determination and differentiation in vertebrates, and is highly expressed in female Mandarin fish. The transient transfection of fish showed that foxl2 not only directly activated the transcription of cyp19a1a and mediated the synthesis of estrogen (Yamaguchi et al., 2007), but it also bound to steroidogenic factor-1 (sf1/nr5a1) to enhance the transcriptional activity of cyp19a1a (Wang et al., 2007). However, no sexual specific expression of nr5a1a or nr5a1b was observed in Mandarin fish, while nr5a2 was highly expressed in the testis. Evidence proves the antagonistic effect of nr5a1 and nr5a2 (Shi et al., 2019), but the role of nr5a2 in this process remains to be explored.

Dmrt1 is a testis-biased gene that is highly conserved across animal phyla and was also highly expressed in male Mandarin fish. Previous studies have demonstrated that dmrt1 not only indirectly down-regulates the expression of cyp19a1a by inhibiting foxl2 (Li et al., 2013; Lindeman et al., 2015), but it also directly binds to the cyp19a1a promoter to down-regulate aromatase activity (Wang et al., 2010). Dmrt1 can promote the expression activity of sox9 (Wei et al., 2019), which has also an antagonistic effect on foxl2 (Nef and Vassalli, 2009; Georges et al., 2014). In Mandarin fish transcriptome data, only one sox9 gene with high testicular expression was found, unlike most teleost fish that have both sox9a and sox9b. In addition, amh, an important downstream gene related to sox9 and dmrt1 (Rodríguez-Marí et al., 2005; Webster et al., 2017), was shown to directly inhibit the expression of cyp19a1a (Rouiller-Fabre et al., 1998), which is also highly expressed in the testes of Mandarin fish.

The combination of existing studies and transcriptome data from sex related genes with differential expression, suggests that Mandarin fish might share a similar regulation pattern for sex-related genes with other teleosts (Fig. 6).

Fig.6 Regulation patterns of sex-related genes in Siniperca chuatsi High expression genes related to females and males are in the red and green boxes, respectively. The arrows show the promoting effect and the horizontal lines show the inhibiting effect.
5 CONCLUSION

In this study, high-throughput RNA sequencing was used to provide transcriptome data of Siniperca chuatsi. The molecular mechanisms of gene function, reproductive regulation and gametogenesis in the HPG axis of Mandarin fish were predicted. A large number of candidate genes for gender preference have been identified, and these should be further studied to better understand their exact functions in these processes. In addition, a large number of SSR and SNP markers were found, providing basic information for the marker-assisted breeding of Siniperca chuatsi.

6 DATA AVAILABILITY STATEMENT

All raw reads of transcriptome sequencing data have been deposited at the NCBI Short Read Archive (SRA) database (SRA accession Nos.: SRR11743000, SRR11743001, SRR11743002, SRR11743003).

Electronic supplementary material

Supplementary material (Supplementary Tables S1–S3 & Figs.S1–S2) is available in the online version of this article at https://doi.org/10.1007/s00343-020-0251-y.

References
Chen S L, Zhang G J, Shao C W, Huang Q F, Liu G, Zhang P, Song W T, An N, Chalopin D, Volff J N, Hong Y H, Li Q Y, Sha Z X, Zhou H L, Xie M S, Yu Q L, Liu Y, Xiang H, Wang N, Wu K, Yang C G, Zhou Q, Liao X L, Yang L F, Hu Q M, Zhang J L, Meng L, Jin L J, Tian Y S, Lian J M, Yang J F, Miao G D, Liu S S, Liang Z, Yan F, Li Y Z, Sun B, Zhang H, Zhang J, Zhu Y, Du M, Zhao Y W, Schartl M, Tang Q S, Wang J. 2014. Whole-genome sequence of a flatfish provides insights into ZW sex chromosome evolution and adaptation to a benthic lifestyle. Nature Genetics, 46(3): 253-260. DOI:10.1038/ng.2890
Chu W Y, Fu G H, Chen J, Chen D G, Meng T, Zhou R X, Xia X J, Zhang J S. 2010. Gene expression profiling in muscle tissues of the commercially important teleost, Siniperca chuatsi L. Aquaculture International, 18(4): 667-678. DOI:10.1007/s10499-009-9289-8
Du X X, Wang B, Liu X M, Liu X B, He Y, Zhang Q Q, Wang X B. 2017. Comparative transcriptome analysis of ovary and testis reveals potential sex-related genes and pathways in spotted knifejaw Oplegnathus punctatus. Gene, 637: 203-210. DOI:10.1016/j.gene.2017.09.055
Fan Z F, You F, Wang L J, Weng S D, Wu Z H, Hu J W, Zou Y X, Tan X G, Zhang P J. 2014. Gonadal transcriptome analysis of male and female olive flounder (Paralichthys olivaceus). BioMed Research International, 2014: 291067. DOI:10.1155/2014/291067
Froese R, Pauly D. 2017. FishBase. World Wide Web electronic publication. www. fishbase. org. Accessed on 2020-5-1.
Georges A, Auguste A, Bessière L, Vanet A, Todeschini A L, Veitia R A. 2014. FOXL2: a central transcription factor of the ovary. Journal of Molecular Endocrinology, 52(1): R17-R33. DOI:10.1530/JME-13-0159
Guiguen Y, Fostier A, Piferrer F, Chang C F. 2010. Ovarian aromatase and estrogens: a pivotal role for gonadal sex differentiation and sex change in fish. General and Comparative Endocrinology, 165(3): 352-366. DOI:10.1016/j.ygcen.2009.03.002
Hagihara S, Yamashita R, Yamamoto S, Ishihara M, Abe T, Ijiri S, Adachi S. 2014. Identification of genes involved in gonadal sex differentiation and the dimorphic expression pattern in undifferentiated gonads of Russian sturgeon Acipenser gueldenstaedtii Brandt & Ratzeburg, 1833. Journal of Applied Ichthyology, 30(6): 1557-1564. DOI:10.1111/jai.12588
Hattori R S, Murai Y, Oura M, Masuda S, Majhi S K, Sakamoto T, Fernandino J I, Somoza G M, Yokota M, Strüssmann C A. 2012. A Y-linked anti-Müllerian hormone duplication takes over a critical role in sex determination. Proceedings of National Academy of Sciences of the United States of America, 109(8): 2955-2959. DOI:10.1073/pnas.1018392109
He S, Liang X F, Sun J, Li L, Yu Y, Huang W, Qu C M, Cao L, Bai X L, Tao Y X. 2013. Insights into food preference in hybrid F1 of Siniperca chuatsi (♀) x Siniperca scherzeri (♂) mandarin fish through transcriptome analysis. BMC Genomics, 14: 601. DOI:10.1186/1471-2164.14-601
He Y H, Li L, Liang X F, He S, Zhao L, Zhang Y P. 2018. Inhibitory neurotransmitter serotonin and excitatory neurotransmitter dopamine both decrease food intake in Chinese perch (Siniperca chuatsi). Fish Physiology and Biochemistry, 44(1): 175-183. DOI:10.1007/s10695-017-0422-8
Hildahl J, Sandvik G K, Edvardsen R B, Fagernes C, Norberg B, Haug T M, Weltzien F A. 2011. Identification and gene expression analysis of three GnRH genes in female Atlantic cod during puberty provides insight into GnRH variant gene loss in fish. General and Comparative Endocrinology, 172(3): 458-467. DOI:10.1016/j.ygcen.2011.04.010
Ikemoto T, Park M K. 2005. Identification and molecular characterization of three GnRH ligands and five GnRH receptors in the spotted green pufferfish. Molecular and Cellular Endocrinology, 242(1-2): 67-79. DOI:10.1016/j.mce.2005.07.004
Janzen F J. 1995. Experimental evidence for the evolutionary significance of temperature-dependent sex determination. Evolution, 49(5): 864-873. DOI:10.1111/j.1558-5646.1995.tb02322.x
Kamiya T, Kai W, Tasumi S, Oka A, Matsunaga T, Mizuno N, Fujita M, Suetake H, Suzuki S, Hosoya S, Tohari S, Brenner S, Miyadai T, Venkatesh B, Suzuki Y, Kikuchi K. 2012. A trans-species missense SNP in Amhr2 is associated with sex determination in the tiger pufferfish, Takifugu rubripes (fugu). PLoS Genetics, 8(7): e1002798. DOI:10.1371/journal.pgen.1002798
Lethimonier C, Madigou T, Muñoz-Cueto J A, Lareyre J J, Kah O. 2004. Evolutionary aspects of GnRHs, GnRH neuronal systems and GnRH receptors in teleost fish. General and Comparative Endocrinology, 135(1): 1-16. DOI:10.1016/j.ygcen.2003.10.007
Li M H, Sun L, Wang D S. 2019a. Roles of estrogens in fish sexual plasticity and sex differentiation. General and Comparative Endocrinology, 277: 9-16. DOI:10.1016/j.ygcen.2018.11.015
Li M H, Sun Y L, Zhao J, Shi H J, Zeng S, Ye K, Jiang D N, Zhou L Y, Sun L N, Tao W J, Nagahama Y, Kocher T D, Wang D S. 2015. A tandem duplicate of anti-müllerian hormone with a missense SNP on the Y chromosome is essential for male sex determination in Nile Tilapia, Oreochromis niloticus. PLoS Genetics, 11(11): e1005678. DOI:10.1371/journal.pgen.1005678
Li M H, Yang H H, Li M R, Sun Y L, Jiang X L, Xie Q P, Wang T R, Shi H J, Sun L N, Zhou L Y, Wang D S. 2013. Antagonistic roles of Dmrt1 and Foxl2 in sex differentiation via estrogen production in tilapia as demonstrated by TALENs. Endocrinology, 154(12): 4814-4825. DOI:10.1210/en.2013-1451
Li Z J, Guo R T, Gu Z Z, Wang X N, Wang Y C, Xu H F, Wang C X, Liu X J. 2019b. Identification of a promoter element mediating kisspeptin-induced increases in GnRH gene expression in sheep. Gene, 699: 1-7. DOI:10.1016/j.gene.2019.03.006
Liang X F, Kiu J K, Huang B Y. 1998. The role of sense organs in the feeding behaviour of Chinese perch. Journal of Fish Biology, 52(5): 1058-1067. DOI:10.1111/j.1095-8649.1998.tb00603.x
Lindeman R E, Gearhart M D, Minkina A, Krentz A D, Bardwell V J, Zarkower D. 2015. Sexual cell-fate reprogramming in the ovary by DMRT1. Current Biology, 25(6): 764-771. DOI:10.1016/j.cub.2015.01.034
Lu S Q, Liu F, Liu Z, Zhang J S, Xie X M. 2008. Comparison in cloning and sequence of growth hormone gene in three species of genus Siniperca. Oceanologia et Limnologia Sinica, 39(4): 354-361. (in Chinese with English abstract)
Matsuda M, Nagahama Y, Shinomiya A, Sato T, Matsuda C, Kobayashi T, Morrey C E, Shibata N, Asakawa S, Shimizu N, Hori H, Hamaguchi S, Sakaizumi M. 2002. DMY is a Y-specific DM-domain gene required for male development in the medaka fish. Nature, 417(6888): 559-563. DOI:10.1038/nature751
Mechaly A S, Viñas J, Piferrer F. 2011. Gene structure analysis of kisspeptin-2 (Kiss2) in the Senegalese sole (Solea senegalensis): characterization of two splice variants of Kiss2, and novel evidence for metabolic regulation of kisspeptin signaling in non-mammalian species. Molecular and Cellular Endocrinology, 339(1-2): 14-24. DOI:10.1016/j.mce.2011.03.004
Miller W L. 2007. Steroidogenic acute regulatory protein (StAR), a novel mitochondrial cholesterol transporter. Biochimica et Biophysica Acta (BBA) — Molecular and Cell Biology of Lipids, 1771(6): 663-676. DOI:10.1016/j.bbalip.2007.02.012
Moncaut N, Somoza G, Power D M, Canario A V M. 2005. Five gonadotrophin-releasing hormone receptors in a teleost fish: isolation, tissue distribution and phylogenetic relationships. Journal of Molecular Endocrinology, 34(3): 767-779. DOI:10.1677/jme.L01757
Myosho T, Otake H, Masuyama H, Matsuda M, Kuroki Y, Fujiyama A, Naruse K, Hamaguchi S, Sakaizumi M. 2012. Tracing the emergence of a novel sex-determining gene in medaka, Oryzias luzonensis. Genetics, 191(1): 163-170. DOI:10.1534/genetics.111.137497
Nanda I, Kondo M, Hornung U, Asakawa S, Winkler C, Shimizu A, Shan Z H, Haaf T, Shimizu N, Shima A, Schmid M, Schartl M. 2002. A duplicated copy of DMRT1 in the sex-determining region of the Y chromosome of the medaka, Oryzias latipes. Proceedings of the National Academy of Sciences of the United States of America, 99(18): 11778-11783. DOI:10.1073/pnas.182314699
Nef S, Vassalli J D. 2009. Complementary pathways in mammalian female sex determination. Journal of Biology, 8(8): 74. DOI:10.1186/jbiol173
Novaira H J, Ng Y, Wolfe A, Radovick S. 2009. Kisspeptin increases GnRH mRNA expression and secretion in GnRH secreting neuronal cell lines. Molecular and Cellular Endocrinology, 311(1-2): 126-134. DOI:10.1016/j.mce.2009.06.011
Ogino Y, Tohyama S, Kohno S, Toyota K, Yamada G, Yatsu R, Kobayashi T, Tatarazako N, Sato T, Matsubara H, Lange A, Tyler C R, Katsu Y, Iguchi T, Miyagawa S. 2018. Functional distinctions associated with the diversity of sex steroid hormone receptors ESR and AR. The Journal of Steroid Biochemistry and Molecular Biology, 184: 38-46. DOI:10.1016/j.jsbmb.2018.06.002
Ospina-Álvarez N, Piferrer F. 2008. Temperature-dependent sex determination in fish revisited: prevalence, a single sex ratio response pattern, and possible effects of climate change. PLoS One, 3(7): e2837. DOI:10.1371/journal.pone.0002837
Reinhart A J, Williams S C, Stocco D M. 1999. Transcriptional regulation of the StAR gene. Molecular and Cellular Endocrinology, 151(1-2): 161-169. DOI:10.1016/S0303-7207(98)00257-3
Robinson-Rechavi M, Laudet V. 2001. Evolutionary rates of duplicate genes in fish and mammals. Molecular Biology and Evolution, 18(4): 681-683. DOI:10.1093/oxfordjournals.molbev.a003849
Robinson-Rechavi M, Marchand O, Escriva H, Bardet P L, Zelus D, Hughes S, Laudet V. 2001. Euteleost fish genomes are characterized by expansion of gene families. Genome Research, 11(5): 781-788. DOI:10.1101/gr.165601
Rodríguez-Marí A, Yan Y L, Bremiller R A, Wilson C, Cañestro C, Postlethwait J H. 2005. Characterization and expression pattern of zebrafish anti-Müllerian hormone (amh) relative to sox9a, sox9b, and cyp19a1a, during gonad development. Gene Expression Patterns, 5(5): 655-667. DOI:10.1016/j.modgep.2005.02.008
Rouiller-Fabre V, Carmona S, Merhi R A, Cate R, Habert R, Vigier B. 1998. Effect of anti-Mullerian hormone on sertoli and leydig cell functions in fetal and immature rats. Endocrinology, 139(3): 1213-1220. DOI:10.1210/endo.139.3.5785
Schartl M, Wilde B, Schlupp I, Parzefall J. 1995. Evolutionary origin of a parthenoform, the Amazon molly Poecilia formosa, on the basis of a molecular genealogy. Evolution, 49(5): 827-835. DOI:10.2307/2410406
Selstam G, Rosberg S, Liljekvist J, Grönquist L, Perklev T, Ahrén K. 1976. Differences in action of LH and FSH on the formation of cyclic amp in the prepubertal rat ovary. European Journal of Endocrinology, 81(1): 150-164. DOI:10.1530/acta.0.0810150
Selvaraj S, Kitano H, Fujinaga Y, Ohga H, Yoneda M, Yamaguchi A, Shimizu A, Matsuyama M. 2010. Molecular characterization, tissue distribution, and mRNA expression profiles of two Kiss genes in the adult male and female chub mackerel (Scomber japonicus) during different gonadal stages. General and Comparative Endocrinology, 169(1): 28-38. DOI:10.1016/).ygcen.2010.07.011
Shi B Y, Lu H J, Zhang L H, Zhang W M. 2019. Nr5a1b promotes and Nr5a2 inhibits transcription of lhb in the orange-spotted grouper, Epinephelus coioides. Biology of Reproduction, 101(4): 800-812. DOI:10.1093/biolre/ioz121
Simpson E R, Mahendroo M S, Means G D, Kilgore M W, Hinshelwood M M, Graham-Lorence S, Amarneh B, Ito Y, Fisher C R, Michael M D, Mendelson C R, Bulun S E. 1994. Aromatase cytochrome P450, the enzyme responsible for estrogen biosynthesis. Endocrine Reviews, 15(3): 342-355. DOI:10.1210/er.15.3.342
Sreenivasan R, Cai M, Bartfai R, Wang X, Christoffels A, Orban L. 2008. Transcriptomic analyses reveal novel genes with sexually dimorphic expression in the zebrafish gonad and brain. PLoS One, 3(3): e1791. DOI:10.1371/journal.pone.0001791
Stocco D M, Wang X J, Jo Y, Manna P R. 2005. Multiple signaling pathways regulating steroidogenesis and steroidogenic acute regulatory protein expression: more complicated than we thought. Molecular Endocrinology, 19(11): 2647-2659. DOI:10.1210/me.2004-0532
Sun F, Liu S K, Gao X Y, Jiang Y L, Perera D, Wang X L, Li C, Sun L Y, Zhang J R, Kaltenboeck L, Dunham R, Liu Z J. 2013. Male-biased genes in catfish as revealed by RNA, Seq analysis of the testis transcriptome. PLoS One, 8(7): e68452. DOI:10.1371/journal.pone.0068452
Tao W J, Yuan J, Zhou L Y, Sun L N, Sun Y L, Yang S J, Li M H, Zeng S, Huang B F, Wang D S. 2013. Characterization of gonadal transcriptomes from Nile tilapia (Oreochromis niloticus) reveals differentially expressed genes. PLoS One, 8(5): e63604. DOI:10.1371/journal.pone.0063604
Tu J G, Tian C X, Zhao P Q, Sun J X, Wang M, Fan Q X, Yuan Y C. 2017. Identification and profiling of growth-related microRNAs in Chinese perch (Siniperca chuatsi). BMC Genomics, 18(1): 489. DOI:10.1186/s12864.017-3851-y
Wang D S, Kobayashi T, Zhou L Y, Paul-Prasanth B, Fumie S, Sakai F, Okubo K, Morohashi K, Nagahama Y. 2007. Foxl2 up-regulates aromatase gene transcription in a female-specific manner by binding to the promoter as well as interacting with ad4 binding protein/steroidogenic factor 1. Molecular Endocrinology, 21(3): 712-725. DOI:10.1210/me.2006-0248
Wang D S, Zhou L Y, Kobayashi T, Matsuda M, Shibata Y, Sakai F, Nagahama Y. 2010. Doublesex, and Mab-3-related transcription factor-1 repression of aromatase transcription, a possible mechanism favoring the male pathway in tilapia. Endocrinology, 151(3): 1331-1340. DOI:10.1210/en.2009-0999
Wang G L, Wang M C, Zhang X W, Chang M X, Xie H X, Nie P. 2017a. Molecular cloning, biological effect, and tissue distribution of interleukin-8 protein in mandarin fish (Siniperca chuasti) upon Flavobacterium columnare infection. Fish & Shellfish Immunology, 66: 112-119. DOI:10.1016/j.fsi.2017.05.016
Wang G, Li J H, Zou P F, Xie H X, Huang B, Nie P, Chang M X. 2012. Expression pattern, promoter activity and bactericidal property of β-defensin from the mandarin fish Siniperca chuatsi. Fish & Shellfish Immunology, 33(3): 522-531. DOI:10.1016/j.fsi.2012.06.003
Wang K Z, Zhu X, Li Y L, Chen D X, Wu P, Chu W Y. 2016. Molecular characterization and expression regulation of Smyd1a and Smyd1b in skeletal muscle of Chinese perch (Siniperca chuatsi). Comparative Biochemistry and Physiology Part B: Biochemistry and Molecular Biology, 194-195: 25-31. DOI:10.1016/j.cbpb.2016.01.004
Wang P F, Zeng S, Xu P, Zhou L, Zeng L, Lu X, Wang H F, Li G F. 2014. Identification and expression analysis of two HSP70 isoforms in mandarin fish Siniperca chuatsi. Fisheries Science, 80(4): 803-817. DOI:10.1007/s12562-014-0747-5
Wang X Q, Li C W, XIE Z G, Fan W J, Zhang J S. 2006. Studies on the growth difference of the male and female Siniperca chuatsi. Freshwater Fisheries, 36(3): 34-37. (in Chinese with English abstract)
Wang Z C, Qiu X M, Kong D, Zhou X X, Guo Z B, Gao C F, Ma S, Hao W W, Jiang Z Q, Liu S C, Zhang T, Meng X S, Wang X L. 2017b. Comparative RNA-Seq analysis of differentially expressed genes in the testis and ovary of Takifugu rubripes. Comparative Biochemistry and Physiology Part D: Genomics and Proteomics, 22: 50-57. DOI:10.1016/j.cbd.2017.02.002
Warner R R. 1984. Mating behavior and hermaphroditism in coral reef fishes: the diverse forms of sexuality found among tropical marine fishes can be viewed as adaptations to their equally diverse mating systems. American Scientist, 72(2): 128-136.
Webster K A, Schach U, Ordaz A, Steinfeld J S, Draper B W, Siegfried K R. 2017. Dmrt1 is necessary for male sexual development in zebrafish. Developmental Biology, 422(1): 33-46. DOI:10.1016/j.ydbio.2016.12.008
Wei L, Li X Y, Li M H, Tang Y H, Wei J, Wang D S. 2019. Dmrt1 directly regulates the transcription of the testis-biased Sox9b gene in Nile tilapia (Oreochromis niloticus). Gene, 687: 109-115. DOI:10.1016/j.gene.2018.11.016
Weltzien F A, Andersson E, Andersen Ø, Shalchian-Tabrizi K, Norberg B. 2004. The brain-pituitary-gonad axis in male teleosts, with special emphasis on flatfish (Pleuronectiformes). Comparative Biochemistry and Physiology Part A: Molecular & Integrative Physiology, 137(3): 447-477. DOI:10.1016/j.cbpb.2003.11.007
Yamaguchi T, Yamaguchi S, Hirai T, Kitano T. 2007. Follicle-stimulating hormone signaling and Foxl2 are involved in transcriptional regulation of aromatase gene during gonadal sex differentiation in Japanese flounder, Paralichthys olivaceus. Biochemical and Biophysical Research Communications, 359(4): 935-940. DOI:10.1016/j.bbrc.2007.05.208
Yano A, Guyomard R, Nicol B, Jouanno E, Quillet E, Klopp C, Cabau C, Bouchez O, Fostier A, Guiguen Y. 2012. An immune-related gene evolved into the master sex, determining gene in rainbow trout, Oncorhynchus mykiss. Current Biology, 22(15): 1423-1428. DOI:10.1016/j.cub.2012.05.045
Zhang G Q, Chu W Y, Hu S N, Meng T, Pan L L, Zhou R X, Liu Z, Zhang J S. 2011. Identification and analysis of muscle-related protein isoforms expressed in the white muscle of the mandarin fish (Siniperca chuatsi). Marine Biotechnology, 13(2): 151-162. DOI:10.1007/s10126-010-9275-1
Zou L J, Gong J, Ji L, Zhuang Y H, Liao H Y, Huang W. 2017. Cloning and expression analysis of CYP19a gene in mandarin fish Siniperca chuatsi. Life Science Research, 21(4): 295-301. (in Chinese with English abstract)