Application of SNP in Genetic Sex Identification and Effect of Estradiol on Gene Expression of Sex-Related Genes in Strongylocentrotus intermedius

Sea urchin (Strongylocentrotus intermedius) is an economically important mariculture species in Asia, and its gonads are the only edible part. The efficiency of genetic breeding in sea urchins is hampered due to the inability to distinguish gender by appearance. In this study, we first identified a sex-associated single nucleotide polymorphism (SNP) by combining type IIB endonuclease restriction site-associated DNA sequencing (2b-RAD-seq) and genome survey. Importantly, this SNP is located within spata4, a gene specifically expressed in male. Knocking down of spata4 by RNA interference (RNAi) in male individuals led to the downregulation of other conserved testis differentiation-related genes and germ cell marker genes. We also revealed that sex ratio in this validated culture population of S. intermedius is not 1:1. Moreover, after a 58-day feeding experiment with estradiol, the expression levels of several conserved genes that are related to testis differentiation, ovary differentiation, and estrogen metabolism were dynamically changed. Taken together, our results will contribute toward improving breeding efficiency, developing sex-controlled breeding, and providing a solid base for understanding sex determination mechanisms in sea urchins.


INTRODUCTION
Sea urchin (Strongylocentrotus intermedius) belongs to Echinodermata and is an important commercial species mainly distributed in several Asian countries, such as China, Korea, and Japan (1). It is well known that gonads of sea urchins are rich in unsaturated fatty acids and are the only edible parts for people (2). It has been reported that gender is one of the strongest influencing factors on the commercial value of sea urchin; for example, increasing rates of body size and gonadal growth in male sea urchins have advantages than females (3). In addition, the immunocompetence of female Paracentrotus lividus urchins might be superior to that of males (4), and the contents of free amino acids and lipid concentrations are different based on gender. As previously described, the gonads of sea urchins are classified into four stages, namely, intergametogenesis and nutritive phagocyte (NP) phagocytosis (stage 1), pregametogenesis and NP renewal (stage 2), gametogenesis and NP utilization (stage 3), and the end of gametogenesis, NP exhaustion, and spawning (stage 4) (5). Recently, three female-specific markers were identified and could be used to identify the genetic sex of Mesocentrotus nudus (6). In addition, the number of genes and micro-RNAs associated with sex determination and sex differentiation have been identified by RNA sequencing (7,8). However, no literature about sex control was reported, and a large part about sex determination and sex differentiation in sea urchin is still unknown.
With the development of next-generation sequencing (NGS) technologies, restriction site-associated DNA sequencing (RADseq) was developed to identify sex-specific markers and discover the mechanisms of sex determination in various species (9)(10)(11). In Mastacembelus armatus, two male-specific DNA fragments and two single nucleotide polymorphisms (SNPs) were identified by RAD-seq, and the Tcl gene was regarded as a candidate sex determination gene (12). In crab species, sex-specific SNPs were identified via RAD-seq, and a WZ/ZZ sex determination system was suggested based on male homogamety and female heterogamety in those sex-specific SNP loci (13). Sea cucumber (Apostichopus japonicus) also belongs to Echinodermata, a malespecific DNA fragment was verified, and a rapid genetic identifying method was developed by sex-specific primer design coupled with PCR amplification (14). Therefore, sexspecific molecular markers obtained by RAD-seq will be helpful in revealing the genetic basis of sex determination and improving breeding efficiency because genetic sex could be identified at early developmental stages.
Vertebrate steroid hormones play important roles in regulating sex differentiation and maintaining secondary sexual characteristics, especially sex steroid hormones including estradiol-17b (E2), testosterone, and progesterone (15,16). In several aquatic vertebrates, the hormonal-induced sex reversal technique had been widely used to achieve monosex culture; thus, the quality, output, and economic value had been dramatically increased (17,18). For example, XY physiological female yellow catfish was obtained by E2 treatment, then YY super-males will be obtained by gynogenesis of XY physiological female fish, and finally all-male yellow catfish will be obtained by crossing YY super-males with normal XX males (19). What is more, the YY super-males, XY males, and XX females can be identified by Y chromosome-linked and X chromosome-linked markers (20). It has been reported that vertebrate-type steroids can be measured in some invertebrates, such as in Scylla serrata (21) and Chlamys farreri (22). However, it is still unclear whether those detected hormones are endogenously synthesized or picked up from the environment (23,24). For example, mollusks can readily absorb vertebrate steroids from the environment and retain them for several weeks (25); thus, whether mollusks contain endogenously synthesized steroids is highly controversial (26). Moreover, steroid androgen exposure has a significant influence on vitellogenin production, ovarian development, and gene expression in some invertebrates, such as Portunus trituberculatus (27) and Saccostrea glomerata (28). In sea urchins, treatment with estradiol dipropionate will change the rate of protein synthesis in oocytes, increase gonad development, accelerate vitellogenesis, and induce maturation of the germ cells (29). Besides, feeding E2 will increase the rate of ovarian growth but inhibited oocyte growth in Lytechinus variegatus (30). However, it is still a mystery whether vertebrate-type steroids can be measured in S. intermedius and whether they have potential effects on gonadal development and gene expression patterns related to sex differentiation.
Spermatogenesis associated 4 gene (spata4) was first identified in human testes and involved in cryptorchidism development (31). Subsequently, the spata4 gene had been cloned from diverse organisms. In Mammalia, SPATA4 is specifically expressed in the testis and may play a significant role in maintaining spermatogenesis ability during adolescence (31,32). In chicken, the spata4 gene is also specifically expressed in the testis, and along with development progress of the testis, its expression level is gradually upregulated (33). In zebrafish and rainbow trout (Oncorhynchus mykiss), the spata4 gene is expressed specifically in the testis and slightly in the ovary (34,35). Moreover, the potential function of spata4 had been investigated by in vitro cell experiment, and it was revealed that regulatory factor X1 (RFX1) could bind the spata4 promoter, regulating its transcriptional activity and endogenous expression, and further regulate the proliferation of Sertoli cells (36). To date, no literature about spata4 was reported in sea urchin.
In the present study, 2b-RAD-seq technology was employed to identify sex-specific tags and SNP loci of sea urchin (S. intermedius). Subsequently, we investigated the 2b-RAD markers by PCR in 100 individuals, and a sex-related SNP locus had been identified. Based on the annotation of the sexrelated SNP, a male specifically expressed gene spata4 had been characterized. Then, the function implication of spata4 in testis development had been investigated by means of RNAi. Furthermore, we measured vertebrate-type steroids in S. intermedius and detected the effect of dietary administration of estradiol on testis development and sex-related gene expression in male individuals. Besides the practical implications for improving breeding efficiency, these findings also provide a molecular basis for understanding sex determination mechanisms in sea urchin.

Sample Preparation and Phenotypic Sex Identification
from markets in Dalian Haibao Fishery Limited Company and used to confirm the sex-specific markers. The tube feet tissues of each sea urchin were removed and stored in absolute ethanol at −20°C for genomic DNA extraction. Genomic DNA was extracted using the protocols of a marine animal tissue genomic DNA extraction kit (Tiangen, DP324).
To determine the phenotypic sex of individuals, the corresponding gonadal tissues were removed and fixed in 4% paraformaldehyde (PFA) at 4°C overnight, and then the gonadal tissues were washed three times with phosphate buffered solution (PBS) at room temperature. Subsequently, the gonadal tissues were embedded in optimal cutting temperature (O.C.T) compound and frozen sections (5 µm) were cut. After stained with hematoxylin/eosin, the gonadal structure was observed in Leica DM4B microscope to determine the phenotypic sex of each sea urchin.

Library Construction and 2b-RAD Sequencing
The 2b-RAD libraries were constructed at Qingdao OE Biotech Co., Ltd. (Qingdao, China) according to a previously report (37). In brief, 1 U BsaXI restriction enzyme was used to digest 200 ng genomic DNA from each individual (10 females and 10 males). The digestion products were ligated to adapter 1 and adapter 2. Then, a set of four primers that introduce individual specific barcodes and sequencing primers was used to amplify the ligation products by PCR. The PCR products were separated in 8% polyacrylamide gels by electrophoresis. After purification, each product was digested with SapI restriction enzyme and purified using streptavidin-coated magnetic beads, and 200 U T4 DNA ligase was added to the supernatant. After purification of the ligation products, barcodes were introduced by PCR with barcode-bearing primers. Then, the PCR products were purified using a MinElute PCR Purification Kit and pooled for sequencing using the Illumina NovaSeq 6000 PE150PE platform.

Data Filtering, Survey Genome, and Screening of Candidate Sex-Specific Markers
The raw reads were merged by Pear software (version 0.9.6) (38), and then the adaptor sequences of the merged reads were trimmed. The terminal 3-bp positions were also removed from each read. Then, low-quality reads that contain ambiguous bases (N) exceeding 8% or without restriction sites were excluded from the raw data and removed from each output read, and clean reads were obtained. The genomic information of S. intermedius is currently unavailable. Also, the 2b-RAD tags we obtained were too short to design PCR primers. One male and one female individual were thus selected to perform short-read de novo sequencing. For the genome survey, DNA was extracted from one female and one male S. intermedius for library construction by using TruSeq DNA LT Sample Prep kit. Briefly, the genomic DNA was sheared into fragments with length~350 bp using S220 Focused-ultrasonicators (Covaris, USA). Adapters were ligated onto the 3′ end of the sheared fragments. After PCR amplification and purification, the final libraries were sequenced on the Illumina NovaSeq platform (Illumina Inc., USA) and 150 bp paired-end reads were generated. The above 150-bp paired-end reads were clustered to build reference sequences for further locus genotyping using ustacks software version 1.34 (39). To identify candidate sex-specific markers, high-quality reads were aligned to the genome reference using SOAP2 (version 2.21) with the following parameters: r = 0, M = 4, v = 2 (40). SNP loci were developed by RAD typing program with default parameter. Finally, candidate sex-specific markers for each individual were generated for further analyses.

Validation of Sex-Specific Markers
PCR amplification was used to validate each candidate sexspecific tag and SNP locus. Firstly, specific primers of candidate sex-specific markers were designed using online software (http://biotools.nubic.northwestern.edu/OligoCalc. html) according to the flank sequences (Supplementary Table  S1). Then, the PCR reactions were carried out with a total volume of 20 ml, containing 50 ng template DNA, 10 µl of 2× Taq buffer (Kang Wei, China), and 0.5 ml of each primer (10 mM), and sterile water was added to the final volume. The PCR programs were as follows: predenaturation at 94°C for 2 min, then 30 cycles of 94°C for 30 s, annealing temperature of primer pair for 30 s, followed by 72°C for 60 s, and a final extension of 10 min at 72°C. The PCR products were separated via 1% agarose gel electrophoresis, and the expected PCR fragments were directly sequenced in both directions by Tsingke Biotechnology (Tianjin, China).

RNA Extraction, Rapid Amplification of cDNA End, and Real-Time Quantitative PCR
Total RNA was extracted from different tissues, including the ovary, testis, tube feet, and intestines, using the SV Total RNA Isolation System (Promega Z3100). Total RNA (1 µg) of the testis was used to establish a cDNA library using the SMARTer ™ Rapid Amplification of cDNA End (RACE) cDNA Amplification Kit (Clontech 634923) according to the recommendations and protocols of the manufacturer. Meanwhile, the first-strand cDNAs of each sample were synthesized according to the protocols of RevertAid First Strand cDNA Synthesis Kit (Fermentas). Real-time quantitative PCR (RT-qPCR) was performed on the LightCycler ® 96 Instrument (Roche) according to a previous description (41), and 18S ribosomal RNA (18S rRNA) was used as the internal control. The relative expression of the target genes was calculated with the 2 −DDCt method. For statistical analysis, one-way ANOVA was calculated with SPSS software and a probability (p) of <0.05 was considered statistically significant. All primers used were designed by online software http://biotools.nubic.northwestern.edu/OligoCalc.html (Supplementary Table S1).

RNA Interference
The gene-specific dsRNA of spata4 was designed by https://www. dkfz.de/signaling/e-rnai3/ and synthesized by the T7 RiboMAX ™ Express RNAi System (Promega) according to the protocol. Firstly, the genetic sex of sea urchins was identified by PCR reaction and Sanger sequencing depended on the sexspecific SNP marker. Then, a total of nine male S. intermedius with a mean shell diameter of 52.5 ± 2.10 mm was selected to inject dsRNA. Meanwhile, the testis development of selected S. intermedius was during gametogenesis and NP utilization (stage 3), and this is a critical stage in spermatogenesis. Before injection, three sea urchins were randomly removed and their gonads were surgically sampled as control samples. Then, spata4-specific dsRNA was heated to 95°C for 3min and cooled to room temperature (25-30°C). Subsequently, sea urchin was injected with 100 µl PBS containing 100 µg dsRNA according to a previous report (5). On day 3 after the first injection, three sea urchins were randomly removed and their gonads were surgically sampled. Then, a second injection was performed on the remaining sea urchins. On day 7 after the first injection, the gonadal tissues of remaining three sea urchins were surgically sampled. Histological examination of the gonads was carried out according to the method above. RT-qPCR was performed to investigate the knockdown effect and the dynamic expression change of other conserved sex-related genes.

Enzyme-Linked Immunosorbent Assay and Dietary Administration of Estradiol
Strongylocentrotus intermedius used for this experiment were purchased from markets in Dalian Haibao Fishery Limited Company, and after being acclimatized to laboratory conditions for 2 weeks, E2 concentration in the gonadal tissues was detected via a commercial enzyme-linked immunosorbent assay (ELISA) kit (Fish estradiol ELISA kit, Mlbio) according to the protocol of the manufacturer. A total of 54 sea urchins individuals with a mean diameter of 31.3 ± 0.43 mm were randomly divided into three groups and cultured in separate seawater tanks at 15°C-17°C. In addition, we made little plastic frames and put them in the seawater tanks, and each sea urchin was cultured and fed on the little plastic frame independently. Among the three groups, two experimental groups had diets with different estradiol contents, one group was fed with food mixed with 50 mg/kg E2 and another group was fed with food mixed with 100 mg/kg E2, while the control group was fed without E2. E2 was purchased from Sigma-Aldrich (E8875-250). Feeding experiments lasted 58 days, and sea urchins were fed with dietary supplementation of E2 every day. On days 3, the excrement of sea urchins was collected to detect E2 level. In addition, on days 15, 30, and 58, six individuals were randomly removed from each group, and their gonads were sampled to analyze the E2 level and for histological examination.
Transcriptome Sequencing, Assembly, Annotation, and Discovery of Differentially Expressed Unigenes On day 58 after the feeding experiments, 12 gonadal samples from the control group and those fed with 100 mg/kg E2 were selected for RNA extraction and transcriptome sequencing. For transcriptome sequencing, every three gonadal samples from the same group were mixed as one sequencing sample and each group had two sequencing samples. The transcriptome sequencing was performed on the Illumina sequencing platform (HiSeqTM 2500 or Illumina HiSeq X Ten) by Qingdao OE Biotech Co., Ltd. (Qingdao, China) according to a previous report (42). The clean reads were obtained after being trimmed off of reads with ploy-N and of low quality. To assemble the clean reads, de novo assemblies were carried out with Trinity (version: trinityrnaseq_r20131110) software. Subsequently, the assembled unigenes were used to perform the BLAST program on the databases of NCBI RefSeq (NR), Swiss-Prot, and KOG protein databases. In order to identify the differentially expressed genes (DEGs), FPKM and read count value of each annotation unigene was calculated using bowtie2 and eXpress, and DEGs were identified using the DESeq functions estimateSizeFactors and nbinomTest with the following parameters: p < 0.05 and fold change >2 or fold change <0.5. Then, DEGs were assigned to different Gene Ontology (GO) terms using the Blast2GO software. Meanwhile, DEGs were blasted against the Kyoto Encyclopedia of Genes and Genomes (KEGG) database by using the automated assignment server (KAAS, http://www. genome.jp/kaas-bin/kaas_main) to identify the biological pathways of these unigenes.

Multiple Sequence Alignment and Phylogenetic Analyses
The amino-acid sequence was predicted by using the DNAMAN software. Multiple sequence alignment was performed using ClustalW. The phylogenetic trees were constructed using the neighbor-joining method in MEGA version 7.0, with 1,000 replicates of bootstrap analysis. The sequences used in this study were downloaded from the NCBI database (Supplementary Table S2).

2b-RAD-Seq Data Analysis
In total, the 2b-RAD-seq produced 140,265,655 raw reads (GenBank: PRJNA761225), and we obtained 131,644,844 clean reads after filtering low-quality reads. In order to cluster and assemble these clean reads, two female and two male samples were selected as references based on read quality and quantity. Following alignment with reference samples, we obtained 979,411 putative tags with an average depth of 33.75× from all 10 female and 10 male samples ( Table 1). In addition, 43,922 putative SNP loci were screened by alignment analysis.

Genome Survey and Validation of Candidate Sex-Specific Tags and SNPs
A total of 102.37G clean reads with a GC content of 37.16% and 99.94G clean reads with a GC content of 36.75% were acquired from the female and male samples, respectively (GenBank: PRJNA762987). The corresponding draft genomes were then assembled for both female and male. In order to identify sex-specific tag candidates, the 2b-RAD tags of each sample were aligned with both the female and male genome survey. We obtained 38 male-specific candidates that present exclusively in 10 male samples and 13 female-specific tags that only exist in 10 female samples (Supplementary Table S3, Additional Fasta File 1). Since a smaller p-value of SNP means a higher correlation with sex, the top 10 sex-specific SNP loci were screened with the lowest p-values for further analysis. Notably, the females appeared to be heterozygous, whereas the males seemed to be homozygous for all 10 sex-specific SNP loci candidates in all samples we analyzed (Supplementary Table S4, Additional Fasta File 2), indicating that S. intermedius may have an ZW/ ZZ sex determination system.
Among the 51 candidate sex-specific tags identified initially, six were too short to be analyzed and were thus excluded. These were female-tag 2, female-tag 5, male-tag 6, male-tag 12, maletag 13, and male-tag 16. To confirm the authenticity of the remaining 45 tags, we designed primers based on the flanking sequences of these tags and performed PCR verification in 10 female and 10 male sea urchins. Unfortunately, we obtained the same PCR products from all 45 tags in both male and female sea urchins. Sanger sequencing analysis of the PCR products suggested that no difference was observed between male and female. Therefore, all these candidate sex-specific tag candidates were excluded from further study. Next, the same method was used to assess the authenticity of the 10 candidate SNP loci. PCR amplification suggested that nine of them did not show any difference between male and female. Hence, these nine sexrelated SNP loci candidates were also excluded for any further analysis. Nonetheless, Sanger sequencing of the PCR products indeed confirmed one candidate sex-specific SNP locus (ref128486). Interestingly, all male individuals we analyzed displayed an apparently homozygous type (G/G), whereas the female individuals showed a G/A heterozygous type ( Figure 1).
Together, these data suggest that this sex-specific SNP may be used to identify the genetic sex of S. intermedius.

Genetic Sex Validation of Strongylocentrotus intermedius in Breeding Populations
To assess whether the above sex-specific SNP locus can be used to characterize the phenotypic sex of sea urchin, we ran PCR reactions with the same primers in 100 sea urchins. Prior to PCR amplification, gonadal histology analyses were performed to determine the phenotypic sex of each sea urchin and 65 females and 35 males were identified (Supplementary Figure S1). Among all 65 phenotypic female individuals, 56 (86.15%) exhibited the G/ A heterozygous genotype (the dominant female genotype) and 9 (16.07%) exhibited the A/A homozygous genotype at the SNP locus. For the rest of the phenotypic male individuals, 33 (94.29%) exhibited the G/G homozygous genotype (the dominant male genotype) and 2 (5.71%) exhibited the G/A heterozygous genotype at the SNP locus ( Table 2). Despite the fact that high consistency (89%) was obtained rapidly from the PCR amplification, the PCR products from 11 individuals (11%) were inconsistent with their phonotypic sex. Importantly, the segregation ratio of female to male was 65:35 (1.86:1) in this culture population. The above data indicate that sex reversal is possible in S. intermedius.

Annotation of Sex-Associated SNP Sequence
The sex-associated SNP locus is located on scaffold ref128486. Ref128486 is 6,149 bp in length. To search for homologous genes on the scaffold, we performed a BLAST program on NCBI databases. Interestingly, a spermatogenesis-associated (SPATA) family gene (spata4) was identified, and it has been reported that SPATA family genes play important roles in spermatogenesis, sperm maturation, and fertilization (43). The spata4 cDNA is  Figure S2). Phylogenetic analysis showed that the topology of clades is basically consistent with the known taxonomic relationship among these analyzed species ( Figure 2B). In addition, the mRNA expression pattern of spata4 in S. intermedius was detected by RT-qPCR in different adult tissues, including the testis, ovary, tube feet, and intestines. The results suggested that spata4 is expressed predominantly in the testis. By contrast, only a few of the spata4 transcripts can be detected in the ovary and intestine (Figure 3).

Knockdown of spata4 by RNAi in Strongylocentrotus intermedius
Considering that spata4 is expressed predominantly in the testis (Figure 3), we further explored its function in spermatogenesis by RNAi, a strategy that has been successfully utilized to study sex-related gene function in sea urchin (5). The genetic sex of nine sea urchins were first identified by PCR amplification and the corresponding products were sequenced. Gonadal tissues from three male sea urchins were removed and sampled as a control group before injection. The specific spata4-targeting double-stranded RNA (dsRNA) was then injected into adult S. intermedius. At 3 and 7 days post injection (dpi), gonadal samples were collected to perform histological analysis and to detect the spata4 mRNA expression with RT-qPCR. We observed that, compared with the control groups, the expression levels of spata4 in the testis were significantly decreased (about 60%) in knockdown samples ( Figure 4A). Moreover, the results of histological examination showed that no obvious apoptotic cells were observed in the dsRNA-treated sea urchins (Supplementary Figure S3). We also measured the expression levels of other well-studied genes involved in testis differentiation, ovary differentiation, and germ cell development. As shown in Figure 4A, the expression of dmrt1 and soxE in the testis differentiation signaling pathway significantly decreased in spata4-knockdown samples compared with the control groups at 3 and 7 dpi. Interestingly, the transcripts of foxl2 and hsd17b8 were similarly suppressed significantly in spata4-knockdown samples at 3 dpi, but the hsd17b8 expression increased about 2.74-fold in the knockdown samples at 7 dpi ( Figure 4B). In addition, the germ cell marker genes, including boule and nanos1 expression, were also significantly decreased after injection in knockdown samples ( Figure 4C).

Dietary Administration of Estradiol in Male Strongylocentrotus intermedius
Invertebrates show a high degree of variability in sex, and phenotypic sex may be affected by environmental factors, such as sex hormone, temperature, and stocking density (44). Because the genetic sex of 11 individuals identified by PCR amplification was inconsistent with their phonotypic sex and the sex ratio of females to males was not 1:1, we speculated that sex reversal is possible in S. intermedius. To test this hypothesis, we first determined the level of E2 in the gonads of S. intermedius by ELISA and found that the concentration of E2 was 17.4 and 13.75 pmol/L in the ovary and testis, respectively. We then characterized male sea urchins with the above SNP and feed them with formulated diet supplemented with 50 or 100 mg/kg E2 (defined as group 1 and group 2, respectively) for 58 days. At 3 days post feeding (dpf), the excrement was collected to detect E2 level. As shown in Figure 5A, the E2 c o n c e n t r a t i o n i n c r e a s e d a l o n g w i t h t h e d i e t a r y supplementation of E2. Gonadal tissues were also collected at 15, 30, and 58 dpf to examine the concentration of E2. As shown in Figure 5B, the E2 concentrations from sea urchins fed with E2 were comparable with those of the control group during a period of 15 days. On day 30, the concentrations of E2 in group 2 were 1.16 times higher than those of the control group. On day 58, the E2 concentrations from group 1 and group 2 were 1.20-and 1.20-fold against the control group, respectively ( Figure 5B). These data suggested that S. intermedius can absorb estradiol from food and accumulate in the gonads.  To explore the dynamic expression changes of sex-related genes after E2 treatment, we performed comparative transcriptome analysis to investigate the expression characteristics and differentially expressed genes (DEGs) in the control group and in the 100-mg/kg E2-treated sea urchins (i.e., group 2). In total, 96.13M and 94.57M clean reads were acquired from group 2 and the control group libraries, respectively (GenBank: PRJNA761244). After assembly with de novo methods and functional annotation, 67,492 unigenes were identified. The expression levels of unigenes were then compared, and 2,355 DEGs with at least two-fold changes (p < 0.05) were identified. Of these DEGs, 1,050 were downregulated and 1,305 were upregulated in the E2 treatment group. In order to categorize the transcripts according to putative function, GO assignment was performed on all DEGs. Not surprisingly, several upregulated genes were enriched significantly in GO terms that are related to estrogen metabolism, spermatogenesis, and oogenesis ( Figure 6A). These include development of primary sexual characteristics (GO:0045137), sperm-egg recognition (GO:0035036), ER-nucleus signaling pathway (GO:0006984), and sexual reproduction (GO:0019953). By contrast, many downregulated genes were significantly enriched in terms of germ plasm (GO:0060293) and positive regulation of fibroblast apoptotic process (GO:2000271) ( Figure 6B). Searching the downregulated unigenes in the KOG database showed significant enrichment of unigenes in pathways previously implicated in ovarian steroidogenesis (ko04913), drug metabolism (ko00982), and cytochrome P450 estrogen signaling pathway (ko04915). Besides, the pathway of oocyte meiosis (ko04114) was recognized as a significantly upregulated term. Together, these data suggested that E2 treatment leads to downregulation of testis-biased genes but upregulation of ovarybiased genes and genes associated with estrogen metabolism pathway in S. intermedius. Lastly, we analyzed the dynamic expression patterns of several sex-related markers that were reported to have been involved in gonad development in the testis of E2-treated S. intermedius by RT-qPCR. As shown in Figure 7A, testis differentiation-associated genes, including dmrt1, soxE, and sp17, were significantly downregulated in the testis samples of E2-treated S. intermedius. Curiously, genes needed for ovary differentiation including hsd17b8, foxl2, and had17b14 were also downregulated ( Figure 7B). Notably, E2 treatment led to a significant decrease in germ cell marker genes including nanos2, boule, seali, sycp3, and vasa in sea urchins ( Figure 7C). In addition, steroidgenic enzyme genes related to active androgen biosynthesis such as cyp17a, cyp450, and nek1 were also downregulated in the testis of E2treated sea urchins ( Figure 7D). As expected, genes associated with estrogen metabolism including ER membrane protein complex subunit 4 (Emc4) were significantly upregulated ( Figure 7E). However, no phenome of sex reversal was observed in the testis of E2-treated S. intermedius according to histological analysis (Supplementary Figure S4).

DISCUSSION
Identification of sex-specific markers is fundamental to understand the sex determination mechanism and is crucial for genetic sex identification and sex-controlled breeding (45,46). Recently, 2b-RAD-seq has been successfully applied to identify sex-specific markers in a great diversity of aquatic organisms (9,10). RNA interference is an effective genetic tool for studying gene function in vivo, especially in animals with low survival rates and long reproductive cycles (5). Currently, a green fluorescent protein (GFP)-derived double-stranded RNA (dsRNA-GFP) is usually used as control in RNAi studies, and no RNAi effect was detected 18S was used as the control. Each bar represents mean ± standard deviation (SD) (n = 3). One-way ANOVA was used to determine statistical analysis. Asterisks (*) indicate significant differences (p ≤ 0.05).
using dsRNA-GFP suggesting target sequence specificity of the RNAi effect (47). However, according to a previous report in Apis mellifera, the expression of nearly 1,400 genes was dynamically changed in response to dsRNA-GFP (48). Hence, it is indeed a big challenge to evaluate target sequence specificity of interference effect.
In pufferfish, a female homozygous (C/C) type and a male heterozygous (C/G) SNP in anti-Mullerian hormone receptor type 2 (amhr2) were identified, and this SNP site can be used to identify the genotypic sex of pufferfish with high accuracy (49). Notably, the missense SNP that resulted in an amino acid mutation (His/Asp) in amhr2 affects its kinase activity, and the combination of the two alleles of amhr2 may determine the sex of Takifugu rubripes (50,51). In the current study, a sex-linked SNP locus was identified in S. intermedius spata4 gene. The SNP locus is located at the 175th codon of the spata4 gene, but the CAC➔CAU change only causes a synonymous mutation in the predicted Spata4 protein. Synonymous mutations are thought to be silent for a long time. Recently, more and more studies reveal that synonymous sites do not evolve neutrally, and synonymous mutations can regulate protein conformation and function by affecting mRNA stability (52,53). However, the effect of synonymous mutations on the expression and function of spata4 gene remains unclear and should be investigated in the future. In most of the species, spata4 is specifically expressed in the testis and may play important roles in maintaining spermatogenesis ability and infertility process (31,43,54). Herein, we revealed that spata4 is expressed predominantly in the testis, and only a few of spata4 transcripts can be detected in the ovary ( Figure 3). As previously reported, dmrt1 and soxE were demonstrated to play critical roles in testis differentiation regulatory pathway (55,56). Meanwhile, boule and nanos1 were previously described as germ cell markers and played conserved roles in regulating germ cell development (57,58). Knocking down spata4 expression in sea urchins led to the downregulation of soxE, dmrt1, boule, and nanos1, indicating a close correlation between spata4 expression and testis development.
Although estradiol can be detected in the gonadal tissues of S. intermedius, we are uncertain whether it is absorbed from the environment or is endogenously synthesized. It has been reported that vertebrate steroids can be absorbed and utilized by invertebrates, especially those living in water, and all reproductively mature fish can release steroids into the environment (23,25,59). Sea urchins (L. variegatus, S. purpuratus, and S. intermedius) are also capable of absorbing vertebrate steroids. A notable example is that estradiol-17b exposure can induce protein synthesis in oocytes of S. intermedius (29).
It has been suggested that sex determination and sex differentiation of invertebrates may be influenced by environment factors such as temperature, nutrition level, and photoperiod (44). Because the ratio between female and male was not around 1:1 in the current study, we speculated that sex reversal might be the case in S. intermedius. After 58 days of dietary administration of estradiol, we did not observe significant sex reversal in male individuals by gonadal histology analysis. However, it should be noted that the expression levels of multiple genes have changed dynamically after estradiol treatment. In particular, many of these genes are known to be involved in testis differentiation, ovary differentiation, and germ cell development. Thus, it is possible that E2 treatment may affect gonadal development in sea urchin.
A vertebrate homologous estrogen receptor (ER) (ADL71443.1) is annotated in transcription database. Nonetheless, its expression level is relatively stable following E2 treatment, implying that the interacting mechanisms between exogenous estrogens and sea urchins may be different from that of vertebrates. In Marisa cornuarietis, estradiol-17b treatment has no obvious impact on the expression levels of ER-like or estrogen-related receptor (ERR) gene. This is not surprising because estradiol-17b does not bind to ER-like or ERR (60,61). Hence, it is attempted to reveal the interactions between estrogen and estrogen receptor in sea urchin.
In conclusion, the current studies identified a sex-related SNP site in S. intermedius and revealed that this SNP is located within the spata4 gene. We confirmed that spata4 is expressed predominantly in the testis, and spata4 may be involved in testis development. In addition, S. intermedius can absorb vertebrate estradiol from dietary food, which caused dynamic changes of the sex-related genes. Taken together, our results provide a powerful tool for understanding the sex-determination mechanisms in sea urchins and have important implications for improving the breeding efficiency of sea urchins.

DATA AVAILABILITY STATEMENT
The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/sra/?term= PRJNA761225, https://www.ncbi.nlm.nih.gov/sra/?term= PRJNA762987, https://www.ncbi.nlm.nih.gov/sra/?term= PRJNA761244. You can also get raw datas by searching the accession numbers (PRJNA761225, PRJNA761244, PRJNA762987) in NCBI. The accession number of spata4 sequence is OK078897. This sequence will be accessible after the indicated release date or the article was published.