Introduction

The genomic mechanisms that enable or limit the evolution of major morphological changes remains one of the great questions of evolutionary biology. Living members of Mollusca represent the broadest morphological disparity of any animal phylum: mollusc body plans encompass squid, worms, living rocks and candy-coloured tree snails, as well as diverse intermediate and additional novel forms known from the fossil record. Early work on chromosome numbers in meiotic division used molluscs, and interest in understanding patterns in chromosome numbers across species also led to fundamental insights in animal polyploidy. Despite playing a key role in early advances, this important phylum lags in terms of quality and taxonomic coverage of whole genome sequence data. Reconstructing the evolution of genome architecture through deep time in diverse molluscs remains critical to understand genome evolution in animals. High quality genomic data for the deeply divergent, morphologically constrained chitons, would be expected to offer an opportunity to explore ancient genetic traits and evolutionary mechanisms preserved across the long span of animal evolution.

Prior work has reasonably assumed that rates of intra-chromosomal gene translocation are constant within major groups (Mackintosh et al., 2023). If this is true, syntenic rearrangement could be a clocklike indicator of divergence times. But rates of genomic rearrangement are not well known in molluscs, nor how they might vary across this vast clade, and higher rates of rearrangement confound reconstruction of ancestral states. In molluscs, an excellent fossil record allows an independent record of divergence times that will provide more insight into the variability (and utility) of rates of rearrangement as a measure of divergence time.

Chitons have long been considered the key taxon to understanding the ancestral molluscan body plan (Wanninger and Wollesen, 2019). Around 1000 living chiton species worldwide all possess an eight-part shell armour that has remained superficially unchanged with relatively little variation since the Carboniferous, over 300 Million years ago (Sigwart, 2017). Chitons are increasingly important for bio-inspired design, which will benefit from genomic tools to understand genetic control of their flexible armour and unique sensory systems (Ampuero et al., 2024; Varney et al., 2021). Chitons are generally conserved, yet the known species richness in chitons is higher than for the more morphologically and behaviourally diverse extant cephalopods. The adaptations that drive speciation in this group remain an open question.

One key issue is whether polyplacophoran molluscs are conservative per se or whether their diagnostic body plan adaptations are so distinctive and strange that this has overshadowed full understanding of additional adaptive traits: fully innervated shells, iron mineralised radula, living at almost all depths and latitudes of the world ocean. While the body plan of chitons has persisted for over 300 My, this is a template for remarkable adaptations that have only recently begun to be appreciated.

Chitons are grazers and most are not highly ecologically specialised (Sigwart and Schwabe, 2017). These animals mostly have separate sexes and all are broadcast spawners, are not migratory, cannot reasonably be subject to strong sexual selection, and many similar species co-occur in parapatric or sympatric radiations (Kelly and Eernisse, 2008). Previously published karyotype data shows closely related species with overlapping ranges can have different numbers of chromosomes, such as Acanthochitona discrepans (1n=8) and A. crinita (1n=9) in the North Atlantic (Certain, 1951). Chromosome rearrangement could be an attractive speculative explanation that could support maintaining species boundaries in this clade.

Here, we sequenced four new reference-quality genomes that cover the three taxonomic orders of living chitons: Deshayesiella sirenkoi (Lepidopleurida), Callochiton septemvalvis (Callochitonida), Acanthochitona discrepans and A. rubrolineata (Chitonida). Callochitonida is sister to Chitonida (Moles et al., 2021) and together these could be considered the Chitonida sensu lato. One aim here was to test the potential genomic differences separating these two orders within Chitonida s.l. We used these four new genomes together with previously published genomic data to reconstruct a phylogeny for the phylum. This allows us to confidently reconstruct ancestral chromosome arrangement of total-group Mollusca, and at different transitions within Polyplacophora.

Results

We sequenced chromosome-level genomes of Deshayesiella sirenkoi (Lepidopleurida) from the Daikoku vent field, Western Pacific Ocean, Callochiton septemvalvis (Chitonida s.l.: Callochitonida) and Acanthochitona discrepans (Chitonida sensu stricto) both from the intertidal of Strangford Lough, N. Ireland, and the congener A. rubrolineata from the intertidal of Qingdao, China. These were sequenced with PacBio HiFi and scaffolded using Hi-C, resulting in high quality assemblies with over 97% BUSCO completeness (94% for Callochiton septemvalvis) and with numbers of chromosomes that differed among genera from 8-13 (Fig. 1, SI Appendix, Fig. S1, Table S1). These species show high levels of heterozygosity, ranging from almost 1% in Deshayesiella to 4.12% in Callochiton (Fig. 1, SI Appendix, Fig. S2).

Fig 1.

Our phylogenetic results confirm the placement of chitons as sister to a monophyletic Conchifera, and we recover the expected topology within Polyplacophora consistent with other recent work using genomic and morphological characters (Fig. 2). Within Conchifera, we confirm the topology of other recent studies (Song et al., 2023); however, our supplementary analyses recovered Scaphopoda sister to Gastropoda but with lower support (SI Appendix, Fig. S3). Comparison with genomes from other molluscan classes shows the molluscan ancestor had a genome composed of 20 linkage groups (Fig. 2), we refer to these as the Molluscan Linkage Groups (MLG) 1-20. Three important fusion events are apparent synapomorphies for Polyplacophora, present in all living chitons and no Conchifera: MLG 4+16+18, MLG 7+10, and MLG 8+9 (Fig. 3, Fig. 4, SI Appendix, Figures S4-S7). There are additional fusions and intra-chromosomal rearrangements that are different in every species sampled.

Fig 2.

Fig 3.

Fig 4.

Discussion

Chiton genomes show frequent and extreme rearrangement

Four new chiton genomes represent the most complete genomes sequenced for the class, adding to several previously published partial or complete genomes (Hong Kong Biodiversity Genomics Consortium et al., 2024; Liu et al., 2023; Varney et al., 2022, 2021) (SI Appendix, Tables S1-S3). Our comparisons of conchiferan and aculiferan (polyplacophoran) linkage groups confirm previous studies that also predicted a haploid karyotype of 20 for the molluscan ancestor based on other metazoans (Simakov et al., 2022).

Previous work demonstrated the variability in chiton karyotypes; species in the clade Chitonida have reported haploid numbers ranging maximally from 6 to 16 (Odierna et al., 2008) with a mode of 11 (Hallinan and Lindberg, 2011). Our data give the first indication for expected chromosome count in a species in Lepidopleurida (1n=11) within this established range. Using five chromosome-level genome assemblies for chitons, we reconstructed the ancestral karyotype for Polyplacophora (more strictly the taxonomic order Neoloricata), and all intermediate phylogenetic nodes to demonstrate the stepwise fusion and rearrangement of gene linkage groups during chiton evolution (Fig. 3).

Chitons demonstrate extreme genome rearrangement, even within a single genus. This represents not only gene order differences but syntenic (co-occurrence) changes in genomic architecture. Species of the genus Acanthochitona have a relatively short divergence time of maximally ∼23 My based on the fossil record (Dell’Angelo et al., 2020). Acanthochitona discrepans and A. rubrolineata each have 8 haploid chromosomes but these result from two different fusions compared to the reconstructed ancestral karyotype of Acanthochitona (Fig. 1, Fig. 4C). Previous karyotype data for A. discrepans were actually based on specimens of A. crinita; the species are very similar but the correct identification can be determined based on the geographic distributions (Vončina et al., 2023). Acanthochitona crinita has 9 haploid chromosomes (Certain, 1951): this is yet another syntenic arrangement within the same genus. There are major changes between congeners in different ocean basins (the Pacific A. rubrolineata) but also between two species in the NE Atlantic (A. discrepans and A. crinita) that are morphologically and ecologically almost indistinguishable.

Living species of Lepidopleurida retain more plesiomorphic morphology and this clade has a deep fossil record extending to the lower Carboniferous (Sigwart, 2017). The largest number of novel fusions among chiton genomes is three, in the lepidopleuran Deshayesiella (Fig. 4A). This in particular contrasts with the naive expectation that gene arrangement in most chiton genomes would be relatively conserved; Lepidopleurida retain plesiomorphic conditions in their morphology (Ampuero et al., 2024; Sigwart, 2017) yet this exemplar shows deviations from the reconstructed ancestral karyotype.

The remaining living chitons (Chitonida s.l.) comprise two sister clades recognised as separate orders: Chitonida and Callochitonida (Moles et al., 2021). Callochiton also has two additional fusion events, and the chitonid Liolophura japonica has a partial genome duplication, with two linkage groups fused apparently at first and then duplicated (Fig. 4B). In our reconstruction of ancestral karyotypes, there is no differences in arrangement between the ancestor of Chitonida sensu stricto or the ancestor of Chitonida sensu lato (Chitonida + Callochitonida). The four species in Chitonida s.l. share a large, fused chromosome (MLG 01+04+16+18) that is notably not well-mixed in Liolophura japonica (i.e., Chr01 in Liolophura japonica, Chr01 in Callochiton septemvalvis, Chr01 in Acanthochitona rubrolineata, and Chr03 in A. discrepans). Part of this pattern has its origin in the ancestral chiton karyotype and is retained in the lepidopleuran Deshayesiella (i.e., MLG 04+16+18). This implies a variable rate of intra-chromosomal rearrangement, with several of the MLGs conserved.

Ancestral karyotypes for Mollusca

Chromosome numbers are not strongly conserved in animals. Changes in chromosome numbers are hypothesized to be an important driver of the diversification of Lepidoptera, with a strong correlation shown between rates of chromosome number changes and speciation (de Vos et al., 2020). The instability in lepidopteran chromosome numbers has mainly been studied from karyotype data but recently confirmed in genomic data (Chen et al., 2019), and changes within genera are known to encompass both neutral and adaptive evolution (Lucek, 2018; Vershinina and Lukhtanov, 2017). However, the level of rearrangement may be less than that in Polyplacophora.

Even within a single species, chromosome fusions are not uncommon, with Roberstonian translocations occurring in roughly 1 out of every 1000 live human births (Wilch and Morton, 2018). A previous study speculated that chromosome loss in related clades of chitons may be the result of Robertsonian translocation (Odierna et al., 2008). Although we still lack information on the telomere, this simple mechanism is not a sufficient explanation for several more complex intrachromosomal rearrangements (e.g. MLG 7+10: Fig. 4). Chromosome fusion does not present any immediate reproductive barrier, yet chromosome rearrangements between sister species can act as Dobzhansky-Muller incompatibilities and generate reproductive isolation (de Vos et al., 2020). Chitons are broadcast spawners, so such barriers may be speculatively advantageous, but so too are bivalves where these syntenic rearrangements are not found.

In order to compare syntenic changes within Polyplacophora to other mollusc clades, we re-analysed all available molluscan genomes in context of the newly identified MLGs. For example, the divergence of the bivalve clade Imparadentia (Lucinida + Venerida) is estimated in the Silurian, ca. 430 Mya (Crouch et al., 2021), or more than 150 My deeper than the divergence of crown group Polyplacophora, yet species within Imparadentia have almost no syntenic rearrangement (SI Appendix, Fig. S7, S8). Bivalves in Pteriomorphia, which have an even deeper origin in the lower Ordovician, ca. 485 Mya, show some rearrangement but members of two orders (Arcida and Pectinida) are highly similar, while two species in Ostreida that are intensively cultivated differ from these others and from each other (SI Appendix, Fig. S8).

Potential adaptive roles may be connected to contrasting mobility of different MLGs. The smallest two MLGs are the most conserved across the five chiton genomes; MLG20 is the only group that remains separate in all bivalve and chiton taxa. This region is also not duplicated in partial genome duplications in Neogastropoda, but is duplicated in whole genome duplication in two hermaphroditic terrestrial gastropods (Achatina fulica and Arion vulgaris) (Fig. 5, SI Appendix, Fig. S5).

Fig 5.

Earlier models based on karyotype data predicted three whole genome duplication events in the evolutionary history of living molluscs: in Neogastropoda, Stylommatophora, and coleoid cephalopods. Available genomes of Stylommatophora confirm whole genome duplication, as well as large scale gene duplications in Neogastropoda (SI Appendix, Fig. S6). In coleoid cephalopod genomes, the co-occurrence of loci from MLGs and intensive fusions (SI Appendix, Fig. S6) might not be easily resolved but may be more likely from the chromosome-disrupted processes reported in previous work (Albertin et al., 2023). By contrast, in Nautilus the linkage group that contains the conserved hox gene sets (Chr11, MLG16) shows no signal of co-occurrence with other MLGs or on other chromosomes (SI Appendix, Fig. S6). New results also show a partial and ancient genome duplication in the chiton Liolophura (Fig. 2, Fig. 4). Given that the MLGs are well conserved in all chitons if not their order (SI Appendix, Fig. S4), the event in Liolophura japonica is most likely a duplication instead of fission. Whole or partial genome duplication is now known from four molluscan classes, with contrasting patterns of ploidy (e.g. gastropods) or tandem (e.g. scaphopods) duplications (Fig. 5).

Reconstructing genome evolution is naturally more difficult for groups with high rates of intra-chromosomal rearrangement (Farré et al., 2019; Muffato et al., 2023). These patterns of co-occurrence of loci on the same chromosomes (synteny) should also persist for longer in evolutionary time, compared to faster rates of change in gene order (Damas et al., 2022; Simakov et al., 2022). Syntenic changes do not necessarily follow overall rates of translocation, which do not differ obviously among molluscs (SI Appendix, Table S4, S5). Confounding effects of rapid gene order or linkage group changes are not limited to issues of phylogenetic distance, when even closely related species have significant syntenic differences. Nonetheless a focus on linkage group arrangements is a promising direction for the many unresolved questions of deep molluscan phylogeny. Recent studies have championed the importance of synteny-based studies on genome evolution, as a basis to understand deep divergences and also the general mechanisms that underlie genomic architecture, fusion, fission, and translocation (Mackintosh et al., 2023). A general trend in insect and vertebrates is for increasing chromosome numbers in derived lineages. Fusion events may be more common in more recently derived chitons or connected to specific adaptations.

Conclusions

The relatively small number of reference quality genomes for most molluscan groups is a temporary limitation to these analyses. High heterozygosity seems to be very common in molluscs, and is clearly a feature of chitons, which causes difficulty for high quality genome assembly. The heterozygosity for Callochiton at 4.12% far exceeds the 2.95% genomic heterozygosity reported as “one of the highest” for Lepidoptera (García-Berro et al., 2023). Large scale analyses based on genomics are equivocal about the drivers of genetic diversity, and further data from diverse clades must be included.

The concept of “diverse groups” is mostly based on perceived variability, typically manifested via morphological, ecological, or genetic differences. As one interesting comparison, Lepidoptera is commonly regarded as super-diverse, but also represent a recognisably constrained form, with variation in striking difference in colours, patterns and wing shapes. While lepidopteran variability is visually conspicuous, chitons are apparently equally variable, as demonstrated through their genetic rearrangement. Chitons, despite exceeding 1000 species in diverse ecological niches, with a wide range of morphological adaptations, are considered a “minor” and “neglected” clade. The group is taxonomically challenging – on the one hand, because of the supposed morphological stasis as a group, on the other, because of their high interspecific variability. Our comparative analyses suggest chromosomal-level changes are a pattern throughout much of chiton evolution, since chromosomal rearrangements are found when comparing congeneric species (Acanthochitona spp., Fig 4C) and also across orders (Fig 4A, 4B).

Chitons exhibit chromosome rearrangements at an almost unprecedented level. This is clear in the present study, the first comparative genomic study for Polyplacophora, compared to orders of magnitude more data for better studied groups. All this poses a more general question on how we define variability, how we perceive it, and also, if we actually understand it at all? Recognition of variability in the genome and phenome is crucial in understanding and re-evaluating unbiased measures of diversity in overlooked groups of organisms. The relationship revealed between the unstable arrangement of chiton genomes and species diversity provides new insights into potential mechanisms for speciation and broader diversification within Mollusca.

Methods

Specimens were collected alive and flash frozen in liquid nitrogen (A. rubrolineata) or frozen at -80°C and stored at -80°C. High molecular weight genomic DNA was extracted from the foot of an individual specimen of each species, following the guidelines of the SDS method, and used for PacBio high fidelity (HiFi) sequencing and Hi-C library preparation. The PacBio Sequel II / IIe instrument was used for sequencing in CCS mode. For transcriptome sequencing, five separate tissues were dissected from the same specimen as DNA extraction, and preserved in RNAlater: foot (F), perinotum (P), radula sac (R), shell edge (S), and visceral mass (V).

Adapters in the reads were checked and removed using HiFiAdapterFilt version 2.0. The genome survey was finished using jellyfish v2.2.10 and Genomescope v2.0 (Ranallo-Benavidez et al., 2020), implemented 21- and 23-mer, which produced the estimated genome size and heterozygosity level for each species (SI Appendix, Fig. S2). The four genomes were assembled de novo based on qualified HiFi reads using hifiasm v0.16.1 (Cheng et al., 2021), with some additional curation steps for C. septemvalvis because of the high heterozygosity (SI Appendix). Potential contamination in the genome was detected and removed using blobtools v1.1.1 (Laetsch & Blaxter, 2017). Duplicated haplotigs and overlaps were removed for D. sirenkoi, A. discrepans and A. rubrolineata using Purge_Dups v1.2.6 (Guan et al., 2020).

The aligner STAR version 2.7.10a or STAR v2.7.3a was employed to map RNA-seq data into the working genome data (Dobin et al., 2012). The resulting alignment file served as an important support in prediction methods. Ab initio gene prediction was performed on the repeat-masked assembly in comparison with other published genomes (SI Appendix). Summary information on assembly, gene model, and annotation is provided in Table S2.

We identified putative orthologous sequences shared among genomes or transcriptomes for 58 molluscs and 4 additional metazoans to generate an alignment for phylogenetic analysis using VEHoP (Li et al., 2024). The alignments were removed if the overlap among them was less than 20 amino acids and there were less than 75% of taxa sampled, then, each alignment was used to construct “approximately maximum likelihood” tree using FastTree version 2.1.11 (Price et al., 2010). Phylogenetic relationships were investigated using IQ-TREE version 2.1.3 (Minh et al., 2020) with the “-MFP” model to compute the best-fit model of each partition and 1000 ultrabootstraps to test the topological support. To test the impact on conchiferan topology we ran a second analysis excluding Solemyida, the earliest diverging clade in Bivalvia (SI Appendix, Fig. S3).

The method for identifying ancient linkage groups followed previously published works. We resconstructed the linkages of the orthologues in molluscs, with the demonstrated commands, and conserved and representative proteins and their presumptive locations (https://github.com/ylify/MLGs). In detail, the ancient and conserved linkage groups were inferred from 7 chromosome-level mollusk assembles, including 1 gastropod (Gibbula magus), 1 bivalve (Mizuhopecten yessonensis), and the 5 chitons (A. rubrolineata, A. discrepans, C. septemvalvis, D. sirenkoi, and Liolophura japonica). The obtained gene set consisted of 4729 homologs, which could be detected in all 7 assembles and located in 20 linkage groups. The ancestral states of nodes were predicted to trace the karyotype evolutionary route in Polyplacophora, from the common ancestor of Mollusca to the genus Acanthochitona. Similarly, the nodes with more than two chromosome-level genomes within Mollusca were also investigated for the chromosome number of their common ancestor. We also compared the translocation rates across Mollusca, based on the non-syntenic rate of change divided by estimated divergence time. In this study, 25 genomes from four classes were selected to calculated the translocation rate at the inter-chromosome level, including 2 scaphopods, 5 chitons, 8 bivalves, and 10 gastropods.

Availability of data and materials

The four chiton genomes projects have been deposited with the NCBI BioProject, with Acanthochitona discrepans in PRJNA1114954, A. rubrolineata in PRJNA1114370, Callochiton septemvalvis in PRJNA1114372, Deshayesiella sirenkoi in PRJNA1114373. In each project, the whole-genome sequencing data, Hi-C data, RNA-seq data have been affiliated. The genome assembly and gene-model predictions are deposited at figshare (10.6084/m9.figshare.27894189). The constructed Molluscan Linkage Groups (MLGs) are available, with the conserved and representative sequences and their predictive corresponding locations (https://github.com/ylify/MLGs/).

The commands used in this study have been deposited on GitHub (https://github.com/ylify/MLGs/), including the genome assembly, repeats and coding regions prediction, phylogenomic inference, and mutual best hits and ancient linkage groups detections.

Acknowledgements

We thank Carola Greve, Damian Baranski, Alexander Ben Hamadou, and Charlotte Gerheim of the Translational Biodiversity Genomics (TBG) project based in the Senckenberg Research Institute and Museum, Frankfurt, for support with lab work and sequencing. We thank the staff of the Queen’s University Marine Laboratory, Portaferry, N Ireland, and Chong Chen (JAMSTEC) for support with field work and specimen collection. This is contribution number 28 of the Senckenberg Ocean Species Alliance. The sampling of Deshayesiella sirenkoi was carried out under the Marine Scientific Research permit number MSR U2022-047 from the United States government, since South Chamorro Seamount is within the Commonwealth of the Northern Mariana Islands.

Additional information

Funding

This study was supported by the Natural Science Foundation of Shandong Province (ZR2023JQ014) and the Fundamental Research Funds for the Central Universities (202172002 and 202241002). This research was further supported by a generous philanthropic donation to the Senckenberg Geselleschaft für Naturforschung that funds the Senckenberg Ocean Species Alliance, and by Leibniz Association project PHENOME (P123/2021).

Author contributions

JDS and JS conceived the study and collected the samples. YL and ZC performed the analyses. JDS and KV led the biological interpretation of the results and review of comparative information from other studies. JDS drafted the manuscript. All authors contributed to interpreting the results and writing the manuscript.

Additional files

Supporting Information