Skip to main content
Download PDF

Comparative transcriptomic analysis revealed dynamic changes of distinct classes of genes during development of the Manila clam (Ruditapes philippinarum)

BMC Genomics volume 23, Article number: 676 (2022) Cite this article

Abstract

Background

The Manila clam Ruditapesphilippinarum is one of the most economically important marine shellfish. However, the molecular mechanisms of early development in Manila clams are largely unknown. In this study, we collected samples from 13 stages of early development in Manila clam and compared the mRNA expression pattern between samples by RNA-seq techniques.

Results

We applied RNA-seq technology to 13 embryonic and larval stages of the Manila clam to identify critical genes and pathways involved in their development and biological characteristics. Important genes associated with different morphologies during the early fertilized egg, cell division, cell differentiation, hatching, and metamorphosis stages were identified. We detected the highest number of differentially expressed genes in the comparison of the pediveliger and single pipe juvenile stages, which is a time when biological characteristics greatly change during metamorphosis. Gene Ontology (GO) enrichment analysis showed that expression levels of microtubule protein-related molecules and Rho genes were upregulated and that GO terms such as ribosome, translation, and organelle were enriched in the early development stages of the Manila clam. Kyoto Encyclopedia of Genes and Genomes pathway analysis showed that the foxo, wnt, and transforming growth factor-beta pathways were significantly enriched during early development. These results provide insights into the molecular mechanisms at work during different periods of early development of Manila clams.

Conclusion

These transcriptomic data provide clues to the molecular mechanisms underlying the development of Manila clam larvae. These results will help to improve Manila clam reproduction and development.

Peer Review reports

Background

Mollusks, which can be aquatic or terrestrial, are the most morphologically diverse group of invertebrate species. However, this phylum exhibits the conserved and characteristic developmental processes of spiral cleavage and trochophore larvae. The early development of mollusks consists of a series of developmental periods during which changes in body form and lifestyle occur [1]. Bivalves are the second largest group in the phylum of Mollusca with regard to species number [2]. Bivalves include almost all economically important shellfish (i.e., oysters, scallops, clams, and mussels). The morphology of the bivalves changesdramatically during early development, and their lifestyle changes from planktonic to benthic or sessile [3,4,5]. Therefore, studies of the molecular processes during early development of bivalves will provide the basis for understanding the evolution and early diversification of animals [6].

The body morphology and lifestyle of bivalves change during development from larvae to adults, and many genes are involved in the transformation from free-swimming to benthic larvae and in the metamorphosis stage [7]. The expression of caveolin, tubulin and tektin provides a reference for the early development of ciliary bands, which are the first organs to emerge during molluscan embryonic differentiation [8]. Tyrosinase may be involved in the biosynthesis of non-calcified shells in Pacific oyster Crassostrea gigas [9]. Dopamine receptor genes may be involved in the larvae shell formation in Pacific oyster C. gigas [10]. In the Meretrix meretrix, mRNA differential expression analysis showed that expression of the MmeFer gene in larvae increased at least 8-fold after the trochophore stage, indicating that it played an important role in the calcification process of the shell in this species [11].

Larval metamorphosis in most bivalves is accompanied by high mortality due to both external and internal factors [12]. A few receptor genes, functional genes, and chemical cues control the metamorphosis of planktonic larvae to benthic larvae [13]. For example, steroid receptors identified by combined multiomics analysis may play an important role during larval metamorphosis in pearl oyster(Pinctada fucata martensii) [14]. In the Pacific oyster, myosin heavy chain knockdown was found to regulate proper myoblast assembly during larval development [15]. Caspase gene expression and in situ hybridization results suggested an important role of this gene in the process of veliger shedding during larval metamorphosis of the Fujian oyster (Crassostreaangulata) [16], and it may also be involved in larval metamorphosis of Meretrix via the apoptotic pathway [17]. The specific expression pattern of molluscan growth and differentiation factor in C. gigas suggested that this factor may play a central role in the metamorphosis process of oyster larvae [18]. The N-methyl-D-aspartate pathway has a recognized regulatory function in the metamorphosis of the Pacific oyster [19], and many neuroactive compounds, including γ-aminobutyric acid, various catecholamines, 5-hydroxytryptamine, and nitric oxide (NO) are associated with metamorphosis in various mollusks [20].

The Manila clam is widely distributed in the sandy mud sediments of tidal flats and shoals along the coasts of China, the northwestern United States, and the Atlantic coast of Europe. It is an economically important marine bivalve, with a production of 4.2 million tons in 2020 [21]. Due to its economic importance and ecological role, many researchers have studied the genetics, breeding, and disease control of this species. Fertilization and hatching rates are important indicators in artificial breeding programs because it has a significant impact on shellfish larval production [22]. Studies have shown that the production of clam juvenile is inextricably linked to hatching rate and larval mortality [23], and it has been demonstrated that functional genes play important roles on the hatching and development process [24]. Understanding the regulatory steps at work during the different stages of early development is critical to artificial breeding and developmental biology research. The growth traits variability and developmental process of the Manila clam has been described in our previous study [25, 26], however, the dynamic changes of distinct classes of genes during embryo and larvae development of the Manila clam is still largely unexplored. In this study, we focused on the molecular processes involved in different developmental stages of Manila clam. Our comprehensive analysis of gene regulation during early development provided a framework for understanding the developmental dynamics of the Manila clam.

The goals of this study were to reveal the molecular mechanisms involved in the early development of the Manila clam and to elucidate the relationship between gene expression and morphological alterations. To investigate the transcriptional processes that occur during clam development, we used RNA-seq to explore the genetic basis of and regulatory genes involved in early development of the Manila clam and to uncover important biological processes, molecular functions and developmental regulation of gene expression during ontogeny.

Materials and methods

Manila clam was collected from the Dalian, Liaoning Province, China, and the gonadal maturity of the population was assessed by regular sampling. After the clams spawned, we mixed sperm with eggs while stirring clockwise to allow full fertilization until fertilized eggs (FE) could be observed (Fig. S1). Samples were also collected and stored for each of the remaining stages. After 0–15 min, the number of first polar body was > 90% (PB1), and specimens were collected and stored in liquid nitrogen. After 20–44 min, the number of second polar body was > 90% (PB2), and specimens were sampled and stored in liquid nitrogen. After 1 h, the eggs entered the two-cell stage (TC), after 84–107 min it entered the eight-cell stage (EC), after 160–180 min it entered the blastula stage (B), and from 3 to 6 h it entered the gastrula stage (G),and then hatched to become trochophore larvae (T), followed by D-pattern larvae (D) after 1 day and umbo-veliger larvae (U) after 2–5 days. After 15 days, the larvae entered the pediveliger larval stage (P). The larvae metamorphosed between 19 and 21 days and grew from a single pipe juvenile (S) to a double plumbed juvenile (J) (Table S1). All methods were carried out in accordance with relevant guidelines and regulations.

Transcriptome sequencing and assembly

Total RNA samples were extracted, and their quality was tested using an Agilent 2100 bioanalyzer. The mRNA was then enriched, and the library was built using the NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina®, operating on 42 sets of samples. After library construction, the libraries (diluted to 1.5 ng/µl) were initially quantified using a Qubit 2.0 Fluorometer and subsequently tested for insert size using the Agilent 2100 bioanalyzer. For insert sizes that met expectations, qRT-PCR was performed to quantify the effective concentration of the library. The libraries were then accurately quantified by qRT-PCR (effective library concentration above 2 nM) to ensure library quality.

Quality control and reads mapping to the reference genome

Analysis of raw image data files generated by the Illumina HiSeq 2500 instrument, where sequenced fragments were converted into sequence data (reads) by the high-throughput sequencer with CASAVA base identification. The raw reads were processed using an internal Perl script [27]. In this step reads with adapters were removed; reads containing N were removed (N indicates that base information could not be determined), and low-quality reads (reads with Qphred < = 20 bases accounted for more than 50% of the entire read length) were removed. After raw data filtering, sequencing error rate checking and GC content distribution checking, clean reads were obtained for subsequent analysis. At the same time, Q20, Q30 and GC content were calculated for the clean data. All downstream analyses were performed using high quality (> Q30) reads. Reference genome and gene model annotation files were obtained from Yan et al. (2019) [28]. Reference genomes were indexed using Hisat2 v2.0.5 and paired clean reads were aligned to the reference genome using Hisat2 v2.0.5 [29]. The new transcripts were assembled by String Tie software [30], using a network flow algorithm together with optional de novo assembly, to assemble these complex data sets into transcripts.

Quantification of gene expression levels

The aims of this process were to count the number of reads for the full range of all genes using the Feature Counts tool in the subread 1.6.4 software [31], to match positions on the reference genome to genes, to filter out reads with low quality, to filter out reads from overlapping regions, and to match non-matched reads to multiple regions. The fragments per kilobase of transcript per million mapped reads (FPKM) is a common metric for estimating gene expression levels, and it is corrected for sequencing depth and sequencing length. Thus, the use of FPKM reduces the effect of gene length and sequencing depth on reads count, and it is generally considered to have a threshold value of 1. The resulting FPKM values were analyzed by principal component analysis (PCA). The expression of differential genes was clustered using the H-cluster method. To visualize the correlation of the samples, we created a correlation coefficient plot (Fig. 1B).

Fig. 1
figure 1

Assembly and information about the transcriptome data of the Manila clam. A Length distribution of unigenes after assembly; the abscissa represents the length range of unigenes, and the ordinate represents the number of unigenes corresponding to the length. B Visualization of Pearson correlation analysis of samples from the 13 developmental stages

Full size image

Differential gene expression analysis

Differential expression analysis of samples within groups was performed using DESeq2 software [32]. It uses a model based on a negative binomial distribution to provide statistical routines for differential expression of gene expression data. The p-values obtained are used to control for false discovery rates to reduce the proportion of false positives using the Benjamini-Hochberg method for p-value correction. P-value (padj) < 0. 05, |log2foldchange| > 0 is generally considered to be the criterion for differential genes.

Analysis of gene enrichment

We performed KEGG pathway enrichment analysis on the differential gene sets using clusterProfile (R 3.4.3) software [33], which is a comprehensive database that integrates genomic, chemical, and phylogenetic information. KEGG pathway enrichment is based on padj < 0. 05 as the threshold for significant enrichment. Gene Ontology (GO) functionally significant enrichment analysis and genomic background comparison of differential expression of genes in GO functionally significant enrichment terms reveals differential expression of genes significantly associated with biological function (http://www.geneontology.org/). The WGCNA (weighted gene co-expression network analysis) algorithm is a typical systems biology algorithm for constructing gene co-expression networks based on high-throughput gene messenger RNA (mRNA) expression data and is widely used in the international biomedical field. GO and KEGG, WGCNA enrichment. The analysis of DEGs is implemented in the cluster Profiler R package (3.4.4), where a correction for gene length bias was applied. GO terms and KEGG pathways were considered significantly enriched for DEGs if padj < 0.05 of DEGs.

Real-time quantitative PCR (RT-qPCR) for DEGs validation

To validate the RNA-Seq data, inhibitor of apoptosis protein (IAP) - Birc7-a, Caspase-7, Caspase-3, EGFR (Epidermal growth factor receptor), Grb2 (Growth factor receptor-bound protein 2) and MAPK signaling pathway gene - Map2k1, a total of six expressed genes (FPKM > 1) were selected for qRT-PCR. Among them, Caspase-7, Caspase-3 and Birc7-a were shown to play important roles in apoptosis [34, 35]. MAPK pathway was reported to have important roles in cell growth and differentiation [36], and EGFR signaling pathway plays important roles in physiological processes such as cell growth, proliferation and differentiation [37], and was detected in this study. Grb2, a growth-related candidate gene, was reported in Charybdis japonica [38]. Total RNA from FE, PB1, PB2, TC, EC, B, G, T, D, U, P, S, and J, samples of 30 mg was isolated using TRIzol with three biological replicates in each group. Primers were designed using Primer5 software [39] (Premier Biosoft International, Palo Alto, CA, USA) (Table S2). β-actin genes were used as an internal control for qRT-PCR analysis. Single-stranded cDNA was synthesized with a reverse transcription system using the reverse transcriptase kit “PrimeScript™ RT Kit” (TaKaRa Bio, Shiga, Japan) according to the manufacturer’s instructions.

qRT-PCR was performed using TB Green Premix ExTaqII (TaKaRa Biologicals, Shiga, Japan). The qRT-PCR was conducted in a final volume of 20 µL, which consisted of 1 µL of cDNA (50 ng/µL), 1 µL of forwarding primer (80 ng/µL), 1 µL of reverse primer (80 ng/µL), 7 µL of ddH2O, and 10 µL of TB Green, with polymerase activation at 94 ℃ for 5 min, followed by 40 cycles at 94 ℃ for 30 s, 60 ℃ for 30 s and 72 ℃ for 30 s. Each sample was processed in triplicate in the Roche LightCycler 480 Real Time PCR System (Roche, Basel, Switzerland). All data were analyzed using the 2−ΔΔCt method.

Results

Transcriptome sequencing

We subjected samples from 13 different developmental stages (FE, PB1, PB2, TC, EC, B, G, T, D, U, P, S, and J) of Manila clams to RNA-seq (Table S3). A total of 1,127,214,960 raw reads and 1,082,461,572 clean reads (96.03% of raw reads) were acquired from the 39 RNA-seq libraries. After raw data filtering, sequencing error rate checking, and GC content distribution analysis, clean reads for subsequent analysis were obtained. This work yielded 16.1 Gb of data with Q20 > 97.11% and Q30 > 91. 65% (Table 1), for a total of 50,285 genes for differentially expressed gene (DEG) analysis. The distribution of sequence lengths indicated that most single genes were more than 2000 base pairs in length; the second highest number of single genes were between 200 and 300 base pairs in length (Fig. 1A). All genes had between 30% and 40% GC content (Fig. S2). The sample correlation matrix based on Pearson’s correlation coefficients (Fig. 1B) showed satisfactory correlation between samples from different developmental periods and between samples from adjacent periods. The transcriptome correlations of the sample collections were analyzed by hierarchical clustering and principal component analysis (PCA). PCA results revealed that the samples (FE, PB1, PB2, EC, TC, P), (S, J), (U, D), and (B, G) differed significantly from each other (Fig. S3).

Table 1 Summary of RNA-seq analysis based on the R. philippinarum transcriptomes
Full size table

DEG analysis

A total of 50,284 genes were obtained after RNA sequencing, including 27,652 genes mapped to the Manila clam genome and 22,632 genes not mapped to the Manila clam genome. The mapped genes were annotated using the Pfam, GO and KEGG databases, resulting in 3,519 genes in the Pfam database, 23,196 genes in the GO database and 960 genes in the KEGG database. The heat map of all genes showed broadly consistent expression between biological duplicate samples, which implied high reliability of the sample (Fig. S4). To investigate the involvement of DEGs in the development of Manila clam, the expression levels of genes were compared between groups at different developmental stages. In total, 26,445 DEGs were obtained with |log2(FC)| ≥0 and p < 0.05. We identified significant differences in the number of DEGs for pair-wise comparisons between two adjacent periods for the 13 developmental stages. Notably, we detected more upregulated genes than downregulated genes in most samples, except PB1 vs. FE, PB2 vs. PB1, and P vs. U (Table S3). Volcano plots of the number of DEGs in the 12 comparison groups (FE vs. PB1, PB2 vs. PB1, TC vs. PB2, EC vs. TC, B vs. EC, G vs. B, T vs. G, D vs. T, U vs. D, P vs. U, S vs. P, and J vs. S) revealed a range of 114 to 13,840 significantly upregulated genes and 314 to 6604 (Fig. 2).

Fig. 2
figure 2

Visualization analysis of expression differences for different comparisons. DEG volcano plots of A FE vs. PB1; B PB1 vs. PB2; C PB2 vs. TC; D TC vs. EC; E EC vs. B; F: B vs. G; G: G vs. T; H: T vs. D; I: D vs. U; J: U vs. P; K: P vs. S; L: S vs. J. One dot in the volcano represents one gene, red dots represent upregulated genes, green dots represent downregulated genes, and blue dots indicate undifferentiated genes. p-value < 0.05, |log2FoldChange| > 0. The abscissa is log2FoldChange and the ordinate is –log10 (p value). The larger the value of –log10 (p value), the more significant the difference of the DEGs

Full size image

We performed GO and KEGG analyses to further evaluate the function of DEGs at different developmental stages. In total, 14,244 genes (28.33% of total genes) were annotated in the GO database, and 3250 genes (9.41% of total genes) were annotated in the KEGG database. The number of DEGs in the GO enrichment analysis was highly variable among the pairwise comparisons (Table 2). The number of items in the 13 samples that passed the GO terminology classification process ranged from 395 to 880. In addition, the KEGG annotation of unigenes from each group reveals that the number of pathways varies between 49 and 117compared with each group (Table 3), and a total of 4734 unigenes assigned to 181 pathways were identified (Table 3), with the number of unigenes in each pathway varying from 1 to 99.

Table 2 Number of DEGs in GO terms
Full size table
Table 3 Number of kegg pathway
Full size table

The key DEGs involved in different developmental stages

The early development of the Manila clam was divided into different period including cell division, cell differentiation, larval growth and development, larval metamorphosis, and shell calcification (Fig. 3). Genes involved in cytokinesis were upregulated during FE-EC, including Chst15, Cdk1, Cdc40, CDC7, Lyar, Cdk11b, Cdc42, SART1, CDC27, CDC5L, DHAR1, and KLHDC3. Most genes involved in cell differentiation were upregulated during the B-G-T period, including Foxj1b, SNRPA, FGF13, SMN1, DNAJB5, Dynein beta, DNAH10, RGS3, Ift22, and Ift88. During this period, the larvae began to develop organs, and the source of nutrition changed from the yolk provided by the egg to exogenous food. During the T-D-U process, the larval shell cuticle gradually calcified until the secondary shell appeared and we found that Aper1, Cht3, Col6a6, FGFR1, Dyrk2, Syncrip, and Serbp1 were upregulated during this process. During the U-P process, the larval veliger organ was shed and the lifestyle changed from planktonic to demersal. The differential genes at work during this change were 17-beta-HSD, Cytochrome b5, Scully-protein, Sperm-flagellar-protein2, IFT27, Gm166Spef2, TR-interacting-protein-3, TGF-beta, Fox-1, EOGT, BI-1, TMBIM4, Tsc6, CARP-1, PDZK1, and Caspase. Egfr, TGFBI1, TGFBI2, wnt8b, WNT5A, IGF2BP2, OGFRL1, KLF1, Klf3, KLF7, and Klf5 were upregulated during the P-J process. Growth and development-related genes, such as EGF1 and TGFBRAP1, were consistently upregulated during the early development of Manila clam larvae.

Fig. 3
figure 3

The heat map shows the expression levels of key DEGs for different molecular functions in samples from 13 different developmental periods. The color intensity from blue to red indicates the magnitude of differential expression. Red color indicates the mostly upregulated genes, and blue color indicates the mostly down regulated genes

Full size image

Figure S5A indicated the important role of chitin and matrix protein-related enzyme genes during Manila clam development, and we found some genes involving formation of protein complex were up-regulated during the middle and late stages of clam larvae development. The heatmap analysis revealed gene expressions of TGF-beta gene family during S and J stages were significantly up-regulated than other periods (Fig. S5B). In addition, apoptosis-related genes and Fox transcription factors, Tyr gene changed in expression during development of Manila clam (Fig. S5C-E).

GO analysis of DEGs

Following the analysis of DEGs, we conducted GO analysis and classified the DEGs into three main categories (Fig. 4): biological processes, cellular components, and molecular functions. Next, the distribution and enrichment of DEGs in GO these functional categories were analysed to identify genomes with differential expression levels (Fig. 4). The numbers of DEGs annotated in the 12 comparison groups FE vs. PB1, PB2 vs. PB1, TC vs. PB2, EC vs. TC, B vs. EC, G vs. B, T vs. G, D vs. T, U vs. D, P vs. U, S vs. P, and J vs. S were 428, 3468, 761, 1134, 10,054, 1301, 8366, 10,753, 852, 8843, 17,319, and 2214, respectively (Fig. 5). The GO terms significantly enriched in each comparison group were Rho protein signal transduction; intracellular membrane-bound organelles; protein-containing complexes; molecular function regulator; membrane-bounded organelles; regulation of nitrogen compound metabolic process; regulation of cellular macromolecule; chitin metabolic process, chitin binding; actin cytoskeleton, regulation of biological quality; aminoglycan metabolic process, oxidoreductase activity; hydrolase activity, acyl-CoA dehydrogenase activity; and proton-transporting, motor activity, respectively. The variety of biological processes identified in the B vs. EC comparison was the largest of all comparison groups, and the cellular components and molecular functions found in S vs. P were the most diverse. The results showed that the numbers of DEGs in the FE vs. PB1 and TC vs. PB2 comparison groups were relatively low, and the numbers of DEGs in the B vs. EC and G vs. B comparison groups were significantly decreased, and the number of DEGs in the P vs. U comparison group was significantly higher than that in the other groups.

Fig. 4
figure 4

The overall GO classification annotation of the unigenes

Full size image
Fig. 5
figure 5

The analysis of GO functional enrichment. The ordinate represents the top 30 enriched GO terms, and the Rich factor is the proportion of the DEGs annotated in the GO term to the total annotated genes in the GO term. The deeper red the color of the Q value, the more significant the enrichment of the GO term. The significant number of DEGs in the GO term is represented by the size of the circle

Full size image

KEGG analysis of DEGs

To understand the pathways at work at different developmental stages, we applied KEGG analysis to the DEGs (Table S4). The largest numbers of DEGs were annotated in the S vs. P comparison group, and relatively few DEGs were annotated in the FE vs. PB1 and TC vs. PB2 comparison groups. The highest numbers of pathways were found in the T vs. G, D vs. T, and S vs. P comparison groups. The lowest number of pathways was observed in the FE vs. PB1 comparison group, and the highest number of significantly enriched pathways was detected in the S vs. P comparison group.

The KEGG enrichment analysis results showed that a number of significantly enriched pathways played important roles in different developmental periods (Fig. 6). The KEGG enrichment analysis disclosed that eleven significantly enriched pathways (Pyruvate metabolism; Glycolysis/Gluconeogenesis; Fatty acid degradation; Oxidative phosphorylation; beta-Alanine metabolism; Phagosome; Carbon metabolism; Citrate cycle; Propanoate metabolism; Valine, leucine and isoleucine degradation; Ribosome) were observed in S vs. P comparison group, which is more than other comparison groups. The PB1 vs. PB2, EC vs. TC, B vs. EC, T vs. G, D vs. T, U vs. D, P vs. U, and J vs. S comparison groups had only four (Carbon metabolism; DNA replication; Endocytosis; mRNA surveillance pathway), one (DNA replication), six (Wnt signaling pathway; Foxo signaling pathway; Spliceosome; Biosynthesis of amino acids; Endocytosis; Protein processing in endoplasmic reticulum; Ribosome), six (Protein processing in endoplasmic reticulum; Ubiquitin mediated proteolysis; TGF-beta signaling pathway; DNA replication; Wnt signaling pathway; Ribosome ), six (Glutathione metabolism; Phagosome; Metabolism of xenobiotics by cytochrome P450; Drug metabolism - cytochrome P450; Oxidative phosphorylation; Ribosome), one(Phagosome), three (Glutathione metabolism; Phagosome; Ribosome) and two (ECM-receptor interaction; Phagosome) significantly enriched pathway, respectively.

Fig. 6
figure 6

Results of KEGG pathway functional enrichment analysis. The ordinate represents the top 30 enriched KEGG pathways, and the Rich factor is the proportion of the DEGs annotated in the KEGG pathway to the total annotated genes in the KEGG pathway. The deeper red the color of the Q value, the more significant the enrichment of the KEGG pathway. The significant number of DEGs in the KEGG pathway is represented by the size of the circle.KEGG is developed by Kanehisa Laboratories and the source was cited [40,41,42]

Full size image

Weighted gene co-expression network analysis

We conducted weighted co-expression network analysis to investigate patterns of gene association between samples from different developmental periods. A soft threshold power of six was selected. Genes with similar expression patterns were grouped into one module among a total of 11 modules (Fig. 7). Based on the module-sample correlation analysis (Fig. S6), most of the genes in the black module are up-regulated in expression from FE to EC stage and down-regulated in expression from B to J stage. Genes in the blue module are down-regulated in expression from FE to G stage and up-regulated between T to U stage and S to J. Interestingly, in the trend of up-regulation of gene expression in the late early developmental stage, these genes are down-regulated in the P stage. In the green module, genes are down-regulated from FE1 to EC and D to J, and up-regulated from B to T. In the red module, genes are down-regulated from FE to EC with P to J and up-regulated from B to U. In the yellow module, genes are down-regulated in the FE to T and P stage and up-regulated in the D to U and S to J stage (Fig. 8).

Fig. 7
figure 7

Clustering of module eigengenes and heat maps of modules in different development stages. The number in parentheses is the number of genes in each module. Blue indicates a negative correlation and red indicates a positive correlation

Full size image
Fig. 8
figure 8

Expression of candidate genes specific to the different developmental stages in the black, blue, green, red, and yellow modules. The figure above shows a heat map of all genes assigned to the black, blue, green, red, and yellow modules at different developmental stages

Full size image

The gene with the highest connectivity in the black module is RAD54B (DNA repair and recombination protein RAD54B, k-value = 306.3). The gene with the highest connectivity in the blue module is RPN2 (Dolichyl-diphosphooligosaccharide-protein glycosyltransferase subunit, k value = 330.02). The gene with high connectivity in the green module is Cherp (Calcium homeostasis endoplasmic reticulum protein, k value = 188.21). The gene with high connectivity in the red module is GTPBP2 (GTP-binding protein 2, k value = 318.69). The gene with high connectivity in the yellow module is ATP5A1 (ATP synthase subunit alpha, k value = 253.

RT-qPCR confirmation of RNA-seq data

To validate the sequencing data obtained from the 13 developmental period samples, we selected six genes that play an important role in the early development of Manila clams, including Mitogen-activated protein kinase (Map2k1), Epidermal growth factor receptor (Egfr), Baculoviral IAP repeat-containing protein 7-A (Birc7-a), Caspase-3 (Casp3), Growth factor receptor-bound protein 2 (Grb2), and Caspase-6 (Casp6). The trends for these genes detected by RT-qPCR were consistent with the RNA-seq data for the 13 development stages (Fig. 9), which confirmed the authenticity of our RNA-seq data.

Fig. 9
figure 9

The relative mRNA expression pattern of different DEGs depending on different molecular functions at different development stages detected by qRT-PCR. The values are given in terms of relative mRNA expression. Data are presented as the means of three replicates ± standard error

Full size image

Discussion

Shell calcification

During Manila clam development, calcification of the shell begins after the D- larval stage and the completion of the secondary shell during early development. A good quality secondary shell provides protection and helps the clams survive. However, shell calcification is subject to a combination of intrinsic and external factors [43]. In this study, we identified several DEGs that were involved in shell matrix formation, such as collagen [44] (a basic component of the extracellular matrix), chitin-binding proteins [45], and shell matrix proteins. Structural domain analysis showed that chitin in samples from different developmental periods contained multiple structural domains, and 37 members of the chitin family were found in Manila clam. In thecurrent study, 10 genes including chitin-related enzymes were significantly upregulated at the D vs. U and S vs. J stages (Fig. S5A). In mussels Mytilus edulis, chitin-related enzymes have been shown to play an important role in biochemical defense processes [46]. In a study on mollusk biomineralization, it was reported that there are many tyrosinase genes in the genome of bivalves and qPCR results showed high expression of tyrosinase genes in the mantle tissue of the Pacific oyster C. gigas [47]. In the present study, 15 tyrosinase genes were upregulated during Manila clam D, U, S, J and shell morphology was altered from D to U in Manila clam (Fig. S5E), indicating tyrosinase genes are likely to play important roles in the larvae shell growth in Manila clam.

Larval metamorphosis

Larval metamorphosis is a complex biological process [48] that involves morphological and lifestyle changes [49]. During the pediveliger stage of clam larvae, the upper disc organs are shed and the lifestyle changes from planktonic to demersal; however, the detailed molecular mechanisms at work during this process are poorly understood. NO is thought to play an important role in larval metamorphosis [50], because habitat change leads to a starvation state, which inhibits NO signaling, which in turn triggers larval metamorphosis [51]. The second messenger NO plays an important role in the regulation of larval metamorphosis in marine invertebrates [52]. NO regulates larval settlement by mediating downstream calcitonin gene-related peptide signaling [53], reducing MAPK expression, and affecting extracellular signal-regulated kinase (ERK signaling), This explains the down-regulation of the MAP2K1 gene, a key gene in the MAPK pathway, during the T-U period in the results of this study (Fig. 9). Additionally, ERK phosphorylation levels are required for transcription of key degenerate genes [54]. In bivalves, such as C. gigas and M. coruscus, NO has been shown to be a negative regulator of larval metamorphosis [55]. An important role for γ-aminobutyric acid in larval metamorphosis has also been reported [56], and γ-aminobutyric acid functions internally as a neurotransmitter to control the initiation of metamorphosis.

Transforming growth factor-beta (TGF-β)

TGFs constitute a superfamily of growth factors with important functions in cell growth, proliferation, differentiation, and apoptosis [57]. Activation of the TGF-β receptor is followed by intracellular effector-Smad protein translocation into the nucleus [58, 59]. In our study, three TGF-βs were significantly upregulated during the B and G stages of Manila clam development, which indicated an important role for these genes during the period of cell differentiation [60]. TGF family analysis revealed that Manila clams contained nine TGF superfamily members, and these genes were expressed in different patterns in the 13 developmental stages (Fig. S5B). An important role of TGF-β in the development of the Pacific oyster gonad has been reported, as its specific expression in gonads peaked when the gonads reached full maturity [61]. Two TGF-β proteins in Drosophila act as heterodimers to stimulate male germ cell proliferation [62]. In Tilapia (Oreochromis mossambicus), TGF-β ligands exhibit tissue specificity in the gonads, and the expression profiles of eight genes in the gonads at different stages of development suggest an important role for TGF-β ligands in sex determination and reproduction [63]. However, the role of TGF-β in the growth of Manila clams requires further study.

Apoptosis

The apoptotic program is a highly conserved process that is widespread in nature [64], and it is thought to be closely related to biological evolution. During sea urchin embryonic development, apoptotic cells are found in the flagellar zone prior to larval metamorphosis [65]. It has been reported that apoptosis is induced in larval epithelial cells during accelerated larval metamorphosis [66]. Caspase is closely associated with apoptosis, and caspases were shown to play a key role in foot loss during larval metamorphosis of Fujian oyster (C. angulata) [67]. In this study, apoptosis-inducing factors and apoptosis protein-related genes were consistently expressed at almost every period (Fig. 9, Fig. S5D), suggesting that apoptosis factors were involved in the early development of Manila clams. Apoptosis is also an important part of the host defense system [68]. However, in the paralytic shellfish toxin induced hemocytes death process, the caspase-dependent pathway was found to induce apoptosis in Pacific oyster immune cells, and the resistance of the oysters to microbial infection diminished [69]. In our study, we identified 47 members of the Manila clam caspase structural domain using hmm search. A study in which pre-metamorphic M. meretrix larvae were treated with Caspase inhibitors showed that apoptosis may be the main mechanism of veliger degeneration during metamorphosis, and in situ labelling showed that caspases were involved in the morphological change process of the larvae [70]. Thus, the role of caspase in the metamorphosis of Manila clam larvae is worthy of further study.

FOXO transcription factor

We found that the FOXO pathway was significantly expressed during the EC to B transition, when the number of upregulated DEGs was 37 and the number of downregulated DEGs was nine. The results indicated that FOXO transcription factors were involved in the EC-B process, during which cell proliferation mainly occurred. FOXO proteins were previously reported to be involved in cellular processes, including cell proliferation and apoptosis [71]. FOXO transcription factors also function as transcriptional activators, and their activity is inhibited by insulin and growth factor signaling [72]. This behaviour explains the observed upregulation of FOXO transcription factor expression during the blastula stage and the gastrula stage. FOXO factors play a role in regulating cell clearance mechanisms, and FOXO3 induces the expression of several autophagy genes involved in the coordination of autophagy-related gene networks [73]. FOXO factors also regulate the ubiquitin-proteasome system [74]. Thus, FOXO is likely to play important roles in the maintenance of cellular homeostasis during the early cell division period in Manila clams.

Glutathione metabolism pathway

Glutathione is a common antioxidant enzyme in bivalves and is used as an important metabolite for protection against oxidative stress in cells [75]. In this study, the glutathione metabolic pathway was significantly enriched during the P. vs. U period versus the T. vs. G. G to T period. Manila clam larvae undergo organogenesis. Some studies have revealed the role of glutathione as more than an antioxidant; the ability of cysteine residues of GSH to be reversibly oxidized in various protein targets acts as a mediator of many processes [76] such as cell proliferation, cell differentiation, and cell death [77, 78]. These studies explain the significant enrichment of the glutathione metabolic pathway during the T to G period, when manila calm undergoes cell differentiation and forms larvae. Glutathione has been found to promote growth in aquatic animals, such as pearl oyster (P. fucata martensii) [79], and white shrimp (Litopenaeus vannamei) [80].

DNA repair and recombination protein RAD54B

Weighted gene co-expression network analysis showed that RAD54B was the gene with the highest k-value in the black module and, in general, genes with the highest connectivity (k-value) ranking can be considered hub genes [81]. Rad54B is a DNA-dependent ATPase that has been described as being involved in recombinant repair of DNA damage [82]. Many RAD homologs for longevity, neurological and immune adaptation have been identified in the lobster genome, including seven homologs of RAD54 [83]. This implies a potential role for RAD genes in biological stress resistance. In the present study, the RAD54B gene was up regulated from FE to EC, indicating that the RAD54B gene has an important role in fertilized egg development in the Manila clam, but the exact function needs further study.

Conclusion

In this study, we analyzed 13 stages of early development in the Manila clam by RNA-seq, including Fertilized egg, 1 st polar body, 2 st polar body, Two-cell, Eight-cell, Blastula, Gastrula, Trochophora, D-larva, Umbo-veliger, Pediveliger, Single pipe juvenile, and Juvenile. We found that a unique gene expression pattern in each stage, and similar expression profile was observed in adjacent stages. The development process of Manila clam is broadly divided into cell division, cell differentiation, larval growth and development, larval metamorphosis, and shell calcification. The rho protein signaling, membrane-bound organelle, were prominent at the early development of Manila clam. KEGG analysis showed that the stage of Manila clam cell differentiation (EC-B) was associated with the Foxo pathway and Manila clam hatching (G-T) was associated with the TGF-beta signaling pathway. In addition, the stage of shell alteration (T-D) is associated with Tyrosine metabolism. These data provide insights into the molecular mechanisms of early Manila clam development and may help to enhance Manila clam reproduction and artificial breeding.

Availability of data and materials

All data generated or analyzed during this study are included in this published article. The raw sequences for R. philippinarum have been deposited in the NCBI PRJNA808620 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA808620).

References

  1. Plazzi F, Ceregato A, Taviani M, et al. A Molecular Phylogeny of Bivalve Mollusks: Ancient Radiations and Divergences as Revealed by Mitochondrial Genes. PLoS ONE. 2011;6(11):e27147. https://doi.org/10.1371/journal.pone.0027147.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. Rdiger B, Mikkelsen P M. Bivalvia–a look at the Branches. Zool J Linnean Soc. 2010;(3):223–35. https://doi.org/10.1111/j.1096-3642.2006.00255.x.

  3. Vermeij S. The Plankton and the Benthos: Origins and Early History of an Evolving Relationship. Paleobiol. 1994;20(3):297–319. https://doi.org/10.1017/S0094837300012793.

    Article  Google Scholar 

  4. Wray G A. Evolution of larvae and developmental modes. 1995.

  5. ERNST, KNIPRATH. Ontogeny of the Molluscan Shell Field: a Review. Zoologica Scripta, 1981. https://doi.org/10.1111/j.1463-6409.1981.tb00485.x.

    Article  Google Scholar 

  6. Jin S, Zhang Y, Thiyagarajan V, et al. Protein expression during the embryonic development of a gastropod. Proteomics. 2010;10(14). https://doi.org/10.1002/pmic.200900846.

  7. Wanninger A, Haszprunar G. Chiton myogenesis: Perspectives for the development and evolution of larval and adult muscle systems in molluscs. J Morphol. 2002;251(2):103-113. https://doi.org/10.1002/jmor.1077.

    Article  PubMed  Google Scholar 

  8. Wanninger A, Haszprunar G. The expression of an engrailed protein during embryonic shell formation of the tusk-shell, Antalisentalis (Mollusca, Scaphopoda). Evol Develop. 2010;3(5):312–21. https://doi.org/10.1046/j.1525-142X.2001.01034.x.

  9. Huan P, Liu G, Wang H, et al. Identification of a tyrosinase gene potentially involved in early larval shell biogenesis of the Pacific oyster Crassostrea gigas. Develop Genes Evol. 2013;223(6):389–94. https://link.springer.com/article/10.1007/s00427-013-0450-z.

  10. Liu Z, Wang L, Yan Y, et al. D1 dopamine receptor is involved in shell formation in larvae of Pacific oyster Crassostrea gigas. Develop Comparative Immunol. 2018;84:337–42. https://doi.org/10.1016/j.dci.2018.03.009.

    Article  CAS  Google Scholar 

  11. Wang X, Liu B, Xiang J. Cloning, characterization and expression of ferritin subunit from clam Meretrix meretrix in different larval stages. Comparative Biochem Physiol Part B Biochem Mol Biol. 2009;154(1):12-16. https://doi.org/10.1016/j.cbpb.2009.04.011.

    Article  CAS  Google Scholar 

  12. Crisp D J. Factors influencing the settlement of marine invertebrate larvae. Chemoreception in Marine Organisms, 1974.https://doi.org/10.1016/0022-0981(91)90068-8.

    Article  Google Scholar 

  13. Clara Sánchez-Lazo, Martinez-Pita I, Young T, et al. Induction of settlement in larvae of the mussel Mytilus galloprovincialis using neuroactive compounds. Aquaculture. 2012;344–349:210–215. https://doi.org/10.1016/j.aquaculture.2012.03.021.

  14. Zhang J, Xiong X, Deng Y, et al. Integrated application of transcriptomics and metabolomics provides insights into the larval metamorphosis of pearl oyster (Pinctada fucatamartensii). Aquaculture. 2021;532:736067. https://doi.org/10.1016/j.aquaculture.2020.736067.

    Article  CAS  Google Scholar 

  15. Li H, Yu H, Li Q. Striated myosin heavy chain gene is a crucial regulator of larval myogenesis in the pacific oyster Crassostrea gigas. Int J Biol Macromolecules. 2021;179(1). https://doi.org/10.1016/j.ijbiomac.2021.03.022.

  16. Zeng Zhen, Ni Jianbin, Han Guo Dong, et al. In situ hybridization of genes related to the metabolism and development of the Portuguese oyster Crassostrea angulata/ Congress and Symposium of the Shellfish Branch of the Chinese Society of Marine Lakes and Marshes, Chinese Zoological Society. 2011.

  17. Huan P, Wang H, Liu B. Transcriptomic Analysis of the Clam Meretrix meretrix on Different Larval Stages. Marine Biotechnol. 2012;14(1):69-78. https://doi.org/10.1007/s10126-011-9389-0.

    Article  CAS  Google Scholar 

  18. Lelong C, Mathieu M, Favrel P. Structure and expression of mGDF, a new member of the transforming growth factor-beta superfamily in the bivalve mollusc Crassostrea gigas. Eur J Biochem. 2010;267(13):3986-3993. https://doi.org/10.1046/j.1432-1327.2000.01432.x.

    Article  Google Scholar 

  19. Susanne V, Penny ME, Li X, et al. First report of a putative involvement of the NMDA pathway in Pacific oyster (Crassostrea gigas) development: Effect of NMDA receptor ligands on oyster metamorphosis with implications for bivalve hatchery management. Aquaculture. 2018;497:140-146. https://doi.org/10.1016/j.aquaculture.2018.07.048.

    Article  CAS  Google Scholar 

  20. Joyce A, Vogeler S. Molluscan bivalve settlement and metamorphosis: Neuroendocrine inducers and morphogenetic responses. Aquaculture. 2018:64–82. https://doi.org/10.1016/j.aquaculture.2018.01.002

  21. FAO (Food & Agriculture Organization), Fishery and Aquaculture Statistics 2017, Food and Agriculture Organization of the United Nations, Rome, 2019.

    Google Scholar 

  22. Watanabe Y, Lo N. Larval Production and Mortality of Pacific Saury, Cololabissaira, in the Northwestern Pacific Ocean. Fishery Bull- Natl Ocean Atmospheric Administration. 1989;87(3):601-613. https://doi.org/10.1093/icesjms/fst089.

    Article  Google Scholar 

  23. Ch’Ng CL, Senoo S. Egg and larval development of a new hybrid grouper, tiger grouper Epinephelusfuscoguttatus × giant grouper E. lanceolatus. Aquaculture Science. 2008;56(4):505-512. https://doi.org/10.11233/aquaculturesci.56.505.

    Article  Google Scholar 

  24. Chong YH, Moriya S, Ogawa S, et al. Deep Brain Photoreceptor (val-opsin) Gene Knockout Using CRISPR/Cas Affects Chorion Formation and Embryonic Hatching in the Zebrafish. Plos One. 2016;11(10):e0165535. https://doi.org/10.1371/journal.pone.0165535.

    Article  CAS  Google Scholar 

  25. Nie H, Zheng M, Wang Z, et al. Transcriptomic analysis provides insights into candidate genes and molecular pathways involved in growth of Manila clam Ruditapesphilippinarum. Funct Integrative Genomics. 2021;21(3–4):1–13. https://doi.org/10.1007/s10142-021-00780-1.

    Article  CAS  Google Scholar 

  26. Nie H, Yan X, Huo Z, et al. Construction of a High-Density Genetic Map and Quantitative Trait Locus Mapping in the Manila clam. Sci Rep. 2017;7(1):229. https://doi.org/10.1038/s41598-017-00246-0.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Simons A. A quality control tool for high throughput sequence data. A quality control tool for high throughput sequence data. 2010. https://doi.org/10.12688/f1000research.21142.2

  28. Yan X, Nie H, Huo Z, et al. Clam genome sequence clarifies the molecular basis of its benthic adaptation and extraordinary shell color diversity. Iscience. 2019;19:1225-1237. https://doi.org/10.1016/j.isci.2019.08.049.

    Article  PubMed  PubMed Central  Google Scholar 

  29. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nature Methods. 2015;12(4):357–60. https://doi.org/10.1038/nmeth.3317.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Mihaela Pertea, Geo M Pertea, Corina M Antonescu1, et al. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5.

  31. Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30(7):923–30. https://doi.org/10.1093/bioinformatics/btt656.

    Article  CAS  PubMed  Google Scholar 

  32. Costa-Silva J, Domingues D, Lopes FM. RNA-Seq differential expression analysis: An extended review and a software tool. PloS one. 2017;12(12):e0190152. https://doi.org/10.1371/journal.pone.0190152.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Chen X, Sun M. Identification of key genes, pathways and potential therapeutic agents for IgA nephropathy using an integrated bioinformatics analysis. J Renin-Angiotensin-Aldosterone Syst. 2020;21(2):1470320320919635. https://doi.org/10.1177/1470320320919635.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Dubrez-Daloz L, Dupoux A, Cartier J. IAPs: more than just inhibitors of apoptosis proteins. Cell Cycle. 2008;7(8):1036–46. https://doi.org/10.4161/cc.7.8.5783.

    Article  CAS  PubMed  Google Scholar 

  35. Wu MH, Jin XK, Yu AQ, et al. Caspase-mediated apoptosis in crustaceans: cloning and functional characterization of EsCaspase-3-like protein from Eriocheir sinensis. Fish Shellfish Immunol. 2014;41(2):625–32. https://doi.org/10.1016/j.fsi.2014.10.017.

    Article  CAS  PubMed  Google Scholar 

  36. Iakovleva NV, Gorbushin AM, Zelck UE. Partial characterization of mitogen-activated protein kinases (MAPK) from haemocytes of the common periwinkle, Littorina littorea (Gastropoda: Prosobranchia). Fish Shellfish Immunol. 2006;20(4):665–8. https://doi.org/10.1016/j.fsi.2005.08.001.

    Article  CAS  PubMed  Google Scholar 

  37. Wang Q, Hao R, Zhao X, Bioscience, et al. Identification of EGFR in pearl oyster (Pinctada fucatamartensii) and correlation analysis of its expression and growth traits. Biotechnol Biochem. 2018;82(7):1073-1080. https://doi.org/10.1080/09168451.2018.1459174.

    Article  CAS  Google Scholar 

  38. Lou F, Yang T, Han Z, et al. Transcriptome analysis for identification of candidate genes related to sex determination and growth in Charybdis japonica. Gene. 2018;677:10–6. https://doi.org/10.1016/j.gene.2018.07.044.

    Article  CAS  PubMed  Google Scholar 

  39. Lalitha S. Primer premier 5. Biotech Software & Internet Report. Comput Software J Sci. 2000;1(6):270–2. https://doi.org/10.1089/152791601750294344.

  40. Kanehisa, M. and Goto, S.; KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 2000, 28, 27–30.

    Article  CAS  Google Scholar 

  41. Kanehisa, M; Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019, 28, 1947–1951.

    Article  CAS  Google Scholar 

  42. Kanehisa M, Furumichi M, Sato Y, Ishiguro-Watanabe M, Tanabe M. KEGG: integrating viruses and cellular organisms. Nucleic Acids Res. 2021;49:D545–51.

  43. Waldbusser GG, Bergschneider H, Green MA. Size-dependent pH effect on calcification in post-larval hard clam Mercenaria spp. Marine Ecol Progress. 2010;417:171–82. https://doi.org/10.3354/meps08809.

    Article  Google Scholar 

  44. Lang R P. Identification of candidate genes for survival and their use in predicting field performance of Pacific oyster Crassostrea gigas families in coastal waters. [D]. Oregon State University, 2008. https://doi.org/10.1002/ece3.719.

    Article  Google Scholar 

  45. Nakayama S, Suzuki M, Endo H, et al. Identification and characterization of a matrix protein (PPP-10) in the periostracum of the pearl oyster, Pinctada fucata. FEBS Open Bio, 2013, 3(1):421–427. https://doi.org/10.1016/j.fob.2013.10.001.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Arivalagan J, Marie B, Chiappetta G, et al. Deciphering shell proteome within different Baltic populations of mytilidmussels illustrates important local variability and potential consequences in the context of changing marine conditions. Sci Total Environment. 2020;745:140878. https://doi.org/10.1016/j.scitotenv.2020.140878.

    Article  CAS  Google Scholar 

  47. Aguilera F. Investigation of gene family evolution and the molecular basis of shell formation in molluscs. 2015. https://doi.org/10.14264/uql.2015.396.

  48. Peters-Didier J, Sewell MA. The role of the hyaline spheres in sea cucumber metamorphosis: lipid storage via transport cells in the blastocoel. EvoDevo. 2019;10(1):1–12. https://doi.org/10.1186/s13227-019-0119-4.

    Article  Google Scholar 

  49. Delalio LJ, Dion SM, Bootes AM, et al. Direct effects of hypoxia and nitric oxide on ecdysone secretion by insect prothoracic glands. J Insect Physiol. 2015;76:56–66. https://doi.org/10.1016/j.jinsphys.2015.02.009.

    Article  CAS  PubMed  Google Scholar 

  50. Vogeler S, Carboni S, Li X, et al. Bivalves are NO different: nitric oxide as negative regulator of metamorphosis in the Pacific oyster, Crassostrea gigas. BMC Develop Biol. 2020;20(1). https://doi.org/10.1186/s12861-020-00232-2

  51. Romero MR, Phuong MA, Bishop C, et al. Nitric oxide signaling differentially affects habitat choice by two larval morphs of the sea slug Alderiawillowi: mechanistic insight into evolutionary transitions in dispersal strategies. J Experiment Biol. 2012;216(6):1114–25. https://doi.org/10.1242/jeb.080747.

    Article  CAS  Google Scholar 

  52. Kitamura H, Kitahara S, Koh HB. Induction of Larval Settlement and Metamorphosis in the Sea Urchins Pseudocentrotusdepressus and Anthocidariscrassispina by Fatty Acids. Fish Sci. 2008;60(3):p311-313. https://doi.org/10.2331/fishsci.60.311.

    Article  Google Scholar 

  53. Zhang Y, He LS, Zhang G, et al. The regulatory role of the NO/cGMP signal transduction cascade during larval attachment and metamorphosis of the barnacle Balanus (= Amphibalanus) amphitrite. J Experiment Biol. 2012;215(Pt 21):3813–22. https://doi.org/10.1242/jeb.070235.

    Article  CAS  Google Scholar 

  54. Immacolata, Castellano, Elena, et al. Nitric Oxide Affects ERK Signaling through Down-Regulation of MAP Kinase Phosphatase Levels during Larval Development of the Ascidian Ciona intestinalis. PLoS ONE. 2014;9(7). https://doi.org/10.1371/journal.pone.0102907.

  55. Zhu Y T, Zhang Y, Liu Y Z, et al. Nitric Oxide Negatively Regulates Larval Metamorphosis in Hard-Shelled Mussel (Mytilus coruscus). Frontiers in Marine Science. 2020;7. https://doi.org/10.3389/fmars.2020.00356.

  56. Biscocho D, Cook JG, Long J, Shah N, Leise EM. GABA is an inhibitory neurotransmitter in the neural circuit regulating metamorphosis in a marine snail. Dev Neurobiol. 2018;78(7):736–53. https://doi.org/10.1002/dneu.22597.

    Article  CAS  PubMed  Google Scholar 

  57. Derynck R, Zhang YE. Smad-dependent and Smad-independent pathways in TGF-beta family signalling. Nature. 2003;425(6958):577–84. https://doi.org/10.1002/jcp.30529.

    Article  CAS  PubMed  Google Scholar 

  58. Shi Y, Joan Massagué. Mechanisms of TGF-beta signaling from cell membrane to the nucleus. Cell. 2003;113(6):685–700. https://doi.org/10.1016/S0092-8674(03)00432-X.

    Article  CAS  PubMed  Google Scholar 

  59. Lin HM, Lee JH, Yadav H, et al. Transforming Growth Factor-β/Smad3 Signaling Regulates Insulin Gene Transcription and Pancreatic Islet β-Cell Function. J Biol Chem. 2009;284(18):12246–57. https://doi.org/10.1074/jbc.M805379200.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Umeda M, Yoshida N, Hisada R, et al. ADAM9 enhances Th17 cell differentiation and autoimmunity by activating TGF-β1. Proceedings of the National Academy of Sciences, 2021, May 4;118(18): e2023230118. https://doi.org/10.1073/pnas.2023230118.

  61. Naimi A, Martinez AS, Specq ML, et al. Molecular cloning and gene expression of Cg-Foxl2 during the development and the adult gametogenetic cycle in the oyster Crassostrea gigas. Comp Biochem Physiol B Biochem Mol Biol. 2009;154(1):134-142. https://doi.org/10.1016/j.cbpb.2009.05.011.

    Article  CAS  PubMed  Google Scholar 

  62. Coon S L, Fitt W K, Bonar D B. Competence and delay of metamorphosis in the Pacific oyster Crassostrea gigas. Marine Biology, 1990, 106(3):379–387.

    Article  Google Scholar 

  63. Shuqing Z, Juan L, Zhilong L, et al. Identification and Evolution of TGF-β Signaling Pathway Members in Twenty-Four Animal Species and Expression in Tilapia. Int J Mol Sci. 2018;19(4):1154. https://doi.org/10.3390/ijms19041154.

    Article  CAS  Google Scholar 

  64. DaVies EL, Fuller MT. Regulation of Self-renewal and Differentiation in Adult Stem Cell Lineages: Lessons from the Drosophila Male Germ Line. Cold Spring HarbSymp Quant Biol. 2008;73:137-145. https://doi.org/10.1101/sqb.2008.73.063.

    Article  CAS  Google Scholar 

  65. Agnello M, Roccheri MC. Apoptosis: focus on sea urchin development. Apoptosis Int J Programmed Cell Death. 2010;15(3):322. https://doi.org/10.1007/s10495-009-0420-0.

    Article  Google Scholar 

  66. Roccheri MC, Tipa C, Bonaventura R, et al. Physiological and induced apoptosis in sea urchin larvae undergoing metamorphosis. Int J Develop Biol. 2002;46(6):801-806. https://doi.org/10.1079/IVP2002312.

    Article  Google Scholar 

  67. Yang B, Li L, Pu F, et al. Molecular cloning of two molluscan caspases and gene functional analysis duringCrassostrea angulate (Fujian oyster) larval metamorphosis. Mol Biol Rep. 2015;42(5):963-975. https://doi.org/10.1007/s11033-014-3833-y.

    Article  CAS  PubMed  Google Scholar 

  68. Li Y, Zhang L, Qu T, et al. Conservation and divergence of mitochondrial apoptosis pathway in the Pacific oyster, Crassostrea gigas. Cell Death Dis. 2017;8(7):e2915. https://doi.org/10.1038/cddis.2017.307.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  69. Abi-Khalil C, Finkelstein DS, Conejero G, et al. The paralytic shellfish toxin, saxitoxin, enters the cytoplasm and induces apoptosis of oyster immune cells through a caspase-dependent pathway. Aquatic Toxicol. 2017;190:133-141. https://doi.org/10.1016/j.aquatox.2017.07.001.

    Article  CAS  Google Scholar 

  70. Wang X, Liu B, Xiang J. Apoptosis and role of caspase in development of Meretrix meretrix larvae. Oceanol Limnol Sin. 2009;40:181–6. https://doi.org/10.3321/j.issn:0029-814X.2009.02.012.

    Article  CAS  Google Scholar 

  71. Mammucari C, Milan G, Romanello V, et al. FoxO3 controls autophagy in skeletal muscle in vivo. Cell Metab. 2007;6(6):458-471. https://doi.org/10.1016/j.cmet.2007.11.001.

    Article  CAS  PubMed  Google Scholar 

  72. Matsumoto M, Han S, Kitamura T, et al. Dual role of transcription factor FoxO1 in controlling hepatic insulin sensitivity and lipid metabolism. J Clin Invest. 2006;116(9):2464–72. https://doi.org/10.1172/JCI27047.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Tzivion G, Dobson M, Ramakrishnan G. FoxO transcription factors; Regulation by AKT and 14-3-3 proteins. Biochimica et Biophysica Acta (BBA)-Mol Cell Res. 2011;1813(11):1938–45. https://doi.org/10.1016/j.bbamcr.2011.06.002

  74. Webb AE, Brunet A. FOXO transcription factors: key regulators of cellular quality control. Trends Biochem Sci. 2014;39(4):159–69. https://doi.org/10.1016/j.tibs.2014.02.003.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Doyotte A, Cossu C, Jacquin MC, et al. Antioxidant enzymes, glutathione and lipid peroxidation as relevant biomarkers of experimental or field exposure in the gills and the digestive gland of the freshwater bivalve Unio tumidus. Aquatic Toxicol. 1997;39(2):93-110. https://doi.org/10.1016/S0166-445X(97)00024-6.

    Article  CAS  Google Scholar 

  76. Hansen J M, Harris C. Glutathione during embryonic development. Biochimica et Biophysica Acta (BBA)-General Subjects. 2015;1850(8):1527–42. https://doi.org/10.1016/j.bbagen.2014.12.001

  77. Laborde E. Glutathione transferases as mediators of signaling pathways involved in cell proliferation and cell death. Cell Death Differentiation. 2010;17(9):1373-80. https://doi.org/10.1038/cdd.2010.80

  78. Ashtiani HRA, Bakhshandi AK, Rahbar M, et al. Glutathione, cell proliferation and differentiation. Afr J Biotechnol. 2011;10(34):6348–63.

    CAS  Google Scholar 

  79. Hao R, Du X, Yang C, et al. Integrated application of transcriptomics and metabolomics provides insights into unsynchronized growth in pearl oyster Pinctada fucatamartensii. Sci Total Environ. 2019;666:46-56. https://doi.org/10.1016/j.scitotenv.2019.02.221.

    Article  CAS  PubMed  Google Scholar 

  80. Wang X, Xu W, Zhou H, et al. Reduced glutathione supplementation in practical diet improves the growth, anti-oxidative capacity, disease resistance and gut morphology of shrimp Litopenaeusvannamei. Fish Shellfish Immunol. 2018;73:152-157. https://doi.org/10.1016/j.fsi.2017.11.043.

    Article  CAS  PubMed  Google Scholar 

  81. Liesecke F, De Craene J O, Besseau S, et al. Improved gene co-expression network quality through expression dataset down-sampling and network aggregation. Sci Rep. 2019;9(1):1–16. https://doi.org/10.1038/s41598-019-50885-8

  82. Murzik U, Hemmerich P, Weidtkamp-Peters S, et al. Rad54B targeting to DNA double-strand break repair sites requires complex formation with S100A11. Mol Biol Cell. 2008;19(7):2926–35. https://doi.org/10.1091/mbc.e07-11-1167.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  83. Polinski JM, Zimin AV, Clark KF, et al. The American lobster genome reveals insights on longevity, neural, and immune adaptations. Sci Adv. 2021;7(26):eabe8290. https://doi.org/10.1126/sciadv.abe8290.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We thank Mr. Ning Li and Mr. Zhengxing Wang for their help in the sample collection.

Funding

This work was supported by the Chinese Ministry of Science and Technology through the National Key Research and Development Program of China (2018YFD0901400, 2019YFD0900704) and supported by China Agriculture Research System of MOF and MARA. The Project is sponsored by “Liaoning BaiQianWan Talents Program”. The work was supported by the Outstanding Chinese and Foreign Youth Exchange Program of China Association for Science and Technology (CAST, 2019), the Dalian Science and Technology Innovation Fund Project (2022), and the Scientific Research funding from Liaoning Provincial Department of Education (LJKZ0701).

Author information

Authors and Affiliations

  1. College of Fisheries and Life Science, Dalian Ocean University, 116023, Dalian, China

    Yanming Zhang, Hongtao Nie, Zhihui Yin & Xiwu Yan

  2. Engineering Research Center of Shellfish Culture and Breeding in Liaoning Province, College of Fisheries and Life Science, Dalian Ocean University, 116023, Dalian, China

    Yanming Zhang, Hongtao Nie, Zhihui Yin & Xiwu Yan

Authors
  1. Yanming Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Hongtao Nie
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Zhihui Yin
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Xiwu Yan
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

HN and XY conceived the experimental design, under which tissue samples featured herein were collected. HN and YZ conceived the transcriptomic study. YZ and ZY carried out the laboratory bench work, analyzed the data and wrote the manuscript. HN and YZ contributed substantially to revised drafts of the manuscript. The author(s) read and approved the final manuscript.

Corresponding author

Correspondence to Hongtao Nie.

Ethics declarations

Ethics approval and consent to participate

The Manila clam is not an endangered or protected species, so no specific permits were required for the study. Ethics approval was not required for this study.

Consent for publication

Not applicable.

Competing interests

We declare that we have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhang, Y., Nie, H., Yin, Z. et al. Comparative transcriptomic analysis revealed dynamic changes of distinct classes of genes during development of the Manila clam (Ruditapes philippinarum). BMC Genomics 23, 676 (2022). https://doi.org/10.1186/s12864-022-08813-0

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-022-08813-0

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us