The Mediterranean mussel Mytilus galloprovincialis: a novel model for developmental studies in mollusks

Historically, our understanding of the origin, mechanisms and evolution of animal development has been predominantly derived from the study of three groups: insects, nematodes and vertebrates (Müller and Grossniklaus, 2010). This resulted in a compendium of available tools, resources and experimental protocols for these organisms that has turned them into extremely amenable and efficient models for the scientific community (Müller and Grossniklaus, 2010; Nigon and Félix, 2017; Stephenson and Metcalfe, 2013). In comparison, the Spiralia, a species-rich and morphologically diverse clade of protostomes characterized by spiral cleavage during early development, have so far only been scantly investigated and still lack powerful experimental models (Martín-Durán and Marlétaz, 2020; Martín-Zamora et al., 2023a). As almost half of all bilaterian phyla are spiralians, members of this clade are of particular interest for understanding the evolution of bilaterian body plans (Henry, 2014; Martín-Durán and Marlétaz, 2020; Martín-Zamora et al., 2023a,b; Piovani et al., 2023). Similarly, morphological similarities between spiralian larvae have been used as taxonomic characteristics, and comparisons between larval types have yielded highly influential hypotheses on animal evolution (Carrier et al., 2018; Martín-Zamora et al., 2023b; Piovani et al., 2023; Wu et al., 2019). The trochophore larva, for example, is widely distributed amongst spiralians, and it has been suggested that it may have already been present in the last common ancestor of all bilaterian animals (Carrier et al., 2018; De Robertis and Tejeda-Muñoz, 2022; Martín-Zamora et al., 2023b; Nielsen, 2005, 2018). However, despite decades of studies, our understanding of the evolutionary origin of the trochophore larva, as well as of the biological processes underpinning its ontogeny and transition to the adult body, remains fragmentary (Carrier et al., 2018; De Robertis and Tejeda-Muñoz, 2022; Piovani et al., 2023; Wada et al., 2020; Wang et al., 2020; Xu et al., 2016).

Within spiralians, mollusks are the largest and anatomically most diverse phylum (Davison and Neiman, 2021; Yang et al., 2020). Yet only a limited number of species is currently available for studying the mechanisms of their development (Davison and Neiman, 2021; Goulding and Lambert, 2016; Lyons and Henry, 2022). This limitation is mainly due to the difficulty of spawning and rearing mollusks in a laboratory setting, and to the related complexity of setting up experimental protocols with only limited access to live material (Davison and Neiman, 2021; Yang et al., 2020). Several mollusk species have nonetheless been highlighted in terms of their potential as ‘models systems’ in diverse fields of biology, including bivalves of the order Mytilida (Davison and Neiman, 2021). Based on this observation, the current work provides an ensemble of updated protocols and new resources to significantly facilitate access to and work with embryos and larvae of a marine mytilid bivalve: the Mediterranean mussel Mytilus galloprovincialis.

M. galloprovincialis is an established model in environmental studies, and a species of great ecological and economic interest with a worldwide distribution (Brzozowska et al., 2012; Daguin and Borsa, 2000; FAO, 2022). Large-scale larval cultures of M. galloprovincialis are easily set up in the laboratory, by following standardized culture protocols from the International Standardization Organization (ISO) and the American Society for Testing Materials (ASTM) (ASTM E724-21, 2021; ISO 17244:2015, 2020). M. galloprovincialis is also the first metazoan with a published pangenome (Gerdol et al., 2020). The M. galloprovincialis genome is in fact characterized by accessory genomic regions encoding up to 20,000 dispensable genes subject to presence-absence variation (Gerdol et al., 2020). This genomic architecture resembles the open pan-genomes found, e.g. in prokaryotes and non-metazoan eukaryotes (Gerdol et al., 2020). This feature of the M. galloprovincialis genome is of significant importance for researchers interested in genome evolution. The presence of detailed culture protocols and of a well-annotated genome make M. galloprovincialis a promising model organism for developmental studies. In addition, basic protocols for pharmacological treatments, quantification of shell biogenesis, quantitative polymerase chain reaction (qPCR) as well as immunohistochemistry and in situ hybridization have recently been described for developing M. galloprovincialis embryos and larvae (Miglioli et al., 2019, 2021a,b).

In this study, we have compiled a detailed developmental transcriptome and have validated it with a standardized qPCR protocol. We further adapted, to M. galloprovincialis embryos and larvae, the previously described protocol for multi-color, fluorescent whole-mount in situ hybridization using hybridization chain reaction (HCR) (Choi et al., 2014; Piovani et al., 2023). We have successfully coupled this approach with immunohistochemistry. The developmental transcriptome of M. galloprovincialis has allowed us to define major developmental periods as well as comparisons with those of another bivalve mollusk: the Pacific oyster Crassostrea gigas. The novel in situ hybridization and immunohistochemistry protocols yield detailed information on the expression of genes marking developmental progression, e.g. in the early shell field, apical sensory cells, muscles, ciliated epithelia and peripheral mantle cells. Taken together, the results of this work represent a valuable resource for future work with M. galloprovincialis embryos and larvae, setting the stage for establishing this species as an emerging model system for studying the developmental biology of mollusks.

The developmental transcriptome of M. galloprovincialis embryos and larvae consists of 15 developmental stages (i.e. a total of 30 libraries from two independent biological replicates). Samples were collected every 4 h under standard culture conditions at 16°C (ASTM E724-21, 2021; ISO, 17244:2015, 2020) from the unfertilized egg to 52 h post-fertilization (hpf), with an additional sampling time point at 72 hpf. The stages sampled for transcriptomic analysis were defined by morphology and acetylated α-tubulin immunoreactivity (Fig. 1A, Fig. S1). As previously described in other mytilids, the M. galloprovincialis zygote undergoes holoblastic unequal cleavage (Kurita et al., 2009). The micromeres form at the animal pole of the embryo at 4 hpf (morula) and subsequent cleavage follows a spiral pattern (Kurita et al., 2009). By 8 hpf (blastula), the embryo reaches the blastula stage. By 12 hpf (gastrula), the macromeres at the vegetal pole start to invaginate, indicating the onset of gastrulation. At 16 hpf, the shell field (SF) gastrula, the stomodeum (or presumptive mouth) has moved anteriorly, and a rosette of dorsal ectodermal cells forms an invagination, indicating the onset of shell field ontogeny (Kniprath, 1980, 1981; Miglioli et al., 2019). At 20 hpf (trochophore 1), the columnar cells of the shell field thicken, creating a deeper invagination. At 24 hpf (trochophore 2), the larva has the stereotypical morphology of a trochophore, with a ventral stomodeum and a dorsal shell field that is still partially invaginated and secretes the organic matrix of the larval shell (Kniprath, 1980; Miglioli et al., 2019). By 28 hpf (trochophore 3) and 32 hpf (trochophore 4), the hinge region of the shell field is completely flat and the prototroch ciliary band surrounds the anterior hemisphere of the larval body. At 36 hpf (early veliger 1), a cavity appears on the dorsal side that connects the former stomodeum (that has now moved anteriorly) and the anal cavity on the posterior side. At 40 hpf (early veliger 2), the D-margin of the shell becomes clearly visible and only the velum (i.e. the ventralized prototroch) of the larva is protruding. By 44 hpf (D-veliger 1), the velum is enclosed in the D-shaped shell. Thereafter, no major morphological changes are observable, except for the progressive outgrowth of the shell with respect to the mantle tissue along the D-border. In summary, the samples used for the transcriptome analysis include representatives of all major developmental stages of early bivalve development.

More than 31 million paired-end clean reads per library were mapped on the published M. galloprovincialis genome (Gerdol et al., 2020), with an average input length of 200 bp (Table S1). More than 85% of the reads were alignable, with 71% of them assigned to a unique locus (Table S1). The output gene count matrix included 52,422 expressed genes (sum of gene raw counts in all libraries>0) out of a total of 60,302 genes in all libraries (Table S2). The minimum number of expressed genes (EGs) was found in unfertilized eggs (≈25,000 at 0 hpf), whereas the highest was found at the late D-veliger stage (≈40,000 at 72 hpf), indicating that at least 15,000 genes (corresponding to about 25% of the genome) are dynamically expressed during early development of M. galloprovincialis (Fig. 1B). The number of EGs increased drastically during the first 12 h of development (i.e. up to the gastrula), with the two replicates yielding similar results. Thereafter, a slow but steady increase in the number of EGs was observed until the last sampling point. Principal component analysis (PCA) of the most variably expressed genes validated sample reproducibility (Fig. 1C). The majority of the variance was explained by the first principal component (78%), which was clearly correlated with developmental timing, lending support to the notion that developmental progression is driven by transcriptomic changes.

Differential expression analysis was performed by pairwise comparison of samples of consecutive stages. The analysis identified 11,996 differentially expressed genes (DEGs) during M. galloprovincialis development (Table S3, Table S11). More than 1500 DEGs were already detected at the morula stage (4 hpf), indicating an early initiation of zygotic gene expression (Fig. 1D). The highest number of DEGs (≈5700) was reached at the blastula stage (8 hpf) and was followed by a progressive decrease as development proceeded. We next annotated the DEGs based on their functional domains (Fig. S2). DEGs were screened for transcription factor domains (Fig. S2A), for domains known to be involved in shell formation (Fig. S2B) and for signaling receptors and membrane transporters (Fig. S2C). In total, 711 transcription factors with 65 different domains were annotated, with the most recurrent (>20 genes) being zinc-finger, homeodomain, helix-loop-helix and fork-head domains (Fig. S2A, Table S4). With respect to shell formation, 779 genes were found to contain domains associated with shell matrix proteins and shell calcification (Fig. S2B) (Ramesh et al., 2019). The most recurrent domains included calcium-binding proteins (EF-hand and calponin), von Willebrand factors (WKA), extracellular matrix proteins (fibronectin), carbohydrate-binding proteins (lectin and concanavalin), chitin-binding proteins (CB), ion channels and transporters (Fig. S2B). Regarding receptors, 225 DEGs were annotated as G protein-coupled receptors of the rhodopsin and secretin families, 92 as ephrin A receptors and 23 belonged to the nuclear receptor superfamily, with other DEGs including scavenger, low density (LD) lipoprotein and TGFβII receptors (Fig. S2C). As many as 125 DEGs were predicted to be membrane transporter-associated domain proteins (Fig. S2C). Altogether, these results confirmed the relevance of the identified DEGs with respect to the major regulatory (i.e. transcription factors) and ontogenetic (e.g. shell matrix proteins) programs controlling early development of bivalve mollusks.

The differential expression analyses were corroborated by qPCR carried out on RNA samples from both biological replicates. Four target genes were selected based on their transcriptomic profiles and transcript abundance. Foxb2 (MGAL_10B093191) and Wnt8a (MGAL_10B085403) were selected because they are characterized by similar trends in their expression in the course of development (both peaking at the SF gastrula), but with significantly different mean expression values (high and low, respectively, for Foxb2 and Wnt8a) (Fig. 1E, Fig. S3). SYTL4 (MGAL_10B044489) and SLC46A3 (MGAL_10B020966) were selected because they have highly dynamic expression profiles (Fig. S3). For all four genes, the qPCR experiments validated the expression profiles obtained from the developmental transcriptome (Fig. 1E, Fig. S3). To allow a spatiotemporal mapping of gene expression, we adapted the hybridization chain reaction (HCR) protocol to M. galloprovincialis embryos and larvae (Fig. S4) (Choi et al., 2014; Piovani et al., 2023). The specificity of the HCR protocol was initially validated by comparison with colorimetric in situ hybridization of genes with known developmental expression patterns: Hox1, Tyrosinase and Tektin (Fig. S4A) (Miglioli et al., 2021b; Piovani et al., 2023). Thereafter, we assessed the efficiency of the HCR protocol using Foxb2 and Wnt8a (Fig. S4B). Although the expression of Foxb2 (higher baseline expression) was conspicuous at all surveyed developmental stages, that of Wnt8a (lower baseline expression) was less distinctive (Fig. S4B). A notable exception was the expression of Wnt8a in the SF gastrula (16 hpf), which is the developmental stage with the highest expression levels of both genes (Fig. S4B).

The HCR protocol proved extremely well-suited to co-expression analyses of developmental genes and thus allowed the parallel mapping of the expression patterns of Foxb2 and Wnt8a: whereas Foxb2 was detectable in the stomodeal and trochal area, Wnt8a was expressed in a discrete group of cells that will give rise to one of the posterior neural ganglions of the larva (Fig. 1F,G, Fig. S4B). These observations were consistent with the expression patterns of these two genes in embryos and larvae of other spiralians (Lartillot et al., 2002; Marlow et al., 2014; Tomer et al., 2010). To go even further, we coupled the HCR protocol with immunocytochemistry, using immunoreactivity against serotonin (5-HT) as a marker for the serotonergic system in the developing nervous system of mussel larvae (Miglioli et al., 2021b; Voronezhskaya et al., 2008). In addition to the identification of serotonergic neurons, co-labeling with serotonin also allowed a standardized orientation of the larvae, which significantly facilitated the interpretation of in situ hybridization signals in trochophore and D-veliger larvae of M. galloprovincialis (Fig. 1F,G). Taken together, our protocol coupling HCR in situ hybridization with immunocytochemistry allows high-resolution mapping of spatiotemporal gene expression dynamics in embryos and larvae of M. galloprovincialis.

Taking advantage of the developmental transcriptome and the HCR protocol, we next carried out a characterization of M. galloprovincialis early development at both the transcriptomic and morphological level. Hierarchical clustering of all transcriptome libraries was performed to identify statistically distinct developmental periods [adjusted-unbiased (au) P-values≥0.99] (Fig. 2A). The analysis divided the libraries in two main developmental clusters in line with developmental progression: one with embryonic libraries from egg to SF gastrula and one regrouping all the larval libraries from trochophore 1 to D-veliger 4 (Fig. 2A). Within the embryonic cluster, the egg and morula libraries were distinct from the blastula, gastrula and SF gastrula libraries, indicating that embryonic development can be subdivided, at the transcriptomic level, into an early and a late embryo period (Fig. 2A). Within the larval cluster, the trochophore libraries formed two sub-clusters, one established by trochophore 1 and 2, and one by trochophore 3 and 4, suggesting a split into an early and a mature trochophore period (Fig. 2A). Likewise, the early veliger and D-veliger libraries were separated into two sub-clusters, establishing the early veliger and D-veliger as two additional periods defining the early development of M. galloprovincialis (Fig. 2A).

To assess potential morphological differences in support of the developmental clusters identified in the transcriptome dataset, HCR in situ hybridization was performed on representative stages of each developmental cluster starting at the SF gastrula (16 hpf), the embryonic stage before the start of larval development (Fig. 2A,B, Fig. S5). Tektin (Tek), Myosin heavy chain (Mhc) and Tyrosinase (Tyr) were selected as tissue markers for, respectively, the ciliated epithelium, the muscular system and the shell field, as previously described for developing bivalve mollusks (Dyachuk and Odintsova, 2009; Miglioli et al., 2019; Piovani et al., 2023; Salamanca-Díaz et al., 2022; Yang et al., 2017). In SF gastrulae (16 hpf), Tek was expressed in the anterior and posterior tufts, with an inconspicuous signal detectable throughout the embryo. In trochophore 1 larvae (20 hpf), Tek expression was still localized in the apical tuft and anterior hemisphere. In trochophore 3 larvae (28 hpf), the Tek signal expanded to the prototroch ciliary band in the anterior hemisphere and to two islets of ciliated epithelium at the posterior pole, in close proximity to the stomodeum and future anus. In early veliger 1 larvae (36 hpf), Tek was detectable in the prototroch on the ventral side of the larva and, in D-veliger 1 larvae (44 hpf), Tek expression marked the velum as well as the oral and anal ciliary tufts, indicating that the prototroch was completely internalized at this stage of larval development (Fig. 2B,C, Fig. S1).

In SF gastrulae (16 hpf), Mhc was expressed in two groups of muscle cells, and this expression subsequently expanded to additional groups of muscle cells in trochophore 1 larvae (20 hpf). In trochophore 3 larvae (28 hpf), Mhc expression indicated a significant expansion of muscle tissue, with three pairs of symmetric larval retractor muscles developing along the prototroch and towards the posterior ciliated epithelium. The muscles were connected dorsally by longitudinal and transversal commissures. In early veliger 1 larvae (36 hpf), the Mhc signal became even more complex, as the larval retractor muscles reached the ciliated epithelium, the larval protractor muscle was formed and elongated towards the prototroch, and the anterior adductor muscle was clearly identifiable. In D-veliger 1 larvae (44 hpf), the expression of Mhc revealed that the anterior protractor muscles were associated with the ciliated epithelium and that the anterior retractors were connected by a muscle ring at the dorsal side of the velum (Fig. 2B,C, Fig. S1). This indicates that the larva has gained the ability to retract the body inside the shell (Dyachuk and Odintsova, 2009).

In SF gastrulae (16 hpf), Tyr expression was detectable in two islets of cells in close proximity to the invaginated shell field, and marked the evagination and expansion of the shell field in trochophore 1 larvae (20 hpf). In the trochophore 3 larva (28 hpf), Tyr expression was concentrated in islets of cells at the D-border of the larval shell, which formed a straight hinge, indicating that the evagination of the shell field was completed. In early veliger 1 larvae (36 hpf), Tyr expression suggested that the larval shell was covering most of the larval body and, in D-veliger 1 larvae (44 hpf), expression of Tyr at the periphery of the larva indicated that the shell was enclosing the complete larval body (Fig. 2B,C, Fig. S1). Taken together, these observations support the developmental staging obtained by hierarchical clustering, as each developmental period was associated with marked changes in larval anatomy: early trochophores were characterized by an invaginating shell field and the differentiation of discrete islets of muscle cells; mature trochophores by a completely evaginated shell field, a differentiated prototroch and developing retractor muscles; early veligers by the shell covering most of the larval body, a ventralized prototroch, and anterior adductor and protractor muscles; and D-veligers by a shell that completely encloses the larval body and muscles that can retract the ciliated epithelium within the larval body (Fig. 2C).

Gene soft clustering was performed on the whole transcriptome to identify potential gene sets and biological processes enriched at specific developmental stages (Fig. 3, Fig. S6). The analysis identified nine gene clusters that, with the exception of early veliger 1 and 2 stages, closely matched the developmental stage clusters defined above (Fig. 3, Table S5). Early embryonic development was subdivided into a maternal and a zygotic gene cluster. For the maternal cluster (i.e. the egg stage), enriched gene ontology terms were associated with protein phosphorylation, microtubule-based movement and signal transduction (Fig. 3, Fig. S6). The zygotic cluster was enriched in genes involved in RNA splicing and DNA replication (Fig. 3, Fig. S6). Late embryonic development, i.e. the blastula and gastrula gene clusters, were both enriched in DNA recombination genes (Fig. 3, Fig. S6). The gene cluster corresponding to the early trochophore was characterized by genes associated with intracellular protein transport, RNA processing and regulation of transcription, whereas the mature trochophore gene cluster was associated with genes involved in signal transduction and apoptosis (Fig. 3, Fig. S6). The gene clusters defining the D-veliger, the mature D-veliger and the late D-veliger were generally enriched in gene ontology terms associated with the G-protein-coupled receptor signaling pathway and transmembrane transport, in addition to translation (for the mature D-veliger) and carbohydrate metabolism (for the late D-veliger) (Fig. 3, Fig. S6). In summary, this analysis corroborates the notion that the early ontogeny of M. galloprovincialis includes six main developmental periods (early embryo, late embryo, early trochophore, mature trochophore, early veliger and D-veliger), each characterized by a distinctive set of anatomical features and transcriptomic signatures.

To assess whether similar transcriptomic dynamics occur during the early development of other bivalve mollusks, gene soft clustering was performed on the developmental transcriptome of the Pacific oyster C. gigas (Liu et al., 2021) and compared with that of M. galloprovincialis (Fig. 4). Nine of the ten soft clusters identified in the C. gigas transcriptome showed similar transcriptomic dynamics to those observed in M. galloprovincialis, covering the main developmental periods of gene expression, from the egg to the late D-veliger stage (Fig. 4A, Table S6). Like our analysis in M. galloprovincialis, that of the C. gigas developmental transcriptome failed to identify a gene cluster for the early veliger libraries. Comparison of the orthologous gene clusters in M. galloprovincialis and C. gigas showed that the D-veliger (40%), late D-veliger (32%), mature D-veliger (31%), zygotic (31%), early trochophore (28%) and mature trochophore (27%) clusters shared the highest percentage of orthologous genes, with the maternal (13%), blastula (11%) and gastrula (9%) clusters being characterized by much lower percentages of orthologous genes (Fig. 4B, Table S7). The orthologous gene sets included several transcription factors (such as fork-head, homeobox and nuclear receptors) and signaling components (of the Wnt, TGF and EGF pathways, for example), suggesting that the embryonic and early larval development of the two bivalves might rely on conserved gene regulatory programs (Table S7). However, only 17% of the M. galloprovincialis genes had an ortholog in C. gigas (i.e. 10,909 out of about 65,000 genes). This might be explained by the overall number of genes annotated in the genomes of the two species, which is significantly higher in M. galloprovincialis than in C. gigas (∼60,000 genes versus ∼30,000 genes) (Gerdol et al., 2020; Peñaloza et al., 2021), but nonetheless implies that the comparative analysis by gene orthology could only be applied to a subset of the genes characterizing the early development of M. galloprovincialis. Taken together, these results suggest that the temporal dynamics of gene expression during early development of M. galloprovincialis and C. gigas are relatively similar. However, the low percentages of orthologous genes between most of the corresponding developmental clusters of the two bivalve species indicates that there might be some degree of divergence in the genetic programs controlling early development of M. galloprovincialis and C. gigas. To shed additional light on conserved aspects of early M. galloprovincialis and C. gigas development, we identified representative marker genes for the early trochophore, mature trochophore and D-veliger periods. Genes were selected based on membership of the same developmental clusters in the two bivalve species and upregulation in the respective mussel libraries. The criteria were matched by 48, 11 and 6 genes for, respectively, the early trochophore, the mature trochophore and the D-veliger (Table S7). We then selected two genes for each developmental period for HCR in situ hybridization (Fig. 4C): Dmbt1 and Dimm for the early trochophore, Cbp2 and NR3A for the mature trochophore, and ColVI-like and EF-hand for the D-veliger (Table 1). Dmbt1 and Dimm were both expressed in the shell field of M. galloprovincialis at the early trochophore stage, Cbp2 and NR3A, respectively, along the hinge axis and throughout the larval shell of the mature trochophore, and ColVI-like and EF-hand, respectively, along the D-border of the shell and in the ciliated epithelium of the velum at the D-veliger stage (Fig. 5). The expression of five of the six genes in tissues secreting the developing shell strongly indicates that the genetic program regulating larval shell formation might be, to some degree, conserved between the two bivalve species.

The results presented in this study provide the first thorough characterization of the early development of M. galloprovincialis at both a transcriptomic and a morphological level. The sampling schedule at 4 h intervals allowed us to closely map the gene expression dynamics underpinning developmental progression. Analysis of the functional domains of DEGs highlighted the power of our transcriptome to delineate early bivalve development. Transcription factors known to act in the regulation of early metazoan development, such as homeobox and fork-head genes, and cohorts of genes involved in early shell formation were thus highly represented among these DEGs (Guo et al., 2023; Lartillot et al., 2002; Morino et al., 2013; Paps et al., 2015; Ramesh et al., 2019; Setiamarga et al., 2013; Zhao et al., 2018; Zheng et al., 2019). In addition, the analysis revealed a strong enrichment of receptors and transporters, which is consistent with previous analyses in bivalve mollusks, indicating that these mediators of intercellular signaling might act in the regulation of larval morphogenesis, neurogenesis and shell biogenesis (Guo et al., 2023; Liu et al., 2020; Marin et al., 2013; Miglioli, 2022; Miglioli et al., 2021a; Song et al., 2016).

To correlate the transcriptional and anatomical changes marking early M. galloprovincialis development, we set up a high-resolution, whole-mount protocol for multicolor in situ hybridization that was compatible with immunostaining approaches. In combination with the developmental transcriptome, this protocol allowed the parallel mapping of tissue marker genes previously described for other bivalve species in M. galloprovincialis embryos and larvae: Tektin (Tek) for the ciliated epithelium; Myosin heavy chain (Mhc) for the muscular system; and Tyrosinase (Tyr) for the tissue forming the shell (Dyachuk and Odintsova, 2009; Miglioli et al., 2019; Piovani et al., 2023; Salamanca-Díaz et al., 2022; Yang et al., 2017; Zhang et al., 2006). We were thus able to sketch the morphogenesis of M. galloprovincialis embryos and larvae, demonstrating that developmental progression at the transcriptomic level is closely matched by anatomical changes, at least at the level of ciliated epithelia, musculature and shell growth. Although not novel per se, the combination of transcriptomic data and fluorescent whole-mount staining protocols has so far found only very limited use in surveys of mollusk development (Krasity et al., 2015; Lopez-Anido et al., 2023 preprint; Piovani et al., 2023). The datasets and protocols presented here therefore represent significant advances that will facilitate future studies on bivalve mollusk development and support the establishment of M. galloprovincialis as a model system in developmental biology.

Hierarchical clustering of libraries and gene soft clustering revealed that the early development of M. galloprovincialis is divided into six different periods, each characterized by distinctive gene expression dynamics (Futschik and Carlisle, 2005; Leclère et al., 2019; Suzuki and Shimodaira, 2006). To assess whether these developmental gene expression dynamics are conserved in other mollusks, we performed gene clustering analyses on the developmental transcriptome of the Pacific oyster C. gigas. Comparisons of the resulting C. gigas gene clusters with those obtained in M. galloprovincialis revealed similar temporal dynamics of gene expression that were, however, not strictly based on expression of orthologous gene sets. This is in line with previous findings indicating that temporal shifts in the use of common gene sets are common in bilaterian larvae and that this process might have facilitated the evolutionary diversification of larval forms and life cycles (Martín-Zamora et al., 2023b). Furthermore, trochophore development in both species was characterized by an early and a mature transcriptomic signature, which we were able to associate, in M. galloprovincialis, with distinctive morphological features and thus specific larval morphotypes. The fact that we were not able to associate the early veliger libraries of both species with specific gene clusters might suggest that this developmental period represents a transition phase from the transcriptomic programs of the trochophore to that of the D-veliger. Taken together, our work identified similar temporal gene expression dynamics in the two bivalve mollusk species that, albeit not driven by strictly orthologous gene sets, corroborated the existence of six conserved developmental periods during bivalve mollusk ontogeny, each characterized by distinctive transcriptomic and anatomical characters: early embryo, late embryo, early trochophore, mature trochophore, early veliger and D-veliger. These results suggest that the staging of bivalve larval development should be carried out with great care to avoid generalized stage assignments that could lead to erroneous comparisons.

In this work, early development of M. galloprovincialis was morphologically characterized using expression of the tissue-specific marker genes Tektin (Tek), Myosin heavy chain (Mhc) and Tyrosinase (Tyr), respectively labeling the ciliated epithelium, musculature and shell-forming tissue. Tek was selected as a marker of ciliated cells based on the single-cell atlas of the Pacific oyster C. gigas (Piovani et al., 2023). Mhc has previously been used to characterize the muscle system in M. coruscus (Dyachuk and Odintsova, 2009) and is expressed in muscles of the quagga mussel Dreissena rostriformis (Salamanca-Díaz et al., 2022). Tyr expression is known to be restricted to the outer edge of growing mantle tissue in both oysters and mussels (Miglioli et al., 2019; Yang et al., 2017). Co-labelling of developing M. galloprovincialis with these three marker genes allowed us to define key morphological events in each larval developmental period and to highlight the main anatomical differences between them. These results suggest that gene expression surveys using Tek, Mhc and Tyr are an effective tool for assessing developmental progression of bivalve mollusks, one that could be generalized to homogenize ontogenetic staging systems within this animal clade.

Detailed comparisons between the developmental transcriptomes of M. galloprovincialis and C. gigas identified sets of highly conserved genes defining the bivalve early trochophore, mature trochophore and D-veliger periods. Of these, we found that two early trochophore genes, Dmbt1 and Dimm, were expressed in the developing shell field. Both genes have previously been described as regulators of terminal differentiation, Dmbt1 in columnar cells and Dimm in neurosecretory cells (Liu et al., 2016; Takito and Al-Awqati, 2004). Expression of the two mature trochophore genes Cbp2 (a chitin binding peritrophin) and NR3A (a nuclear receptor) was also associated with the developing shell, so was that of the D-veliger gene ColVI-like, which was detectable along the edge of the shell-forming tissue. Only the D-veliger gene EF-hand was not expressed in developing shell tissues but was found in the ciliated epithelium instead.

The expression dynamics of the conserved genes identified in our comparative analysis closely follow the stereotypical phases of shell formation in bivalve mollusks: terminal differentiation of shell field cells in the early trochophore, organic matrix growth in the mature trochophore and calcification in the D-veliger (Kniprath, 1980, 1981). Cbp2 proteins, for example, are required for the assembly of the organic layer of the larval shell in the mature trochophore, while ColVI-like proteins are necessary for the successive formation of the nacre in the D-veliger (Marin et al., 2013; Zhao et al., 2018; Zheng et al., 2017). These results strongly suggest that the genes identified by our comparative analysis are part of the molecular toolbox for biomineralization in bivalve mollusks (Davison and Neiman, 2021; Marin et al., 2013; McDougall and Degnan, 2018) and indicate that the gene networks mediating shell formation are among the most conserved genetic components of early bivalve mollusk development. In addition, the implication of some of these genes (Dimm and NR3A) in neurogenesis in ecdysozoans and deuterostomes (Coumailleau et al., 2015; Liu et al., 2016) lends support to the notion that genes involved in nervous system development might have been co-opted for shell formation during the early evolutionary diversification of mollusks (McDougall and Degnan, 2018; Wollesen et al., 2017).

The results reported in this study represent the first thorough transcriptomic analysis of the development of the Mediterranean mussel M. galloprovincialis. This work further provides a collection of protocols, for qPCR, immunohistochemistry and HCR in situ hybridization, that are essential for the validation and exploitation of transcriptomic data. We found that the combination of fluorescent immunohistochemistry and HCR in situ hybridization is particularly powerful for assessing the colocalization of developmental factors in developing mussel embryos and larvae. Potential applications of this transcriptome, and associated experimental protocols, include developmental studies, cross-species comparisons and ecotoxicological surveys. Taking advantage of the developmental transcriptome and the novel protocols, we are the first to document the expression of a suite of genes illustrating the developmental progression of M. galloprovincialis, hence defining a set of markers that can be used in other bivalve mollusk systems to characterize early development. We are confident that this represents the first step for M. galloprovincialis to become a reference for developmental studies of bivalve mollusks. The next steps need to include the establishment of protocols for transgenesis to allow functional analyses, which are currently limited to pharmacological approaches (Davison and Neiman, 2021). Tackling this challenge should be a priority for the scientific community to further improve the experimental amenability of this promising animal model.

Sexually mature M. galloprovincialis adults were harvested from the natural population in the bay of Villefranche-sur-Mer, France (43.682°N, 7.319°E) during the spawning seasons of 2021, 2022 and 2023 (January through April). Animals were acclimated to and maintained in laboratory conditions, as previously described (Miglioli et al., 2021a). To induce gamete release, 30 adult mussels were exposed to a heat shock. Briefly, epibionts, dirt and algae were scrubbed off the shells under running filtered seawater. Byssus and threads were carefully cut off. Clean mussels were placed on ice for 10 min and then immersed, in individual 200 ml containers and on a rocking shaker, in filtered seawater kept at 28°C until gamete release. Spawning mussels were immediately moved to filtered seawater at 16°C. Before fertilization, gametes were collected, washed and maintained in Millipore-filtered seawater (using a 0.2 µm Millipore filter). Fertilization was carried out as previously described with a 1:10 ratio of eggs to sperm (Miglioli et al., 2021a). After 30 min, the fertilization success rate (i.e. the ratio of the number of fertilized eggs over the total number eggs×100) was established by microscopic observation. Parental pairs with a high fertilization success rate (>90%) were selected for the following experiments. Two experimental replicates were set up, each consisting of a homogeneous pool of fertilized eggs from two independent parental pairs. Embryos from each replicate were distributed into 200 ml flasks for suspended cultures and brought to a density of 100 embryos per ml with Millipore-filtered seawater. A total of 14 flasks per replicate were prepared, corresponding to the pre-defined sampling strategy: in addition to sampling before fertilization, sample collections were carried out every 4 h from 4 hpf to 52 hpf, with an additional sampling timepoint at 72 hpf. Embryo cultures were maintained at 16°C throughout the experiment (ASTM E724-21, 2021; ISO, 17244:2015, 2020; Miglioli et al., 2021b). Samples for immunohistochemistry and in situ hybridization were produced following the same culture protocol. After filtration, larvae were fixed in 4% paraformaldehyde in 1×PBS overnight at 4°C, washed three times for 15 min in 1×PBS with 0.01% Tween 20, and then stored in 100% methanol at −20°C (Miglioli et al., 2021b).

RNA extraction and purification were carried out with the RNAqueous-Micro Kit (Invitrogen) and the RNeasy MinElute Cleanup Kit (Qiagen) following published protocols (Leclère et al., 2019). Unfertilized eggs were directly collected from spawning females, and embryos and larvae were sampled with a 60 µm mesh filter and concentrated in RNAse-free 1.5 ml tubes. The tubes were then placed on ice and gently centrifuged multiple times to remove the seawater from the samples. The pellet was immediately resuspended in lysis buffer, homogenized by energetic pipetting and flash-frozen in liquid nitrogen. For each sample, 10 µl of the 1.5 ml filtrate were placed in a separate tube and fixed overnight in 4% paraformaldehyde in 1× PBS at 4°C (Miglioli et al., 2021b). The fixed larvae were stained with Hoechst (Invitrogen) and imaged with a SP8 confocal microscope (Leica Microsystems) to validate proper developmental progression (Miglioli et al., 2021b). Lysed samples were stored at −80°C until RNA extraction and purification (Leclère et al., 2019). RNA concentration and quality was assessed with a Nanodrop (Thermo Scientific) and a 2100 Bioanalyzer (Agilent). Samples were shipped on dry ice to BGI Genomics (Hong Kong) for high-throughput sequencing. At least thirty million clean paired-end reads of 200 bp each were produced using the BGISEQ-500 Platform. Read cleaning and trimming was performed by BGI Genomics. After sequencing, the raw reads were filtered. Data filtering included removal of adapter sequences, contaminations and low-quality reads from the raw data. The entire read was deleted if more than 28% matched the adapter sequence, if more than 40% had a quality value lower than 20 or if there were more than 3% of unidentifiable nucleotides in the read. This step was completed by the SOAP nuke software developed by BGI  Genomics. Libraries were not subjected to further trimming (Dobin et al., 2013).

Sequencing quality of the libraries was assessed with FastQC (v. 0.11.7) (Andrews, 2010). The libraries were mapped on the indexed M. galloprovincialis reference genome with STAR (v. 2.7.10a) using default parameters (Dobin et al., 2013; Gerdol et al., 2020). Read counts were obtained using StringTie (v.2.2.1) in counting mode with the genome annotation file (Pertea et al., 2015, 2016). The resulting gene count matrix was exported using prepDE.py and processed in Rstudio (v. 4.2.2) with the DESeq2 (v. 1.24.0) package using default parameters and a design file, indicating replicate number and developmental time for each library (Love et al., 2014). The dataset was then pre-filtered by keeping all genes with at least 10 reads in at least two replicated libraries to reduce the potential batch effect derived from the variability in gene presence-absence characterizing the species (Gerdol et al., 2020). The resulting matrix was normalized using the regularized logarithms of counts (DESeq2::rlog) and used to perform principal component analyses and heatmaps (Kolde, 2019). Stage-by-stage differential expression analysis was performed with DESeq2 using the pairwise comparison function, with the reference stage for comparisons re-levelled for each contrast (Love et al., 2014). Pairwise comparisons were performed on samples of consecutive stages (0 hpf versus 4 hpf, 4 hpf versus 8 hpf, etc.). DEGs were filtered by adjusted P-value (Padj.cutoff<0.05) and log2 fold change (lfc.cutoff>0.58).

To discriminate stages of M. galloprovincialis development from a transcriptomic perspective, bootstrapped hierarchical clustering with pvclust (v. 2.2-0) default parameters was performed on regularized logarithms (rlog) of the DESeq2 count matrix using the mean expression values of replicated libraries (Leclère et al., 2019; Love et al., 2014; Suzuki and Shimodaira, 2006). Soft clustering with the Mfuzz (v. 2.56.0) package was performed on mean expression values of replicated libraries normalized by transcripts per million to identify representative gene clusters of different developmental phases and potential transcriptomic waves (Futschik and Carlisle, 2005). The number of soft clusters was determined using the elbow method to the minimum centroid distance (Futschik and Carlisle, 2005). Genes were assigned to the cluster with the highest membership probability.

The predicted protein sequences of the DEGs were annotated using BLASTp (v. 2.11.0) against the UniProt/SwissProt database, conserved domains were identified with hmmerscan from HMMER 3.3 against the Pfam database and transcription factor domains were screened with AnimalTFDB (v2.0) (Finn et al., 2015, 2017; Jones et al., 2014; Mahram and Herbordt, 2015; Mistry et al., 2021; Sayers et al., 2021; Shen et al., 2023). Ontology enrichment analysis was performed for each gene soft cluster. The protein sequence corresponding to genes of interest were extracted from the M. galloprovincialis proteome and annotated with Gene Ontology (GO) identifiers using Interproscan (v. 95.0) (Finn et al., 2017; Gerdol et al., 2020; Jones et al., 2014). Enrichment analysis of Biological Processes (BP) was conducted with topGO R (v.3.17). The significance of the enrichment was computed with Fisher's exact test (Alexa and Rahnenführer, 2009).

The transcripts per million matrix of the C. gigas developmental transcriptome from the egg to the D-veliger stage was downloaded from MolluscDB and gene soft clustering was performed as previously described (Futschik and Carlisle, 2005; Liu et al., 2021). Orthologous genes were defined with DIAMOND (v.0.9.22) using the rbhXpress script (Buchfink et al., 2014; El Hilali and Copley, 2023 preprint). Representative genes for different periods of early development conserved between C. gigas and M. galloprovincialis were selected based on membership to corresponding gene clusters and differential expression dynamics in the respective M. galloprovincialis transcriptome libraries.

Complementary first strand DNAs (cDNAs) were synthesized from total RNA extracts of both biological replicates at all sample timepoints using the SuperScript VILO Kit (Invitrogen). The synthesis of cDNA was performed using random primers following the manufacturer's instructions: 10 min at 25°C, 60 min at 42°C and 5 min at 85°C. To validate the differential gene expression analysis and the cluster analysis, four target genes were selected: Foxb2 (MGAL_10B093191), Wnt8a (MGAL_10B085403), SYTL4 (MGAL_10B044489) and SLC46A3 (MGAL_10B020966). EF-α1 was used for data normalization as one of the best-performing reference gene controls (Balbi et al., 2016). Primer pairs were designed within the open reading frame of each gene using Primer3 (Koressaar and Remm, 2007). Sequences, efficiency and amplicon sizes are available in Table S8. The qPCR experiments were carried out for each set of primers, in experimental triplicates. The amplification was followed in real time using LightCycler 480 SYBR Green 1 master mix (Roche) in 45 cycles (initial 5 min at 95°C, followed by 45 cycles of 20 s at 95°C, 20 s at 60°C and 30 s at 72°C) in a LightCycler 480 (Roche). Quantification was carried out in parallel for each of the selected genes and the reference gene. Quantification of expression at a given stage was determined in relation to the expression of the reference gene using the comparative CT method (Schmittgen and Livak, 2008). To compare qPCR and transcriptome results, expression of the target genes from both datasets were normalized with respect to their maximum value. Statistical differences in qPCR results were assessed using an unpaired Student's t-test (de Winter, 2013).

For hybridization chain reaction (HCR)-based in situ hybridization experiments, hybridization probe sets were designed using a user-friendly Python interface (https://github.com/rwnull/insitu_probe_generator) and synthesized by OligoPools (Twist Bioscience). Amplifiers (B2: 546 nm, B1: 647 nm, B3: 514 nm) were purchased from Molecular Instruments. Probe sets with relative amplifier pairings are summarized in Table S9. Larval samples stored in methanol were rehydrated with four washes of 5 min each followed by a 10 min incubation in 1×PBS with 0.01% Triton X-100 at room temperature. Permeabilization was carried out by a 30-min incubation in a detergent solution [1% SDS, 0.005% Tween 20, 50 mM HCl pH7.5, 1 mM EDTA (pH8.0) and 0.15 M NaCl] at room temperature. After two washes of 1 min each and two washes of 5 min each in 1×PBS with 0.01% Triton X-100 at room temperature, the specimens were placed in pre-hybridization for 30 min at 37°C with gentle shaking in hybridization buffer (HB) [30% formamide, 5×SSC, 9 mM citric acid (pH 6.0), 0.01% Tween 20, 2 mg heparin, 1×Denhardt's solution and 10% dextran sulphate]. After addition of the probes (at a final concentration of 10 nM), hybridization was carried out overnight at 37°C. After hybridization, 1.5 volumes of washing solution [30% formamide, 5×SSC, 9 mM citric acid (pH6.0), 0.01% Tween 20 and 2 mg heparin] with respect to the volume of the hybridization solution were added to the samples, and four washes of 15 min each were subsequently performed with washing solution at 37°C. Specimens were then washed, at room temperature, twice for 5 min in 5×SSC with 0.01% Tween 20 and once for 5 min in 1×SSC with 0.01% Tween 20, and then incubated for 30 min in amplification buffer (5×SSC, 0.01% Tween 20 and 10% dextran sulphate) before amplifiers were added (at a final concentration of 60 nM) to start the amplification chain reaction. Amplifiers were prepared as follows: the required volume of hairpins from each amplifier were placed in a 0.6 ml tube and incubated for 90 s at 95°C. After the heat shock, amplifiers were placed in the dark at room temperature for 30 min. The amplifiers were then resuspended in amplification buffer and distributed to the samples. Amplification was carried out at room temperature in the dark overnight. Excess hairpins were removed in two washes of 5 min each and two washes of 30 min each in 1×SSC with 0.01% Tween 20, followed by two washes of 5 min in 1×PBS with 0.01% Tween 20.

Co-labeling by immunocytochemistry with anti-serotonin (5-HT) rabbit antibody (1:10,000; 20080, Immunostar) and goat anti-rabbit Alexa Fluor 488-conjugated secondary antibody (1:500; 4412S, Cell Signaling Technology) was performed immediately after HCR amplification using a previously described protocol (Miglioli et al., 2021b). The same protocol was used to carry out labelling with anti-acetylated α-tubulin mouse antibody (1:10,000; T6793, Sigma-Aldrich) and goat anti-mouse Alexa Fluor 488-conjugated secondary antibody (1:500; A-11001, Thermo Fisher Scientific). Specimens were then prepared for imaging by incubation in 1 µg/l Hoechst (Invitrogen) in 1×PBS for 30 min at room temperature. After three washes of 5 min in 1×PBS with 0.01% Tween 20, specimens were mounted with the antifading agent CitiFluor AF1 (Agar Scientific) and imaged on a SP8 confocal microscope (Leica Microsystems), as previously described (Miglioli et al., 2021a). Images were analyzed with ImageJ (Schindelin et al., 2012) and fluorescent signals were corrected as follows: autofluorescence from tissue and shell was subtracted with ‘math’, the signal was homogenized with the ‘despeckle’ or ‘smooth’ functions and z-stacks were merged by either ‘maximum intensity’ or ‘sum slices’. Fluorescent HCR in situ hybridization signals were validated with colorimetric in situ hybridization experiments using genes with known developmental expression as positive controls and with no-probe HCR in situ hybridization assays as negative controls (Table S10). Colorimetric in situ hybridization was performed as previously described for M. galloprovincialis embryos and larvae (Miglioli et al., 2019).

The authors are indebted to Jenifer Croce, Guy Lhomond and Lucas Leclère for technical assistance and fruitful discussions. The authors thank Sébastien Schaub of the Plateforme d'Imagerie par Microscopie (PIM) as well as Laurent Gilletta and Axel Duchene of the Service Aquariologie (SA) of the Institut de la Mer de Villefranche (IMEV), France, which is supported by EMBRC-France and funded by the Agence Nationale de la Recherche (ANR) (ANR-10-INBS-02). We thank the members of the Ascidian BioCell and EvoInSiDe teams at the Laboratoire de Biologie du Développement de Villefranche-sur-Mer, France, as well as the members of the team of Environmental Physiology at the Università Degli Studi Di Genova, Italy, for their support of this work.

The peer review history is available online at https://journals.biologists.com/dev/lookup/doi/10.1242/dev.202256.reviewer-comments.pdf