Expression of calcification‐related ion transporters during blue mussel larval development

Abstract The physiological processes driving the rapid rates of calcification in larval bivalves are poorly understood. Here, we use a calcification substrate‐limited approach (low dissolved inorganic carbon, C T) and mRNA sequencing to identify proteins involved in bicarbonate acquisition during shell formation. As a secondary approach, we examined expression of ion transport and shell matrix proteins (SMPs) over the course of larval development and shell formation. We reared four families of Mytilus edulis under ambient (ca. 1865 µmol/kg) and low C T (ca. 941 µmol/kg) conditions and compared expression patterns at six developmental time points. Larvae reared under low C T exhibited a developmental delay, and a small subset of contigs was differentially regulated between ambient and low C T conditions. Of particular note was the identification of one contig encoding an anion transporter (SLC26) which was strongly upregulated (2.3–2.9 fold) under low C T conditions. By analyzing gene expression profiles over the course of larval development, we are able to isolate sequences encoding ion transport and SMPs to enhance our understanding of cellular pathways underlying larval calcification processes. In particular, we observe the differential expression of contigs encoding SLC4 family members (sodium bicarbonate cotransporters, anion exchangers), calcium‐transporting ATPases, sodium/calcium exchangers, and SMPs such as nacrein, tyrosinase, and transcripts related to chitin production. With a range of candidate genes, this work identifies ion transport pathways in bivalve larvae and by applying comparative genomics to investigate temporal expression patterns, provides a foundation for further studies to functionally characterize the proteins involved in larval calcification.


| INTRODUC TI ON
During bivalve ontogeny, major developmental steps occur rapidly within the first days of life to produce the larval shell, prodissoconch I (PD I). Successful deposition of the PD I shell is essential for subsequent development as the calcified shell provides structural support for swimming and feeding (Galtsoff, 1964).
Calcification begins as the precipitation of a calcium carbonate shell onto an organic template during the trochophore larval stage which is characterized by a free-swimming, ciliated larva (Bayne, 1976). The supply of calcification substrates (Ca 2+ and HCO 3 − ) to the site of calcification by means of transepithelial transport is crucial (Ramesh, Hu, Thomsen, Bleich, & Melzner, 2017). In addition, organic matrix components such as shell matrix proteins and carbohydrates (up to 5% of shell mass, Simkiss & Wilbur, 2012) are synthesized and incorporated into the growing shell, providing the organic template for Ca 2+ deposition. Although comprising a small proportion of the shell content, shell matrix proteins are a diverse assemblage of proteins that are important for polymorph control, crystal nucleation, deposition kinetics, and structure (Marin, Luquet, Marie, & Medakovic, 2008;Miyamoto et al., 2013). Within 2 days of development, bivalve larvae precipitate an inorganic aragonitic shell that is almost equivalent to their somatic mass and covers the entire larval body (Waldbusser et al., 2013).
Accompanying this transition into PD I larval stage is a 250-fold increase in larval calcium content within a few hours (Ramesh et al., 2017). In bivalve larvae, the substrates for calcification (Ca 2+ and HCO 3 − ) are not accumulated prior to calcification but are continually supplied during larval shell formation, most likely by means of transepithelial transport (Ramesh et al., 2017).
The transport of calcium and bicarbonate ions to, and removal of protons from, the calcification space in bivalves is regulated by the expression and activity of membrane-bound ion transport proteins (Ramesh et al., 2017;Sillanpaa, Sundh, & Sundell, 2018). The close relationship between calcification and acid-base-related ion transport processes makes it difficult to distinguish between these processes as intracellular pH regulation is also achieved via movement of protons and bicarbonate across the cell membrane. Primary active ion transporters such as V-type H + ATPases and H + /K + -ATPases utilize energy (ATP) for proton translocation. Secondarily active proton transport is driven by the electrochemical gradient provided by other membrane-bound transport proteins, often the ubiquitous Na + /K + ATPase (NKA). These secondarily active proteins are important for pH i regulation and are involved in the translocation of protons and bicarbonate via Na + -coupled pathways (Na + /H + exchangers SLC9, Na + HCO 3 − cotransporters, SLC4), anion exchangers (Cl − /HCO 3 − exchangers, SLC4), or voltage-gated proton channels.
The cellular pathways elicited to regulate intracellular pH are fundamentally conserved in eukaryotes, and the presence of all relevant transporter families in bivalves has been confirmed by genome and transcriptome sequencing projects in the last decade (Li, Sun, et al., 2016;Murgarella et al., 2016;Takeuchi et al., 2012;Zhang et al., 2012). Several of these membrane-bound proteins have been functionally characterized to be involved in mollusc intracellular pH regulation (Boron & DeWeer, 1976;Ellington, 1993;Sillanpaa et al., 2018;Zange, Grieshaber, & Hans, 1990 Goncalves et al., 2016;Hüning et al., 2013). On the other hand, a number of studies have successfully used RNAi to demonstrate functional roles of specific shell matrix proteins and have shown that their presence is critical to "normal" shell development (Fang et al., 2011;Funabara et al., 2014;Suzuki et al., 2009). Previous studies on mollusc larval development have revealed several shell matrix proteins and genes associated with shell production Liu et al., 2015) and indicate that bivalve larvae utilize notably different shell matrix proteins when compared with adults (Zhao et al., 2018).
The absence of information on the fundamental mechanisms of bivalve larval calcification physiology limits our ability to predict how these organisms can respond and adapt to environmental change. One way to address hypotheses regarding the role of various ion transporters related to mineral formation and organic deposition in the larval calcification pathway is by challenging calcification in a substrate-limited environment. In this study, we used RNA-Seq techniques to identify ion transport and shell matrix proteins involved in the different developmental stages of calcifying larvae of the blue mussel, Mytilus edulis. We hypothesized that the genes involved in calcification in M. edulis would exhibit severe changes in expression related to the rapid rates of PD I calcification in these organisms. Based on previous empirical data on larval mussels which demonstrated that limiting conditions of dissolved inorganic carbon (C T ) elicited strong reductions in larval calcification (Thomsen, Haynert, Wegner, & Melzner, 2015), we used a substrate-limited approach (low C T ) to gain insight into the role of these genes, primarily those involved in inorganic carbon acquisition and crystal formation.  (Stuckas, Stoof, Quesada, & Tiedemann, 2009). We will refer to them as Baltic

| MATERIAL S AND ME THODS
Mytilus edulis-like according to Stuckas et al. (2017). Spawning was induced by exposing the adults to rapidly elevated water temperature between 18 and 25°C. Spawning individuals were separated, and gametes were collected in individual beakers filled with 0.2 µm filtered seawater (FSW). Eggs were fertilized with sperm, and fertilization success was estimated by the presence of polar bodies/cell cleavage. Cleaving embryos were reared in 10-L Duran glass bottles at a density of 10 embryos ml −1 bubbled with pressurized air through plastic tubing. All experiments were performed at 17°C. A total of four separate fertilizations were obtained to conduct replicate experiments (n = 4).  Dickson (1990) and Roy et al. (1993), respectively. Cleaving embryos were added to treated water once pH had increased to stable values (ca. 8.14).

| Experimental treatments and sample collection
Temperature, salinity, and carbonate chemistry parameters of experimental conditions are shown in Table S1.
For each experiment, embryos from single fertilizations were added to both control (FSW) and treatment (low C T ) bottles and sampled at six developmental time points as determined by hours postfertilization (hpf). Samples were collected at 20, 22, 24, 27, 30, and 35 hpf from control bottles, based on time points that were previously identified to be critical for calcification (Ramesh et al., 2017).
To correct for developmental delay in treatment bottles, samples were collected at similar developmental stages to those in control bottles, as determined by frequent microscopic observation of percentage shell cover of the larval body ( Figure S1). Stage 1 occurred prior to calcification while at Stage 2, the onset of calcification was observed, and larvae exhibited a dorsal flattening at the region of the shell field. At Stages 3 and 4, the trochophore shell was observed by the presence of a shell that covered ca. 10% and 49% of the larval body, respectively. Finally, at Stages 5 and 6, larvae had secreted a shell that covered ca. 76% and 100% of the larval body, respectively.
Samples were quickly concentrated on a mesh (55 µm), transferred to 1.5-ml Eppendorf tubes and centrifuged at 10,000 g to form a larval pellet (~6,000 larvae). Seawater was removed using a pipette, and samples were flash frozen in liquid nitrogen and stored at −80°C.
For each sample time point, ca. 100 larvae for photographs were fixed in 4% paraformaldehyde prepared in FSW, buffered to pH 8.2 using 5 mM NaOH. Samples were photographed using a Zeiss Axio Scope A1 microscope equipped with a ProgRes CF Jenoptik camera and ProgRes Capture Pro software (v. 2.9.0.1).

| Bioinformatics analysis
All bioinformatics analyses were carried out using default parameters, unless otherwise specified. Illumina adapter trimming of the reads was performed using Trimmomatic v.0.33 (Bolger, Lohse, & Usadel, 2014), and the reads were further trimmed based on quality and length using Fastq-mcf v.1.04.636 (Aronesty, 2011), setting the Phred quality score to 30 and minimum read length to 60 bp.
A published mantle transcriptome of Baltic M. edulis-like individuals (Yarra, 2018, PRJNA494236), collected from the same geographic coordinates as the animals in this study, was used for mapping reads. The cleaned reads were aligned to the Baltic M. edulis-like mantle transcriptome (Yarra, 2018) using Bowtie v.1.1.1 (Langmead, Trapnell, Pop, & Salzberg, 2009), and the digital measure of transcript abundance was calculated using RSEM (RNA-Seq by Expectation- (Li & Dewey, 2011). All contigs with digital expression levels less than 2 counts per million (CPM) at the Trinity "gene" level, in at least half the libraries, were filtered out before analysis for differential expression. Preliminary analysis of the data revealed mislabeling of four samples, and the mislabeling was corrected as discussed in Appendix S1.
Contigs from the mantle transcriptome were annotated with a few different databases (Yarra, 2018). Sequence similarity searches of the transcript sequences were performed using BLAST (Altschul, Gish, Miller, Myers, & Lipman, 1990, blastx) with an E-value cut off of 1e −10 against public databases SwissProt (accessed 08 January, where at least 40% of the query sequence was aligned to, were considered to reflect strong sequence similarity. Transdecoder (part of the Trinity pipeline) was used to translate contigs into putative protein sequences of at least 20 codons. Translated protein sequences were mined for domain and family information using Interproscan (Jones et al., 2014), and Gene Ontology (GO) terms for contigs were assigned based on the Interpro database (Finn et al., 2017).
Differentially expressed contigs between developmental stages and treatments were identified using edgeR 3.20 (Robinson, McCarthy, & Smyth, 2010). contigs, and a cystic fibrosis transmembrane conductance regulator from zebrafish as an outlier. All accession ID's for protein sequences used in the tree are provided in Table S2.

| Larval development
Manipulation of seawater carbonate chemistry by the addition of 1 M HCl resulted in a reduction of C T from 1865.5 ± 26.2 µmol/ kg seawater under control conditions to 941.7 ± 51.3 µmol/kg seawater. In addition, a reduction in bicarbonate availability from 1840.8 ± 23.2 to 888.3 ± 47.5 and carbonate availability from 108 ± 3.7 to 43.7 ± 5 was observed (Table S1). Further, seawater C T reductions were associated with a decrease in pCO 2 from 423.4 ± 7.2 µatm under control conditions to 244.6 ± 23.7 µatm and Ω aragonite from 1.7 ± 0.03 to 0.6 ± 0.08. Development at reduced C T resulted in a developmental delay starting at 22 hpf, corresponding to the onset of calcification (Ramesh et al., 2017, Figure S1, Table 1). The mean developmental delay was 1.71 ± 1.38 hr, and in one family, the delay was observed to go up to 6 hr ( Table 1). The variability in developmental delay across replicates may be attributed to biological variability in larval energy budgets or differences in the degree of M. edulis × M trossulus hybridization between the four replicate families. Morphologically distinct developmental stages were ascribed to Stages 1-6 for further analyses.

| Quality control of sequencing reads
Sequencing of 48 larval libraries yielded a total of 590 million reads, with 541 million reads remaining after filtering based on quality and length. Cleaned reads were aligned to the Baltic M. edulis-like mantle transcriptome (Yarra, 2018), and mapping rates of approximately TA B L E 1 Morphological stages at which Baltic Mytilus edulis-like larvae were sampled during the experiment

Control families
Reared under low C T conditions Early trochophore Precalcification 20 20 20 20 20 Trochophore Larvae exhibit dorsal flattening at the region of the shell field which marks the onset of calcification

| Gene expression analysis
Three types of differential gene expression analysis were conducted for this dataset: pairwise comparisons between treatment and control groups at each developmental stage, pairwise comparisons between developmental stages, and time course comparisons between the developmental stages using an auto-regressive hidden Markov statistical model.
In the pairwise comparisons between treatment and control libraries, very few contigs (53) were found to be differentially expressed at each developmental stage (Table S3).
However, multiple contigs were found to be differentially expressed between the pairwise comparison of developmental stages. A large number of contigs (22,564) were differentially expressed between Stage 1 and each subsequent developmental stage, and multiple contigs were also found to be differentially expressed between each consecutive developmental stage (Table S4). Two hundred and forty-five contigs were differentially expressed following the onset of larval shell deposition (Stage 2) in comparison to the precalcifying ontogenetic stage (Stage 1; at 20 hpf, Table S5).
The time series-based differential gene expression analysis revealed two expression paths to be the most prevalent during PD I development in the Baltic M. edulis-like trochophore stage (

| Identification of transport pathways involved in calcification
The primary objective of this study was to identify candidates of ion transporter proteins potentially involved in providing substrates  induced a small set of contigs to be differentially regulated (Table 4,   (Table S3).

| Putative homology of bicarbonate ion transporters
There are multiple families of bicarbonate transporters present in eukaryotes (Alper & Sharma, 2013;Pushkin & Kurtz, 2006). To further characterize the sequence similarity of differentially expressed contigs with bicarbonate-transporting domains, a phylogenetic tree from multiple sequence alignments of translated protein sequences was assessed ( Figure 2). The phylogenetic analysis presented here indicates that like most eukaryotes, Baltic M. edulis-like larvae possess several cellular HCO 3 − -transporting proteins. Two putative mussel HCO 3 − -transporting proteins, which share high sequence similarity with membrane-bound SLC4 and SLC26 proteins, are the most interesting in the context of the present study ( Figure 3).
These sequences have been associated with calcification processes as they were found to be upregulated during adult mussel shell repair (Yarra, 2018).
In contrast to the small number of contigs exhibiting differential expression in response to substrate limitation, several contigs putatively encoding ion transport proteins corresponding to solute carrier families SLC4, SLC9, and SLC26 were differentially expressed during the course of larval development and shell deposition ( Figure 1). Among these SLC families, several contigs exhibited progressive increases in expression during the course of development (Table S3). These sequences encoded proteins such as sarco/endoplasmic reticulum Ca 2+ -ATPase, sodium/calcium exchangers (NCX), and the sodium/potassium ATPase. The putative ion transport pathways involved in larval calcification based on expression patterns for contigs of interest are presented schematically in Figure 3.

| Shell matrix proteins
Multiple genes that encode shell matrix proteins previously identified in the shell matrices of adult Mytilus spp. and expressed by the adult mantle tissue, particularly during shell repair, were found to be differentially expressed during larval shell development.
Approximately, 33% of the contigs annotated with SMP domains displayed an increasing expression profile starting from the trochophore stage (Table S6). A few shell matrix proteins (α-carbonic anhydrase, β-lactamase, concanavalin A, and cyclophilin PPIase) displayed decreasing expression levels as the initial shell was completed.

| D ISCUSS I ON
In this study, we employed a two-stage analysis. First, we used a

| Substrate limitation approach
We used low dissolved inorganic carbon, C T to challenge larval acidbase regulatory systems and identify candidate contigs putatively involved in calcification. In comparison to controls, substrate limitation treatments (mean C T = 941.7 µmol/kg seawater) were characterized by a ca. 51% and ca. 59% decrease in bicarbonate and carbonate availability, respectively. Such modifications in seawater carbonate chemistry were accompanied by a developmental delay in shell accretion that was consistent with previous empirical data ( Thomsen et al., 2015). These observations in developmental delay demonstrated that low C T conditions were correlated with decreases in larval calcification rates in Mytilus edulis (Thomsen et al., 2015 Figure S2), indicating that more subtle changes in gene expression were not detected with the present experimental design. Phylogenetic tree analyses revealed similarity between the M. edulis SLC26 contig and the human SLC26A11 sulfate/anion transporter ( Figure 3). Recently, the function of SLC26A11 transporters as sodium-independent sulfate transporters has been critically reviewed based on observations of their function as a chloride channel in mice neurons using electrophysiological techniques (Rahmati et al., 2013). Furthermore, transfection of the SLC26A11 transporter in mouse neuronal cells hints toward the activation of V-type H + -ATPases by SLC26A11 transporters inducing proton translocation (Rahmati et al., 2013). Alternatively, the upregulation of this trans-  (Schmidt et al., 2005). Interestingly, the function of transcription factor Sox 8 is strongly linked to the expression of runt-related transcription factor 2 (Runx2; Schmidt et al., 2005), where Runx2 (upregulated under low C T at Stages 3 and 5) is an important regulator of mammalian bone formation (Fowlkes et al., 2008;Franceschi & Xiao, 2002;Takarada et al., 2013) and arterial calcification (Ruffenach et al., 2016). Two other transcripts putatively encoding tumor necrosis factor α-inducing proteins and in-terleukin1 receptor-associated kinase 4 (IRAK4) that were observed to be upregulated under low C T conditions in M. edulis larvae have also been associated with osteoblast differentiation and bone resorption (Katsuyama et al., 2014;Tintut, Patel, Parhami, & Demer, 2000). Tumor necrosis factor α is involved in the activation of nuclear factor-kappa B (NF-κB, upregulated under low C T at Stage 5), a protein whose activation is linked to arterial calcification in humans (Zhao et al., 2012). Substrate limitation also induced upregulation of a contig with sequence similarity to a zinc transporter during Stage 5 of M. edulis larval development. In mammals, the knockout of zinc transporters has been demonstrated to result in skeletal disorders (Fukada et al., 2008) and the role of zinc is established in bone mineralization (Yamaguchi, 1998) and mollusc shell deposition (Tan & Mai, 2001). Finally, the expression of a contig encoding putative Ctype lectins was also upregulated under low C T in M. edulis larvae. Ctype lectins form important components of the shell organic matrix in molluscs (Mann, Edsinger-Gonzales, & Mann, 2012). For example, the shell matrix protein perlucin that is expressed in M. edulis larvae during PD I formation contains a C-type lectin domain. Therefore, our substrate (C T ) limitation technique elicited an expression change in several contigs that can be related to M. edulis larval development.
The low number of differentially expressed contigs with respect to substrate limitation was certainly linked to the high variability in gene expression, as observed before (Hüning et al., 2013;Yarra, 2018), but was also in line with previous studies on calcifying larvae that observe no significant changes in gene expression in response to simulated ocean acidification (Evans, Chan, Menge, & Hofmann, 2013; Kelly et al., 2016). In contrast, adult bivalves have been observed to exhibit differential regulation of genes related to ion and acid-base regulation in response to elevated seawater pCO 2 . Differences in larval and adult transcriptomic responses to CO 2 in other studies may reflect differences in acclimation pCO 2 , experimental design and individual variability. Alternatively, acid-base regulation may be ensued via posttranslational mechanisms, for example, the translocation of membrane-bound transport proteins to compensate for increased transport of calcification substrates (Roa, Munévar, & Tresguerres, 2014;Tresguerres, Parks, Wood, & Goss, 2007) and phosphorylation of ion transport proteins (Flemmer et al., 2010;Levitan, 1994) that can induce their rapid activation (Ramnanan & Storey, 2006;Rapoport & Murad, 1983
Among all the contigs putatively encoding ion transport and shell matrix proteins investigated in this study, a sequence encoding an NBC exhibited high transcript abundance during the larval development of M. edulis with a peak in expression during early calcification F I G U R E 2 An ion transport model depicting the cellular transport processes of membrane-bound ion transport proteins exhibiting elevated expression during Mytilus edulis-like larval development. The Na + /K + ATPase (NKA) provides the electrochemical gradient for secondary ion transport via proteins such as the Na + /H + exchanger (NHE, SLC9), sodium bicarbonate cotransporters (NBC, SLC4), and Na + / Ca 2+ exchangers (NCX, SLC8). Additionally, septate junctions may regulate the permeability of ions via the paracellular pathway (Jonusaite, Kelly, Donini, 2017). The putative precipitation of calcium carbonate in an extracellular calcification space is described in gray. The precise cellular location (apical/basolateral) and distribution of these ion transport proteins are unknown. In situ hybridization studies are necessary to ascertain the expression of these transporters calcifying epithelia F I G U R E 3 Phylogenetic tree depicting relationships between bicarbonate transporter families in Homo sapiens (human), Stylophora pistillata (Spi), Crassostrea gigas (Cgi), Strongylocentrotus purpuratus (Spu), and larval mussels (TRINITY). All sequences, along with accession IDs, are provided in Table S2. Starred sequences were differentially expressed in adult mussels during shell regeneration (Yarra, 2018) (Figure 1k, TRINITY_DN175808_c1_g5). The peak in expression of NBC encoding contigs is accompanied by the onset of shell formation. Following early trochophore development, expression levels of this NBC sequence rapidly decreased. Within the contigs encoding HCO 3 − transport in the transcriptome that exhibited differential expression during larval development, two contigs clustering with SLC4 (TRINITY_DN167998_c0_g2) and SLC26 (TRINITY_DN173725_c0_ g1) families were also observed to be upregulated during induced shell repair in adult Mytilus edulis, further supporting the role of these transcripts in substrate acquisition for calcification (Yarra, 2018).

| Ca 2+ transport
Prior to the onset of calcification at the trochophore larval stage, Ca 2+ is not accumulated and stored by mussel larvae (Ramesh et al., 2017). Rapid calcification of the PD I shell in mussels is accompanied by a tremendous uptake of calcium by larvae within a few hours (Ramesh et al., 2017). In contrast to larval sea urchins (Vidavsky et al., 2016;Vidavsky, Masic, Schertel, Weiner, & Addadi, 2015) the acquisition of the calcification substrates Ca 2+ and HCO 3 − from seawater via endocytotic transport does not seem to be a major pathway for calcium acquisition in larval mussels, suggesting that uptake of Ca 2+ likely occurs via transepithelial pathways (Ramesh et al., 2017 within the sarco/endoplasmic reticulum (Arruda & Hotamisligil, 2015). In adult bivalves, SERCA has been suggested to play a role in biomineralization due to its high expression (Truebano et al., 2010) and localized expression of one SERCA isoform (Fan et al., 2007) in mantle tissue. Aside from SERCA, organisms may also employ calcium-binding proteins to reduce free Ca 2+ concentrations intracellularly. Expression of contigs for one such calcium-binding protein, calbindin, was observed to increase during the course of M. edulis larval development with a peak in expression at PD I stage. Sodium/calcium exchangers (SLC8) is a group of membranebound transport proteins that facilitate the reversible exchange of three sodium ions (Na + ) for one calcium ion and has an established role in mammalian osteoblast (bone) calcification and avian eggshell mineralization (Cheidde, Viera, Lima, Saad, & Heilberg, 2003;Sosnoski & Gay, 2007

| Na + transport
The Na + /K + ATPase (NKA) protein is an active membrane-bound pump present on the basolateral membrane. It is critical for maintenance of cell membrane potential and generates the electrochemical gradient necessary to facilitate the subsequent transport of ions by secondary transport proteins (Boron & Boulpaep, 2009 (Figure 1u-y). In accordance with the upregulation of contigs encoding NKA during early shell formation, NKA activity has also been recorded to peak during early shell formation in oyster larvae (Frieder, Applebaum, Pan, Hedgecock, & Manahan, 2017). Similar transport processes are present in the primary mesenchymal cells in sea urchin larvae which are responsible for calcification and skeletogenesis where amiloride-sensitive ion transport proteins such as the NHE have been demonstrated to be significant for cellular pH regulation (Stumpp et al., 2012). In addition to elevated expression of contigs encoding NHEs at the onset of larval calcification, we also observed a peak in expression of gene encoding Na + /H + exchange regulatory factor (NHERF). Na + /H + exchange regulatory factor proteins are involved in regulating the function of NHE and have a pivotal role in bone formation, where their regulation of NHEs is crucial for osteoblast differentiation and strength (Liu et al., 2012). and the observed increases in pH at the site of calcification in larval mussels (Ramesh et al., 2017).

| Transport of other ions
Parallel to the acquisition of substrates for calcification (Ca 2+ and HCO 3 − ) and removal of proton by-products of calcification, there are ion transport proteins that are essential for maintenance of cellular electrogenic gradients, cell volume, etc. During the ontogenetic development of M. edulis, expression of several contigs encoding chloride (Cl − ) and potassium (K + ) channels was observed.
In particular, elevated expression for various chloride channels during PD I larval stage was detected ( Figure S3). Efflux of Cl − from the basolateral membrane due to acquisition of HCO Expression of several K + channels was seen in the larval transcriptome. However, significant changes in the expression of only one K + channel was observed during the course of development, the inward rectifier K + channel (K ir , Figure S3). Interestingly, this group of membrane-bound transport proteins was also found to be upregulated during shell repair in adult M. edulis (Yarra, 2018).
K ir channels are important in sustaining electrochemical gradients and cell resting potential by recycling K + ions (Weber, Cunningham, & Schulte, 2001). In addition, their absence in mammalian osteoblasts has been observed to inhibit osteoblastgenesis due to a decreased efficiency in production of an extracellular matrix (Sacco et al., 2015).

| Shell matrix proteins (SMPs)
It has long been known that SMPs play a critical role in calcium deposition and shell development (Weiner & Traub, 1984;Wheeler & Sikes, 1984), but these are poorly characterized in larvae.
Interestingly, approximately, 65% of the SMPs expressed in the adult Baltic M. edulis-like mantle transcriptome (Yarra, 2018) and ca. 65% of the SMPs extracted from adult shells were expressed during larval development. Similar to M. edulis-like larvae, three SMPs (nacrein, EGF-like, and tyrosinase) have also been observed to be expressed in other bivalve larvae (Fang et al., 2011;Liu et al., 2015). Multiple contigs with sequence similarity to other adult shell SMPs, but not yet identified in larval shells, were found to be differentially expressed throughout PD I development. Most of the differentially expressed contigs encoded domains involved in structuring the shell or tissue, such as β-hexosaminidase, glycoside hydrolase, chitin synthase, chitin binding, von Willebrand factor A, and Fibronectin type III. These were all found to have an increasing expression profile as the shell field was expanding over the surface of the larvae. Contigs containing copper-binding domains such as amine oxidase, dopamine-β-hydroxylase, and tyrosinase were also shown to increase in expression as the PD I shell was formed. Of particular interest were contigs containing tyrosinase domains that showed a decrease in expression right before the end of PD I formation. Tyrosinase proteins are involved in periostracum formation (Zhang, Xie, Huang, Chen, & Zhang, 2006), and the drop in expression of contigs containing tyrosinase domains may reflect the completion of periostracum formation in PD I.
Other SMPs linked to mineral deposition and crystallographic control include proteases and protease inhibitors that behave antagonistically, where the former are known to be important for crystal nucleation (Hershey et al., 2016;Tiaden et al., 2012). During Although bivalve larval shells are composed of only one calcium carbonate polymorph, aragonite (Kudo et al., 2010;Yokoo et al., 2011), several SMPs (tyrosinase, gigasin-like and alveoline-like proteins) that were previously identified in the calcitic fibrous prism structures of adult mytilid shells Liao et al., 2015) were also expressed in the Baltic M. edulis-like developing larvae.
This emphasizes the need for further investigation into larval stages over a longer developmental time period to determine whether there is specific SMP partitioning with development.

| CON CLUS IONS
By rearing mussel larvae under conditions of substrate limitation for calcification and analyzing differential gene expression patterns, we were able to identify a membrane-bound transport protein potentially involved in HCO 3 − acquisition, belonging to the SLC26 family of anion transporters and other candidate genes previously identified in human biomineralization. Interestingly, the present study identifies only a small subset of contigs to be differentially expressed under substrate limitation. Although this may be improved by a stronger C T treatment (>50% decrease relative to control), higher sequencing depth, and level of replication, bivalve larvae may possess a fixed capacity to modify their transcriptomic developmental program. This is consistent with previous studies on bivalve larval development that observed no significant changes in gene expression in response to induced acid-base stress (Kelly et al., 2016). Our data demonstrated an increased expression of contigs encoding for Ca 2+ and HCO 3 − -transporting proteins during larval development, in particular, once larval shell formation started.
In particular, the dynamic expression patterns and high expression levels of contigs encoding SERCA, NCX, and NBC hint toward the role of these ion transport pathways in bivalve larval calcification.
Similarly, the analyses of SMP expression patterns revealed several proteins with hypothesized roles in shell structure, crystallographic control, and periostracum deposition to be upregulated during larval development. To date, functional analyses using RNA interference techniques have been limited to SMPs (Fang et al., 2011;Funabara et al., 2014;Suzuki et al., 2009). However, knock out/knock down techniques are required to establish the role of candidate ion transporters in larval calcification. The identification of candidate biomineralization genes in this study paves the way for future in depth investigations. Heterologous expression techniques in conjunction with electrophysiology techniques should be used to characterize the substrate specificity and stoichiometry of these ion transport genes (in particular, the SLC26A11 ortholog) and their complex interactions with SMPs to produce a robust larval shell. We thank Ulrike Panknin and Nancy Kühne for help with the library preparation and sequencing.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
K.R. and F.M. designed the study. K.R. conducted the larval experiments, RNA extractions, and sequencing library preparations. T.Y conducted the bioinformatics analysis. K.R. and T.Y. analyzed the data and wrote the manuscript. All authors contributed to paper revisions.