Clam Genome Sequence Clarifies the Molecular Basis of Its Benthic Adaptation and Extraordinary Shell Color Diversity.

Ruditapes philippinarum is an economically important bivalve with remarkable diversity in its shell coloration patterns. In this study, we sequenced the whole genome of the Manila clam and investigated the molecular basis of its adaptation to hypoxia, acidification, and parasite stress with transcriptome sequencing and an RNA sequence analysis of different tissues and developmental stages to clarify these major issues. A number of immune-related gene families are expanded in the R. philippinarum genome, such as TEP, C3, C1qDC, Hsp70, SABL, and lysozyme, which are potentially important for its stress resistance and adaptation to a coastal benthic life. The transcriptome analyses demonstrated the dynamic and orchestrated specific expression of numerous innate immune-related genes in response to experimental challenge with pathogens. These findings suggest that the expansion of immune- and stress-related genes may play vital roles in resistance to adverse environments and has a profound effect on the clam's adaptation to benthic life.

Introduction

Bivalves are represented by approximately 20,000 species, distributed throughout aquatic habitats (Daniel, 2013). The Manila clam, Ruditapes philippinarum, is a bivalve mollusk with a worldwide distribution and important commercial, culinary, and ecological value. The production of R. philippinarum reached over 4.0 million tons, equivalent to 3.7 million USD, in 2015 (FAO, 2017). As a filter-feeding animal living buried in sediment, R. philippinarum is an important “sentinel” species used to assess the quality of the marine environment (Milan et al., 2011). Benthic bivalves play important roles in food chains, including as food for humans, and most benthic bivalves are critical in bioturbation, bioirrigation, and the breakdown of organic matter (Norkko and Shumway, 2011). However, massive mortality in the Manila clam has been reported arising from parasitic infections, which has severely affected the Manila clam industry (Choi et al., 2002, Pretto et al., 2014). Moreover, ocean acidification, increased water temperatures, and hypoxic events are increasing worldwide problems (Steckbauer et al., 2015). The molecular mechanisms underlying the tolerance of benthic organisms for diverse biotic and abiotic factors, such as parasites, ocean acidification, and hypoxia, have received much attention in recent years (Lim et al., 2006, Shimokawa et al., 2010, Zhao et al., 2018, Kodama et al., 2018).

Bivalve mollusks have a biphasic life cycle, in which planktonic larvae (trochophores and veligers) settle (as pediveliger larvae) and undergo metamorphosis (to the juvenile stage). The transition of benthic bivalves from a pelagic to a benthic lifestyle during their early development, their burrowing to avoid predation, and their tolerance of hypoxic environments have been attractive characteristics for research (Tamai, 1993, Kodama et al., 2018). However, knowledge of the underlying molecular mechanisms that regulate these processes is still very limited (Findlay and Battin, 2016). The Manila clam is an excellent marine model system, with a robust and multifaceted immune system (Dyachuk, 2016). The immune system of bivalves relies on the recognition of conserved pathogen-associated molecular patterns shared by broad classes of microorganisms, which allows the bivalves to successfully defend themselves against infection (Gestal et al., 2008, Ki et al., 2006, Song et al., 2010). Sediment is a specific habitat of bivalves, with characteristic microbial diversity, low input of sunlight, and a low concentration of oxygen (Wang et al., 2012). Benthic bivalves live in sediment and play important roles in natural biochemical cycles and in the material exchanges between water bodies and sediments (Vaughn and Hakenkamp, 2001). How benthic bivalves adapt to and survive in such hypoxic and pathogen-rich environments is of fundamental interest.

Despite its burrowing lifestyle, R. philippinarum displays a large range of shell colors and patterns (Protas and Patel, 2008, Zhang and Yan, 2010). The Manila clam has several different shell colors, including white, zebra, and orange in its natural habitat (Figure S1), and it has been selected over several generations for high stress resistance and distinctive shell color with selective breeding programs (Zhang and Yan, 2010). Previous experimental crosses have shown that shell color has a genetic basis in the Manila clam (Peignon et al., 1995). Color is one of the most conspicuous phenotypic traits in R. philippinarum, and the diverse colors of molluscan shells are generally believed to be determined by the presence of biological pigments (Mann and Jackson, 2014, Lemer et al., 2015, Feng et al., 2018). Several studies have been conducted of the physiological, biochemical, and molecular responses to different types of biotic and abiotic stress in different shell color strains of the Manila clam (Zhang and Yan, 2010, Nie et al., 2017a, Nie et al., 2017b). In recent years, RNA sequencing (RNA-seq) has been used to identify the genes differentially expressed in the differently colored shells of several bivalves (Sun et al., 2015, Feng et al., 2015, Yue et al., 2015). However, the genes potentially responsible for shell color formation and variation in strains of R. philippinarum are still largely unknown. The expression of genes for tyrosinase (Tyr) and melanogenesis-associated transcription factor (Mitf) is responsible for the production of melanin from tyrosine, promoting melanin production and stimulating pigmentation (Krumholz et al., 2011). In mollusks, tyrosinase is secreted and appears to contribute to shell pigmentation and melanin synthesis, which can also physically encapsulate pathogens, and is therefore, an important component of the immune system (Cerenius et al., 2008).

The main goals of this study were to provide a robust reference genome for the Manila clam and to analyze its transcriptomic responses to various stimuli and different ontogenetic stages. Here, we report the whole genome sequence of the Manila clam and its de novo assembly and a comparative analysis with a Korean strain of the Manila clam (Mun et al., 2017). We also examined the expression of pigmentation-related genes (Tyr and Mitf), immune pathways, and defense mechanisms in the Manila clam. To better understand these biological issues, RNA-seq analyses of hypoxia-challenged, acidification-exposed, and parasite-infected clams were performed, as were RNA-seq analyses of different shell color strains and developmental stages of the Manila clam. This study provides clues to the molecular mechanisms underlying the shell pigmentation, immune defenses, and resistance to different stresses in this clam. Comparative genomic analyses of gene expansion, contraction, and positive selection in R. philippinarum and Lottia gigantea (Simakov et al., 2013), Crassostrea gigas (Zhang et al., 2012), Pinctada fucata (Takeuchi et al., 2012), and Patinopecten yessoensis (Wang et al., 2017) also identify important genetic differences among these species and provide insights into the molecular basis of adaptation to a benthic lifestyle in clams.

Results

Genome Characterization

The genome of R. philippinarum is shown in Figure 1A. The Manila clam genome has a high level of heterozygosity and an estimated size of 1.32 Gb (Table S1 and Figure S2). To reduce the heterozygosity, genomic DNA isolated from an inbred individual (six generations, 2.0% estimated heterozygosity) (Figure S1D) was used for whole-genome shotgun sequencing (286.1 Gb of high-quality paired-end data and 135.1 Gb of mate-pair data; Table S2 and Figure S3). The R. philippinarum assembly comprised 30,811 scaffolds with an N50 length of 345 kb and a total length of 1.12 Gb (Table S3). The heterozygosity of the assembly was reduced to 1.038%, compared with that of a wild-type individual (1.69%) (Table S4). Evaluation of the genome showed the completeness of the assembly (Tables S5–S7). A total of 233 CEGMA (Core Eukaryotic Genes Mapping Approach) core genes with a ratio of completeness of 93.95% (Table S8), together with 91.0% complete and 3.9% fragmented Metazoa BUSCO (Benchmarking Universal Single-Copy Orthologs) orthologs (Table S9), were identified in the assembled genome, indicating the high degree of completeness of the gene regions. Repetitive regions comprised ~38.30% of the R. philippinarum assembly (Table S10, and Figure S4), and combined transposable elements were predicted to constitute 35.8% of the R. philippinarum genome (Table S11). Transcriptome-based, ab initio, and homology-based predictions identified 27,652 protein-coding genes (Table S12), a similar number to other Lophotrochozoa genomes (Table S13 and Figure S5). Of these, 24,570 (88.85%) were annotated based on various public databases (Table S14), with the best BLASTp hits in the SwissProt, TrEMBL, and NCBI nonredundant (NR) protein databases. Gene domain annotations were made by searching the InterPro database. All genes were functionally annotated by their alignment against Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) terms. Noncoding RNA was also annotated (Table S15).

Figure 1

Genome Landscape and Phylogenetic Analysis of R. philippinarum

(A) From outer to inner circles: a marker distribution on 19 chromosomes at the Mb scale; b gene distribution on each chromosome; blue lines indicate genes on the forward direction strand, and yellow lines indicate genes on the reverse direction strand; c GC content within a 1-Mb sliding window; d repeat content within a 1-Mb sliding window; e SNP density of the inbreed offspring; f SNP density of the wild sample.

(B) Phylogenetic tree among 14 species. The split of R. philippinarum was estimated at 476.9 million years ago.

Comparative Analysis of R. philippinarum (Rph) and Korean Clam (Hrph) Genomes

We compared the genomic statistics for Rph and Hrph (Table 1). The contig N50 for Rph was 28,111 bp, whereas that for Hrph was 6,520 bp. The complete BUSCO (C) ratio of the Rph genome was 92.2%, which was much higher than that of the Hrph genome (69.5%). The number of protein-coding genes in Rph was 27,652, whereas that in Hrph was 108,034, and the complete BUSCO (C) ratios of the genes were 91.0% and 81.5%, respectively (Table 1). There was about 67% collinearity between the Rph and Hrph genomes. To analyze the gene families, the genes of Hrph were filtered from 108,034 to 19,776. A gene family cluster analysis was performed between Rph, Hrph, and six other species, including C. gigas, L. gigantea, P. fucata, Patinopecten yessoensis, Octopus bimaculoides, and Lingula anatina, and 436 orthologous gene families were shared among these species. Rph shared 1,266 orthologous gene families with the other species, excluding Hrph. Of these 1,266 genes, 1,240 were not annotated and 26 genes were annotated in more than one copy in the Hrph genome, and 1,246 were annotated in the collinear regions of the Rph and Hrph genomes. Hrph shared 56 orthologous gene families with other species, excluding Rph. Of the 56 genes, 16 genes were not annotated and 40 genes were annotated in more than one copy in the Rph genome, and 48 genes were annotated in the collinear regions of the Rph and Hrph genomes. Compared with Hrph, 7 expanded and 78 contracted gene families were detected in Rph with Fisher's exact test (Table S16). A GO enrichment analysis of the expanded gene families in Rph, relative to Hrph, revealed that active transmembrane transporter activity, neurotransmitter transport, hydrolase activity, and scavenger receptor activity are involved in their immune functions and ecological adaptation (Table S16).

Table 1

Comparative Analysis between the Current Genome of Manila Clam and the Genome of Korean Manila Clam

Type	Rph	Hrph
Estimated genome size (Gb)	1.32	1.37
Assembly total length (bp)	1,122,973,377	1,078,771,101
Sequencing depth	320.63	74.2
Contig N50 (bp)	28,111	6,520
Scaffold N50 (bp)	56,467,786	119,518
Average contig length (bp)	12,041	3,186
Max contig length (bp)	249,659	80,451
Average scaffold length (bp)	56,456,578	80,439
Max scaffold length (bp)	204,629,219	1,050,406
Number of contig > 2 Kb	61,395	121,896
Number of scaffolds > 2 Kb	19	13,318
Genome BUSCO assessment	C:92.2%[S:90.3%,D:1.9%],F:1.6%,M:6.2%,n:978	C:69.5%[S:66.6%,D:2.9%],F:11.2%,M:19.3%,n:978
Repeat percentage	38.29	26.38
Gene models number	27,652	108,034
Average gene length (bp)	12,875	5,117
Average number of exons per gene	7.30	4.17
Average exon length (bp)	200.30	232
Average number of introns per gene	6.29	3.17
Average intron length (bp)	1,697.30	1,230
Gene BUSCO assessment	C:91.0%[S:89.3%,D:1.7%],F:3.9%,M:5.1%,n:978	C:81.5%[S:58.9%,D:22.6%],F:13.2%,M:5.3%,n:978

Rph, the current genome of Manila clam; Hrph, the genome of Korean Manila clam.

Comparative Genome Analysis of the Gene Families

The results of a gene family cluster analysis of R. philippinarum and other 13 species are shown in Figure S6. Compared with four other Mollusca species (C. gigas, L. gigantea, P. fucata, and Patinopecten yessoensis), R. philippinarum has 1,582 unique gene families and shares 6,975 gene families with other species (Figure S7, Tables S17 and S18). A phylogenetic tree was constructed using 288 shared single-copy orthologues of 14 species (Figure 1B). As the closest bivalve clade to the class Gastropoda (L. gigantean), R. philippinarum was estimated to have split from it about 476.9 million years ago (Figure 1B). Compared with the four other Mollusca species, 2,197 expanded and 3,943 contracted gene families were identified in R. philippinarum with Fisher's exact test (Tables S19 and S20–S22), and 29 positively selected genes were detected in the R. philippinarum genome as the foreground branches and genes of four other Mollusca species as background branches (Tables S23 and S24). GO and KEGG enrichment analyses of the expanded genes revealed that they are involved a number of immune-related pathways in R. philippinarum, including in the defense response to scavenger receptor activity, cytokine activity, apoptotic process, G-protein-coupled receptor activity, sulfotransferase activity, tumor necrosis factor (TNF) receptor binding, neuropeptide signaling pathway (Table S25), complement and coagulation cascades, phagosome, Toll-like receptor signaling pathway, NF-κB signaling pathway, TNF signaling pathway, and PI3K–AKT signaling pathway (Table S26).

Differential Expression of Tyr and Mitf Genes in Manila Clam

Manila clams are capable of producing an extraordinary diversity of shell color patterns. We examined the molecular processes of the melanogenesis and tyrosine metabolism pathways, which are involved in pigment synthesis. Tyrosinase (encoded by Tyr) is a crucial enzyme in the process of melanin synthesis, and the Tyr genes are regulated by a melanogenesis-associated transcription factor (encoded by Mitf). In this study, transcriptomic and gene expression analyses of the Tyr and Mitf genes in four shell color strains of R. philippinarum were conducted with RNA-seq and quantitative PCR (qPCR) (Figure 2). In total, 21 Tyr genes and 12 homologues of the Mitf gene were identified in the R. philippinarum genome. Most of the Mitf mRNAs were upregulated in the orange strain (O), and three of them were uniquely and strongly expressed in O compared with the other three color strains (Figure 2A). To confirm the differential expression of Tyr and Mitf in the different shell color strains, six Tyr and seven Mitf genes were further evaluated with RT-qPCR in four different shell color strains of R. philippinarum, and all of them were differentially expressed in the four shell color strains (Figure 2). All seven Mitf genes were upregulated in the mantle of the O strain, compared with the zebra strain, and three Mitf genes (Mitf1, Mitf2, and Mitf5) were upregulated in the mantle of the white zebra (WZ) strain (Figure 2A). An RNA-seq analysis was conducted in the different developmental stages of R. philippinarum, including the trochophore, D-shaped larval, umbo-veliger larval, spat, early juvenile before shell pigmentation, and late juvenile after shell pigmentation stages. The expression of four Mitf genes was upregulated between the early juvenile stage before shell pigmentation and the juvenile stage after shell pigmentation (Figure 2A). Interestingly, the Tyr gene family showed a distinctive expression pattern and stage-specific expression in pelagic larvae (trochophore, D-shaped larva, and umbo-veliger larva), benthic juveniles (early juvenile and late juvenile), and adult clams (Figure 2B).

Figure 2

The Expression of Key Genes Impacting Shell Color Patterns in R. philippinarum and Evolutionary Analysis of Tyr Genes

(A) RNA-seq hierarchical clustering of Mitf genes at four different developmental stages, at noncolored and colored developmental stages in two different shell color strains, and in four different shell color strains, as well as RT-qPCR expression of seven Mitf genes in four different shell color strains.

(B) RNA-seq hierarchical clustering of Tyr genes at four different developmental stages, at noncolored and colored developmental stages in two different shell color strains, and in four different shell color strains, as well as RT-qPCR expression of six Tyr genes in four different shell color strains. The abbreviations include trochophora (A11), D-shaped larvae (A12), umbo larvae (A23), pediveliger larvae (z3), uncolored zebra clam (Za), colored zebra clam (Zb), uncolored orange clam (Oa), colored orange clam (Ob), Orange strain (O),White strain (W), White Zebra strain (WZ), and Zebra strain (Z).

Expansion of Immune-Related Genes in R. philippinarum

Compared with four other Mollusca species (C. gigas, P. fucata, L. gigantea, and Patinopecten yessoensis), a number of immune-related genes are expanded in R. philippinarum, including the thioester-containing proteins (TEP), complement component C3 (C3), C1q domain-containing proteins (C1qDC), lysozymes, sialic acid-binding lectin (SABL), and heat shock protein 70 (Hsp70) (Figure 3). These gene families expanded in R. philippinarum are enriched in complement and coagulation cascades, phagosome, and lytic pathways (Figure 3). Pathogen recognition in R. philippinarum can be attributed to pathogen recognition receptors C1qDCs, Bf, and carbohydrate-binding lectins, such as SABLs and mannan-binding lectin, through the complement pathways (Figure 3A).

Figure 3

Analysis of Immunity in R. philippinarum

(A) A schematic diagram depicting the complement system in R. philippinarum. Hexagons with bold black borders indicate gene families (HSP70s, C3, SABLs, and lysozymes) expanded in R. philippinarum. The pathogen recognition in R. philippinarum can be equipped with pathogen recognition receptors C1qDCs, Bf, and carbohydrate-binding lectins such as SABLs and MBL through three well-defined activation pathways—the classical, alternative, and lectin pathways. Then, the central complement component C3 was proteolytically activated by converting into fragments C3a and C3b. Upon proteolytic activation, C3b bonds covalently to the surface molecules of microbes using its intrachain thioester bond and initiates the formation of the membrane-attack complex (MAC), which can eliminate the invading pathogens by the lytic pathway. On the other hand, C3a and degraded C3b can be recognized by complement C3a receptor 1 (C3AR1) and complement component (3b/4b) receptor 1 (CR1) on phagocytes, respectively, which may induce the chemotaxis and opsonization of target cells. C1qDCs, C1q domain-containing proteins; SABLs, sialic acid-binding lectins; MBL, mannose-binding lectin; MASPs, Mannan-binding protein associated serine proteases; Bf, Factor B; C3, complement component 3; HSP70s, heat shock 70 kDa proteins; C8, complement component 8; SRs, scavenger receptor.

(B) Phylogenetic tree of TEP superfamily from selected organisms. The tree shows the three classic subgroups of the TEP superfamily, the C3 subgroup in purple, the alpha 2-macroglobulin protease inhibitors (A2M) in red, and the CD109/iTEPs subgroup in blue. The pink blocks, green blocks, yellow blocks, orange blocks and gray blocks indicate the genes of R. philippinarum, P. yessoensis, P. fucata, C. gigas, and L. gigantea. R. philippinarum has more genes of C3 and CD109/iTEPs subgroups (seven and seventeen, respectively) in comparison with C. gigas, L. gigantea, P. fucata, and P. yessoensis, which suggested a species-specific gene expansion in R. philippinarum. The abbreviations of species are: ACA, Aplysia californica; ADI, Acropora digitifera; AGA, Anopheles gambiae; AJA, Apostichopus japonicas; BBE, Branchiostoma belcheri; BFL, Branchiostoma floridae; CFA, Chlamys farreri; CGI, Crassostrea gigas; DLI, Diadumene lineata; DME, Drosophila melanogaster; DRE, Danio rerio; EBU, Eptatretus burger; ESC, Euprymna scolopes; ESI, Eriocheir sinensis; ETA, Euphaedusa tau; FCH, Fenneropenaeus chinensis; HAD, Hasarius adansoni; HAS, Homo sapiens; HCU, Hyriopsis cumingii; ISC, Ixodes scapularis; ICA, Lethenteron camtschaticum; LGI, Lottia gigantea; LJA, Lethenteron japonicum; LSP, Limulus sp.; PFU, Pinctada fucata; rde, Ruditapes decussatus; RPH, Ruditapes philippinarum; SEX, Swiftia exserta; SPU, Strongylocentrotus purpuratus; TTR, Tachypleus tridentatus.

(C) Phylogenetic tree of HSP70 family genes in R. philippinarum, C. gigas, L. gigantea, P. fucata, and P. yessoensis. The pink blocks, green blocks, yellow blocks, orange blocks, and gray blocks indicate the genes of R. philippinarum, P. yessoensis, P. fucata, C. gigas, and L. gigantea. All bivalve HSPA12 genes from R. philippinarum, C. gigas, P. fucata, and P. yessoensis are clustered with red blocks. Gene names in purple belong to R. philippinarum. The other clusters consisting of HSPA1, HSPA4, HSPA5, HSPA8, HSPA9, HSPA13, HSPA14, HSPH1, HSP70B2, and HSP68 are clustered with blue blocks.

Ruditapes philippinarum has the largest number of genes in the C3 and CD109/iTEP subgroups (7 and 17 members, respectively), compared with C. gigas, L. gigantea, P. fucata, and Patinopecten yessoensis (Table S27), suggesting species-specific gene family expansion in R. philippinarum. Phylogenetic analyses also indicated extensive expansion of the Hsp70 genes from the Hspa12 subfamily in R. philippinarum and also in the bivalve lineages (Figure 3C). The expansion of immune-related genes (SABL, C3, C1qDC, lysozyme, and Hsp70) provides evidence of an association between gene duplication and the organism's adaptation to pathogen-rich environments (Figure 4). The characteristics and expression of five complement C3 genes (xfsc0000076.27, xfsc0000154.2, xfsc0000773.3, xfsc0000773.2, and xfsc0001170.6) were investigated in the blood and digestive gland of R. philippinarum at 3, 6, 12, 24, 48, 72, and 96 h after treatment with lipopolysaccharide (LPS), peptidoglycan (PGN), or polyinosinic–polycytidylic acid (poly[I:C]). All these complement C3 genes were significantly upregulated in the immune responses of R. philippinarum to different pathogenic stimuli (LPS, PGN, and poly[I:C]) (Figure S8), indicating that the complement C3 genes play critical roles in the immunity of the Manila clam.

Figure 4

The Heatmap of the Immune-Related Genes Number in L. gigantea, R. philippinarum, P. yessoensis, C. gigas, and P. fucata

RNA-Seq Analyses at Different Ontogenic Stages

To identify the key genes involved in the development and metamorphosis processes, RNA-seq analyses were conducted in the different ontogenic stages of R. philippinarum. The Hox gene clusters and their expression at different ontogenic stages are shown in Figure 5. The expression of 11 Hox genes was detected in 17 different stages, and the Post1 and Lox5 genes were expressed in the gastrula stage (Figure 5B). Most Hox genes were strongly expressed in the early trochophore stage, when the larvae begin to secrete the shell to coat and protect the embryo. The Hox1, Post2, Hox3, and Lox4 genes were mainly expressed in trochophore and veliger larvae (Figure S9). In contrast, the Hox2 and Hox4 genes were strongly expressed in the umbo larvae (Figure 5B). The Post2 and Lox4 genes were strongly expressed in the trochophore, veliger, and umbo larvae. The Antp and Lox2 genes were strongly expressed in the umbo larvae, spats, and juveniles, whereas Lox5 was strongly expressed in the gastrula stage (Figure 5B). The expression of immune-related genes TEP and C3 was negligible throughout the different planktonic developmental stages (trochophore, D-shaped larva, and umbo-veliger larva) of R. philippinarum, but was induced during metamorphosis and increased significantly during the transition from the umbo-veliger larva to the spat stage (Figure S10). However, eight Hsp70 genes were strongly expressed in the planktonic larval stages (trochophore, D-shaped larva, and umbo-veliger larva) compared with their expression in the benthic spat stage, whereas 40 Hsp70 genes were strongly expressed in the spat stage during metamorphosis (Figure S11).

Figure 5

Hox Genes Clusters and the Expressions at Different Developmental Stages in R. philippinarum

(A) Clustering of Hox genes in R. philippinarum (Rph), C. gigas (Cgi), L. gigantea (Lgi), P. fucata (Pfu), and P. yessoensis (Pye) genomes. The relative position and orientation of the genes are indicated.

(B) RT-qPCR hierarchical clustering of Hox genes at seventeen different developmental stages in R. philippinarum. Egg, eggs; Zyg, zygotes; Pol1, first polar body; Pol2, second polar body; C2, 2 cell; C4, 4 cell; C8, 8 cell; C16, 16 cell; C32, 32 cell; Gas, gastrula; ET, early trochophora; LT, late trochophora; Vel, veliger larvae; UL, umbo larvae; Ped, pediveliger larvae; Juv, juvenile; HC, have color juvenile.

Transcriptomic Responses to Hypoxia, Acidification, and Parasites

We hypothesized that the expansion of the complement and immune-related genes plays a critical role in the resistance of the Manila clam to different stresses and in its benthic adaptation. The transcriptomic responses of the Manila clam to hypoxia, acidification, and parasites were analyzed to identify the complement- and immune-responsive genes. An RNA-seq analysis identified 102 differentially expressed genes (DEGs), including 75 upregulated genes and 27 downregulated genes, corresponding to 12 major physiological functions in the gills of R. philippinarum after exposure to hypoxia for 2 days (Table S28). The expression of hypoxia-responsive genes encoding Hsp70, C-type lectin, complement factor B-1 (CFB), Ikappa, melatonin receptor type 1A (MRTA), serine/threonine-protein phosphatase (STP), growth arrest, DNA damage-inducible protein (GADD), phosphoenolpyruvate carboxykinase (PECKG), and Tyr (xfSc0000104.14) was upregulated (Figure S12), whereas the expression of genes encoding TNF, NFX1, GTPase IMAP (IMAP), E3, and hemicentin was downregulated when analyzed with RT-qPCR after exposure to hypoxia for 0, 2, 5, or 8 days (Figure S12). These genes are implicated in different physiological pathways, including respiration and energy metabolism (carbonic anhydrase, SLC2A1) and antioxidant and immune defense (CFB, C-type lectins, Hsp70, lysozyme, Tyr, TNF, ANKK1, IκB, GPx).

To confirm the DEGs identified in clams under acidic conditions, 11 genes (six upregulated and five downregulated) were selected from the DEGs in several canonical pathways of interest and analyzed with RT-qPCR (Figure S13). A number of vital complement- and immune-responsive genes were significantly differentially expressed, including complement C1q, SABL, C-type lectin, Hsp70, and CYP450. An RNA-seq analysis of clams infected with digenetic trematodes identified 460 DEGs, including 416 upregulated and 44 downregulated genes (Table S29). Eleven candidate genes associated with immune, antioxidant, and molecular chaperone functions were chosen from the DEGs to confirm their expression levels with qPCR (Figure S14). In the infected clams, the expression of complement C1q and C-type lectin was significantly upregulated in the mantle and gonad. Expression of the CYP450 gene was significantly downregulated in the gills but significantly upregulated in the mantle and gonad. The expression of inhibitor of apoptosis (IAP) and glutathione S-transferase (GST) was significantly elevated in the three tissues (Figure S14). The expression of Hsp70 was significantly elevated in the gills of R. philippinarum after 2 and 5 days of hypoxia and after 12 h of acidification stress (both p < 0.05). A transcriptomic analysis of clams exposed to the parasite revealed that the expression of Hsp70 was significantly increased (p < 0.05) in three tissues tested, including the gill, mantle, and gonad (Figure S14). Overall, the RNA-seq analysis revealed that the expanded Hsp70 genes are involved in multistress resistance in R. philippinarum.

Expansion of the G-Protein-Coupled Receptor-Related Genes

G-protein-coupled receptors (GPCRs) were enriched in the unique and expanded genes of R. philippinarum. The numbers of genes containing the seven-transmembrane-associated domain in R. philippinarum, Patinopecten yessoensis, P. fucata, C. gigas, and L. gigantea were 858, 463, 533, 606, and 385, respectively (Table S30). Compared with the four other Mollusca species (C. gigas, P. fucata, L. gigantea, and Patinopecten yessoensis), the number of GPCR-related genes was expanded in R. philippinarum, including the chemokine receptor, D dopamine receptor, FMRFamide receptor, G-protein-coupled receptor GRL101, neuropeptide receptor, probable G-protein-coupled receptor 139, probable G-protein-coupled receptor Mth, thyrotropin-releasing hormone receptor, and cysteinyl leukotriene receptor (Table S30, Figure S15). These expanded GPCR gene families in R. philippinarum may participate in the regulation of various physiological functions, especially immune functions and the inflammatory response, which are potentially important for stress resistance and the clam's adaptation to benthic life.

Discussion

Notably, the contig N50 of R. philippinarum was 28.1 kb, which is similar to those of other species in the bivalve lineage. Moreover, the complete BUSCOs in the genome of R. philippinarum was 92.2%, which was higher than that in Modiolus philippinarum, Limnoperna fortunei, Argopecten purpuratus, C. gigas, O. bimaculoides, and Saccostrea glomerata (Table S31 and Figure S16). In summary, the genome assembly index of R. philippinarum was similar to that of other species of Mollusca.

Ruditapes philippinarum has a variety of colors and a distinct set of genes that are either expanded or highly expressed and facilitate its adaptation to extreme environments. Many bivalves express multiple tyrosinases, and the expansion of the tyrosinase genes appears to be a common feature of bivalves, with more than 20 gene family members present in R. philippinarum, C. gigas, P. fucata, and Patinopecten yessoensis (Table S32). In the present study, qPCR analyses of the Mitf and Tyr genes showed that the Tyr genes have a distinctive expression pattern in different developmental stages and that Mitf gene expression increases significantly from the early development stages that lack pigmentation to the juvenile stages with colored shells. This indicates the potential roles of these genes in shell biomineralization and shell color pigmentation. Three Mitf genes are uniquely expressed in the O strain, suggesting that these Mitf proteins are involved in the pigmentation responsible for its orange color. The distinctive expression patterns of the Mitf and Tyr genes that were demonstrated with RNA-seq and qPCR analyses shed light on the molecular basis of shell pigmentation and the color diversity in the Manila clam.

Here, we have provided an overview of the immune-related genes of the Manila clam that were identified with genomic, transcriptomic, and comparative genomic analyses. The ability of the immune system to detect pathogen-derived molecules is mediated by pattern recognition receptors (PRRs) (Janeway and Medzhitov, 2002). Several groups of invertebrate PRRs have been characterized in the Manila clam, such as TEPs, Toll-like receptors, peptidoglycan recognition proteins, fibrinogen-related proteins, galectins, C-type lectins, and complement-related proteins. In this study, the TEP superfamily genes were shown to be expanded in the R. philippinarum genome. TEPs bind to the surfaces of exogenous microorganisms with conserved thioester bonds, triggering a series of immune responses, including the activation of the complement system and phagocytosis (Fujita, 2002).

Complement component 3 (C3) is central to the complement system, and plays an important role in immune defense, immune regulation, and immune pathology. C3 is the core member of the complement system because it is the intersection of several activation pathways, and depends on a positive feedback loop from C3b. Upon its proteolytic activation, C3b binds covalently to the surface molecules of microbes via intrachain thioester bonds and initiates the formation of the membrane-attack complex (MAC), which eliminates invading pathogens via the lytic pathway (Anderluh and Gilbert, 2014). A C3a-like peptide has been shown to mediate the chemotaxis of hemocytes in Ciona intestinalis (Pinto et al., 2003) and Pyura stolonifera (Raftos et al., 2003). The activated fragments of C3 (C3a and C3b) are also involved in the activation of macrophages and mastocytes (Gao et al., 2013).

C3 is a gene conserved throughout evolution, from cnidarians to humans (Sahu and Lambris, 2001). The number of C3 gene isoforms in highly evolved vertebrates is small, but there is usually a large number of C3 genes in primitive vertebrates and invertebrates, although the number varies between species (Huan et al., 2012). C3 is an essential molecule in the complement system, acting as an opsonin in the phagocytosis of microorganisms. Therefore, as an important part of the immune system, C3 is particularly crucial to shellfish health because most shellfish lack an adaptive immune system. The number of genes in the TEP superfamily (C3, CD109, and A2M) is higher in the Manila clam than in other bivalve species (Table S27), suggesting that TEP expansion is related to the immune defense system and benthic adaptation of the Manila clam. The complement pathway genes (SABL, C1qDCs, C3, C3a, C3b) are also expanded in the Manila clam, which is evidence of an association between gene duplication and pathogen-rich ecological conditions. The pathogen recognition mechanisms of R. philippinarum include the pathogen-recognition receptors C1qDCs, Bf, and carbohydrate-binding lectins. The expansion of the C1q genes in the Manila clam has also been reported by Mun et al. (2017). Zhang et al. (2015) demonstrated that the C1qDC genes are expanded in the Pacific oyster and that C1qDC expression is upregulated by both biotic and abiotic challenges.

Transcriptomic analyses and qPCR confirmation showed that Hsp70 expression increased significantly in R. philippinarum after hypoxia stress for 2 or 5 days. The expression of Hsp70 also increased significantly after exposure to a parasite or acidification. Therefore, the Hsp70 genes are considered to be involved in the multistress resistance of R. philippinarum. Hierarchical clustering showed that the Hsp70 genes have different expression patterns at different developmental stages (Figure S11). The first cluster of eight Hsp70 mRNAs is only expressed in the pelagic larval stages, whereas other Hsp70 genes are mainly expressed in the pediveliger larval stages, implying that Hsp70 has potential roles in the immune responses in the different larval stages of the Manila clam. Overall, the expansion of the immune-related gene families (HSP70, C3, SABL, and lysozyme) is vital to the resistance to extreme conditions and the benthic adaptation of R. philippinarum.

The Hox gene family directs morphological development and tissue differentiation along all the principal axes of an embryo. Many Hox genes are expressed during early embryogenesis, suggesting their roles in development (Krumlauf et al., 1987, Utset et al., 1987, Dolecki et al., 1988). The overall form assumed during embryonic development is controlled by the actions of the Hox proteins, their cofactors, and ultimately by the DNA regulatory elements that control their expression (Lufkin, 2005, Hejnol and Martindale, 2009, Hinman et al., 2003). In Gibbula varia, the strong expression of Gva-Hox2 was observed in trochophore and D-shaped larvae, whereas Gva-Post2 was expressed in cells around the shell field in the early trochophore, which form the future mantle edge (Samadi and Steiner, 2009). In the present study, Post2 was strongly expressed in the trochophore, veliger, and umbo larvae, demonstrating its significant role in larval development. Lox5 was mainly expressed in the umbo larva, suggesting that Lox5 is involved in regulating morphogenesis. In G. varia, the expression of Gva-Lox5 is lost in the velum of the posttorsional veliger, although the gene is still expressed in the cerebral ganglia and their commissures (Samadi and Steiner, 2009). High levels of Hox1 and Hox4 expression were detected in the D-shaped larvae of both the oyster and scallop (Zhang et al., 2012, Wang et al., 2017). The Hox proteins play pivotal roles in molluscan larval development and shell formation (Lufkin, 2005).

In conclusion, the expansion of the immune-related genes and complement pathways plays an important role in the benthic adaptation of the Manila clam. The ecology of a benthic lifestyle implies that these animals interact with high concentrations of microorganisms at the low oxygen levels in sediments. However, the Manila clam has well-adapted defense strategies to protect itself from pathogens and is recognized as a robust species that is extremely resistant to adverse environmental conditions.

Limitations of the Study

Here, we have provided a comprehensive framework for understanding the genetic adaptations of R. philippinarum to the benthic environment. Although the quality and contiguity of the R. philippinarum genome assembly was carefully validated and was generally reliable for the current study, the contig/scaffold N50 did not reach the optimal level for molluscan genome assemblies. With the development of the high-throughput sequencing technology, the quality of the assembly will be improved in future studies. The asymmetry in the quality of the two clam genomes compared means that the findings from this comparison must be considered with caution. Functional experimental assays are also required to confirm the expansion of these gene families and to identify the targets involved in the adaptation of R. philippinarum to its hypoxic, pathogen-rich environment.

Methods

All methods can be found in the accompanying Transparent Methods supplemental file.