The published genome (taxid: 307972) data of the sea cucumber at the scaffold level (Zhang et al., 2017) was downloaded from the NCBI database (https://www.ncbi.nlm.nih.gov/). A BLASTP search was then performed using the common conserved domain protein sequence of the DM domain (XM_030995783.1) and the High Mobility Group box (HMG box) (XP_786809.2) to screen the homologous genes of the Dmrt family and SoxE, respectively. The E-value was set to ≤ e− 5. The suspected members of the Dmrt gene family and SoxE gene were validated by blasting against the NCBI database, in order to assess the reliability of the analysis.

The coding sequence of the Dmrt gene was predicted by ORF Finder online software (https://www.ncbi.nlm.nih.gov/orffinder). BioEdit software was used to perform multiple sequence alignments and visual analysis of the Dmrt proteins in the sea cucumber. The sequences used in this study were downloaded from the NCBI database (Supplementary Table S1). Phylogenetic analysis was performed using MEGA-7 software with bootstrapping (1,000 replicates) and the neighbour-joining (NJ) method.

Adult sea cucumbers were collected from the coastal areas of Dalian, China. No endangered or protected species were involved in this study. The expression patterns of the Dmrt genes and SoxE gene in different tissues were analysed by RT-qPCR. Briefly, total RNA from different tissues, including the tube feet, testes, intestines, ovaries, stomach, longitudinal muscle, coelomocytes and respiratory tree, were extracted using the Silica Membrane-based Vacuum Pump (SV) Total RNA Isolation System (Promega Z3100). The middle part of the intestine was used for gene expression analysis. According to previously reports (Sun et al., 2021), the gonads of A. japonicus are classified into four stages, including the early growing stage (Stage 1), growing stage (Stage 2), mature stage (Stage 3) and post-spawning stage (Stage 4). After gonadal histology analysis, gonadal tissues at different development stages were collected to examine the dynamic expression changes in Dmrt1 and SoxE. Gonadal tissues at Stage 2 were used to detect tissue expression specificity.

RT-qPCR experiments were performed in 20-μl reactions containing 10 µL of Fast Start Essential DNA Green Master (Roche, Mannheim, Germany), 6.4 µL of sterile water, 0.8 µL of each 10 mM primer, and 2 µL of cDNA. The protocol was as follows: 95°C (10 min) for heat denaturing; then, 40 cycles of 95°C (15 s) and 60°C (1 min). The efficiency of primers discovery through standard curve formulation reached 98%. The housekeeping gene NADH was selected as the internal reference gene (Sun et al., 2021). Data were performed from three independent experiments. Each sample was analyzed in triplicates, and the relative expression levels of target genes were calculated with the 2^—ΔΔCt method. For statistical analysis, one-way ANOVA was calculated with SPSS software after the normal distribution and homogeneity of variance test (SPSS Inc.), and a probability (p) of ≤0.05 was considered statistically significant. The primers used in this study were designed using online software (http://biotools.nubic.northwestern.edu/OligoCalc.html) (Supplementary Table S2).

Gene-specific dsRNAs for the Dmrt1 and SoxE genes were designed using an online platform (https://www.dkfz.de/signaling/e-rnai3/). DsRNAs were synthesised using a T7 RiboMAX™ Express RNAi System (Promega), according to the manufacturer’s protocol. Three different dsRNAs were designed for each gene to ensure interference efficiency. The genetic sex of each sea cucumber was determined according to a previous report (Wei et al., 2021). In briefly, the tube feet tissues were used to extract genomic DNA by alkaline lysis method, then male-specific primers were used to amplify the male specific DNA marker, one specific band was amplified in males, but not in females. For the Dmrt1 RNAi experiment, 24 male sea cucumbers with an average weight of 100 ± 20 g were selected and were randomly divided into two groups: Dmrt1-knockdown (n = 12) and control (n = 12). For the SoxE RNAi experiment, 24 male sea cucumbers with an average weight of 60 ± 10 g were selected. The 24 males were divided into two groups: SoxE-knockdown (n = 12) and control group (n = 12). The injection experiments were performed as previously reported with slight modification (Sun et al., 2021). Specifically, before injecting, three gene-specific dsRNAs were mixed, then 100 μg of gene-specific dsRNA per 50 g body weight was injected into each sea cucumber in the knockdown group. Meanwhile, sea cucumbers in the control group were injected with GFP dsRNA at the same concentration. For the Dmrt1 RNAi experiment, three sea cucumbers from each group were removed randomly on the third day, seventh day and 10th day following the first injection. The gonads were surgically sampled for future research. In addition, repeat injections were performed on the third, sixth and ninth days following the first injection. For the SoxE RNAi experiment, three sea cucumbers from each group were removed randomly on the third, 7th and 23rd day following the first injection. The gonads were surgically sampled for future research. In addition, repeat injections were performed every 3 days after the first injection.

Gonadal tissues were surgically removed and fixed in 4% paraformaldehyde (PFA) at 4°C overnight. Then, they were washed three times in phosphate-buffered solution (PBS) and balanced in 30% sucrose at room temperature for 2 h. Next, the treated tissues were embedded in optimal cutting temperature (O.C.T) compound and frozen sections (5 μm) were cut and stained with haematoxylin/eosin. Digital photos were acquired under a Leica DM4B microscope

Genomic information for A. japonicus is available from NCBI (assembly ASM275485v1) (Zhang et al., 2017). To identify the Dmrt and SoxE genes in A. japonicus, the protein sequences of the conserved DM domain and HMG box of S. purpuratus (XM_030995783.1 and XP_786809.2) were used as the queries to perform a BLASTP search, respectively. Five predicted Dmrt genes (PIK43860.1, PIK34621.1, PIK44057.1, PIK33706.1 and PIK41536.1) were identified. Each Dmrt gene contained a complete coding DNA sequence (CDS) with lengths of 1,029, 651, 945, 507 and 849 bp, respectively. All identified Dmrt proteins had a DM domain (Supplementary Figure S1A), but had a different number of exons. PIK41536.1, PIK43860.1 and PIK34621.1 contained just a single intron, PIK44057.1 had two introns and PIK33706.1 had no introns (Supplementary Figure S1B). The identity between Dmrt cDNA sequences is about 32 %. Likewise, a predicted SoxE gene encoding a protein of 459 amino acids was identified AjSoxE. As expected, there was a 71-amino acid HMG box domain in the protein sequence (Supplementary Figure S1C).

To define these five Dmrt genes, 40 Dmrt genes from 17 species were downloaded from the NCBI database and multiple protein sequence alignments and phylogenetic analyses were performed. The protein sequences of the Dmrt genes exhibiting low similarity between species (4.94%–57.01%). However, the protein sequence of DM domain was highly conserved such that the similarity with other species was generally more than 60% (Supplementary Figures S2–S5). The phylogenetic analyses indicated that the Dmrt family genes can be classified into five branches (Dmrt1, -2, -3, -4 and -5). PIK41536.1 and PIK43860.1 were clustered into the Dmrt1 and Dmrt2 branches, and were named AjDmrt1 and AjDmrt2, respectively. PIK34621.1 and PIK44057.1 were both clustered with Dmrt3 homologs from other species, so were named AjDmrt3a and AjDmrt3b, respectively. The deduced amino acid sequences of PIK33706.1 were clustered into the Dmrt5 branch and were defined as AjDmrt5 (Figure 1A). Surprisingly, none of the AjDmrt genes clustered with Dmrt4 homologs from other species’ branches. The topology of clades of Dmrt1, Dmrt2, Dmrt3 and Dmrt5 were basically consistent with the known taxonomic relationships among these species. Moreover, the phylogenetic analyses showed that AjSoxE was firstly grouped with SoxE in invertebrate species, then clustered into Sox8, Sox9 and Sox10 branches in jawed vertebrates (Figure 1B). This phylogenetic tree further confirms that the ancestral SoxE gene is replicated in vertebrates and has produced at least three members (Sox8, Sox9, and Sox10).

In order to reveal the expression differences, the transcripts of all five AjDmrt genes and the SoxE gene in adult tissues were examined by RT-qPCR. Specific primers of five AjDmrt genes were designed in the parts with the greatest difference in the nucleic acid sequence. AjDmrt1 was expressed predominantly in the testes, only low AjDmrt1 expression was detected in other tissues, including the ovaries, stomach, tube feet, longitudinal muscle, respiratory tree, and coelomocytes (Figure 2A). In contrast, AjDmrt2 was higher expressed in the tube feet (Figure 2B). AjDmrt3a and AjDmrt3b were expressed predominantly in the coelomocytes and tube feet, with low expression of AjDmrt3a and AjDmrt3b detected in other tissues, including the intestines, stomach, longitudinal muscle, ovaries, and respiratory tree (Figures 2C,D). The expression pattern of AjDmrt5 in adult tissue was similar to that of AjDmrt2, with predominant expression in the tube feet and slight expression in the other analysed tissues (Figure 2E). Strikingly, the SoxE gene had the highest expression in the tube feet. Its expression differed significantly between the male and female gonads, such that a large number of transcripts were detected in the ovaries while SoxE was nearly undetectable in the testes (Figure 2F).

Considering that the expression of Dmrt1 and SoxE differed significantly between the ovaries and testes in A. japonicus, gonadal tissues at different development stages were collected to examine the dynamic expression changes in Dmrt1 and SoxE. After gonadal histology analysis, gonadal tissues of four different development stages (i.e., the early growing stage (Stage 1), growing stage (Stage 2), mature stage (Stage 3) and post-spawning stage (Stage 4)) were obtained from male and female individuals, respectively. The expression of Dmrt1 was continuously maintained at a high level during the spermatogenesis process and was significantly downregulated during Stage 2 in the testes (Figure 3A). SoxE expression in the testes at the corresponding stage was much lower than that of Dmrt1, and there were no significant differences between Stage 1, Stage 2 and Stage 3. Along with the occurrence of spermatogenesis, the expression level of SoxE increased up to 6.3-fold during Stage 4 against to Stage1 (Figure 3B). Meanwhile, a large amount of SoxE transcript was detected in Stage 1 ovary tissue, but the relative expression level continuously decreased together with the oogenesis development process, reaching the lowest in Stage 4 (Figure 3C).

RNAi was further utilised to examine the underlying function of Dmrt1 in testes development in male A. japonicus. As shown in Figure 4A, the expression levels of Dmrt1 significantly decreased, down 40.85% and 47.80% in the knockdown group compared to the control group at 3 days post-injection (dpi) and 7 dpi, respectively. On the 10 dpi, Dmrt1 expression was only 3.31% of that in the control group. Compared to the control group, the transcripts of SoxE were reduced in the knockdown group by 88.65% at 3 dpi, 45.52% at 7 dpi and 77.10% at 10 dpi (Figure 4B). Moreover, the dynamic expression changes in foxl2 and piwi, which have been reported to be involved in ovary differentiation and germ cell development, respectively (Ottolenghi et al., 2005; Gonzalez et al., 2021), were investigated. Interestingly, the expression level of foxl2 was significantly downregulated (about 46.17%) in the knockdown group at 3 dpi, but sharply increased 1.67–2.17-fold at 7 dpi and 10 dpi (Figure 4C). The expression level of piwi did not significantly differ between the knockdown and control groups at 3 dpi and 7 dpi, but was significantly reduced by 58.36% at 10 dpi in knockdown adults (Figure 4D). Although the expression levels of genes involved in testis differentiation, ovary differentiation and germ cell development showed dynamic changes after knockdown of Dmrt1 by RNAi, no obvious histological changes and no apoptosis cells were observed in the knockdown testes (Supplementary Figure S6).

Since SoxE has been reported to play a conserved role in sexual differentiation, and it shows a sexually dimorphic expression pattern in the gonads of A. japonicus, its function during testis development was investigated by RNAi in an adult male individual. As shown in Figure 5A, the transcript level of SoxE was significantly decreased by 79.46% with SoxE knockdown compared to the control at 3 dpi. Its expression was continuously downregulated at 7 dpi and 23 dpi. In addition, the expression of Dmrt1 was sharply decreased to almost an undetectable level in the knockdown testes (Figure 5B). Conversely, the expression level of foxl2 was increased about 8.09-fold and 7.34-fold in the knockdown sample compared to the control sample at 3 dpi and 23 dpi, respectively. However, the expression level of foxl2 was not significantly different compared to the control at 7 dpi (Figure 5C). The transcripts of piwi were significantly decreased by 79.42%, 84.60% and 48.87% at 3 dpi, 7 dpi and 23 dpi, respectively (Figure 5D ). However, no obvious histological changes or apoptosis cells were observed in the knockdown testes (Supplementary Figure S7).