Discovery of Sex-Determining Genes in Portunus trituberculatus: A Comparison of Male and Female Transcriptomes During Early Developmental Stages

Portunus trituberculatus is one of the main mariculture crabs of high economic value. To identify genes involved in sex determination, we first performed sex-specific transcriptome sequencing at six larval development stages using a DNA/RNA co-extraction method. A total of 907,952,938 and 828,774,880 reads were obtained from female and male crabs, respectively. 2,379 differentially expressed genes (DEGs) were found between females and males, and were mainly enriched in DNA replication, folate biosynthesis, and retinol metabolism pathways. Furthermore, transcription patterns of genes in the sex-determining region (SD) were analyzed based on the transcriptome data, and one Dmrt gene (PtDMY) was found to be exclusively expressed in males during early developmental stages. Notably, some known sex-related genes, including IAG, Dmrt11E, DmrtB1, and DmrtC2 were significantly down-regulated after knocking down PtDMY. Our results suggested that PtDMY is involved in sex determination and may be one of the key upstream regulators of the sex determination pathway. In addition, the massive volume of transcriptome data obtained in this study provided an important basis for the systematic study of sex determination mechanisms in P. trituberculatus.


INTRODUCTION
Crustaceans are the second largest class of arthropods (Alwi et al., 2015). There are multiple sex determination systems in crustaceans that have been identified by karyotype analyses, including XY/ZW systems (Cui et al., 2015;Fang et al., 2020). It is an important material to study the evolution of sex staining and the mechanisms of sex determination (Torrecilla et al., 2017;Toyota et al., 2021). In recent years, with the development of sequencing technologies, sex-related quantitative trait locus (QTL) have been mapped and sex markers have been identified in many crustaceans of economic significance, including Litopenaeus vannamei, Eriocheir sinensis, Scylla paramamosain, and Exopalaemon carinicauda (Cui et al., 2015;Yu et al., 2015;Shi et al., 2018;Fang et al., 2020;Lv et al., 2020a). Multiple studies have contributed to the identification of sex determination genes and have helped clarify the mechanisms of sex determination in crustaceans.
Sex-determining genes, such as sex-determining region of Y-chromosome (Sry) and double-sex and mab-3 related transcription factor (Dmrt), have been identified in mammals, fish, birds, and reptiles (Masaru et al., 2002;Prokop et al., 2020;Han et al., 2021;Naoki and Makoto, 2021). Most sexdetermining genes are transcription factors located on the sex chromosomes. For example, Dmrt in Oryzias latipes is located on the Y chromosome (Feng et al., 2021), and Dmrt in Xenopus laevis is located on the W chromosome (Yoshimoto et al., 2008).
Studies have shown that the androgenic gland (AG) is an important organ for sex determination, and the Insulin-like androgenic gland hormone (IAG) secreted by the AG is necessary for male genesis (Sun and Li, 2021). Notably, the knockdown of IAG can lead to sex reversal in Macrobrachium rosenbergii (Lezer et al., 2015). Recently, a Y-linked Dmrt was also identified in Sagmariasus verreauxi (Chandler et al., 2017), which shed new light on the prevalence of Dmrt genes as master sex-determinants in crustaceans. However, studies on sex determination genes in crustaceans are relatively backward, and the regulatory networks of sex determination are not clear, thus, limiting the identification of the mechanisms of sex determination in crustaceans.
P. trituberculatus, known as the swimming crab, is the main mariculture crab in China with important economic value and an annual production of ∼100,000 tons. The growth traits of P. trituberculatus are significantly different between males and females, and therefore, it is of great economic significance to detail the mechanisms of sex determination and establish sexcontrol technologies. We located the sex QTL and identified sex markers in P. trituberculatus (Lv et al., 2018), and established sex identification technologies based those markers . We also confirmed the XY sex determination system in P. trituberculatus according to the differences in sex markers between males and females. However, due to the lack of genomic information, unambiguous sex-determining genes have not been identified.
In this study, we first performed sex-specific transcriptome sequencing of six larval development stages using a DNA/RNA co-extraction method. Through comparative transcriptome analyses, DEGs between males and females were identified and enriched. In addition, we analyzed the transcription patterns of genes in the SD during the early development stages and identified a key candidate gene for sex determination. Our results contribute to the interpretation of sex determination mechanisms in P. trituberculatus.

Experimental Materials
P. trituberculatus used in this study was cultured in Changyi Haifeng Aquaculture Co., Ltd., Weifang City, Shandong Province, China. There are 12 developmental stages from embryonic development to adult crab (Figure 1), including the fertilized egg (multicellular stage, Mc; blastula stage, B; gastrulation stage, G; egg-nauplius stage, En; Egg-zoea stage, Ez), zoea (zoea 1, Z1; zoea 2, Z2; zoea 3, Z3; zoea 4, Z4), megalopa (M), and juvenile crab (juvenile crab 1, J1; juvenile crab 2, J2). Samples at each stage were separately picked and placed into RNase-free centrifuge tubes (more than nine individuals per tube and three biological replicates for each period) and immediately placed into liquid nitrogen. Samples were then used to prepare cDNA for qRT-PCR. An additional sample was collected during Z2-J2 (more than 30 individuals per stage), and was immediately soaked in absolute ethanol and placed in liquid nitrogen. Samples were then used for DNA/RNA co-extraction of each individual.

DNA/RNA Co-extraction and Sex Identification
Individual DNA and RNA were extracted simultaneously using the SteadyPure Universal RNA Extraction kit (Accurate Biology, Hunan). The DNA has been tested by agarose gel electrophoresis to determine the fragment size. And the OD of DNA was measured by uv spectrophotometer, the OD260/OD280 was about 1.8. DNA was used for sex identification using previously established sex identification methods . RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, United States).

Library Construction and Transcriptome Sequencing
After sex identification, the RNA of three or four individuals of the same sex were mixed into one sample, and three parallel biological samples were prepared for males and females at each developmental stage. To select cDNA fragments that were 370-420 bp in length, the library fragments were purified using the AMPure XP system (Beckman Coulter, Beverly, United States). Clustering of the index-coded samples was performed on a cBot Cluster Generation System using the TruSeq PE Cluster Kit, v3-cBot-HS (Illumina), according to the manufacturer's instructions. Following cluster generation, the libraries were sequenced on an Illumina Novaseq platform with 150 bp paired-end read datasets generated.

Quality Control
Clean data were obtained by removing low quality reads, and reads containing adapters and ploy-N sequences. The Q20, Q30, and GC contents of the clean data were calculated. The index of the reference genome was built using Hisat2 v2.0.5, and paired-end clean reads were aligned to the reference genome using Hisat2 v2.0.5.

Quantification of Gene Transcription
The mapped reads of each sample were assembled using StringTie (v1.3.3b) (Pertea et al., 2016) in a reference-based approach. FeatureCounts v1.5.0-p3 was used to count the number of reads mapped to each gene. The fragments per kilobase of transcript sequence per million base pairs sequenced (FPKM) of each gene was then calculated based on the length of the gene and read counts mapped to the gene. The FPKM considers the effect  of the sequencing depth and gene length for the read counts at the same time.

Differential Transcription Analysis
Differential gene transcription analyses of males and females were performed using the DESeq2 R package (1.20.0) (Love et al., 2014). DESeq2 provides statistical routines to determine the differential transcription of digital gene transcription data using a model based on a negative binomial distribution. Genes with an adjusted p-value ≤ 0.05 and a | log2FoldChange| ≥ 1.0 were considered differentially expressed.

Quantitative Real-Time PCR
Primers were designed using Primer 5.0 software . The volume of the qRT-PCR was 10 µl, containing 5 µL of 2 × SYBR Master Mix (RR420A, Takara Bio, Japan), 0.2 µL FIGURE 2 | Differential gene transcription analyses in the early stages of growth.   of each primer and ROX Reference Dye II, 1 µL of cDNA template, and 3.4 µl of RNase free dH 2 O. qRT-PCR was carried out in a FAST-7500 system (ABI-7500, Thermo Fisher Scientific, Singapore) with the following parameters: 95 • C for 30 s; 40 cycles of 95 • C for 5 s, 60 • C for 30 s, and 72 • C for 30 s. The relative quantification of qRT-PCR data was calculated using the 2 − Ct method (Livak and Schmittgen, 2001).

RNA Interference
Z1 stage is the first period when eggs hatch into larvae, and the expression of PtDMY was relatively high in the Z1 stage. Therefore, Z1 stage samples were used for RNAi. The PtDMY gene sequence was sent to Shanghai Sangon Biology Co., Ltd. for primer design and preparation of RNAi reagents (G3-Sense, GCGCAUCCAUACUAACAAATT; G3-Antisense, UUUGUUAGUAUGGAUGCGCTT; NC-Sense, UUCUCCGAACGUGUCACGUTT; NC-Antisense, ACGUGACACGUUCGGAGAATT). RNAi was performed by feeding Rotifers soaked in RNAi reagent to crabs. The concentration of RNAi reagent for soaking was 3.3 µg/mL, and feeding was initiated after soaking for 2-3 h. The Rotifers were fed 3 times per day. Crabs at Z1 stage were sampled after feeding for 48 h. RNA extraction and reverse transcription were performed on the collected samples, and the effect of RNAi was detected using qRT-PCR.

Transcriptome Sequencing
Male and female samples from six developmental stages were sequenced on an Illumina platform and 86.84 Gb of clean bases were produced. In total, 907,952,938 and 828,774,880 reads were obtained from females and males, respectively. The mean values of Q20 and Q30 were 96.355 and 91.6175, respectively. The average percentage of genome mapping was 84.21% ( Table 1). The raw data were deposited in the Short Read Archive of the National Center for Biotechnology Information (NCBI) with accession number PRJNA631920.

Differentially Expressed Genes Analysis
Through comparative analyses, a total of 2,379 significant DEGs were found between females and males (p-value ≤ 0.05 and | log2FoldChange| ≥ 1.0). The number of DEGs were highest in the Z4 and J2 periods, accounting for 20.22 and 33.12%, respectively (Figure 2). Gene Ontology (GO) enrichment analyses showed that chitin metabolic and glucosaminecontaining compounds metabolic, among others, were the most significantly enriched (Figure 3). Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses revealed that 20 pathways were enriched. The pathways for DNA replication, retinol metabolism, and neuroactive ligand-receptor interactions were the most significantly enriched (Figure 4).

Transcription Patterns of Genes in the Sex-Determining Region
In our previous study, we located the SD region on two contigs in P. trituberculatus (Lv et al., 2020b). Based on the genomic annotations, 41 genes were found in the SD  Table 2), among which, eight genes were significantly differentially expressed between females and males. Further analyses showed that the transcription of three genes in males was significantly higher than that in females at all stages of development (Figure 5). In particular, Contig 475.13 was expressed almost exclusively in males, and belonged to the Dmrt family (named PtDMY).

Expression Analysis of PtDMY
To further clarify the expression pattern of PtDMY in early development, qRT-PCR was performed using samples from 12 developmental periods from fertilized egg to juvenile crab, including Z1-J2 which were sex-specific. The results showed that PtDMY was almost undetectable in early embryonic development (Mc-G) and was up-regulated from the En phase to a peak at the Ez-Z1, before being down-regulated. It is worth noting that the expression of PtDMY in Z1-J2 was only detected in males, with almost no expression in females. Thus, those results were consistent with the transcriptome results (Figure 6).

RNAi for PtDMY
In this study, we successfully knocked down the expression of PtDMY in zoea larva after 48 h of RNAi (Figure 7). Furthermore, we analyzed the expression of known sex-related genes after PtDMY knockdown, including IAG, Dmrt, and SPATA. The results showed that most genes were significantly down-regulated, among which, IAG decreased by 78.48% and Dmrt decreased by an average of 66.49% (Figure 8).

DISCUSSION
Sex determination and differentiation are among the most interesting topics in biology, however, such processes remain unclear for most crustaceans . P. trituberculatus is one of the main mariculture crabs worldwide, with high economic value and differences in values between males and females. Thus, the identification of sex determination mechanisms in P. trituberculatus will facilitate sex control technologies for monosex breeding during production.
Sex determination in crustaceans generally occurs during early developmental stages (Aflalo et al., 2006). Therefore, sexspecific transcriptome sequencing at early developmental stage will provide the necessary information for the identification of genes related to sex determination. In fact, the sex determining genes in the silkworm moth, Bombyx mori, have been successfully identified using this approach (Takashi et al., 2014). Beside, DEGs between female and male L. vannamei at the Z1 stage have also been obtained using transcriptome sequencing . However, there have been no reports on sex-specific transcriptome sequencing at that early developmental stage in crabs, and thus, research on genetic sex markers is currently lacking. For the first time, this study sequenced and analyzed the sex-specific transcriptomes of crabs at early developmental stages. Compared to previous work in shrimp , our study included six developmental stages and obtained more sequencing data (272 G), which allowed us to identify sex-related genes more robustly.
A total of 2,379 DEGs were identified between female and male transcriptomes from Z2 to J2, and suggested that many genes may be involved in sex determination or differentiation. GO enrichment results showed that genes involved in chitin metabolic processes (GO:0006030), glucosamine-containing compound metabolic processes (GO:1901071), and amino sugar metabolic processes (GO:0006040) were the most significantly enriched in biological process (BP), which were related to molting, development, and growth (Kimura et al., 2008;Rocha et al., 2012). Crustacean growth requires molting, and the growth rate of male and female crabs differs between young crabs to adult crabs . DEGs were enriched in processes related to molting and development, which indicated that the differentiation of growth and development of male and female crabs may have begun in the early stage of larval development.
The Dmrt, which is characterized by one or more highly conserved DNA-binding (DM) domains, was involved in sex determination (Meng et al., 1999). In Medaka (XX/XY sexdetermination system), the Y-linked DMY gene was characterized as sex-determining genes, leading to male sexual development (Masaru et al., 2002). In X. laevis (ZZ/ZW sex-determination system), one W-linked DM-W gene was identified as a likely sex-determining gene (Yoshimoto et al., 2008). In Chicken and Flatfish (ZZ/ZW sex-determination system), Dmrt was located on the Z chromosome, and the level of expression may induce male development (Nanda et al., 2000;Chen et al., 2014). Recently, an invertebrate sex-linked Dmrt was found in Sagmariasus verreauxi (Chandler et al., 2017), suggesting that it also plays an important role in sex determination in crustaceans. In this work, one Dmrt gene in the SD showed obvious male-specific transcription, which was consistent with the characteristics of closely linked with males (Lv et al., 2020b). We named that gene PtDMY and hypothesized that it played a key role in sex determination.
To verify the function of PtDMY, we successfully knocked down its expression during larval development using RNAi technology. We found multiple known sex-related genes were significantly down-regulated after PtDMY knockdown, including IAG and Dmrt (Figure 8). Previous studies showed that IAG secreted by the AG is necessary for male genesis (Sun and Li, 2021), and is regulated by Dmrt in P. trituberculatus . It should be noted that three Dmrt genes were significantly down regulated after PtDMY knockdown, which strongly suggested that PtDMY may be a key upstream factor in the sex determination pathway. Furthermore, we systematically analyzed the expression pattern of PtDMY from the multicellular stage of crab to the juvenile stage. The results showed that its expression peaked in the Ez to Z1 stages (Figure 6), indicating that this period may be key for sex determination.
Interestingly, no genes were continuously highly expressed during female development. There are several possibilities for this observation: (1) The low expression leads to the female related genes have not been found; (2) the effect of genes related to female genesis is short, and there is no need for long-term high expression; or (3) there are no genes in P. trituberculatus that determine female sex. Combining the results of this and previous studies (Takashi et al., 2014;Wang et al., 2020), a model for sex determination and differentiation in P. trituberculatus was preliminarily outlined (Figure 9). The species belongs to the male heterogametic XY system (Lv et al., 2019), where PtDMY is the key gene for sex determination, and is expressed almost exclusively in males during larval development (Z1-J2). Dmrt and IAG are important sex differentiation-related genes, and are regulated by PtDMY. Notably, multiple aspects of the sex determination mechanism in P. trituberculatus remain unknown, and our research is just begin. In follow-up studies, further research is required to more fully detail the sex-determination pathway in crab.

CONCLUSION
In this work, sex-specific transcriptome sequencing at six larval development stages was performed and 2,379 DEGs were found between females and males. one Dmrt gene located in SD was exclusively expressed in males during early developmental stages, which was considered to be a key upstream regulators of the sex determination pathway. Our results here provided an important basis for clarifying the sex determination mechanisms of P. trituberculatus.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI (accession: PRJNA631920).

ETHICS STATEMENT
The animal study was reviewed and approved by the Yellow Sea Fisheries Research Institute. Written informed consent was obtained from the owners for the participation of their animals in this study.

AUTHOR CONTRIBUTIONS
WZ conducted experiments, analyzed data and results, and manuscript writing. JL was experiment design, analyzed data and results, and manuscript revising. WL and BG conducted experiments. PL was experiment design and review manuscript. All authors contributed to the article and approved the submitted version.