Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

The Scaly-foot Snail genome and implications for the origins of biomineralised armour

The Scaly-foot Snail genome and implications for the origins of biomineralised armour ARTICLE https://doi.org/10.1038/s41467-020-15522-3 OPEN The Scaly-foot Snail genome and implications for the origins of biomineralised armour 1,10 2,10 2,10 3 4,5 3 Jin Sun , Chong Chen , Norio Miyamoto , Runsheng Li , Julia D. Sigwart , Ting Xu , 1 1 3 1 1 6 Yanan Sun , Wai Chuen Wong , Jack C. H. Ip , Weipeng Zhang , Yi Lan , Dass Bissessur , 2,7 2 2 8 8 Tomo-o Watsuji , Hiromi Kayama Watanabe , Yoshihiro Takaki , Kazuho Ikeo , Nobuyuki Fujii , ✉ ✉ 9 3 2 1 Kazutoshi Yoshitake , Jian-Wen Qiu , Ken Takai & Pei-Yuan Qian The Scaly-foot Snail, Chrysomallon squamiferum, presents a combination of biomineralised features, reminiscent of enigmatic early fossil taxa with complex shells and sclerites such as sachtids, but in a recently-diverged living species which even has iron-infused hard parts. Thus the Scaly-foot Snail is an ideal model to study the genomic mechanisms underlying the evolutionary diversification of biomineralised armour. Here, we present a high-quality whole- genome assembly and tissue-specific transcriptomic data, and show that scale and shell formation in the Scaly-foot Snail employ independent subsets of 25 highly-expressed tran- scription factors. Comparisons with other lophotrochozoan genomes imply that this biomi- neralisation toolkit is ancient, though expression patterns differ across major lineages. We suggest that the ability of lophotrochozoan lineages to generate a wide range of hard parts, exemplified by the remarkable morphological disparity in Mollusca, draws on a capacity for dynamic modification of the expression and positioning of toolkit elements across the genome. Department of Ocean Science, Division of Life Science and Hong Kong Branch of the Southern Marine Science and Engineering Guangdong Laboratory (Guanzhou), The Hong Kong University of Science and Technology, Hong Kong, China. X-STAR, Japan Agency for Marine-Earth Science and Technology (JAMSTEC), 2-15 Natsushima-cho, Yokosuka, Kanagawa 237-0061, Japan. Department of Biology, Hong Kong Baptist University, Hong Kong, China. 4 5 6 Marine Laboratory, Queen’s University Belfast, Portaferry, N. Ireland. Senckenberg Museum, Frankfurt, Germany. Department for Continental Shelf, Maritime Zones Administration & Exploration, Ministry of Defence and Rodrigues, 2nd Floor, Belmont House, 12 Intendance Street, Port-Louis 11328, 7 8 Mauritius. Department of Food and Nutrition, Higashi-Chikushi Junior College, Kitakyusyu, Japan. National Institute of Genetics, 1111 Yata, Mishima, Shizuoka, Japan. Graduate School of Agricultural and Life Sciences, The University of Tokyo, 1-1-1, Yayoi, Bunkyo-ku, Tokyo, Japan. These authors contributed equally: Jin Sun, Chong Chen, Norio Miyamoto. email: kent@jamstec.go.jp; boqianpy@ust.hk NATURE COMMUNICATIONS | (2020)11:1657 | https://doi.org/10.1038/s41467-020-15522-3 | www.nature.com/naturecommunications 1 1234567890():,; 10 MB 10 MB 20 MB 10 MB 10 MB 30 MB 10 MB 20 MB 10 MB 40 MB 30 MB 30 MB 20 MB 20 MB 10 MB 10 MB ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15522-3 he appearance of biomineralised skeletons in the Cambrian Living hydrothermal vent taxa represent Cenozoic radiations 9,10 precipitated an evolutionary arms race and the original (≤66 million years ago, Mya), including peltospirid gastropods Texplosive diversification of modern animal forms . Since such as the Scaly-foot Snail; the Scaly-foot Snail represents a then, a long history of body plan modifications, reversals, and recent evolution of a complex scleritome. On its discovery in the convergences has obscured the relationships among some animal iron-rich Kairei hydrothermal vent field in the Indian Ocean in lineages . Understanding the genomic toolkit that enabled inno- 2001, the complex armature of the Scaly-foot Snail was compared vations in skeletons and armour is critical to reconstructing the with apparently plesiomorphic scleritomes found in aculiferan 3,4 early radiation of major clades . In particular, Lophotrochozoa molluscs and important early fossil taxa such as Halkieria and (including annelids, brachiopods, molluscs and others) presents a other sachtids . Later, a second population lacking iron sulfide significant and unresolved problem for phylogenetic reconstruc- mineralisation was discovered in the Solitaire hydrothermal site tion, and controversies over the affinities of key fossils . The characterised by low concentrations of iron (1/58 of that found at Scaly-foot Snail, Chrysomallon squamiferum, is unique among Kairei ), which presents an opportunity to understand physio- gastropod molluscs in having dense, imbricating chitinous logical and molecular mechanisms behind its unique miner- sclerites covering the whole distal surface of the soft foot , alisation pattern (also see Supplementary Note 1). Recent results forming a dermal scale armour in addition to a solid calcium indicated that the Scaly-foot Snail mediates its scale biominer- carbonate coiled shell (Fig. 1). These hard parts, including alisation by supplying sulfur through nano-scale channel-like sclerites and the shell, are often mineralised with iron sulfide, columns in the scales, which reacts with iron ions diffusing making it the only known metazoan using iron as a significant inwards from the vent fluid to make iron sulfide nanoparticles . component of skeleton construction . Many sachtids also had tissue-filled canals within their sclerites, Total Mantle (shell) Scales Novel genes Fig. 1 Key features of the Scaly-foot Snail Chrysomallon squamiferum genome. Circos plot showing the 15 pseudo-chromosomal linkage groups; with the Scaly-foot Snail at centre. The outer ring (red peaks) indicates gene density in each pseudo-chromosome, and the inner rings shows the normalized density of the highly expressed genes in the shell-secreting mantle (black peaks) and scale-secreting epithelium (white peaks, semi-transparent overlaid on black mantle peaks) as well as the density of novel genes (grey). The expression level is the average fold change of the target tissue versus the other four types of tissues (n = 5); the sliding window size is 10 kb. MB, megabases. Source data are provided as a Source Data file. 2 NATURE COMMUNICATIONS | (2020)11:1657 | https://doi.org/10.1038/s41467-020-15522-3 | www.nature.com/naturecommunications 0.6 0.4 0.2 0.6 0.4 0.2 0.8 0.6 0.4 0.2 30 MB 30 MB 20 MB 20 MB 10 MB 10 MB 20 MB 40 MB 10 MB 30 MB 20 MB 20 MB 10 MB 10 MB 10 MB 10 MB 0 NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15522-3 ARTICLE ab c Scaly-foot Snail Other Mollusca 100% Lophotrochozoa 80% 60% 4 40% 20% 85% 90% 95% 100% Scales Mantle Oes. Gland Gill 85% 90% 95% 100% BUSCO BUSCO Annotated Novel genes (Genome completeness) (Genome completeness) Fig. 2 Summary of the Scaly-foot Snail genome. a, b Quality comparisons between the Scaly-foot Snail genome and other available lophotrochozoan genomes: a Contig N50 vs BUSCO (Benchmarking Universal Single-Copy Orthologs); b Number of scaffolds vs BUSCO (for full comparison see Supplementary Table 2). Red dot indicates the Scaly-foot Snail genome, blue dots indicate other molluscan genomes, orange dots indicate other non- molluscan lophotrochozoan genomes. c Proportion of annotated genes (pink) and novel genes (grey) in four key tissue types of the Scaly-foot Snail, including scale-secreting epithelium, shell-secreting mantle, oesophageal gland (oes. gland), and gill. Source data are provided as a Source Data file. although entirely different in proportion and arrangement. The The Scaly-foot Snail genome represents the most complete channels within the scales of the Scaly-foot Snail are likely linked and continuous genome among the assembled mollusc an or to another key adaptation, in that it hosts sulfur-oxidising bac- lophotrochozoan genomes to date (Fig. 2), with a metazoan teria within cells of a highly vascularised, hypertrophied oeso- BUSCO (Benchmarking Universal Single-Copy Orthologs) score 13,14 phageal gland , and the sulfur may originate as metabolites of 96.6% for the genome assembly and 87.5% for the predicted from the endosymbionts . Taking these observations together, it transcripts. This genome assembly represents a successful is unclear whether the evolution of the sclerites in the Scaly-foot application of methods originally developed for model organ- Snail should be interpreted as a recurring ancestral phenome, or a isms to a little-studied taxon, providing a benchmark in quality recently derived adaptive novelty. compared with other published lophotrochozoan genomes The search for a biomineralisation toolkit underlying hard part confounded by heterozygosity and high repeat contents (Fig. 2a, evolution in molluscs and lophotrochozoans has previously b; Supplementary Table 2). focused on molluscan mantle gene expression and shell forma- tion . Molluscs have repeatedly ‘invented’ additional sclerite-like hard structures, in chitons, aplacophorans, Chrysomallon, and Gene family analyses. The Scaly-foot Snail has fewer novel gene also in other gastropods, cephalopods, and bivalves . Here, we families than other lophotrochozoans with high-quality genomes use a complete genome assembly of the Scaly-foot Snail and available (such as Pomacea canaliculata, Mizuhopecten yessoensis, tissue-specific transcriptomics to compare with data from other Lingula anatina, Achatina fulica, Sinonovacula constricta and lophotrochozoans in order to test whether there is indeed a Capitella teleta: Fig. 3a; Supplementary Table 2), which appears to universal biomineralisation toolkit in Mollusca or Lopho- conflict with an expectation that evolutionary novelties are trochozoa, and to identify specific genomic tools that enable the associated with novel genes (Fig. 3a). To corroborate this, only Scaly-foot Snail to repeatedly modify and duplicate hard parts. 11% of the Scaly-foot Snail gene families are not found among four other lophotrochozoan genomes, while in the other taxa, which apparently lack the dramatic morphological novelties of Results the Scaly-foot Snail, this figure is over 20% and up to 35% Genome assembly. The genome of a single specimen of Chry- (Fig. 3a). In addition, only 4.8% of gene families in the Scaly-foot somallon squamiferum (Fig. 1, Supplementary Fig. 1, Supple- Snail genome were found to be novel to gastropods, while as mentary Table 1) collected from the Kairei hydrothermal vent much as 87.0% of them may be unique to Lophotrochozoa or field was sequenced with a combination of Oxford Nanopore even pre-date the origin of Lophotrochozoa (Fig. 3b). Technologies and Illumina platforms (Supplementary Note 2; Comparative analyses among available lophotrochozoan gen- Supplementary Table 2). Sulfur-oxidising endosymbionts within omes (n = 15; Supplementary Table 4) showed that 351 gene the body are a source for significant potential contamination families were significantly expanded in the Scaly-foot Snail. A GO and therefore bacterial endosymbiont genomes were removed enrichment analysis on the expanded gene families revealed 30 from downstream analyses (Supplementary Note 2). The genome overrepresented GO categories in the Scaly-foot Snail genome of the Scaly-foot Snail is relatively compact for Mollusca (444.4 (Supplementary Fig. 2, Supplementary Data 1). These expanded Mb; Supplementary Table 1) and is highly heterozygous (1.38%) gene families appear to be involved in the secretion process of but with relatively low repeat content (25.2%). With additional proteinaceous materials (e.g. scavenger receptor activity, carbo- Hi-C data using a second specimen from the same population, hydrate binding, chitin-related metabolic process and chitin 1032 contigs (N50 = 1.88 Mb) were anchored to 15 pseudo- binding) and symbiosis (regulation of innate immune response chromosomal linkage groups (Fig. 1). A total of 16,917 gene and symbiont-containing vacuoles), suggesting their contribu- models were predicted (85.7% comparatively annotated; 2415 tions to both biomineralisation and regulation of endosymbiosis. additional novel genes) with evidence from transcripts, homo- The expanded genes were biased in distribution across the logue proteins, and ab initio methods (Supplementary Note 2). chromosomes, with Chr11 and Chr12 being especially enriched The number of genes in the Scaly-foot Snail genome appears to be (Fold change > 2 and FDR < 0.01, Supplementary Note 3). low, and the fact that 97.3% of the de novo assembled tran- We recovered the complete Antennapedia (ANTP) Hox gene scriptome could be mapped to the genome indicates that the complex containing 11 Hox genes on Chr11 (Fig. 4a) as in Lottia 19 20 genome assembly has high completeness. gigantea (a gastropod) and Mizuhopecten yessoensis (a NATURE COMMUNICATIONS | (2020)11:1657 | https://doi.org/10.1038/s41467-020-15522-3 | www.nature.com/naturecommunications 3 Contig N50 Scaffolds s 10 K 12 K 8 K 10 K 6 K 8 K 4 K 6 K 2 K 4 K 0 K 2 K 8 K 0 K 16 K 6 K ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15522-3 a b Eumetazoa Bilateria Lophotrochozoa Mollusca Gastropoda (11%) 76 73 202 45 2708 85 (26%) (21%) 106 127 380 277 2614 364 (35%) (20%) Fig. 3 Gene family analyses of lophotrochozoan genomes. a Venn diagram showing the number of shared and unique gene families among five lophotrochozoan genomes. b Circos plot showing the proportion and origin of shared gene families across five lophotrochozoan taxa. Arc values correspond to the number of gene families. Genes shared across Eumetazoa are indicated by grey lines, Bilateria by green lines, Lophotrochozoa by orange lines, Mollusca by blue lines, and Gastropoda by purple lines. K, kilo (1000). Source data are provided as a Source Data file. Mantle Scales Gill OG DMBT1 (scavenger domain) x 65 Gbx Evx Hox1 Hox2 Hox3 Hox4 Hox5 Lox5 Antp Lox4 Lox2 Post2 Post1 En En En Mox Mox Mnx DLx Chr11, 17.5M Gbx Evx Hox1 Hox2 Hox3 Hox4 Hox5 Lox5 Antp Lox4 Lox2 Post2 Post1 En En En Mox Mox Mnx DLx Mantle Scales Gill OG Hox4 Hox1 Evx Hox3 Fig. 4 Tissue-specific gene expression in the Scaly-foot Snail. a Arrangement of DMBT1 gene tandem repeats (grey), Hox genes and Hox-like genes (red) on the pseudo-chromosome Chr11 and their expression patterns in four types of tissues including mantle, scale-secreting epithelium, gill, oesophageal gland (OG); darker shades indicate higher levels of expression. SRCR (scavenger receptor cysteine-rich) domain, scavenger receptor cysteine-rich domain. Source data are provided as a Source Data file. b Tissue-specific expression patterns of transcriptome factors shown with in situ hybridisation, including Hox4 and Hox1 in the mantle edge (arrowheads; additional tissue is visible in the background including scales also expressing Hox1) and tissue sections of the scale-secreting epithelium expressing Evx and Hox3. Scale bars = 500 μm for mantle and 100 μm for scales. In situ hybridisation experiments were repeated independently on three individuals (twice on each individual for Hox1 and Hox3, once for Hox4 and Evx) with similar results. 4 NATURE COMMUNICATIONS | (2020)11:1657 | https://doi.org/10.1038/s41467-020-15522-3 | www.nature.com/naturecommunications 14 K 4 K 12 K 2 K 0 K 10 K 10 K 8 K 8 K 6 K 6 K 4 K 4 K 2 K 2 K 0 K 0 K –1 0 1 2 3 Z score P NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15522-3 ARTICLE bivalve); among molluscs, the intact ordered cluster has only been ‘aculiferan’ molluscs (i.e. chitons and the worm-like Caudofoveata recovered in one gastropod (Lottia), and two bivalve (Mizuho- and Solenogastres). The Scaly-foot Snail possesses an expanded 21 4 pecten, and Azumapecten farreri) genomes , but are critical to suite of mineral-secreting tissues , in the foot epidermis as well as understanding body plan evolution . The derived Scaly-foot the shell-forming mantle tissue. Many distinctive morphological Snail homeobox genes show reorganisation by intra- adaptations in Chrysomallon are connected to supporting endo- chromosomal rearrangement (Supplementary Figs. 3, 4). Syntenic symbiotic microbes housed in an enlarged oesophageal gland, a comparisons reveal both extensive inter- and intra-chromosomal condition resulting from rapid evolution . Indeed, even the rearrangements between the Scaly-foot Snail and other molluscs scales of the Scaly-foot Snail are apparently related to supporting (Supplementary Fig. 3). its endosymbionts . Gene expression levels in different tissue types of the Scaly-foot To assess tissue-specific functions, particularly in relation to Snail were assessed by transcriptome sequencing (five specimen biomineralisation, a paired test in DESeq2 was applied to replicates) using samples fixed in situ in the snail’s native deep- distinguish tissue highly expressed genes (n = 5, Figs. 1, 4a, and sea vent habitat to avoid the confounding effect of environment Supplementary Data 5–7), further confirmed with real-time PCR and hydrostatic pressure changes on transcription. The five (two replicates, five tissue types, Supplementary Note 5) and primary tissue types targeted included the shell-secreting mantle, in situ hybridisation (ISH; for mantle and scale-secreting scale-secreting epithelium, oesophageal gland where the bacter- epithelium; Figs. 4b, 5 and Supplementary Fig. 6). Given the iocytes housing the endosymbionts are located, the gill, and the high quality of the Scaly-foot Snail genome, genes related to foot musculature. biomineralisation can be accurately linked to their organisation in Among the expanded genes on Chr11 (262 genes), one that the genome and their respective expression levels in relevant controls protein binding to proteoglycan (for instance, chitin) tissues, i.e. shell-secreting mantle and scale-secreting epithelium. known as DMBT1 with scavenger receptor cysteine-rich (SRCR) Visualisation of expression levels by genome organisation domains was found to be highly expressed in both the scale- revealed similar expression patterns for scale-secreting epithelium secreting epithelium and shell-secreting mantle (Fig. 4a) with up to and mantle at chromosomal level across the genome, also 65 tandemly duplicated paralogues compared with one or two revealing a pattern where the pseudo-chromosome Chr11, which copies in other molluscs . Proteins with SRCR domains are contains the abovementioned ANTP Hox cluster as well as the commonly detected in shell matrix proteins . These gene copies DMBT genes, contained a high density of highly expressed genes were highly expressed in the mantle, but were also especially highly in these two biomineralising tissues (Fig. 1). The expression expressed in the scale-secreting epithelium (Fig. 4a). The Scaly-foot profiles of biomineralisation genes in the scale-secreting epithe- Snail DMBT1 gene copies had an average of 2.3 SRCR domains, lium was most similar to those of the shell-secreting mantle tissue but the number varied greatly among individual paralogues (0–5; (81 shared highly expressed genes, 11.2% of the scale-secreting Supplementary Data 2). As SRCR domains are known to evolve epithelium or 6.8% of mantle highly expressed genes); both rapidly to generate high gene diversity, the variation of the domain expressed common genes involved in biomineralisation (Supple- numbers in each paralogue likely reflects rapid evolution and mentary Data 5, 6), in separate tissues with different and complex expansion of this gene to coordinate protein secretion of dermal functions. scales and the shell periostracum. A gene tree of the Scaly-foot We identified 25 highly expressed transcription factors in the Snail DMBT1 shows that the paralogues segregate into two clades, scale-secreting epithelium (12) and mantle (13), often confirmed with only one copy from each clade being present in the by ISH (Figs. 4b, 6 and Supplementary Data 5, 6); the two vetigastropod Haliotis genome (Supplementary Fig. 5). biomineralisation tissues did not share any highly expressed Molluscan shell matrix proteins are commonly rich in transcription factors (Fig. 3c). Several of these were only detected repetitive low-complexity domains (RLCDs) , which have been as highly expressed in the Scaly-foot Snail hard parts among the attributed to a number of functions such as chitin or calcium available lophotrochozoan tissue-specific data (Fig. 3). However, binding and structural support, but such domains are lacking in all transcription factors reported to date from other lophotro- DMBT1 paralogues in the Scaly-foot Snail (Supplementary chozoan skeletal elements (brachiopod mantle and chaetae, and Data 3). Analyses of RLCD-rich genes across the Scaly-foot Snail annelid chaetae, n = 10 genes) were found to be also highly genome revealed that the ratio of RLCDs in genes highly expressed in the Scaly-foot Snail, in either the mantle (n = 6 expressed in the scale-secreting epithelium is similar to the overall genes) or scale-secreting epithelium (n = 4). Some of these co- ratio of the whole genome, whereas a higher ratio is seen in the opted transcription factors, and others, have also been reported mantle. This indicates that RLCDs likely contribute to shell from other molluscs including aculiferan mantle and sclerite- biomineralisation like known for other molluscs, but not the secreting tissue, gastropod or bivalve mantle tissue, or tissues 21,22 scales; suggesting that RLCDs may be linked with calcium secreting the radula, beak, or operculum (Fig. 6b and carbonate biomineralisation which is lacking in the scales Supplementary Table 4). The involvement of these transcription (Supplementary Table 3). factors in the formation of various hard parts across different A further 79 gene families were found to be significantly lophotrochozoan groups ranging widely from plesiomorphic (e.g. contracted in the Scaly-foot Snail genome relative to other brachiopod shell) to more recently evolved armature (e.g. Scaly- lophotrochozoans, including several enzymes that play critical foot Snail scales) signify that they together comprise the roles in steroid and amino acid synthesis, in keeping with similar biomineralisation toolkit preserved from the ancestral lopho- gene reductions in other chemosymbiont-dependent holobionts trochozoan genome. (Supplementary Data 4). Numerous enriched GO terms (Supple- Taxon coverage of high quality lophotrochozoan genomes is mentary Note 4) in highly-expressed genes in the oesophageal currently too limited to make significant phylogenetic inferences. gland are relevant to symbiosis, including oxidation-reduction Internal nodes within our time calibrated phlylogenetic recon- process compensating the reactive oxygen species generated by struction confirm and refine some previous divergence estimates. the endosymbiont . For example, Pteriomorphia, the best-sampled lineage within Mollusca, is recovered with a divergence estimate around 470 Mya (483.5–465.1 Mya) which corresponds to the early fossil Biomineralization toolkit. Until the early 2000s, the only living record for the group . Dating of earlier nodes (origins of molluscs known to regularly possess imbricating sclerites were Gastropoda, Bivalvia and Lophotrochozoa), and the topology, are NATURE COMMUNICATIONS | (2020)11:1657 | https://doi.org/10.1038/s41467-020-15522-3 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15522-3 a bc Azumapecten farreri Chitin binding Peritrophin-A 32,768 Crassotrea gigas Pinctada fucata Lingula anatina Pomacea canaliculata 100 Chrysomallon squamiferum CsqKR_Scaf21_18.24 Chrysomallon squamiferum CsqKR_Scaf21_18.25 Haliotis discus hannai Aplysia californica Scale Mantle Gill Foot OG Radix auricularia Tissue Radix auricularia Aplysia californica Pif protein Pomacea canaliculata de Lanistes nyassanus Chrysomallon squamiferum CsqKR_Scaf193_7.1 7 Lottia gigantea Bathymodiolus platifrons 94 Modiolus philippinarum Crassotrea gigas 100 66 Pinctada fucata Azumapecten farreri Mizuhopecten yessoensis 0.03 Modiolus philippinarum Scale Mantle Gill Foot OG Chitin binding Peritrophin-A Octopus bimaculoides von Willebrand factor typeA Tissue 0.3 Fig. 5 Characterisation of chitin binding peritrophin-A gene and pif gene in Chrysomallon squamiferum.a Phylogenetic analysis of chitin binding peritrophin-A and pif genes with names in red indicating genes from the Scaly-foot Snail, and domain architecture of the specific genes displayed on the right of names (purple indicates chitin binding peritrophin-A and blue indicates von Willebrand factor type A). b In situ hybridization of chitin binding Peritrophin-A in scales. c Boxplot showing the expression level of chitin binding peritrophin-A in five different tissue types including mantle, scale-secreting epithelium, gill, oesophageal gland (OG), and foot. The box depicts first to third quartile with whiskers indicating maximum and minimum expression levels, and the centre line refers to the median. Source data are provided as a Source Data file. d In situ hybridization showing the expression of pif gene in mantle. e Boxplot showing the expression level of pif in five different tissues types. Within the Boxplot, the centre line refers to the median, and the boxplot depicts the first to the third quartile, with the whiskers indicating maximum and minimum expression levels. likely confounded by limitations of the available data. Within genes (i.e, Scaly-foot Snail synapomorphies), or shared with, but gastropods, we recovered Neomphaliones (represented by as yet undiscovered in, other related taxa. Chrysomallon) sister to Vetigastropoda, and this clade sister to Genes known to produce proteins involved in molluscan shell Patellogastropoda (Fig. 6a). One recent phylotranscriptomics formation, such as pif, chitin-binding peritrophin-A domain 29 32 study similarly found Vetigastropoda sister to Patellogastro- , were also conserved and highly gene, and chitin synthase poda, but that analysis did not include any examples of expressed in the shell-secreting mantle and scale-secreting Neomphaliones. Neither our genome tree nor any phylotran- epithelium of the Scaly-foot Snail (Fig. 5 and Supplementary scriptomic analyses to date have included all gastropod subclasses Fig. 6). A GO enrichment analysis on highly expressed genes in as no whole genome is available for Neritimorpha . Early multi- the mantle and scale-secreting epithelium revealed that the gene phylogenetic analyses repeatedly encountered severe long categories integral component of membrane, scavenger receptor branch attraction artefacts particularly from patellogastropod activity, and cell-matrix adhesion were enriched in both of these sequences . A recent mitogenome analyses including all biomineralisation tissues (Supplementary Note 4), further gastropod subclasses found Patellogastropoda as the earliest- indicating that they elements of a shared biomineralisation 10,30 branching lineage within living gastropods , which is the toolkit. topology predicted by morphology . The debate on internal These possible downstream biomineralisation genes may be relationships among gastropods is likely to continue until more controlled by the transcription factors in the biomineralisation representatives from all subclasses become available for robust toolkit, and in a number of cases their specific expression in genome-based phylogeny. the scale-secreting epithelium were confirmed by ISH (Fig. 4b). Highly expressed ‘novel’ or taxon-specific genes were dis- The gene exhibiting the highest expression level among the scale- tributed across all tissue types in the Scaly-foot Snail. Although secreting epithelium highly expressed genes was chitin-binding there were relatively slightly more novel genes in biomineralising peritrophin-A. This gene (CsqKR_Scaf21_18.25), together with tissues, the presence of novel genes was not restricted to areas of its recent duplicated paralogue (CsqKR_Scaf21_18.24), were also obvious morphological adaptations. There were fewer novel genes highly expressed in the shell-secreting mantle but at lower levels in the scale-secreting epithelium than the shell-secreting mantle, than in the scale-secreting epithelium. They also show sequence and fewer still in the symbiont-hosting oesophageal gland similarity with the well-known shell matrix protein pif , which (Fig. 2c). Although these novel genes do not correspond to other by contrast was highly expressed in the mantle in the Scaly-foot previously annotated genes, and were newly discovered within the Snail C. squamiferum, similar to other Mollusca, but not in the Scaly-foot Snail genome, they may yet be present in other taxa. scale-secreting epithelium. A similar pattern is seen in the chitin This is only the first genome in the subclass Neomphaliones, and synthase gene family, as the Scaly-foot Snail genome possesses the poor taxon coverage and limited completeness of existing diverse subgroups of this family, and different paralogues were molluscan genomes do not provide a reliable reference framework found to be highly expressed in the mantle or scale-secreting to infer whether these ‘novel’ genes are truly lineage-specific epithelium (Supplementary Fig. 6). These chitin synthase genes 6 NATURE COMMUNICATIONS | (2020)11:1657 | https://doi.org/10.1038/s41467-020-15522-3 | www.nature.com/naturecommunications TPM TPM NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15522-3 ARTICLE Octopus bimaculoides 28.91 Mizuhopecten yessoensis Azumapecten farreri 470.06 548.49 Crassostrea gigas 334.95 Pinctada fucata 431.95 Modiolus philippinarum 100.08 Bathymodiolus platifrons 529.62 385.70 Haliotis discus hannai 453.12 Chrysomallon squamiferum Lottia gigantea 613.35 483.94 92.18 Lanistes nyassanus Pomacea canaliculata 404.51 Aplysia californica 196.54 116.84 Radix auricularia Biomphalaria glabrata 509.50 Capitella teleta 567.27 Notospermus geniculatus 503.51 Lingula anatina Phoronis australis C O S D C P T J K Pg N 600 500 400 300 200 100 0 Ma Mantle Mantle Mantle Mantle Fig. 6 Genomic and transcriptomic comparisons across Lophotrochozoa and the biomineralisation toolkit. a Genome-based phylogeny of selected taxa showing the position of the Scaly-foot Snail among lophotrochozoans and divergence times among molluscan lineages. Error bars indicate 95% confidence levels. The Scaly-foot Snail is highlighted in red, Gastropoda in pink, Mollusca in blue, and Lophotrochozoa in orange. b Transcription factors shown to be involved in armature secretion in the Scaly-foot Snail, compared with conchiferan molluscs (bivalves, gastropods, cephalopods), aculiferan molluscs (chitons), brachiopod, and annelids. Transcription factors only verified for the Scaly-foot Snail are highlighted in red, those known to be shared across Mollusca in blue, and those shared across Lophotrochozoa in orange. The top row for each group (darker shading) shows records of significant expression in shell-secreting mantle and bottom row (lighter shading) shows expression for other hard parts such as scales. Ma, million years ago. Source data are provided as a Source Data file. are clearly involved in the biomineralisation process, and may second population from the iron-poor Solitaire vent field that have been co-opted in the evolution of the sclerites. naturally lack any iron sulfide mineral coating . In the Kairei Although some highly expressed genes were different in the population, with iron sulfide mineralisation, transmembrane two biomineralising tissues, positions on the genome of transporter activity was comparatively enriched in highly transcription factors involved in the biomineralization toolkit expressed genes (Supplementary Data 8), supporting the were found to be synchronous with expression patterns seen in hypothesis the iron sulfide precipitation is indeed mediated by the both mantle and scale-secreting epithelium. We propose that they gastropod through transportation and precipitation of sulfur control downstream novel and existing biomineralization genes species . Another gene, metal tolerance protein (MTP) 9, in different ways to fabricate different kinds of armour or skeletal exhibited over 27-fold increased expression level (Supplementary elements by upstream activation. Fig. 7) in the population with iron sulfide mineralisation com- pared with the one without. This gene is widely found in inver- tebrates, protists, and even plants . Although poorly studied in invertebrates, functional assays of MTP in plants revealed Iron sulfide biomineralisation. We compared gene expression potential pathways for enhanced tolerance of metal ions and levels in the scale-secreting epithelium and shell-secreting mantle maintaining intracellular homoeostasis . of the Scaly-foot Snail from the iron-rich Kairei vent field with a NATURE COMMUNICATIONS | (2020)11:1657 | https://doi.org/10.1038/s41467-020-15522-3 | www.nature.com/naturecommunications 7 Hox4 Hox1 En Post1 Hox5 Arx Brachyury ETS Dlx Gsc Msx Six3/6 Hox2 Antp Zic Evx Mox grh Hox3 Gbx Sox14 Lox4 Lox5 Pax3/7 Pax6 Mollusca Lophotrochozoa ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15522-3 The Scaly-foot Snail biomineralises iron sulfide nanoparticles 50 Gb of Illumina Novaseq reads with the paired-end mode and a read length of 150 bp. by allowing sulfur that it actively recruits or deposits in the scales to react with iron ions diffusing in from its highly iron-enriched Oxford Nanopore Technologies (ONT) library preparation. A total of 2–3 μg environment , and the MTP9 gene likely helps the snail tolerate HMW genomic DNA in 10 mM Tris-HCl (pH 8.0) were used for each library such high levels of iron in its surrounding environment. This preparation. All libraries were prepared using the Ligation Sequencing Kit 1D would allow the Scaly-foot Snail to gain finer-scale control of (SQK-LSK108, ONT, Oxford, UK). The standard protocols [1D gDNA selecting for nanomaterial-scale production, which is completed at much long reads (SQK-LSK108) Protocol] from Oxford Nanopore Technologies were modified and performed as follows. lower temperatures than can be currently controlled in laboratory 12 The optional DNA repair step was not performed. End repair and dA-tailing settings . Many other species of molluscs incorporate iron into were performed using the NEBNext Ultra II End-Repair/dA-tailing Module (NEB hard parts, particularly the radula . Molluscs have been high- E7546). A total volume of 120 μl reaction mixture included 14 μl Ultra II End-Prep lighted as a model for generating biogenic nanomaterials at low buffer, 6 μl Ultra II End-Prep enzyme mix, and 100 μl genomic DNA. The reaction temperatures , but the genomic tools to open this part of the mixture was incubated at 20 °C for 30 min and 65 °C for 20 min using a thermal cycler. Clean-up was performed using a 0.4× volume AMPure XP beads (48 μl), biomineralisation toolkit for other applications were previously incubated at room temperature with gentle stirring for 5 min, washed twice with unavailable. 200 μl fresh 70% ethanol, and briefly air-dried for 1 min to obtain the pellet. DNA was eluted by adding 31 μl of nuclease-free water (NFW), resuspending the beads, and incubating for 10 min at 37 °C. One microlitre of aliquot was quantified by Discussion Qubit to ensure that ≥1.5 μg of DNA were retained. By comparing genomic elements and tissue-specific gene Ligation was then performed by gently mixing 20 μl Adaptor Mix (AMX1D, expression patterns in the Scaly-foot Snail scale- and shell- SQK-LSK108, ONT), 50 μl NEB Blunt/TA Master Mix (NEB M0367), and 30 μl dA-tailed DNA while incubating at room temperature for 30 min. The adaptor- secreting tissues, as well as other biomineralising tissues in ligated reaction was cleaned up with a 0.6 × volume (60 μl) of AMPure XP beads, lophotrochozoans, we revealed an ancient biomineralisation incubated for 5 min at room temperature, and followed by resuspending the pellet toolkit comprising at least 25 transcription factors that contribute in 500 μl Adapter Bead Binding buffer (ABB, SQK-LSK108, ONT). The purified- to biomineralisation across all lophotrochozoan hard parts ligated DNA was eluted using 15 μl of Elution Buffer (ELB, SQK-LSK108, ONT), resuspending beads, and incubating for 10 min at 37 °C. One microlitre of aliquot investigated to date. Gene families in the Scaly-foot Snail genome was quantified by Qubit to ensure that ≥500 ng of DNA were retained. The aliquot have predominantly ancient origins, as seen in other lopho- of the adapted and tethered DNA (the pre-sequencing Mix) was used for loading trochozoans (Fig. 3b), but their distribution and duplications into MinION Flow Cell. across various lineages are nonsynchronous with phylogenetic positions (Fig. 3a), underlining rapid modifications. Although MinION sequencing. MinION sequencing was performed as per manufacturer’s novel genes do appear to play important roles in downstream guidelines using R9 flow cells (FLO-MIN106, ONT). Priming of the SpotON Flow production of the hard parts, the hard parts themselves arise by Cell was preformed following the standard protocol (1D gDNA selecting for long reads (SQK-LSK108) Protocol). Directly after priming, 75 μl of the prepared library deploying elements from the conserved biomineralisation toolkit. mixed with 12 μl of the pre-sequencing Mix (adapted and tethered DNA library), Comparison between sclerites from unrelated groups of molluscs 25.5 μl of Library Loading Bead (LLB, ONT), 35 μl of Running Buffer Fuel Mix —the Scaly-foot Snail, and aculiferan molluscs, or Cambrian (RBF, ONT), and 2.5 μl of NFW were loaded through the SportON sample port in fossils—may underline the longevity of these gene families. a dropwise fashion. MinION sequencing was operated with MinKNOW 1.3.23 and fastq files were base-called with Albacore v2.3.4 (ONT) with the default setting. The true power of the biomineralisation toolkit lies in the Reads <3 kb were discarded. capacity for dynamic combination of the components being switched on or off, expanded or reduced, and relocated within the Genome assembly. The Illumina reads were trimmed with Trimmomatic v0.33 . genome, which creates compounding changes in phenome. The They were further assembled by Platanus v1.2.4 using the following settings: –k31 extreme morphological disparity of mollusc biomineralisation –u 0.2 –t29 –s 10. The genome size was estimated to be 444.4 Mb using the 17-mer underpins their successful diversification. Similarly, the re-use, re- histogram generated (Supplementary Note 1). The histogram was also submitted to arrangement and re-deployment of conserved genomic elements GenomeScope (http://qb.cshl.edu/genomescope/); and the genome heterozygosity 37 was estimated to be 1.38%. This indicates that the C. squamiferum genome is explains both our challenges over more than 540 million years relatively compact in Mollusca and also highly heterozygous. The total size, N50, in obtaining genomic and phylogenetic resolution, and their and mean length of the assembled genome from Platanus were 469.0 Mb, 18.2 kb evolutionary success. and 1.8 kb, respectively, indicating high fragmentation. We then used a range of bioinformatics pipelines to assemble the genome with ONT reads, including the ONT-only approaches (i.e. [canu version 1.7 ], [canu Methods 40 version 1.7 + smartdenovo (https://github.com/ruanjue/smartdenovo) ], and Sample collection. Specimens of the Scaly-foot Snail Chrysomallon squamiferum [minimap2 version 2.17-r943-dirty + miniasm version 0.3-r179], etc.) and the used in the present study were collected by the manned submersible Shinkai 6500 Illumina + ONT hybrid approach (e.g. MaSuRCA version 3.2.6). The detailed on-board multiple deep-sea research expeditions of R/V Yokosuka. For genome codes and settings of each assembly pipeline are detailed in Supplementary Note 2. sequencing, a specimen collected from the Kairei hydrothermal vent field (25° Following comparison of assembly statistics of different pipelines 19.23ʹS, 70°02.42ʹE, 2415 m depth, cruise YK16-02E, Feb 2016) and immediately (Supplementary Note 2), the genome assembled from the canu + smartdenovo placed into −80 °C upon recovery was used. For gene expression analyses, speci- pipeline was the best one and therefore used for downstream analyses. The mens fixed in situ using RNA stabilising agent from both the Kairei field (25° assembly was firstly error corrected five times with the ONT fastq file by Racon 19.23ʹS, 70°02.42ʹE, 2415 m depth, cruise YK13-03, Mar 2013) and the Solitaire 41 42 v1.4 and sequentially corrected twice with Illumina reads using Pilon v1.13 . field (19°33.41ʹS, 65°50.89ʹE, 2606 m depth, cruise YK13-02, Feb 2013) were used. The resultant error-corrected assembled genome had a total size, N50 and mean Specimens for in situ hybridisation were collected at the same time from the size of 407.8 Mb, 1.91 Mb and 393.7 kb, respectively. Analysis of the genome Solitaire field, but was fixed in 4% paraformaldehyde (PFA) solution upon recovery assembly completeness with metazoan Benchmarking Universal Single-Copy and later transferred to 80% ethanol. 43 Orthologs (BUSCO) v3.0 revealed that the completeness of the genome is 96.6%, only 2.5% of the reported metazoan genes were missing, and confirmed the presence of single-copy metazoan genes. High molecular weight DNA extraction and quantification. The foot and mantle of a Kairei field Scaly-foot Snail (specimen code E02B1) were used for DNA extraction. High Molecular Weight (HMW) DNA was extracted using the Microbial sequence contamination removal. Genome binning was performed by MagAttract HMW DNA Kit (Qiagen, Hilden, Germany) according to the manu- using MetaBAT 2 with the assembled contigs as input to check the microbial facturer’s instructions. The HMW DNA was further purified and concentrated sequence contamination. The resulting output genomes were examined for com- using the Genomic DNA Clean & Concentrator (gDCC-10) kit (ZYMO Research, pleteness and potential contamination using CheckM v1.0.7 , based on the pre- Irvine, CA, USA). DNA quality was assessed by running 1 μl of the sample on a sence of particular marker gens. Open reading frames (ORFs) were predicted using BioDrop µLITE (BioDrop, Holliston, MA, USA) to ensure the purified of DNA Prodigal v2.6.3 and only ORFs with closed end were retained. Phylogenetic with the OD 260/280 of 1.8 and the OD 260/230 of 2.0–2.2. Concentration of DNA analysis was performed based on 16 bacterial single-copy marker genes (frr, rplB, was assessed using the dsDNA HS assay on a Qubit fluorometer v3.0 (Thermo rplC, rplD, rplE, rplF, rplN, rplS, rpmA, rpsB, rpsC, rpsE, rpsJ, rpsS, smpB and tsf) Fisher Scientific, Singapore). A total of 1 μg DNA was used to obtain approximate share by all the genomes included for analyses. Protein sequence of the 16 genes 8 NATURE COMMUNICATIONS | (2020)11:1657 | https://doi.org/10.1038/s41467-020-15522-3 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15522-3 ARTICLE 47 61 were aligned individually by ClustalW in MEGA 6 , and then linked together to BLASTp results were used to assign the OGs by OrthoMCL v2.0.9 pipeline . OGs construct a maximum likelihood tree. For HGT analysis, all predicted protein from selected lophotrochozoan taxa (Fig. 6a) were used for the phylogenomic sequences from the two symbionts were BLASTP searched against the NCBI NR analysis. Only single-copy genes in each OG and genes that can be found in at least database (https://www.ncbi.nlm.nih.gov/refseq/) using an e-value 1e−5 with 50 60% of taxa (i.e. at least 11 species) were retained for downstream phylogenomic best hits. The BLASTP output files were queried into MEGAN 6 for taxonomy analysis, resulting in 1375 OGs. Gene sequences within each OGs were aligned by analysis in a latent class analysis (LCA) model. Summary of removed microbial MUSCLE, and the bona fide alignments were kept after trimming by TrimAL . sequences can be found in Supplementary Note 2. These alignments were concatenated, and the total alignments which contained missing sequences included 435,071 distinct alignment patterns across 19 species. The phylogenetic tree was conducted by RAxML v8.2.4 with the partition Hi-C sequencing and genome scaffolding. Another individual of the Scaly-foot information of each orthologue gene and GTR + gamma model. The final tree file Snail from the same locality, the Kairei field, and the same collection lot (specimen was viewed by FigTree v1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/). All boot- code E02B2) was used for Hi-C library preparation. For Hi-C library preparation strap values were 100, indicating full support. from the Chrysomallon squamiferum specimen E02B2 from Kairei field, the fol- MCMCtree , part of Phylogenetic Analysis by Maximum Likelihood (PAML) lowing methods were used, modified from Lieberman-Aiden et al. . Approxi- v4.9 (http://abacus.gene.ucl.ac.uk/software/paml.html), was used to predict the mately 2 g wet weight of foot tissue stored in −80 °C freezer was thawed on ice and divergence time among molluscs, with nodes constrained by fossil records and resuspended with 37% formaldehyde in serum-free DMEM for animal chromatin geographic events . The following fossil records and geographic events were used cross-linked. Then, the suspended tissues were homogenized and incubated at to constrain the nodes in the MCMC tree: A hard max time-point of 150 Ma for L. room temperature (RT) for 5 min, and glycine was added to a final concentration nyassanus and P. canaliculata, which correspond to the split of South America and of 0.25 M. The solution was incubated at RT for another 5 min and transferred on Africa ; minimum = 168.6 Ma and soft maximum = 473.4 Ma for A. californica ice for 15 min. The cells were further lysed in a pre-chilled lysis buffer which and B. glabrata ; hard minimum bound = 390 Ma for Caenogastropoda and includes 10 mM NaCl, 0.2% IGEPAL CA-630 (Sigma-Aldrich), 1X protease Heterobranchia ; minimum = 470.2 Ma and soft maximum = 531.5 Ma for A. inhibitor in 10 mM PH = 8.0 Tris-HCL buffer. The chromatin digestion, labelling californica (or B. glabrata) and L. gigantea ; and minimum = 532 Ma and soft and ligation with biotin steps followed Lieberman-Aiden et al. . The protein and maximum = 549 Ma for the first appearance of molluscs ; hard minimum = biotinylated free-ends were removed, and DNA was purified and sequenced by 465.0 Ma for the first appearance of Pteriomorpha ; and minimum = 550.25 Ma Illumina Novaseq platform with the paired-end mode and a read length of 150 bp. and soft maximum = 636.1 Ma for the first appearance of Lophotrochozoan . The The Hi-C raw Illumina reads were trimmed with Trimmomatic v0.33 and best protein substitution model was LG + I + G and was employed to each site. then assessed by HiC-Pro v2.10.0 . Only the valid reads generated by HiC-Pro The burn-in, sample frequency, and number of samples were set as 10,000,000, were further processed by the Juicer 1.5.6 pipeline . Then, genomic scaffolding 1000, and 10,000, respectively. was conducted with the 3D de novo assembly pipeline using the default diploid parameters. Several manual corrections were done in Juicebox to ensure the scaffolds within the same pseudo-chromosomal linkage groups met the Hi-C Gene family analyses. The OGs deduced from OrthoMCL v2.0.9 were also used for linkage characteristics. In this process, three contigs were split into two parts and 69,70 the gene family expansion and contraction analysis by CAFÉ v3.1 pipeline .Only anchored to different chromosomes. In total, 1025 contigs were scaffolded into 15 those with gene family wide P value <0.01 and a taxon-specificViterbi P value <0.01 pseudo-chromosomal linkage groups, and only seven contigs were not anchored were considered as an event of expansion/contraction. due to insufficient Hi-C linkage found on them. The Hox gene cluster found in two To determine the gene age of the OGs, the OrthoMCL result was also used to contigs in the pre-Hi-C version was linked into an intact one, suggesting a good deduce the time of origin, after removing taxon-specific OGs. OGs with genes only performance of the Hi-C scaffolding method. The final genome assembly statistics found in gastropods were classified as Gastropoda-specific. The similar approach can be found in Supplementary Table 1. was applied to assign OGs as Mollusca-specific, Lophotrochozoa-specific, Bilaterian-specific and Eumetazoan-specific. A binary (present/absent) matrix was also deduced from the OrthoMCL results, and taxon-specific OGs were further Genome quality check and repeats identification. Quast v4.0 was used to check excluded. The binary matrix was used to plot a correlation matrix heatmap using genome assembly quality with ONT reads, for assembly assessment report see the corrplot library in R 3.5.2 . Supplementary Note 2. Repeats and transposable elements were annotated using In order to examine the chromosomal distribution of the expanded genes or the the RepeatModeler 2.0 and RepeatMasker 4.0.8 pipeline with the searching genes that are highly expressed in each tissue, a hypergeometric test was applied to programme of NCBI RMblast v2.9.0. The species-specific repeat library was assess whether there was any bias in any particular chromosome, with the annotated with RepeatModeler. Afterwards, RepeatMasker was run twice, one assumption that genes are randomly distributed in each chromosome. The using the species-specific repeat library and another one using the repeats in distribution in the expanded gene families is summarised in Supplementary Note 3. Repbase (https://www.girinst.org/repbase/). All the results were summarized and classified with the perl script buildSummary.pl in the RepeatMasker package, for summary table see Supplementary Note 2. Transcriptome sequencing. The Scaly-foot Snail used in genome sequencing (specimen code E02B1) was also dissected into a number of tissues/organs, namely digestive gland, scales, foot muscle, ctenidium (gill), mantle, nerve cord, nephri- Genome annotations. Two versions of transcriptome assemblies were generated: dium, endosymbiont-containing oesophageal gland, testis and ventricle, following a (1) the transcriptome reads (see Transcriptome sequencing section below) were de previously published account of the Scaly-foot Snail anatomy . This dataset was novo assembled by Trinity v2.6.5; (2) the transcriptome reads were aligned to the only used for gene prediction, as the sample was not fixed in situ. Meanwhile, for genome using histat2 v2.1.0 , and the aligned .bam file was assembled by Trinity comparative purposes, only in situ fixed samples were considered; two individuals under the genome-guided mode. These two versions of transcriptome were merged from the Kairei field and another three from the Solitaire field were used, targeting by PASA pipeline v2.2.0 with the aligners of BLAT and gmap and were further tissues of interest including scales, foot, gill, mantle and oesophageal gland. clustered with cd-hit-est v4.6 with a minimum sequence identity of 0.95 . RNA was extracted using TRIzol reagents and further sequenced using the Maker v3.0 was initially ran with the transcript evidence alone, and only gene Illumina Novaseq platform with the approximate output of 5 Gb for each tissue, models with an AED score <0.01 were retained. Gene models with <3 exons, with read length of 150 bp and paired-end mode. The raw reads were cleaned with incomplete open reading frame, and with an inter-genic region <3 kb were Trimmomatic v0.36 . Gene expression level in each tissue was quantified by removed. The rest of bona fide gene models were used to train Augustus v3.1, a de 58 Kallisto v0.44.0 with sequence based bias correction. Differentially expressed novo gene predictor . Then, the gene model prediction was performed by using genes were determined by DESeq2 using the normalization method of Loess, a Maker again, but with transcript evidence, protein evidence, Augustus gene minimum read count of 10, and paired test (n = 5). predictions, as well as an automatic annotation integration of these data into a Tissue-specific genes for the tissues were determined based on their expression consensus annotation based on their evidence-based weights. levels compared across all tissue types. To minimise the batch effect from pooling Gene functions were determined by searching the NCBI non-redundant different sample collection events and localities, a paired test method was applied. database (https://www.ncbi.nlm.nih.gov/refseq/) and SwissProt database via Only genes that were overexpressed with a fold change of over 2 and FDR < 0.05 UniProt (https://www.uniprot.org/) with the settings of: -evalue 1e-5 -word_size 3- against other tissue types were classified as highly expressed. Differentially num_alignments 20 -max_hsps 20. The Gene Ontology (GO) functional categories 59 expressed genes in the scales and mantle between the Kairei and Solitaire fields were deduced from the BLAST2GO pro v.4.1.9 software. Kyoto Encyclopedia of were also compared with shed light on the iron sulfide biomineralisation. Genes and Genomes (KEGG) annotation was performed using the KEGG Dominant functions of these target genes were further assessed with GO Automatic Annotation Server (https://www.genome.jp/kegg/kaas/) with the bi- 73 59 enrichment analyses using GOEAST v.1.30 or in the BLAST2GO v.4.1.9 directional best hit (BBH) method. A sensitive HMM scanning method on the package (see Supplementary Note 4 for results). known pfam functional domains with an e-value of 0.05 was also be used to classify the gene families. The low complexity protein was predicted with online XSTREAM v1.73 (https://amnewmanlab.stanford.edu/xstream/). Real-time PCR validation. Real-time PCR was employed to validate gene expression patterns in the two main tissues of interest, scales and mantle, as well as Phylogenetic analyses. The orthologue groups (OGs) were determined by a three selected other tissue types for comparison, including foot muscle, oesophageal BLASTp search against protein sequences of other available high-quality mollus- gland, and gill. The primers for real-time PCR were designed with the on-line can, lophotrochozoan, and metazoan genomes (see Supplementary Note 6). The NCBI Primer-BLAST tool (for a list of primers see Supplementary Note 5). PCR NATURE COMMUNICATIONS | (2020)11:1657 | https://doi.org/10.1038/s41467-020-15522-3 | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15522-3 underlying Figs. 1, 2c, 3b, 4a, 5c and 6b, and Supplementary Figs. 1, 2, 3, 6b, 7c, 8, 9, 10, product length was set within the range of 100–200 bp, and optimal melting 11, 12 and 13 are provided as a Source Data file. Although specimens of Chrysomallon temperature was set as 60.0 °C. Only Primer-pair with the least possibilities of self- squamiferum are available from the authors at a reasonable request, the numbers complementarity and self 3ʹ complementarity was selected for each gene. The available are very limited. elongation factor 1-alpha (EF1a) was selected as the internal standard gene, as its expression level remain almost constant across all the samples. Total RNA was extracted from each type of tissue by Trizol method and trace Code availability amount of DNA was removed with TURBO DNA-free kit (Thermo Fisher All codes, commands, and intermediate files for the bioinformatics analyses carried out Scientific), and the first strand cDNA was synthesized by using High Capacity (using freely or commercially available software, as listed in the Methods section) in the cDNA Reverse Transcription Kit (Applied Biosystems). Real-Time PCR was present study are contained in Supplementary Data 9. performed with SYBR Green RT-PCR Reagents Kit (Applied Biosystems) on LightCycler 480 II (Roche) with the following procedures: (1) polymerase activation at 95 °C for 10 min, and (2) annealing and extending at 57 °C for 1 min Received: 21 March 2019; Accepted: 13 March 2020; with a total of 40 cycles. The specificity of primer pairs for the PCR amplification was checked by the melting curve method. Triplicates were applied for each gene, ΔΔCt 74 and the relative gene expression level was calculated based on the 2 method . In situ hybridisation. To localize the expression of genes involved in scale and mantle formation, in situ hybridisation was performed on scale-secreting epidermis References and the mantle of Scaly-foot Snails collected from the Solitaire field, fixed in 4% 1. Vermeij, G. J. The origin of skeletons. Palaios 4, 585–589 (1989). PFA solution and stored in 80% ethanol. In situ hybridisation was performed 2. Marlétaz, F., Peijnenburg, K. T. C. A., Goto, T., Satoh, N. & Rokhsar, D. S. A according to methods detailed in Miyamoto et al. with slight modifications, new spiralian phylogeny places the enigmatic arrow worms among detailed as follows. gnathiferans. Curr. Biol. 29, 312–318.e313 (2019). Whole-mount in situ hybridisation was carried out for the mantle. Samples 3. Zhuravlev, A. Y. & Wood, R. A. The two phases of the Cambrian Explosion. were rehydrated and washed with PBST (i.e. PBS containing 0.1% Tween 20) three Sci. Rep. 8, 16656 (2018). times. The samples were digested with 10 µg/ml proteinase K/PBST for 30 min at 4. McDougall, C. & Degnan, B. M. The evolution of mollusc shells. Wiley 37 °C. After a brief wash with PBST, the samples were post-fixed in 4% PFA/PBST Interdiscip. Rev.: Developmental Biol. 7, e313 (2018). for 10 min RT (20–25 °C) and washed three times in PBST. Samples were 5. Kocot, K. M. et al. Phylogenomics of Lophotrochozoa with consideration of prehybridised in a prehybridisation solution (50% formamide, 5 × SSC, 5 × systematic error. Syst. Biol. 66, 256–282 (2016). Denhardt’s solution, 100 µg/ml yeast RNA, 0.1% Tween 20) at 60 °C for 4 h and 6. Chen, C., Copley, J. T., Linse, K., Rogers, A. D. & Sigwart, J. How the mollusc then hybridised with a hybridisation solution containing a digoxigenin (DIG)- got its scales: convergent evolution of the molluscan scleritome. Biol. J. Linn. labelled RNA probe at 60 °C, for 3 days. Hybridised samples were washed twice in a Soc. 114, 949–954 (2015). solution of 50% formamide, 4 × SSC, and 0.1% Tween 20 for 30 min each; then 7. Warén, A., Bengtson, S., Goffredi, S. K. & Van Dover, C. L. A hot-vent twice in 50% formamide, 2 × SSC, and 0.1% Tween 20 for 30 min twice; 2 × SSC, gastropod with iron sulfide dermal sclerites. Science 302, 1007–1007 and 0.1% Tween 20 for 30 min each time; and twice in 0.2 × SSC and 0.1% Tween (2003). 20 for 30 min each, at 60 °C. These were then washed with MABT (i.e. maleic acid 8. Vrijenhoek, R. C. On the instability and evolutionary age of deep-sea buffer containing 0.1% Tween 20) three times for 30 min at RT, blocked in 2% chemosynthetic communities. Deep Sea Res. Pt. II 92, 189–200 (2013). blocking reagent (Roche) in MABT for 2 h at RT, and incubated overnight with a 1:1500 dilution of antDIG-AP antibody (Roche) in the blocking buffer at 4 °C. 9. Kaim, A., Jenkins, R. G., Tanabe, K. & Kiel, S. Mollusks from late Mesozoic seep deposits, chiefly in California. Zootaxa 3861, 401–440 (2014). Samples were then further washed six times with MABT for 60 min each on a rocker and then transferred into TNT buffer (100 mM Tris pH 9.5, 100 mM NaCl, 10. Stöger, I. et al. The continuing debate on deep molluscan phylogeny: evidence 0.1% Tween 20). A chromogenic reaction was performed using nitro blue for Serialia (Mollusca, Monoplacophora + Polyplacophora). BioMed. Res. Int. tetrazolium chloride/5-bromo-4-chloro-3-indoly-phosphate (NBT/BCIP; Roche) in 2013, 18 (2013). AP buffer (100 mM Tris pH 9.5, 100 mM NaCl, 50 mM MgCl , 0.1% Tween 20 and 11. Nakamura, K. et al. Discovery of new hydrothermal activity and 2% polyvinyl alcohol) until signals were visible. The reaction was stopped in PBST, chemosynthetic fauna on the Central Indian Ridge at 18°–20°S. PLoS ONE 7, post-fixed in 10% formalin/PBST, rewashed with PBST, mounted with 40% e32965 (2012). glycerol, and observed under a light microscope (IX71, Olympus). 12. Okada, S. et al. The making of natural iron sulfide nanoparticles in a hot vent In situ hybridisation of the scale-secreting epidermis was carried out with snail. Proc. Natl. Acad. Sci. USA, https://doi.org/10.1073/pnas.1908533116 sections. Samples were washed three times with PBS and mounted in Tissue-Tek O. (2019). C.T. compound to use in frozen sectioning. Frozen sections were air-dried, washed 13. Goffredi, S. K., Waren, A., Orphan, V. J., Van Dover, C. L. & Vrijenhoek, R. C. with PBST, and fixed in 4% PFA/PBS for 10 min at RT. The slides were washed Novel forms of structural integration between microbes and a hydrothermal with PBS and digested with 1 µg/mL proteinase K in PBS for 10 min at RT. After a vent gastropod from the Indian Ocean. Appl Environ. Microbiol. 70, brief wash with PBS, the samples were post-fixed in 4% PFA/PBS for 10 min at RT. 3082–3090 (2004). The slides were washed with PBS three times. The samples were prehybridized for 14. Chen, C., Copley, J. T., Linse, K., Rogers, A. D. & Sigwart, J. D. The heart of a at least 1 h in prehybridisation solution (50% formamide, 5 × SSC, 5 × Denhardt’s dragon: 3D anatomical reconstruction of the ‘scaly-foot gastropod’ (Mollusca: solution, 50 µg/mL yeast RNA) at 60 °C and hybridized with a DIG-labelled RNA Gastropoda: Neomphalina) reveals its extraordinary circulatory system. Front probe at 60 °C for 3 days. The slides were washed with a solution of 50% Zool. 12, 13 (2015). formamide and 2 × SSC for 30 min; twice in 2 × SSC for 30 min each; 0.2 × SSC for 15. Sigwart, J. D. Molluscs all beneath the sun, one shell, two shells, more, or 30 min twice at 60 °C. They were further rinsed three times with MAB, blocked in none. Curr. Biol. 27, R708–R710 (2017). 2% blocking reagent/MAB for 2 h at RT and incubated overnight at 4 °C with a 16. Nakagawa, S. et al. Allying with armored snails: the complete genome of 1:1500 dilution of anti-DIG-AP antibody (Roche) in blocking buffer. Finally, they gammaproteobacterial endosymbiont. ISME J. 8,40–51 (2014). were washed six times with MAB for 30 min each and transferred into AP buffer 17. Jain, M. et al. Nanopore sequencing and assembly of a human genome with (100 mM Tris pH 9.5, 100 mM NaCl, 50 mM MgCl and 2% polyvinyl alcohol). A ultra-long reads. Nat. Biotechnol. 36, 338 (2018). chromogenic reaction was performed using NBT/BCIP in AP buffer until a signal 18. Aguilera, F., McDougall, C. & Degnan, B. M. Co-option and de novo gene was visible. The reaction was stopped in PBS, postfixed in 4% PFA/PBS, washed evolution underlie molluscan shell diversity. Mol. Biol. Evol. 34, 779–792 with PBS and mounted with 80% glycerol. The hybridised samples were observed (2017). under a light microscope (IX71, Olympus). 19. Simakov, O. et al. Insights into bilaterian evolution from three spiralian genomes. Nature 493, 526 (2012). Reporting summary. Further information on research design is available in the 20. Wang, S. et al. Scallop genome provides insights into evolution of bilaterian Nature Research Reporting Summary linked to this article. karyotype and development. Nat. Ecol. Evol. 1, 0120 (2017). 21. Wanninger, A. & Wollesen, T. The evolution of molluscs. Biol. Rev. 94, 102–115 (2019). Data availability 22. Schiemann, S. M. et al. Clustered brachiopod Hox genes are not expressed The Chrysomallon squamiferum genome that support the findings of this study have been collinearly and are associated with lophotrochozoan novelties. Proc. Natl deposited in the NCBI Sequence Read Archive under the BioProject number Acad. Sci. USA 114, E1913–E1922 (2017). PRJNA523462, all raw sequencing data, including Illumina and Nanopore reads, are also 23. Shinzato, C. et al. Using the Acropora digitifera genome to understand coral deposited under the same BioProject number. The assembled genome, transcriptome, 76 responses to environmental change. Nature 476, 320 (2011). predicted transcripts, proteins have been deposited in Dryad . Publicly available datasets 24. Mollenhauer, J. et al. DMBT1, a new member of the SRCR superfamily, on used in the study include the following: NCBI NR database (https://www.ncbi.nlm.nih. chromosome 10q25.3–26.1 is deleted in malignant brain tumours. Nat. Genet. gov/refseq/), Repbase (https://www.girinst.org/repbase/), SwissProt database via UniProt 17,32–39 (1997). (https://www.uniprot.org/), and KEGG (https://www.genome.jp/kegg/). The source data 10 NATURE COMMUNICATIONS | (2020)11:1657 | https://doi.org/10.1038/s41467-020-15522-3 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15522-3 ARTICLE 25. Nam, B.-H. et al. Genome sequence of pacific abalone (Haliotis discus hannai): 56. Li, W. & Godzik, A. Cd-hit: a fast program for clustering and comparing large the first draft genome in family Haliotidae. GigaScience 6, gix014 (2017). sets of protein or nucleotide sequences. Bioinformatics 22, 1658–1659 (2006). 26. Kocot, K. M., Aguilera, F., McDougall, C., Jackson, D. J. & Degnan, B. M. Sea 57. Cantarel, B. L. et al. MAKER: an easy-to-use annotation pipeline designed for shell diversity and rapidly evolving secretomes: insights into the evolution of emerging model organism genomes. Genome Res. 18, 188–196 (2008). biomineralization. Front Zool. 13, 23 (2016). 58. Stanke, M. & Morgenstern, B. AUGUSTUS: a web server for gene prediction 27. Chen, C., Uematsu, K., Linse, K. & Sigwart, J. D. By more ways than one: rapid in eukaryotes that allows user-defined constraints. Nucleic Acids Res. 33, convergence at hydrothermal vents shown by 3D anatomical reconstruction of W465–W467 (2005). Gigantopelta (Mollusca: Neomphalina). BMC Evol. Biol. 17, 62 (2017). 59. Götz, S. et al. High-throughput functional annotation and data mining with 28. Cope, J. & Babin, C. Diversification of bivalves in the Ordovician. Geobios 32, the Blast2GO suite. Nucleic Acids Re.s 36, 3420–3435 (2008). 175–185 (1999). 60. Kanehisa, M. & Goto, S. KEGG: Kyoto encyclopedia of genes and genomes. 29. Cunha, T. J. & Giribet, G. A congruent topology for deep gastropod Nucleic Acids Res. 28,27–30 (2000). relationships. Proc. R. Soc. B: Biol. Sci. 286, 20182776 (2019). 61. Li, L., Stoeckert, C. J. Jr. & Roos, D. S. OrthoMCL: identification of ortholog 30. Uribe, J. E., Irisarri, I., Templado, J. & Zardoya, R. New patellogastropod groups for eukaryotic genomes. Genome Res. 13, 2178–2189 (2003). mitogenomes help counteracting long-branch attraction in the deep 62. Capella-Gutierrez, S., Silla-Martinez, J. M. & Gabaldon, T. trimAl: a tool for phylogeny of gastropod mollusks. Mol. Phylogenet. Evol. 133,12–23 (2019). automated alignment trimming in large-scale phylogenetic analyses. 31. Ponder, W. F. & Lindberg, D. R. Towards a phylogeny of gastropod molluscs: Bioinformatics 25, 1972–1973 (2009). an analysis using morphological characters. Zool. J. Linn. Soc. 119,83–265 63. Stamatakis, A., Ludwig, T. & Meier, H. RAxML-III: a fast program for (1997). maximum likelihood-based inference of large phylogenetic trees. 32. Zakrzewski, A.-C. et al. Early divergence, broad distribution, and high Bioinformatics 21, 456–463 (2005). diversity of animal chitin synthases. Genome Biol. Evol. 6, 316–325 (2014). 64. dos Reis, M. & Yang, Z. Approximate likelihood calculation on a phylogeny for 33. Suzuki, M. et al. An acidic matrix protein, Pif, is a key macromolecule for Bayesian estimation of divergence times. Mol. Biol. Evol. 28, 2161–2172 (2011). nacre formation. Science 325, 1388–1390 (2009). 65. Benton, M. J. et al. Constraints on the timescale of animal evolutionary 34. Ricachenevsky, F., Menguer, P., Sperotto, R., Williams, L. & Fett, J. Roles of history. Palaeontol. Electron 18,1–106 (2015). plant metal tolerance proteins (MTP) in metal storage and potential use in 66. Hayes, K. A. et al. Molluscan models in evolutionary biology: apple snails biofortification strategies. Front. Plant Sci. 4, 144 (2013). (Gastropoda: Ampullariidae) as a system for addressing fundamental 35. Hilgers, L., Hofreiter, M., Hartmann, S. & von Rintelen, T. Novel genes, questions. Am. Malacol. Bull. 27,47–58 (2009). ancient genes, and gene co-option contributed to the genetic basis of the 67. Benton, M. J., Donoghue, P. C. J. & Asher, R. J. in The Timetree of Life (eds S. radula, a molluscan innovation. Mol. Biol. Evol. 35, 1638–1652 (2018). Blair Hedges, S. & Kumar, S.) 35–86 (Oxford University Press, 2009). 36. Grunenfelder, L. et al. Stress and damage mitigation from oriented 68. Jörger, K. M. et al. On the origin of Acochlidia and other enigmatic nanostructures within the radular teeth of Cryptochiton stelleri. Adv. Funct. euthyneuran gastropods, with implications for the systematics of Mater. 24, 6093–6104 (2014). Heterobranchia. BMC Evol. Biol. 10, 323 (2010). 37. Kocot, K. M., McDougall, C. & Degnan, B. M. in Physiology of Molluscs, A 69. Han, M. V., Thomas, G. W., Lugo-Martinez, J. & Hahn, M. W. Estimating Collection of Selected Reviews Vol. 1 (eds Saleuddin, S. & Mukai, S.) (Apple gene gain and loss rates in the presence of error in genome assembly and Academic Press, 2017). annotation using CAFE 3. Mol. Biol. Evol. 30, 1987–1997 (2013). 38. Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for 70. Sun, J. et al. Adaptation to deep-sea chemosynthetic environments as revealed Illumina sequence data. Bioinformatics 30, 2114–2120 (2014). by mussel genomes. Nat. Ecol. Evol. 1, 0121 (2017). 39. Koren, S. et al. Canu: scalable and accurate long-read assembly via adaptive k- 71. R Development Core Team. R: A Language and Environment Forstatistical mer weighting and repeat separation. Genome Res. 27, 722–736 (2017). Computing (R Foundation for Statistical Computing, 2008). 40. Schmidt, M. H.-W. et al. De Novo assembly of a new Solanum pennellii 72. Bray, N. L., Pimentel, H., Melsted, P. & Pachter, L. Near-optimal probabilistic accession using nanopore sequencing. Plant Cell 29, 2336–2348 (2017). RNA-seq quantification. Nat. Biotechnol. 34, 525 (2016). 41. Vaser, R., Sovic, I., Nagarajan, N. & Sikic, M. Fast and accurate de novo 73. Zheng, Q. & Wang, X. J. GOEAST: a web-based software toolkit for Gene genome assembly from long uncorrected reads. Genome Res. 27, 737–746 Ontology enrichment analysis. Nucleic Acids Res. 36, W358–W363 (2008). (2017). 74. Livak, K. J. & Schmittgen, T. D. Analysis of relative gene expression data using 42. Walker, B. J. et al. Pilon: An integrated tool for comprehensive microbial real-time quantitative PCR and the 2(-Delta Delta C(T)) method. Methods 25, variant detection and genome assembly improvement. PLoS ONE 9, e112963 402–408 (2001). (2014). 75. Miyamoto, N., Yoshida, M.-a, Koga, H. & Fujiwara, Y. Genetic mechanisms of 43. Simao, F. A., Waterhouse, R. M., Ioannidis, P., Kriventseva, E. V. & Zdobnov, bone digestion and nutrient absorption in the bone-eating worm Osedax E. M. BUSCO: assessing genome assembly and annotation completeness with japonicus inferred from transcriptome and gene expression analyses. BMC single-copy orthologs. Bioinformatics 31, 3210–3212 (2015). Evol. Biol. 17, 17 (2017). 44. Kang, D. D., Froula, J., Egan, R. & Wang, Z. MetaBAT, an efficient tool for 76. Sun, J. et al. Data from: The scaly-foot snail genome and implications for the accurately reconstructing single genomes from complex microbial ancient origins of biomineralised armour. Dryad, https://doi.org/10.5061/ communities. PeerJ 3, e1165 (2015). dryad.24053dn (2020). 45. Parks, D. H., Imelfort, M., Skennerton, C. T., Hugenholtz, P. & Tyson, G. W. CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res. 25, 1043–1055 (2015). Acknowledgements 46. Hyatt, D. et al. Prodigal: prokaryotic gene recognition and translation The authors thank the captain and crew of R/V Yokosuka during cruises YK13-02 initiation site identification. BMC Bioinforma. 11, 119 (2010). (principal scientist: Manabu Nishizawa), YK13-03 (principal scientist: Kentaro Naka- 47. Tamura, K., Stecher, G., Peterson, D., Filipski, A. & Kumar, S. MEGA6: mura) and YK16-02E (principal scientist: Ken Takai), and pilots of DSV Shinkai 6500. molecular evolutionary genetics analysis version 6.0. Mol. Biol. Evol. 30, Takuro Nunoura (JAMSTEC) is gratefully acknowledged for bringing together colla- 2725–2729 (2013). borative efforts. This research was financially supported by the China Ocean Mineral 48. Huson, D. H. et al. MEGAN community edition - interactive exploration and Resource Research and Development Association (DY135-E2-1-03 to P.-Y.Q.), the Hong analysis of large-scale microbiome sequencing data. PLoS Comput. Biol. 12, Kong Branch of South Marine Science and Engineering Guangdong Laboratory e1004957 (2016). (SMSEGL20Sc01 to P.-Y.Q.), the Research Grants Council of Hong Kong (GRF grant No. 49. Lieberman-Aiden, E. et al. Comprehensive mapping of long-range interactions 16101219 to J.S., C.C., and P.-Y.Q.), and a Japan Society for the Promotion of Science reveals folding principles of the human genome. Science 326, 289–293 (2009). Grant-in-Aid for Scientific Research (18K06401 to C.C. and H.K.W.). Illumina 50. Servant, N. et al. HiC-Pro: an optimized and flexible pipeline for Hi-C data sequencing was performed by Novogene (Beijing, China). Sampling in the Mauritian processing. Genome Biol. 16, 259 (2015). EEZ was approved by Ministry of Foreign Affairs, Regional Integration, and Interna- 51. Durand, N. C. et al. Juicer provides a one-click system for analyzing loop- tional Trade, Mauritian Government (Ref. 29/2014; 50/38/24 V2). resolution Hi-C experiments. Cell Syst. 3,95–98 (2016). 52. Dudchenko, O. et al. De novo assembly of the Aedes aegypti genome using Hi- C yields chromosome-length scaffolds. Science 356,92–95 (2017). Author contributions 53. Smit, A. F. A. & Hubley, R. RepeatModeler/RepeatModeler, http://www. P.-Y.Q., K.T., C.C., J.S., Y.T., J.-W.Q. and J.D.S. conceived and designed the project. C.C., repeatmasker.org/ (2008). T.W., H.K.W., K.T., J.D.S., K.I., N.F., K.Y., D.B. and N.M. collected and fixed the samples. 54. Kim, D., Langmead, B. & Salzberg, S. L. HISAT: a fast spliced aligner with low C.C., J.S., T.W. and N.M. dissected the samples. J.S., Y.S., T.X. and Y.L. extracted RNA, memory requirements. Nat. Methods 12, 357–360 (2015). DNA, and performed the ONT sequencing. W.Z. checked bacteria contamination. N.M. 55. Haas, B. J. et al. Improving the Arabidopsis genome annotation using performed in situ hybridisation. R.L. and J.S. performed synteny analysis. J.S. and J.C.H.I. maximal transcript alignment assemblies. Nucleic Acids Res. 31, 5654–5666 performed gene family analyses. J.S. assembled the genome, annotated the genes, and (2003). performed bioinformatics analyses. W.C.W. performed the real-time PCR. C.C., J.D.S., NATURE COMMUNICATIONS | (2020)11:1657 | https://doi.org/10.1038/s41467-020-15522-3 | www.nature.com/naturecommunications 11 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15522-3 J.S. and N.M. interpreted the data and drafted the paper. All authors contributed to the Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in final paper. published maps and institutional affiliations. Competing interests Open Access This article is licensed under a Creative Commons The authors declare no competing interests. Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give Additional information appropriate credit to the original author(s) and the source, provide a link to the Creative Supplementary information is available for this paper at https://doi.org/10.1038/s41467- Commons license, and indicate if changes were made. The images or other third party 020-15522-3. material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the Correspondence and requests for materials should be addressed to K.T. or P.-Y.Q. article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from Peer review information Nature Communications thanks Bernard Degnan, Kevin Kocot the copyright holder. To view a copy of this license, visit http://creativecommons.org/ and David Lindberg for their contribution to the peer review of this work. Peer reviewer licenses/by/4.0/. reports are available. Reprints and permission information is available at http://www.nature.com/reprints © The Author(s) 2020 12 NATURE COMMUNICATIONS | (2020)11:1657 | https://doi.org/10.1038/s41467-020-15522-3 | www.nature.com/naturecommunications http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Nature Communications Springer Journals

Loading next page...
 
/lp/springer-journals/the-scaly-foot-snail-genome-and-implications-for-the-origins-of-HdUnMQjiH4

References (85)

Publisher
Springer Journals
Copyright
Copyright © The Author(s) 2020
eISSN
2041-1723
DOI
10.1038/s41467-020-15522-3
Publisher site
See Article on Publisher Site

Abstract

ARTICLE https://doi.org/10.1038/s41467-020-15522-3 OPEN The Scaly-foot Snail genome and implications for the origins of biomineralised armour 1,10 2,10 2,10 3 4,5 3 Jin Sun , Chong Chen , Norio Miyamoto , Runsheng Li , Julia D. Sigwart , Ting Xu , 1 1 3 1 1 6 Yanan Sun , Wai Chuen Wong , Jack C. H. Ip , Weipeng Zhang , Yi Lan , Dass Bissessur , 2,7 2 2 8 8 Tomo-o Watsuji , Hiromi Kayama Watanabe , Yoshihiro Takaki , Kazuho Ikeo , Nobuyuki Fujii , ✉ ✉ 9 3 2 1 Kazutoshi Yoshitake , Jian-Wen Qiu , Ken Takai & Pei-Yuan Qian The Scaly-foot Snail, Chrysomallon squamiferum, presents a combination of biomineralised features, reminiscent of enigmatic early fossil taxa with complex shells and sclerites such as sachtids, but in a recently-diverged living species which even has iron-infused hard parts. Thus the Scaly-foot Snail is an ideal model to study the genomic mechanisms underlying the evolutionary diversification of biomineralised armour. Here, we present a high-quality whole- genome assembly and tissue-specific transcriptomic data, and show that scale and shell formation in the Scaly-foot Snail employ independent subsets of 25 highly-expressed tran- scription factors. Comparisons with other lophotrochozoan genomes imply that this biomi- neralisation toolkit is ancient, though expression patterns differ across major lineages. We suggest that the ability of lophotrochozoan lineages to generate a wide range of hard parts, exemplified by the remarkable morphological disparity in Mollusca, draws on a capacity for dynamic modification of the expression and positioning of toolkit elements across the genome. Department of Ocean Science, Division of Life Science and Hong Kong Branch of the Southern Marine Science and Engineering Guangdong Laboratory (Guanzhou), The Hong Kong University of Science and Technology, Hong Kong, China. X-STAR, Japan Agency for Marine-Earth Science and Technology (JAMSTEC), 2-15 Natsushima-cho, Yokosuka, Kanagawa 237-0061, Japan. Department of Biology, Hong Kong Baptist University, Hong Kong, China. 4 5 6 Marine Laboratory, Queen’s University Belfast, Portaferry, N. Ireland. Senckenberg Museum, Frankfurt, Germany. Department for Continental Shelf, Maritime Zones Administration & Exploration, Ministry of Defence and Rodrigues, 2nd Floor, Belmont House, 12 Intendance Street, Port-Louis 11328, 7 8 Mauritius. Department of Food and Nutrition, Higashi-Chikushi Junior College, Kitakyusyu, Japan. National Institute of Genetics, 1111 Yata, Mishima, Shizuoka, Japan. Graduate School of Agricultural and Life Sciences, The University of Tokyo, 1-1-1, Yayoi, Bunkyo-ku, Tokyo, Japan. These authors contributed equally: Jin Sun, Chong Chen, Norio Miyamoto. email: kent@jamstec.go.jp; boqianpy@ust.hk NATURE COMMUNICATIONS | (2020)11:1657 | https://doi.org/10.1038/s41467-020-15522-3 | www.nature.com/naturecommunications 1 1234567890():,; 10 MB 10 MB 20 MB 10 MB 10 MB 30 MB 10 MB 20 MB 10 MB 40 MB 30 MB 30 MB 20 MB 20 MB 10 MB 10 MB ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15522-3 he appearance of biomineralised skeletons in the Cambrian Living hydrothermal vent taxa represent Cenozoic radiations 9,10 precipitated an evolutionary arms race and the original (≤66 million years ago, Mya), including peltospirid gastropods Texplosive diversification of modern animal forms . Since such as the Scaly-foot Snail; the Scaly-foot Snail represents a then, a long history of body plan modifications, reversals, and recent evolution of a complex scleritome. On its discovery in the convergences has obscured the relationships among some animal iron-rich Kairei hydrothermal vent field in the Indian Ocean in lineages . Understanding the genomic toolkit that enabled inno- 2001, the complex armature of the Scaly-foot Snail was compared vations in skeletons and armour is critical to reconstructing the with apparently plesiomorphic scleritomes found in aculiferan 3,4 early radiation of major clades . In particular, Lophotrochozoa molluscs and important early fossil taxa such as Halkieria and (including annelids, brachiopods, molluscs and others) presents a other sachtids . Later, a second population lacking iron sulfide significant and unresolved problem for phylogenetic reconstruc- mineralisation was discovered in the Solitaire hydrothermal site tion, and controversies over the affinities of key fossils . The characterised by low concentrations of iron (1/58 of that found at Scaly-foot Snail, Chrysomallon squamiferum, is unique among Kairei ), which presents an opportunity to understand physio- gastropod molluscs in having dense, imbricating chitinous logical and molecular mechanisms behind its unique miner- sclerites covering the whole distal surface of the soft foot , alisation pattern (also see Supplementary Note 1). Recent results forming a dermal scale armour in addition to a solid calcium indicated that the Scaly-foot Snail mediates its scale biominer- carbonate coiled shell (Fig. 1). These hard parts, including alisation by supplying sulfur through nano-scale channel-like sclerites and the shell, are often mineralised with iron sulfide, columns in the scales, which reacts with iron ions diffusing making it the only known metazoan using iron as a significant inwards from the vent fluid to make iron sulfide nanoparticles . component of skeleton construction . Many sachtids also had tissue-filled canals within their sclerites, Total Mantle (shell) Scales Novel genes Fig. 1 Key features of the Scaly-foot Snail Chrysomallon squamiferum genome. Circos plot showing the 15 pseudo-chromosomal linkage groups; with the Scaly-foot Snail at centre. The outer ring (red peaks) indicates gene density in each pseudo-chromosome, and the inner rings shows the normalized density of the highly expressed genes in the shell-secreting mantle (black peaks) and scale-secreting epithelium (white peaks, semi-transparent overlaid on black mantle peaks) as well as the density of novel genes (grey). The expression level is the average fold change of the target tissue versus the other four types of tissues (n = 5); the sliding window size is 10 kb. MB, megabases. Source data are provided as a Source Data file. 2 NATURE COMMUNICATIONS | (2020)11:1657 | https://doi.org/10.1038/s41467-020-15522-3 | www.nature.com/naturecommunications 0.6 0.4 0.2 0.6 0.4 0.2 0.8 0.6 0.4 0.2 30 MB 30 MB 20 MB 20 MB 10 MB 10 MB 20 MB 40 MB 10 MB 30 MB 20 MB 20 MB 10 MB 10 MB 10 MB 10 MB 0 NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15522-3 ARTICLE ab c Scaly-foot Snail Other Mollusca 100% Lophotrochozoa 80% 60% 4 40% 20% 85% 90% 95% 100% Scales Mantle Oes. Gland Gill 85% 90% 95% 100% BUSCO BUSCO Annotated Novel genes (Genome completeness) (Genome completeness) Fig. 2 Summary of the Scaly-foot Snail genome. a, b Quality comparisons between the Scaly-foot Snail genome and other available lophotrochozoan genomes: a Contig N50 vs BUSCO (Benchmarking Universal Single-Copy Orthologs); b Number of scaffolds vs BUSCO (for full comparison see Supplementary Table 2). Red dot indicates the Scaly-foot Snail genome, blue dots indicate other molluscan genomes, orange dots indicate other non- molluscan lophotrochozoan genomes. c Proportion of annotated genes (pink) and novel genes (grey) in four key tissue types of the Scaly-foot Snail, including scale-secreting epithelium, shell-secreting mantle, oesophageal gland (oes. gland), and gill. Source data are provided as a Source Data file. although entirely different in proportion and arrangement. The The Scaly-foot Snail genome represents the most complete channels within the scales of the Scaly-foot Snail are likely linked and continuous genome among the assembled mollusc an or to another key adaptation, in that it hosts sulfur-oxidising bac- lophotrochozoan genomes to date (Fig. 2), with a metazoan teria within cells of a highly vascularised, hypertrophied oeso- BUSCO (Benchmarking Universal Single-Copy Orthologs) score 13,14 phageal gland , and the sulfur may originate as metabolites of 96.6% for the genome assembly and 87.5% for the predicted from the endosymbionts . Taking these observations together, it transcripts. This genome assembly represents a successful is unclear whether the evolution of the sclerites in the Scaly-foot application of methods originally developed for model organ- Snail should be interpreted as a recurring ancestral phenome, or a isms to a little-studied taxon, providing a benchmark in quality recently derived adaptive novelty. compared with other published lophotrochozoan genomes The search for a biomineralisation toolkit underlying hard part confounded by heterozygosity and high repeat contents (Fig. 2a, evolution in molluscs and lophotrochozoans has previously b; Supplementary Table 2). focused on molluscan mantle gene expression and shell forma- tion . Molluscs have repeatedly ‘invented’ additional sclerite-like hard structures, in chitons, aplacophorans, Chrysomallon, and Gene family analyses. The Scaly-foot Snail has fewer novel gene also in other gastropods, cephalopods, and bivalves . Here, we families than other lophotrochozoans with high-quality genomes use a complete genome assembly of the Scaly-foot Snail and available (such as Pomacea canaliculata, Mizuhopecten yessoensis, tissue-specific transcriptomics to compare with data from other Lingula anatina, Achatina fulica, Sinonovacula constricta and lophotrochozoans in order to test whether there is indeed a Capitella teleta: Fig. 3a; Supplementary Table 2), which appears to universal biomineralisation toolkit in Mollusca or Lopho- conflict with an expectation that evolutionary novelties are trochozoa, and to identify specific genomic tools that enable the associated with novel genes (Fig. 3a). To corroborate this, only Scaly-foot Snail to repeatedly modify and duplicate hard parts. 11% of the Scaly-foot Snail gene families are not found among four other lophotrochozoan genomes, while in the other taxa, which apparently lack the dramatic morphological novelties of Results the Scaly-foot Snail, this figure is over 20% and up to 35% Genome assembly. The genome of a single specimen of Chry- (Fig. 3a). In addition, only 4.8% of gene families in the Scaly-foot somallon squamiferum (Fig. 1, Supplementary Fig. 1, Supple- Snail genome were found to be novel to gastropods, while as mentary Table 1) collected from the Kairei hydrothermal vent much as 87.0% of them may be unique to Lophotrochozoa or field was sequenced with a combination of Oxford Nanopore even pre-date the origin of Lophotrochozoa (Fig. 3b). Technologies and Illumina platforms (Supplementary Note 2; Comparative analyses among available lophotrochozoan gen- Supplementary Table 2). Sulfur-oxidising endosymbionts within omes (n = 15; Supplementary Table 4) showed that 351 gene the body are a source for significant potential contamination families were significantly expanded in the Scaly-foot Snail. A GO and therefore bacterial endosymbiont genomes were removed enrichment analysis on the expanded gene families revealed 30 from downstream analyses (Supplementary Note 2). The genome overrepresented GO categories in the Scaly-foot Snail genome of the Scaly-foot Snail is relatively compact for Mollusca (444.4 (Supplementary Fig. 2, Supplementary Data 1). These expanded Mb; Supplementary Table 1) and is highly heterozygous (1.38%) gene families appear to be involved in the secretion process of but with relatively low repeat content (25.2%). With additional proteinaceous materials (e.g. scavenger receptor activity, carbo- Hi-C data using a second specimen from the same population, hydrate binding, chitin-related metabolic process and chitin 1032 contigs (N50 = 1.88 Mb) were anchored to 15 pseudo- binding) and symbiosis (regulation of innate immune response chromosomal linkage groups (Fig. 1). A total of 16,917 gene and symbiont-containing vacuoles), suggesting their contribu- models were predicted (85.7% comparatively annotated; 2415 tions to both biomineralisation and regulation of endosymbiosis. additional novel genes) with evidence from transcripts, homo- The expanded genes were biased in distribution across the logue proteins, and ab initio methods (Supplementary Note 2). chromosomes, with Chr11 and Chr12 being especially enriched The number of genes in the Scaly-foot Snail genome appears to be (Fold change > 2 and FDR < 0.01, Supplementary Note 3). low, and the fact that 97.3% of the de novo assembled tran- We recovered the complete Antennapedia (ANTP) Hox gene scriptome could be mapped to the genome indicates that the complex containing 11 Hox genes on Chr11 (Fig. 4a) as in Lottia 19 20 genome assembly has high completeness. gigantea (a gastropod) and Mizuhopecten yessoensis (a NATURE COMMUNICATIONS | (2020)11:1657 | https://doi.org/10.1038/s41467-020-15522-3 | www.nature.com/naturecommunications 3 Contig N50 Scaffolds s 10 K 12 K 8 K 10 K 6 K 8 K 4 K 6 K 2 K 4 K 0 K 2 K 8 K 0 K 16 K 6 K ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15522-3 a b Eumetazoa Bilateria Lophotrochozoa Mollusca Gastropoda (11%) 76 73 202 45 2708 85 (26%) (21%) 106 127 380 277 2614 364 (35%) (20%) Fig. 3 Gene family analyses of lophotrochozoan genomes. a Venn diagram showing the number of shared and unique gene families among five lophotrochozoan genomes. b Circos plot showing the proportion and origin of shared gene families across five lophotrochozoan taxa. Arc values correspond to the number of gene families. Genes shared across Eumetazoa are indicated by grey lines, Bilateria by green lines, Lophotrochozoa by orange lines, Mollusca by blue lines, and Gastropoda by purple lines. K, kilo (1000). Source data are provided as a Source Data file. Mantle Scales Gill OG DMBT1 (scavenger domain) x 65 Gbx Evx Hox1 Hox2 Hox3 Hox4 Hox5 Lox5 Antp Lox4 Lox2 Post2 Post1 En En En Mox Mox Mnx DLx Chr11, 17.5M Gbx Evx Hox1 Hox2 Hox3 Hox4 Hox5 Lox5 Antp Lox4 Lox2 Post2 Post1 En En En Mox Mox Mnx DLx Mantle Scales Gill OG Hox4 Hox1 Evx Hox3 Fig. 4 Tissue-specific gene expression in the Scaly-foot Snail. a Arrangement of DMBT1 gene tandem repeats (grey), Hox genes and Hox-like genes (red) on the pseudo-chromosome Chr11 and their expression patterns in four types of tissues including mantle, scale-secreting epithelium, gill, oesophageal gland (OG); darker shades indicate higher levels of expression. SRCR (scavenger receptor cysteine-rich) domain, scavenger receptor cysteine-rich domain. Source data are provided as a Source Data file. b Tissue-specific expression patterns of transcriptome factors shown with in situ hybridisation, including Hox4 and Hox1 in the mantle edge (arrowheads; additional tissue is visible in the background including scales also expressing Hox1) and tissue sections of the scale-secreting epithelium expressing Evx and Hox3. Scale bars = 500 μm for mantle and 100 μm for scales. In situ hybridisation experiments were repeated independently on three individuals (twice on each individual for Hox1 and Hox3, once for Hox4 and Evx) with similar results. 4 NATURE COMMUNICATIONS | (2020)11:1657 | https://doi.org/10.1038/s41467-020-15522-3 | www.nature.com/naturecommunications 14 K 4 K 12 K 2 K 0 K 10 K 10 K 8 K 8 K 6 K 6 K 4 K 4 K 2 K 2 K 0 K 0 K –1 0 1 2 3 Z score P NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15522-3 ARTICLE bivalve); among molluscs, the intact ordered cluster has only been ‘aculiferan’ molluscs (i.e. chitons and the worm-like Caudofoveata recovered in one gastropod (Lottia), and two bivalve (Mizuho- and Solenogastres). The Scaly-foot Snail possesses an expanded 21 4 pecten, and Azumapecten farreri) genomes , but are critical to suite of mineral-secreting tissues , in the foot epidermis as well as understanding body plan evolution . The derived Scaly-foot the shell-forming mantle tissue. Many distinctive morphological Snail homeobox genes show reorganisation by intra- adaptations in Chrysomallon are connected to supporting endo- chromosomal rearrangement (Supplementary Figs. 3, 4). Syntenic symbiotic microbes housed in an enlarged oesophageal gland, a comparisons reveal both extensive inter- and intra-chromosomal condition resulting from rapid evolution . Indeed, even the rearrangements between the Scaly-foot Snail and other molluscs scales of the Scaly-foot Snail are apparently related to supporting (Supplementary Fig. 3). its endosymbionts . Gene expression levels in different tissue types of the Scaly-foot To assess tissue-specific functions, particularly in relation to Snail were assessed by transcriptome sequencing (five specimen biomineralisation, a paired test in DESeq2 was applied to replicates) using samples fixed in situ in the snail’s native deep- distinguish tissue highly expressed genes (n = 5, Figs. 1, 4a, and sea vent habitat to avoid the confounding effect of environment Supplementary Data 5–7), further confirmed with real-time PCR and hydrostatic pressure changes on transcription. The five (two replicates, five tissue types, Supplementary Note 5) and primary tissue types targeted included the shell-secreting mantle, in situ hybridisation (ISH; for mantle and scale-secreting scale-secreting epithelium, oesophageal gland where the bacter- epithelium; Figs. 4b, 5 and Supplementary Fig. 6). Given the iocytes housing the endosymbionts are located, the gill, and the high quality of the Scaly-foot Snail genome, genes related to foot musculature. biomineralisation can be accurately linked to their organisation in Among the expanded genes on Chr11 (262 genes), one that the genome and their respective expression levels in relevant controls protein binding to proteoglycan (for instance, chitin) tissues, i.e. shell-secreting mantle and scale-secreting epithelium. known as DMBT1 with scavenger receptor cysteine-rich (SRCR) Visualisation of expression levels by genome organisation domains was found to be highly expressed in both the scale- revealed similar expression patterns for scale-secreting epithelium secreting epithelium and shell-secreting mantle (Fig. 4a) with up to and mantle at chromosomal level across the genome, also 65 tandemly duplicated paralogues compared with one or two revealing a pattern where the pseudo-chromosome Chr11, which copies in other molluscs . Proteins with SRCR domains are contains the abovementioned ANTP Hox cluster as well as the commonly detected in shell matrix proteins . These gene copies DMBT genes, contained a high density of highly expressed genes were highly expressed in the mantle, but were also especially highly in these two biomineralising tissues (Fig. 1). The expression expressed in the scale-secreting epithelium (Fig. 4a). The Scaly-foot profiles of biomineralisation genes in the scale-secreting epithe- Snail DMBT1 gene copies had an average of 2.3 SRCR domains, lium was most similar to those of the shell-secreting mantle tissue but the number varied greatly among individual paralogues (0–5; (81 shared highly expressed genes, 11.2% of the scale-secreting Supplementary Data 2). As SRCR domains are known to evolve epithelium or 6.8% of mantle highly expressed genes); both rapidly to generate high gene diversity, the variation of the domain expressed common genes involved in biomineralisation (Supple- numbers in each paralogue likely reflects rapid evolution and mentary Data 5, 6), in separate tissues with different and complex expansion of this gene to coordinate protein secretion of dermal functions. scales and the shell periostracum. A gene tree of the Scaly-foot We identified 25 highly expressed transcription factors in the Snail DMBT1 shows that the paralogues segregate into two clades, scale-secreting epithelium (12) and mantle (13), often confirmed with only one copy from each clade being present in the by ISH (Figs. 4b, 6 and Supplementary Data 5, 6); the two vetigastropod Haliotis genome (Supplementary Fig. 5). biomineralisation tissues did not share any highly expressed Molluscan shell matrix proteins are commonly rich in transcription factors (Fig. 3c). Several of these were only detected repetitive low-complexity domains (RLCDs) , which have been as highly expressed in the Scaly-foot Snail hard parts among the attributed to a number of functions such as chitin or calcium available lophotrochozoan tissue-specific data (Fig. 3). However, binding and structural support, but such domains are lacking in all transcription factors reported to date from other lophotro- DMBT1 paralogues in the Scaly-foot Snail (Supplementary chozoan skeletal elements (brachiopod mantle and chaetae, and Data 3). Analyses of RLCD-rich genes across the Scaly-foot Snail annelid chaetae, n = 10 genes) were found to be also highly genome revealed that the ratio of RLCDs in genes highly expressed in the Scaly-foot Snail, in either the mantle (n = 6 expressed in the scale-secreting epithelium is similar to the overall genes) or scale-secreting epithelium (n = 4). Some of these co- ratio of the whole genome, whereas a higher ratio is seen in the opted transcription factors, and others, have also been reported mantle. This indicates that RLCDs likely contribute to shell from other molluscs including aculiferan mantle and sclerite- biomineralisation like known for other molluscs, but not the secreting tissue, gastropod or bivalve mantle tissue, or tissues 21,22 scales; suggesting that RLCDs may be linked with calcium secreting the radula, beak, or operculum (Fig. 6b and carbonate biomineralisation which is lacking in the scales Supplementary Table 4). The involvement of these transcription (Supplementary Table 3). factors in the formation of various hard parts across different A further 79 gene families were found to be significantly lophotrochozoan groups ranging widely from plesiomorphic (e.g. contracted in the Scaly-foot Snail genome relative to other brachiopod shell) to more recently evolved armature (e.g. Scaly- lophotrochozoans, including several enzymes that play critical foot Snail scales) signify that they together comprise the roles in steroid and amino acid synthesis, in keeping with similar biomineralisation toolkit preserved from the ancestral lopho- gene reductions in other chemosymbiont-dependent holobionts trochozoan genome. (Supplementary Data 4). Numerous enriched GO terms (Supple- Taxon coverage of high quality lophotrochozoan genomes is mentary Note 4) in highly-expressed genes in the oesophageal currently too limited to make significant phylogenetic inferences. gland are relevant to symbiosis, including oxidation-reduction Internal nodes within our time calibrated phlylogenetic recon- process compensating the reactive oxygen species generated by struction confirm and refine some previous divergence estimates. the endosymbiont . For example, Pteriomorphia, the best-sampled lineage within Mollusca, is recovered with a divergence estimate around 470 Mya (483.5–465.1 Mya) which corresponds to the early fossil Biomineralization toolkit. Until the early 2000s, the only living record for the group . Dating of earlier nodes (origins of molluscs known to regularly possess imbricating sclerites were Gastropoda, Bivalvia and Lophotrochozoa), and the topology, are NATURE COMMUNICATIONS | (2020)11:1657 | https://doi.org/10.1038/s41467-020-15522-3 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15522-3 a bc Azumapecten farreri Chitin binding Peritrophin-A 32,768 Crassotrea gigas Pinctada fucata Lingula anatina Pomacea canaliculata 100 Chrysomallon squamiferum CsqKR_Scaf21_18.24 Chrysomallon squamiferum CsqKR_Scaf21_18.25 Haliotis discus hannai Aplysia californica Scale Mantle Gill Foot OG Radix auricularia Tissue Radix auricularia Aplysia californica Pif protein Pomacea canaliculata de Lanistes nyassanus Chrysomallon squamiferum CsqKR_Scaf193_7.1 7 Lottia gigantea Bathymodiolus platifrons 94 Modiolus philippinarum Crassotrea gigas 100 66 Pinctada fucata Azumapecten farreri Mizuhopecten yessoensis 0.03 Modiolus philippinarum Scale Mantle Gill Foot OG Chitin binding Peritrophin-A Octopus bimaculoides von Willebrand factor typeA Tissue 0.3 Fig. 5 Characterisation of chitin binding peritrophin-A gene and pif gene in Chrysomallon squamiferum.a Phylogenetic analysis of chitin binding peritrophin-A and pif genes with names in red indicating genes from the Scaly-foot Snail, and domain architecture of the specific genes displayed on the right of names (purple indicates chitin binding peritrophin-A and blue indicates von Willebrand factor type A). b In situ hybridization of chitin binding Peritrophin-A in scales. c Boxplot showing the expression level of chitin binding peritrophin-A in five different tissue types including mantle, scale-secreting epithelium, gill, oesophageal gland (OG), and foot. The box depicts first to third quartile with whiskers indicating maximum and minimum expression levels, and the centre line refers to the median. Source data are provided as a Source Data file. d In situ hybridization showing the expression of pif gene in mantle. e Boxplot showing the expression level of pif in five different tissues types. Within the Boxplot, the centre line refers to the median, and the boxplot depicts the first to the third quartile, with the whiskers indicating maximum and minimum expression levels. likely confounded by limitations of the available data. Within genes (i.e, Scaly-foot Snail synapomorphies), or shared with, but gastropods, we recovered Neomphaliones (represented by as yet undiscovered in, other related taxa. Chrysomallon) sister to Vetigastropoda, and this clade sister to Genes known to produce proteins involved in molluscan shell Patellogastropoda (Fig. 6a). One recent phylotranscriptomics formation, such as pif, chitin-binding peritrophin-A domain 29 32 study similarly found Vetigastropoda sister to Patellogastro- , were also conserved and highly gene, and chitin synthase poda, but that analysis did not include any examples of expressed in the shell-secreting mantle and scale-secreting Neomphaliones. Neither our genome tree nor any phylotran- epithelium of the Scaly-foot Snail (Fig. 5 and Supplementary scriptomic analyses to date have included all gastropod subclasses Fig. 6). A GO enrichment analysis on highly expressed genes in as no whole genome is available for Neritimorpha . Early multi- the mantle and scale-secreting epithelium revealed that the gene phylogenetic analyses repeatedly encountered severe long categories integral component of membrane, scavenger receptor branch attraction artefacts particularly from patellogastropod activity, and cell-matrix adhesion were enriched in both of these sequences . A recent mitogenome analyses including all biomineralisation tissues (Supplementary Note 4), further gastropod subclasses found Patellogastropoda as the earliest- indicating that they elements of a shared biomineralisation 10,30 branching lineage within living gastropods , which is the toolkit. topology predicted by morphology . The debate on internal These possible downstream biomineralisation genes may be relationships among gastropods is likely to continue until more controlled by the transcription factors in the biomineralisation representatives from all subclasses become available for robust toolkit, and in a number of cases their specific expression in genome-based phylogeny. the scale-secreting epithelium were confirmed by ISH (Fig. 4b). Highly expressed ‘novel’ or taxon-specific genes were dis- The gene exhibiting the highest expression level among the scale- tributed across all tissue types in the Scaly-foot Snail. Although secreting epithelium highly expressed genes was chitin-binding there were relatively slightly more novel genes in biomineralising peritrophin-A. This gene (CsqKR_Scaf21_18.25), together with tissues, the presence of novel genes was not restricted to areas of its recent duplicated paralogue (CsqKR_Scaf21_18.24), were also obvious morphological adaptations. There were fewer novel genes highly expressed in the shell-secreting mantle but at lower levels in the scale-secreting epithelium than the shell-secreting mantle, than in the scale-secreting epithelium. They also show sequence and fewer still in the symbiont-hosting oesophageal gland similarity with the well-known shell matrix protein pif , which (Fig. 2c). Although these novel genes do not correspond to other by contrast was highly expressed in the mantle in the Scaly-foot previously annotated genes, and were newly discovered within the Snail C. squamiferum, similar to other Mollusca, but not in the Scaly-foot Snail genome, they may yet be present in other taxa. scale-secreting epithelium. A similar pattern is seen in the chitin This is only the first genome in the subclass Neomphaliones, and synthase gene family, as the Scaly-foot Snail genome possesses the poor taxon coverage and limited completeness of existing diverse subgroups of this family, and different paralogues were molluscan genomes do not provide a reliable reference framework found to be highly expressed in the mantle or scale-secreting to infer whether these ‘novel’ genes are truly lineage-specific epithelium (Supplementary Fig. 6). These chitin synthase genes 6 NATURE COMMUNICATIONS | (2020)11:1657 | https://doi.org/10.1038/s41467-020-15522-3 | www.nature.com/naturecommunications TPM TPM NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15522-3 ARTICLE Octopus bimaculoides 28.91 Mizuhopecten yessoensis Azumapecten farreri 470.06 548.49 Crassostrea gigas 334.95 Pinctada fucata 431.95 Modiolus philippinarum 100.08 Bathymodiolus platifrons 529.62 385.70 Haliotis discus hannai 453.12 Chrysomallon squamiferum Lottia gigantea 613.35 483.94 92.18 Lanistes nyassanus Pomacea canaliculata 404.51 Aplysia californica 196.54 116.84 Radix auricularia Biomphalaria glabrata 509.50 Capitella teleta 567.27 Notospermus geniculatus 503.51 Lingula anatina Phoronis australis C O S D C P T J K Pg N 600 500 400 300 200 100 0 Ma Mantle Mantle Mantle Mantle Fig. 6 Genomic and transcriptomic comparisons across Lophotrochozoa and the biomineralisation toolkit. a Genome-based phylogeny of selected taxa showing the position of the Scaly-foot Snail among lophotrochozoans and divergence times among molluscan lineages. Error bars indicate 95% confidence levels. The Scaly-foot Snail is highlighted in red, Gastropoda in pink, Mollusca in blue, and Lophotrochozoa in orange. b Transcription factors shown to be involved in armature secretion in the Scaly-foot Snail, compared with conchiferan molluscs (bivalves, gastropods, cephalopods), aculiferan molluscs (chitons), brachiopod, and annelids. Transcription factors only verified for the Scaly-foot Snail are highlighted in red, those known to be shared across Mollusca in blue, and those shared across Lophotrochozoa in orange. The top row for each group (darker shading) shows records of significant expression in shell-secreting mantle and bottom row (lighter shading) shows expression for other hard parts such as scales. Ma, million years ago. Source data are provided as a Source Data file. are clearly involved in the biomineralisation process, and may second population from the iron-poor Solitaire vent field that have been co-opted in the evolution of the sclerites. naturally lack any iron sulfide mineral coating . In the Kairei Although some highly expressed genes were different in the population, with iron sulfide mineralisation, transmembrane two biomineralising tissues, positions on the genome of transporter activity was comparatively enriched in highly transcription factors involved in the biomineralization toolkit expressed genes (Supplementary Data 8), supporting the were found to be synchronous with expression patterns seen in hypothesis the iron sulfide precipitation is indeed mediated by the both mantle and scale-secreting epithelium. We propose that they gastropod through transportation and precipitation of sulfur control downstream novel and existing biomineralization genes species . Another gene, metal tolerance protein (MTP) 9, in different ways to fabricate different kinds of armour or skeletal exhibited over 27-fold increased expression level (Supplementary elements by upstream activation. Fig. 7) in the population with iron sulfide mineralisation com- pared with the one without. This gene is widely found in inver- tebrates, protists, and even plants . Although poorly studied in invertebrates, functional assays of MTP in plants revealed Iron sulfide biomineralisation. We compared gene expression potential pathways for enhanced tolerance of metal ions and levels in the scale-secreting epithelium and shell-secreting mantle maintaining intracellular homoeostasis . of the Scaly-foot Snail from the iron-rich Kairei vent field with a NATURE COMMUNICATIONS | (2020)11:1657 | https://doi.org/10.1038/s41467-020-15522-3 | www.nature.com/naturecommunications 7 Hox4 Hox1 En Post1 Hox5 Arx Brachyury ETS Dlx Gsc Msx Six3/6 Hox2 Antp Zic Evx Mox grh Hox3 Gbx Sox14 Lox4 Lox5 Pax3/7 Pax6 Mollusca Lophotrochozoa ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15522-3 The Scaly-foot Snail biomineralises iron sulfide nanoparticles 50 Gb of Illumina Novaseq reads with the paired-end mode and a read length of 150 bp. by allowing sulfur that it actively recruits or deposits in the scales to react with iron ions diffusing in from its highly iron-enriched Oxford Nanopore Technologies (ONT) library preparation. A total of 2–3 μg environment , and the MTP9 gene likely helps the snail tolerate HMW genomic DNA in 10 mM Tris-HCl (pH 8.0) were used for each library such high levels of iron in its surrounding environment. This preparation. All libraries were prepared using the Ligation Sequencing Kit 1D would allow the Scaly-foot Snail to gain finer-scale control of (SQK-LSK108, ONT, Oxford, UK). The standard protocols [1D gDNA selecting for nanomaterial-scale production, which is completed at much long reads (SQK-LSK108) Protocol] from Oxford Nanopore Technologies were modified and performed as follows. lower temperatures than can be currently controlled in laboratory 12 The optional DNA repair step was not performed. End repair and dA-tailing settings . Many other species of molluscs incorporate iron into were performed using the NEBNext Ultra II End-Repair/dA-tailing Module (NEB hard parts, particularly the radula . Molluscs have been high- E7546). A total volume of 120 μl reaction mixture included 14 μl Ultra II End-Prep lighted as a model for generating biogenic nanomaterials at low buffer, 6 μl Ultra II End-Prep enzyme mix, and 100 μl genomic DNA. The reaction temperatures , but the genomic tools to open this part of the mixture was incubated at 20 °C for 30 min and 65 °C for 20 min using a thermal cycler. Clean-up was performed using a 0.4× volume AMPure XP beads (48 μl), biomineralisation toolkit for other applications were previously incubated at room temperature with gentle stirring for 5 min, washed twice with unavailable. 200 μl fresh 70% ethanol, and briefly air-dried for 1 min to obtain the pellet. DNA was eluted by adding 31 μl of nuclease-free water (NFW), resuspending the beads, and incubating for 10 min at 37 °C. One microlitre of aliquot was quantified by Discussion Qubit to ensure that ≥1.5 μg of DNA were retained. By comparing genomic elements and tissue-specific gene Ligation was then performed by gently mixing 20 μl Adaptor Mix (AMX1D, expression patterns in the Scaly-foot Snail scale- and shell- SQK-LSK108, ONT), 50 μl NEB Blunt/TA Master Mix (NEB M0367), and 30 μl dA-tailed DNA while incubating at room temperature for 30 min. The adaptor- secreting tissues, as well as other biomineralising tissues in ligated reaction was cleaned up with a 0.6 × volume (60 μl) of AMPure XP beads, lophotrochozoans, we revealed an ancient biomineralisation incubated for 5 min at room temperature, and followed by resuspending the pellet toolkit comprising at least 25 transcription factors that contribute in 500 μl Adapter Bead Binding buffer (ABB, SQK-LSK108, ONT). The purified- to biomineralisation across all lophotrochozoan hard parts ligated DNA was eluted using 15 μl of Elution Buffer (ELB, SQK-LSK108, ONT), resuspending beads, and incubating for 10 min at 37 °C. One microlitre of aliquot investigated to date. Gene families in the Scaly-foot Snail genome was quantified by Qubit to ensure that ≥500 ng of DNA were retained. The aliquot have predominantly ancient origins, as seen in other lopho- of the adapted and tethered DNA (the pre-sequencing Mix) was used for loading trochozoans (Fig. 3b), but their distribution and duplications into MinION Flow Cell. across various lineages are nonsynchronous with phylogenetic positions (Fig. 3a), underlining rapid modifications. Although MinION sequencing. MinION sequencing was performed as per manufacturer’s novel genes do appear to play important roles in downstream guidelines using R9 flow cells (FLO-MIN106, ONT). Priming of the SpotON Flow production of the hard parts, the hard parts themselves arise by Cell was preformed following the standard protocol (1D gDNA selecting for long reads (SQK-LSK108) Protocol). Directly after priming, 75 μl of the prepared library deploying elements from the conserved biomineralisation toolkit. mixed with 12 μl of the pre-sequencing Mix (adapted and tethered DNA library), Comparison between sclerites from unrelated groups of molluscs 25.5 μl of Library Loading Bead (LLB, ONT), 35 μl of Running Buffer Fuel Mix —the Scaly-foot Snail, and aculiferan molluscs, or Cambrian (RBF, ONT), and 2.5 μl of NFW were loaded through the SportON sample port in fossils—may underline the longevity of these gene families. a dropwise fashion. MinION sequencing was operated with MinKNOW 1.3.23 and fastq files were base-called with Albacore v2.3.4 (ONT) with the default setting. The true power of the biomineralisation toolkit lies in the Reads <3 kb were discarded. capacity for dynamic combination of the components being switched on or off, expanded or reduced, and relocated within the Genome assembly. The Illumina reads were trimmed with Trimmomatic v0.33 . genome, which creates compounding changes in phenome. The They were further assembled by Platanus v1.2.4 using the following settings: –k31 extreme morphological disparity of mollusc biomineralisation –u 0.2 –t29 –s 10. The genome size was estimated to be 444.4 Mb using the 17-mer underpins their successful diversification. Similarly, the re-use, re- histogram generated (Supplementary Note 1). The histogram was also submitted to arrangement and re-deployment of conserved genomic elements GenomeScope (http://qb.cshl.edu/genomescope/); and the genome heterozygosity 37 was estimated to be 1.38%. This indicates that the C. squamiferum genome is explains both our challenges over more than 540 million years relatively compact in Mollusca and also highly heterozygous. The total size, N50, in obtaining genomic and phylogenetic resolution, and their and mean length of the assembled genome from Platanus were 469.0 Mb, 18.2 kb evolutionary success. and 1.8 kb, respectively, indicating high fragmentation. We then used a range of bioinformatics pipelines to assemble the genome with ONT reads, including the ONT-only approaches (i.e. [canu version 1.7 ], [canu Methods 40 version 1.7 + smartdenovo (https://github.com/ruanjue/smartdenovo) ], and Sample collection. Specimens of the Scaly-foot Snail Chrysomallon squamiferum [minimap2 version 2.17-r943-dirty + miniasm version 0.3-r179], etc.) and the used in the present study were collected by the manned submersible Shinkai 6500 Illumina + ONT hybrid approach (e.g. MaSuRCA version 3.2.6). The detailed on-board multiple deep-sea research expeditions of R/V Yokosuka. For genome codes and settings of each assembly pipeline are detailed in Supplementary Note 2. sequencing, a specimen collected from the Kairei hydrothermal vent field (25° Following comparison of assembly statistics of different pipelines 19.23ʹS, 70°02.42ʹE, 2415 m depth, cruise YK16-02E, Feb 2016) and immediately (Supplementary Note 2), the genome assembled from the canu + smartdenovo placed into −80 °C upon recovery was used. For gene expression analyses, speci- pipeline was the best one and therefore used for downstream analyses. The mens fixed in situ using RNA stabilising agent from both the Kairei field (25° assembly was firstly error corrected five times with the ONT fastq file by Racon 19.23ʹS, 70°02.42ʹE, 2415 m depth, cruise YK13-03, Mar 2013) and the Solitaire 41 42 v1.4 and sequentially corrected twice with Illumina reads using Pilon v1.13 . field (19°33.41ʹS, 65°50.89ʹE, 2606 m depth, cruise YK13-02, Feb 2013) were used. The resultant error-corrected assembled genome had a total size, N50 and mean Specimens for in situ hybridisation were collected at the same time from the size of 407.8 Mb, 1.91 Mb and 393.7 kb, respectively. Analysis of the genome Solitaire field, but was fixed in 4% paraformaldehyde (PFA) solution upon recovery assembly completeness with metazoan Benchmarking Universal Single-Copy and later transferred to 80% ethanol. 43 Orthologs (BUSCO) v3.0 revealed that the completeness of the genome is 96.6%, only 2.5% of the reported metazoan genes were missing, and confirmed the presence of single-copy metazoan genes. High molecular weight DNA extraction and quantification. The foot and mantle of a Kairei field Scaly-foot Snail (specimen code E02B1) were used for DNA extraction. High Molecular Weight (HMW) DNA was extracted using the Microbial sequence contamination removal. Genome binning was performed by MagAttract HMW DNA Kit (Qiagen, Hilden, Germany) according to the manu- using MetaBAT 2 with the assembled contigs as input to check the microbial facturer’s instructions. The HMW DNA was further purified and concentrated sequence contamination. The resulting output genomes were examined for com- using the Genomic DNA Clean & Concentrator (gDCC-10) kit (ZYMO Research, pleteness and potential contamination using CheckM v1.0.7 , based on the pre- Irvine, CA, USA). DNA quality was assessed by running 1 μl of the sample on a sence of particular marker gens. Open reading frames (ORFs) were predicted using BioDrop µLITE (BioDrop, Holliston, MA, USA) to ensure the purified of DNA Prodigal v2.6.3 and only ORFs with closed end were retained. Phylogenetic with the OD 260/280 of 1.8 and the OD 260/230 of 2.0–2.2. Concentration of DNA analysis was performed based on 16 bacterial single-copy marker genes (frr, rplB, was assessed using the dsDNA HS assay on a Qubit fluorometer v3.0 (Thermo rplC, rplD, rplE, rplF, rplN, rplS, rpmA, rpsB, rpsC, rpsE, rpsJ, rpsS, smpB and tsf) Fisher Scientific, Singapore). A total of 1 μg DNA was used to obtain approximate share by all the genomes included for analyses. Protein sequence of the 16 genes 8 NATURE COMMUNICATIONS | (2020)11:1657 | https://doi.org/10.1038/s41467-020-15522-3 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15522-3 ARTICLE 47 61 were aligned individually by ClustalW in MEGA 6 , and then linked together to BLASTp results were used to assign the OGs by OrthoMCL v2.0.9 pipeline . OGs construct a maximum likelihood tree. For HGT analysis, all predicted protein from selected lophotrochozoan taxa (Fig. 6a) were used for the phylogenomic sequences from the two symbionts were BLASTP searched against the NCBI NR analysis. Only single-copy genes in each OG and genes that can be found in at least database (https://www.ncbi.nlm.nih.gov/refseq/) using an e-value 1e−5 with 50 60% of taxa (i.e. at least 11 species) were retained for downstream phylogenomic best hits. The BLASTP output files were queried into MEGAN 6 for taxonomy analysis, resulting in 1375 OGs. Gene sequences within each OGs were aligned by analysis in a latent class analysis (LCA) model. Summary of removed microbial MUSCLE, and the bona fide alignments were kept after trimming by TrimAL . sequences can be found in Supplementary Note 2. These alignments were concatenated, and the total alignments which contained missing sequences included 435,071 distinct alignment patterns across 19 species. The phylogenetic tree was conducted by RAxML v8.2.4 with the partition Hi-C sequencing and genome scaffolding. Another individual of the Scaly-foot information of each orthologue gene and GTR + gamma model. The final tree file Snail from the same locality, the Kairei field, and the same collection lot (specimen was viewed by FigTree v1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/). All boot- code E02B2) was used for Hi-C library preparation. For Hi-C library preparation strap values were 100, indicating full support. from the Chrysomallon squamiferum specimen E02B2 from Kairei field, the fol- MCMCtree , part of Phylogenetic Analysis by Maximum Likelihood (PAML) lowing methods were used, modified from Lieberman-Aiden et al. . Approxi- v4.9 (http://abacus.gene.ucl.ac.uk/software/paml.html), was used to predict the mately 2 g wet weight of foot tissue stored in −80 °C freezer was thawed on ice and divergence time among molluscs, with nodes constrained by fossil records and resuspended with 37% formaldehyde in serum-free DMEM for animal chromatin geographic events . The following fossil records and geographic events were used cross-linked. Then, the suspended tissues were homogenized and incubated at to constrain the nodes in the MCMC tree: A hard max time-point of 150 Ma for L. room temperature (RT) for 5 min, and glycine was added to a final concentration nyassanus and P. canaliculata, which correspond to the split of South America and of 0.25 M. The solution was incubated at RT for another 5 min and transferred on Africa ; minimum = 168.6 Ma and soft maximum = 473.4 Ma for A. californica ice for 15 min. The cells were further lysed in a pre-chilled lysis buffer which and B. glabrata ; hard minimum bound = 390 Ma for Caenogastropoda and includes 10 mM NaCl, 0.2% IGEPAL CA-630 (Sigma-Aldrich), 1X protease Heterobranchia ; minimum = 470.2 Ma and soft maximum = 531.5 Ma for A. inhibitor in 10 mM PH = 8.0 Tris-HCL buffer. The chromatin digestion, labelling californica (or B. glabrata) and L. gigantea ; and minimum = 532 Ma and soft and ligation with biotin steps followed Lieberman-Aiden et al. . The protein and maximum = 549 Ma for the first appearance of molluscs ; hard minimum = biotinylated free-ends were removed, and DNA was purified and sequenced by 465.0 Ma for the first appearance of Pteriomorpha ; and minimum = 550.25 Ma Illumina Novaseq platform with the paired-end mode and a read length of 150 bp. and soft maximum = 636.1 Ma for the first appearance of Lophotrochozoan . The The Hi-C raw Illumina reads were trimmed with Trimmomatic v0.33 and best protein substitution model was LG + I + G and was employed to each site. then assessed by HiC-Pro v2.10.0 . Only the valid reads generated by HiC-Pro The burn-in, sample frequency, and number of samples were set as 10,000,000, were further processed by the Juicer 1.5.6 pipeline . Then, genomic scaffolding 1000, and 10,000, respectively. was conducted with the 3D de novo assembly pipeline using the default diploid parameters. Several manual corrections were done in Juicebox to ensure the scaffolds within the same pseudo-chromosomal linkage groups met the Hi-C Gene family analyses. The OGs deduced from OrthoMCL v2.0.9 were also used for linkage characteristics. In this process, three contigs were split into two parts and 69,70 the gene family expansion and contraction analysis by CAFÉ v3.1 pipeline .Only anchored to different chromosomes. In total, 1025 contigs were scaffolded into 15 those with gene family wide P value <0.01 and a taxon-specificViterbi P value <0.01 pseudo-chromosomal linkage groups, and only seven contigs were not anchored were considered as an event of expansion/contraction. due to insufficient Hi-C linkage found on them. The Hox gene cluster found in two To determine the gene age of the OGs, the OrthoMCL result was also used to contigs in the pre-Hi-C version was linked into an intact one, suggesting a good deduce the time of origin, after removing taxon-specific OGs. OGs with genes only performance of the Hi-C scaffolding method. The final genome assembly statistics found in gastropods were classified as Gastropoda-specific. The similar approach can be found in Supplementary Table 1. was applied to assign OGs as Mollusca-specific, Lophotrochozoa-specific, Bilaterian-specific and Eumetazoan-specific. A binary (present/absent) matrix was also deduced from the OrthoMCL results, and taxon-specific OGs were further Genome quality check and repeats identification. Quast v4.0 was used to check excluded. The binary matrix was used to plot a correlation matrix heatmap using genome assembly quality with ONT reads, for assembly assessment report see the corrplot library in R 3.5.2 . Supplementary Note 2. Repeats and transposable elements were annotated using In order to examine the chromosomal distribution of the expanded genes or the the RepeatModeler 2.0 and RepeatMasker 4.0.8 pipeline with the searching genes that are highly expressed in each tissue, a hypergeometric test was applied to programme of NCBI RMblast v2.9.0. The species-specific repeat library was assess whether there was any bias in any particular chromosome, with the annotated with RepeatModeler. Afterwards, RepeatMasker was run twice, one assumption that genes are randomly distributed in each chromosome. The using the species-specific repeat library and another one using the repeats in distribution in the expanded gene families is summarised in Supplementary Note 3. Repbase (https://www.girinst.org/repbase/). All the results were summarized and classified with the perl script buildSummary.pl in the RepeatMasker package, for summary table see Supplementary Note 2. Transcriptome sequencing. The Scaly-foot Snail used in genome sequencing (specimen code E02B1) was also dissected into a number of tissues/organs, namely digestive gland, scales, foot muscle, ctenidium (gill), mantle, nerve cord, nephri- Genome annotations. Two versions of transcriptome assemblies were generated: dium, endosymbiont-containing oesophageal gland, testis and ventricle, following a (1) the transcriptome reads (see Transcriptome sequencing section below) were de previously published account of the Scaly-foot Snail anatomy . This dataset was novo assembled by Trinity v2.6.5; (2) the transcriptome reads were aligned to the only used for gene prediction, as the sample was not fixed in situ. Meanwhile, for genome using histat2 v2.1.0 , and the aligned .bam file was assembled by Trinity comparative purposes, only in situ fixed samples were considered; two individuals under the genome-guided mode. These two versions of transcriptome were merged from the Kairei field and another three from the Solitaire field were used, targeting by PASA pipeline v2.2.0 with the aligners of BLAT and gmap and were further tissues of interest including scales, foot, gill, mantle and oesophageal gland. clustered with cd-hit-est v4.6 with a minimum sequence identity of 0.95 . RNA was extracted using TRIzol reagents and further sequenced using the Maker v3.0 was initially ran with the transcript evidence alone, and only gene Illumina Novaseq platform with the approximate output of 5 Gb for each tissue, models with an AED score <0.01 were retained. Gene models with <3 exons, with read length of 150 bp and paired-end mode. The raw reads were cleaned with incomplete open reading frame, and with an inter-genic region <3 kb were Trimmomatic v0.36 . Gene expression level in each tissue was quantified by removed. The rest of bona fide gene models were used to train Augustus v3.1, a de 58 Kallisto v0.44.0 with sequence based bias correction. Differentially expressed novo gene predictor . Then, the gene model prediction was performed by using genes were determined by DESeq2 using the normalization method of Loess, a Maker again, but with transcript evidence, protein evidence, Augustus gene minimum read count of 10, and paired test (n = 5). predictions, as well as an automatic annotation integration of these data into a Tissue-specific genes for the tissues were determined based on their expression consensus annotation based on their evidence-based weights. levels compared across all tissue types. To minimise the batch effect from pooling Gene functions were determined by searching the NCBI non-redundant different sample collection events and localities, a paired test method was applied. database (https://www.ncbi.nlm.nih.gov/refseq/) and SwissProt database via Only genes that were overexpressed with a fold change of over 2 and FDR < 0.05 UniProt (https://www.uniprot.org/) with the settings of: -evalue 1e-5 -word_size 3- against other tissue types were classified as highly expressed. Differentially num_alignments 20 -max_hsps 20. The Gene Ontology (GO) functional categories 59 expressed genes in the scales and mantle between the Kairei and Solitaire fields were deduced from the BLAST2GO pro v.4.1.9 software. Kyoto Encyclopedia of were also compared with shed light on the iron sulfide biomineralisation. Genes and Genomes (KEGG) annotation was performed using the KEGG Dominant functions of these target genes were further assessed with GO Automatic Annotation Server (https://www.genome.jp/kegg/kaas/) with the bi- 73 59 enrichment analyses using GOEAST v.1.30 or in the BLAST2GO v.4.1.9 directional best hit (BBH) method. A sensitive HMM scanning method on the package (see Supplementary Note 4 for results). known pfam functional domains with an e-value of 0.05 was also be used to classify the gene families. The low complexity protein was predicted with online XSTREAM v1.73 (https://amnewmanlab.stanford.edu/xstream/). Real-time PCR validation. Real-time PCR was employed to validate gene expression patterns in the two main tissues of interest, scales and mantle, as well as Phylogenetic analyses. The orthologue groups (OGs) were determined by a three selected other tissue types for comparison, including foot muscle, oesophageal BLASTp search against protein sequences of other available high-quality mollus- gland, and gill. The primers for real-time PCR were designed with the on-line can, lophotrochozoan, and metazoan genomes (see Supplementary Note 6). The NCBI Primer-BLAST tool (for a list of primers see Supplementary Note 5). PCR NATURE COMMUNICATIONS | (2020)11:1657 | https://doi.org/10.1038/s41467-020-15522-3 | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15522-3 underlying Figs. 1, 2c, 3b, 4a, 5c and 6b, and Supplementary Figs. 1, 2, 3, 6b, 7c, 8, 9, 10, product length was set within the range of 100–200 bp, and optimal melting 11, 12 and 13 are provided as a Source Data file. Although specimens of Chrysomallon temperature was set as 60.0 °C. Only Primer-pair with the least possibilities of self- squamiferum are available from the authors at a reasonable request, the numbers complementarity and self 3ʹ complementarity was selected for each gene. The available are very limited. elongation factor 1-alpha (EF1a) was selected as the internal standard gene, as its expression level remain almost constant across all the samples. Total RNA was extracted from each type of tissue by Trizol method and trace Code availability amount of DNA was removed with TURBO DNA-free kit (Thermo Fisher All codes, commands, and intermediate files for the bioinformatics analyses carried out Scientific), and the first strand cDNA was synthesized by using High Capacity (using freely or commercially available software, as listed in the Methods section) in the cDNA Reverse Transcription Kit (Applied Biosystems). Real-Time PCR was present study are contained in Supplementary Data 9. performed with SYBR Green RT-PCR Reagents Kit (Applied Biosystems) on LightCycler 480 II (Roche) with the following procedures: (1) polymerase activation at 95 °C for 10 min, and (2) annealing and extending at 57 °C for 1 min Received: 21 March 2019; Accepted: 13 March 2020; with a total of 40 cycles. The specificity of primer pairs for the PCR amplification was checked by the melting curve method. Triplicates were applied for each gene, ΔΔCt 74 and the relative gene expression level was calculated based on the 2 method . In situ hybridisation. To localize the expression of genes involved in scale and mantle formation, in situ hybridisation was performed on scale-secreting epidermis References and the mantle of Scaly-foot Snails collected from the Solitaire field, fixed in 4% 1. Vermeij, G. J. The origin of skeletons. Palaios 4, 585–589 (1989). PFA solution and stored in 80% ethanol. In situ hybridisation was performed 2. Marlétaz, F., Peijnenburg, K. T. C. A., Goto, T., Satoh, N. & Rokhsar, D. S. A according to methods detailed in Miyamoto et al. with slight modifications, new spiralian phylogeny places the enigmatic arrow worms among detailed as follows. gnathiferans. Curr. Biol. 29, 312–318.e313 (2019). Whole-mount in situ hybridisation was carried out for the mantle. Samples 3. Zhuravlev, A. Y. & Wood, R. A. The two phases of the Cambrian Explosion. were rehydrated and washed with PBST (i.e. PBS containing 0.1% Tween 20) three Sci. Rep. 8, 16656 (2018). times. The samples were digested with 10 µg/ml proteinase K/PBST for 30 min at 4. McDougall, C. & Degnan, B. M. The evolution of mollusc shells. Wiley 37 °C. After a brief wash with PBST, the samples were post-fixed in 4% PFA/PBST Interdiscip. Rev.: Developmental Biol. 7, e313 (2018). for 10 min RT (20–25 °C) and washed three times in PBST. Samples were 5. Kocot, K. M. et al. Phylogenomics of Lophotrochozoa with consideration of prehybridised in a prehybridisation solution (50% formamide, 5 × SSC, 5 × systematic error. Syst. Biol. 66, 256–282 (2016). Denhardt’s solution, 100 µg/ml yeast RNA, 0.1% Tween 20) at 60 °C for 4 h and 6. Chen, C., Copley, J. T., Linse, K., Rogers, A. D. & Sigwart, J. How the mollusc then hybridised with a hybridisation solution containing a digoxigenin (DIG)- got its scales: convergent evolution of the molluscan scleritome. Biol. J. Linn. labelled RNA probe at 60 °C, for 3 days. Hybridised samples were washed twice in a Soc. 114, 949–954 (2015). solution of 50% formamide, 4 × SSC, and 0.1% Tween 20 for 30 min each; then 7. Warén, A., Bengtson, S., Goffredi, S. K. & Van Dover, C. L. A hot-vent twice in 50% formamide, 2 × SSC, and 0.1% Tween 20 for 30 min twice; 2 × SSC, gastropod with iron sulfide dermal sclerites. Science 302, 1007–1007 and 0.1% Tween 20 for 30 min each time; and twice in 0.2 × SSC and 0.1% Tween (2003). 20 for 30 min each, at 60 °C. These were then washed with MABT (i.e. maleic acid 8. Vrijenhoek, R. C. On the instability and evolutionary age of deep-sea buffer containing 0.1% Tween 20) three times for 30 min at RT, blocked in 2% chemosynthetic communities. Deep Sea Res. Pt. II 92, 189–200 (2013). blocking reagent (Roche) in MABT for 2 h at RT, and incubated overnight with a 1:1500 dilution of antDIG-AP antibody (Roche) in the blocking buffer at 4 °C. 9. Kaim, A., Jenkins, R. G., Tanabe, K. & Kiel, S. Mollusks from late Mesozoic seep deposits, chiefly in California. Zootaxa 3861, 401–440 (2014). Samples were then further washed six times with MABT for 60 min each on a rocker and then transferred into TNT buffer (100 mM Tris pH 9.5, 100 mM NaCl, 10. Stöger, I. et al. The continuing debate on deep molluscan phylogeny: evidence 0.1% Tween 20). A chromogenic reaction was performed using nitro blue for Serialia (Mollusca, Monoplacophora + Polyplacophora). BioMed. Res. Int. tetrazolium chloride/5-bromo-4-chloro-3-indoly-phosphate (NBT/BCIP; Roche) in 2013, 18 (2013). AP buffer (100 mM Tris pH 9.5, 100 mM NaCl, 50 mM MgCl , 0.1% Tween 20 and 11. Nakamura, K. et al. Discovery of new hydrothermal activity and 2% polyvinyl alcohol) until signals were visible. The reaction was stopped in PBST, chemosynthetic fauna on the Central Indian Ridge at 18°–20°S. PLoS ONE 7, post-fixed in 10% formalin/PBST, rewashed with PBST, mounted with 40% e32965 (2012). glycerol, and observed under a light microscope (IX71, Olympus). 12. Okada, S. et al. The making of natural iron sulfide nanoparticles in a hot vent In situ hybridisation of the scale-secreting epidermis was carried out with snail. Proc. Natl. Acad. Sci. USA, https://doi.org/10.1073/pnas.1908533116 sections. Samples were washed three times with PBS and mounted in Tissue-Tek O. (2019). C.T. compound to use in frozen sectioning. Frozen sections were air-dried, washed 13. Goffredi, S. K., Waren, A., Orphan, V. J., Van Dover, C. L. & Vrijenhoek, R. C. with PBST, and fixed in 4% PFA/PBS for 10 min at RT. The slides were washed Novel forms of structural integration between microbes and a hydrothermal with PBS and digested with 1 µg/mL proteinase K in PBS for 10 min at RT. After a vent gastropod from the Indian Ocean. Appl Environ. Microbiol. 70, brief wash with PBS, the samples were post-fixed in 4% PFA/PBS for 10 min at RT. 3082–3090 (2004). The slides were washed with PBS three times. The samples were prehybridized for 14. Chen, C., Copley, J. T., Linse, K., Rogers, A. D. & Sigwart, J. D. The heart of a at least 1 h in prehybridisation solution (50% formamide, 5 × SSC, 5 × Denhardt’s dragon: 3D anatomical reconstruction of the ‘scaly-foot gastropod’ (Mollusca: solution, 50 µg/mL yeast RNA) at 60 °C and hybridized with a DIG-labelled RNA Gastropoda: Neomphalina) reveals its extraordinary circulatory system. Front probe at 60 °C for 3 days. The slides were washed with a solution of 50% Zool. 12, 13 (2015). formamide and 2 × SSC for 30 min; twice in 2 × SSC for 30 min each; 0.2 × SSC for 15. Sigwart, J. D. Molluscs all beneath the sun, one shell, two shells, more, or 30 min twice at 60 °C. They were further rinsed three times with MAB, blocked in none. Curr. Biol. 27, R708–R710 (2017). 2% blocking reagent/MAB for 2 h at RT and incubated overnight at 4 °C with a 16. Nakagawa, S. et al. Allying with armored snails: the complete genome of 1:1500 dilution of anti-DIG-AP antibody (Roche) in blocking buffer. Finally, they gammaproteobacterial endosymbiont. ISME J. 8,40–51 (2014). were washed six times with MAB for 30 min each and transferred into AP buffer 17. Jain, M. et al. Nanopore sequencing and assembly of a human genome with (100 mM Tris pH 9.5, 100 mM NaCl, 50 mM MgCl and 2% polyvinyl alcohol). A ultra-long reads. Nat. Biotechnol. 36, 338 (2018). chromogenic reaction was performed using NBT/BCIP in AP buffer until a signal 18. Aguilera, F., McDougall, C. & Degnan, B. M. Co-option and de novo gene was visible. The reaction was stopped in PBS, postfixed in 4% PFA/PBS, washed evolution underlie molluscan shell diversity. Mol. Biol. Evol. 34, 779–792 with PBS and mounted with 80% glycerol. The hybridised samples were observed (2017). under a light microscope (IX71, Olympus). 19. Simakov, O. et al. Insights into bilaterian evolution from three spiralian genomes. Nature 493, 526 (2012). Reporting summary. Further information on research design is available in the 20. Wang, S. et al. Scallop genome provides insights into evolution of bilaterian Nature Research Reporting Summary linked to this article. karyotype and development. Nat. Ecol. Evol. 1, 0120 (2017). 21. Wanninger, A. & Wollesen, T. The evolution of molluscs. Biol. Rev. 94, 102–115 (2019). Data availability 22. Schiemann, S. M. et al. Clustered brachiopod Hox genes are not expressed The Chrysomallon squamiferum genome that support the findings of this study have been collinearly and are associated with lophotrochozoan novelties. Proc. Natl deposited in the NCBI Sequence Read Archive under the BioProject number Acad. Sci. USA 114, E1913–E1922 (2017). PRJNA523462, all raw sequencing data, including Illumina and Nanopore reads, are also 23. Shinzato, C. et al. Using the Acropora digitifera genome to understand coral deposited under the same BioProject number. The assembled genome, transcriptome, 76 responses to environmental change. Nature 476, 320 (2011). predicted transcripts, proteins have been deposited in Dryad . Publicly available datasets 24. Mollenhauer, J. et al. DMBT1, a new member of the SRCR superfamily, on used in the study include the following: NCBI NR database (https://www.ncbi.nlm.nih. chromosome 10q25.3–26.1 is deleted in malignant brain tumours. Nat. Genet. gov/refseq/), Repbase (https://www.girinst.org/repbase/), SwissProt database via UniProt 17,32–39 (1997). (https://www.uniprot.org/), and KEGG (https://www.genome.jp/kegg/). The source data 10 NATURE COMMUNICATIONS | (2020)11:1657 | https://doi.org/10.1038/s41467-020-15522-3 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15522-3 ARTICLE 25. Nam, B.-H. et al. Genome sequence of pacific abalone (Haliotis discus hannai): 56. Li, W. & Godzik, A. Cd-hit: a fast program for clustering and comparing large the first draft genome in family Haliotidae. GigaScience 6, gix014 (2017). sets of protein or nucleotide sequences. Bioinformatics 22, 1658–1659 (2006). 26. Kocot, K. M., Aguilera, F., McDougall, C., Jackson, D. J. & Degnan, B. M. Sea 57. Cantarel, B. L. et al. MAKER: an easy-to-use annotation pipeline designed for shell diversity and rapidly evolving secretomes: insights into the evolution of emerging model organism genomes. Genome Res. 18, 188–196 (2008). biomineralization. Front Zool. 13, 23 (2016). 58. Stanke, M. & Morgenstern, B. AUGUSTUS: a web server for gene prediction 27. Chen, C., Uematsu, K., Linse, K. & Sigwart, J. D. By more ways than one: rapid in eukaryotes that allows user-defined constraints. Nucleic Acids Res. 33, convergence at hydrothermal vents shown by 3D anatomical reconstruction of W465–W467 (2005). Gigantopelta (Mollusca: Neomphalina). BMC Evol. Biol. 17, 62 (2017). 59. Götz, S. et al. High-throughput functional annotation and data mining with 28. Cope, J. & Babin, C. Diversification of bivalves in the Ordovician. Geobios 32, the Blast2GO suite. Nucleic Acids Re.s 36, 3420–3435 (2008). 175–185 (1999). 60. Kanehisa, M. & Goto, S. KEGG: Kyoto encyclopedia of genes and genomes. 29. Cunha, T. J. & Giribet, G. A congruent topology for deep gastropod Nucleic Acids Res. 28,27–30 (2000). relationships. Proc. R. Soc. B: Biol. Sci. 286, 20182776 (2019). 61. Li, L., Stoeckert, C. J. Jr. & Roos, D. S. OrthoMCL: identification of ortholog 30. Uribe, J. E., Irisarri, I., Templado, J. & Zardoya, R. New patellogastropod groups for eukaryotic genomes. Genome Res. 13, 2178–2189 (2003). mitogenomes help counteracting long-branch attraction in the deep 62. Capella-Gutierrez, S., Silla-Martinez, J. M. & Gabaldon, T. trimAl: a tool for phylogeny of gastropod mollusks. Mol. Phylogenet. Evol. 133,12–23 (2019). automated alignment trimming in large-scale phylogenetic analyses. 31. Ponder, W. F. & Lindberg, D. R. Towards a phylogeny of gastropod molluscs: Bioinformatics 25, 1972–1973 (2009). an analysis using morphological characters. Zool. J. Linn. Soc. 119,83–265 63. Stamatakis, A., Ludwig, T. & Meier, H. RAxML-III: a fast program for (1997). maximum likelihood-based inference of large phylogenetic trees. 32. Zakrzewski, A.-C. et al. Early divergence, broad distribution, and high Bioinformatics 21, 456–463 (2005). diversity of animal chitin synthases. Genome Biol. Evol. 6, 316–325 (2014). 64. dos Reis, M. & Yang, Z. Approximate likelihood calculation on a phylogeny for 33. Suzuki, M. et al. An acidic matrix protein, Pif, is a key macromolecule for Bayesian estimation of divergence times. Mol. Biol. Evol. 28, 2161–2172 (2011). nacre formation. Science 325, 1388–1390 (2009). 65. Benton, M. J. et al. Constraints on the timescale of animal evolutionary 34. Ricachenevsky, F., Menguer, P., Sperotto, R., Williams, L. & Fett, J. Roles of history. Palaeontol. Electron 18,1–106 (2015). plant metal tolerance proteins (MTP) in metal storage and potential use in 66. Hayes, K. A. et al. Molluscan models in evolutionary biology: apple snails biofortification strategies. Front. Plant Sci. 4, 144 (2013). (Gastropoda: Ampullariidae) as a system for addressing fundamental 35. Hilgers, L., Hofreiter, M., Hartmann, S. & von Rintelen, T. Novel genes, questions. Am. Malacol. Bull. 27,47–58 (2009). ancient genes, and gene co-option contributed to the genetic basis of the 67. Benton, M. J., Donoghue, P. C. J. & Asher, R. J. in The Timetree of Life (eds S. radula, a molluscan innovation. Mol. Biol. Evol. 35, 1638–1652 (2018). Blair Hedges, S. & Kumar, S.) 35–86 (Oxford University Press, 2009). 36. Grunenfelder, L. et al. Stress and damage mitigation from oriented 68. Jörger, K. M. et al. On the origin of Acochlidia and other enigmatic nanostructures within the radular teeth of Cryptochiton stelleri. Adv. Funct. euthyneuran gastropods, with implications for the systematics of Mater. 24, 6093–6104 (2014). Heterobranchia. BMC Evol. Biol. 10, 323 (2010). 37. Kocot, K. M., McDougall, C. & Degnan, B. M. in Physiology of Molluscs, A 69. Han, M. V., Thomas, G. W., Lugo-Martinez, J. & Hahn, M. W. Estimating Collection of Selected Reviews Vol. 1 (eds Saleuddin, S. & Mukai, S.) (Apple gene gain and loss rates in the presence of error in genome assembly and Academic Press, 2017). annotation using CAFE 3. Mol. Biol. Evol. 30, 1987–1997 (2013). 38. Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for 70. Sun, J. et al. Adaptation to deep-sea chemosynthetic environments as revealed Illumina sequence data. Bioinformatics 30, 2114–2120 (2014). by mussel genomes. Nat. Ecol. Evol. 1, 0121 (2017). 39. Koren, S. et al. Canu: scalable and accurate long-read assembly via adaptive k- 71. R Development Core Team. R: A Language and Environment Forstatistical mer weighting and repeat separation. Genome Res. 27, 722–736 (2017). Computing (R Foundation for Statistical Computing, 2008). 40. Schmidt, M. H.-W. et al. De Novo assembly of a new Solanum pennellii 72. Bray, N. L., Pimentel, H., Melsted, P. & Pachter, L. Near-optimal probabilistic accession using nanopore sequencing. Plant Cell 29, 2336–2348 (2017). RNA-seq quantification. Nat. Biotechnol. 34, 525 (2016). 41. Vaser, R., Sovic, I., Nagarajan, N. & Sikic, M. Fast and accurate de novo 73. Zheng, Q. & Wang, X. J. GOEAST: a web-based software toolkit for Gene genome assembly from long uncorrected reads. Genome Res. 27, 737–746 Ontology enrichment analysis. Nucleic Acids Res. 36, W358–W363 (2008). (2017). 74. Livak, K. J. & Schmittgen, T. D. Analysis of relative gene expression data using 42. Walker, B. J. et al. Pilon: An integrated tool for comprehensive microbial real-time quantitative PCR and the 2(-Delta Delta C(T)) method. Methods 25, variant detection and genome assembly improvement. PLoS ONE 9, e112963 402–408 (2001). (2014). 75. Miyamoto, N., Yoshida, M.-a, Koga, H. & Fujiwara, Y. Genetic mechanisms of 43. Simao, F. A., Waterhouse, R. M., Ioannidis, P., Kriventseva, E. V. & Zdobnov, bone digestion and nutrient absorption in the bone-eating worm Osedax E. M. BUSCO: assessing genome assembly and annotation completeness with japonicus inferred from transcriptome and gene expression analyses. BMC single-copy orthologs. Bioinformatics 31, 3210–3212 (2015). Evol. Biol. 17, 17 (2017). 44. Kang, D. D., Froula, J., Egan, R. & Wang, Z. MetaBAT, an efficient tool for 76. Sun, J. et al. Data from: The scaly-foot snail genome and implications for the accurately reconstructing single genomes from complex microbial ancient origins of biomineralised armour. Dryad, https://doi.org/10.5061/ communities. PeerJ 3, e1165 (2015). dryad.24053dn (2020). 45. Parks, D. H., Imelfort, M., Skennerton, C. T., Hugenholtz, P. & Tyson, G. W. CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res. 25, 1043–1055 (2015). Acknowledgements 46. Hyatt, D. et al. Prodigal: prokaryotic gene recognition and translation The authors thank the captain and crew of R/V Yokosuka during cruises YK13-02 initiation site identification. BMC Bioinforma. 11, 119 (2010). (principal scientist: Manabu Nishizawa), YK13-03 (principal scientist: Kentaro Naka- 47. Tamura, K., Stecher, G., Peterson, D., Filipski, A. & Kumar, S. MEGA6: mura) and YK16-02E (principal scientist: Ken Takai), and pilots of DSV Shinkai 6500. molecular evolutionary genetics analysis version 6.0. Mol. Biol. Evol. 30, Takuro Nunoura (JAMSTEC) is gratefully acknowledged for bringing together colla- 2725–2729 (2013). borative efforts. This research was financially supported by the China Ocean Mineral 48. Huson, D. H. et al. MEGAN community edition - interactive exploration and Resource Research and Development Association (DY135-E2-1-03 to P.-Y.Q.), the Hong analysis of large-scale microbiome sequencing data. PLoS Comput. Biol. 12, Kong Branch of South Marine Science and Engineering Guangdong Laboratory e1004957 (2016). (SMSEGL20Sc01 to P.-Y.Q.), the Research Grants Council of Hong Kong (GRF grant No. 49. Lieberman-Aiden, E. et al. Comprehensive mapping of long-range interactions 16101219 to J.S., C.C., and P.-Y.Q.), and a Japan Society for the Promotion of Science reveals folding principles of the human genome. Science 326, 289–293 (2009). Grant-in-Aid for Scientific Research (18K06401 to C.C. and H.K.W.). Illumina 50. Servant, N. et al. HiC-Pro: an optimized and flexible pipeline for Hi-C data sequencing was performed by Novogene (Beijing, China). Sampling in the Mauritian processing. Genome Biol. 16, 259 (2015). EEZ was approved by Ministry of Foreign Affairs, Regional Integration, and Interna- 51. Durand, N. C. et al. Juicer provides a one-click system for analyzing loop- tional Trade, Mauritian Government (Ref. 29/2014; 50/38/24 V2). resolution Hi-C experiments. Cell Syst. 3,95–98 (2016). 52. Dudchenko, O. et al. De novo assembly of the Aedes aegypti genome using Hi- C yields chromosome-length scaffolds. Science 356,92–95 (2017). Author contributions 53. Smit, A. F. A. & Hubley, R. RepeatModeler/RepeatModeler, http://www. P.-Y.Q., K.T., C.C., J.S., Y.T., J.-W.Q. and J.D.S. conceived and designed the project. C.C., repeatmasker.org/ (2008). T.W., H.K.W., K.T., J.D.S., K.I., N.F., K.Y., D.B. and N.M. collected and fixed the samples. 54. Kim, D., Langmead, B. & Salzberg, S. L. HISAT: a fast spliced aligner with low C.C., J.S., T.W. and N.M. dissected the samples. J.S., Y.S., T.X. and Y.L. extracted RNA, memory requirements. Nat. Methods 12, 357–360 (2015). DNA, and performed the ONT sequencing. W.Z. checked bacteria contamination. N.M. 55. Haas, B. J. et al. Improving the Arabidopsis genome annotation using performed in situ hybridisation. R.L. and J.S. performed synteny analysis. J.S. and J.C.H.I. maximal transcript alignment assemblies. Nucleic Acids Res. 31, 5654–5666 performed gene family analyses. J.S. assembled the genome, annotated the genes, and (2003). performed bioinformatics analyses. W.C.W. performed the real-time PCR. C.C., J.D.S., NATURE COMMUNICATIONS | (2020)11:1657 | https://doi.org/10.1038/s41467-020-15522-3 | www.nature.com/naturecommunications 11 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15522-3 J.S. and N.M. interpreted the data and drafted the paper. All authors contributed to the Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in final paper. published maps and institutional affiliations. Competing interests Open Access This article is licensed under a Creative Commons The authors declare no competing interests. Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give Additional information appropriate credit to the original author(s) and the source, provide a link to the Creative Supplementary information is available for this paper at https://doi.org/10.1038/s41467- Commons license, and indicate if changes were made. The images or other third party 020-15522-3. material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the Correspondence and requests for materials should be addressed to K.T. or P.-Y.Q. article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from Peer review information Nature Communications thanks Bernard Degnan, Kevin Kocot the copyright holder. To view a copy of this license, visit http://creativecommons.org/ and David Lindberg for their contribution to the peer review of this work. Peer reviewer licenses/by/4.0/. reports are available. Reprints and permission information is available at http://www.nature.com/reprints © The Author(s) 2020 12 NATURE COMMUNICATIONS | (2020)11:1657 | https://doi.org/10.1038/s41467-020-15522-3 | www.nature.com/naturecommunications

Journal

Nature CommunicationsSpringer Journals

Published: Apr 8, 2020

There are no references for this article.