Identification of Sex-Specific Markers Through 2b-RAD Sequencing in the Sea Urchin (Mesocentrotus nudus)

Sex-specific markers play an important role in revealing sex-determination mechanism. Sea urchin (Mesocentrotus nudus) is an economically important mariculture species in several Asian countries and its gonads are the sole edible parts for people. However, growth rate and immunocompetence differ by sex in this species, sex-specific markers have not been identified, and the sex-determination mechanism of sea urchin remains undetermined. In this study, type IIB endonuclease restriction-site associated DNA sequencing (2b-RAD-seq) and a genome survey of M. nudus were performed, and three female-specific markers and three female heterogametic single nucleotide polymorphism (SNP) loci were identified. We validated these sex-specific markers via PCR amplification in a large number of individuals, including wild and artificially bred populations. Several open reading frames (ORFs) were predicted, although there are no potential genes known for sex determination and sex differentiation within the scaffold in which the sex-specific markers are located. Importantly, the female-specific sequences and female heterozygous SNP loci indicate that a female heterogametic and male homogametic ZW/ZZ sex-determination system should exist in M. nudus. The results provide a solid basis for revealing the sex-determination mechanism of this species, and open up new possibilities for developing sex-control breeding in sea urchin.

Identification of Sex-Specific Markers Through 2b-RAD Sequencing in the Sea Urchin (Mesocentrotus nudus) INTRODUCTION Sexual dimorphism is common in aquaculture animals, and growth rate, immunocompetence, and body size generally differ significantly by sex (Mei and Gui, 2015;Jiang et al., 2017). Sea urchin is widely distributed and commercially harvested all over the world (Inomata et al., 2016), and the gonads are the only edible parts for people. In addition, edible sea urchins are important mariculture species in several Asian countries. This group of animals shows gender-specific differences in immune response, for example, the immunocompetence of female Paracentrotus lividus urchins might be superior to that of males (Arizza et al., 2013). Moreover, the gonadal growth rate in Abatus cavernosus urchins is faster in males than in females (Gil et al., 2020), and the types and the contents of free amino acids, as well as lipid concentrations, are also closely linked to sexual dimorphism in urchins (Osako et al., 2007;Martínez-Pita et al., 2010;Díaz de Vivar et al., 2019). However, most sea urchins lack conspicuous sexually dimorphic phenotypes, making it impossible to distinguish the phenotypic sex of living animals. This dramatically increases the costs and reduces the efficiency of genetic breeding. In addition to commercial value, sea urchins are excellent organisms for the study of early embryonic development (Juliano et al., 2010), interactions among genes (Yuh et al., 1998;Oliveri et al., 2002), evolution (Kober and Bernardi, 2013), fertilization (Chun et al., 2014), and immunity (Branco et al., 2014). Therefore, it is important to reveal the genetic mechanism of sex determination in this species. To do so, it is critical to develop reliable genetic sex-identification markers (Kobayashi et al., 2013;Liu et al., 2013;Mei and Gui, 2015;Chen et al., 2018). A successful system would improve breeding efficiency and make it possible to implement sex-control breeding in these echinoderms (Wang et al., 2009;Martínez et al., 2014).
Several methods have been applied to identify sex-specific molecular markers, including traditional methods and highthroughput sequencing methods. Amplified fragment length polymorphism (AFLP) is a traditional method to score random SNP markers across an entire genome. However, this approach is expensive and requires complex experiments (Griffiths and Orr, 1999;Flachowsky et al., 2001). With the development of nextgeneration sequencing, several methods have been developed to obtain thousands of SNPs in a single experiment, such as double-digest restriction-site-associated DNA sequencing (dd-RAD-seq), sequence-based genotyping (SBG) and type IIB endonucleases restriction-site-associated DNA sequencing (2b-RAD-seq). Dd-RAD-seq technology can generate number of SNPs allowing for the construction of phylogenetic trees, detection of genetic variation of a population and accurate genetic maps (Zhou et al., 2014;Puchta et al., 2018). Due to the multiplexing a large number of samples, dd-RAD-seq requires the appropriate tools for analysis (Peterson et al., 2012). Using SBG can achieve genome-wide SNP discovery and genotyping of large populations without the information of a reference genome sequence, and this approach had been successfully used in the plants (Truong et al., 2012). Compared to the other traditional RAD-seq methods, 2b-RAD-seq is cost-effective and simple method, and allows for nearly every restriction site in the genome to be screened in parallel (Wang et al., 2012). Some studies have already successfully applied the 2b-RAD approach to explore sexspecific DNA fragments and SNPs in a great diversity of aquatic organisms (Liu et al., 2018;Zhou et al., 2019;Zhu et al., 2021). In addition to quickly identify genetic sex, sexspecific markers also provide insights into the sex-determination system in various animals (Liu et al., 2018). In several species of Scylla crab, sex-specific SNPs have been identified by 2b-RAD-seq, among which female are heterogamety and male are homogamety, while indicating that a ZW/ZZ sexdetermination system may exist in the above organisms (Shi et al., 2018). In Mastacembelus armatus, two male-specific 2b-RAD tags and two SNPs were identified through 2b-RAD-seq, and reported that the Tcl gene was a candidate sex-determination gene (Xue et al., 2020). However, few studies have identified sex-specific markers in sea urchins.
Mesocentrotus nudus belongs to the Strongylocentrotidae family of Echinodermata, and is mainly distributed and farmed in coastal areas of several Asian countries, such as China, South Korea, and Japan (Agatsuma, 2020). To further understand the sex-determination mechanism in this species, we used 2b-RADseq to identify sex-specific markers. Subsequently, we verified the candidate markers in a large number of samples, including laboratory breeding population, cultured population, and wild populations, and perfect accordance was observed between sex-specific makers and sexual phenotype. Our results can be used to significantly improve breeding efficiency, and provide useful information for understanding the sex-determination mechanisms in sea urchins. To determine the gender of each sea urchin, the corresponding gonad tissues were sampled in 4% paraformaldehyde (PFA) for histological examination.

DNA Extraction
Genomic DNA was extracted following the protocols of a marine animal tissue genomic DNA extraction kit (Tiangen, DP324). The quality and concentration of DNA was determined by 1% agarose gel electrophoresis and using a SimpliNano spectrophotometer (Biochrom, United Kingdom).

Gonadal Histology and Identifying of Gender
Samples of gonad tissues were fixed in PFA at 4 • C, then washed three times with phosphate buffered solution (PBS) at room temperature. After incubating the samples in 30% sucrose for 3 h, they were embedded in optimal cutting temperature (OCT) compound and sectioned at 3 µm. Then they were stained with hematoxylin/eosin. The gender of each sea urchin was determined by observing the gonadal structure with a microscope (Leica DM4B microscope).

Library Construction and 2b-RAD Sequencing
The genomic DNA of 20 individuals, including 10 males and 10 females, was used for 2b-RAD-seq. The 2b-RAD libraries were constructed at Qingdao OE Biotech Co., Ltd. (Qingdao, China)  following previous methods with slight modification (Wang et al., 2012). Briefly, 200 ng genomic DNA from each individual was digested with 1 U BsaXI restriction enzyme, then 10 µL digestion product was ligated with adaptor. Next, adaptor primers were used to amplify the ligation products by PCR. The quality and concentration of each PCR product was determined on an 8% polyacrylamide gel. Purified PCR products were digested with 2U SapI restriction enzyme at 37 • C for 30 min, and the digested products were further purified using streptavidincoated magnetic beads. The supernatant was incubated with T4 DNA ligase at 16 • C for 45 min, followed by purification of the ligation products. Barcodes were introduced by PCR with barcode-bearing primers, the PCR products were purified before sequencing. Finally, each library was pooled for sequencing using the Illumina Nova PE150PE platform.

Data Filtering, Survey Genome, and Screening of Candidate Sex-Specific Markers
Paired-end reads were merged using Pear software (version 0.9.6) (Zhang et al., 2014), and terminal tag sequences were excluded from each output read. Reads with ambiguous bases (N) exceeding 8%, poor-quality reads, and those without restriction sites were excluded from the raw data. Ultimately, clean reads    (Catchen et al., 2011). For genome survey, K-mer distribution was estimated by using jellyfish1 (version: 2.2.6) 1 with parameters "-m 17 -C." The genome size and heterozygosity ratio were estimated by the GenomeScope2. To identify candidate sex-specific markers, high-quality 33 bp reads were aligned against the genome survey using SOAP2 software (Li et al., 2009), and sequences with fewer than three reads were precluded (Li et al., 2009). Finally, candidate sex-specific markers for each individual were generated for further analyses.

Validation of Sex Specific Sequences
According to the flank sequence of sex-specific markers, specific primers were designed using online software 2 according to the flank sequences of the sex-specific marker. Then PCR reactions were carried out in a total volume of 20 µL, containing 10 µL 2 X Taq buffer (KangWei), 1 µL forward and reverse primer (10 mM), and 50-100 ng template DNA, adding sterile water to reach the final volume. The PCR amplification program was as follows: predenaturation at 94 • C for 2 min, then 30 cycles of 94 • C for 30 s, 58 • C for 30 s, and a final extension of 2 min at 72 • C. The PCR products were determined via 1% agarose gel electrophoresis and sequenced by Tsingke Biotechnology (Tianjin, China).

2b-RAD-Seq
In total, the 2b-RAD-seq produced 137,684,505 raw reads, 76,661,950 for females and 61,022,555 for males. After filtering low-quality reads, 72,125,643 clean reads for females and 56,991,085 clean reads for males were obtained. Based on the read quality and quantity, two female and two male samples were selected as references, clean reads were clustered and assembled. Subsequently, the reads of 20 individuals were aligned with reference samples, with an average alignment rate of 72.04%. After alignment analysis, 2,384,989 tags with an average depth of 36.78 × were obtained from the 20 individuals, and 43,922 polymorphic SNP markers were also screened ( Table 1).
FIGURE 1 | Validation of three candidate female-specific 2b-RAD-tags in 10 females and 9 males. Label M represents the DL2,000 DNA marker.

Genome Survey
Because the 2b-RAD tags were only 27 bp in length, too short to design PCR primers, and the genomic information for M. nudus is unavailable, a male and a female individual were subjected to short-read de novo sequencing. This generated 648,230,824 raw reads from the female; after filtering low-quality reads, 620,114,248 clean reads were acquired with a GC content of 37.06%. Additionally, a total of 669,375,890 raw reads were obtained from male individual, and 639,002,820 clean reads with a GC content of 37.23% were produced after filtering.  The estimated genome size using SOAP de novo K-mer module (K = 17) was 626,279,328 bp with high heterozygosity (1.3%). The SRA raw reads were deposited in the GenBank public database with the accession number of PRJNA741812. Using these clean reads, a draft genome of each sex was assembled (1,582,413 scaffolds in female and 5,953,609 scaffolds in male). The longest scaffold was 52,656 bp with N50 lengths of 1,579 bp in the female genome, and 10,758 bp with N50 lengths of 150 bp in the male genome. Moreover, there were 78,550 (4.96%) and 943 (0.01%) scaffolds with lengths of more than 2 kb in the female and male draft genome, respectively ( Table 2).

Identification of Candidate Sex-Specific Tags and SNPs
To identify putative sex-specific makers, the 2b-RAD tags and SNPs were aligned against the genome survey. After data filtering, we obtained 13 putative female-specific tags, which were present  in 10 female samples (Table 3). In addition, 20,201 putative SNP loci were also obtained for further analysis. A smaller p-value of SNP means a higher correlation with sex; hence, we screened the top 10 sex-specific SNP locus with the lowest p-values for further analysis. Notably, the females appeared heterozygous whereas the males appeared homozygous for all 10-candidate sex-specific SNPs (Table 4), these data suggest that a ZW/ZZ sex-determination system exists in M. nudus.

Validation of Candidate Sex-Specific Tags and SNPs
To confirm the authenticity of the candidate sex-specific tags and SNPs, primers were designed according to the flanking sequences of the 2b-RAD tags and SNPs from the genome survey (Table 5).
Then, PCR reactions were performed using specific primers for 10 female and 9 male laboratory breeding sea urchins. As shown in Figure 1, the PCR products from 3 of the 13 female tags (female tags 3, 7, and 8) amplified one expected sizes band in females but not in males. Next, we purified these PCR products and sequenced them via Sanger sequencing. The nucleotide sequences were consistent with the female genomic sequences, further confirming the female-specific segments. Considering that the other 10 tags amplified products from both sexes (data not shown), these 10 tags were excluded from further study. In addition, the 10 candidate SNP loci were validated via PCR assays using specific primers from 20 individuals, and all PCR products were sequenced. Three of the candidate sex-specific SNPs were confirmed (Figure 2), and the others were excluded.

Genetic Sex Identification in M. nudus Using Sex-Specific Tags
To further validate the accuracy of the three female-specific tags, we performed PCR reactions on 95 sea urchins from breeding population outside the laboratory and wild population. Initially, the gender of each sea urchin was determined by gonadal histology, and 38 females (Figure 3) and 57 males (Figure 4) FIGURE 7 | Electrophoretogram of the amplified products using female-tag 7 for genotypic sex identification in M. nudus. Label M represents the DL2,000 DNA marker (500 and 250 bp DNA ladder are shown). were identified in breeding population, 10 females and 9 males were identified in wild population (Figure 5). Then the three female-specific tags, including female-tag 3, female-tag 7, femaletag 8, were used to identify the genetic sex of these individuals. As expected, all three female-specific tags were only detected in the females, but not in males (Figures 6-8). These data suggest that the three female-specific markers could be used to identify the genetic sex of M. nudus with a very high rate of accuracy.

Sex-Specific Sequence Annotation
The three female specific tags (female-tag 3, female-tag 7, femaletag 8) locate at scaffold 2838205, 655827, 439353, respectively. The scaffold 2838205 is 717 bp in length, and 655827, 439353 are 4,146 and 5,496 bp, respectively. To better understand these scaffolds, we predicated homologous genes by BLAST program. The databases of Softberry and NCBI, genomes of Strongylocentrotus purpuratus were used for annotation of the scaffolds. Although several open reading frames (ORFs) were predicted, there are no potential genes known for sex determination or differentiation within any scaffold (Figure 9). In addition, potential ORFs fail to be identified in the scaffold 2838205 due to a short sequence.

DISCUSSION
Mesocentrotus nudus is an economically important species in several Asian countries, but its benefits differ by sex and the gender of animals cannot be easily ascertained. In addition, sex-specific markers have not been identified, and the sexdetermination mechanism of sea urchin remains undetermined. According to previous reports, the chromosomes of sea urchins are too small to be analyzed accurately, and heteromorphic sex chromosomes cannot be observed in most species (Saotome, 1987;Saotome et al., 2002;Eno et al., 2009). As an exception, a heteromorphic chromosome has been observed in Paracentrotus lividus sea urchins, for which an XX/XY sex-determination system has also been deduced (Lipani et al., 1996). In many animals, sex determination system has been deduced using sexspecific markers which developed by 2b-RAD-seq. In redtail catfish (Mystus wyckioides), some Y-specific fragments and X-specific fragments have been identified by 2b-RAD-seq, and XX/XY sex determination system was deduced because of the X-specific fragment showing dosage effect between females and males (Zhou et al., 2019). In this study, we revealed that all of the identified sex-specific markers were female-specific, and the identified SNP markers were heterozygous in females but homozygous in males (Figures 1, 2), our data suggest that a ZW/ZZ sex-determination system exists in M. nudus. However, it has been reported that the diploid chromosome number of M. nudus is 42 (2n = 42), and no heteromorphic sex chromosomes had not been observed (Saotome, 1987).
In the future, we will use these female-specific sequences for chromosome fluorescence in situ hybridization (FISH) and attempt to locate sex chromosomes. Our results may help elucidate the origin of sex-determination systems.
Although it is difficult to observe heteromorphic sex chromosomes, sex-determining regions have been reported in several aquaculture animals (Wen et al., 2020;Xue et al., 2020). Liu et al. (2018) obtained a male-specific sequence 8,661 bp in length from both bighead carp (Hypophthalmichthys nobilis) and silver carp (Hypophthalmichthys molitrix) via 2b-RAD-seq, and suggested that it was located on the Y chromosome or in a nascent Y-linked region. Similarly, using male-specific 2b-RAD tags, Kirubakaran et al. (2019) identified a potential sexdetermining region of 9 kb in Atlantic cod, and suggested that the zkY gene located in this region was a candidate sex-determining gene. In the present study, we identified three female-specific regions, providing a solid basis to reveal the sex-determination mechanisms in sea urchin. After homologous analysis by BLAST, several ORFs were predicted in these female-specific scaffolds, but no potential genes known for sex determination or differentiation within any scaffold. This phenome is similar to most previously reports that even though the sex-specific markers were extended into long segments within hundreds of genes, no sex related genes have been obtained (Pan et al., 2015;Shi et al., 2018;Zhou et al., 2019;Fang et al., 2020). In the future, we will use genome walking to analyze the differences in sequences between male and female individuals. The complete genome should also be sequenced to further elucidate the linkages among the three sex-specific sequences, as well as to reveal the sex-determination mechanism in sea urchin.
The sex-specific markers identified in our study were used to identify gender of more than one hundred of sea urchins with 100% accuracy (Figures 6-8). However, the need for DNA sample preparation limits the efficiency of this method. A direct PCR method would significantly reduce the identification time (Wang et al., 1993;Satya et al., 2013;Gao et al., 2020). This method uses NaOH-incubated products rather than purified DNA as template, and DNA can be extracted in 30 min from the tube feet of sea urchins. Furthermore, we compared the effectiveness of PCR using two different DNA extraction methods (using a DNA extraction kit and the alkaline extraction DNA procedure) and perfect consistency was observed between the two (data not shown). However, most aquaculture farms do not have PCR equipment and other required devices, which limits the application of this method. The loop-mediated isothermal amplification (LAMP) method can quickly amplify DNA with high specificity and efficiency under isothermal conditions, and LAMP equipment is not expensive (Notomi et al., 2000;Nagamine et al., 2002). Therefore, it would be of great value to develop a rapid sex-identification system that combines LAMP technology with female-specific markers.

CONCLUSION
In conclusion, we identified sex-specific markers and SNPs via 2b-RAD-seq and a genome survey in M. nudus. These were found to belong to three female-specific scaffolds and were successfully used to accurately identify genetic sex of M. nudus. In addition, a ZW/ZZ sex-determination system was deduced for M. nudus, based on the presence of female-specific markers and the sexspecific SNPs in females appeared heterozygous but homozygous in males. Taken together, our results will help improve breeding efficiency, and provide a powerful tool for uncovering sexdetermination mechanisms in sea urchins.

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: https://www.ncbi.nlm. nih.gov/genbank/, PRJNA741812.

ETHICS STATEMENT
Ethical review and approval was not required for the animal study because sea urchin belongs to invertebrates, all animal experiments in this study are carried out with the permission of the local government.

AUTHOR CONTRIBUTIONS
ZS and YC: conceptualization and resources. ZC and JZ: methodology. ZC, JZ, and BL: investigation. ZC and CZ: data curation. ZS, ZC, and JZ: writing-original draft preparation.
ZS: writing-review and editing, funding acquisition. All authors have read and agreed to the published version of the manuscript.

FUNDING
This research was funded by the National Natural Science Foundation of China (Grant No. 31802276).