{'fraction_non_alphanumeric': 0.040144459228283594, 'fraction_numerical': 0.024310967496673636, 'mean_word_length': 5.069566220581449, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Clam Genome Sequence Clarifies the Molecular Basis of Its Benthic Adaptation and Extraordinary Shell Color Diversity

Summary 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., 2017aNie 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, acidificationexposed, 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.

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). 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-kB 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).

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  (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).
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.

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 . 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. 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).

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, IkB, 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 RNAseq 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  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.

DECLARATION OF INTERESTS
The authors declare no competing interests.

Lorem ipsum
Missing Complete

Sampling
The fifth generation of orange shell color strains of R. philippinarum was collected from the Zhangzidao aquaculture farm in Dalian, Liaoning Province, China. The clams with a similar shell length averaging at 20.90±1.34 mm, respectively, were used in this study. After being transported from the field to the laboratory, the clams were acclimated in aerated 70 L plastic tanks, containing water at 20 °C with salinity of 30. All the samples were acclimated in the laboratory for one week before the sampling. The Manila clam is not an endangered or protected species, so no specific permits were required for the study. One of the clams was sacrificed for DNA extraction.

The present study was approved by College of Fisheries and Life Science, Dalian Ocean
University. Two adult individuals of orange strain of Manila clam Ruditapes philippinarum were selected for whole genome sequencing and genomic DNA was extracted from the whole individuals for short and long libraries of the whole-genome shotgun sequencing. One of the F6 Orange strains of Manila clam with distinctive orange shell color and highly inbred individual was selected for whole genome sequence (Fig. S1D).

Genome sequencing
Genomic DNA was broken into random fragments according to the whole genome shotgun Gb raw data produced by MP libraries.
To reduce the effect of sequencing error to the assembly, the following type of reads would be filtered: reads with ≥10% unidentified nucleotides for libraries, reads from short-insert libraries having more than 50% bases with quality score less than 5, reads having the sequencing adapters, read1 and read2 of two paired-end reads that were completely identical (considered to be products of PCR duplication). After filtering out the adapter sequences, low quality and duplicated reads, a total of 361.4 Gb (~274-fold coverage) data were remained as clean data for the whole-genome assembly.
To characterize the distribution of single-nucleotide polymorphism (SNP) and estimate the level of heterozygosity in the R. philippinarum genome, short pair-ended reads from the inbred and wild samples were aligned to the assembled genome. We identified 33,830,021 and 11,758,346 heterozygous SNPs in the assembled (inbred) and resequenced (wild) genome with the polymorphism rate of 1.038% and 1.690%, respectively (Table S4). As expected, inbreeding has significantly reduced its polymorphisms. The level of polymorphism in R. philippinarum is similar to that in Lingula genome (1.6%), but higher than that found in the scallop Patinopecten yessoensis (1.04%) and the Pacific oyster Crassostrea gigas (1.30%), indicating that the R.

Estimation of genome size
The distribution of K-mer frequency is widely used for the estimation of genome size. Total of 58 Gb high-quality short-insert reads (230 bp) were used to estimate the Manila clam genome size by analyzing the K-mer frequency (defined as 17 bp here), the coverage sequencing depth of which follows Poisson distribution . The formula 'Genome size = total_kmer_num/kmer_depth' (Liu et al., 2013), where total_kmer_num is the total number of K-mer and kmer_depth is the peak sequencing depth on the K-mer frequency distribution map, was used to estimate the genome size (Fig. S2). And the genome size of R. philippinarum were determined to be 1.32 Gb with high ratio of heterozygosity (Table S1).

Genome assembly
It is a challenging task to obtain a high-quality genome of R. philippinarum due to the high levels of heterozygosity (~2.0%). To tackle this challenge, an innovation assembly software NOVOheter2.0, developed by Novogene Bioinformatics Institute (Beijing, China) on basis of SOAPdenovo package , was applied for the Manila clam genome assembly. Firstly, the heterozygous K-mers were determined from coverage information obtained from de Bruijn graph in R. philippinarum genome. All high-quality reads from short-insert size libraries were used to split and counted to generate 19-mer frequency distribution, and 19-mers with sequencing depth ranged from 23 to 124 were regarded as the initial heterozygous 19-mer set. Next, a de Bruijn graph was constructed by splitting the high-quality reads from short-insert size libraries into 45-mers with SOAPdenovo setting parameters as '-e 1 -M 1 -R', and contigs with sequencing depth ranged from 20 to 90 were regarded as heterozygous contigs (Fig. S3).
Secondly, PE and MP reads were mapped to the heterozygous K-mer set to obtain the coverage relationship, and an effective link between two heterozygous contigs was defined if there were at least three read pairs mapped. Two contigs in a bubble represents two potential haploids in a diploid genome. The heterozygous contigs were clustered together as bubble contigs considering the order and distance between two heterozygous contigs. After clustering and classifying heterozygous contigs, there were two big groups of bubble contigs waiting to be removed in the next step.
The heterozygous contigs were matched back to the corresponding edges in de Bruijn graph constructed by 45-mers, one of the heterozygous contig groups on the edges were removed to solve the bubble structures or intricate structures in the graph caused by heterozygosity of diploid R. philippinarum genome. After merging bubble contigs, the contig sequences were obtained finally.
After obtaining contig sequences, all high-quality short-insert PE reads and large-insert MP reads were realigned onto these contigs by Platanus assembler (v1.2.1) (Kajitani et al., 2014) with scaffold protocol (default parameters). And the shared paired-end relationships between each pair of adjacent contigs were used to construct scaffolds.
Gaps in the constructed scaffolds were filled by GapCloser (English et al., 2012) (http://sourceforge.net/projects/soapdenovo2/files/GapCloser) with parameters setting as '-l 150 -p 30' based on the paired-end information of the read pairs that had one end mapped to the unique contig and the others located in the gap region.
After removing the heterozygosity in R. philippinarum genome by HaploMerger (Huang et al., 2012) with default parameters, a total of 30,811 scaffolds were assembled covering total length over 1,123 Mb. And a scaffold N50 length of 345,005 bp, contig N50 length of 28,111 bp and total length of 1.12 Gb were gained in the final Manila clam genome assembly (Table S3).

Genome evaluation
The high-quality short insert PE reads were aligned onto the assembled genome to evaluate the accuracy by using BWA (Li and Durbin, 2009) (http://bio-bwa.sourceforge.net/) with parameters setting as '-o 1 -i 15'. The distribution of the sequencing depth at each position was calculated to measure the completeness of the genome assembly. The assembly with GC content equal to 31.89% (Table S5) was evaluated through mapping reads with approximately 93.32% of mapping rate and 99.83% of coverage (Table S6), explaining the completeness of the genome assembly.
Raw data of transcriptomes from the adductor muscle (5.84 Gb) were sequenced by Illumina sequencing platform. Transcripts were assembled by Trinity (Grabherr et al., 2011) with parameters setting as '-ss 0.5 -jc 0 -minkmercov 2 -minglue 2' to be as expressed sequence tags (ESTs). The 78,505 random ESTs were aligned to the scaffolds by Blat (Kent, 2002) (http://genome.ucsc.edu/goldenpath/help/blatSpec.html) with default parameters, 81.52% of which (> 200 bp) were more than 90% covered by a scaffold (Table S7). The final assembly version covered almost of the genome and the region of coding genes.
CEGMA (Parra et al., 2007) (Core Eukaryotic Genes Mapping Approach, http://korflab.ucdavis.edu/dataseda/cegma/) defined a set of conserved protein families that occur in a wide range of eukaryotes, and present a mapping procedure that accurately identifies their exon-intron structures in a novel genomic sequence. Through mapping to the 248 ultra conserved eukaryotic genes, a total of 233 genes with a ratio of 93.95% were found (Table S8), supporting a high ratio of core conserved gene families, which indicated that our assembly reached a better integrity of the core genes.

Annotation of repeats
After the genome assembly, repeat annotation was carried out for the R. philippinarum genome.
Repetitive sequences include transposable elements (TEs) and tandem repeats. Two approaches were used to discover TEs. RepeatMasker (Smit, Chen, 2004) (version 3.3.0) found TEs in a integrate repeat library which was derived from known repeat library (Repbase 15.02) and the de novo repeat library, built by RepeatModeler (Smit, 2017) (Vision 1.0.5), RepeatScout (Price et al, 2005), and LTR_FINDER (Xu and Wang, 2007). RepeatProteinMask (Chen, 2004) was performed to detect TEs in R. philippinarum genome by comparing the TE protein database. Tandem repeats were ascertained in the genome using Tandem Repeats Finder (Benson, 1999). Repetitive sequences comprised ~38.30% of the R. philippinarum genome, and the results were shown in detail in Table S10, Table S11 and Fig. S4.
There were two ways to assemble the RNA-seq data into the unique sequences of transcripts. One was mapping the RNA-seq data to the R. philippinarum genome using Tophat (Trapnell et al., 2009) (version 2.0.8) and using cufflinks (Trapnell et al., 2010) (version 2.1.1) (http://cufflinks.cbcb.umd.edu/) for transcript assembly. The other was applying Trinity (Grabherr et al, 2011) to assemble the RNA-seq data, and then PASA (Haas et al., 2003) software (http://pasapipeline.github.io/) improved the gene structures. A w eighted and non-redundant gene set was generated by EVidenceModeler (EVM) (Haas et al., 2008) which merged all genes models predicted by the above three approaches. Combining with transcript assembly, PASA adjusted the gene models generated by EVM. At last, the gene set were filtered according to the following standards: coding region lengths of ≤ 50 amino acids, supported only by de novo methods and with FPKM < 5. The final reference gene set contained 27,652 protein coding genes. (Table S12,   Table S13 and Fig. S5)

Functional annotation
Functional annotation of protein coding genes was obtained according to best BLAST hit by BLASTP (E-value ≤ 1E-05) against SwissProt (Bairoch and Apweiler, 2000), TrEMBL (Bairoch and Apweiler, 2000) and NCBI non-redundant (NR) protein databases. Motifs and domains were annotated by using InterProScan (Mulder and Apweiler, 2007) (Kanehisa and Goto, 2000) pathway to identify the best match classification for each gene.

Non-coding RNA annotation
The tRNA genes were predicted by tRNAscan-SE software (Lowe and Eddy, 1997). The miRNA and snRNA fragments were identified by INFERNAL (Nawrocki et al, 2009) (Table S15)

Gene family cluster
Gene families were generated by OrthoMCL   program, the longest transcript was selected from alternative splicing transcripts belonging to one gene and genes with ≤ 50 amino acids were removed. The alignments with high-scoring segment pairs (Hsps) were conjoined for each gene pair by solar (Yu et al., 2006). To identify homologous gene-pairs, more than 30% coverage of the aligned regions in both homologous genes was required. Finally, the alignments were clustered into gene families using OrthoMCL with 1.5 inflation index. After clustering, 33,076 gene families and 288 single-copy orthologs were detected across R. philippinarum and 13 other species (Fig. S6 and Fig. S7). GO enrichment and KEGG enrichment of unique gene families in R. philippinarum were presented in Table S16 and   Table S17.

Phylogenetic tree construction and divergence time estimation
The above mentioned 288 shared single-copy orthologs were utilized to construct the phylogenetic tree. Protein sequences of these orthologs were aligned by MUSCLE (Edgar, 2004), then protein alignments were transformed to CDS alignments. Using this CDS alignments, the phylogenetic tree was constructed by the ML (maximum likelihood) TREE algorithm in RAxML software (Stamatakis, 2006;Stamatakis et al., 2008)  sample-number=100,000, and sample-frequency=2. Four calibration points were selected from Simakov et al., (2015) and TimeTree (Hedges et al., 2015) website ( were presented in Table S19-S22.

Positively selected genes
The protein alignments of 4,494 single-copy gene families in R.philippinarum, C.gigas, L.gigantea P.fucata and P.yessoensis were generated using MUSCLE. Gblocks (Castresana, 2000) was applied to filter poorly aligned positions and divergent regions of the protein alignments before transformed to CDS alignments. With the R. philippinarum as foreground branch, positive selection sites were detected based on branch-site models (Zhang et al., 2005) of PAML using the CDS alignments. P-values were computed using the χ2 statistic and adjusted by FDR method. In addition, in order to ensure the accuracy of the positive selection results, the genes with low alignment quality were filtered with the following criteria (PMID:25362486): gene length more than 300 bp, and no gap within 3 amino acids around the positive selected site in the five species.
Finally, 29 positively selected genes were found in R. philippinarum genome. GO enrichment and KEGG enrichment of positively selected genes in R. philippinarum were presented in Table S23 and Table S24.

Comparative genome analysis between R. philippinarum and Korean Manila clam genome
The Korean Manila clam genome sequencing data (Mun et al., 2017) previously published was used to conduct the comparative analysis between the two draft Manila clam genome. Since the gene prediction of the genome of Korean clams is the genomic data (V1) without deleting heterozygosity effectively, there is a large redundancy in the resulting of 108,034 gene models.
The software of HaploMerger removes the hybrid segment without interrupting the original genomic version (Huang et al., 2012). Based on which, we mapped the genome data of the last version (V2) of Korean clams to the genome of V1 version using blast. The sequences would be removed except the sequences were matched perfectly without "N". The number of genes finally obtained was 19,776 according to the gene models in the matched sequences.

Linkage map-based chromosome anchoring
A total of 20,581marker sequences with physical distance from 19 linkage group were obtained from a linkage map of R. philippinarum reported by Nie et al. (2017c). Only markers with a unique location were used for anchoring and orienting scaffolds to different linkage groups.
Scaffolds in conflict with the genetic map (such as markers from a different linkage on the same scaffold) were checked manually with 10 kb mate-paired reads. At last, 3,888 markers were used to construct the chromosome.

RNA-seq of different developmental stages of R. philippinarum
Healthy R. philippinarum (shell length: 3.2±0.3 cm) were collected from Zhangzidao (Dalian, Liaoning province, China) and maintained in 70 L tanks containing aerated sand-filtered seawater (salinity: 31±1 ppt, pH 8.1±0.1) at 22±1°C for one week prior to experimentation. Artificial breeding larvae were cultivated and collected at different stages including trochophore (A11), D-shape larva (A12), umbo larva (A23) and pediveliger (Z3). The early phase of juvenile without shell color and late phase of juvenile with shell color were also sampled for RNA-seq in the orange strain and zebra strain, respectively. The collected samples were stored at -80°C before RNA extraction.

Transcriptome response to acidification, parasite and hypoxia stress
After preliminary acclimation, one hundred clams dividing into two groups were cultured on two 70 L tank for 3 months at pH 8.0 (control) and 7.4 (treated). For water acidification, CO 2 gas was diffused directly into seawater contained in each tank to achieve testing pH conditions. Seawater pH was continuously monitored and controlled us ing a pH meter. Total alkalinity (TA), pH, temperature and salinity, were measured daily to determine CO 2 partial pressure ( philippinarum, collected from Zhuanghe, Dalian, China, were divided into two groups of 50 animals. After a 7-day acclimatization period in tanks supplied with aerated filtered seawater, clams were exposed to hypoxia (Dissolved oxygen, DO: 2 mg/L) and normoxia (DO: 8 mg/L) as control groups. At day 0, the start of the experiment, the DO in the water was decreased to 2 mg/L using an oxygen depletion system through a fine bubble aerator where nitrogen was injected.
Oxygen removal was controlled by nitrogen flow. Surface gas exchange in the rearing tank was limited by setting the water inflow under the water surface. The O2 concentration in the tank was monitored using YSI oxymeter and adjusted when necessary to keep hypoxia level constant all along the experiment. Normoxia was obtained by equilibrating seawater with air. The vitality and motility of the clams suffered no adverse affects when DO level changed and no mortality were observed either in the control or in the hypoxia-exposed clams.
Three biological replicates per each time and condition were collected. Total RNA was extracted from the gills of 12 controls and 12 hypoxia exposed clams at 0, 2, 5, 8 days of exposure using RNAprep pure Tissue Kit (TianGene, Beijing, China), according to the manufacturer's protocol.
RNA purity was determined with a NanoPhotometer® spectrophotometer (Implen, CA, USA). Sequencing reads were aligned to the R. philippinarum genome using hisat2 with default parameters. Gene expression levels in terms of FPKM were estimated by custom Perl scripts.
Differentially expressed gene (DEG) analysis was carried out using edgeR with three biological replicates, and genes with a fold change value ≥ 2 and adjusted p-value < 0.05 were defined as significant DEGs.

Identification of the immune-related genes
The immune-related protein sets from NCBI and Uniprot (reference proteins) were searched against the five protein sets (R.philippinarum, C. gigas, L. gigantea, P. fucata, P. yessoensis) using BLASTP with the parameter "-e 1e-5". The genes were identified as candidates following criterion: 1) the genes have a similarity of more than 40% with the reference proteins; 2) the genes have the same domain structures with the reference proteins; 3) the annotation of the genes are same reference proteins. The motifs and domains identification of the genes were mainly performed with the Pfam and SMRT database. The genes were functionally annotated based on homologous searches against databases of NR database (NCBI).

Identification of the G protein-coupled receptor (GPCR) related genes
The protein-coding genes of R. philippinarum, P. yessoensis, P. fucata, C. gigas and L. gigantea were annotated based on homologous searches against the database of PFAM, NR and SwissProt.
The genes with the 7 transmembrane (TM) domain were considered as the GPCR-associated family genes, that were classified into different function according the NR and SwissProt annotation. The P-values were obtained from Chi-square tests for the significant expansion gene function in the R. philippinarum and the total number of each function was used as background.

R. philippinarum and immune challenge
Manila clams were collected from a local commercial farm and transported to the laboratory.

Real-time RT-PCR confirmation of Illumina sequencing data
To validate our Illumina sequencing data, differentially expressed genes were selected for quantitative RT-PCR (RT-qPCR) analys is, using the same treated groups and control groups RNA samples as were used for transcriptome profiling. All data were analyzed using the 2 -△△ Ct method.

Data availability
This genome project has been registered in NCBI under the BioProject accession PRJNA479743.
The R. philippinarum assembly has been deposited with the accession number QUSP00000000 under PRJNA479743. All sequence data of R. philippinarum genome have been deposited in NCBI Sequence Read Archive under the accession numbers of SRR7716263-SRR7716297.   Note: The recovered matches were classified as 'complete' (C) if their lengths were within the expectation of the BUSCO profile match lengths. If these were found more than once they were classified as 'duplicated'(D). The matches that were only partially recovered were classified as 'fragmented'(F), and BUSCO groups for which there were no matches that pass the tests of orthology were classified as 'missing'(M).The number of total BUSCO groups (n) was 843. Note: Repeats in R. philippinarum genome were consisted with tendam repeats and interpersed repeats. The TRF (Tandem repeats finder) was used to predict tandem repeats, RepeatMasker and RepeatProteinMask were used to predict interpersed repeats. Note: The genomic abundance of TEs were predicted by Repeatmasker through ab-initio method and homology based method. The TE proteins were predicted based on TE protein database by RepeatProteinMask. The final results of TEs in R. philippinarum genome were integrated of both methods without redundancy. "Other" refer to the repeats that can be classified by RepeatMasker, but not included by the classes above; "Unknown" refer to the repeats that can't be classified by RepeatMasker. Note: # Statistics calculated from the gene set predicted form each method. $ Statistics calculated from the gene set predicted by homolog proteins from each species. * Statistics calculated from the R.philippinarum genome.
Table S13 General statistics of protein-coding genes of homolog species. Related to Figure 1.

Species Number
Table S17 GO enrichment of unique gene families in R. philippinarum compared with four other Mollusca species (C. gigas, L. gigantea, P. fucata and P. yessoensis). Related to Figure 1.