Clam focal and systemic immune responses to QPX infection revealed by RNA-seq technology

The hard clam Mercenaria mercenaria is an important seafood species widely exploited along the eastern coasts of the United States and play a crucial role in coastal ecology and economy. Severe hard clam mortalities have been associated with the protistan parasite QPX (Quahog Parasite Unknown). QPX infection establishes in pallial organs with the lesions typically characterized as nodules, which represent inflammatory masses formed by hemocyte infiltration and encapsulation of parasites. QPX infection is known to induce host changes on both the whole-organism level and at specific lesion areas, which imply systemic and focal defense responses, respectively. However, little is known about the molecular mechanisms underlying these alterations. RNA-seq was performed using Illumina Hiseq 2000 (641 Million 100 bp reads) to characterize M. mercenaria focal and systemic immune responses to QPX. Transcripts were assembled and the expression levels were compared between nodule and healthy tissues from infected clams, and between these and tissues from healthy clams. De novo assembly reconstructed a consensus transcriptome of 62,980 sequences that was functionally-annotated. A total of 3,131 transcripts were identified as differentially expressed in different tissues. Results allowed the identification of host immune factors implicated in the systemic and focal responses against QPX and unraveled the pathways involved in parasite neutralization. Among transcripts significantly modulated upon host-pathogen interactions, those involved in non-self recognition, signal transduction and defense response were over-represented. Alterations in pathways regulating hemocyte focal adhesion, migration and apoptosis were also demonstrated. Our study is the first attempt to thoroughly characterize M. mercenaria transcriptome and identify molecular features associated with QPX infection. It is also one of the first studies contrasting focal and systemic responses to infections in invertebrates using high-throughput sequencing. Results identified the molecular signatures of clam systemic and focal defense responses, to collectively mediate immune processes such as hemocyte recruitment and local inflammation. These investigations improve our understanding of bivalve immunity and provide molecular targets for probing the biological bases of clam resistance towards QPX.


Background
The hard clam, Mercenaria mercenaria, is an ecologically-and economically-important marine bivalve species that thrives along the northeastern coasts of the United States and Maritime Canada. In the past few decades, the hard clam industry has been severely impacted by a protistan parasite called QPX (Quahog Parasite Unknown), which is responsible for mortality episodes in both wild and cultured clam populations [1][2][3][4][5][6][7]. QPX is believed to be an opportunistic pathogen and has been detected in a wide variety of environmental substrates including sediments, marine aggregates and other organic matrices [8][9][10]. Interestingly, previous reports highlight the ability of QPX to sustain very low abundance in clams without causing disease outbreaks until it encounters hosts with reduced immunity or following shifts of environmental conditions that add to the virulence of the parasite, under which conditions QPX can take advantage to establish infection sometimes leading to severe clam mortality events [4,11].
Lesions caused by QPX, usually associated with the presence of nodules, are commonly found in pallial tissues, such as alongside the inner rim of the mantle or at the base of the siphon [1,2]. These places are widely considered as the portal of entry for QPX cells acquired from the surrounding environment during suspensionfeeding [5,7]. The QPX nodules represent inflammatory masses containing both parasite cells and abundant clam hemocytes, resulting from a series of comprehensive host immune responses induced by the infection that leads to massive focal hemocyte infiltration, parasite encapsulation and partial necrosis of the affected area [1]. Like other invertebrates, the hard clam lacks the specific immune responses and their defense mechanism mainly relies on the effectors of innate immunity, which is mediated by circulating hemocytes and highly diversified humoral antimicrobial factors. These cellular and humoral immune components work in a synergistic way to initiate the recognition, segregation and ultimately elimination of pathogens and other non-self entities [12,13]. The launching of innate immune responses involves myriad cellular and humoral events modulated not only at the infection sites (focally) but also at a larger, whole-organism scale (systemically). In general, the focal response represents the alterations driven by direct host-pathogen interactions at the infection sites where direct cell-cell (e.g., molecular patterns) interactions mediate the response; while the systemic response reflects overall modifications within the host as a result of the ongoing infection and is mainly associated with dynamic changes of circulating hemocytes and their secreted immune mediators.
Most of the previous investigations have solely focused on the systemic response of M. mercenaria against QPX during the infection events, where changes in cellular and humoral immune parameters (e.g., anti-QPX activity and lysozyme activity in clam plasma, hemocyte phagocytic activity, reactive oxygen species (ROS) production, etc.) as well as expression of a limited number of immune-related genes in tissues and circulating hemocytes were assessed [11,[14][15][16]. In contrast, no previous studies have focused on the characterization of clam focal response at the infection sites. Given the fact that QPX disease is usually focal with formation of welldelimited lesions, the study of clam immune responses at the infection site in the lesions per se is of specific value as it provides insights to better characterize cellular interactions between the hard clam and QPX upon their contact. In this framework, QPX disease in clams offer a unique opportunity to contrast focal and systemic responses against microbial diseases in invertebrates allowing for a more comprehensive understanding of defense strategies used by these animals to fend microbial attacks.
Our study aimed to characterize the gene regulation features of M. mercenaria during QPX infection by profiling the transcripts at the infection lesion and compare focal clam responses with systemic responses detected in healthy tissues from infected clams (in addition to a parallel comparison with tissues from healthy clams). This study allowed the identification of factors involved in the interactions with the parasite as well as molecular pathways activated by the host to neutralize QPX.

Results and discussion
Illumina sequencing and de novo transcriptome assembly The main objective of this study was to identify molecular features of M. mercenaria in response to QPX infection and to compare the immune-related pathways involved in the lesion-specific focal response with the whole-organism scale systemic response. The highthroughput Illumina RNA sequencing and de novo assembly employed in this investigation allowed the construction of the transcriptome in the absence of M. mercenaria genome information. A total of 640,596,320 of 100 bp raw reads were generated from the Illumina paired-end sequencing with about 27 to 48 Millions paired-end reads generated from each of the 9 sequenced libraries (Table 1, Fig. 1a). The short read sequences generated from this RNAseq project have been deposited at the NCBI short Read Archive database under the SRA accession number SRP068241. Trimming and filtering procedures yielded 606,021,407 clean reads that were used for the de novo assembly of consensus transcriptome based on all sequenced RNA libraries in order to maximize the diversity of transcripts. This allowed 90.61 to 92.20 % of the reads from the 9 libraries be used for the transcriptome assembly. A total of 62,980 transcripts ranging from 201 to 23,103 bp with average size of 1297.59 bp and median size of 835 bp were produced from the assembly after low FPKM and rare isoforms filtering. The size distribution of all the de novo assembled transcripts is shown in Fig. 1b. Once the transcriptome was constructed, the 9 libraries were individually remapped to the 62,980 transcripts and resulted with 85.27 to 89.05 % of reads remapping. Theses counting data were then used for DE analysis.

Transcriptome functional annotation
The transcriptome annotation performed using Blas-t2GO returned a total of 19,107 transcripts (30.3 %) with significant BlastX homology matches to other sequences in NCBI nr database (E-value < 10E-03) (Fig. 1). Not surprisingly, the top 3 species that had the most similarity to M. mercenaria sequences were mollusks with available genomes and included the Pacific oyster, Crassostrea gigas (7,467), followed by the owl limpet Lottia  . 1c). KEGG Orthology (KO) terms were assigned to 6,425 sequences and reference pathways were mapped to the KEGG database based on the assigned KO terms (Fig. 1a, Additional file 1). A total of 29,815 sequences were identified to match to at least one conserved protein domain in the InterPro database (Fig. 1a, Additional file 1). Gene ontology (GO) assignments were used to classify functions of the predicted clam proteins. Based on sequence similarity (E-value of 10E-03), 13,584 sequences were assigned to at least one GO annotation (Fig. 1a, Additional file 1). As summarized in Fig. 2, a total of 8,168, 4,600 and 4,231 sequences were respectively categorized into the three main categories: biological process, cellular component, and molecular function at the second functional annotation level. The most dominant terms presented in the three categories are the "cellular process", "metabolic process", "binding", "catalytic activity", "cell", and "organelle". Very few transcripts were clustered into "rhythmic process", "cell killing", "protein tag", "channel regulator activity", "nucleoid" or "virion". It is noticeable that a good fraction of transcripts were clustered into the immunerelated categories of response to stimulus (503), immune system process (43) and biological adhesion (38). Those transcripts were of special interest given that they might be involved in the M. mercenaria defense and resistance toward QPX infection.
A significant portion (69.7 %) of M. mercenaria transcripts did not match any BlastX hit in NCBI nr database, in agreement with previous transcriptomic studies in mollusks [17][18][19][20]. Most of the unannotated transcripts may represent transcripts spanning untranslated mRNA regions, or transcripts containing only nonconserved protein domains [21,22].

Identification of differentially expressed transcripts
The generated transcriptome was used as a reference for downstream investigations of global gene expression in the three different tissues of interest (nodule, nonnodule and healthy) to identify genes associated with M. mercenaria's focal and systemic immune response against QPX. A gene-isoform relationship was estimated using RSEM over Trinity output isoforms. Results showed that about 43 % (27,021) of all the transcripts had 1 isoform, 19 % (12,307) had 2 isoforms and 38 % (23,652) had 3 isoforms, suggesting extensive isoform diversity in M. mercenaria transcriptome. Transcript isoform variation could affect mRNA stability, localization and translation, as well as the production of protein variants that differ in localization or function [23].
By comparing the number of transcripts expressed in each sample, the contribution of specific samples to the analysis can be estimated. The highest number of expressed transcripts was found in the nodules of infected clams, which were closely followed by that found in the healthy clam samples (Fig. 3a). The lowest number of expressed transcripts came from non-nodule samples of QPX infected clams, with about 1,500 less transcripts expressed than the other two samples. Read coverage, which is critical in accurate determination of fold change, averaged 477, 422 and 509 reads per transcript for nodule, non-nodule and healthy tissue samples, respectively (Fig. 3b).
Statistical analysis by DEseq identified 3,131 differentially expressed (DE) transcripts from the pair-wise comparisons (|log 2 (fold change)| >2, adjusted p-value <0.001) between clam tissue samples (Fig. 4). In nodules, a total of 829 transcripts, including 408 over-expressed and 421 under-expressed transcripts were identified as compared to non-nodule samples of QPX infected clams. Compared to tissues from healthy clams, 1,591 DE transcripts were identified in nodules with 864 overexpressed and 727 under-expressed transcripts. Similarly, 1,681 DE transcripts were obtained from the comparison between healthy and non-nodule clam tissues, of which 513 and 1,168 were over-and under-expressed, respectively (Fig. 4, Additional file 3). A total of 1,694 of these DE transcripts had protein homologs found in NCBI nr database by Blastx searches (e-value <10E-03), which were further examined for their putative functions during  Table 1 for more details). The Y-axis indicates the number of transcripts expressed in the samples. b Summary statistics of the reads coverage in each sample M. mercenaria immune response toward QPX. Annotated transcripts were subsequently grouped into curated categories according to their biological functions based on the gene ontology (GO) terms and literature searches highlighting immune functions. DE transcript sets were further examined in reference to the assigned KO terms for the analysis of pathways regulation. Significantly enriched KEGG pathways were identified via the Fisher's exact test (P < 0.01) [24].  Table 1 for more details). For B and C: the cross-shaded areas inside each bar represent the number of transcripts with higher (B) or lower (C) expression levels in a condition as compared to the other two (e.g., the expression levels of 38 focally over-expressed transcripts in B were higher in nodule than in non-nodule tissues and healthy clams) Here we specifically focus on DE transcripts drawn from comparisons between "nodule vs. non-nodule" and "non-nodule vs. healthy" tissue samples, which were respectively considered to reflect the transcriptomic changes caused by "focal" and "systemic" immune response of M. mercenaria toward QPX infection, respectively. The overview of DE transcripts drawn from these two responses is presented in Fig. 4 and Additional files 2 and 3.

Differential expression of immune related transcripts Focal response
Clam focal response reflected the alterations caused by direct clam-QPX interaction at the infection site (Additional file 2). QPX nodules are inflammatory masses resulting from massive hemocytes infiltration and encapsulation of parasite cells [1,2]. This process largely relies on the motility and adhesion properties of hemocytes, thus allowing these cells to migrate throughout the circulatory system and recruit to the infection site. Hemocytes can sense stimuli in host tissues through an array of cell surface receptors, and use these cues to adjust their behavior accordingly [25]. The activation of hemocytes requires the binding of specific ligands to the cell surface receptors, which subsequently initiate the transduction of extracellular signals into the cytoplasm via a variety of signaling pathways, thus inducing a series of hemocytemediated immune response such as phagocytosis, encapsulation (prominent response against QPX in clams), ROS production, as well as secretion of immune effectors and cytokines [26][27][28]. A collection of DE transcripts involved in these defense processes was identified during focal response, suggesting that strong and comprehensive hostpathogen interactions were taking places inside the QPX lesions ( Table 2, Additional file 2). A large fraction of the DE transcripts of focal response was annotated as receptors or molecules with receptor activities, which putatively contribute to the host defense against QPX as (1) cell surface receptors expressed on hemocytes that mediate the recognition and phagocytosis or encapsulation of foreign entities through microbe associated molecular patterns (MAMP); (2) signaling receptors activating intracellular signaling cascades or (3) the soluble bridging molecules mediating the linkage between MAMPs and hemocytes [29]. Among those receptors, most are identified as pathogen pattern recognition receptors (PRRs), which include the C-type lectins (CTLs), the scavenger receptors (SRs) and the toll-like receptors (TLRs).
The C-type mannose receptor-2 (MRC2) identified during the focal response (Table 2) is a member of the C-type lectins (CTLs) superfamily, a large group of Ca 2+ -dependent carbohydrate-binding proteins that play crucial roles in innate immunity. CTLs recognize pathogens and facilitate their phagocytosis [30,31] or encapsulation [32][33][34]. MRCs are also key regulators of inflammatory responses and contribute to the removal of harmful inflammatory agents [35,36]. The 16-fold over-expression of MRC2 during the focal response suggested that active hemocyte encapsulation and local inflammation was induced by QPX at the infection lesions, which is consistent with the results of histopathological observations [1]. Another over-expressed CTL member, the perlucin-like protein, has been previously shown to trigger immune response in Manila clams during microbial infection [37].
Scavenger receptors (SRs) were also among the strongly over-expressed transcripts in infection foci ( Table 2). These included somatomedin-b and thrombospondin type-1 domain-containing (RPE spondin), insulin-related peptide receptor, hemicentin-1, lysyl oxidase-like protein 2 and mam domain-containing glycosylphosphatidylinositol anchor protein 1. SRs are structurally diverse PRRs that share the common function of recognizing oxidized or acetylated low-density lipoprotein (LDL) [38]. They contribute to innate immunity by recognizing MAMPs and mediating non-opsonic phagocytosis [29,39]. They are extensively found on immune cells and are able to interact with both modified-host components and exogenous ligands, which makes SRs a key component in host defense, apoptosis, inflammation and lipoprotein homeostasis [40][41][42]. For example, scallop SRs bind not only to acetylated LDL but also to MAMP including lipopolysaccharides (LPS), peptidoglycans (PGN), mannan and zymosan particles [39]. The sea urchin genome encodes approximately 150 genes consisting of one or more scavenger receptor cysteine-rich (SRCR) domains [43], and the members of this gene family exhibit dynamic shifts in transcription after immune challenge [40,44].
Our results also show an over-expression of TLR-1 and Toll-8/tollo in nodules (Table 2), which is in agreement with previous investigations showing up-regulation of TLRs in M. mercenaria mantle following QPX challenge [15]. TLRs are among the most ancient and conserved PRRs. They are expressed by immune cells and interact with a large variety of MAMPs. Bivalve TLRs have been characterized in the oyster C. gigas and the scallop C. farreri where they exhibited significant response to LPS stimulation [45,46]. Transcriptional modulation of TLRs has also been reported in Ruditapes philippinarum and Mya arenaria following MAMPs stimulation and bacterial challenge [47,48]. Interestingly, Toll-8 (Tollo) has been shown to participate in Drosophila epithelial immunity where it mediates host cells communication that subsequently activates systemic immune responses [49]. This suggests that the Toll pathway could be one of the crucial pivoting links that allow coordination between focal and systemic immune components during infection.
The QPX nodules are formed as the result of granulomatous inflammation, which is a chronic inflammatory Table 2 Transcripts with annotated functions (GO terms) related to immune recognition, signaling and regulation that were differentially expressed during M. mercenaria focal response against QPX. Additional information on these transcripts is given in Additional file 2. "Inf" designates an infinite fold change calculated for focal response as the expression of that transcript in non-nodule tissue was equal to 0   [50,51]. The formation of granuloma requires local recruitment of hemocytes at the site of infection to execute extracellular defense processes around the invaders [50]. An array of transcripts associated with cell migration, adhesion and proliferation was regulated in nodules, including Gprotein coupled receptors (GPCRs) and integrins families (Table 2). GPCRs regulate inflammatory response via binding to chemokines and chemoattractants, thus activating pathways mediating hemocyte migration and adhesion [52]. They also activate transcription factors in immune cells, thus modulating the synthesis and secretion of certain pro-or anti-inflammatory substances [53]. On the other hand, integrins represent a major group of cell adhesion mediators [54]. They not only modulate the cell-cell and cell-extracellular matrix adhesion, but also affect multiple signal transduction cascades regulating cell survival and proliferation [54].
Overexpression of GPCRs and integrins in nodules suggests their role in hemocytes adhesion and aggregation associated with the formation of granuloma [50,55]. Several enzymes regulating ROS production were also over-expressed during focal response (Table 3). These included a dual oxidase, which is a key component mediating host-microbe interactions in mucosa [56,57]. Dual oxidase regulates oxidative burst and ROS production in the gill muscosa of the shrimp Marsupenaeus japonicus, favoring shrimp survivorship during viral infections [58]. Interestingly, transcripts of dual oxidase were only expressed in nodules, suggesting this enzyme was induced upon direct clam-QPX interactions as a part of the mantle mucosa-related immune response. Other transcripts associated with oxidation-reduction processes also exhibited somewhat nodule-exclusive pattern, including the allene oxide synthase-lipoxygenase (AOSL), lysyl oxidase-like protein (LOXL), ww domaincontaining oxidoreductase (WWOX), c-terminal binding protein (CtBP), isocitrate dehydrogenase (ICD) and methylenetetrahydrofolate reductase (MTHFR). These molecules are important for maintaining the redox homeostasis of extracellular environment as they are key regulators for oxi-reduction reactions. Over-expression of these transcripts in nodules suggests the need for the host to timely balance out excessive ROS and other toxic intermediates produced during interaction with QPX. In addition to redox-regulation, many of these molecules also take part in the immune modulation indirectly. For example, AOSL play a role in coral immunity by controlling the production of the inflammation regulator arachidonic acid during apoptosis [59], and LOXL acts both as a scavenger receptor and regulator for extracellular matrix remodeling that initiate hemocyte migration and tissue regeneration [60], while WWOX was shown to promote proliferation of immune cells through inhibition of their apoptosis [61,62]. In addition, ICD, MTHFR and cytochrome p450 are major detoxification enzymes [63,64]. In fact, immune cells and their secreted effectors require the proper redox state in the extracellular environments to exert their immune functions, which makes the maintenance of redox homeostasis essential for persistent and effective host defense [65,66]. This is particularly true in the case of QPX disease where the neutralization of parasites depends on extracellular killing pathways [26].
Apoptosis is an essential host mechanism to effectively remove damaged and infected cells without causing inflammatory destructions to surrounding tissues [67,68]. Interestingly, apoptosis seems to be largely inhibited during M. mercenaria focal response, as shown by the under-expression of pro-apoptosis transcripts ( Table 2). For example, the tumor necrosis factor (TNF)-like protein and 3-hydroxy-3-methylglutaryl-coenzyme A (HMG-CoA) reductase-like protein were significantly under-expressed in nodules. Similarly, the pro-apoptotic p53-induced protein and solute carrier family 25 member protein were also under-expressed in nodules. Meanwhile, inhibitor of apoptosis protein (IAP) was overexpressed in nodules. IAPs regulate immune cell expansion and survival in highly inflammatory environments in mammals [69] and they may share similar function in clams by preventing hemocytes from death during Table 2 Transcripts with annotated functions (GO terms) related to immune recognition, signaling and regulation that were differentially expressed during M. mercenaria focal response against QPX. Additional information on these transcripts is given in Additional file 2. "Inf" designates an infinite fold change calculated for focal response as the expression of that transcript in non-nodule tissue was equal to 0 (Continued)  interaction with QPX. In fact, ROS production during parasite killing may trigger apoptotic cell death in molluscs [67], and proper control of apoptosis mechanisms is required to maintain cellular homeostasis during immune response. This suspected inhibition of host apoptosis is supported by the above-mentioned overexpression of integrins, as these were shown to protect cells from apoptosis and induce anti-apoptotic pathways during cell adhesion and spreading in the snails Lymnaea stagnalis [70] and B. glabrata [71]. Infection and tissue injury trigger host immune responses via immune signaling pathways [72], by activating transcription factors and initiating the production of immune effectors and regulators. Immune signaling pathways identified in mollusks include Toll, MAPK/ JNK and JAK/STAT signaling pathways [27,47,73,74]. During M. mercenaria focal response to QPX, a variety of transcripts encoding kinases and phosphatases were over-expressed (Table 2), suggesting the involvement of MAPK and other kinase-mediated cascades in regulating the focal inflammatory response [75], whereas the under-expression of EF-hand domain containing protein and calmodulin may indicate the suppression of calcium-regulated pathways [76]. Overexpression of E3 ubiquitin-protein ligase and its upstream regulator COP9 signalosome suggests the activation of the damage surveillance ubiquitin/proteasome pathway [77][78][79]. In parallel, the over-expression of rho GTPase and rho kinase suggests the induction of anti-apoptotic Rhomediated signaling pathway [80,81] and reinforce the idea that apoptotic inhibition is extensively initiated by M. mercenaria to help fight QPX. However, a very limited number of DE transcripts was detected in relation to those conventional signaling pathways of innate immunity, such as the complement pathway and the Toll/ TLR pathway. Only 1 transcript encoding macrophageexpressed gene protein 1 (MPEG1), a putative member of the complement pathway [82,83], was differentially expressed in QPX nodules. As for the Toll/TLR pathway, only a few receptors were identified ( Table 2) but none of their downstream components.
Interestingly, components of the Notch and the Wnt signaling pathways were over-expressed during focal response ( Table 2). These included two putative Notch family members, the neurogenic locus Notch protein and the mediator protein nicastrin, and the tyrosin- Table 3 Transcripts with putative functions (GO terms) related to immune effectors that were differentially expressed during M. mercenaria focal response against QPX. Additional information on these transcripts is given in Additional file 2. "Inf" designates an infinite fold change calculated for focal response as the expression of that transcript in non-nodule tissue was equal to 0 (Continued) protein kinase RYK which belongs to the Wnt pathway. Wnt signaling pathway regulates many cellular immune processes and is evolutionarily conserved across taxa [84,85]. Wnt signaling has been intensively exploited for its regulatory functions during wound healing and tissue regeneration [86,87], so its over-expression may be related to wound healing to repair damage resulting from tissue digestion by the parasite or tissue necrosis. The Notch signaling network was scarcely explored in bivalves even though it has been reported to be associated with several aspects of immune response in mammals [88], especially in regulating granulomatous reactions to foreign bodies [89]. As an evolutionarily conserved pathway involved in modulating the intercellular signaling, the Notch pathway presumably shares an equally important role in M. mercenaria by modulating the formation of granuloma. Notch signaling triggers macrophage expression of genes involved in pro-inflammatory responses [90], but can suppress inflammation responses triggered by canonical TLR cascade [91], in agreement with our observations. These tightly regulated mechanisms ensure tailored immune responses against different pathogens and are crucial for the host to achieve high immune efficiency while avoiding excessive immune activation and self-inflicted damages. Several proteases (A.K.A proteinases, peptidases) were also differentially regulated during M. mercenaria focal response against QPX. These proteases mostly belong to the serine and metallo protease families and were generally over-expressed in QPX nodules (Table 3). Proteases serve as key immune modulators partially through their ability to digest and remodel the extracellular matrix and tissues associated with hemocyte activation [92]. Commonly associated with lysosomes and granules of inflammatory cells, serine proteases participate in immune regulation either directly by degradation of pathogens or indirectly through activation of cell surface receptors and signal molecules [92][93][94]. The function of metalloproteases in immune regulation is even more diverse, acting as immune effectors, signal transducers, and mediators of immune cell development and migration [95]. Metalloproteases are also known to be involved in many pro-inflammatory pathways, particularly in the Notch pathway where they act as a type of downstream element to Notch [96].
At the same time, immune effectors with known universal protease inhibitor activities, such as alpha2macroglobulin (α 2 M) and thioester-containing protein (TEP, a subfamily of α 2 M), were also collectively overexpressed in QPX nodules ( Table 3). The α 2 M superfamily inhibits peptidases of diverse origins [97,98]. The simultaneous over-expression of proteases and protease inhibitors may reflect a finely adjusted defense response of M. mercenaria to maintain homeostasis and regulate self-and pathogen-derived proteases, as shown in other host-pathogen systems [92,99], including bivalves [100]. Proteases have been identified as major virulence factors of QPX [101,102], and are thought to degrade host proteinaceous immune effectors and hydrolyze host tissues to fulfill nutritional requirements. Therefore, inhibition of pathogen proteases contribute to host protection, and was shown to represent a determinant factor for resistance against infectious diseases in bivalves [103][104][105][106]. In parallel, TEPs have been extensively studied in mollusk immunity [99,107], and beside their function as protease inhibitors, they also play a role as PRRs or opsonins to facilitate microbial phagocytosis and encapsulation. Consistent over-expression of TEPs was noted in this study in agreement with findings following experimental infection with QPX [15]. These results support a critical role of TEPs in clam immune response against QPX either via the protease-inhibitor activity of these proteins, or by mediating parasite encapsulation, or both.
In addition, several metal ion transporters were overexpressed in nodules (Table 3), including the putative copper ion binding protein ceruloplasmin precursor, the transferrin enzyme ferric-chelate reductase and the divalent metal transporter (A.K.A natural resistanceassociated macrophage protein; 328 fold increase). These molecules contribute to host defense by controlling the supply of essential micronutrients in the vicinity of infection sites thus reducing parasite survival [108,109] and favoring the production of antimicrobial factors [110]. It is noteworthy to point out that some of the focally over-expressed transcripts might be partly driven by the dramatic increase of hemocyte proportion within tissues in the vicinity of infection foci as compared to the surrounding host tissues [26,111]. This would be especially the case for transcripts known to be highly expressed in hemocytes, such as the cell surface PRRs, secreted humoral immune effectors, cell signal transducers and enzymes associated with ROS production.

Systemic response
Significant transcriptomic regulations were observed during M. mercenaria systemic immune response against QPX, with a total of 1,681 DE transcripts, which is about two times the number of focal DE transcripts (829). However, only about one third of the DE transcripts (513) were over-expressed in response to QPX infection, the larger remaining part (1,168) represented significantly under-expressed transcripts (Fig. 4, Additional file 3), possibly due to the chronic stress imposed by the infection. Transcriptome-wide depression has been demonstrated in many marine invertebrates as the result of pathogenic or environmental stress [18,59,[112][113][114]. A considerable number of systemically underexpressed transcripts was related to metabolism and biomineralization, which could be the result of host resource allocation during on-going infection. The chronic inflammation induced by infection likely created extra energy demands, which require resources being allocated from other physiological processes, such as growth and reproduction, to immune processes which is critical for the survivorship of the host. In fact, slow growth and lower tissue conditions are frequently observed in QPXinfected M. mercenaria [2], and similar energy trade-offs existed between immune defense and other energy expenditure pathways [115][116][117][118]. Interestingly, a suite of transcripts over-expressed during focal response was significantly under-expressed in non-nodule tissues as compared to healthy clams (Additional file 3), which included several immune effectors and mediators associated with nodule formation and focal inflammation (e.g., integrins, notch proteins and peroxidases). In fact, maintaining high levels of these focally-induced molecules could be costly and dangerous as some are toxic to both the parasite and the host, so their production must be restrained within areas where they can directly exert defense function, and reduced outside the infection foci to minimize risks of undesirable effects on the host [116,117].
On the other hand, systemically over-expressed transcripts included stress proteins and other soluble immune factors such as lysozyme (c-type lysozyme 2), lectins (C1q domain containing protein, macrophage mannose receptor 1, low affinity immunoglobulin epsilon fc receptor), AMP (hemocyte defensin), proteases (cathepsin K, calpain 11, isoaspartyl peptidase/L-asparaginase, ASRGL, counting factor associated protein d) and ferric-chelate reductase (Table 4, Additional file 3). Over-expression of host stress proteins, such as heat shock proteins (HSP 70, HSP 90) and universal stress protein (USP) was also noted, in agreement with observations made during infection in other bivalve species [119][120][121]. Increased levels of stress proteins provide host cells with protection against incorrect protein folding caused by infection, inflammation, oxidative stress and other destructive events [122,123]. The systemic over-expression of soluble immune effectors (e.g., humoral proteins) may help maintaining comparatively high immune capacity to prevent the spread of QPX (or secondary pathogens) throughout the host. In addition, transcripts of antiapoptotic factors (IL17, deoxyguanosine, baculoviral map repeat-containing proteins) were also overexpressed during the systemic response, indicating that anti-apoptotic processes noted during focal response are not limited to the infection foci.
Pathway alterations during M. mercenaria's response to QPX Transcriptomic alterations during both focal and systemic response discussed above were also highlighted in the pathway enrichment analysis of the DE transcripts. This analysis aims at extracting an overview of phenotypic changes on the underlying functional level, to reduce the complexity of biological information given by the long lists of DE genes/transcripts [124]. The KEGG pathways of focal adhesion (04510), ECM-receptor interaction (04512), Notch signaling pathway (04330) and apoptosis (0421) were significantly over-represented during both focal and systemic response (Fig. 5), even though fold enrichment were generally higher during the focal response. Other immune-related pathways particularly enriched during the focal response included regulation of actin cytoskeleton (04810), cell adhesion molecules (CAMs, 04514), the leukocyte transendothelial migration (04670), complement and coagulation cascade (04610) and Wnt signaling pathway (04310). These pathways are critically involved in the immune cells activation during migration, attachment and parasite encapsulation, which serve as the underlying mechanisms for nodule formation and QPX killing. On the other hand, basic metabolic pathways such as the citrate cycle (TCA) and pyruvate metabolism were specially enriched during the systemic response. These alterations were largely in accordance with the under-expression of metabolism-associated DE transcripts in infected tissue compared to the healthy tissue (Fig. 4), possibly reflecting changes in the energy allocation strategy during infection as discussed above.
Interdependence of KEGG pathways widely exist and most of these are interrelated with each other via shared components, forming a signaling network to allow for pathway crosstalk. To investigate these interactions, we extracted the DE transcripts shared by multiple enriched pathways and constructed a sketch of the hypothetical pathways network that are significantly altered by QPX infection (Fig. 6). In this framework, M. mercenaria response to QPX infection was initiated upon the sensing of danger signals via cell membrane receptors. The signals subsequently transmitted down through the MAPK, Wnt and Notch pathways and triggered the production of a series of host defense factors as the end results. In parallel, activation of pathways regulating actin cytoskeleton and leukocyte transendothelial migration facilitated the recruitment of hemocytes to the infection area to build a barrier of cellular defense against the parasite. Recruited hemocytes then attached to and encapsulated QPX cells as suggested by the modulation of focal adhesion and ECM receptor interaction pathways. These cellular activities were performed under a tight regulation of the apoptosis pathway to determine cell fates, resulting in either the survival or death of M. mercenaria cells.

Distinctive transcriptomic pattern of healthy clams
A suite of transcripts (407, Fig. 4, Additional file 4 and Table 5) exhibited higher transcription levels in healthy clams as compared to diseased clams (considering both nodule and non-nodule tissues). A considerable fraction of these transcripts were related to metabolic processes, nucleic acids binding and transcriptional regulation. The over-expressed immune-related transcripts identified in healthy clams are of particular interest as they may be involved in M. mercenaria resistance towards QPX. For example, an antimicrobial protein (aplysianin A [125,126]) was exclusively identified in healthy clams with almost no detection in diseased clams. In addition, the highest expression levels of a serine protease inhibitor were also observed in healthy clams. A serine protease inhibitor has been linked to oyster (Crassostrea virginica) resistance against the protozoan parasite Perkinsus marinus [103][104][105][106]127]. Therefore, the high expression of the serine protease inhibitor in healthy clams supports its involvement in M. mercenaria's resistance against QPX, likely by inhibiting the activity of parasite proteases. Moreover, a pathogen recognition protein (ctype lectin domain family 10 member A-like) was significantly higher in healthy clams as compared to diseased animals, which may also contribute to clam resistance against QPX by promoting microbial recognition and encapsulation [128]. Previous studies demonstrated that clam genetic background affects M. mercenaria resistance toward QPX [4,129,130], therefore immune- Table 4 Selected transcripts with annotated functions (GO terms) related to immune response that were over-expressed during M. mercenaria systemic response against QPX. Additional information on these transcripts is given in Additional file 3. "Inf" designates an infinite fold change calculated for focal response as the expression of that transcript in non-nodule tissue was equal to 0  Fig. 4 and Table 5) exhibited lowest expression levels in healthy clams. These include stress protein HSP 70 and USP, as well as the detoxification molecule cytochrome p450, which together highlight the stress experienced by clams as the result of QPX infection [63,[131][132][133][134][135]. In addition, other immune related transcripts including a protease (isoaspartyl peptidase/L-asparaginase, ASRGL), a protease inhibitor (GTP-binding protein yptV4), a Ctype lectin (MRC1), and molecules involved in tissue regeneration and cell signaling, were also underexpressed in healthy clams as compared to infected individuals, suggesting their role in fighting the infection.

Conclusions
This is one of the first studies contrasting focal and systemic immune responses to infections in invertebrates using high-throughput sequencing. Resulting transcriptome represents a significant addition to the so far limited public genomic information available for M. mercenaria. The transcriptomic profiles of healthy and infected clams reflected complex interactions between the host immune system and the pathogen leading to molecular changes at both the infection foci and the mercenaria reflected prevailing transcriptomic suppression accompanied with a contrasting over-expression of stress proteins and soluble antimicrobial effectors; whereas the focal response highlighted cell-cell interactions between hemocytes and the parasite that typically result in local inflammation, extracellular degradation, encapsulation, granuloma formation, and wound repair. What needs to be kept in mind is that the regulation of these genes can be the result of either an effective immune response or a symptom of a future death. In fact, several apoptotic and anti-apoptotic genes were regulated. This highlights a fine adjustment of M. mercenaria defense mechanisms to precisely adapt to the infection development (e.g., through the modulation of energy allocation, apoptotic and antiapoptotic processes and mobilization of different signaling pathways). In addition, the identification of immune-related transcripts that were particularly associated with healthy clams offered new perspectives on the molecular features putatively involved in clam resistance against QPX.

Clam tissue and RNA samples preparation
Adult hard clams (54 ± 5 mm in length, mean ± standard deviation) were collected from a QPX enzootic area in Massachusetts. Clams were shucked and grossly examined for the presence of nodules along the mantle rim. Nodule tissues were dissected and a small piece of each nodule was microscopically examined (fresh mount) for the identification of QPX cells. Aliquot tissues from positive nodules were submitted to RNA extraction. Meanwhile, a piece of seemingly healthy tissue that is anatomically symmetrical to the nodule was collected from the same clam (e.g., "healthy" tissue from a diseased clam) and divided into 2 aliquots, with the first aliquot used for QPX detection by qPCR [10] and the second used for RNA extraction. Mantle tissues were also collected from "seemingly" healthy clams (no visible nodules) and used for QPX detection and RNA extraction. Following confirmation of disease status by qPCR, the samples were divided into the following 3 categories: (1) infected tissue from a diseased clam, (2) "healthy"(non-nodule) tissue from a diseased clam, and (3) "healthy" tissue from a healthy clam. Total RNA was individually extracted using TRIzol ® Reagent (Invitrogen, Carlsbad, CA, USA). Further RNA clean-up and oncolumn DNase digestion were performed with RNeasy Mini Kit (Qiagen Ltd., Crawley, UK) according to the manufacturer's guideline. RNA quantity and quality were analyzed on Nanodrop ND-1000 (Thermo Scientific, Wilmington, USA). Only RNA samples with absorption ratios of A260/A280 close to 2.0 were used for RNA-seq analysis. A total amount of 3 μg RNAs per sample was used and pooled into a total of 3 pools (3 clams per pool) representing each tissue category (9 pools in all). The pooling strategy is shown in Table 1.

RNA sequencing, de novo assembly and annotation
The sequencing of each pooled RNA sample as a paired end (PE) reads library (100 pb) was performed on Illumina HiSeq 2000 platform at the McGill University and Genome Quebec Innovation Center (Montreal, Canada), producing from 27.1 up to 46.8 millions of PE reads per sample (Table 1). Raw reads were filtered and trimmed according to length and quality score (min length 60 nt, end trimming quality 25, min quality filtering: 20 on 75 % of the read length) using the FASTX-Toolkit software v 0.0.13 (http://hannonlab.cshl.edu/fastx_toolkit/). rRNA cleaning was performed using the riboPicker software v 4.0.3 (http://ribopicker.sourceforge.net) [136] against SILVA database v111. And finally pair retrieval was performed using a homemade python script accessible through GitHub (https://github.com/ppericard/bioinfotoolkit). High quality filtered sequence reads from all libraries were combined and subsequently used for de novo assembly. Assembly was based on the de Bruijn graph assembler Trinity (http://trinityrnaseq.sourceforge.net/) [137] Table 5 Selected transcripts with annotated functions (GO terms) related to immune response that were differentially expressed in naïve M. mercenaria as compared to QPX infected individuals. Additional information on these transcripts is given in Additional file 4. "Inf" designates an infinite fold change calculated for focal response as the expression of that transcript in non-nodule tissue was equal to 0 using the default parameters. Assembly quality was controlled by re-mapping raw reads back to the transcripts using bowtie 2 and RSEM within Trinity package scripts. The assembled transcripts with sequence length longer than 200 bp, re-mapping FPKM (fragments per kilobase of transcript per million mapped reads) greater than 1 and isoform discovery level greater than 1 % were then considered for the following annotation and transcript abundance quantitation. Annotation of this de novo assembled transcriptome was performed using Blast2GO (http:// www.blast2go.com) with a semi-automated functional annotation based on sequence homology search. Putative gene identities was obtained by Blastx search against National Center for Biotechnology Information (NCBI) nonredundant sequences (nr) database with the E-value threshold setting at 1E-03. Putative gene functions were predicted by sequence similarity search against Gene Ontology (GO, http://www.geneontology.org/) database and assigning GO annotation terms to each mapped transcript. Protein domain search and enzyme annotation were also performed using InterPro scan and the Kyoto Encyclopedia of Genes and Genomes (KEGG). KEGG Orthology (KO) terms and KEGG pathways were also assigned to the assembled sequences using the online KEGG Automatic Annotation Server (KAAS, http://www.genome.jp/kegg/kaas/) using bidirectional best-hit methods [138]. This server provides KO annotation and pathway mapping.

Differential gene expression analysis
The consensus transcriptome generated in the previous steps was used as a reference for transcript abundance analysis. Reads from each single library were separately mapped to this reference transcriptome using the "RNA-Seq by Expectation-Maximization (RSEM)" method that is bundled within the Trinity package. The expression level of each transcript was determined as the total mapped reads count. The differences in gene expression between clam tissue samples (nodule, non-nodule and healthy tissues) were estimated using the DESeq Bioconductor package (https://github.com/Bioconductor-mirror/SESeq/tree/release-3.2) [139] in R statistical software (R Development Core Team, 2010; http://www.R-project.org). The threshold for defining significant differentially expressed (DE) transcripts between two different conditions (3 replications in each condition) was set as adjusted p-value smaller than 0.001 and absolute log 2 (fold change) values greater than 2. Expression patterns of DE transcripts were also analyzed by a K-means clustering method using Euclidean distance based on expression levels over all input samples. For further analysis, only those DE transcripts with annotation were considered as candidates of interest and were subsequently divided into over-and under-expressed groups.

Availability of supporting data
The data sets supporting the results of this study are available at the NCBI short Read Archive database under the SRA accession number SRP068241.