Single cell RNA sequencing data from Dreissena rostriformis that had previously been generated (Salamanca-Díaz et al., 2022) were used for the assessment of unknown genes expressed in the shell field as well as for further analyses. In the following, a summary of all major steps from animal acquisition through the in silico analyses performed herein is provided.

Sexually mature individuals of Dreissena rostriformis were collected from the Danube River in Vienna, Austria (N 48°14′45.812″, O 16°23′38.145″). Collection took place between April and September 2019. Adults were gathered from underneath stones and transferred to the laboratory where they were cleaned and maintained in aquaria with filtered river water (FRW) at 19°C.

Spawning of animals was induced by incubating sexually mature specimens in a 10⁻³ M solution of serotonin for 15 min (Sigma-Aldrich, Darmstadt, Germany) in FRW, followed by one wash and subsequent maintenance in FRW. Individuals were kept isolated in FRW in 50 ml glass beakers and after approximately 30 min, up to 50% of the treated specimens started to spawn. Fertilization occurred when three to four drops of sperm-containing water were added to 50 ml glass beakers with oocytes. After fertilization, water was changed every half an hour for the first 3 h and then every 6 h to remove excess sperm and avoid bacterial and fungal growth. The embryos were cultured at 23°C.

To discard false positives while screening for novel genes in the shell field, the orthology assessment was made using the whole genome of Dreissena rostriformis. After that, we analyzed the genes that are exclusively part of the transcriptomic signature from the shell field of the trochophore larva. Around one third of all orthogroups (31.4%) predicted from 30 different metazoan species contain Dreissena rostriformis (“DRERO”) genes (cf. Supplementary Table S3). In addition, we identified D. rostriformis lineage-specific orthogroups with non-annotated genes, meaning there is a noteworthy number of genes that have no known match with any other animal sampled. However, all other species used in our analysis show similar low percentages of genes that can be assigned to known orthogroups (Supplementary Table S3), corroborating the common notion of the vital role of lineage specific genes or families during animal genome evolution (Fernández and Gabaldón, 2020). In D. rostriformis, such genes identified from the genome mount up to 19.8% (7469 genes; see Supplementary Tables S2, S3).

Using the outputs from the OrthoFinder and Single-cell seq pipelines, we characterized the shell field-specific genes and their orthogroup correspondence (Supplementary Table S2). In total, we analyzed 357 genes differentially expressed in shell field cells from the trochophore stage of Dreissena rostriformis. Gene ontology terms of these genes showed enrichment in shell formation-associated processes such as vesicle-mediated transport, phospholipid metabolic processing, integrin-mediated signaling pathway, and positive regulation of cell cycle G2/M phase progress (Salamanca-Díaz et al., 2022; Supplementary Table S2). Tertiary structure analysis using the Phyre2 webserver coincide with and thus confirm the results from SignalP and Interproscan (Supplementary Tables S2, S5, S6 and Supplementary File S1). In addition, expression dynamics of shell matrix genes during development of the trochophore of D. rostriformis were analyzed and compared with shell field-specific genes from pre-metamorphosis stages of the oyster Crassostrea gigas using previously published RNA-seq data (Figure 2, Supplementary Tables S7, S8) (Zhao et al., 2018; Calcino et al., 2019). Expression of most of these genes starts early in development, i.e., shortly after fertilization, likely by maternal transcripts. High normalized peaks of transcription are seen throughout the late gastrula and trochophore stages (during which the protoconch I is established) and continue in the veliger stages (continuous protoconch II formation) i.e., between 13 and 48 hpf. In situ hybridization experiments of some of these genes have previously shown a high level of expression in stages of embryonic shell (protoconch I) formation (e.g., Hox 1, hic31) (Salamanca-Díaz et al., 2021; Salamanca-Díaz et al., 2022). Furthermore, numerous orthogroups that contain genes that are specific to the embryonic shell field are shared across Metazoa (Supplementary Tables S2, S4. This demonstrates multiple cooption events of these genes into various functions in the respective metazoan lineages (Supplementary Table S1). However, we also found a group of genes (68 genes in 61 orthogroups; e.g., engrailed, cyclin-A2, carbonic anhydrase, tyrosinase homologs) active in D. rostriformis shell field formation that are also involved in shell formation of other mollusks (Nederbragt et al., 2002; Iijima et al., 2008; Kin et al., 2009; Samadi and Steiner, 2009; Huan et al., 2020; Zhao et al., 2020). These genes are also present in all other metazoans screened for herein and are commonly known to be related to body plan specification, cell cycle, and metalloenzyme activity (Figure 1 and Supplementary Table S2).

A closer analysis of specific taxonomic orthogroups (e.g., Protostomia, Lophotrochozoa) revealed that the majority of the genes (48.7%) that are differentially expressed in the shell field in D. rostriformis have orthologs in other taxa. However, their distribution between the given taxa is highly variable (Figure 1, Supplementary Tables S2, S4). While 19% of the shell field-specific genes are shared with other metazoan taxa, only 2.24% of the total number of shell field-specific genes are shared with other protostome species (8 genes in 6 orthogroups), and 1.1% of the same total number of genes are restricted to lophotrochozoans (4 genes in 4 orthogroups). Within molluscs, 1.1% of the shell field-specific genes are shared with other conchiferans, another 1.1% were only found in the sampled bivalves and gastropods, and 2.5% are possibly bivalve-specific. A quarter of the shell field-specific genes were only found in D. rostriformis and are not shared with other taxa. These may either genus- or species-specific genes, however their evolutionary history needs to be assessed in more depth once more bivalve datasets become available. The majority of genes in all orthogroups do not have a match in the InterPro database but show low level similarities with human orthologs (e-values higher than 1). Among the few genes with annotations in these groups, there is a member of the Claudin protein family, a keratin ortholog, epidermal growth factor domains, heat shock 70 kDa protein, and an endonuclease 2 ortholog. Genes expressed in the shell field which are restricted to Mollusca lack confident annotations (Supplementary Table S2). Human blast hits show only few domain commonalities to genes with a role in protein modification, DNA repair, nuclear envelope component and ion exchange, i.e., DDB1 and CUL4 associated factor 4, DNA repair protein XRCC1, nuclear envelope integral membrane protein 2, and sodium-driven bicarbonate exchanger.

The number of hitherto non-annotated genes shows a tendency to decrease when analyzing the different lineages inside Mollusca. This suddenly changes in the branch leading to Dreissena rostriformis, where the number of shell field-specific genes notably increases (Figure 1). Within Conchifera, a putative Bivalvia + Gastropoda clade shows 2 hitherto undescribed genes in 2 separate orthogroups and Bivalvia alone presents a unique set of 9 genes in 9 orthogroups (Figure 1, Supplementary Table S2). In this gene set, there are low e-values and little similarity to human as well as Crassostrea gigas orthologs, with blast hits to genes associated to antioxidant reactions, cell migration, cell attachment, and cellular proliferation, i.e., superoxide dismutase, myomegalin, laminin subunit beta-4, and ETS domain-containing transcription factor ERF. Additionally, from the genes expressed in the D. rostriformis embryonic shell field which were not assigned to any orthogroup or have a specific identity, and thus are considered here for Dreissena to be lineage-specific (86 genes in total, 24% of all shell field genes), 39 have transmembrane domains and 41 have signal peptides. This suggests that almost half of this gene subset is probably crucial for cell signaling since it has domains that interact directly with the outside of the cell membrane (Supplementary Tables S2, S5, S6). Altogether, our results show that, while there is a core gene set expressed in the embryonic shell field which is present throughout Metazoa, there is also strong indication of novel gene emergence that is specific to the embryonic shell field of the Dreissena trochophore.

Previous studies have characterized the shell secretomes from larval and adult stages of two marine bivalves, Pinctada fucata and Crassostrea gigas. They found that, despite having some common gene expression signatures (e.g., carbonic anhydrase, chitin binding protein, and von Willebrand factor type A), they also show distinct expression patterns of larval shell matrix proteins depending on species and developmental stages. One significant subset of the genes (around 90 out of 156 genes) involved in shell secretion is expressed in trochophore stages, while the other genes are expressed during the later D-shape veliger stages, suggesting different molecular signatures underlying embryonic versus larval shell formation (Zhao et al., 2018). This calls into question the homology of embryonic and larval shells in Bivalvia. Since solid data on the genes involved in bivalve teleoconch formation are still lacking, evolutionary relationships between the adult and the two transitory protoconch shell types currently remain unknown. This underlines that more in-depth comparative studies are needed to assess the decades-old question of (ontogenetic) homology of conchiferan embryonic, larval, and adult shells within the respective sublineages (particularly bivalves, gastropods, and scaphopods).

Our study shows that 24% (86) of the genes differentially expressed in the shell field of the trochophore of D. rostriformis could not be assigned to any orthogroup and may thus be genus- or species-specific (Figure 1A). From these, 13 unassigned genes have only incomplete annotations in specific regions of each gene sequence, 41 have low e-value similarity with human or Crassostrea orthologs, and 32 of these genes have no known annotation or ortholog match (Supplementary Table S2). Similar trends are also known from other mollusks, where unassigned and undescribed genes expressed in shell- and plate-forming cells appear to be highly taxon-specific. For example, secretomes from adult gastropods, bivalves, polyplacophorans, and a nautiloid cephalopod show considerable levels of lineage-specific orphan genes (Jackson et al., 2006, 2009; Immel et al., 2016; Kocot et al., 2016; Marin, 2020; Setiamarga et al., 2021). It has been argued previously that the rapid evolutionary rate of these genes may be a possible reason for the lack of orthology detection of these shell matrix toolbox genes (Aguilera et al., 2017). This, in turn, could be the result of evolutionary responses to the widely varying ecological conditions shell-bearing mollusks are exposed to, since most of the gene products in the shell field are in direct contact with the environment. Interestingly, almost half of these lineage-specific orphan genes have transmembrane domains and/or signaling peptides (Supplementary Tables S2, S5, S6). This suggests that genes expressed in the shell field at the trochophore stage might be significantly influenced by the ecology of the larva. Previous studies have found that molluscan shell proteomes drastically change when ecological factors such as the pH or the temperature are altered, but combined experimental and transcriptomic studies are currently too scarce for robust conclusions on an evolutionary level (Timmins-Schiffman et al., 2014; Wei et al., 2015). However, since the environmental conditions during protoconch I and protoconch II formation are identical in D. rostriformis, this might hint towards an independent evolutionary origin (and thus argue against ontogenetic homology) of these shell types. This is further supported by the fact that, after shell field formation, there is a fluctuation of gene expression throughout development, i.e., Chitin binding domain ortholog (Gene.49769) and voltage-dependent calcium channel subunit alpha-2/delta-4 human ortholog (Gene.25093) (Figure 2A), demonstrating putatively different expression dynamics during protoconch I and protoconch II formation, respectively. A similar tendency emerges when comparing temporal expression dynamics of shell-specific genes of D. rostriformis with C. gigas (Figure 2B). The Pacific oyster seems to have different sets of genes with alternate expression throughout developmental stages where shell field formation is active, just as in Dreissena, thus calling into question the homology of bivalve ontogenetic shell types (cf. Zhao et al., 2018). However, further comparative studies employing different developmental stages of the same as well as similar developmental stages of different species are needed to further assess this assumption.

Our single-cell RNAseq and OrthoFinder analyses grouped 271 (76%) out of 357 identified genes that are differentially expressed in shell field cells from the trochophore stage of Dreissena rostriformis into orthogroups shared with different taxa. From these resulting orthogroups, there is a fraction of Dreissena trochophore shell field-specific genes that are shared with the rest of the sampled metazoans (68 genes in 61 orthogroups, 19%) (Figure 1, Supplementary Table S2). It has been shown previously that mantle secretomes in other bivalves and gastropods also possess a wide range of gene families that originated prior to the emergence of the conchiferan clade (Kocot et al., 2016; Aguilera et al., 2017). Among this, a set of genes from the shell field of Dreissena, which are present in other metazoans, is known to be involved in extracellular matrix formation, such as orthologs of laminin, C-type lectin domains, and immunoglobulins (Supplementary Table S2). Additionally, among this set there are genes containing leucine-rich repeat domains and semaphorins, which are also found throughout metazoans (Supplementary Table S2). Moreover, genes from this subset of orthogroups were coopted into embryonic shell formation in conchiferan mollusks such as Dreissena, e.g., Hox 1 (Gene.152834), Hox 4 (Gene.66474), Lox 4 (Gene.142102), and engrailed (Gene.126286) (Figure 1, Supplementary Table S2) as previously described for other species of mollusks (Samadi & Steiner, 2009; Fritsch et al., 2015; Wollesen et al., 2018; Huan et al., 2020; Liu et al., 2020; Salamanca-Díaz et al., 2021). Furthermore, in these orthogroups there are genes that have been found to be also involved in biomineralization processes in echinoderms and vertebrates, e.g., Cyclophilin-type (Gene.103270) and Carbonic anhydrase (Gene.82229) (Livingston et al., 2006; Mann et al., 2008; Mann et al., 2010; Mann and Edsinger, 2014). Such an organic matrix is formed prior to secretion of the mineralized part of the shell and is thus of crucial importance for conchiferan mollusks, but the respective factors involved are also present in other metazoans that lack a shell (Supplementary Table S2) (Marie et al., 2010; Marie et al., 2011a; Marie et al., 2011b; Marie et al., 2012; Marin et al., 2014). Altogether, our data point towards a shared “molecular biomineralization toolbox” across Metazoa, but a broader taxon sampling especially from key invertebrate phyla are required for deeper evolutionary insights. Given the fact that numerous animal phyla contain taxa with mineralized hard parts, including accessible representatives such as annelids, brachiopods, other lophotrochozoans, as well as numerous arthropods, this hypothesis can be tested by comparative studies using single-cell RNA transcriptomic approaches.

Taken together, the quagga mussel Dreissena rostriformis shows a mosaic of co-option of known metazoan genes and de novo recruitment of genes with hitherto unknown function or ortholog match into embryonic (protoconch I) shell formation. Our data suggest that not only adult but also embryonic bivalve shells are highly plastic in the gene repertoire that underlie their ontogeny, which may be indicative of non-homology of bivalve—and possibly conchiferan - ontogenetic shell types.