Introduction

The genomic mechanisms that enable or limit the evolution of major 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 (1). Despite playing a key role in early advances, this important phylum lags in terms of quality and taxonomic coverage of whole genome sequence data (2). Reconstructing the evolution of genome architecture through deep time in diverse molluscs remains critical to understand genome evolution in animals.

Prior work has reasonably assumed that rates of intra-chromosomal gene translocation are constant within major groups(3). 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(4). 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(5). Chitons are increasingly important for bio-inspired design, which will benefit from genomic tools to understand genetic control of their flexible armour, unique sensory systems, and iron mineralised radula(6). Chitons have separate sexes and most reproduce by broadcast spawning, the animals are grazers and most are not highly ecologically specialised(7). Chitons are generally conserved, yet the known species richness in chitons is higher than for the more morphologically and behaviourally diverse cephalopods. The adaptations that drive speciation in this group remain an open question.

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). We use these 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 & 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(6, 8–10)(SI Appendix, Tables S1-S3). We sequenced chromosome-level genomes of Deshayesiella sirenkoi (Lepidopleurida) from the Daikoku vent field, Western Pacific Ocean, Callochiton septemvalvis (Callochitonida) and Acanthochitona discrepans (Chitonida) 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(11); 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. This confirms previous studies that also predicted a haploid karyotype of 20 for the molluscan ancestor based on other metazoans(12). 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.

Previous work demonstrated the variability in chiton karyotypes; species in the clade Chitonida have reported haploid numbers ranging maximally from 6 to 16(13) with a mode of 11(1). 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(14). 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(15). Acanathochitona crinita has 9 haploid chromosomes(16): 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 Carboinferous(5). The largest number of novel fusions among chiton genomes is three, in the lepidopleuran Deshayesiella (Fig. 4A). The remaining living chitons form two sister clades recognised as separate orders: Chitonida and Callochitonida(17). 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 or the ancestor of (Chitonida + Callochitonida). The four species in Chitonida (sensu lato, including Callochiton) 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(18). The instability in lepidopteran chromosome numbers has mainly been studied from karyotype data but recently confirmed in genomic data(19), and changes within genera are known to encompass both neutral and adaptive evolution(20, 21). 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(22). A previous study speculated that chromosome loss in related clades of chitons may be the result of Robertsonian translocation(13). 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(18). 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(23), 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 (1). Available genomes of Stylommatophora confirm whole genome duplication (24), 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 (25). By contrast, in Nautilus the linkage group that contains the conserved hox gene sets (Chr11, MLG16) (26) 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(27, 28). 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(12, 29). 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(3). A general trend in insect and vertebrates is for increasing chromosome numbers in derived lineages (29, 30). Fusion events may be more common in more recently derived chitons or connected to specific adaptations.

Understanding diversification requires diverse clades

Chitons 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(31). Chromosome rearrangement is an attractive speculative explanation that could support maintaining species boundaries in this clade, for example Acanthochitona discrepans (1n=8) and A. crinita (1n=9) (16). 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(32). Large scale analyses based on genomics are equivocal about the drivers of genetic diversity (33, 34), and further data from diverse clades must be included.

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.

The concept of “diverse groups” is mostly based on perceived variability, typically manifested via morphological, ecological, or genetic differences. Lepidopterans are commonly regarded as super-diverse, but their striking colours and different wing shapes are variations within a clearly constrained theme. 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 (14). Chitons exhibit chromosome rearrangements, even at the genus level within Acanthochitona, at an almost unprecedented level. This is clear in the first comparative genomic study, 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? While lepidopteran variability is conspicuous, chitons are equally variable, clearly demonstrated through their genetic rearrangement. Recognition of variability in the genome and phenome is crucial in understanding and re-evaluating unbiased measures of diversity in overlooked groups of organisms.

Materials and 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. The qualified HiFi reads (i.e., adapter-free) were processed for de novo genome assembly using hifiasm v0.16.1-r375 (35). 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 (36), 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 (35), 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 (37). Duplicated haplotigs and overlaps were removed for D. sirenkoi, A. discrepans and A. rubrolineata using Purge_Dups v1.2.6 (38).

The aligner STAR version 2.7.10a or STAR v2.7.3a (39) was employed to map RNA-seq data into the working genome data. 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. 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 (40). Phylogenetic relationships were investigated using IQ-TREE version 2.1.3 (41) 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 (12, 42, 43). 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.

Data availability

The four chiton genomes projects have been deposited with the NCBI BioProject, with Acanthochitona discrepans in PRJNA1114954, A. rubrolineatain 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 also available on GitHub, as well as the constructed Molluscan Linkage Groups (MLGs) is also available, with the conserved and representative sequences and their predictive corresponding locations (https://github.com/ylify/MLGs/).

Code availability

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 antient 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.

This study was supported by the Natural Science Foundation of Shandong Province (ZR2023IQ014) 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).