Between‐population differences in constitutive and infection‐induced gene expression in threespine stickleback

Abstract Vertebrate immunity is a complex system consisting of a mix of constitutive and inducible defences. Furthermore, host immunity is subject to selective pressure from a range of parasites and pathogens which can produce variation in these defences across populations. As populations evolve immune responses to parasites, they may adapt via a combination of (1) constitutive differences, (2) shared inducible responses, or (3) divergent inducible responses. Here, we leverage a powerful natural host‐parasite model system (Gasterosteus aculeatus and Schistochephalus solidus) to tease apart the relative contributions of these three types of adaptations to among‐population divergence in response to parasites. Gene expression analyses revealed limited evidence of significant divergence in constitutive expression of immune defence, and strong signatures of conserved inducible responses to the parasite. Furthermore, our results highlight a handful of immune‐related genes which show divergent inducible responses which may contribute disproportionately to functional differences in infection success or failure. In addition to investigating variation in evolutionary adaptation to parasite selection, we also leverage this unique data set to improve understanding of cellular mechanisms underlying a putative resistance phenotype (fibrosis). Combined, our results provide a case study in evolutionary immunology showing that a very small number of genes may contribute to genotype differences in infection response.

resistance (Brunner et al., 2013;Zhang et al., 2019), or might adopt entirely different strategies (e.g., resistance vs. tolerance ;Sagonas et al., 2016). Thus, host populations often evolve divergent immune traits, whether they are faced with geographically varying parasite communities, or even a single widespread parasite.
Vertebrate immunity is a highly derived system consisting of a complex mix of defences and responses adapted to defend hosts against diverse pathogens and parasites (Boehm, 2012). These include constitutively expressed defences (i.e., physical barriers, antimicrobial peptides, phagocytosis, etc Paludan et al., 2021), as well as inducible responses, which are only initiated upon exposure to a potential pathogen (Frost, 1999;Paludan et al., 2021). Inducible responses can further be split into two somewhat overlapping categories: innate responses, which rapidly target broad pathogen categories (e.g., Toll-like receptor signalling; Kumar et al., 2011), and adaptive responses, which are slower to start but more fine-tuned to specific pathogens and involve retained pathogen memory (e.g., T and B cell mediated immunity; Mirzaei, 2020). All of these complex responses are potential targets of selection as a result of hostparasite interactions.
As populations evolve to adapt to a given parasite, to what extent does this evolution lead to between-population differences in constitutive defences, versus infection-induced responses? A simple strategy to study this question is to rear different host populations in a common garden environment, then experimentally infect some individuals and measure the resulting phenotypic variation. In particular, transcriptomic measures of gene relative expression are an increasingly common method to evaluate highly multivariate immune traits (Dheilly et al., 2014). One can then statistically partition the resulting variation in gene expression into relative contributions of genotype (G), environment (E, in particular, infection state), and genotype × environment interaction (G × E), to address the above question.
If the populations differ in heritable constitutive defences, then different genotypes raised in a common garden setting may exhibit fixed differences in expression (a genotype effect, G), insensitive to infection. This can happen, for instance, if one population consistently has higher numbers of certain types of immune cells, leading to higher tissue-level expression of genes unique to that cell type.
Numerous studies have documented variation in constitutive immunity as a driver of variation in response to parasites or pathogens (Ali et al., 2012;Evison et al., 2016;Kamiya et al., 2016;Schmitt et al., 2013). For example, genotypes of potato with higher resistance to the pathogen Phytophthora are characterized by increased constitutive expression of immune defences (Ali et al., 2012).
In contrast to constitutive variation, an environmental effect, in this case parasite infection, may trigger consistent inducible immune responses involving a cascade of changes in cell population abundance, signalling, and activation state (i.e., E). For example, social immune responses to mites in honeybees are conserved across multiple independent populations (Oddie et al., 2018).
Often these changes can be observed via gene expression (Ingham et al., 2008).
Induced responses to parasite infections may also be variable across populations (genotype by infection interaction, a particular case of more general G × E interactions). In such cases infection may change gene expression more strongly in one population than in another, or even in opposing directions. Plastic or induced responses are known to contribute to variation in parasite resistance (Reeson et al., 1998;Wakelin & Donachie, 1983). In Drosophila melanogaster, strains of flies resistant to a pathogenic bacteria are able to induce stronger and quicker innate immune responses than their susceptible counterparts (Okado et al., 2009).
Here, we present an estimate of the relative contributions of genotype, infection, and genotype × infection interactions to between-population differences in an immune organ's transcriptome, in threespine stickleback, Gasterosteus aculeatus. Sticklebacks' evolutionary history has yielded extensive natural variation in hostparasite interactions (Bolnick et al., 2020a(Bolnick et al., ,2020bPoulin et al., 2011;Young & Maccoll, 2017). Ancestral marine stickleback repeatedly colonized freshwater environments during the Pleistocene deglaciation. These new freshwater populations are exposed to a novel helminth cestode parasite, Schistocephalus solidus, which is absent in marine environments (McKinnon & Rundle, 2002;Simmonds & Barber, 2016) and rare in anadromous stickleback. Consequently, numerous genetically isolated freshwater populations of stickleback have been independently evolving in response to this freshwaterassociated parasite for thousands of generations (Weber, Kalbe, et al., 2017). The resulting replicated adaptation to S. solidus resulted in replicated evolution of greater resistance compared to marine fish (Weber, Kalbe, et al., 2017). But adaptation has not been entirely parallel: some populations evolved to drastically suppress cestode growth whereas others evolved to be relatively tolerant, reproducing despite the presence of large cestodes . Preliminary findings indicate that the growth suppression in these lakes is caused, in part, by the formation of peritoneal fibrotic scar tissue that traps the parasite (Lohman et al., 2017;Weber et al., 2021), and can sometimes kill the parasite. This among-lake variation in parasite resistance and fibrosis presents a powerful system for investigating relative contributions of G, E, and G × E, adaptations in driving between-population differences in immune phenotypes (e.g., resistance vs. tolerance).
Specifically, we bred F2 and backcross hybrids between a resistant population (Roberts Lake) and a tolerant population (Gosling Lake) in the laboratory, and experimentally exposed them to S. solidus. We then assayed infection outcome and stickleback gene expression in an immune organ (pronephros, a.k.a. head kidney). The head kidney is an important site of hematopoiesis and immune cell development in stickleback (Kum & Sekki, 2011). Transcriptomic analyses reveal a mix of G, E, and G × E effects, but expression variation was dominated by constitutive genetic variation (G) with very few G × E interactions. Furthermore, we highlight the transcriptomic signatures and putative immunological mechanisms which underly the putative resistance phenotype, fibrosis peritoneal fibrosis (De Lisle & Bolnick, 2020;Hund et al., 2020;Lohman et al., 2017;Weber et al., 2021). Combining these avenues of research, we provide new evidence suggesting diverse roles of dynamic adaptative responses of hosts to parasites, contributing to improved understanding of host-parasite evolutionary dynamics.

| Experimental design
We used minnow traps to capture reproductively mature stickleback from Roberts Lake and Gosling Lake, on Vancouver Island in British Columbia. These populations represent two ends of the natural spectrum of parasite prevalence: high parasite load in Gosling, low in Roberts . Furthermore, data suggest that Roberts Lake fish have evolved a fibrosis-based immune response to suppress parasite growth, which is not present in Gosling Lake fish (Hund et al., 2020;Lohman et al., 2017;Weber et al., 2021).
Wild-caught gravid females were stripped of eggs, which we fertilized using sperm obtained from macerated testes of males from the same lake (within-population crosses, denoted ROB or GOS) or the other lake (F1 hybrids, RG or GR depending on cross direction). Fish were collected with permission from the Ministry of Forests, Lands, and Natural Resource Operations of British Columbia (Scientific Fish Collection permit NA12-77018 and NA12-84188). The resulting eggs were shipped back to Austin, Texas, hatched, and reared to maturity. A subset of these first-generation laboratory-raised adults were experimentally infected with Schistocephalus solidus cestodes, or sham-exposed as a control. The resulting infection rates, cestode growth rates, and host immune traits are reported in Weber, Kalbe, et al., 2017, and host immune gene expression is described in Lohman et al. (2017). The remaining laboratory-reared F1 adults were artificially crossed among families (outbred) to generate F2 hybrids, including both intercrosses (F1 × F1 hybrids), and reciprocal backcrosses (ROB × F1 or GOS × F1). In contrast to Lohman et al. (2017), the fish examined here are second generation laboratoryraised individuals, allowing us to more confidently exclude maternal effects from wild-caught parents. We experimentally exposed 711 one-year-old F2 hybrids to S. solidus cestodes, following standard procedures (Weber, Kalbe, et al., 2017;. Briefly, we obtained mature cestodes from wild-caught stickleback from Gosling Lake or Echo Lake (Roberts Lake fish do not carry mature cestodes). We obtained the cestodes by dissecting freshly euthanized fish, then paired the cestodes by mass to mate them in nylon biopsy bags in artificial media, mimicking bird intestines where the cestodes typically mate (Wedekind et al., 1998). We collected the resulting eggs, and stored these at 4°C for up to one year. We hatched the eggs and fed them to Macrocyclops albidus copepods. The copepods were screened for successful infections after 14 days; then 5 infected copepods were fed to individually-isolated stickleback, as described in (Weber, Kalbe, et al., 2017;. We filtered the water after the exposure trial to ensure the copepods had been consumed. All F2 hybrid stickleback used in this trial were exposed to S. solidus (no sham exposures), to maximize infection rate for QTL mapping that has been described elsewhere (Ling et al., 2020).
However, only a subset of fish were successfully infected, providing a contrast between infected versus uninfected fish. Prior transcriptomic and flow cytometry data suggest that at the chosen time point (42 days post-exposure), fish with failed infections are phenotypically similar to sham exposed fish (Lohman et al., 2017). The experimentally infected fish were maintained for 42 days post-exposure, then euthanized with MS-222 and dissected to obtain (1) one head kidney (pronephros) for flow cytometry; (2) one head kidney for gene expression analysis, preserved in RNAlater at -80°C; (3) fish mass and length and sex; (4) the mass and number of successfully established cestodes, and (5) the presence or absence of fibrosis.
All fish handling was approved by the University of Texas IACUC (AUP-2010-00024).

| Flow cytometry
Flow cytometry data on head kidney cell population ratios (granulocytes versus lymphocytes) and activity (baseline ROS and oxidative burst) were generated following methods described by Weber, Kalbe, et al. (2017)), ). Data were analysed using flowjo software (Treestar). Populations of granulocytes and lymphocytes were separated by linear forward scatter (FSC) and side scatter (SSC), providing counts of the relative abundance of each cell type. ROS production by granulocytes was measured following protocols for PMA stimulation described in Weber, Kalbe, et al. (2017) and .

| RNA extraction and transcriptome sequencing
We extracted RNA from one head kidney using the Ambion MagMAX-96 Total RNA Isolation Kit, following a modified version of the manufacturer's protocol (see Supporting Information). Each head kidney (hereafter "sample") was separately placed in lysis/binding solution and homogenized using a motorized pestle. After initial purification using magnetic beads provided by the kit, DNA was removed by adding TURBO DNAse and a second purification with Serapure magnetic beads, leaving only RNA. The RNA yield of each sample was quantified using a Tecan NanoQuant Plate.
RNAseq libraries were constructed using TagSeq methodologies detailed in Lohman et al. (2016) with modifications. After fragmentation of the RNA in a magnesium buffer (NEB Next RNA fragmentation buffer), the RNA fragments were purified using Agencourt RNAClean XP beads. A poly-dT primer (3ILL-30TV) was annealed to the poly-A tail of mRNAs, after which the first cDNA strand was synthesized, which was amplified in a second PCR reaction. The PCR products were purified with Serapure magnetic beads, quantified (Quant-IT PicoGreen) and normalized (1 ng/µl), after which all libraries were PCR-barcoded using Illumina i5 and i7 indexes.
Fragment size selection occurred via automated gel extraction and final quantification was performed using qubit 2.0. The libraries were sequenced using a hiseq 2500 at the Genomics Sequencing and Analysis Facility of the University of Texas at Austin.
Samples with less than 500,000 mapped reads were removed from subsequent analyses, resulting in a final n = 390. Finally, we annotated transcripts with a blastx comparison to the uniprotkb database (http:// www.unipr ot.org/help/unipr otkb) with the parameters: max target seqs = 10; evalue = 1e −5 . Results were filtered to obtain the match with the highest evalue and bit score for each transcript.

| Analysis with DESeq2
Raw data, including read count matrixes and metadata, as well as code for all of the following analyses can be found on GitHub (https:// github.com/lfues s/TagSeqMS). To test for differential expression, we used the r package DESeq2. Transcripts were filtered to remove those that were not expressed in more than 195 samples (approximately half of the sample set). The remaining 15,354 sequences were tested for differential expression using the following model: where Y ij is the count of transcript I in individual j, β Room is a fixed effect with two levels corresponding to the room in which fish were reared, β Cross is a fixed effect with three levels: F2, GBC, and RBC, β Infection is a fixed effect with two levels: infected or uninfected, β Fibrosis is a fixed effect with two levels: fibrotic or nonfibrotic, β ROS is a continuous factor corresponding to measure reactive oxygen production per sample, and β Sex is a fixed effect with three levels: male, female, or unknown (for the few samples where sex could not be identified with confidence). β Batch is a random effect corresponding to the lane of sequence sampling.
Fibrosis and infection effects were tested independently from one another (i.e., not all fibrotic fish were labelled infected). All p-values were multiple test corrected using a 10% FDR (Benjamini & Hochberg, 1995).

| Expression pathway and upstream regulator analyses
Differentially activated biological pathways and upstream regulators were assessed for model factors of interest using the ingenuity pathway analysis software (IPA; Qiagen Inc., https://www.qiage nbioi nform atics.com/produ cts/ingen uity-pathw ay-analysis). For each factor of interest, logfold change, unadjusted p-value, and annotation (spID) per transcript was input into the software. Transcripts were then filtered to retain only those with unadjusted p-values < .05. Transcripts with duplicated IDs were averaged for analysis. Pathway analyses were used to identify pathways significantly affected by each factor, and generate corresponding p-values and activation scores (z-scores).
Factors tested using IPA were infection, fibrosis, and cross.

| Experimental infection outcomes
Of the 711 fish exposed to cestodes, 268 were infected with S. solidus after 42 days. Additionally, peritoneal fibrosis was observed in

| Constitutive variation among populations (G)
To test for signatures of variation in constitutive immune defences among our populations, we considered baseline variation in gene expression among our three crosses (main effect of host genotype, G, irrespective of infection status). Constitutive gene expression varied considerably. Of the 15,354 genes tested, 11,321 varied constitutively between at least two given cross types: 7745, 10,601, and 1202 transcripts were differentially expressed between F2 versus GBC, F2 versus RBC, and RBC versus GBC, respectively (p adj < .10, 10% FDR; Figure 1; Supporting Information S1). A total of 1,216 of these differentially expressed genes have immunological functions (Table 1). However this is not more immune-related genes than expected from chance (chi-squared test; p = .2627, χ 2 = 1.2543, df = 1).
In addition to broad variation among cross-types in gene expression, numerous pathways and upstream regulators, including many involved in immunity, were also significantly differentially expressed between crosses (p adj < .10, 10% FDR), independent of infection status. A total of 240, 313, and 94 pathways, representative of many different immune components, varied in activation state when comparing F2 versus GBC, F2 versus RBC, and RBC versus GBC respectively (Table 2; Figure S1; Supporting Information S2).

| General response to infection (E)
Host response to infection by a common cestode parasite, S. solidus, involved multiple genes and pathways. Comparing 158 infected versus 232 uninfected fish (all three cross types), we found 2,369 differentially expressed transcripts (p adj < .10, 10% FDR; Supporting Information S1), 2,223 of which were annotated. Of these, 341 transcripts were annotated to encode for proteins involved in immunity. Biological pathway and upstream regulator analyses also indicated broad effects of infection on hosts. A total of 169 pathways were significantly activated as a result of infection (p adj < .10, Supporting Information S2). Thirty-one of these pathways were linked to immunity ( Figure 2). All of these significant immune pathways were suppressed (lower relative transcript abundance) in infected fish relative to uninfected fish. Upstream regulator analysis identified 121 regulators which were differentially activated as a result of infection, 27 of which had roles in immunity ( Figure S2, Supporting Information S3). Similar to pathway results, regulators displayed broad patterns of immune suppression.
A large fraction of the infection-responsive genes (1,812 of the 2,369 differentially expressed genes associated with infection, or 76%) also varied in expression between two or more cross types (e.g., genes that exhibit both G and E effects). However, the overlap between the infection-responsive and constitutively-different genes was no greater than random expectations (p = .2821, χ 2 = 1.1569, df = 1). A particularly noteworthy result is that genes with higher average expression in GBC versus RBC (G effect) are more likely to exhibit a negative response to infection (decreased expression) in both genotypes. Specifically, we found an inverse correlation between directional effect sizes of genotype versus infection main effects for the subset of genes significantly differentially expressed both in infected fish and between RBC versus GBC fish (p = .0005439, Our results showed little overlap with those from previous transcriptomic study of the response to S. solidus infection in pure-cross families (as opposed to hybrids) of these same stickleback populations (Lohman et al., 2017). Fourteen genes were shared between the studies; all but four responded in similar directions (Table 3). Additionally, 76 of the genes that were differentially expressed here were also F I G U R E 1 Venn diagram of overlap in significantly differentially expressed genes among all three cross type comparisons TA B L E 1 Example list of genes differentially expressed among cross types with putative functions in immunity. A full list of differentially expressed genes for each contrast can be found in Supporting Information S1

| Variation in induced immune responses between genotypes (G × E)
A total of 569 genes exhibited interactions between cross type and infection; however, the vast majority of these were only significant when considering differences in response to infec-

| Transcriptomic signatures of fibrosis
We also leveraged these data to identify changes in expression associated with a putative resistance phenotype: fibrosis. We identified strong transcriptomic signatures of fibrosis: 5826 genes were differentially expressed in fibrotic fish compared to those not displaying the fibrosis phenotype, 840 of which have putative roles in immunity and fibrosis. Similar to the effects of infection, the majority of these genes were downregulated in fibrotic fish. Pathway and upstream regulator analyses also revealed significant downregulation of immune-related processes in fibrotic fish (Figure 4; Figure   S2). Many (257) biological pathways were significantly differentially activated in fibrotic fish, 53 of which are related to immunity, defence, or fibrosis responses. Most of these immune-related pathways were also significantly differentially activated in fish infected with S. solidus. However, these shared pathways were almost always TA B L E 3 Comparison of infection-associated significantly differentially expressed genes to results from a previous study (Lohman et al., 2017) using the same two source populations  Figure 5).

| Variation in constitutive expression of immune defences across genotypes (G)
Here we document significant variation in gene expression as a result of both genetic variation between two populations, and cestode infection (an environmental effect), and their interaction. By far the dominant source of variation in gene expression arose from heritable differences between genotypes (Roberts Lake backcrosses, F2 hybrids, and Gosling Lake backcrosses), irrespective of infection status.
Given the large number of heritable differences in gene expression between the populations, it is challenging to interpret these differences in a functional framework. The expression differences may be adaptive or neutral. If adaptive they may serve any of a vast number of functions given the many phenotypes that differ between these lakes including skeletal morphology, size, behaviour, male breeding colours, diet, microbiota, diverse species of parasites besides S. solidus (Bolnick et al., 2020b;Ling et al., 2020;Snowberg et al., 2015). Constitutive expression differences do contain many immune genes (  (Duan et al., 2021;Xie et al., 2021). When considering pathways and upstream regulators, B cell receptor signalling was activated in GBC compared to RBC fish whereas RBC fish activated both IL-7 and IL-3 signalling pathways compared to GBC fish.
B cells are an essential component of fish immunity (Sunyer, 2012), including serving roles in the defence against parasite-induced proliferative kidney disease in salmonids (Abos et al., 2018;Bailey et al., 2017). Additionally, localized expression of interleukins is known to play a role in fish-parasite interactions (Perez-Cordon et al., 2014).
In vertebrates, IL-3 specifically plays roles in development of essential antiparasitic immune cells such as mast cells and basophils (Lantz et al., 1998). Finally, RBC fish had higher activation of upstream regulator TGF-beta 1 compared to GBC fish. TGF-beta 1 is an important regulator of fibrosis (Gressner et al., 2002), which is an important parasite resistance phenotype in stickleback (see later sections; Weber et al., 2021). Consequently, while not statistically overrepresented as a gene ontology category, some of the observed differences in constitutive expression of immune genes may indeed contribute to divergence in parasite responses.

| Conserved responses to a cestode parasite (E)
The next most important factor in our data was a main effect of infection (an environmental effect) irrespective of genotype. These genes represent a multigenic and multipathway response to infection representative of a diversity of immune components. This type of multifaceted response is common across systems: host response to parasite infection often involves multiple arms of the immune system (Anthony et al., 2007;McRae et al., 2015;Medzhitov, 2007;Nakada-Tsukui & Nozaki, 2016). were also suppressed including Fcy receptor-mediated phagocytosis in macrophages and monocytes, phagosome formation, and phagosome maturation. Upstream regulator pathway analysis confirmed these patterns of broadscale reduction in immunity in infected fish.
Antifibrotic regulator hepatic nuclear factor 4-alpha, HNF4A (Ji et al., 2018;Yue et al., 2010) was significantly activated as a result of infection. Intriguingly, HNF4A is also a major target of natural selection driven differences between Roberts and Gosling Lake, falling within the quantitative trait locus for production of granulomas that encase and frequently kill the cestode (Weber et al., 2021). Several important adaptive immune components such as interleukin-4 (IL4; (Heeb et al., 2020) and IgE complex (Galli & Tsai, 2012) were also suppressed in infected fish. Together, gene expression and pathway analyses reveal wide patterns of immune suppression associated with infection of G. aculeatus with S. solidus. This suppression in infected fish could arise from three processes: initially immunesuppressed fish might have been more vulnerable to infection, fish observed to be infected might be those that have adapted tolerance rather than resistance strategies, or successfully established parasites may be actively suppressing host immunity. Time-series analyses of individual host expression could distinguish these alternatives, but is not practical because fish must be euthanized to obtain head kidney samples.
Comparison of the results of our study to previous transcriptomic analyses of stickleback response to S. solidus revealed limited overlap in infection-responsive genes. When comparing our results to those from a previous study of fish from pure-cross families generated from the same populations (Lohman et al., 2017) we found limited shared genes, though the patterns of expression within this small group were mostly congruent (Table 3). Two shared genes of particular interest were interleukin 8 and fibronectin (also a genomic target of selection within one of the resistance QTL mapped by Weber et al. (2021). Infection was associated with higher expression

| Variation in response to inducible responses across genotypes (G × E)
A third potential source of variation in response to parasite expo- through regulation of CD4 + T cell functioning (Crotty, 2011;Crotty et al., 2003). Furthermore, SH2D1A may affect NKT cell processes (Nichols et al., 2005;Wu et al., 2016) and differentiation of T H 2 cells (Wu et al., 2001). Interestingly, expression of this key component

| Expression changes associated with fibrosis
In addition to evaluating the roles of variation in G, E, G × E responses in divergent parasite responses, our findings also shed light on a putative resistance phenotype in this system: fibrosis. Fibrosis is a common immune pathology across vertebrates (Sgalla et al., 2016;Vrtílek & Bolnick, 2020;Wick et al., 2010), frequently associated with parasitic infections (Kamdem et al., 2018;Niu et al., 2019;Wilson et al., 2007). Often excessive fibrotic responses can cause health issues, including in humans (Henderson et al., 2020;Todd et al., 2012;Wynn, 2008). In G. aculeatus, recent study has demonstrated that peritoneal fibrosis is an induced response to S. solidus in some stickleback populations (Hund et al., 2020), and is associated with reduced cestode growth (Weber et al., 2021). More broadly, teleost fish in general are almost all susceptible to peritoneal fibrosis in response to an immune challenge, however little mechanistic knowledge exists regarding the cellular processes controlling these F I G U R E 6 Summary of major changes in gene expression, pathways, and upstream regulators that may affect the resistance phenotype, fibrosis. Lines above each box demonstrate up or down regulation of each component in fibrotic (blue) or infected (red) fish. Lines from each box show the predicted relationship between a given component and the fibrosis phenotype (based on existing literature). In each case arrows indicate upregulation/activation whereas blunt end lines indicate downregulation/inhibition responses (Vrtílek & Bolnick, 2020 (Liu et al., 2018;Xavier et al., 2017). Similar to the effects of infection, the majority of these genes were downregulated in fibrotic fish. We also observed 36 collagen genes that were differentially expressed in fibrotic fish, almost all of which were also downregulated. Fibrotic tissue is formed via the deposition of collagen and other extracellular matrix components (Henderson et al., 2020;Wynn, 2008). The large-scale downregulation of these complement and collagen genes, as well as other immune and fibrosisrelated transcripts, in fibrotic fish may be an artefact of the time of sampling: 42 days after exposure. It is possible that these genes are activated earlier on in the infection time course, but have since been downregulated, potentially as a mechanism for attenuation of fibrosis (Hund et al., 2020). Future time-course expression studies are needed to clarify the temporal dynamics of gene expression during fibrosis.
Pathway and upstream regulator analysis also revealed downregulation of immune-related process in fibrotic fish, however this downregulation was substantially less than downregulation of the same pathways in infected fish. Many of these pathways differentially expressed in both fibrotic and infected fish have potential ties to fibrotic processes, including hepatic fibrosis signalling pathway and IGF-1 signalling. IGF-1 stimulates the differentiation of fibroblasts to promote a potent fibrosis response (Hung et al., 2013). Furthermore, antifibrotic regulator hepatic nuclear factor 4-alpha, HNF4A (Ji et al., 2018;Yue et al., 2010) was significantly activated as a result of infection, but suppressed in fibrotic fish. Combined these results suggest that while infection is marked by a general reduction in immune pathways, this suppression is weaker or absent in fibrotic fish that more effectively reduce cestode growth ( Figure 6).
The association between infection and immunosuppression is particularly interesting as immunosuppression is widely known in helminths in general (Maizels et al., 2004(Maizels et al., , 2018Maizels & McSorley, 2016), and has been documented for S. solidus in particular (Scharsack et al., 2004). It is possible the described patterns are indicative of cestode suppression of host immune response during infection. In contrast, fibrotic fish show weaker signatures of immune reduction, and in some instances demonstrate patterns of expression opposite to those displayed by infected, nonfibrotic fish. Furthermore, these specific pathways and genes that demonstrate disparate patterns between fibrotic and infected-but-nonfibrotic fish (e.x. HNF4A), are also highly variable among cross-types. Combined, these analyses suggest that cestodes may act to suppress immune responses in their host, but that this immunosuppression differs between host genotypes, and between fibrotic and nonfibrotic fish. These findings suggest that evolution of resistance may also be dependent on the acquisition of traits to overcome or avoid parasite manipulation of host immunity.

| CON CLUS IONS
Vertebrate immunity is a highly derived, complex system of constitutive and inducible immune responses. Consequently, the effects of selection induced by novel parasites and pathogens can be stochastic and difficult to quantify. Here we provide evidence that divergence in response to parasite selection is a mix of conserved responses and variable inducible responses to the parasite. While we document strong variation in general constitutive expression among cross types, immune defence genes are not statistically overrepresented among these constitutive genetic differences. In contrast we observe a small number of genes closely linked to parasite resistance and immunity which vary in expression in response to infection among our cross types (G × E). The ratio of significantly differentially expressed genes with G, E, and G × E effects is 11,321:2369:569, showing that constitutive genetic divergence is the dominant contribution to expression variation in this system. However, when considering only genes with documented immunological roles (i.e., genes with immune annotations), the ratio is 1216:342:65.
Comparing this ratio to the G:E:G × E ratios for all gene annotations, inducible responses (while still a minority) are statistically overrepresented. From this observation we conclude that, at least in this system, inducible responses, as opposed to constitutive defences, contributes disproportionately to variation in immune gene expression. Genotype-specific immune responses affected a very small minority of genes, though these may nevertheless be functionally quite important. In addition to providing further clarification regarding mechanisms of divergent adaptation to parasites, we also identify transcriptomic signatures contributing to a putative resistance phenotype, fibrosis. Not only do we identify key genes and pathways associated with this response, but we also provide evidence that fibrotic (i.e., resistant) fish apparently are refractory to immune suppression associated with infection in other fish. This suggests host evolution of counter mechanisms may also be key in the evolution of resistance.

DATA AVA I L A B I L I T Y S TAT E M E N T
Raw data, including read count matrixes and metadata, and code can be found on GitHub (https://github.com/lfues s/TagSeqMS).