With their incredible diversity in morphology, life history and ecology, molluscs are of great interest to diverse fields of biology ranging from ecology to evolutionary developmental biology to biomechanics. To date there are high-quality genomes available for six of the eight molluscan classes publicly available, although sampling is heavily biased toward the economically and ecologically important members of the clade Conchifera (e.g. Gastropoda, Bivalvia and Cephalopoda; reviewed by Gomes-dos-Santos et al. 2020 and Sigwart et al. 2021). As the sister taxon to all other molluscs, the clade Aculifera (Polyplacophora + Aplacophora) is important to our understanding of molluscan evolution and innovation. Despite their importance, only three aculiferan genomes have been published to date, two species of chitons (Polyplacophora; Varney et al. 2021; Varney and Yap-Chiongco et al. 2022) and one species of Caudofoveata ( Wang et al. 2024). Therefore, we sequenced the genomes of two species of Solenogastres – Neomenia megatrapezata Salvini-Plawen & Paar-Gausch 2004 and Epimenia babai Salvini-Plawen, 1997 – to facilitate further study of this evolutionarily important taxon.

The two clades of aplacophoran molluscs are united by their distinctive vermiform shape, loss or reduction of the foot, and a chitinous cuticle covered by calcareous integumental spines and scales (sclerites). Solenogastres are distinctive with the presence of a ciliated foot restricted to a ventral groove, an undifferentiated gut, and the lack of true ctenidia, while members of Caudofoveata lack a foot entirely and have the presence of a specialized oral shield used for burrowing in mud ( Todt 2013). The current taxonomy of Solenogastres divides the group into 24 families and four orders (Pholidoskepia, Cavibelonia, Sterrofustia, and Neomeniamorpha) although recent phylogenomic analyses have indicated the need for significant taxonomic revision within this group ( Kocot et al. 2019; Yap-Chiongco et al. 2024).

Whereas most solenogasters are just a few millimeters in length ( Todt 2013) both N. megatrapezata and E. babai are large enough to obtain high-quality genomic and transcriptomic data from a single individual. E. babai is a species of Cavibelonia that is relatively easily collected at SCUBA-accessible depths where it feeds on the soft coral Scleronephthya gracillima ( Okusu 2002). Because of its relatively large-body size (up to 20 cm in length) and reliable spawning within the laboratory, E. babai remains one of very few species of Solenogastres in which development has been described ( Baba 1938, 1940, 1951; Okusu, 2002; Todt and Wanninger 2010). Sequencing of its genome coupled with characterized development and reliability of spawning make E. babai a desirable species to serve as a developmental model within Solenogastres. N. megatrapezata (up to 18 cm in length) is a species of Neomeniamorpha found in the Southern Ocean off Antarctica ( Salvini-Plawen & Paar-Gausch 2004; Kocot et al. 2019). Neomenia is interesting to study trait evolution within Solenogastres as species within this genus have a complex accessory genital apparatus and lack both ventral foregut glands and a radula ( García-Álvarez & Salvini-Plawen 2007). Together, the sequencing of these two species expands our ability to explore solenogaster evolution and, more broadly, provides an important resource in identifying large-scale patterns in molluscan evolution.

For N. megatrapezata, Illumina genome sequencing yielded 443.02 M reads. Illumina transcriptome sequencing of the same individual of N. megatrapezata used for genome sequencing (Ap227) yielded 29.72 M reads and 32.93 M reads for a second individual of N. megatrapezata (Ap259). Transcriptomes of other Neomenia species resulted in 110.68 M reads for N. aff. herwigi and 24.40 M reads for N. permagna. GenomeScope analysis estimated the genome size of N. megatrapezata to be 294 Mbp and a heterozygosity of 0.142% based on the trimmed PE reads. Following adapter trimming, two flowcells of Oxford Nanopore sequencing yielded 101.31 M reads. Assembly with MaSuRCA resulted in an assembly of 9,870 contigs totaling 347 Mbp (N50 = 119 Kbp, L50 = 1,041) with a BUSCO completeness score of 87.1% (83.9% single, 3.2% duplicated, 8.0% fragmented, and 4.9% missing). After polishing and purging redundant haplotigs, the final N. megatrapezata assembly was reduced to 6,168 contigs totaling 412 Mbp (N50 = 132 Kbp, L50 = 881) with a BUSCO completeness score of 85.1% (81.7% single, 3.4% duplicated, 8.1% fragmented, and 6.8% missing). BRAKER predicted 65,328 gene models with 90.6% of BUSCOs detected (80.5% single, 10.1% duplicated, 4.5% fragmented, and 4.9% missing). Removal of gene models not supported by transcriptome or protein evidence resulted in a final gene set for N. megatrapezata of 22,463 with a BUSCO completeness score of 90.2% (81.0% single, 9.2% duplicated, 4.7% fragmented, and 5.1% missing).

For E. babai, Illumina PE sequencing yielded 693.24 M reads. Illumina transcriptome sequencing of E. babai and Epimenia sp. yielded 70.06 M and 131.47 M reads, respectively. GenomeScope analysis of PE reads failed to converge for E. babai with k-mer sizes of 17 and 21. Therefore, we used the trimmed nanopore data in the program KMC 3 ( Kokot et al. 2017) with the options -k21 -fm -t16 -m500 -ci1 -cs10000 to obtain k-mer statistics followed by kmc-tools with the transform option to create a histogram for input into GenomeScope2. Oxford Nanopore sequencing of three flowcells data yielded 14.79M reads post adapter trimming with PoreChop. The resulting analysis based on Nanopore data estimated a genome size for E. babai of 543 Mbp and heterozygosity of 1.16%. MaSuRCA yielded an assembly of 11,695 contigs totaling 572 Mbp (N50 = 271 Kbp, L50 = 618) with a BUSCO completeness score of 91.3% (77.7% single, 13.6% duplicated, 5.2% fragmented, and 3.5% missing). The final polished and purged Epimenia assembly was reduced to 5,965 contigs totaling 628 Mbp (N50 = 413 Kbp, L50 = 370) with a BUSCO completeness score of 90.1% (82.0% single, 8.1% duplicated, 6.0% fragmented, and 3.9% missing). BRAKER predicted 123,904 gene models with a BUSCO completeness score of 93.3% (80.9% single, 12.4% duplicated, 4.7% fragmented, and 2.0% missing). Removal of genes not supported by transcriptome or protein evidence resulted in a final gene set for E. babai of 25,393 with a BUSCO completeness score of 92.4% (80.5% single, 11.9% duplicated, 4.9% fragmented, and 2.7% missing).

Molluscan genomes sequenced to date range in size from 359 Mbp ( Simakov et al. 2013) to over 6 Gbp ( Song et al. 2023) and genome size estimates by Adachi et al. (2021) for 141 diverse molluscs (Polyplacophora, Gastropoda, Bivalvia, and Cephalopoda) based on flow cytometry ranged from 469.4 Mbp (measured as a C-value of 0.48) to 5.31 Gbp (C-value of 5.43) with a mean size of 2.61 Gbp (C-value of 2.67) among the sampled molluscan taxa. Kocot et al. (2016) utilized Feulgen image analysis densitometry (FAID) to estimate genome sizes in select species of aplacophorans, including the same individual of Epimenia babai sequenced here. Genome size of E. babai was estimated to be 750 Mbp and that of Neomenia permagna was estimated to be 293 Mbp ( Kocot et al. 2016). Genome size for E. babai inferred by GenomeScope based on Oxford Nanopore data is smaller than expected based on FIAD at 543 Mbp. However, we were unable to obtain an estimate based on more accurate PE data as the analysis failed to converge. The estimated genome size for N. megatrapezata of 294 Mbp inferred by GenomeScope is comparable to the FIAD estimate of 293 Mbp for the closely related species N. permagna. With final assembly sizes of 628 Mbp and 412 Mbp for E. babai and N. megatrapezata, respectively, solenogaster genomes are fairly small with respect to other molluscs. Based on computational and FAID genome size estimates, genome size within Aculifera is quite variable ranging from as small as 293 Mbp in N. permagna to 2.84 Gbp in Cryptochiton stelleri ( Kocot et al. 2016; Hinegardner, 1974). Comparing the aplacophoran classes, the estimated genome size of Chaetoderma sp. (Caudofoveata) is quite large (2.45 Gbp) in comparison to both E. babai and N. megatrapezata ( Table 1). Although in some groups it has been shown that genome size may co-evolve with life-history traits, the observed high level of variation in genome size within Aculifera raises questions about the underlying mechanisms driving this heterogeneity in genome size ( Ritchie et al. 2017; Beaudreau et al. 2021). Further comparisons between solenogaster, caudofoveate, and chiton genomes will aid in our understanding of the evolutionary history of Aplacophora and Aculifera and the potential interplay between genome size and, e.g., life history.

Comparison of the full set of Epimenia and Neomenia gene models to the gene models from 21 other lophotrochozoans in OrthoFinder resulted in 251,468 groups of homologous sequences. Our pipeline selected 7,561 single-copy genes sampled for at least 18 of the 23 taxa. The final PhyloPyPruner output resulted in a supermatrix of 3,577 genes with a total of 1,169,375 amino acids and 20.3% missing data. Of these, E. babai was sampled for 2,319 and N. megatrapezata was sampled for 2,232.The resulting tree ( Figure 1) is strongly supported overall with maximal support for Solenogastres and for Aculifera as the sister taxon to all other sampled molluscs, consistent with recent studies (e.g., Kocot et al. 2020; Song et al. 2023).

In summary, the genomes of N. megatrapezata and E. babai are relatively complete and comparable to other published aculiferan genomes ( Table 1) with BUSCO completeness scores of 85.1% and 90.1%, respectively. These represent the first sequenced genomes for Solenogastres and are a valuable resource for future studies of molluscan evolution.