The Iron-Responsive Genome of the Chiton Acanthopleura granulata

Abstract Molluscs biomineralize structures that vary in composition, form, and function, prompting questions about the genetic mechanisms responsible for their production and the evolution of these mechanisms. Chitons (Mollusca, Polyplacophora) are a promising system for studies of biomineralization because they build a range of calcified structures including shell plates and spine- or scale-like sclerites. Chitons also harden the calcified teeth of their rasp-like radula with a coat of iron (as magnetite). Here we present the genome of the West Indian fuzzy chiton Acanthopleura granulata, the first from any aculiferan mollusc. The A. granulata genome contains homologs of many genes associated with biomineralization in conchiferan molluscs. We expected chitons to lack genes previously identified from pathways conchiferans use to make biominerals like calcite and nacre because chitons do not use these materials in their shells. Surprisingly, the A. granulata genome has homologs of many of these genes, suggesting that the ancestral mollusc may have had a more diverse biomineralization toolkit than expected. The A. granulata genome has features that may be specialized for iron biomineralization, including a higher proportion of genes regulated directly by iron than other molluscs. A. granulata also produces two isoforms of soma-like ferritin: one is regulated by iron and similar in sequence to the soma-like ferritins of other molluscs, and the other is constitutively translated and is not found in other molluscs. The A. granulata genome is a resource for future studies of molluscan evolution and biomineralization.


SIGNIFICANCE STATEMENT 23
Chitons are molluscs that make shell plates, spine-or scale-like sclerites, and iron-coated teeth. 24 Currently, all molluscs with sequenced genomes lie within one major clade (Conchifera). Sequencing the 25 genome of a representative from the other major clade (Aculifera) helps us learn about the origins and 26 evolution of molluscan traits. The genome of the West Indian Fuzzy Chiton, Acanthopleura granulata, 27 reveals chitons have homologs of many genes other molluscs use to make shells, suggesting all molluscs 28 share some shell-making pathways. The genome of A. granulata has more genes that may be regulated 29 directly by iron than other molluscs, and chitons produce a unique isoform of a major iron-transport 30 protein (ferritin), suggesting that chitons have genomic specializations that contribute to their 31 production of iron-coated teeth. 32

INTRODUCTION 34
Animals construct hardened structures by combining organic and inorganic components, a process 35 termed biomineralization. To do so, they secrete proteins that initiate and guide the crystallization of 36 inorganic molecules. Animals also incorporate proteins into biomineralized structures, enhancing their 37 strength and flexibility (Cölfen 2010). Molluscs have long been models for studying the genetic 38 mechanisms associated with biomineralization because they craft a wide range of materials into shells, 39 spines, scales, and teeth (McDougall & Degnan 2018). The ability of molluscs to produce diverse 40 biomineralized structures likely contributes to their remarkable morphological and ecological diversity. 41 Chitons (Polyplacophora, Figure 1A Chitons biomineralize their teeth ( Figure 1D) from a unique combination of materials. Most molluscs 53 have a feeding organ, the radula, that bears rows of teeth built from chitin and, in many species, 54 hardened with minerals such as calcium carbonate or silica. Chitons instead harden the cores of their 55 teeth with calcium phosphate (as apatite), and then reinforce their cutting edges with iron (as 56 magnetite) (Lowenstam 1962). These iron coatings allow chitons to scrape algae from rocks without 57 rapidly dulling or damaging their teeth. Chitons produce new teeth throughout their lives, making new 58 rows within days (Shaw et al. 2002;Joester & Brooker 2016). To make new teeth, chitons continuously 59 sequester and transport high concentrations of iron (Kim et al. 1989;Shaw et al. 2002Shaw et al. , 2010. 60 Continuous iron biomineralization may thus present a physiological challenge to chitons because free 61 iron causes oxidative stress (Dixon & Stockwell 2014). 62 To date, most investigations of biomineralization in molluscs have focused on species from the classes 63 Bivalvia and Gastropoda. These, together with Monoplacophora, Cephalopoda, and Scaphopoda make 64 up the clade Conchifera. The sister clade to Conchifera is Aculifera, which is made up of Polyplacophora 65 and Aplacophora. Conchifera and Aculifera diverged approximately 550mya (Vinther et To advance the study of molluscan evolution and to better  76  understand the genetic mechanisms of biomineralization, we sequenced the genome of the West Indian  77 fuzzy chiton Acanthopleura granulata. Exploring the A. granulata genome allowed us to: 1) identify 78 genes chitons may use to build their shell plates, sclerites, and teeth; 2) seek genomic signatures 79 associated with the biomineralization of iron and the mitigation of iron-induced oxidative stress; and 3) 80 better understand the origin and evolution of biomineralization in molluscs. 81

RESULTS AND DISCUSSION 82
We sequenced the genome of a single specimen of A. granulata. We combined reads from one lane of 83 Illumina HiSeq X paired-end sequencing (124 Gb of 2 X 150 bp reads,~204X coverage) with reads from 84 four Oxford Nanopore flowcells run on the GridION platform (22.87 Gb, 37X coverage). Using the hybrid 85 assembler MaSuRCA and optical mapping, we produced a haploid genome assembly for A. granulata 86 that is 606.9 Mbp, slightly smaller than the 743 Mbp haploid genome size estimated by flow cytometry 87 (Roebuck 2017 Figure 2). 95 We generated gene models for A. granulata by 1) sequencing transcriptomes from eight different 96 tissues from the same specimen used for genome sequencing, 2) combining these transcriptomes into a 97 single assembly and aligning the combined transcriptome to the genome, and 3) training de novo gene 98 predictors using both our combined transcriptome and protein sequences predicted from the 99 transcriptomes of other aculiferans. Following these steps, we produced a set of 81,691 gene models 100 that is 96.9% complete according to a BUSCO transcriptomic analysis. This score is similar to the 101 completeness score of the A. granulata genome, so it is likely this set of gene models missed few genes, 102 if any, in the genome assembly. However, of the BUSCO genes expected to be single-copy in all animals, 103 17.2% were represented by multiple gene models. Using Markov clustering to eliminate redundant 104 isoforms, we generated a reduced set of 20,470 gene models that is 94.7% complete. In this smaller set 105 of gene models, only 0.5% of the BUSCO genes have multiple copies, supporting Markov clustering as an 106 effective method for reducing the redundancy of gene models. To characterize proteins based on shared 107 functional domains and sequence similarity, we analyzed the set of 20,470 gene models with 108 InterProScan. We identified at least one GO term for 12,301 genes and a Pfam match for 15,710 genes. 109 We also conducted a KEGG analysis and identified 7,341 proteins that could be assigned to putative 110 molecular pathways. 111 To provide a robust dataset for phylogenetic analysis and gene family evolution analyses, we identified 112 homologous genes shared between A. granulata and other molluscs. We used the larger set of gene 113 models from A. granulata to ensure a more complete input dataset, knowing that any duplicate gene 114 models for the same locus would cluster within the same orthologous group. We compared gene 115 4 models from the A. granulata genome to those from the genomes of nineteen other lophotrochozoans, 116 including fourteen molluscs, two annelids, one brachiopod, one phoronid, and one nemertean. This 117 resulted in 59,276 groups of homologous sequences including 3,379 found in all 20 genomes. 118 We used a tree-based approach to identify orthologous genes shared among all 20 taxa and 119 reconstructed molluscan phylogeny using the 2,593 orthologs present in at least 17 of the 20 genomes 120 we searched. This dataset totaled 950,322 amino acid positions with 16.2% missing data. We recovered 121 A. granulata as the sister taxon of all other molluscs with sequenced genomes ( Figure 1F). We 122 conducted an additional phylogenetic analysis that included more taxa by using transcriptomes in 123 addition to genomes and recovered Acanthopleura within the family Chitonidae in the order Chitonida, have high heterozygosity because this species is a broadcast spawner with a wide geographic range 133 (Glynn 1970). To compare heterozygosity across molluscs, we selected a set of high-quality molluscan 134 genomes for which short-read data are available (Supplementary Table 1). Using k-mer based analysis, 135 we found the highest heterozygosity among the seven genomes we analyzed was 3.15% in the blood 136 clam S. broughtonii, and the other genomes had heterozygosities between those of A. granulata and S. 137 broughtonii. Our findings indicate that heterozygosity may be influenced by more than an animal's 138 reproductive mode, larval type, and geographic range (Supplementary Table 2), and that molluscan 139 genomes should not be assumed to have high heterozygosity. 140 The A. granulata genome is arranged differently than other molluscan genomes and has fewer repetitive 141 elements. Compared to a non-molluscan lophotrochozoan, Lingula anatina (a brachiopod), A. granulata 142 has more repetitive elements of certain types in its genome. Conversely, A. granulata has fewer 143 repetitive elements in its genome than any conchiferan mollusc (Supplementary Table 2). This suggests 144 multiple proliferations of repetitive elements during molluscan evolution. Repetitive elements 145 contribute to structural changes in genomes by providing breakpoints that increase the likelihood of 146 chromosomal rearrangements (Weckselblatt & Rudd 2015). Consistent with this prediction, synteny is 147 lower between A. granulata and all conchiferan molluscs we examined than it is between any two of 148 these conchiferans, and the genomes of conchiferans and A. granulata have little synteny with the 149 genome of L. anatina (Supplementary Figure 5). A recent study found greater conservation of bilaterian 150 ancestral linkage groups (ALGs) between a scallop and non-molluscs than between other bivalves and 151 non-molluscs, suggesting substantial genomic rearrangements occurred within bivalves (Wang et al. 152 2017). Molluscan genomes appear to rearrange frequently across evolutionary time, and perhaps 153 rearrange more frequently in conchiferans due to the proliferation of repetitive elements. However, 154 limited synteny between A. granulata and the scallop genome suggests that the chiton genome has also 155 5 undergone significant rearrangement relative to the hypothetical ancestral bilaterian genome 156 organization. 157 The Hox cluster is a widely conserved set of regulatory genes that contribute to the patterning of the 158 anterior-posterior axes in bilaterian animals. In lophotrochozoans, the genes are typically collinear, 159 beginning with Hox1 and ending with Post1. Although several gastropods and bivalves possess intact 160 Hox clusters, this cluster is dispersed in some bivalves and cephalopods ( structures. 172

Acanthopleura granulata shares many biomineralization genes with conchiferan molluscs 173
We expected chitons to lack many genes previously identified in molluscan biomineralization pathways 174 because their shell plates and sclerites lack both calcite and nacre. We were surprised to find homologs 175 in the A. granulata genome of many biomineralization genes known from conchiferans (Supplementary  176  Table 4). For example, we found an ortholog in A. granulata for Pif. In pterid bivalves, the Pif mRNA 177 encodes a protein that is cleaved into two peptides, PIF97 and PIF80 (Suzuki et al. 2009). These peptides 178 have different roles in biomineralization: PIF80 binds nacre and aids in nacre formation (Suzuki et al. 179 2009), whereas PIF97 binds to chitin and guides the growth of calcium carbonate crystals (Suzuki et al. 180 2013). A recent study in the gastropod Lymnaea stagnalis identified a PIF-like protein that does not 181 promote nacre production, but is expressed in tissues that suggest this PIF protein still plays a role in 182 shell biomineralization (Ishikawa et al. 2020). We found that A. granulata possesses a Pif homolog, but 183 appears to only produce PIF97 rather than two separate peptides. The expression of Pif mRNA was 184 highest in girdle tissue in A. granulata and lowest in the radula, suggesting that PIF peptides may play a 185 role in sclerite formation in chitons (Supplementary Table 4). We hypothesize that the last common 186 ancestor of extant molluscs used PIF97 to help build mineralized structures, and that production of 187 PIF80 is novel to bivalves. 188 The ancestral mollusc likely produced mineralized structures, but whether the ancestral mollusc had a 189 single shell, multiple shell plates, or sclerites remains a matter of debate (Scherholz et  matrices. Chitin is a major component of the extracellular matrices of all molluscan shells and radulae, 196 6 and the A. granulata genome contains genes for chitin synthase, chitinase, and chitin-binding proteins. 197 We also found homologs of lustrin and dermatopontin, two proteins expressed in the extracellular 198 matrices of conchiferans that increase the elasticity and flexibility of their shells (Gaume et al. 2014); 199 Supplementary Table 5).  Table 6). We found 27 of these 31 genes code for proteins with 209 signal peptides, indicating they may be secreted as part of the extracellular matrix during 210 biomineralization (Supplementary Table 6). Among these genes, we found three collagens, one 211 chitinase, and one carbonic anhydrase, all possible contributors to shell formation and repair (Patel 212 2004); Supplementary Table 6). Several of the genes encoding proteins with silk-like domains are highly 213 expressed in the girdle of A. granulata, suggesting a role in the mineralization of sclerites 214 (Supplementary Figure 6). 215

A. granulata has more genes with iron response elements (IREs) than other molluscs 216
Chitons have more iron in their hemolymph than any other animal studied to date (Kim et al. 1988). Iron 217 presents physiological challenges to animals because iron can cause oxidative stress. We hypothesize 218 that the ability of chitons to biomineralize iron requires them to respond quickly to changes in 219 concentration of this potentially toxic metal. To assess the iron-responsiveness of the A. granulata 220 genome, we searched it for iron response elements (IREs), three-dimensional hairpin structures that 221 form in either the 3' or 5' untranslated regions (UTRs) of mRNA molecules and control translation via 222 binding by iron regulatory protein (IRP; Supplementary Figure 7). We also examined IREs in several high-223 quality molluscan genomes that include UTRs as part of their available annotation data. All of the 224 molluscan genomes we examined had similar proportions of 3' to 5' IREs ( Figure 3A). Despite having the 225 fewest gene models, the genome of A. granulata has more IREs than the genomes of any other mollusc 226 we examined. We predicted 271 IREs in the A. granulata genome, compared with an average of 119 IREs 227 across other molluscan genomes (Supplementary Table 7). The highest number of predicted IREs in a 228 conchiferan came from the genome of the blood clam Scapharca broughtonii, which had 201. The blood 229 clam is so named because it is one of relatively few molluscs that produces hemoglobin for use as a 230 respiratory pigment (Kawamoto 1928 because they must absorb and transport larger amounts of iron to produce iron-coated teeth and 233 hemoglobin, respectively. 234 We next asked if the expression of genes with IREs in A. granulata differs between iron-rich and iron-235 poor tissues. In the absence of iron, IRPs bind IREs. When this occurs in the 5' UTR of an mRNA, 236 ribosomes are blocked from translating the protein; thus, mRNAs with 5' IREs will not be translated in 237 the absence of free iron. When IRPs bind IREs in the 3' UTR of an mRNA, they block endonucleases from 238 degrading the mRNA, allowing multiple translations from a single mRNA molecule; thus, the translation 239 of mRNAs with 3' IREs is higher in the absence of free iron (Supplementary Figure 7). We compared the 240 expression of genes with IREs between transcriptomes sequenced from the foot, girdle, ctenidia, and 241 four developmentally-distinct regions of the radula from the same specimen of A. granulata we used for 242 genome sequencing, and found expression of genes with IREs in all tissue types ( Figure 3B). We 243 quantified the expression of genes that contained 3' IREs and 5' IREs separately and found many genes 244 with 5' IREs had higher expression in the three anterior, iron-rich regions of the radula than the rest of 245 the body (Supplemental Figure 8). Genes with 5' IREs that are expressed in iron-rich tissues are likely 246 transcribed into proteins in only those tissues because elsewhere in the body IRPs will bind to IREs and 247 block translation. Consequently, the proteins encoded by these genes may be unique to the radula. 248 After identifying genes in A. granulata that contain 5' IREs and are expressed at relatively high levels in 249 the iron-rich anterior regions of the radula, we asked if these genes might have roles in the 250 biomineralization of the radula. We used GO analysis to compare the molecular functions of the protein 251 sequences coded by the genes with 5' IREs to the protein sequences coded by the full set of genes from 252 the A. granulata genome. We found that genes with a 5' IRE that are highly expressed in the anterior of 253 the radula are more likely to be associated with the molecular functions 'response to inorganic 254 substance', 'response to calcium ion', and 'response to metal ion' (Supplementary Figure 9). This 255 suggests that genes with a 5' IRE that are highly expressed in the radula may be involved in the 256 biomineralization of the apatite (calcium) cores of teeth and their magnetite (iron) caps. A previous 257 study by Nemoto et al. identified a novel biomineralization protein (RTMP1) in the radula of another 258 chiton (Cryptochiton stelleri), and proposed that RTMP1 played a role in iron biomineralization (Nemoto 259 et al. 2019). We examined the mRNA of RTMP1 in C. stelleri and did not detect an IRE in either its 5' or 3' 260 UTR. Thus there are genes that may be important to biomineralization in the chiton radula whose 261 expression levels are not influenced by IREs. 262 Two isoforms of ferritin may provide chitons with tissue-specific protection from oxidative stress 263 All metazoans require iron. However, free iron poses a threat to animals because it catalyzes the 264 production of reactive oxygen species, which inflict damage on DNA and tissues (Dixon & Stockwell 265 2014). To transport iron safely, metazoans use the iron-binding protein ferritin. Previous work suggests 266 that chitons use ferritin to transport iron to their radula (Kim et al. 1988). An iron response element 267 (IRE) is present in the 5' UTR of the heavy chain (or soma-like) ferritin that is expressed by all metazoans 268 (Piccinelli & Samuelsson 2007-7). We found two isoforms of heavy chain ferritin in our gene models for 269 A. granulata: a first isoform (isoform 1) that contains the conserved 5' IRE, and a second isoform 270 (isoform 2) that does not ( Figure 4A). 271 Isoform 1 of ferritin from A. granulata contains an IRE in the 5' UTR, allowing this isoform to be 272 translated only in the presence of free iron. By regulating the translation of ferritin, cells can transcribe 273 ferritin mRNA continuously so they are primed to produce large quantities of ferritin protein rapidly if 274 conditions require it. If no free iron is present, IRP will bind to the IRE and block translation. We found 275 isoform 1 of ferritin is expressed at similar levels in all the transcriptomes we sequenced for A. 276 granulata, including those for the foot, girdle, gonad, ctenidia, and all four regions of the radula ( Figure  277 4B). Thus, when A. granulata needs to bind excess iron, it may be able to rapidly produce isoform 1 of 278 8 ferritin protein throughout its body. We examined other mollusc genomes and transcriptomes and 279 found a ferritin isoform present in all of them that is similar to A. granulata isoform 1 and contains the 280 5' IRE . 281 Isoform 2 of ferritin in A. granulata lacks the 5' IRE present in isoform 1. We identified an alternative 282 transcription initiation site downstream of ferritin exon 1 in the A. granulata genome. Isoform 2 of 283 ferritin, initiated at this downstream site, contains a different exon 1 than isoform 1 of ferritin, but 284 shares exons 2-4 with isoform 1 (Figure 4). We examined other molluscan genomes and transcriptomes 285 and did not find evidence for expression of a ferritin isoform similar to A. granulata isoform 2 (data 286 available on Dryad). In A. granulata, transcripts of isoform 2 are expressed at a lower level than isoform 287 1 throughout all body tissues (foot, girdle, gonad, ctenidia) and in the posterior region of the radula that 288 lacks iron mineralization ( Figure 4B). Expression of isoform 2 is almost undetectable in the iron-rich 289 regions of the radula. Without the 5' IRE, translation of the mRNA of isoform 2 is not blocked in the 290 absence of free iron. The 5' IRE in ferritin is an important regulatory mechanism for protein production. 291 In rats, for example, the expression of ferritin mRNAs is relatively constant across tissues but protein 292 levels vary (Rogers & Munro 1987). Further, mutations in the 5' IRE of ferritin cause hyperferritinaemia 293 in mammals, an iron-related medical condition caused by an overproduction of ferritin protein 294 (Thomson et al. 1999). We hypothesize that chitons use isoform 2 of ferritin to produce a low level of 295 ferritin protein constitutively in tissues outside their radula as protection from the high concentrations 296 of iron circulating through their bodies. 297

Conclusions 298
The A. granulata genome is the first available genome for any chiton or any aculiferan. The information 299 it provides improves our understanding of the evolution of biomineralization across Mollusca as well as 300 lineage-specific innovations within chitons. Chitons are a valuable system for investigating 301 biomineralization because they produce shell plates, sclerites, and iron-clad teeth. The unique 302 combination of structures produced by chitons makes the A. granulata genome a resource for future 303 studies of biomineralization. Although many genes involved in molluscan shell secretion are rapidly 304 evolving (Jackson et al. 2006;Kocot et al. 2016), we were able to identify homologs of many conchiferan 305 biomineralization genes in the A. granulata genome. The expression of several genes associated with 306 conchiferan shell secretion in the girdle of A. granulata suggests these genes may function in sclerite 307 biomineralization in chitons. This suggests a common underlying biomineralization mechanism for 308 conchiferan shells and aculiferan sclerites, structures known to share some developmental pathways 309 even though they arise via different cell lineages (Wollesen et al. 2017). 310 All metazoans require iron, but they must balance iron use against potential oxidative damage. 311 Regulating iron is a particular concern for chitons because they biomineralize their teeth with 312 magnetite. The genome of A. granulata contains more genes with iron response elements (IREs) than 313 the genome of any other mollusc examined to date, indicating it has a larger proportion of genes 314 regulated directly by iron. We identified two isoforms of ferritin in A. granulata, one that is iron-315 responsive and a second that is constitutively translated. We propose the second isoform of ferritin 316 protects tissues outside the radula from oxidative stress by binding free iron. Chitons are an emerging 317 model for studies of both biomineralization and iron homeostasis. The A. granulata genome will aid 318 9 future studies by suggesting specific proteins and pathways to target with comparative studies of gene 319 expression and gene manipulation. 320

Specimen collection and preservation 322
We collected a single male specimen of Acanthopleura granulata from Harry Harris State Park in the 323 Florida Keys (Special Activity License #SAL-17-1983-SR). We cut the majority of the foot into ~1 mm 2 324 cubes and froze them at -80°C. We froze additional pieces of foot, girdle (dissected such that the tissue 325 sample would not contain shell secretory tissue), ctenidia, gonad, and radula in RNAlater and stored 326 them at -80°C as well. 327

Genome and transcriptome sequencing 328
We extracted high molecular weight DNA from frozen samples of foot tissue from A. granulata using a 329 CTAB-phenol chloroform method. We cleaned DNA for short read generation with the Zymo Clean and 330 Concentrator Blobtools identified a large proportion of scaffolds as chordate. We identified contaminants as 362 sequences that differed from the majority of scaffolds in both GC content and coverage and used BLAST 363 to verify these sequences as bacterial before removing them from the assembly. 364 We scaffolded this reduced assembly with one lane of Bionano SAPHYR optical mapping, using two 365 enzymes (BssSI and DLE1) and Bionano Solve v3.4's scaffolding software, which resulted in 87 scaffolds. 366 We ran REAPR v. 1.0.18 (Hunt et al. 2013), which map short read data and collect mapping statistics 367 simultaneously, to determine accuracy of the assembly overall relative to all short-read data generated, 368 and found despite reducing heterozygosity in the final assembly, 85.31% of paired-end reads map 369 perfectly back to the genome assembly, indicating a complete genome assembly relative to the paired-370 end data. 371 To assess our genome assembly, we ran QUAST v. separately to cluster isoforms. We also generated a composite transcriptome of all eight tissues (eight 378 total transcriptomes including four separate radula regions) by combining reads and then following the 379 same process described above. We used this composite transcriptome for annotation. 380

Genome annotation 381
To annotate the A. granulata genome, we first generated a custom repeat library with RepeatModeler v. protein sequences from several other species of chitons that were generated previously (see 385 Supplementary File 1). Using the highest quality gene models from the first as a maker-input gff3 (AED 386 <0.5), we ran a second round of MAKER. From these resulting gene models, we used those with an AED 387 <0

Phylogenetic analyses 446
For species tree reconstruction, in cases where two or more sequences were present for any taxon in a 447 single-gene alignment, we used PhyloPyPruner 0.9.5 (https://pypi.org/project/phylopypruner/) to 448 reduce the alignment to a set of strict orthologs. This tool uses single-gene trees to screen putative 449 orthogroups for paralogy. To build single-gene trees based on orthologs, we trimmed alignments with 450 BMGE v1.12.2 (Criscuolo & Gribaldo 2010) and constructed approximately maximum likelihood trees for 451 each alignment with FastTree2(Price et al. 2010) using the "slow" and "gamma" options. We then used 452 these alignments in PhyloPyPruner with the following settings: --min-len 100 --min-support 0.75 --mask 453 pdist --trim-lb 3 --trim-divergent 0.75 --min-pdist 0.01 --trim-freq-paralogs 3 --prune MI. For datasets 1 454 ("genomes") and 3 ("all_taxa"), only orthogroups sampled for at least 85% of the total number of taxa 455 were retained for concatenation. For dataset 2 ("biomin_subset"), only orthogroups sampled for all 456 eight taxa were retained. Phylogenetic analyses were conducted on the supermatrix produced by 457 PhyloPyPruner of sequences from A. granulata that matched the biomineralization genes of interest, and allowed us to 470 narrow down our list of potential biomineralization genes present in A. granulata. 471 We used the above set of gene queries from other molluscs to identify the ortholog group from the 472 above OrthoFinder2 on the subset of genomes selected as biomineralization representatives across 473 Mollusca. We used the complete CDS or longest mRNA for each gene as a nucleotide query to search 474 our orthogroups, again with an e-value cutoff of 1e-8 to identify the orthogroup(s) likely contained that 475 14 particular biomineralization gene of interest. This produced a list of orthogroups that contained 476 sequences with high similarity to the query, often multiple orthogroups per gene (Supplementary Table  477 4). This was expected due to clustering within OrthoFinder2. We used NCBI BLAST to verify the identity 478 of the orthologous gene sequences by verifying that the top hits for each in BLAST matched the 479 biomineralization gene of interest. We then examined these orthogroups to locate the previously 480 identified A. granulata gene model that matched to each biomineralization protein. The query 481 sequences for each gene sought in A. granulata are available in Supplementary Table 8. to reduce redundancy in the annotations of the A. granulata genome). We then ran SIRE on each of 496 these sets of predicted transcripts. We only accepted predicted IREs scored as "high quality" according 497 to the SIRE metric (indicating both sequence and structural characteristics of a functional IRE). We 498 pulled chiton genes containing a high quality IRE from the eight different tissue transcriptomes 499 generated for genome annotation and assessed expression by mapping each back to the genome with 500 Salmon v. 0.11.3 (Patro et al. 2017-4) to generate quantifications of reads per transcript, and running 501 these quantifications through edgeR (Robinson et al. 2010) to account for transcript length (TPM) and 502 permit direct comparisons of gene expression. We also separated 3' and 5' IREs by subsetting the high 503 quality IREs based on whether the IRE was located at the beginning or end of the sequence. We made 504 heatmaps with log-transformed data to compensate for outliers in expression levels with R package 505 prettyheatmap (Kolde 2012  most anterior region, contains teeth used for feeding; R2 contains teeth that are developed but are not 828 yet used for feeding; R3 contains developing teeth that contain iron oxide; and R4, the most posterior 829 region, contains developing teeth that have yet to be coated with iron. We found lower expression of 830 most IRE-containing genes in the anterior regions of the radula. (LogTPM used to allow both data ranges 831 to appear legibly on the same graph) 832 Relative expression of both isoforms of ferritin across A. granulata tissues. The radula is divided into four 836 developmentally distinct regions as in Figure 3. Isoform 1 is transcribed more highly throughout the 837 body than isoform 2. Isoform 2 is transcribed at lower levels in the anterior (iron-rich) regions of the 838 radula than in other tissues. 839     A. granulata has more IREs than all other molluscs examined, but the relative proportion of 5' and 3' IREs appears consistent across molluscan genomes. (b) The relative expression [log10(TPM)] of transcripts containing IREs in the different tissues of A. granulata. The radula is divided into four developmentally distinct regions: R1,the most anterior region, contains teeth used for feeding; R2 contains teeth that are developed but are not yet used for feeding; R3 contains developing teeth that contain iron oxide; and R4, the most posterior region, contains developing teeth that have yet to be coated with iron. We found lower expression of most IRE-containing genes in the anterior regions of the radula.  . (a) The locations of the transcription initiation sites and exons of isoform 1 of ferritin (orange, above) and isoform 2 of ferritin (blue, below). A 5' IRE (red) is present in the 5' untranslated region of isoform 1, but not in isoform 2. (b) Relative expression of both isoforms of ferritin across A. granulata tissues. The radula is divided into four developmentally distinct regions as in Figure 3. Isoform 1 is transcribed more highly throughout the body than isoform 2. Isoform 2 is transcribed at lower levels in the anterior (iron-rich) regions of the radula than in other tissues.