The sea cucumber, E. quinquesemita, was collected by diving at the depth of about 3 m in Tokyo Bay, near the pier of Hakkeijima Sea Paradise in Yokohama, Kanagawa Prefecture, or at the depth of about 5–10 m near Aquamarine Fukushima in Onahama, Fukushima Prefecture. Animals were kept in tanks until used for experiments. Methods for induction of evisceration and regeneration were based on Okada and Kondo, 2019. The day of evisceration was designated as 0 days post evisceration (dpe). Briefly, evisceration was induced by injection of approximately 100 μL of 0.45 M KCl into the coelom and the eviscerated animals were kept in aquaria with sea water at 13–18°C without feeding, up to 26 dpe. Embryos were collected after spontaneous spawning in laboratory tanks and kept in 2 L of filtered sea water in a plastic container at 20°C without feeding.

To identify MET-related genes, Hox and Parahox genes and neuronal patterning genes in the focal sea cucumber species, RNA sequencing was performed. Total RNAs were extracted from 2 embryonic stages (early and late gastrula) and juveniles of E. quinquesemita using QIAzol lysis reagent (Qiagen), and subsequently the RNA samples were pooled together. The pooled sample was then subjected to library preparation for RNA sequencing and sequenced on a Hiseq 4000 platform (Illumina) at a commercial service provider, Eurofins Genomics (Tokyo). The resulting 100-bp (base pairs) paired-end reads were deposited in the Sequence Read Archive of the DNA Data Bank of Japan (DDBJ) under the accession number DRA017482. Adapter and low-quality sequences were removed from the paired-end reads using Trimmomatic v0.39 (Bolger et al., 2014) with the following options: ILLUMINACLIP : TruSeq3-PE.fa:2:30:10, LEADING:20, TRAILING:20, SLIDINGWINDOW:4:20, MINLEN:25. Using the filtered reads, de novo transcriptome assembly was performed using Trinity v2.8.5 (Grabherr et al., 2011) with default options. To obtain the target orthologous genes, BLAST searches against the assembly were carried out with target gene orthologs in Strongylocentrotus purpuratus and Apostichopus japonicus (Tu et al., 2012; Zhang et al., 2017). The top-hit sequences from the E. quinquesemita transcriptome database were defined as putative orthologs of the target genes. To confirm the orthologs, reciprocal BLAST searches and phylogenetic analyses were performed ( Supplementary Figures 1 – 7 ).

To analyze gene expression patterns by quantitative PCR, total RNAs were firstly extracted from the stages during regeneration after the evisceration in the focal sea cucumber species ( Table 1 ). Since the regeneration requires about 3 weeks, samples were prepared at 0, 4, 7, 14, and 26 dpe. Regenerating digestive tracts were isolated by dissection and tissues from 2–4 individuals were pooled as each sample to obtain enough total RNAs. Separately from anterior and posterior digestive tracts, total RNAs were extracted by using QIAzol lysis reagent (Qiagen) and the RNA purification was performed with Agencourt AMPure XP (Beckman Coulter). Intact digestive tracts from mature individuals were also used for comparison. The extracted RNAs were reverse-transcribed using a High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems). Quantifications of target transcript levels were performed using Fast SYBR Green Master Mix and an ABI Prism 7500 instrument (Applied Biosystems). To evaluate endogenous control levels of constitutive expressions of putative reference genes, i.e., elongation factor 1 alpha (EF1α), tubulin beta chain (TUBB), NADH dehydrogenase (NADH) and ribosomal protein S18 (RPS18), were evaluated using geNorm (Vandesompele et al., 2002) and Normfinder (Andersen et al., 2004). The results indicated that EF1α was the most appropriate reference gene for comparisons among stages and body parts. Primers for qPCR ( Supplementary Table 1 ) were designed using Primer Express software (ver. 3.0.0, Applied Biosystems). Data acquisition and analyses were performed using ABI Prism 7500 software ver. 2.0.4 (Applied Biosystems) with the relative standard curve method. For statistics, biological triplicates were subjected to one-way analysis of variance (ANOVA) followed by Tukey’s multiple comparisons tests (p < 0.05) using R ver. 4.0.2.

To investigate the localization of target genes, whole mount in situ hybridizations (WISHs) were carried out, based on previous studies with minor modifications (Mashanov et al., 2010; Omori et al., 2011; Kikuchi et al., 2015). The PCR products were used as templates to transcribe riboprobes with a DIG RNA Labeling Kit (Roche). After dissection, the tissue samples were fixed with 4% paraformaldehyde in 0.5 M NaCl, 0.1 M 3-(N-morpholino) propanesulfonic acid (MOPS pH 7.0) overnight at 4°C. The samples were washed with phosphate buffered saline (PBS), and then decalcified with 0.5 M ethylenediaminetetraacetic acid (EDTA) (pH 8.0) in PBS for 3–4 days at room temperature. After decalcification, the samples were kept in 99.5% ethanol at –20°C until use. The samples were washed in PBS containing 0.1% Tween-20 (PBST), treated with 1/50 volume of proteinase K (TaKaRa) in PBST for 20 min at 37°C, acetylated sequentially in 0.25% and 0.5% acetic anhydride in 0.1 M triethanolamine, 5 min each. After the samples were washed twice with PBST for 5 min each, prehybridization was performed at 58°C for 2 h or longer in hybridization buffer containing 50% formamide, 5× saline-sodium citrate (SSC), 100 μg/mL yeast RNA, 5× Denhardt’s solution, 0.1% Tween 20. The riboprobes were diluted in hybridization buffer at 58°C to a final concentration of about 400 ng/ml and denatured at 80°C for 5 min. The hybridization was carried out at 58°C overnight. After hybridization, the samples were washed with 50% formamide in 5× SSC at 58°C for 20 min, 5× SSC at 58°C for 50 min, 2× SSC at 58°C for 50 min, 0.1× SSC at 58°C for 15 min and then PBST at room temperature for 15 min twice. Subsequently, the samples were incubated in 0.1% blocking reagent (Roche) in PBST (blocking buffer) at room temperature for 30 min, followed by an incubation in 1/2000 volume of anti-DIG-AP (Roche) in blocking buffer at 4°C overnight. Following a wash with PBST (15 min, 8 times), immunodetection was performed using BM purple (Roche) at room temperature. After detection, the samples were washed with PBST and then kept in 10% formalin/PBS.

Animals were anesthetized in 72 g/L MgCl₂ in sea water for approximately 15 minutes to 1 hour before fixation. An individual body was dissected along the anterior-posterior (oral-aboral) axis of the body at the right ventral interambulacral zone, exposing the whole-body cavity, and dissected into anterior and posterior parts. The dissected bodies were fixed in Bouin’s Fixative or 4% paraformaldehyde in 0.5 M NaCl, 0.1 M 3-(N-morpholino) propanesulfonic acid (MOPS pH 7.0) overnight at 4°C. The specimens were washed with PBS, and then decalcified with 0. 5M EDTA (pH 8.0) in PBS for 3–4 days at room temperature. After decalcification, the specimens were kept in 99.5% ethanol at –20°C until embedding in paraffin (Paraplast X-TRA, SIGMA) for histological observation. Serial sections of 6- μm thickness were produced and were stained with hematoxylin-eosin (HE) or toluidine blue (TB) and observed under a light microscope (ECLIPSE TE300, Nikon). Sections were mostly made perpendicular to the anterior-posterior axis of the body. In total, 4 animals were used for histological observations.

RNA-sequencing analysis of the transcriptome derived from gastrula embryos and juveniles of E. quinquesemita yielded 28,686,525 pairs of 100-bp reads. After filtering, 26,497,155 paired reads were retained. De novo transcriptome assembly of the sequence reads generated 350,742 contigs with a total of 296,553,944 bases and N50 of 1,513 bp. Using the transcriptome data, a gene database was constructed to search target genes, i.e., MET-related genes, Hox genes and neural patterning genes. As a result, the gene orthologs from the focal sea cucumber species were obtained by reciprocal BLAST searches and by constructing phylogenetic trees ( Supplementary Figures 1 – 3 ). Based on the obtained gene sequences from the RNA-sequencing result, primers for qPCR were designed ( Supplementary Table 1 ).

Firstly, the expression patterns of MET-related genes, i.e., twist and snail, during the regeneration after anterior evisceration in E. quinquesemita were investigated. The results clearly showed that both of these genes were highly expressed in the anterior tissues of the regenerating digestive tract, in comparison with the posterior ones ( Figure 2B ). At 0 dpe, the expression level of twist in the anterior part was still low, but the twist expression became highly upregulated at 4 dpe, and was gradually decreased at the later stages. In the case of snail, the expression level was gradually upregulated, reaching the highest expression level at around 14–26 dpe, at the time when the twist expression started to decrease. The results of whole mount in situ hybridization showed that twist was expressed in anterior mesenteries ( Figure 2C ). On the other hand, in situ signals were not detected for snail.

Anterior Hox genes such as hox1, hox5, and hox8 were highly expressed in the anterior tissues, while posterior Hox genes such as hox9/10, hox11/13/a, and hox11/13b were highly expressed in the posterior ones ( Figure 3A ). No difference between anterior and posterior tissues was detected in the expression levels of hox7. Temporal expression peaks varied among the Hox genes. hox5 showed the highest expression level at 0 dpe in anterior tissue, while hox 8 showed it at 4 dpe in anterior tissue. For hox1, hox9/10, hox11/13a, and hox11/13b, the temporal expression peak occurred at 14–26 dpe.

Like Hox genes, Parahox genes, i.e., gsx, xlox, and cdx, also showed colinear expression patterns: xlox was expressed in anterior tissues, while cdx was expressed exclusively in posterior tissues ( Figure 3B ). These genes were expressed at a later stage of regeneration, after 14 dpe. Among the three examined Parahox genes, gsx showed a distinctive pattern, in which no difference was seen between anterior and posterior tissues, and the expression level was highest at 7 to 14 dpe.

Histological observations clearly revealed the process of nerve ring regeneration ( Figures 4B–E ). At the time when anterior evisceration occurred (0 dpe), the anterior end of the sea cucumber body shrank to close the wound, so the radial nerve tissues came close to each other. At the early stage (4 dpe), neural tissues showing histological characteristics of cell proliferation extended from wounded radial nerves in the body wall, entering the regenerating tissues at the center of the body trunk ( Figure 4B ). At 6 and 12 dpe, the apical ends of regenerating nerves were elongated in the perpendicular direction, forming ring structures by connecting to each other ( Figures 4C, D ). At 12 dpe, an almost complete ring-like structure was formed by the elongated neural tissues connecting to each other ( Figure 4D ). At 21 dpe, the formation of the nerve ring was completed ( Figure 4E ). In addition, at 12 dpe, tentacular canals of the water vascular system were apparent ( Figure 4D ) and at 21 dpe, regenerating tentacles were observed ( Figure 4E ).

Orthologous genes of 6 neural patterning genes, i.e., otx, six3/6, pax6, elav, musashi, and neuroD, were obtained from the gene database in E. quinquesemita ( Supplementary Figures 4 – 7 ). Based on the sequences of the gene orthologs, primers for qPCR were designed ( Supplementary Table 1 ). F or otx, six3/6, and pax6, the expression levels were higher in the anterior tissues, and almost no expression was detected in the posterior tissues ( Figure 5 ). For otx and six3/6, the expression peaks were seen at 4 dpe, and then the expression levels were decreased. For pax6, the expression peak was seen at 26 dpe. A neural marker gene, elav, showed a strong expression level immediately after the evisceration (0 dpe) in the anterior tissue ( Figures 3 , 4 ). The expression level of musashi, which is involved in neural differentiation was slightly higher in the anterior tissues, but it was relatively constant over the entire regeneration period. For neuroD, no difference in expression levels was detected between anterior and posterior tissues, or among regeneration stages ( Figure 5 ).

In general, twist and snail suppress cadherin expression (Foroni et al., 2012) and the downregulation of cadherin genes leads to the maintenance of mesenchyme ( Figure 2A , Kang and Massagué, 2004; Barrallo-Gimeno and Nieto, 2005; Lamouille et al., 2014; Pei et al., 2019). Thus, our results suggest that the expression level of cadherin could be low in anterior regenerating tissues, that is, the number of mesenchyme cells should be increased, along with the upregulation of twist and snail. This process is consistent with the histological observations showing that mesenchymal cells proliferated at the earlier stages of digestive-tract regeneration in this species (Okada and Kondo, 2019). In addition, the timing of twist expression was earlier than that of snail expression ( Figure 2B ), which is congruent with a report showing that snail expression is regulated downstream of twist (Smit et al., 2009). Therefore, this time difference between the expressions of these two factors is suggested to reflect the epistatic relationship between the two factors.

Recent RNA-seq analysis in the related species Eupentacta Fraudatrix, which also shows anterior evisceration, revealed the upregulation of EMT-related factors such as snai2 and id2 during regeneration after evisceration (Boyko et al., 2020). It was also suggested in a posterior eviscerating species that EMT was involved in the intestinal regeneration (García-Arrarás et al., 2011). Therefore, it is possible that induction of EMT at the onset of intestinal regeneration via regulators such as twist and snail is widely shared among holothurians.

On the other hand, histological observations suggest that MET is involved in the digestive tract regeneration of anterior tissues of E. quinquesemita, in which epithelial tissues appeared among mesenchymal tissues (Okada and Kondo, 2019). Although the qPCR results in this study suggest that EMT occurs in anterior regenerating tissues as a whole, down regulation of twist and snail should occur locally at the tissues where the epithelial formation occurs after 7 dpe. Histological investigation suggested that MET does not occur in the anterior regenerating tissues of E. fraudatrix (Mashanov and Dolmatov, 2001, Mashanov et al., 2005), so the dynamics of gene expression related to MET might be different between species.

The results of this study revealed a tendency for anterior Hox genes to be expressed anteriorly while posterior ones were expressed posteriorly in regenerating tissues ( Figure 3A ). It is known that, in sea urchin larvae, xlox is expressed at the border between stomach and gut (Arnone et al., 2006). This is consistent with our result in which xlox was expressed in the anterior region, since the border in E. quinquesemita is located at a relatively anterior portion of the digestive tract (Byrne, 2001; Okada and Kondo, 2019). As shown by histological observations, the developmental speed in posterior regenerating tissues is a bit slower than that in anterior tissues (Okada and Kondo, 2019). Therefore, the expression of cdx, which is required for formation of the posterior digestive tract, was also seen at later stages of the regeneration process ( Figure 3B ). The expression of another Parahox gene, gsx, is known to be localized at a part of the neural tube in the development of amphioxus (Brooke et al., 1998; Arnone et al., 2006). It is generally believed that, in ancestral animals, gsx is expressed at the anterior (or oral) side of the digestive tract (Samadi and Steiner, 2010), although there are many animal species that do not exhibit gsx expression in their digestive tracts (Wollesen et al., 2015). In this study, gsx expression was shown to peak during the middle period of regeneration ( Figure 3B ), suggesting that this gene shows diverse expression patterns depending on the animal species.

Overall, our results suggest that Hox and Parahox genes provide spatial information along the A-P axis for digestive-tract differentiation. Time differences of the expression peaks among Hox/Parahox genes indicate that the timepoints for determination and differentiation vary among digestive-tract parts. Upregulation of Hox genes in the regenerating intestine has been reported in several holothurians (e.g., Méndez et al., 2000), warranting further investigation of the conservation and diversification of detailed gene expression patterns and histology among species, especially anterior eviscerating species and posterior eviscerating species.

Our histological observations focusing on the nervous tissue regeneration showed that the nerves were regenerated from the anterior ends of the remaining radial nerve cords, and later connected to each other to form the ring structure ( Figure 4F ). A similar regeneration pattern is known in the closely related species E. fraudatrix (Dolmatov, 1992). Since seasonal evisceration is known in E. quinquesemita (Byrne, 1985), the complete renewal of the nerve ring observed in these species might play a key role in their lifecycle.

Gene expression analyses by real time qPCR showed that the neural patterning genes otx, six3/6 and pax6 were upregulated in the anterior body part during the regeneration process. In particular, otx and six3/6 were upregulated at a relatively early stage of regeneration (4 dpe, Figure 5 ). The timing of upregulation corresponds to the timing when the regenerating nerve tissues start to enter the regenerating tissues composed of mesenchymal cells that later form the digestive tract. This may suggest that otx and six3/6 are required for the differentiation and elongation of nerve tissues around the digestive tract. It was previously shown that, in the echinoderm species Holopneustes purpurescens (Echinoidea), the otx gene is expressed in the nerve ring during larval stages (Morris et al., 2004; Morris and Byrne, 2005). In a feather star, Anneissia japonica, otx and six3 are expressed in oral nervous tissues (Omori et al., 2020). Therefore, it is suggested that neural patterning genes such as otx and six3/6 are co-opted from normal development to the re-formation of the nerve ring after anterior evisceration.

On the other hand, pax6 was upregulated later in the regeneration process (14–26 dpe). Around at 12 dpe, histological observations showed that the morphological characters of tentacles became apparent ( Figures 4D, E ). As it is known that pax6 is expressed in tentacles in adults of the sea cucumber Apostichopus japonicus (Liu et al., 2020), the upregulation of pax6 associated with the tentacular formation is consistent with previous knowledge. In other holothurians, RNA-seq analyses suggest that conserved paracrine factors involved in embryonic development such as wnt6, wnt9, and bmp1 are upregulated during intestinal regeneration (Ortiz-Pineda et al., 2009; Sun et al., 2013; Yuan et al., 2019; Auger et al., 2023) and notch is upregulated in nerve cord regeneration (Mashanov et al., 2014). It is known that, in regeneration processes in vertebrates, expressions of transcription factors such as pax, six and hox genes are regulated downstream of the signals such as the bmp and wnt pathways (e.g., Grogg et al., 2005; Patel et al., 2022), so similar epistatic relationships are also suggested in the regeneration in echinoderms.