Dual transcriptomics reveals co‐evolutionary mechanisms of intestinal parasite infections in blue mussels Mytilus edulis

On theoretical grounds, antagonistic co‐evolution between hosts and their parasites should be a widespread phenomenon but only received little empirical support so far. Consequently, the underlying molecular mechanisms and evolutionary steps remain elusive, especially in nonmodel systems. Here, we utilized the natural history of invasive parasites to document the molecular underpinnings of co‐evolutionary trajectories. We applied a dual‐species transcriptomics approach to experimental cross‐infections of blue mussel Mytilus edulis hosts and their invasive parasitic copepods Mytilicola intestinalis from two invasion fronts in the Wadden Sea. We identified differentially regulated genes from an experimental infection contrast for hosts (infected vs. control) and a sympatry contrast (sympatric vs. allopatric combinations) for both hosts and parasites. The damage incurred by Mytilicola infection and the following immune response of the host were mainly reflected in cell division processes, wound healing, apoptosis and the production of reactive oxygen species (ROS). Furthermore, the functional coupling of host and parasite sympatry contrasts revealed the concerted regulation of chitin digestion by a Chitotriosidase 1 homolog in hosts with several cuticle proteins in the parasite. Together with the coupled regulation of ROS producers and antagonists, these genes represent candidates that mediate the different evolutionary trajectories within the parasite's invasion. The host–parasite combination‐specific coupling of these effector mechanisms suggests that underlying recognition mechanisms create specificity and local adaptation. In this way, our study demonstrates the use of invasive species’ natural history to elucidate molecular mechanisms of host–parasite co‐evolution in the wild.

Shine, 2012; Weclawski et al., 2013), without the need to raise multiple generations in the laboratory to use in time-shift experiments (i.e., in Brockhurst et al., 2007;Schulte, Makus, Hasert, Michiels, & Schulenburg, 2010). Unlike time-shift experiments, however, hosts at invasion fronts are at the start of their co-evolutionary interaction, whereas parasites have been moving across the host population for some time near the invasion origin. It is nevertheless possible to test the mode of evolution along invasion pathways or to test for local adaptation, by taking hosts and parasites from separate locations along the historical invasion front, thus creating a sympatry vs. allopatry contrast. For instance, the same mode of parasite evolution after invasion took place repeatedly in rabbits infected with the Myxoma virus in Australia and the United Kingdom, leading to parallel changes in viral virulence across locations (Kerr et al., 2015). Experiments with invasive species can thus address questions about the adaptive evolution of phenotypes in nature, especially for nonmodel organisms with long generation times that prohibit experiments over multiple generations (Colautti & Barrett, 2013).
The invasion of the parasitic copepod Mytilicola intestinalis into the Wadden Sea unites all aspects of an invasion to study host-parasite co-evolution (Feis et al., 2016). Native in the Mediterranean and Adriatic Seas, M. intestinalis invaded mussel beds in the late 1930s (Caspers, 1939) and has spread southwest towards the island of Texel (The Netherlands) (Drinkwaard, 1993) and northward to Sylt (Germany) (Dethlefsen, 1972), arriving at both islands around the same time in the early 1970s. The invasion routes are thus replicates of a newly begun co-evolutionary relationship between native blue mussel Mytilus edulis populations and the parasite M. intestinalis. Although the parasitic lifestyle of M. intestinalis was sometimes debated (Davey, Gee, Bayne, & Moore, 1977), we could already show that M. intestinalis and mussels have evolved different evolutionary trajectories along their invasion routes. Specifically, sympatric combinations of M. intestinalis and mussel hosts caused more damage to the host than allopatric combinations (Feis et al., 2016), probably owing to fact that Mytilicola parasites mainly consume host tissue (Goedknegt et al., 2017). Mussels from the southwest invasion front on Texel, for example, were more resistant than those from Sylt, while mussels from Sylt tended to be more tolerant to being parasitized (Feis et al., 2016).
To understand the evolutionary process leading to the phenotypic differences, the molecular mechanisms in the mussel-Mytilicola relationship need to be deciphered. The major immune pathways are reasonably well characterized in bivalves and gastropods (Gerdol, 2017). For mussels (Mytilus spp.), these were mainly derived from transcriptional and metabolite studies in response to microbial pathogens, for example, Gram-positive and Gram-negative bacteria (Liu et al., 2014;Toubiana et al., 2014;Venier et al., 2011), yeast and filamentous fungi (Sonthi et al., 2012). Because -omics tools were essential in the description of important pathways in mussel immunity (Gerdol & Venier, 2015), this approach can also be useful to characterize the immune response towards macroparasitic infections. Macroparasites have been targeted only by very few mussel studies (e.g., trematode infections, Gorbushin & Iakovleva, 2011), despite a high parasite richness of this widespread epibenthic filter feeder (Thieltges, Engelsma, Wendling, & Wegner, 2013). Methodologically, macroparasite infections have the advantage that host and parasite transcriptomes can be sequenced separately, extending a dual transcriptomics approach to a dual-species transcriptomics approach. In microbial infections, hosts and parasites often cannot be separated physically and classical dual transcriptomics relies on bioinformatic separation of host and parasite transcripts leading to heavy bias towards the host with only few parasite transcripts (Greenwood, Ezquerra, Behrens, Branca, & Mallet, 2016). Nevertheless, a dual transcriptomics approach can be very helpful to identify interacting genes (Schulze, Schleicher, Guthke, & Linde, 2016;Westermann, Gorski, & Vogel, 2012) and has been applied to several model organisms (e.g., Choi, Aliota, Mayhew, Erickson, & Christensen, 2014;Pittman, Aliota, & Knoll, 2014;Rosani et al., 2015;Tierney et al., 2012). For bivalves, the interaction between the Pacific oyster Magallana gigas (previously known as Crassostrea gigas, see Salvi & Mariottini, 2017) and Ostreid Herpesvirus type 1 (OsHV-1) has been characterized using this approach, revealing that OsHV-1 interfered with active host defence by suppression of cytokine signalling through interferon-related pathways (Rosani et al., 2015).
So far, these dual transcriptomics studies have mainly focused on the infection process as such, but have not incorporated the geographic aspects of co-evolutionary rapid local adaptation or contrasted sympatric and allopatric combinations of hosts and parasites.
Here, we incorporate these aspects by combining the natural history of the Mytilicola invasion with a dual-species RNA seq analysis of a reciprocal infection experiment. This approach can uncover the molecular underpinnings of the mussel-Mytilicola relationship at the transcriptional level and thus help explain the different evolutionary trajectories observed along the invasion route of Mytilicola (Feis et al., 2016). We characterize this interaction on three different levels representing our explicit experimental contrasts. First, we ask which genes and processes are involved in regulatory responses of hosts infected by Mytilicola by contrasting experimentally infected hosts with control hosts, representing the infection contrast. Second, we identify host genes that are differentially regulated in infections with locally co-evolved parasites compared to infections with parasites that lack a recent shared co-evolutionary history, representing the sympatry-allopatry contrast. Lastly, our dual-species transcriptomics approach allows us to use the same sympatry-allopatry contrast to identify differentially regulated genes in the parasite, which can be matched to the host sympatry-allopatry contrast, and couple phenotypic evolutionary trajectories to co-evolved gene regulatory processes in antagonistically interacting species.

| Experimental procedures
The details of the underlying cross-infection experiment are described in Feis et al. (2016). Briefly, mussels M. edulis and parasitic copepods M. intestinalis were collected at two locations, representing the two invasion pathways of the parasite. These were at Vlakte van Kerken on Texel, The Netherlands, representing the southwestern pathway and in K€ onigshafen on Sylt, Germany, representing the northern pathway. The mussels used in the experiment were treated against any previous Mytilicola spp. infections according to the method of Blateau, Le Coguic, Mailhe, and Grizel (1992) and were left to recover after the anti-Mytilicola treatment for at least 2 weeks before the experiment started. The copepodites used in the experiment to infect mussels were hatched from egg sacks of gravid M. intestinalis from the two locations. Mussels were infected with 24 copepodites from one or two mothers. Therefore, copepodites within one infection were related; however, we do not know whether multiple males sired the offspring from one female. Assuming that relatedness within the population will be 0, we can be sure that relatedness was higher within a single infection than between infections. In a full-factorial design, we generated different experimental contrasts, that is, infected vs. uninfected control mussels, sympatric (Texel host 9 Texel parasite, Sylt host 9 Sylt parasite) and allopatric (Texel host 9 Sylt parasite and vice versa) host and parasite combinations. In the control treatment, we did not expose mussels to any copepodites. Each mussel was kept in a temperaturecontrolled room at 18°C in its own bottle that was supplied with its own flow-through inlet distributing filtered seawater equally (salinity of the water input was 29-30 psu). Bottles with mussels were randomly relocated within and between containers twice per week.
Mussels were dissected under a stereomicroscope 80 days after the experiment had started. The intestinal tissue of each mussel, which is in direct contact with the parasite, and M. intestinalis specimens were put in RNAlater (Qiagen, Hilden, Germany) separately and were frozen at À80°C until further processing.
2.2 | RNA extraction, library preparation and sequencing RNA extractions of mussels and parasites were done with the RNeasy Plant Mini Kit (QIAGEN, Hilden, Germany). Since extraction of single female parasites did not yield enough RNA for sequencing, we pooled three female M. intestinalis per sample. We had three biological replicates for each treatment combination (see Table 1), except for the Texel host 9 Texel parasite combination, for which we had two biological replicates as not enough material was present for pooling a third replicate. Total RNA concentrations were mea-    UniRef, only retaining hits with an e-value cut-off of <1 9 10 À5 .

| Mytilicola intestinalis: assembly and annotation
We used TransDecoder to find the most probable ORFs (Haas et al., 2013) and InterProScan (Jones et al., 2014) to gain more insight into the potential function of differentially expressed M. intestinalis contigs without a BLAST annotation. After removing contigs that had a bivalve (host) or phytoplankton (food) annotation, 44,616 contigs remained (see Table 2 for assembly statistics). The quality of the parasite transcriptome backbone was assessed through BUSCO (Simao, Waterhouse, Ioannidis, Kriventseva, & Zdobnov, 2015) (Table S1).
The same gene in different organisms may have different functions, but belong to the same gene regulatory network (e.g., Hinman & Davidson, 2007). Understanding the potential pitfalls of inferring function from annotations based on sequence similarity without further experimental evidence on each specific gene for mussel and M. intestinalis, we will assume homologous function for the differentially expressed genes we report.

| Differential expression analyses
Differential expression analyses and GO enrichment analysis were performed in the R statistical environment version 3.4.1 (R Core Team 2017), after the reads from each library were mapped back to the respective transcriptome backbone, using the RNA seq tool in CLC Genomics Workbench 8.5.1 (CLCbio, Denmark), with the length and similarity fraction parameters set to 0.5 and 0.8, respectively.
The resulting count data were further analysed using the DESeq2 package (Love, Huber, & Anders, 2014) from BioConductor (Huber et al., 2015), where we selected the 18 transcriptomes of mussels and 11 of parasites that belong to the experimental groups described in this article.
To address the different questions for mussel hosts, we first analysed the complete mussel RNAseq data set to examine the difference between all infected vs. all control mussels and then focused on a subset including only the infected mussels to differentiate between sympatric and allopatric infections. For the parasite, we only contrasted sympatric and allopatric infections. The DESeq experiment design was full rank for all data sets. To include only meaningfully expressed contigs, we applied a mapping cut-off filter by excluding contigs with fewer than five hits per sample on average. Main effects and interaction effects were tested by Wald tests. For each main and interaction effect, we used the fdrtool for false discovery rate correction with the Benjamin-Hochberg method (Strimmer, 2008a(Strimmer, , 2008b. We consider contigs differentially expressed if the adjusted p-value (p adj ) < .05 and the absolute log2fold change (LFC) ≥1.5.
To gain further insights into the functions of differentially expressed contigs for mussels, we performed GO enrichment analyses for biological processes of the main and interaction terms based on the GO annotations from Trinotate using the R-package topGO (Alexa, Rahnenf€ uhrer, & Lengauer, 2006). A similar approach for M. intestinalis was not successful because of the low number of differentially expressed contigs with GO annotations. Due to their functional relevance in the interaction with M. intestinalis, we examined the differentially expressed mussel genes with the GO annotations (including offspring) in the biological processes of "response to stimulus" (GO:0050896), "immune system processes" (GO:0002376) and "death" (GO:0016265) closer, and in addition, we manually curated this DEG set by including contigs with the same or similar protein name that lacked a GO annotation.

| Mussel and Mytilicola transcriptome backbone quality
We obtained a total of 348 million cleaned reads for the mussel M. edulis and 108 million cleaned reads for the parasite M. intestinalis, resulting in assemblies of 163,270 and 45,645 contigs, respectively (Table 2). Applying our mapping cut-off filter resulted in 63,577 contigs for the infection vs. control differential expression analyses and in 64,759 mussel and 19,546 parasite contigs for the sympatry vs. allopatry contrast.
3.2 | Differential expression analysis in mussels: infection vs. control We found a total of 471 differentially expressed contigs (p adj < .05 and LFC ≥ 1.5, 233 had an annotation, Table S2): 226 in the infection contrast, 176 in the host source contrast and 184 for the interaction infection by host source (Figure 1a,   (Neubauer & Pitson, 2013), was also downregulated in M. intestinalis-infested mussels.
3.3 | Differential expression analysis in mussels: sympatry vs. allopatry In comparison with the infection contrast, we found more than twice as many differentially expressed contigs (946 at p adj < .05 and LFC ≥ 1.5, of which 262 were annotated, Table S3). The majority (

| Differential expression analysis in Mytilicola parasites
We found 53 differentially expressed contigs (p adj < .05 and LFC ≥ 1.5). In contrast to the host sympatry contrast, where most genes were differentially regulated depending on parasite source, in the parasite transcriptomes most differentially regulated genes were found between host sources (30 genes, 57%), while only 13 (25%) originated from differences between the parasites themselves. The interaction term representing regulation patterns for specific host 9 parasite combinations yielded 16 genes (30%). About half (26 genes) of the 53 differentially expressed genes had a significant BLAST or InterProScan hit, leaving 27 genes without putative functional annotation ( Figure 3, Table S4).

| 1511
Differences between the host sources were mainly due to downregulation of genes in the parasites infecting Sylt hosts (Figure 3).
Where homologies could be determined, these mainly included genes involved in development of the cuticle (four cuticle proteins, one chitin-binding domain, LFC range À1.8 to À4.4; Figure S1), which could indicate that the external physical barrier of the parasite suffered less damage in the guts of Sylt mussels-a pattern matching the downregulation of Chitotriosidase 1 in Sylt mussels (Figures 2, 4 and S2). Demann, unpublished data). In this way, we probably captured the chronic effects of the parasite rather than the acute ones, but such chronic effects can also be expected for the natural situation. At the functional level, we found that with 13 out of 37 differentially expressed annotated immune genes ROS pathways were important in the interaction between M. intestinalis and its mussel host.
Because the ROS response is a fairly general and unspecific response to stress and diseases (Naviaux, 2012), there have to be gene regulatory mechanisms that create the observed specificity reflecting local adaptation by coupling of the matching host and parasite genes.
Here, the higher expression of ROS-producing genes in sympatric hosts was coupled to ROS antagonist production in the corresponding parasites demonstrating this specificity (Figures 4, S1, and S2).
This link between host and parasite transcriptomes was also reflected in other functionally related interactions. For example, chitin digestion of the hosts was coupled to cuticle protein expression in the parasite and might be a core mechanism underlying the resistance phenotype differences previously observed along the invasion fronts (Feis et al., 2016). Such coupled gene expression pattern not only demonstrates the usefulness of dual transcriptomics approaches (Greenwood et al., 2016), but, more importantly, it shows that regulation of functionally important differences evolved rapidly along the two invasion fronts, supporting the usefulness of invasive species' natural history to understand host-parasite co-evolution in the wild.

M. intestinalis
Mussels that were experimentally infected with M. intestinalis parasites regulated the expression of genes with functions in cell cycling and cell division and associated with wound healing (Figure 1a).
Mytilicola intestinalis is known to obstruct the mussel's intestine and perforate the intestinal wall with its appendages and hooks (Figueras, Jardon, & Cladras, 1991;Moore, Lowe, & Gee, 1977;Robledo, Santar em, & Figueras, 1994), and this wounding is reflected in the transcriptome of the host. Increased cell division suggests that wounds in the gut wall are being closed and the intestinal epithelium is being replaced after disruption to guarantee proper function of the gut. In addition, downregulation of anticlotting agents (e.g., a homolog of fibrinolytic enzyme isozyme C) in M. intestinalis-infested mussels indicates an increased rate of blood clotting in infected mussels. Wounding and repair of physical damage therefore represent the major regulators of host gene expression.
Apart from the primary response to physical damage, host transcriptomes in the sympatry-allopatry contrast were also enriched in genes linked to an immune response. Similar to other studies investigating immune responses in mussels (e.g., Liu et al., 2014), we found the upregulation of genes involved in the production of ROS in M. intestinalis-infested mussels. This was accompanied by upregulation of genes providing antioxidants, likely to limit the negative effects of ROS on the host's own tissues. This mechanism of defence is indeed general and does not only combat pathogen infections, but also plays a role in defence against infection by macroparasites (for instance, see Coustau et al., 2015). However, we additionally found significant differences in ROS response between sympatric and allopatric combinations. Therefore, although the ROS response may be a general reaction to infections, there is likely a certain degree of specificity in the regulation of these processes leading to local adaptation.
We further found a differential regulation of the homolog maternal embryonic leucine zipper kinase, an enzyme that can be involved in apoptotic processes (Jiang & Zhang, 2013). Apoptosis may in part be a consequence of ROS production (Mat es & S anchez-Jim enez, 2000). Apoptotic processes play an important role in the immune reaction to foreign substances and were predicted to be important as first line of defence in the digestive tissues of mussels to macroparasites such as M. intestinalis (Romero, Novoa, & Figueras, 2015). Consequently, apart from the direct response to wounding, the upregulation of genes associated with cell division processes may also be linked to apoptosis by replacing apoptotic cells.
Other immune genes that were differentially expressed in infected mussels (Figure 2a 071,41,162,50,761,5,027,106,242 and 14,158 for ROS;seven contigs numbered 1,728,38,961,14,694,116,190,90,969,61,914,105,058 for ROS antagonists) and parasites (two contigs numbered 8,928 and 8,821). The middle shows attack on cuticle with the read counts of Chitotriosidase 1 of the hosts (both contigs) and cuticle proteins of the parasite (five contigs numbered 27,351, 14,629, 9,376, 18,979 and 9,312). The lower panel shows proteolytic activity of the parasite (read counts of tick legumain) as a potential virulence factor. Interaction plots of the individual differentially regulated genes can be found in Figures S1 and S2. Note different scales on y-axes used to illustrate the relative differences within rather than the absolute differences between species [Colour figure can be viewed at wileyonlinelibrary.com] for parasite feeding activity on host tissue. Legumain is expressed in the midgut in blood digesting ticks, where it serves as a haemoglobinase, processing haemoglobin and activating other proteases (Dall & Brandstetter, 2016;Sojka, Francischetti, Calvo, & Kotsyfakis, 2011).
In other invertebrate groups (Fuzita et al., 2015 and references therein), but also in the protochordate Branchiostoma belcheri (Teng, Wada, & Zhang, 2009), legumain plays a role in food digestion: it breaks down macromolecules in the acidic gut environment. Thus, the higher legumain expression levels in sympatric combinations could be a direct indicator of lower host condition (Feis et al., 2016).
Furthermore, other genes potentially involved in feeding (e.g., MAM and a low-density lipoprotein [LDL]-receptor class A [LDLRA] domain-containing protein) showed similar host-parasite combination-specific expression patterns, indicating once more a stronger exploitation of sympatric hosts. Arthropods have to acquire cholesterol from their diet as they are incapable of its de novo biosynthesis (Hassett, 2004;Perner et al., 2016), LDL is a rich source of cholesterol (Perner et al., 2016). Consequently, LDL receptors are also upregulated in scorpion digestion (Fuzita et al., 2015). The increased transcription of these genes in sympatric combinations might therefore indicate increased digestion of host tissue by Mytilicola, leading to higher energy drain from the host, which could further contribute to lower body condition in sympatric combinations (Feis et al., 2016).
Other genes with host-parasite combination-specific expression patterns in the sympatry contrast were involved in the parasite's immune response or defence against the host's immune response.
The upregulation of a C-type lectin homolog in sympatric combination ( Figure S1) suggests an increased immune activity of the parasite itself. C-type lectins are widespread pattern recognition molecules that are involved in carbohydrate recognition often associated with bacterial lipopolysaccharide in crustaceans (Kawabata & Iwanaga, 1999), suggesting that the parasite's immunological surveillance may be fine-tuned towards bacteria encountered in sympatric host guts. Since immune resources of the host are diverted to combat the parasite, less immune resources might be available to maintain bacterial homeostasis in the gut. Disturbed host microbiota can contain higher number of secondary opportunists (Lokmer & Wegner, 2015) that could potentially also harm the parasite. This effect might be particularly strong in the already weakened sympatric hosts therefore calling for an increased level of immunological surveillance on the parasite's side as well.
A strong relative increase in host ROS-producer transcription after infection (Figure 4) was one of the main defence mechanisms.
Here, the pattern of host ROS-producer transcription supports the hypothesis of a higher degree of self-damage in sympatric combinations (Feis et al., 2016). Host expression patterns of ROS antagonists, on the other hand, displayed opposite patterns, with especially high levels of host ROS antagonist expression in Texel hosts 9 Texel parasites combinations ( Figure 4). However, the parasite expression profiles revealed that the parasites in Texel host 9 Texel parasite combinations hardly express any ROS antagonists. Assuming that transcription will positively correlate with ROS production this may indicate that the total amount of ROS scavengers might be lower and more ROS is freely available. A qualitatively similar pattern could be derived for the combination of Sylt hosts 9 Sylt parasites, where the increase in host ROS production expression was not associated with a remarkable increase in either host or parasite ROS antagonist expression, which might also result in a higher relative amount of freely available ROS in the mussel gut. The high level of ROS antagonist transcription in Sylt parasites in general might therefore explain the small reduction in host condition following infection with parasites from this location (Feis et al., 2016). The link between ROS production and ROS antagonist expression couples the expression patterns of parasites to those of the host and may be associated with to the phenotype of lowered body condition potentially resulting from oxidative stress. Further studies should therefore measure ROS production and oxidative stress directly to test this link on a functional level (Griendling et al., 2016).
Apart from the expression patterns of ROS pathways in the interaction between host and parasite, the differential expression of Chitotriosidase 1 homologs (Chit1) in mussels and cuticle proteins in M. intestinalis could be another example of a direct coupling of gene expression between hosts and parasites. As an essential component of the innate immunity against chitin-coated pathogens, Chit1 also regulates inflammatory processes (Elmonem, van den Heuvel, & Levtchenko, 2016;Van Eijk et al., 2005). However, due to its digestive function in mussels (Birkbeck & McHenery, 1984;Lesser & Macmanes, 2016), oysters (Yang et al., 2015) and insectivorous scorpions (Fuzita et al., 2015), it is also another example of the dual roles of genes in digestion and immunity.
Interestingly, the upregulation of Chit1 corresponded to the upregulation of a number of cuticle proteins in the parasites infecting Texel hosts (Figures 4, S1, and S2), suggesting that direct targeting of the parasite's chitin cuticle might have led to the parasite's response of strengthening the weakened exoskeleton with increased cuticle production. Since we only captured successful infections of parasites that were able to respond to the host attack, higher resistance in Texel mussels and higher infectivity of Texel parasites could be explained by successful expulsion of parasites that failed to counter the host chitinase immune effector by upregulation of their cuticle production. Hence, the coupling of host Chit1 and parasite cuticle proteins might characterize another direct interaction and help to explain how different evolutionary trajectories were realized along both invasion fronts (Feis et al., 2016).

| CONCLUSION
In this study, our dual-species transcriptomics approach elucidated some molecular underpinnings of the young co-evolutionary relationship between mussels M. edulis and their intestinal macroparasites M. intestinalis. The mechanistic interaction pathways described by several candidate genes partly matched the different