The genome of an apodid holothuroid (Chiridota heheva) provides insights into its adaptation to a deep-sea reducing environment

Cold seeps and hydrothermal vents are deep-sea reducing environments that are characterized by lacking oxygen and photosynthesis-derived nutrients. Most animals acquire nutrition in cold seeps or hydrothermal vents by maintaining epi- or endosymbiotic relationship with chemoautotrophic microorganisms. Although several seep- and vent-dwelling animals hosting symbiotic microbes have been well-studied, the genomic basis of adaptation to deep-sea reducing environment in nonsymbiotic animals is still lacking. Here, we report a high-quality genome of Chiridota heheva Pawson & Vance, 2004, which thrives by extracting organic components from sediment detritus and suspended material, as a reference for nonsymbiotic animal’s adaptation to deep-sea reducing environments. The expansion of the aerolysin-like protein family in C. heheva compared with other echinoderms might be involved in the disintegration of microbes during digestion. Moreover, several hypoxia-related genes (Pyruvate Kinase M2, PKM2; Phospholysine Phosphohistidine Inorganic Pyrophosphate Phosphatase, LHPP; Poly(A)-specific Ribonuclease Subunit PAN2, PAN2; and Ribosomal RNA Processing 9, RRP9) were subject to positive selection in the genome of C. heheva, which contributes to their adaptation to hypoxic environments.

E chinodermata is a phylum of marine animals comprising 5 extant classes, including Holothuroidea (feather star, subphylum Pelmatozoa), Asteriodea and Ophiuroidea (starfish and brittle star, subphylum Asterozoa), and Echinoidea and Holothuroidea (sea urchin and sea cucumber, subphylum Echinozoa) 1 . Adult echinoderms are characterized by having a body showing pentameral symmetry, a water vascular system with external tube feet (podia), and an endoskeleton consisting of calcareous ossicles 2 . Echinoderms exhibit a high divergence in morphology, from the star-like architecture in Asteroidea to the worm-like architecture in Holothuroidea 3,4 .
Compared with other echinoderms, holothurians have a unique body architecture and evolutionary history. The worm-like body of the holothurian preserves the pentameral symmetry structurally along the oral-aboral axis 5 . In addition, holothurians have a soft and stretchable body, in which the ossicles are greatly reduced in size 2 . The order Apodida is a group of holothurians that are found in both shallow-water and deep-sea environments 6 . Phylogenetic analyses showed that Apodida is sister to other orders of Holothuroidea 7,8 . Apodid holothurians lack tube feet and complex respiratory trees, making them morphologically distinct from other holothurians 2 . In contrast to other classes of Echinodermata, which experienced a severe evolutionary bottleneck during the Permian-Triassic mass extinction interval, Holothuroidea did not experience family-level extinction through the interval. The deposit-feeding lifestyle of holothurians conferred a selective advantage during the primary productivity collapse of the Permian-Triassic massextinction 9 . As the genomes of only two shallow-water holothurians (Apostichopus japonicus and Parastichopus parvimensis) have been assembled and analyzed [10][11][12] , it is critical to study the genomes of more holothurians to dissect their special morphological characteristics and evolutionary history.
Cold seeps are areas where methane, hydrogen sulfide, and other hydrocarbons seep or emanate as gas from deep geologic sources 13 . Hydrocarbon-fluid seepage from cold seeps is completely devoid of O 2 and comprises high levels of sulfides. After reacting with sulfides contained in the fluid, any free O 2 is removed from the deep-sea water. Thus, cold seeps are characterized by high hydrostatic pressure, low temperature, lack of oxygen, and photosynthesis-derived nutrients, and high concentrations of reducing chemicals 14 . Chemosynthetic microbes oxidize the reduced chemicals contained in the hydrocarbon fluids to produce energy and fix carbon into organic matter, which supports a large amount of invertebrates, including tubeworms, mussels, clams, and gastropods 15 . Most of these macrobenthos depend on the epi-or endosymbiotic relationships with chemoautotrophic microorganisms for nutrition 14,16,17 . Recent genomic analyses have revealed the genetic basis of adaptation in several seep-and vent-dwelling macrobenthos hosting symbiotic bacteria [18][19][20][21] . However, the genomic basis of nutrient acquisition and hypoxic adaptation in cold seep-dwelling nonsymbiotic animals is still lacking with only one reported genome 22 .
Echinoderms are a rare component of deep-sea chemosynthetic ecosystems 23 . Chiridota heheva Pawson & Vance, 2004 (Apodida: Chiridotidae) is one of the few echinoderms that occupies all three types of chemosynthetic ecosystems (hydrothermal vent, cold seep, and organic fall) 24 . This suggests that the species is well adapted to deep-sea reducing environments. Unlike most cold seep-and hydrothermal vent-dwelling species, C. heheva does not host chemosynthetic bacteria 6 . It derives nutrients from a variety of sources, extracting organic components from sediment detritus, suspended material, and wood fragments when available 6,25 . The cosmopolitan distribution and special lifestyle of C. heheva make it an ideal model to study adaptation to deep-sea reducing environments in nonsymbiotic animals.
Here, we assembled and annotated a high-quality genome of C. heheva collected from the Haima cold seep in the South China Sea. Evolutionary analysis revealed that the ancestor of C. heheva diverged from the ancestors of two shallow-water holothurians (A. japonicus and P. parvimensis) approximately 429 Ma ago. Additionally, demographic analysis suggested that C. heheva might have colonized the current habitat in the early Miocene. Comparative genomic analysis showed that the aerolysin-like protein (ALP) family was significantly expanded in C. heheva compared with other echinoderms. The expansion of the ALP family might be involved in the disintegration of microbes during digestion, which in turn facilitated its adaptation to cold seep environments. Moreover, several hypoxia-related genes were subject to positive selection in the genome of C. heheva, which contributes to their adaptation to hypoxic environments.

Results and discussion
Characterization and genome assembly of C. heheva. The sequenced sample was collected at a depth of 1385 meters using manned submersible Shenhaiyongshi from the Haima cold seep in the South China Sea (16°73.228′ N, 110°46.143′ E) (Fig. 1 Table 3). The C. heheva genome was assembled into 4609 contigs, with a total size of 1.107 Gb and contig N50 of 1.22 Mb (Table 1). We determined the completeness of the assembled genome by running benchmarking universal single-copy orthologs (BUSCO) and sequencing quality assessment tool (SQUAT) software. BUSCO analysis with metazoan (obd10) gene set showed that the assembled C. heheva genome contained 89.6% complete single-copy orthologs (Supplementary Table 4). Additionally, 91.1% of Illumina reads could be aligned to the assembled genome with high confidence in SQUAT assessment (Supplementary Table 5). These results indicate the high integrity of our assembled genome.
Genome annotation. Repetitive elements represented 624. 38 Mb in the C. heheva genome assembly (Supplementary Table 6). Long interspersed nuclear elements (LINEs) were the largest class of annotated transposable elements (TEs), making up 9.72% of the genome. DNA transposons, which were the second largest class of TEs, represented 33.59 Mb (3.03%) of the genome. Additionally, the C. heheva genome comprised a large proportion (38.39%) of unclassified interspersed repeats. Comparative genomic analysis among C. heheva and other echinoderms revealed that the C. heheva genome consisted of the largest number of TEs (Fig. 2a, b; Supplementary Table 7). Repetitive elements constituted 56.40% of the C. heheva genome, and they accounted for 26.68% and 25.02% of the genomes of A. japonicus and P. parvimensis, respectively. The differences in the repeat content were close to the size differences between the genomes of C. heheva and the other two holothurians. This suggests that repeats contributed to the size differences among the genomes of these three holothurians. Notably, the proportion of LINEs in the C. heheva genome was substantially higher than that in the genomes of other echinoderms (Fig. 2b). Kimura distance-based copy-divergence analysis identified a recent expansion of LINEs in the C. heheva genome (Fig. 2c). Protein-coding genes were identified in the genome of C. heheva through a combination of ab initio and homology-based protein-prediction approaches. In total, we derived 36,527 gene models in the C. heheva genome. The structure of predicted genes in C. heheva is slightly different to that of other previously sequenced echinoderm genomes. With longer exons and introns, genes in C. heheva are longer than the ones in A. japonicus (Supplementary Table 8).
Phylogenomic analysis and demographic inference. With more than 1400 extant species, Holothuroidea is one of the largest classes in the phylum Echinodermata 1 . In addition, holothurians are well adapted to diverse marine environments, with habitats ranging from shallow intertidal areas to hadal trenches 27,28 . However, due to the lack of body fossils, evolutionary study of Holothuroidea is more difficult than other classes of Echinodermata. To investigate the evolutionary history of C. heheva, a maximum-likelihood (ML) phylogenetic tree was reconstructed using single-copy orthologs of C. heheva and 16 other deuterostomes ( Supplementary Fig. 2). Chiridota heheva appeared sister to two other holothurians. In addition, divergence times were determined among 7 echinoderms that had whole-genome sequences (Fig. 3a). The divergence time of A. japonica and other echinoderms was estimated to be approximately 539 million years (Ma), which is generally consistent with the fossil records 29, 30 . The ancestor of Chiridota heheva diverged from the ancestors of two shallow-water holothurians (A. japonicus and P. parvimensis) approximately 429 Ma ago. As Apodida is the basal taxon in Holothuroidea, these results support the view that holothurians had evolved by the Early Ordovician 31 .
To better investigate the evolution of holothurians, we inferred the histories of ancestral-population sizes of C. heheva and A. japonicus using the pairwise sequential Markovian coalescent (PSMC) method (Fig. 3b). Chiridota heheva experienced a decline in population size approximately 21 Ma ago. Ocean temperature increased slowly between the late Oligocene and early Miocene (21-27 Ma ago) after long-term cooling from the end of the Eocene 32,33 . Furthermore, species diversity within Echinodermata started to increase in the early Miocene 34,35 . These results indicate that C. heheva might have colonized the current habitat in the early Miocene. A decline in ancestral-population size in A. japonicus started in the late Miocene (approximately 8 Ma ago). Chiridota heheva also experienced a moderate decline in population size in the early Pliocene. Additionally, the oceans experienced a decrease in temperature during the late Miocene (7-5.4 Ma ago) 36 . These results suggest that global cooling and environmental changes in the late Miocene were an important driver of demographic changes in both shallow-water and deepsea holothurians.
Hox/ParaHox gene clusters. Apodida do not have tube feet or complex respiratory trees, which are commonly found in other holothurians 37 . It has been demonstrated that Hox genes play a critical role in embryonic development 38 . In addition, previous studies proposed that the presence/absence and expression pattern of Hox genes might contribute to morphological patterning and embryonic development in echinoderms 10,11 . Therefore, to determine whether Hox genes contribute to morphological divergence in Holothuroidea, we identified Hox gene clusters and their evolutionary sister complex, the ParaHox gene cluster, in the genomes of C. heheva and 6 other echinoderms ( Supplementary Fig. 3). A Hox cluster and a ParaHox cluster could be identified in the genomes of all 7 species. The gene composition and arrangement of both Hox and ParaHox clusters were highly consistent between the genomes of C. heheva and A. japonicus, suggesting that Hox/ ParaHox genes do not control the development of tube feet and respiratory trees in Apodida. Hox4 and Hox6 were missing in the genomes of both C. heheva and A. japonicus, which is inconsistent with the view that the loss of Hox4 or Hox6 in echinoderms is a lineage-specific event 5 .  Table 9). Taken together, these results indicate that the C. heheva NLR repertoire shows different abundances and structural complexities than the sea urchin. We performed phylogenetic analysis of C. heheva NLRs and other representative NLRs of metazoans, including humans, Amphimedon queenslandica, S. purpuratus, Acropora digitifera, Nematostella vectensis, Pinctada fucata, Capitella teleta, mollusks, and arthropods 43 . We found that the majority of C. heheva NLRs form a monophyletic lineage with sea-urchin NLRs (Fig. 4), supporting the lineage-specific evolution of NLRs in Echinodermata 44 . Given that human IPAF (ice proteaseactivating factor) and NAIP (neuronal apoptosis-inhibitory protein) proteins were reported to have originated before the evolution of vertebrates 44 , one C. heheva NLR clustering with these two proteins indicates that this NLR may have an ancient independent origin (Fig. 4).
Gene-family evolution. We performed gene-family analysis based on the phylogenetic tree of 7 echinoderms (Fig. 3a). Compared with other echinoderms, 66 gene families were expanded, and 25 gene families were contracted in C. heheva (P < 0.05) (Supplementary Data 1 and Supplementary Table 10). Several significantly expanded gene families are involved in the processes of cell cycle progression, protein folding, and ribosome assembly. As high hydrostatic pressure causes cell cycle delay and affects protein folding 45,46 , expansion of these families may have contributed to the adaptation of C. heheva to cold seep environments.
Aerolysins, which are pore-forming toxins (PFTs), were first characterized as virulence factors in the pathogenic bacterium Aeromonas hydrophyla 47,48 . As typical pore-forming proteins, aerolysin and related proteins are found in a large variety of species and possess diverse functions 49 . ALPs in eukaryotes originated from recurrent horizontal gene transfer (HGT) 50 . ALPs of the same origin have similar functions, while the ones of different origins possess diverse functions 50 . The ALPs were   significantly expanded in the genome of C. heheva (7 copies) compared with other echinoderms (0 or 1 copy) (P < 0.05) (Supplementary Data 1). To investigate the possible origin and function of C. heheva ALPs, we reconstructed the phylogenetic tree of ALPs in echinoderms and diverse species. Chiridota heheva ALPs did not cluster with ALPs from other echinoderms. Additionally, these two groups of ALPs were clustered with aerolysins from distinct groups of bacteria (Fig. 5). This suggests that ALPs from C. heheva and other echinoderms have different origins. Chiridota heheva ALPs form a clade with ALPs from stony corals (Stylophora pistillata, Pocillopora damicornis, and Orbicella faveolata) and sea anemones (Nematostella vectensis and Ecaiptasia diaphana). This indicates that ALPs from C. heheva, stony corals, and sea anemones might have the same origin and similar biological functions. It was shown that ALPs from hydra and sea anemones (N. vectensis) are involved in prey disintegration after predation by lysing cells through pore formation on membranes 50,51 . The microbial communities of cold seeps are very different from those of other seafloor ecosystems 52 . Moreover, some of these microbes have unique cellular structures that might be difficult to disintegrate 53 , which impedes nutrient acquisition of C. heheva from free-living microbes of cold seeps. Therefore, the expansion of the ALP family might have contributed to microbe digestion in C. heheva, which in turn facilitated its adaptation to cold seep environments.

Percentage (%) C . h e h e v a A . j a p o n i c u s P . p a r v i m e n s i s A . j a p o n i c a A . p l a n c i L . v a r i e g a t u s S . p u r p u r a t u s c
Positively selected genes. To better understand the genetic basis of its adaptation to a deep-sea reducing environment, we searched for positively selected genes (PSGs) in C. heheva. Compared with 6 other echinoderms, 27 PSGs were identified in the C. heheva genome (Supplementary Table 11). Four hypoxia-related genes (pyruvate kinase M2, PKM2; phospholysine phosphohistidine inorganic pyrophosphate phosphatase, LHPP; poly(A)-specific ribonuclease subunit PAN2, PAN2; and ribosomal RNA processing 9, RRP9) were identified as PSGs in C. heheva [54][55][56][57] . PKM2 promotes transactivation of HIF-1 target genes by directly interacting with the HIF-1α subunit. In addition, the transcription of the PKM2 gene is activated by HIF-1. This positive-feedback loop increases glycolysis and lactate production and decreases oxygen consumption under hypoxic conditions 54 . LHPP interacts with PKM2 to induce ubiquitin-mediated degradation of PKM2 and impede the glycolysis and respiration under hypoxia 55 . Thus, selection against these two interacting genes (PKM2 and LHPP) might play a key role in the hypoxic adaptation in C. heheva. Interestingly, the LHPP was also subject to positive selection in cetaceans, which are hypoxia-tolerant mammals 58 . Furthermore, both C. heheva and cetaceans have the same amino acid substitution at position 118 of the LHPP protein (Fig. 6), which indicates a possible convergent evolution in the LHPP during the adaptation of cetaceans and C. heheva to hypoxic environments. A positively charged amino acid (histidine, H) in two shallow-water  Fig. 3 Evolutionary history of C. heheva. a A species tree of 7 echinoderm species. In total, 988 single-copy orthologs were used to reconstruct the phylogenetic tree. The divergence time between species pairs was listed above each node, and 95% confidence interval of the estimated divergence time was denoted as blue bar. The numbers of protein families that were significantly expanded (red) and contracted (blue) (P < 0.05) in each species are denoted beside the species names. b Demographic history of C. heheva (blue) and A. japonicus (red). The changes of ancestral-population size of C. heheva and A. japonicus were inferred using the PSMC method. Time in history was estimated by assuming a generation time of 3 years and a mutation rate of 1.0 × 10 −8 . holothurians is substituted to a negatively charged one (aspartic acid, D) in C. heheva at this position, which might cause a conformation change that contributes to the hypoxic adaptation in C. heheva. To study the potential underlying structural effects of this substitution, we predicted the three-dimensional structures of LHPP from echinoderms ( Supplementary Fig. 4). The substitution, which is located in an α-helix, does not change the conformation of LHPP. The effect of the substitution of the LHPP protein needs to be further investigated. A large number of metazoans reside in cold seeps and hydrothermal vents, which are challenging environments with high concentration of toxic compounds and chronic hypoxia 59,60 . Several physiological and molecular modifications of the respiratory system have been identified to cope with hypoxia in organisms living in these environments 60 . The concentration of hemoglobin, which is an oxygen-binding protein, is higher in seep-and vent-dwelling species than closely related species living in well-oxygenated environments 60 . In addition, hemoglobins from deep-sea organisms have higher affinity for oxygen than hemoglobins from shallowwater relatives 61 . This facilitates seep-and vent-dwelling species to thrive in the extreme environments by improving the efficiency of oxygen transportation. Four hypoxia-related genes (PKM2, PAN2, LHPP, and RRP9) were identified to be positively selected in C. heheva. PKM2 increases glycolysis and decreases oxygen consumption by promoting transactivation of HIF-1 target genes through directly interacting with the HIF-1α subunit under hypoxic conditions 54 . This suggests that animals living in deep-sea chemosynthetic environments might also adapt to hypoxic conditions through reprogramming glucose metabolism. Intriguingly, LHPP gene was subjected to positive selection in both C. heheva and cetaceans. This indicates a possible convergent evolution, in which echinoderms and mammals utilize similar strategies to cope with hypoxic challenges.     A q u 1 _ 2 0 9 6 9 6 ad i_v 1_ 15 68 9 Ch iHe h_g ene 119 99 S P U _ 0 0 3 5 3 9

Methods
A q u 1 _ 2 0 8 1 8 9 Ch iH eh _g en e7 54 5 Ch iH eh _g en e1 59 03   To construct Nanopore sequencing library, high-molecular-weight genomic DNA was prepared by the CTAB method. The quality and quantity of the DNA were measured via standard agarose-gel electrophoresis and with a Qubit 4.0 Fluorometer (Invitrogen). Sequencing library was constructed and sequenced by Nanopore PromethION platform (Oxford Nanopore Technologies). Additionally, DNA was extracted to construct Illumina sequencing library. The quality and quantity of the DNA were measured via standard agarose-gel electrophoresis and with a Qubit 2.0 Fluorometer (Invitrogen). Sequencing library was constructed and sequenced by Illumina Novaseq platform (Illumina).

Capitella teleta Homo sapiens P in c ta d a fu c a ta
Mitochondrial and nuclear genome assembly. Low-quality (reads with ≥10% unidentified nucleotide and/or ≥ 50% nucleotides having phred score < 5) and sequencing-adapter-contaminated Illumina reads were filtered and trimmed with Fastp (v0.21.0) 62 to obtain high-quality Illumina reads, which were used in the following analyses. Mitochondrial genome of C. heheva was assembled using the two-step mode of mitoZ (v2.4) 63 with the high-quality Illumina reads. The assembled genome was annotated using mitoZ (v2.4) with parameter "-clade Echinodermata".
The size and heterozygosity of C. heheva genome were estimated using highquality Illumina reads by k-mer frequency-distribution method. The number of kmers and the peak depth of k-mer sizes at 17 was obtained using Jellyfish (v2.3.0) 64 with the -C setting. Genome size was estimated based on the k-mer analysis as described previously 65 . The heterozygosity of C. heheva genome was determined by fitting the k-mer distribution of Arabidopsis thaliana using Kmerfreq implemented in SOAPdenovo2 (r242) 66 .
Low-quality Nanopore reads were filtered using custom Python script. Two draft-genome assemblies were generated using filtered Nanopore reads with Shasta (v0.4.0) 67 and WTDBG2 (v2.5) 68 , respectively. The contigs of the two draft assemblies were subject to error correction using filtered Nanopore reads with Racon (v1.4.16) 69 three times. The corrected contigs were then polished with high-quality Illumina reads with Pilon (v1.23) 70 three times. The error-corrected contigs of Shasta assembly and WTDBG2 assembly were assembled into longer sequences using quickmerge (v0.3) 71 . The merged contigs were subject to error correction using filtered Nanopore reads with Racon three times, and then using high-quality Illumina reads with Pilon three times. As the heterozygosity of C. heheva genome is high, haplotypic duplications in the assembled genome were identified and removed using purge_dups (v1.2.3) 72 . The completeness and quality of the assembly was evaluated using BUSCO (v4.0.5) 73 against the conserved Metazoa dataset (obd10), and SQUAT with high-quality Illumina reads 74 .
Genome annotation. Repetitive elements in the assembly were identified by de novo predictions using RepeatMasker (v4.1.0) (https://www.repeatmasker.org/). A de novo repeat library for C. heheva was built using RepeatModeler (v2.0.1) 75 . To identify repetitive elements, sequences from the C. heheva assembly were aligned to the de novo repeat library using RepeatMasker (v4.1.0). Additionally, repetitive elements in C. heheva genome assembly were identified by homology searches against known repeat databases using RepeatMasker (v4.1.0). A repeat landscape of C. heheva genome was obtained using an R script that was modified from https:// github.com/ValentinaBoP/TransposableElements. To compare the proportion and composition of repetitive elements among the genomes of echinoderms, genome sequences of Strongylocentrotus purpuratus (GCA_000002235.3), Lytechinus variegatus (GCA_000239495.2), Acanthaster planci (GCA_001949165.1), and Anneissia japonica (GCA_011630105.1) were downloaded from NCBI. Genome sequence of Parastichopus parvimensis was downloaded from echinobase (http:// bouzouki.bio.cs.cmu.edu/Echinobase/PpDownloads, retrieved September 2021). Repetitive elements in the genomes of these species were identified by homology searches against known repeat databases using RepeatMasker (v4.1.0). The proportion and composition of repetitive elements of Apostichopus japonicus was obtained from Li et al. (2018) 11 .  51 10 0   We applied a combination of homology-based and de novo predication methods to build consensus-gene models for the C. heheva genome assembly. For homology-based gene prediction, protein sequences of Helobdella robusta, Lytechinus variegatus, Strongylocentrotus purpuratus, Dimorphilus gyrociliatus, Apostichopus japonicus and Acanthaster planci were aligned to the C. heheva genome assembly using tblastn. The exon-intron structures then were determined according to the alignment results using GenomeThreader (v1.7.0) 76 . In addition, de novo gene prediction was performed using Augustus (v3.3.2) 77 , with the parameters obtained by training the software with protein sequences of Drosophila melanogaster and Parasteatoda tepidariorum. Two sets of gene models were integrated into a nonredundant consensus-gene set using EvidenceModeler (v1.1.1) 78 . To identify functions of the predicted proteins, we aligned the C. heheva protein models against NCBI NR, trEMBL, and SwissProt database using blastp (Evalue threshold: 10 −5 ), and against eggNOG database 79 using eggNOG-Mapper 80 . In addition, KEGG annotation of the protein models was performed using GhostKOALA 81 .
To estimate the divergent time among echinoderms, single-copy orthologs were identified among A. japonica, A. planci, A. japonicus, P. parvimensis, C. heheva, L. variegatus, and S. purpuratus after running OrthoMCL pipeline as mentioned above. Single-copy orthologs were aligned using MAFFT (v7.310), trimmed using ClipKit (v1.1.3) with 'gappy' mode, and concatenated using PhyloSuite (v1.2.2) 88 . Divergent time among 7 echinoderms were estimated using the concatenated alignment with MCMCtree module of the PAML package (v4.9) 89 . MCMCtree analysis was performed using the maximum-likelihood tree that was reconstructed by IQ-TREE2 as a guide tree and calibrated with the divergent time obtained from TimeTree database (minimum = 193 million years and soft maximum = 350 million years between L. variegatus and S. purpuratus) 90 .
Demographic inference of C. heheva and A. japonicus. Paired-end Illumina reads of A. japonicus were downloaded from NCBI SRA database 11 . The reads of A. japonicus were filtered and trimmed with fastp (v0.21.0). The Illumina clean reads of C. heheva and A. japonicus were aligned to the respective reference-genome assembly using BWA (v0.7.17) 91 with "mem" function. Genetic variants were identified using Samtools (v1.9) 92 . Whole-genome consensus sequence was generated with the genetic variants using Samtools (v 1.9). PSMC (v0.6.5) 93 was used to infer the demographic history of C. heheva and A. japonicus using the whole genome consensus sequences. The substitution mutations rate and generation time of C. heheva and A. japonicus was set to 1.0 × 10 −8 and 3 years according to the previous study of A. planci 94 .
Homeobox gene analysis. Homeobox genes in C. heheva genome were identified by following the procedure described previously 95 . Homeodomain sequences, which were retrieved from HomeoDB database (http://homeodb.zoo.ox.ac.uk) 96 , were aligned to C. heheva genome assembly using tbalstn. Sequences of the candidate homeobox genes were extracted based on the alignment results. The extracted sequences were aligned against NCBI NR and HomeoDB database to classify the homeobox genes.
Identification of NOD-like receptors (NLRs) in C. heheva. We used HMMER (v3.1) 97 to search against the proteome of C. heheva with the HMM profile of NACHT domain (PF05729) retrieved from Pfam 34.0 as the query and an e cutoff value of 0.01. Proteins identified by the HMM search were retrieved from the proteome and aligned with 964 representative proteins from eukaryotes and prokaryotes 98 , and other representative metazoan NLRs 43 using hmmalign method implemented in HMMER (v3.1) based on the STAND NTPase domain. The alignment was refined by manual editing. The large-scale phylogenetic analysis was performed using an approximate maximum likelihood method implemented in FastTree 99 . Representative SWACOS and MalT NTPases were used as outgroups 98 . Significant hit clustering with metazoan NLRs was regarded as NLRs, and proteindomain organizations were annotated through hmmscan method implemented in HMMER (v3.1).
Phylogenetic analysis of Chiridota NLRs and other representative metazoan NLRs. To explore the evolutionary relationships among C. heheva NLRs and other representative metazoan NLRs, we reconstructed the phylogenetic tree of NLRs. The NACHT domains of C. heheva NLRs and representative metazoan NLRs were aligned using MAFFT (v7.310), and then refined by manual editing. The representative metazoan NLRs were chosen from literature 43 . The phylogenetic tree was reconstructed using a maximum-likelihood method implemented in IQ-TREE2 (v2.1.2). The best-fit substitution model was selected by using ModelFinder   Fig. 6 A possible amino acid substitution of LHPP that contributed to hypoxic adaptation in C. heheva and cetaceans. The maximum-likelihood phylogenetic tree of cetaceans, C. heheva, and other echinoderms was reconstructed using 598 single-copy orthologs. C. heheva and cetaceans, which are tolerant to hypoxia, have the same amino acid substitution at position 118 of LHPP protein.
algorithm. Branch supports were assessed using the UFBoot approach with 1000 replicates.
Gene-family expansion and contraction analysis. r8s (v1.7) 100 was applied to obtain the ultrametric tree of 7 echinoderm species, which is calibrated with the divergent time between A. planci and S. purpuratus (541 mya) obtained from TimeTree database. CAFÉ(v5) 101 was applied to determine the significance of gene-family expansion and contraction among 7 echinoderm species based on the ultrametric tree and the gene clusters determined by OrthoMCL (v2.0.9). The divergence time reported by TimeTree database might not be precise as it is a consensus of divergence times estimated in previous studies. Therefore, we repeated the analysis twice, in which the divergence time between A. planci and S. purpuratus was set to 461 mya and 495 mya according to the previous studies, respectively 102,103 . All the three analyses had the same result.
We used HMMER (v3.3.2) to search against NCBI nonredundant protein database (accessed on July 2021) with the HMM profile of aerolysin domain (PF01117) retrieved from Pfam 34.0 as the query and a e cutoff value of 0.001. Proteins identified by the HMM search were retrieved and filtered for the ones that have less than 75 residues. The filtered proteins were aligned with aerolysin-like proteins (ALPs) from C. heheva, A. japonicas, and P. parvimensis using MAFFT (v7.310) and trimmed using ClipKit (v1.1.3) with 'gappy' mode. The phylogenetic tree was reconstructed using a maximum-likelihood method implemented in IQ-TREE2 (v2.1.2). The best-fit substitution model was selected by using ModelFinder algorithm. Branch supports were assessed using the UFBoot approach with 1000 replicates.
Identification and analysis of positively selected genes. Branch-site models implemented in the codeml module of the PAML package is widely used to identify positively selected genes (PSGs). Thus, we identified PSGs in the C. heheva genome within the single-copy orthologs among 7 echinoderm species, based on the branch-site models using GWideCodeML (v1.1) 104 . C. heheva was set as the 'foreground' phylogeny, and the other species were set as the 'background' phylogeny. An alternative branch site model (Model = 2, NSsites = 2, and fix_omega = 0) and a neutral branch site model (Model = 2, NSsites = 2, fix_omega = 1, and omega = 1) were tested. Genes with Bayesian empirical Bayes (BEB) sites > 90% and a corrected P-value < 0.1 were identified to have been subject to positive selection.
Statistics and reproducibility. Alpha levels of 0.05 were regarded as statistically significant throughout the study, unless otherwise specified.

Data availability
Raw reads and genome assembly are accessible in NCBI under BioProject number PRJNA752986. Assembled genome sequences are accessible under Whole Genome Shotgun project number JAIGNY000000000. Raw reads and genome assembly are also available at the CNGB Sequence Archive (CNSA) of China National GeneBank DataBase (CNGBdb) with accession number CNP0002134. The genome assembly, related annotation files, and source files for generating figures can be accessed through

Code availability
Custom script used in this study is available at Figshare 106 (https://doi.org/10.6084/ m9.figshare.15302229). Versions and parameters for other software packages used in this study are described in the reporting summary and elsewhere in the "Methods." Received: 2 October 2021; Accepted: 16 February 2022;