Sexually mature M. lateralis adults were selected from the culture system of MOE Key Laboratory of Marine Genetic and Breeding, and spawning was induced with drying in the shade followed by 26°C thermal shock method. Fertilized eggs of M. lateralis were incubated in a 15-L tank with filtered seawater at 22°C. After removal of abnormal embryos and larvae, samples of four developmental stages, blastula (3 h post fertilization (hpf)), gastrula (5 hpf), trochophore larva (8 hpf), and D-shaped larva (14 hpf) were collected ( Figure S1 ), immediately frozen in liquid nitrogen, and preserved at −80°C until RNA extraction.

The total RNA of 12 samples, with three parallel of each stage, was extracted individually by phenol–chloroform extraction (Liu et al., 2015). The extracted RNAs were analyzed with 1.2% agarose gel electrophoresis and NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA). RNA concentration was measured using a Qubit 2 fluorometer (Invitrogen, Carlsbad, CA, USA), and integrity was assessed using the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). Then, the quantified RNA samples were subsequently used for cDNA library construction.

Equal amounts of RNAs from these four developmental stages of M. lateralis were pooled together for the construction of PacBio isoform sequencing (Iso-seq) library. The library was prepared according to the Isoform Sequencing protocol (Iso-Seq) using the Clontech SMARTer PCR cDNA Synthesis Kit (Clontech, Mountain View, CA, USA) and the BluePippin Size Selection System (Sage Science, Beverly, MA, USA) as described by Pacific Biosciences (PN 100-092-800-03). PacBio library was sequenced on the PacBio Sequel platform (Pacific Biosciences, Menlo Park, CA, USA).

Sequence data were processed using the SMRTlink 5.1 software (Hon et al., 2020). The subreads were processed by the Circular Consensus Sequence (CCS) algorithm to obtain the CCS (min_length 50, max_drop_fraction 0.8, no_polish TRUE, min_zscore −9,999.0, min_passes 2, min_predicted_accuracy 0.8, max_length 15,000). CCSs were then classified into full-length and non-full-length reads using pbclassify.py (ignore polyA false and minSeqLength 200). All fasta files produced were then fed into the cluster step for isoform-level clustering (ICE), followed by final Arrow polishing (hq_quiver_min_accuracy 0.99, bin_by_primer false, bin_size_kb 1, qv_trim_5p 100, and qv_trim_3p 30) (Cao et al., 2020). To further improve the sequencing accuracy, additional nucleotide errors in consensus reads were corrected using the Illumina RNA-Seq data of the present study with the software LoRDEC (Salmela and Rivals, 2014). Any redundancy in consensus reads was removed by CD-HIT (-c 0.95, -T 6, -G 0, -aL 0.00, -aS 0.99) (Fu et al., 2012) to obtain the final transcripts, which were then considered as the reference transcriptome (designated as ref) for the subsequent analysis.

The obtained non-redundant transcript sequences were functionally annotated using the following databases: National Center for Biotechnology Information (NCBI) non-redundant protein (NR), EuKaryotic Orthologous Groups (KOG), Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG), NCBI nucleotide sequences (NT), Protein Family (Pfam), and Gene Ontology (GO). The first four database annotations were performed using diamond v0.8.36 with an e-value threshold of 1e−5 (Buchfink et al., 2015). The ncbi-blast-2.7.1+ with an e-value threshold of 1e−5 was utilized for NT database annotation, and the Pfam databases annotations were performed using Hmmscan of the HMMER 3.1 package (Nguyen et al., 2016). Script of protein annotation results based on Pfam database was used for GO annotations (Ashburner et al., 2000).

Predictive analysis of coding sequence (CDS) was performed using ANGEL v. 2.4 software, and bivalve species confident protein sequences were used for ANGEL training (Shimizu et al., 2006). Transcription factors were predicted using the AnimalTFDB 2.0 database with iTAK software (Tatusov et al., 2003). Coding–Non-coding index (CNCI) (Sun et al., 2013), Coding Potential Calculator (CPC) (Kong et al., 2007), Pfam-scan (Finn et al., 2016), Predictor of Long Non-coding RNAs, and Messenger RNAs based on an Improved k-mer scheme (PLEK) (Li et al., 2014) were integrated to identify long non-coding RNAs (lncRNAs) in the transcripts. Transcripts longer than 200 nt without coding potential were selected as lncRNA candidates.

Simple sequence repeats (SSRs) of the transcriptome sequences were identified using MISA 1.0 with default parameters (Beier et al., 2017). The minimum repeat time for core repeat motifs was set as 10 for mononucleotides, 6 for dinucleotides, and 5 for trinucleotides, tetranucleotides, pentanucleotides, and hexanucleotides.

Twelve RNA samples of M. lateralis embryos/larvae from four developmental stages were used for Illumina cDNA library construction using NEBNext Ultra™ RNA Library Prep Kit for Illumina (NEB, Ipswich, MA, USA) following the manufacturer’s instructions. Qualified libraries were applied to next-generation sequencing (NGS) using an Illumina Hiseq X Ten platform (Illumina, USA) to generate 150-bp paired-end sequence reads. Raw reads were firstly preprocessed for quality control to obtain the clean reads by removing reads containing adapters, ambiguous N nucleotides, and low-quality reads. Then, all clean reads obtained were mapped onto the PacBio Iso-Seq ref with bowtie2 (-end-to-end, -sensitive) (Langmead and Salzberg, 2012).

The mapping results from bowtie2 were calculated by RSEM (Li and Dewey, 2011), and the read count values of each gene from each sample were converted into FPKM (expected number of fragments per kilobase of transcription sequence per mills base pairs sequence) to determine the expression level of each gene. Both heatmap and hierarchical cluster analyses were carried out using the R package “pheatmap”.

DESeq 1.10.1 package (Anders and Huber, 2010) was used for differential expression analysis of genes between two sequential developmental stages. Unigenes with false discovery rate (FDR)< 0.01 and |log₂(FoldChange)| > 1 were considered to be the differentially expressed genes (DEGs).

Gene Ontology enrichment analysis of DEGs between different stages was implemented by the R package “GOseq” (Young et al., 2010), in which gene length bias was corrected. GO terms with a Benjamini and Hochberg (BH)-corrected p-value of less than 0.05 were considered significantly enriched. The enrichment analysis of the KEGG pathway was performed using KOBAS (2.0) (Xie et al., 2011). The p-value correction method and the threshold value were consistent with the above.

The R package “WGCNA” was utilized to construct the co-expression network (Langfelder and Horvath, 2008), and 25,664 annotated genes with FPKM values greater than 0.1 were selected. Firstly, samples were clustered to assess the presence of any obvious outliers. Then, the R function pickSoftThreshold was used to calculate the soft thresholding power β, and the connectivity between genes met a scale-free network distribution (scale-free R ² = 0.9) when the value of soft thresholding power β was set to 6. Hierarchical clustering and the dynamic tree cut function were used to detect modules, and gene significance (GS) and module membership (MM) were calculated to relate modules with clinical traits. The minimum module size was set to 30 (CutHeight = 0.99, minModuleSize = 30), and the modules with a higher correlation were merged (r< 0.25). The corresponding gene information of each module was extracted for further analysis. Lastly, Cytoscape 3.7.0 software was used to visualize the network of eigengenes.

Six unigenes related to larval shell formation (tyr, chitin synthase, CA, EF-hand, a calmodulin coding gene (designated as C-2442), and solute carrier family 4 member 8 (slc4a8)) were selected to confirm the RNA-Seq results with quantitative RT-PCR (RT-qPCR) according to their expression difference and functional annotation. RNA samples were isolated from M. lateralis of different developmental stages, and first-strand cDNAs were synthesized using PrimeScript™ RT reagent Kit (Takara, Dalian, China). Primers were designed with Primer Premier5 ( Table S1 ), and Ef-1b was used as a reference gene (Li et al., 2019). PCRs were conducted using SYBR Green Real-Time PCR Master Mix (Takara, China) and the LightCycler 480 Real-time fluorescence quantitative PCR instrument (Roche, Basel, Switzerland). Data were analyzed using the 2^−ΔΔCt method (Livak and Schmittgen, 2001) and presented as mean ± SD. All reactions were performed with three biological replicates and three technical replicates. Statistical analysis was tested using one-way ANOVA followed by Duncan’s test (SPSS 22.0, USA), and a statistically significant difference was set as p< 0.05.

Through PacBio Iso-Seq sequencing, a total of 6,623,650 subreads (15.11-Gb nucleotides) were obtained after initial quality control by removal of the adaptor reads and subreads<50 bp ( Table 1 ). The subreads were then self-corrected, and a total of 337,922 CCSs were produced, of which 238,919 sequences were identified as full-length non-chimeric (FLNC) reads. After iterative clustering for error correction (ICE) and polishing using Arrow, 184,784 high-quality full-length consensus transcripts were generated. After the raw reads generated were trimmed and filtered by the Illumina RNA-Seq, more than 54,095,792 high-quality clean reads ( Table 2 ) were generated for each sample with Q30 values greater than 92.18%. To build the full-length transcripts with higher accuracy and lower error rate, all full-length CCSs were corrected with Illumina RNA-Seq data by LoRDEC. At last, 121,424 unigenes, as the gene reference sequences, were obtained after the removal of redundant sequences by CD-HIT ( Table 1 ). Then, the clean reads of each sample obtained by Illumina sequencing were mapped to the ref, and the mapping rate ranged from 64.64% to 84.75% ( Table 2 ).

This study combined the PacBio Iso-Seq and Illumina RNA-Seq technologies to generate a full-length transcriptome for M. lateralis. After correction, the average length of full-length consensus transcripts was 3,267 bp with an N50 of 3,879 bp ( Table S2 ), which were much longer than that of the de novo assembled transcripts (Hu et al., 2019; Nie et al., 2021; Núñez-Acuña et al., 2022) and comparative with the full-length transcriptome of bivalve species recently reported (Liao et al., 2022; Zeng et al., 2022). The average mapping rate of Illumina sequences was greater than 75%, indicating that a good reference dataset was established. The results suggested the full-length transcriptome obtained in the present provided good genetic resources for M. lateralis.

The identified 121,424 unigenes were annotated by searching against seven databases, and a total of 73,367 (60.42%) ( Table S3 ) were annotated in at least one database, as well as 10,036 were annotated in all seven databases ( Table 3 , Figure 1A ). According to the prediction by NR database, 62,786 unigenes were annotated in 642 homologous species, where the top 20 species were shown in Figure 1C . C. gigas, whose genome was the first to be described as a marine bivalve (Qi et al., 2021), was found to contain the highest number of homologous sequences (29,485, 46.96%) with M. lateralis. The second and third species were Lottia gigantea (5628, 4.63%) and Lingula anatina (3576, 2.95%), respectively. More than 60% of unigenes have been annotated in at least one database, which provided important information for studying the shell formation as well as other development processes of M. lateralis.

There were 37,777 (31.11%) unigenes assigned to the three main GO categories and then categorized into 55 GO terms. The three most abundant terms in the biological process were “cellular process”, “metabolic process”, and “single-organism process”. Within the cellular component category, “cell part”, “cell”, and “organelle” were the prominent terms. “Binding”, “catalytic activity”, and “transporter activity” were dominant in the molecular function category ( Figure 1D ). The term “extracellular matrix” associated with the cellular component category is generally supposed to be related to shell formation genes (Song et al., 2022), and collagen and calcium-binding protein were recognized in this term.

A total of 59,418 (48.93%) transcripts were annotated in 1,280 pathways in the KEGG database. According to the KEGG classification, the unigenes were assigned to six level 1 KEGG categories (cellular processes, environmental information processing, genetic information processing, metabolism, organismal systems, and human diseases) and 45 level 2 subcategories ( Figure 1B ). The top five prominent subcategories were “Signal transduction”, “Cancers: Overview”, “Transport and catabolism”, “Endocrine system”, and “Cell growth and death”. Notably, it was reported that the “calcium signaling pathway”, “Wnt signaling pathway”, and “tyrosine metabolism pathway” probably play significant roles in the biomineralization process (Shi M. et al., 2013; Liu et al., 2015). Among these three pathways, calmodulin, tyr, β-catenin, TCF, and wnt5 are supposed to be related to calcium signal transduction and the shell formation process (Shi Y. et al., 2013; Gao et al., 2016; Zhang et al., 2021).

The KOG analysis demonstrated that 40,132 (33.05%) unigenes were classified into 26 functional clusters ( Figure 1E ). In M. lateralis, the most numerous five categories were “General function prediction only (R)”, “Signal transduction mechanisms (T)”, “Posttranslational modification, protein turnover, chaperones (O)”, “transcription (K)”, and “Intracellular trafficking, secretion, and vesicular transport (U)”. A total of 1,071 unigenes were assigned to “Inorganic ion transport and metabolism (P)”, 39 genes of which were annotated as members of the SLC4 family, a HCO 3 − transporter family that can supply calcified substrates ( Ca 2   + and HCO 3 − ) to calcified sites through cross-epithelial transport (Ramesh et al., 2017). Moreover, 6,636 unigenes were assigned to the category “Signal transduction mechanisms (T)”, within which tyr, CA, calmodulin, and chitin synthase were reported to be involved in shell formation (Miyamoto et al., 1996; Miglioli et al., 2019; Wang et al., 2022).

In this study, animalTFDB 2.0 was used to identify animal transcription factors in M. lateralis. A total of 3,015 transcripts were predicted as animal transcription factors ( Table S4 ). The first three transcription factor families were the zf-C2HC transcription factor family (1177), Homeobox transcription factor family (249), and ZBTB transcription factor family (238) ( Figure 2A ).

SSR is widely and evenly distributed in the genome of eukaryotes. The software MISA was used to search the SSR profiles in the unigenes of M. lateralis. A total of 39,049 SSR-containing unigenes were detected. In the four developmental stages, mono-nucleotide repeat motifs (9–12) were the most abundant, followed by di-nucleotide repeats and tri-nucleotide repeats (5–8) ( Figure 2B ).

Four computational approaches (CPC, CNCI, CPAT, and Pfam) were combined to identify lncRNAs. As is shown in Figure 2C , 36,386 unigenes were predicted as lncRNAs by all four methods, accounting for 29.97% of the total predicted lncRNAs.

To investigate the expression patterns of unigenes in M. lateralis, the Illumina clean reads were mapped back to the ref to determine their expression levels. The FPKM values of more than 95% of genes were ≤5, and approximately 1% of the genes were >15 ( Table S5 ). A Venn diagram was drawn to exhibit the unigenes with FPKM > 0.1 in all four stages (30,178 unigenes) ( Figure 3A ), and the number of genes uniquely expressed in each stage was 776 (blastula), 1,974 (gastrula), 4,002 (trochophore larva), and 3,072 (D-shaped larva). There were 327 (blastula), 366 (gastrula), 413 (trochophore larva), and 485 (D-shaped larva) unigenes with high expression (FPKM > 60) in each of the four developmental stages. Specific to these highly expressed unigenes (FPKM > 60,619 genes), a hierarchical clustering analysis was performed with the normalized value of log₁₀(FPKM + 1) ( Figure 3B ). The expression patterns of genes between the blastula and gastrula were similar, as most of the genes highly expressed in the blastula were also highly expressed in the gastrula. However, the genes highly expressed in the D-shaped larva were very low or even not expressed in the other three development stages ( Figure 3B ).

Pairwise comparisons between the two adjacent developmental stages (designated as b&g for the blastula and gastrula, g&t for the gastrula and trochophore larva, and t&d for the trochophore larva and D-shaped larva) were performed to identify the DEGs ( Table S6 ). The maximal DEGs were observed between the trochophore larva and D-shaped larva (17,829), while the fewest existed between the blastula and gastrula, (4512). Moreover, 10,637 DEGs were identified between gastrula and trochophore larva ( Table 4 ). Among those genes, tyr, CA, insoluble shell matrix protein-encoding gene 1 (ISMP1), and EF-hand, which play key roles in embryo shell development in bivalve molluscs (Zhang and He, 2011; Zhao et al., 2018), were identified to have significant differential expression between the trochophore larva and the D-shaped larva in M. lateralis.

GO and KEGG enrichment analyses were used to mine DEGs related to shell formation. The most representative level 2 GO terms of DEGs from the three groups were “binding”, “cellular process”, and “metabolic process” ( Table S7 ). In the “biological process” category, the top three GO terms for b&g were “organic cyclic compound metabolic process”, “cellular aromatic compound metabolic process”, and “heterocycle metabolic process”, which indicated that compound metabolic process may be very different between the blastula and gastrula. The DEGs in g&t were mainly assigned to “single-organism process”, “biological regulation”, and “regulation of biological process”. The most enriched GO terms for t&d were “single-organism process”, “single-organism cellular process”, and “localization”. In the “molecular function” category, more DEGs in g&t were enriched in “transition metal ion binding”, while DEGs of t&d were significantly enriched in “chitin binding” categories, as well as in “transcription factor activity”. In addition, 80, 172, and 254 DEGs in the b&g, g&t, and t&d groups, respectively, were involved in the “calcium ion binding function” GO category, some of which were involved in shell formation (EF-hand and calmodulin). The upregulated DEGs in the D-shaped larva such as collagen and calcium-binding proteins were involved in the “extracellular matrix”.

KEGG enrichment analysis showed that DEGs in b&g were mostly enriched in the pathways of “Huntington’s disease”, “Pathogenic Escherichia coli infection”, and “Phagosome” ( Table S8 , Figure 4A ). “Metabolic pathways”, “Biosynthesis of secondary metabolites”, and “Protein processing in endoplasmic reticulum” were the top KEGG pathways in the g&t group ( Table S8 ; Figure 4B ). For the t&d group, the top three pathways were “Metabolic pathways”, “Cell cycle”, and “Biosynthesis of secondary metabolites” ( Table S8 ; Figure 4C ). Moreover, more DEGs in the t&d group were enriched into the “calcium signaling pathway” and “endocrine and other factor-regulated calcium reabsorption” pathways than those in the g&t group. Thirty DEGs enriched in the “calcium signaling pathway” were annotated as calcium/calcium-like protein-encoding genes, which are involved in the regulation of calcium transport and secretion (Li et al., 2004). A total of nine DEGs enriched in “endocrine and other factor-regulated calcium reabsorption” were annotated as solute carrier family 8 (slc8a), which encode proteins mediating Ca²⁺ extrusion in most cell types (Brini and Carafoli, 2011). The increasing genes involved in the calcium reabsorption pathway during larval development suggested that the trochophore larva was the key stage for the individuals to absorb calcium ions, which was coincident with the discoveries in mussels (Ramesh et al., 2017).

To construct and analyze the gene co-expression network, 25,664 annotated genes with FPKM values greater than 0.1 were utilized. With the soft thresholding power β of 6 and the minimum module size of 30, a total of 14 modules were identified through hierarchical clustering ( Figure 5A ). To explain the gene expression variation, module eigengenes (MEs) were calculated to represent each module. Relationships between the identified modules were mapped ( Figure 5B ), and the results indicated that gene expressions were relatively independent between modules ( Figure 5D ). With the connectivity analysis of eigengenes, these 14 modules could be clustered into two clusters ( Figure 5C ).

The developmental stage (blastula, gastrula, trochophore larva, and D-shaped larva) was used as the traits to select the trait-related modules, and the association between the modules and traits was explored by measuring the correlation between ME values and features. The results indicated that pink and salmon modules were most positively correlated with D-shaped larva (r = 0.99, p = 4e–09; r = 0.89, p = 5e−04), and negative correlations were observed with the other three traits ( Figure 6A ). Meanwhile, a high degree of connectivity between pink and salmon modules was revealed ( Figure 5C ). In the pink module, insoluble matrix shell protein, tyr, and chitin-binding protein genes could be identified. Moreover, EF-hand and C-2442 were found in the salmon module. GO analysis of genes in the modules pink and salmon was conducted, and the results showed that the genes within the pink module were significantly revealed to have multiple functions associated with biomineralization, such as calcium ion binding and transmembrane transporter activity ( Figure 6B ), while the genes that belonged to the salmon module were significantly enriched in translation and protein metabolic process ( Figure 6C ).

To further analyze the functions represented by the pink module, which exhibited the highest concentration in the D-shaped larva, we use a |GS| over 0.2 and an |MM| over 0.8 as cutoff criteria and identified 948 hub genes. According to the recommendation of the weighted gene correlation network analysis (WGCNA) official website, the top 85 genes with the highest weight value were selected, and the co-expression network was visualized using Cytoscape ( Figure 6D ). Among these hub genes, slc4a8, a bicarbonate transporter from the slc4 family, was identified, which is involved in calcium transport and shell repair, and calcification during blue muscle larval development (Ramesh et al., 2019). We found that rnf41 (E3 ubiquitin-protein ligase NRDP1) was highly connected with slc4a8 (weight = 0.6144), and these two genes may regulate the calcification process during the early shell formation (Fang et al., 2012).

To focus on the analysis of larval shell formation in M. lateralis, all DEGs were retrieved, and 12 candidate genes were obtained to be related to biomineralization. Among these 12 DEGs, 8 unigenes function in shell matrix deposition (ISMP1, ISMP2, ISMP5, chitin synthase, tyr, chitin-binding protein, collagen, and pu14), and 4 unigenes were related to ion transportation (CA, slc4a8, C-2442, and EF-hand).

In order to verify the genes that were putatively involved in the biomineralization process, we selected six unigenes, including two genes involved in matrix deposition (chitin synthase and tyr), two HCO 3 − transporter genes (CA and slc4a8), and two calcium transporter genes (EF-hand and C-2442) for RT-qPCR analysis. As expected, the expression profiles of these genes showed similar trends with the RNA-Seq data ( Figure 7 ).

CA and chitin synthase exhibited similar expression patterns that considerably increased to the peak in trochophore larva and then dropped quickly in D-shaped larva. It is worth noting that the expression of CA dramatically increased to 503.8-fold that in the blastula. CA plays several essential roles in multiple physiological processes such as pH regulation, electrolyte balance, ionic transportation, carboxylation or decarboxylation reactions, and biocalcification, which is considered one of the most important enzymes in CaCO₃ biomineralization due to its ability to catalyze the hydration of carbon dioxide and provides bicarbonate ions that react with Ca²⁺ to form CaCO₃ (Supuran, 2011; Ramesh et al., 2017). It is proposed that molluscs usually absorb calcium first and CA catalyzes the hydration of CO₂ reacting with calcium to facilitate its deposition at the stratum corneum formed by the shell matrix after calcium accumulates to a certain extent in the body (Evans, 2019; Sharker et al., 2021). Chitin synthases are crucial enzymes for the metabolism of chitin. Arthropod chitin forms a major part of the exoskeleton, and the supporting layer of the tubular trachea system and chitin synthase is a key factor to regulate most insect growth and molting (Jiang et al., 2021). In shellfish, chitin synthase has been proven to be involved in shell biosynthesis of trochophore larva chitin shell, which is then covered by the calcified shell at early D-shaped larva (Zhang et al., 2012; Liu et al., 2018). We observed that during the development of M. lateralis, the shell formation process was completed within 1 h. The obvious expression fluctuation suggested that CA and chitin synthase played a key role during the larval shell formation and calcium accumulation process in such a short time window.

slc4a8 and C-2442 began to upregulate in the trochophore larva and reached the peak in D-shaped larva. SLC4 family bicarbonate transporters mediate the co-transport of sodium and bicarbonate. In sea urchin embryos, the SLC4 family member regulates intra- and extracellular exchange of HCO 3 − with Na⁺ and is fundamentally linked to the biological precipitation of CaCO₃ (Hu et al., 2018). In molluscs, some SLC4A family members have been identified in mantle transcriptome data of C. gigas (Zhang et al., 2012), Pinctada fucata (Shi Y. et al., 2013), and Tridacna squamosa (Ip et al., 2017), but the function of these transporters and a description of the Ca²⁺ transporting mechanisms still need to be explored. C-2442 codes calmodulin, and the protein has been found to mediate the enhanced calcium deposition induced by CO₂ exposure (Wang et al., 2022). The expression characteristics of both genes suggested that ion transportation mainly began at the trochophore larva and continuously provided a substrate for biomineralization.

EF-hand also encodes calmodulin, which is structurally conserved and functionally preserved throughout the animal kingdoms (Xu et al., 2017). It has been proven to be served as the calcium ion transporter during shell formation in a variety of shellfish (Huang et al., 2007; Sun et al., 2015; Jain et al., 2017). In molluscs, tyr is strongly expressed in the edge of the mantle, which is the main biomineralization tissue, and tyrosinase proteins have also been detected in the proteome of the shell prismatic layer and nacreous layer, indicating its key role in the shell matrix deposition (Huan et al., 2013; Miglioli et al., 2019; Ren et al., 2020; Zhu et al., 2021). Moreover, tyrosinase was confirmed to participate in the embryo shell ontogeny of the organic matrix in M. galloprovincialis and Hyriopsis cumingii (Miglioli et al., 2019; Ren et al., 2020). In M. lateralis, the expression pattern of these two genes was obviously upregulated in D-shaped larva. The results indicated that these genes may be involved in the formation of larval calcified shells. The characteristics of these genes suggested that the process of ion transport occurs earlier than the deposition of the shell matrix, and the trochophore larva and D-shaped larva were the key stages of larval shell formation. Further studies are needed to functionally validate the regulatory mechanisms of larval shell formation, and M. lateralis will provide a promising platform to facilitate the related research.