Parasite manipulation of host phenotypes inferred from transcriptional analyses in a trematode‐amphipod system

Manipulation of host phenotypes by parasites is hypothesized to be an adaptive strategy enhancing parasite transmission across hosts and generations. Characterizing the molecular mechanisms of manipulation is important to advance our understanding of host–parasite coevolution. The trematode (Levinseniella byrdi) is known to alter the colour and behaviour of its amphipod host (Orchestia grillus) presumably increasing predation of amphipods which enhances trematode transmission through its life cycle. We sampled 24 infected and 24 uninfected amphipods from a salt marsh in Massachusetts to perform differential gene expression analysis. In addition, we constructed novel genomic tools for O. grillus including a de novo genome and transcriptome. We discovered that trematode infection results in upregulation of amphipod transcripts associated with pigmentation and detection of external stimuli, and downregulation of multiple amphipod transcripts implicated in invertebrate immune responses, such as vacuolar ATPase genes. We hypothesize that suppression of immune genes and the altered expression of genes associated with coloration and behaviour may allow the trematode to persist in the amphipod and engage in further biochemical manipulation that promotes transmission. The genomic tools and transcriptomic analyses reported provide new opportunities to discover how parasites alter diverse pathways underlying host phenotypic changes in natural populations.


| INTRODUC TI ON
Parasites often form tight ecological associations with their hosts, leading to evolutionary arms races (Poulin, 2007).From the perspective of the parasite, these complex co-evolutionary dynamics often promote parasitic strategies to avoid host detection or, in some cases, modify the host's behaviour (Poulin, 2013).This is prominent for parasites that depend on multiple hosts to complete their life cycle.In these cases, parasites may modify the behaviour of intermediary hosts to ensure successful trophic transmission (Moore, 2002).These putative manipulations range from physiological alterations to behavioural changes and can have significant impacts on the community structure of their hosts (Hasik et al., 2023;Wood et al., 2007).For example, modification of host phenotypes by parasites is often documented among helminths with complex life cycles.In these contexts, parasite-induced changes in host phenotype often increase the vulnerability of intermediate hosts to predation by definitive hosts (Thomas et al., 2005).The term manipulation invokes a change in host phenotype that is induced by parasite-encoded gene products and has fitness benefits for the parasite by facilitating transmission.These adaptive scenarios have been attributed to an extended phenotype of parasite evolution (Dawkins, 1982;Thomas et al., 2012).A classic example of this phenomenon is the trematode Dicrocoelium dendriticum, which causes infected ants to ascend vegetation for probable herbivory by sheep (Carney, 1969;Moore, 2002).In other instances, multiple host phenotypes may be manipulated.For example, the avian trematode Leucochloridium induces changes in snail crypsis, such that infected hosts display an altered coloration and pulsating antennae, which mimics a caterpillar (Moore, 2002).Similarly, ants infected by an avian nematode develop berry-mimicking abdominal gasters and are driven to conspicuously perch in view of frugivorous birds (Yanoviak et al., 2008).The argument for manipulation is strengthened when multiple host phenotypes are altered in ways that enhance parasite transmission, but evidence for actual enhanced predation is needed to distinguish true manipulation from parasite-induced phenotypic alterations ('PIPA', (Fayard et al., 2020)).
While previous research on the parasitic manipulation of host behaviour has focused on its evolution and benefits to the parasite, there are relatively few studies elucidating the specific molecular mechanisms underlying this phenomenon (Carr et al., 2021;Hughes & Libersat, 2019;Ponton et al., 2006;Stoldt et al., 2021;Thomas et al., 2005;Will et al., 2020).This is in part due to the difficulty of studying the genetic basis of behaviour in a dynamic system.
While experimental manipulations coupled with advances in nextgeneration sequencing technologies have greatly enhanced our understanding of host-parasite evolution in model systems (Carr et al., 2021;Greenwood et al., 2016;Hebert et al., 2015;Poulin & Maure, 2015;Wijayawardena et al., 2016;Will et al., 2020), these studies are dominated by a limited collection of taxa: 34% of empirical studies involve two acanthocephalan parasite genera Polymorphus and Pomphorhynchus and 39% of studies focused on the single amphipod genus Gammarus (Poulin & Maure, 2015).As a result, knowledge of the specific molecular mechanisms underpinning host-parasite manipulation is lacking (Cezilly et al., 2013;Herbison et al., 2018).Progress in this area has been made with a marine amphipod host Echinogammarus marinus and acanthocephalan parasites that influence its photo-and geotactic behaviours (Guler & Ford, 2010).In a subsequent study, this host shows altered expression of neurological genes related to serotonin activity in response to a trematode parasite (Guler et al., 2015).These examples offer promising approaches to unravel the complex molecular interactions that may govern the evolution of host parasite coevolution.
The present study uses natural populations of a trematodeamphipod system to advance the mechanistic understanding of parasitic manipulation of host behaviour.We investigate the transcriptional basis of parasite-driven alterations to the amphipod Orchestia grillus, the second intermediate host of the avian trematode Levinseniella byrdi (Bousfield & Heard, 1986;Heard, 1968;Johnson & Heard, 2017).Notably, this system provides a novel perspective complementing the many examples of other amphipods that have been studied in the context of host-parasite dynamics (Giari et al., 2020;Guler et al., 2015;Poulin & Maure, 2015).The life cycle of L. byrdi begins when infected marsh-dwelling birds shed faeces containing thousands of developing eggs.Ova are then passively consumed by the first intermediate host, hydrobiid snails, feeding on marsh substrate.In turn, snails shed free-swimming larval cercariae capable of penetrating the exoskeleton of its second intermediate host, O. grillus.Approximately 25-30 days post-infection, cercaria encyst into large metacercarial cysts within the amphipod body cavity, and dramatically change host phenotypes (see Figure 1).First, infected amphipods change from their natural brown or grey colour to bright orange, and second, infected individuals appear to lose their negative phototaxis increasing the probability that they spend more time in exposed areas of the salt marsh than uninfected individuals (Johnson, 2022).These altered phenotypes may increase avian biochemical manipulation that promotes transmission.The genomic tools and transcriptomic analyses reported provide new opportunities to discover how parasites alter diverse pathways underlying host phenotypic changes in natural populations.

K E Y W O R D S
amphipod, differential expression, ecological genomics, host-parasite co-evolution, infection response, Orchestia grillus, parasite manipulation, population genetics, trematode predation risk and facilitate trophic transmission to definitive avian hosts (Johnson & Heard, 2017).
While the obvious phenotypic consequences of L. byrdi infection are well-described (Bousfield & Heard, 1986;Johnson, 2011Johnson, , 2022;;Johnson et al., 2009;Johnson & Heard, 2017), the molecular signals and mechanisms underpinning drastic changes in O. grillus coloration and behaviour remain unexplored.In this paper, we capitalize on the biology of the amphipod O. grillus and the trematode L. byrdi to elucidate transcriptomic responses to parasitism and its behavioural consequences.To this end, we assembled a de novo draft genome for the amphipod and sequenced full transcriptomes of 24 infected and 24 uninfected individuals collected from a single natural population.We sought to characterize which transcripts show differential abundance in response to trematode infection.Gene ontology analyses of the differentially expressed genes show a striking correspondence to the primary phenotypic changes observed in infected amphipods.Overall, our findings provide insights into transcriptomic responses of host manipulation in natural populations and identify candidate genes and pathways that likely govern this host-parasite interaction.

| Amphipod datasets and draft genome
We generated two datasets for this study.First, two whole-genome sequences were obtained from DNA extracted from a single leg of a male and a single leg of a female amphipod, both sequenced to deep coverage using Illumina paired end reads.Leg tissue is very unlikely to be infected with trematode DNA as the cysts form inside the body cavity.Thus, this assembled genome should contain little or no trematode sequences.Second, a panel of 48 individual transcriptomes was sequenced from 24 infected and 24 uninfected amphipods using paired-end Illumina RNAseq.Both these panels were generated from individuals collected at the banks of the Club Head Creek in the Plum Island Estuary, Massachusetts, USA (Figure 1a).
For the transcriptomic panel, we classified individuals as infected or uninfected based on coloration and presence/absence of cysts (Figure 1b,c).The L. byrdi cysts are relatively large (0.4 mm diameter) and easy to find in the body cavity, and infected amphipods have 1-4 cysts per individual, with one cyst per infected amphipod being the most common (Johnson & Heard, 2017).These two criteria (colour and presence of at least one cyst) simplified the classification of infection status.We chose to extract RNA from whole individual animals given that three potentially different amphipod phenotypes are altered under infected conditions (red colour change, locomotor or phototactic responses, and immune responses).Whole body RNA extractions would dilute the signal from any one of these traits, but would ensure that each of the transcriptional factors anticipated to underly the altered phenotypes would be sampled.
Since no genomic tools were available for the species, we generated a de novo draft genome using the reads from the female amphipod leg DNA (see Data S1).We assembled the draft genome using a two-step approach: a first pass done with the SparseAssembler (Ye et al., 2012), and a second "heterozygote corrected" assembly and scaffolding done with Redundans (Pryszcz & Gabaldon, 2016).
The first pass assembly produced a highly fragmented genome with 195,956 contigs (>500 bp) and an N50 of 4695 bp.The total length of genome assembled was 812,068,892 bp.Subsequent use of the Redundans pipeline reduced the number of contigs to 87,264 and increased the N50 to 17,506.The draft genome generated is available in NCBI (GenBank assembly accession: GCA_014899125.1;GenBank: VOMN00000000.1;Genome name = Ogril1).This draft genome was used as the mapping reference for all genomic and transcriptomic datasets (see additional details in Data S1).

| Predicting amphipod genes
Gene prediction was done using the ab initio approach Augustus (Stanke et al., 2006) using RNA-seq mapping evidence and a model trained on the well-annotated genome of Drosophila melanogaster to take advantage of gene ontology information for that genome.
The RNA-seq evidence was obtained from mapping the pooled RNA reads of all 24 uninfected amphipod individuals to the reference (we used a similar code to (Nunez et al., 2020)).This method predicted the existence of ~21,000 putative protein coding genes.In order to gain further evidence regarding the identity of these loci, we ran pairwise BLAST comparisons of the predicted Ogril1 amino acid sequences relative to the transcriptome of Drosophila melanogaster (assembly release 6 plus ISO1 MT; Dmel6; NCBI ID: 47).Of all the predicted structures, only ~10,000 sequences had BLAST hits with e-values <10 −10 , which were retained for further functional annotation analyses.

| Genomic and transcriptomic analysis
We first generated a panel of single nucleotide polymorphisms (SNPs) by mapping RNA reads from all 48 individuals to the reference genome.A subset of scaffolds showed unexpected levels of heterozygosity and population genetic structure that are indicative of assembly errors or linkage to sex chromosomes, and were therefore excluded from subsequent transcriptomic analyses (Data S1).A principal component analysis (PCA) built using the expression levels of all transcripts does not show any clustering patterns in the data related to infection status (Figure 2a).We conducted the RNA-seq analysis using three parallel approaches: EBSeq (Leng et al., 2013), edgeR (Robinson et al., 2010), and DESeq2 (Anders & Huber, 2010;Love et al., 2014).While each of these approaches showed statistically significant patterns of differential expression between the 24 infected and the 24 uninfected transcriptomes, the analyses differ in the number of differentially expressed (DE) genes detected (DESeq2 = 754 DE genes, EBSeq = 557 DE genes, edgeR = 1067 DE genes; see Table S1 and our Reproducible Analyses in GitHub (https://github.com/DavidRandL ab/Amphi pod-host-paras itegenomics)).In all three analyses there were more down-regulated genes than upregulated genes in infected individuals relative to uninfected (Figures 2b and S1).Notably, when the PCA is repeated using the top 100 differentially expressed genes, individuals cluster by infection status (Figure 2c).

| Functional analyses
We used Gene Ontology (GO) information available for D. melanogaster as an aid to interpreting the patterns of differentially expressed genes.The GO categories from the DE genes identified by DESeq2 and edgeR are similar, while those from EBSeq are different.We note that EBSeq prioritizes the roles of differentially spliced transcripts in the identification of DE genes, and this may be the source of different results.Nevertheless there are 610 genes shared across at least two of the three analyses that show significant differential expression (Figure 3a).The GO categories for up-and downregulated genes for the 610 shared genes are shown in Figure 3b.
Details of all three individual GO analyses are reported in Table S2 and Figure S2.GO analyses for significant DE genes without reference to the direction of expression are shown in Figure S3 and Table S2.Some categories are seen in both the directional and nondirectional GO analyses, while others are unique.We have focused on the up-and down-regulated GO analyses based on the directionality of some phenotypic effects: more orange pigment, suppression of immune responses.Because we annotated Ogril1 based on homology to the genome of D. melanogaster, all candidate genes have a corresponding gene ontology (GO) term from FlyBase (https:// flyba se.org/).We considered adding other crustacean genomes to aid in the annotation (e.g., Daphnia or Hyalella), but the goal of the functional annotation is to associate the orthologs we have identified in Orchestia with particular phenotypic effects of known genes.This effort requires a model system where mutational studies have allowed the gene to be causally associated with a particular phenotype or pathway and Drosophila offers the most complete reference of this type.
Using these FlyBase GO terms, we conducted enrichment analyses for the concordant targets.We partitioned this analysis on two groups: downregulated and upregulated targets in infected individuals (Figure 3b).For genes that are upregulated in the infected group, we observe a conspicuous enrichment of genes related to detection of external stimuli as well as pigmentation (Figure 3b in red).In the case of downregulated genes, we observe enrichment for a variety of metabolic functions including pH, transmembrane transport, and ion regulation (Figure 3b in blue).
In addition to down-regulation of individual genes involved in lipid and carbohydrate metabolism, there is a notable enrichment of genes encoding subunits of the Vacuolar-type H + -ATPase complex (V-ATPase; Figure 4; see Table S1).V-ATPase is an important proton pump in eukaryotic cells (Beyenbach & Wieczorek, 2006;Dietz et al., 2001;Harvey & Wieczorek, 1997;Nelson et al., 2000).
It is composed of two functional domains V 0 and V 1 .The V 1 domain is a peripheral complex with ATP hydrolysis capacity.It is composed of 8 subunits A-H (14-70 kDa); V 0 on the other hand is tasked with proton translocation across the membrane and it is composed of 5 subunits (a, d, c, c', and c"; 17-100 kDa) (Beyenbach & Wieczorek, 2006).Our analysis reveals a general pattern of downregulation in infected individuals of 13 amphipod genes with homology to V-ATPases.These genes correspond to 6 homologues of subunit a, 2 homologues of subunit H, and one homologue of subunits d, B, C, D, and c' (Figure 4).

| DISCUSS ION
Manipulation of host phenotypes is a common strategy used by parasites to increase their survival as well as transmission.Despite the fascinating biology encapsulated in this process, we still lack a fundamental understanding of the specific molecular mechanisms underpinning host-parasite manipulation (Cezilly et al., 2013;Herbison et al., 2018).Progress in this area is being made in certain ant-fungal systems and amphipod-trematode systems (De Bekker et al., 2015;Elya et al., 2018;Guler et al., 2015;Will et al., 2020).While each of these host-parasite systems shows evidence for parasite manipulation of host phenotypes, the mechanisms driving each interaction are likely to be distinct at the cellular and molecular level.An open question is whether there are common pathways that emerge across different host-parasite systems (Guler et al., 2015;Ponton et al., 2006;Will et al., 2020).In this paper, we investigated the transcriptomic basis of parasite-driven alterations by the avian trematode L. byrdi to its second intermediate host, the amphipod O. grillus.
We sequenced the transcriptomes of 24 infected and 24 uninfected amphipods collected from the same natural population.Our results suggest that hundreds of gene transcripts vary in abundance in response to infection, and over 600 transcripts are shared across three different analysis pipelines (DESeq2, EBSeq and edgeR).We observe F I G U R E 2 Principal components analysis (PCA) and differential expression.(a) PCA built using expression data analysed in DESeq2 for all transcripts detected in the 24 infected and 24 uninfected amphipods.There is no clear genome-wide transcriptional differentiation of infected (blue) an uninfected (red) amphipods.(b) Volcano plot of differential expression with fold change on the x-axis and FDR-corrected significance on the y-axis.The sign of the values on the x-axis is based on a contrast of (infected -uninfected), thus down-regulated transcripts with negative fold change tend to occur in infected amphipods.(c) PCA built using the top 100 differentially expressed transcripts between infected versus uninfected amphipods, showing a clear differentiation across PC1 for these transcripts.upregulation of pigmentation-related genes and genes involved in detection of external stimuli in infected individuals.These are intuitive results given the clear changes in pigmentation of amphipods as a consequence of infection, and their increased tendency to be found in exposed areas of salt marsh habitat (Figure 1b,c).In addition, we observe that infected individuals show many downregulated genes involved in carbohydrate metabolism, transmembrane transport, sterol transport, ion transport, and pH regulation (e.g., vacuolar-ATPases).While we acknowledge that predicting organismal behaviour from transcriptional profiles spans multiple important cause-effect relationships, the phenotype-genotype associations identified here highlight potential gene targets driving host manipulation in the wild.These motivate testable mechanistic hypotheses about the genomic underpinnings of host-parasite dynamics.

| Do upregulated genes underlie amphipod phenotypes that enhance trematode transmission?
Our data show evidence for the upregulation of pigmentation and photoreception genes in infected, relative to uninfected, individuals.Notably, the bright red, infected amphipods not only spend more time in open exposed areas in the marsh substrate, they are sluggish and do not scurry under vegetation when disturbed ( (Johnson, 2022) and personal observation).Although field observations do not indicate that the amphipods become positively phototactic upon infection, their altered behaviour results in them being found in open areas of the salt marsh substrate more often than uninfected individuals.This may be due to their movements being more sluggish or negative phototaxis being diminished, or both.We note that our transcriptional data were obtained from infected and uninfected amphipods that were housed together in 'common garden' containers with moist paper towels and a 12:12 light:dark cycle for one month after collection from the field, before RNA was extracted.It is possible that the upregulation of photoreceptor genes in infected individuals could be due to more exposure to light in the field, but it seems unlikely that this expression pattern would persist for one month of laboratory acclimation in the same containers and light environment as the uninfected amphipods.We do not believe the differential expression of pigment and photoreception genes is L. byrdi in the amphipod (Johnson, 2022;Johnson & Heard, 2017).
These kinds of phenotypic changes, particularly in intermediate hosts in multi-host systems, are interpreted as an adaptive induction of the host phenotypes to increase visibility and enhance predation, facilitating transmission of the parasite to the next host in the life cycle (Poulin, 2013).From the perspective of the trematode L. byrdi in our study system, transmission from the amphipod to the final bird host requires that birds eat amphipods.As a result, phenotypes that enhance predation of the amphipods by birds are likely to be targeted by natural selection to increase trematode transmission and fitness.Positive demonstration of enhanced predation is a key criterion for the adaptive nature of the parasitism (Poulin & Maure, 2015), and direct evidence for enhanced predation is still lacking in this O. grillus -L.byrdi system, but observational and experimental studies do show that L. byrdi infections reduce locomotor activity of amphipods resulting in them spending more time in exposed patches of the salt marsh (Johnson, 2022;Johnson et al., 2009;Johnson & Heard, 2017).This is likely to increase opportunities for predation.
The genes identified in the significant Gene Ontology categories "eye pigmentation" and "pigment biosynthetic process" are homologues of Drosophila genes Culd, Pu, v, Hn, Rab32 (see Table S2 for DESeq2) that alter fly eye or body colour or phototransduction (see Flybase (www.flybase.org)).These are strong candidates for mechanisms related to increased orange pigmentation in infected amphipods as per the gene descriptions in Flybase, summarized here.Vermillion (v) is a classic eye colour mutant in flies also affecting body colour (Phillips & Forrest, 1980).Culd encodes a photoreceptor-cell enriched transmembrane protein that localizes to endocytic vesicles and is required for endocytic trafficking of the products of ninaE ('neither inactivation nor afterpotential E', the rhodopsin expressed in the largest class of photoreceptors) and trpl ('transient receptor potential-like' encodes a plasma membrane cation channel that is highly eye-enriched and contributes to the electrical response to light in photoreceptors (https://flyba se.org/repor ts/FBgn0 002940)).Culd mutants affect sensory detection in photoreceptors (Xu & Wang, 2016).Punch (Pu) encodes a GTP cyclohydrolase and mutants show defective body and eye pigmentation (Mackay & O'donnell, 1983).Henna (Hn) encodes a tryptophan phenylalanine hydroxylase and is a dual function enzyme also involved in pteridine synthesis affecting eye and body pigments (Ferre et al., 1985), as well as locomotor function.Rab32 encodes a small GTPase that contributes to vesicle trafficking regulation and is involved in eye development, autophagy, and lipid storage (Ma et al., 2004).
The significant Gene Ontology category "detection of external stimulus" and "detection of abiotic stimulus" includes homologues of Drosophila genes Arr1, Pkd2, trp, cac, inaE, and nompC (see Table S2), which are involved in photopigment and phototransduction.The annotation of these genes, as reported in Flybase (www.flyba se.org) and summarized here, also provides strong candidates for mechanisms related to the amphipod behaviour of altered phototaxis.Arrestin 1 (Arr1) encodes a protein involved in rhodopsin inactivation that contributes to photoreceptor maintenance (Dolph et al., 1993).Polycystic kidney disease 2 (Pkd2) encodes a calcium permeable cation channel of the transient receptor potential (TRP) family.Its roles include sperm motility, smooth muscle contraction, larval feeding behaviour, detection of mechanical stimuli, and cold nociception.Transient receptor potential (trp) encodes another nonselective plasma membrane cation channel that, like trpl, is highly eye-enriched and contributes to the electrical response to light in photoreceptors (Minke & Parnas, 2006).Cacophony (cac) encodes the primary structural subunit of a voltage-gated calcium channel, which is located at presynaptic active zones and functions in evoked neurotransmitter release at neuromuscular synapses.It contributes to male courtship behaviour and a wide range of neurophysiological processes (Smith et al., 1998).The gene 'Inactivation no afterpotential E' (inaE) encodes a diacylglycerol lipase involved in phototransduction (Leung et al., 2012).'No mechanoreceptor potential C' (nompC) encodes a pore-forming subunit for a mechanosensitive non-selective cation channel.It belongs to the TRP channel family and is expressed in peripheral sensory neurons.It senses gentle touch and regulates locomotion in larval body wall neurons (Yan et al., 2013).These annotations from Flybase and original references provide remarkable concordance with the physiological and behavioural processes related to the most overt phenotypic changes in infected amphipods (Johnson, 2022).
It is strongly suggestive that these genes for pigmentation are upregulated in infected amphipods, consistent in direction with increased pigmentation.For the behaviour phenotypes, however, increased expression of genes affecting 'detection of stimuli' may suggest 'more response' to stimuli, but could involve genes that reduce negative phototaxis leaving infected amphipods less likely to move away from exposed areas, a behaviour observed in both field and laboratory conditions (Johnson, 2022).We acknowledge important caveats regarding causality of the behavioural traits, as we have not established if the transcriptional changes are the drivers of the behaviours, or if the behaviours mediated by infection status lead to altered gene expression.Moreover, the direction of expression (up-vs.down-regulation of specific genes) cannot be equated with 'increased' or 'decreased' activities of certain behaviours.Studies in the amphipod Echinogammarus marinus show that infection with trematode parasites alters photo-and geotaxis of the amphipods in a manner that would increase predation, and involves serotonergic activity (Guler & Ford, 2010).Notably, the transcriptional analyses of neurological genes associated with serotonin function in Echinogammarus marinus show that while individual genes are altered in expression as a result of trematode infection, the direction of altered expression varies: five genes showed increases and five showed decreases (Guler et al., 2015).The specific neurobiological bases of the genes showing altered expression in our O. grillus amphipod study need further investigation, but the genes identified through homology and annotation in Drosophila show striking correspondence to the most evident phenotypic changes in infected amphipods.It is further noted that pigments that affect coloration can be critical in neurological processes (Yamaguchi & Hearing, 2014) and the detection of light (Lebhardt & Desplan, 2017), so the GO categories 'pigmentation' and 'detection of stimuli' can be functionally related.

| Do down-regulated genes underlie manipulation of amphipod physiology to enhance trematode fitness?
Gene ontology analyses show strong down regulation of genes involved in carbohydrate metabolism, transmembrane transport, sterol transport, ion transport and pH regulation (e.g., vacuolar-ATPases).A possible explanation for this pattern is that decreased expression of genes involved in energy metabolism could result in lower rates of catabolism in the host amphipods, allowing more resources to be available for trematodes to engage in anabolic processes for growth and encystation in their intermediate hosts.
Genes related to amylase, fumarase, B-galactosidase and chitinases are all down regulated; all are familiar genes encoding catabolic processes of various carbohydrates, including the chitinous exoskeleton of arthropods.Moulting, a process that involves chitin metabolism, has been associated with host defence against parasite infection in Daphnia (Duneau & Ebert, 2012).Down-regulation of chitinase(s) could reflect a defence mechanism of the host to avoid further infection, or an induced response by the trematode to facilitate persistence.Resolving these effects would require a more fine-grained temporal analysis of the infection cycle than we have provided here.
A particularly notable down regulated gene is Adenosine deaminase-related growth factor A (Adgf-A), which encodes an extracellular deaminase that regulates the level of extracellular adenosine by converting adenosine to inosine.This activity is important for regulation of systemic metabolism, especially during stress and infection (Bajgar & Dolezal, 2018).Down regulation of this gene could mask the physiological evidence for an infected state, allowing the trematode to remain less detected by the amphipod innate immune or infection-surveillance pathways.It is notable that infected amphipods are distinctly less active than uninfected amphipods, appearing almost sluggish.This difference may be associated with the systemic down-regulated of metabolic genes with infection.However, it remains to be determined whether this is an adaptive component of the trematode's biology to ensure transmission, rather than a nonadaptive consequence of simply being infected.
It remains possible that the transcriptional basis of particular parasite induced phenotypic alterations ('PIPA' sensu (Fayard et al., 2020)) may involve combinations of up-and down-regulated genes.We performed GO analyses of all DE genes without reference to the direction of expression (Figure S3).Non-directional GO enrichment analyses recover many of the categories in the up-and down-regulated genes, but the directional analyses reveal that certain GO terms are only found in up-or down-regulated analyses even if they persist in the non-directional analyses (compare Figure 2 with Figures S2 and S3).It may be that the impact of directional DE genes are 'strong enough' to remain significant in non-directional analyses, or the non-directional analyses capture up-and down-regulated genes that contribute to the same biological pathway.But the categories associated with the most evident amphipod phenotypes are more apparent in the directional analyses.

| V-ATPase downregulation: An immunosuppression mechanism?
The general down regulation of vacuolar-ATPase subunits in infected individuals provides another line of evidence that trematodes induce some form of immunosuppression in infected individuals (Figure 4).V-ATPases are a family of ATP-dependent proton pumps which display high levels of evolutionary conservation in eukaryotes (Beyenbach & Wieczorek, 2006;Dietz et al., 2001;Harvey & Wieczorek, 1997;Nelson et al., 2000).As proton pumps, these multi-subunit structures are known to play a variety of important functions in cells including responses to infection and cancer progression (Hackam et al., 1998;Natale & Mccullough, 1998;Pamarthy et al., 2018).For example, when mice are challenged with infection by the helminth Fasciola hepatica, the parasites secrete a 'helminth defense molecule' (HDM) that reduces the efficacy of the immune response via inhibition of lysosomal Vacuolar-ATPase (Robinson et al., 2011(Robinson et al., , 2012)).This can contribute to the persistence of the pathogen in the host.Comparable evidence has shown that V-ATPase also plays a role in mediating parasite invasion in mosquitoes (Huang et al., 2006).Overall, these studies suggest that V-ATPase pH regulation function is pivotal for a cell's capacity to mount parasitic immune responses in both vertebrate and invertebrate systems.As such, we hypothesize that the downregulation of V-ATPases in infected amphipods is a mechanism used by L. byrdi trematodes to suppress the host amphipod's immune response.In the mouse-helminth model, a HDM was identified that is secreted by the helminth, proteolytically processed by the host lysosome, and prevents vacuolar acidification by inhibiting the activity of V-ATPase (Robinson et al., 2012).This mammalian system involves components of both the innate and adaptive immune systems, while our amphipod-trematode system involves only innate immune pathways.The parallel findings in our amphipod-trematode system and the mouse-helminth system may point to conserved mechanisms by which pathogens can avoid detection by suppressing host immune functions.We have not identified a comparable 'trematode defense molecule' that directly inhibits the V-ATPase, but our results suggest that efforts with paired RNAseq and proteomic analyses in the amphipod and trematode might uncover similar targets as an outcome of trematode infection.Experimental manipulation combined with next-generation sequencing approaches has been done effectively in the ant-fungal systems (Elya et al., 2018;Will et al., 2020).The difficulty of culturing the trematode in our salt marsh, three-host system makes experimental manipulations a challenge, but the new findings do point to candidate targets that may be identified with other 'omics approaches.
Although we observed a general pattern of suppression of V-ATPase in O. grillus infected with L. byrdi, we should point out possible caveats.At least three of the infected O. grillus individuals have gene expression profiles of V-ATPase subunits that are qualitatively more similar to uninfected individuals (Figure 4).
One possibility is that the infected individuals we sampled were at different stages of disease progression, and that the hypothesized effect of L. byrdi on host immunosuppression occurs at a later stage than colour phenotype manipulation.Alternatively, this result may reflect idiosyncratic infection outcomes, where hosts may undergo colour shift and yet resist immunosuppression.Both of these possibilities -staggered effects of parasites on host gene expression throughout disease progression (Will et al., 2020), and population genetic variability in infection outcomes (Poulin, 2013;Wang et al., 2020) --are documented in other host-parasite systems characterized by host manipulation.
Parasite load may also contribute to variable gene expression and infection outcomes (Amorim et al., 2019).While we were able to bioinformatically extract L. byrdi reads from host cyst samples, this was not adequate to accurately estimate their abundance.
All infected amphipods had at least one L. byrdi cyst, and a few had several, but the range of cyst number was low, and a single cyst was the most common condition.Moreover, it is difficult to infer the physiological impact of single vs. multiple cysts as they differed in size.But it only takes one metacercaria to induce the primary phenotypic effects observed in infected amphipods (Johnson, 2022).
At the outset of this study we expected to identify many upregulated immune genes in the infected amphipods.Multiple studies in Drosophila have shown that immune responses vary in response to the type of parasite (Hoffmann, 2003;Roxstrom-Lindquist et al., 2004).For example, when Drosophila are infected with a protozoan parasite, immune responses normally include upregulation of GNBP genes, lysozyme genes, and lectin genes among others (Roxstrom-Lindquist et al., 2004).We did not observe upregulation of these genes among our candidates but it is possible that such genes are either missing from our assembly or are misannotated.An alternative hypothesis is related to our V-ATPase findings: L. byrdi suppression of V-ATPase may prevent the initiation of a canonical immune response in the amphipod host.We classified individuals as 'infected' based on bright coloration and the presence of an encysted trematode inside an amphipod.Individuals were placed in a common lab environment for ~one month prior to RNA extraction and these individuals could have matured past an early-infectionstage immune response after trematode larva(e) entered the amphipod.This may have reduced our ability to detect an initial immune gene expression response.

| Host-parasite transcriptomics and reverse ecology
In this paper we have provided evidence for the transcriptomic underpinnings of host manipulation in a natural population of amphipods.These findings showcase the power of "reverse ecology" (Ellison et al., 2011;Levy & Borenstein, 2012;Li et al., 2008), an ecological version of 'bulk-segregant analysis' which leverages in-depth knowledge of distinct phenotypic characteristics of an organism's natural history to learn the molecular drivers of ecologically important phenotypes.Reverse ecology is a powerful approach in the era of genomics since it envisions non-model systems as "natural laboratories".However, there are clear drawbacks-and room for improvement-which must be acknowledged.As previously mentioned, all our analyses depend on the quality of our genomic tools, and like most natural systems this remains a fundamental bottleneck in ecological genomics research programs.
We further acknowledge that the inferences from DESeq2, EBSeq and edgeR are not completely concordant, with EBseq results for DE genes and GO terms being distinct from DESeq2 and edgeR.
As noted above, EBSeq highlights different splice variants as an informative component of 'differential expression', which may account for the distinct results.Nevertheless, there are a core set of ~600 genes shared by these three methods that lend support for the main inferences we discuss.Despite these caveats, the additional insights provided by RNAseq analyses lend support for the notion that L. byrdi parasites are manipulating their O. grillus hosts to aid transmission.The evidence presented strongly suggests adaptive manipulation because: (A) we see at least three phenotypic changes: (1) body colour, (2) neutralized phototaxis, and (3) reduced fleeing from potential predators (Johnson, 2022;Johnson et al., 2009;Johnson & Heard, 2017); (B) we provide novel support for a parasite-induced suppression of the host immune system without incapacitating the amphipods with an overwhelming parasite burden; and (C) 99% of the amphipods found in the exposed habitats at low tide -and therefore susceptible to bird predation -are infected with L. byrdi (Johnson et al., 2009).
The suite of phenotypes and transcriptional profiles reported here -spanning alterations affecting pigmentation and crypsis, neurological functions and immune effects -provide an opportunity to test the neuropsychoimmune hypothesis of host manipulation (Fayard et al., 2020).That each of these broad phenotypic categories can be implicated in parasite persistence or increased risk of predation, strengthens the support for manipulation.A finer time course of infection across developmental stages and tissue types would be highly informative in assessing the integration across the pigmentneurological-immune axes.
The candidate genes uncovered in this study may help extend the breadth of inference in this reverse ecology approach.As the sequence reads for both our genomic and transcriptomic data are freely available, we hope these resources will enable others to provide new insights into the co-evolutionary signatures of hostparasite interactions.

| Nucleic acid extractions and sequencing
DNA was extracted separately from a single leg of a male and a female O. grillus using the Qiagen DNeasy Blood and Tissue Kit and tissue disruption, via TissueLyser.To extract cyst DNA, four cysts were dissected out from one infected amphipod and pooled together for extraction.Samples were quantified by Nanodrop and diluted to equimolar concentrations prior to sequencing.O. grillus targeted for RNA extraction were preserved at −80°C in RNAlater-ICE (ThermoFisher Scientific) and dissected in chilled phosphate buffered saline (PBS).Amphipods were classified as "infected" if they were orange and carried at least one cyst, or as "uninfected" if they were brown and carried no cysts.A total of 24 infected and 24 uninfected samples were homogenized individually using a TissueLyser.
RNA was extracted using a Qiagen RNeasy Mini Kit and diluted to equimolar concentrations prior to sequencing.DNA and RNA libraries were sequenced by Genewiz (Cambridge, MA).Three DNA samples (uninfected male leg, uninfected female leg, and pooled cysts) were multiplexed and sequenced using 2 × 125 paired-end bp sequencing on two lanes of an Illumina HiSeq2500.Forty-eight RNA samples stemming from 24 infected and 24 uninfected amphipods were multiplexed across six lanes of a HiSeq2500 using 2 × 125 bp paired-end sequencing.

| Genome assembly and annotation
The O. grillus genome was assembled de novo using short reads and a de Bruijn graph approach implemented as described by (Ye et al., 2012) (k-mer size used = 51).Reads were quality controlled (QC) prior to assembly using the QC tools implemented in the String Graph Aligner (SGA) suite (complexity and homopolymer filter) (Simpson & Durbin, 2010) and BBmap (Bushnell, 2014).
Because the O. grillus individuals used for assembly are diploids collected from a natural population, we expected a high degree of heterozygosity in the assembly.We conducted genome scaffolding accounting for this high degree of variation using the algorithm Redundans implemented by (Pryszcz & Gabaldon, 2016).
Redundans uses information contained in the contigs and the paired-end reads to collapse heterozygous contigs into homozygous scaffolds.We estimated summary statistics of the assembly using Quast 5.0.0 (Gurevich et al., 2013).We benchmarked genome completion using Benchmarking Universal Single-Copy Ortholog (BUSCO 3.0.2) gene metrics (Simao et al., 2015).We predicted genetic features using the software Augustus 3.3 (Stanke et al., 2006).To improve gene prediction, we mapped all pooled amphipod reads from the 24 uninfected individuals to the de novo genome using STAR 2.6.1 (Dobin et al., 2013) and subsequently extracted intron-exon boundary hints.We chose to use only uninfected individuals to avoid the risk of trematode RNA reads introducing noise into the prediction.We utilized these hints to inform an Augustus analysis using a model trained on Drosophila melanogaster.This generated a gene feature file (GFF) containing all genomic features predicted in the genome.In order to identify homologous loci for all genome and transcriptome features, we conducted local BLAST runs against the whole translated transcriptome of Drosophila melanogaster (assembly release 6 plus ISO1 MT; Dmel6; NCBI ID: 47).We accomplished this by converting Dmel6's translated transcriptome into a BLAST-able database using the makeblastdb command in BLAST 2.7.1.BLAST searches were done using the blastp command enforcing a minimum e-value of 10 −10 .

| Read mapping, counting, and post-processing
Reads from each individual (either infected or infected) were mapped onto the de novo genome using STAR 2.6.1.Reads were mapped onto the de novo transcriptome using BWA-MEM 0.7.15 (Li, 2013), and were post-processed using samtools 1.9 (Li et al., 2009).The R package Rsubread 1.32.4 (Liao et al. 2019) was used to generate count tables for differential expression analysis, filtering for reads mapping to exon features contained in a GFF file produced by the Augustus gene prediction.For the population genetic analyses, we generated variant call files (vcf) using bamtools 1.9 (Barnett et al., 2011).Prior to vcf generation, we conducted stringent filtering on the bam files.We set a minimum mapping quality of 35 in all samples.We applied the following samflag filters to the bam: -f 0 × 0002 (keep reads mapped in proper pairs), -F 0 × 0004 (remove unmapped reads), -F 0 × 0008 (remove reads with unmapped mates).PCR duplicates were removed and bam files sorted by coordinate using picardtools (https://broad insti tute.github.io/picard/).Mapping RNA reads onto a reference genome results in gapped alignments and this may introduce problems when performing SNP calling at the intron-exon boundaries.To correct for this, we utilized the SplitNCigarReads function in GATK 4.0.9(https://gatk.broadinsti tute.org/).Variant calling was performed with bcftools (Li, 2020) using the multiallelic caller algorithm (−m) and a mutation rate prior of 2.8 × 10 −9 .

| Bioinformatic analysis: Differential expression and population analysis
Differential expression analysis was conducted using three approaches in R: DESeq2 (Anders & Huber, 2010), EBSeq (Leng et al., 2013), and edgeR (Robinson et al., 2010).The same read count table and gene annotation file was used across all three analyses.These files and the scripts for all differential expression analyses are available in the GitHub site listed below.Population genetic analyses were conducted in vcftools 0.1.16(Danecek et al., 2011).Statistical analyses were conducted in R using tools from the tidyverse suite 1.2.1 (Wickham et al., 2019) and adegenet 2.1.1 (Jombart, 2008).
VCF files were analysed using vcfR 1.8.0 (Knaus & Grunwald, 2017) and VariantAnnotation 1.28.2 (Obenchain et al., 2014).DSJ and TJF dedicate this manuscript to him and thank him for his mentorship and friendship these past decades.

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors have no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The genomic and transcriptomic data are available through NCBI Bioproject number PRJNA557566.The analyses are available at the

F
Field sampling and bioinformatic analyses of the amphipod-trematode system.(a) Samples were collected from The Plum Island Estuary -Long Term Ecological Reserve in Ipswich, Massachusetts USA (red box; see Methods for details).(b) The life cycle of the trematode Levinseniella byrdi through its three hosts: eggs are laid by avian hosts, ingested by the first intermediate host Hydrobiid snails, which pass the larval stages to the second intermediate host, the amphipod Orchestia grillus, which are then preyed upon by avian hosts.(c, d) Amphipods infected by L. byrdi change colour from light grey or brown to orange and move into more exposed areas of the marsh substrate, which may increase rates of predation (b-d from (Johnson & Heard, 2017)).(e) The experimental protocol for obtaining genomic and transcriptomic data, and the bioinformatic pipelines used in the data analyses.1365294x,2023, 18, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.17093by Mbl Whoi Library, Wiley Online Library on [23/07/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 5031 RAND et al.
an artefact of microhabitat effects influencing gene expression of the field-collected individuals.These patterns of differential gene expression are consistent with the known pigmentation and behavioural changes induced by F I G U R E 3 Differential expression of gene functional categories.(a) Venn diagram of the significantly differentially expressed transcripts from the three different analysis pipelines edgeR, EBseq, and DEseq2.There are 610 transcripts shared by at least two of the three pipelines, and 100's of transcripts unique to each pipeline.(b) Gene Ontology (GO) terms associated with significant up-or down-regulated transcripts among the 610 FDR-corrected differentially expressed genes shared by at least two of the three analysis pipelines.The GO categories identified using DESeq2 and edgeR are similar, while those from EBseq differ somewhat (see TableS2for details).GO annotations were based on homology of amphipod transcripts with Drosophila melanogaster gene sequence based on BLAST values of <10 −10 (see Methods).The notable categories are upregulation of 'detection of external stimuli', 'eye pigmentation' and 'pigment biosynthetic pathway', and down regulation of general metabolic processes and 'ion membrane transport'.These are associated with phenotypic effects of infection in adult amphipods.

F
Infected amphipods show downregulation in V-type proton ATPases.The heatmap plot shows down-or up-regulated transcripts across various subunits of V-type proton ATPases.Colours above columns indicate infection status.Colours on rows indicate different types of subunits, as noted in the legend at right.Parasite-induced down regulation of host V-ATPases has been demonstrated in other systems as playing a role in the inhibition of host immune responses against parasite infection (Robinson et al., 2012), implying a common role for these proteins in the mechanisms of host-parasite interactions.1365294x, 2023, 18, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.17093by Mbl Whoi Library, Wiley Online Library on [23/07/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 1365294x, 2023, 18, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.17093by Mbl Whoi Library, Wiley Online Library on [23/07/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 1365294x, 2023, 18, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.17093by Mbl Whoi Library, Wiley Online Library on [23/07/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Sample collectionOrchestia grillus amphipods were collected from Plum Island Estuary in northeastern Massachusetts, USA during October 2013.Specifically, we sampled the salt marshes adjacent to Club Head Creek (42°44′25.4″N 70°50′15.0″W) by manually sorting through vegetation and detritus for visually obvious amphipods.Individuals identified as putatively infected (i.e., orange) and uninfected (i.e., brown) O. grillus were collected in sterile microtubes and transported to Brown University.Amphipods were common gardened in the laboratory for 30 days in the same container with wet paper towels, in an effort to reduce any transcriptional effects associated with local field collecting conditions.Infection status was later confirmed by checking for cysts (see below).Samples were snap-frozen in liquid nitrogen and stored at −80°C until nucleic acid extraction.
1365294x, 2023, 18, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.17093by Mbl Whoi Library, Wiley Online Library on [23/07/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License The project was initiated by Cohort 3 of the NSF IGERT program in Reverse Ecology: Computational Integration of Genomes, Organisms and Environments (NSF DGE 0966060, PI: DMR, co-PI: ZGC).DMR, ZGC, AEG, KBN, ANS, WM, KK, DSJ conceived of the project and participated in the field work.KBN, ANS, WM maintained the animals and prepared RNA and DNA.TJF advised on trematode biology.JCBN, SBW, SR, JTB, WM, NS performed genome and transcriptome assembly and annotation.JCBN, SBW, SR, JTB, WM, YR performed bioinformatic analyses.DMR, KBN JCBN, SR, ANS wrote initial drafts of the manuscript.DAF, MGZ, AL, NO, DMM, KK, BRPB, AEG, ZGC, participated in project development and manuscript editing.DMR, DSJ, ZGC, JCBN, YR wrote or edited the submitted manuscript.ACK N O WLE D G E M ENTS This project was supported by an NSF IGERT award in Reverse Ecology: Computational Integration of Genomes, Organisms and Environments (NSF DGE 0966060), by NSF support for The Plum Island Ecosystems Long-Term Ecological Research site (NSF OCE 1637630).DMR acknowledges the support of NIH awards R01GM067862, 1R35GM139607 and P20GM109035, and DSJ acknowledges the support of NSF DEB 1902712.We thank Richard Heard for parasitological advice.Dr. Heard passed away in 2023.