Functional immunogenetic variation, rather than local adaptation, predicts ectoparasite infection intensity in a model fish species.

Natural host populations differ in their susceptibility to infection by parasites, and these intrapopulation differences are still an incompletely understood component of host-parasite dynamics. In this study, we used controlled infection experiments with wild-caught guppies (Poecilia reticulata) and their ectoparasite Gyrodactylus turnbulli to investigate the roles of local adaptation and host genetic composition (immunogenetic and neutral) in explaining differences in susceptibility to infection. We found differences between our four study host populations that were consistent between two parasite source populations, with no indication of local adaptation by either host or parasite at two tested spatial scales. Greater values of host population genetic variability metrics broadly aligned with lower population mean infection intensity, with the best alignments associated with major histocompatibility complex (MHC) "supertypes". Controlling for intrapopulation differences and potential inbreeding variance, we found a significant negative relationship between individual-level functional MHC variability and infection: fish carrying more MHC supertypes experienced infections of lower severity, with limited evidence for supertype-specific effects. We conclude that population-level differences in host infection susceptibility probably reflect variation in parasite selective pressure and/or host evolutionary potential, underpinned by functional immunogenetic variation.

INTRODUCTION

While hosts are under selection to combat pathogens, pathogens are under concurrent selection to overcome host defences (Schmid‐Hempel, 2011; Woolhouse et al., 2002). This interaction is ubiquitous (e.g., Poulin, 1999) and has consequences ranging across the evolution of host immune systems (Frank, 2002), epidemiology and the emergence of new infectious diseases (Ebert, 1994; Lively, 2016), adaptive radiation (Karvonen & Seehausen, 2012), the maintenance of sex (Morran et al., 2011), and conservation biology (Altizer et al., 2003). Hosts within a population typically vary in their susceptibility to infection (Woolhouse et al., 2002), and, in spatially heterogeneous host‐pathogen systems, host populations often differ in their average observed susceptibility (Brunner et al., 2017; Eizaguirre & Lenz, 2010; Lively & Dybdahl, 2000). Integrating these within‐ and between‐population interactions is important for understanding the evolutionary and epidemiological consequences of these dynamic processes (Brunner et al., 2017; Carlsson‐Granér & Thrall, 2002; El Nagar & MacColl, 2016; Hess, 1996; Penczykowski et al., 2016; Schneider et al., 2017; Smith et al., 2003; Soubeyrand et al., 2009; Thompson, 2005). Host adaptation to local parasites, for example, may constitute a barrier to gene flow between host populations, facilitating speciation (El Nagar & MacColl, 2016), or increase host susceptibility to parasites transmitted from distant populations or different species (Daszak et al., 2000).

Less genetically variable host populations are often reported to be more susceptible to infection (e.g., Gibson & Nguyen, 2020). Such associations may arise if genetically homogenous hosts are easier for parasites to adapt to (reviewed in refs King & Lively, 2012; Radwan et al., 2010) and/or because individuals with higher homozygosity – genome‐wide or at immunity genes – are more susceptible (e.g., Acevedo‐Whitehouse et al., 2003; Luikart et al., 2008; Ortego et al., 2007). Differences between hosts and parasites in the processes contributing to genetic diversity may thus play a critical role in local coevolutionary outcomes. In turn, this may help explain why, while shorter generation times and larger populations should give parasites the edge over hosts in local adaptation arms races (Gandon & Michalakis, 2002; Price, 1980), the majority of reciprocal infection experiments report no significant local adaptation by either parasite or host (approximately 56%; Greischar & Koskella, 2007). Complementary tests of local reciprocal adaptation alongside data on immunogenetic markers thus offer a potentially useful, but underutilised, approach for understanding patterns of infection within and among host populations.

Here, we performed a controlled infection test for local adaptation using a model fish‐ectoparasite system, complemented with the study of two sets of highly polymorphic genetic markers, one presumed to be neutral and the other known to be under intense selection from parasites. For the neutral marker set, we used microsatellites: a well‐characterised and well‐utilised set is available for our host species, and their polyallelic nature makes them useful in direct comparisons with our marker under selection, the major histocompatibility complex (MHC). Genes of the MHC encode molecules involved in immune responses in vertebrates (Klein, 1986), and decades of research has been devoted to the complex suite of selection pressures that maintain and promote the gene family's extreme polymorphism, which includes parasite‐mediated selection (Radwan et al., 2020; Snell, 1968; Spurgin & Richardson, 2010), sexual selection (Ejsmond et al., 2014; Penn & Potts, 1999), and selection acting on the MHC‐linked sheltered load (van Oosterhout, 2009). The ecological pertinence of the MHC is well established, including numerous studies reporting associations between MHC alleles and resistance/susceptibility to parasites in wild systems (e.g., Buczek et al., 2016; Fraser & Neff, 2010; Kaufmann et al., 2017; Schad et al., 2005). The role of MHC genes and MHC variability in causing differences in host resistance/susceptibility between populations is less well understood, with the hypothesis that populations with more MHC variants have lower parasite burdens supported by observational evidence (Meyer‐Lucht & Sommer, 2009), mesocosm experiments (e.g., Eizaguirre et al., 2012), and wild cage experiments (e.g., Bolnick & Stutz, 2017), but limited exposure‐controlled experimental testing (but see Smallbone et al., 2021). An important concept in the study of MHC evolution is that of “supertypes” (STs), groups of MHC alleles that encode peptides with similar antigen‐binding properties. STs may better capture the functional breadth of host defence than alleles or phylogenetic groupings. In ecological MHC studies, STs are usually assigned by statistical clustering (see Materials and Methods), but the concept is founded in laboratory immunology (Doytchinova & Flower, 2005; Sandberg et al., 1998; Sidney et al., 1996).

Guppies (Poecilia reticulata) are tropical freshwater fish native to northern South America and the Caribbean, and have been an important model species in elucidating processes as diverse as sexual selection, predator‐prey interactions, ecological competition, and, most relevant to the present study, host‐parasite dynamics. Furthermore, their MHC has been well‐characterised and well‐studied (e.g., Fraser & Neff, 2010; Smallbone et al., 2021). Monogenean ectoparasites in the genus Gyrodactylus are widespread across wild guppy populations, but their prevalence varies greatly between and within populations, and through time (Dargent et al., 2013; Mohammed et al., 2020; van Oosterhout et al., 2003; Stephenson et al., 2015, 2017). The known pathogenicity of the parasites, coupled with the relative ease with which they can be maintained in a laboratory and used in exposure‐controlled infection trials, make the guppy‐Gyrodactylus system an excellent model for studying a wide range of host‐parasite interactions, including the effects of parasitism on some of the processes described above (e.g., sexual selection, predator‐prey; reviewed by Bakke et al., 2007). In the present study, we use this highly tractable system to investigate (i) whether wild host populations show consistent differences in resistance across parasite strains, (ii) whether parasites are adapted to their local hosts and vice versa, (iii) what role interpopulation variation in MHC traits may play in this dynamic, and (iv) whether MHC genotypes predict infection intensity at the individual level and if this varies between host populations/parasite lineages.

MATERIALS AND METHODS

Host collection and rearing

We collected (hand seine) juvenile guppies (standard length 5–12 mm) from four wild populations in March 2016, two on Trinidad (Lopinot,”Lop”; Santa Cruz, “SC”) and two on Tobago (Scarborough Health Centre,”HC”; Roxborough, “Rox”; Table S1.1). Surveys undertaken as preliminary work for Phillips et al. (2018) had shown that all sites had populations of gyrodactylids, and previous population genetic analyses have shown significant neutral and MHC differentiation between all host population pairs, with differentiation stronger between islands than within islands (Herdegen‐Radwan et al., 2020; Phillips et al., 2018). At our field station in Tobago, we treated all fish with salt water (15 ppt, 5 min) to kill any gyrodactylids (Schelkle et al., 2011), confirmed by briefly anaesthetising all fish (0.02% tricaine methanesulphonate; “MS‐222”) and screening them under a dissecting microscope with cold illumination multiple times over several days, according to Schelkle et al. (2009). We saw no signs of other ectoparasites or disease at any point in the study. Each population was then reared in a separate aquarium (80 L, 50–100 fish per aquarium) and fed daily with live Artemia nauplii and generic, pet‐shop fish flakes (Aquarian).

Parasite collection and rearing

In June 2016 we returned to each of the four sites and collected 50–60 guppies to act as gyrodactylid donors for our experiment. The prevalence and intensity of gyrodactylid infections on these fish were too low for our preferred protocol of infecting experimental fish with parasites straight from the wild, and for sourcing parasites from all four populations, so we cultured Gyrodactylus turnbulli from populations Lop and HC using fish from a gyrodactylid‐free captive population. Details of gyrodactylid species identification are in Appendix S1 (see also Cable & van Oosterhout, 2007a, 2007b; King et al., 2009; Xavier et al., 2015). This “farm” host population, maintained in an 800L mesocosm, had been founded 18 months earlier by crossing captive‐reared virgin females from a Tobagonian population with males from a Trinidadian population (Appendix S1), and had been verified as gyrodactylid‐free at P, F1 and F2 generations. Neither founding population of the farm stock features in the present study. Farm gyrodactylid lineages were established by briefly anaesthetising both a wild donor and a recipient fish, and, under a dissecting microscope with cold illumination, bringing together their caudal fins (tail) until a single gyrodactylid moved from the donor to the recipient. Any extra gyrodactylids that jumped were removed with watchmaker's forceps. Donor and recipient were then separated and revived, and the recipient was moved to a 500 ml isolation container. After six days the procedure was repeated, using infected farm hosts to make single‐gyrodactylid infections on a fresh batch of farm recipients. Gyrodactylid cultures were subsequently maintained by keeping 1–3 parasite‐naïve recipients in an isolation container with an infected donor for 3–4 h, then moving each new recipient to its own isolation container (after Stewart et al., 2017). Farm hosts were fed fish flakes daily, with water changed every other day. All gyrodactylids used in experimental infections could be traced to their original wild founder.

Experimental design

We performed exposure‐controlled gyrodactylid infection trials (Cable & van Oosterhout, 2007a, 2007b) on the fish captured as juveniles in March 2016. Experimental infections were established by briefly anaesthetising parasite donor (infected farm fish) and recipient and allowing two gyrodactylids to move to the recipient (any extras were removed – see above). Recipients were measured (standard length) before infection. All infected fish were females, with length ≥15.0 mm. Each experimental host was kept in its own 500 ml isolation container at ambient shade temperature and fed with fish flakes every other day. The day after infection (day 1), we anaesthetised each experimental host and counted the number of gyrodactylids it carried and repeated this every other day thereafter until day 17 or until the fish was observed to be gyrodactylid‐free for five days. A typical gyrodactylid infection in a parasite‐naïve wild guppy will start with exponential growth, peak after 7–11 days, and rapidly fall away to single figures of parasites (e.g., Phillips et al., 2018; see also Figures S11.1, S11.2). Each fish received a water change at the time of screening. Every 12 h, we checked whether all fish were alive, increasing this to every 4 h if a host's infection intensity rose above 70 gyrodactylids. Fish found dead were promptly preserved in 1 ml absolute ethanol (changed as in Appendix S1). We preserved a fin clip (caudal fin, 2–4 mm2; 0.3 ml absolute ethanol) from all fish that survived the experiment. Any fish that cleared its infection within the first seven days was reinfected 4–6 days after first being recorded as “clear”, as rapid clearance may be a stochastic effect of gyrodactylid quality (further details below). We initiated experimental infections in two blocks four days apart, balanced by fish and gyrodactylid population. All reinfections were initiated six days after the second block.

Genetic analyses

We extracted genomic DNA (20–100 ng per sample) from guppy fin clips using MagJET Genomic DNA kits (Thermo Scientific). We then PCR amplified a 217 bp fragment of the MHC class II second exon, which codes for the highly polymorphic β‐chain of the MHC molecule's antigen binding groove, using the primers of Phillips et al. (2018) designed from mRNA fragments conserved among poecilids (PCR conditions as Phillips et al., 2018). The PCR included the fusion primers required for DNA sequencing with an IonTorrent Personal Genome Machine (PGM; Life Technologies), as well as a unique combination of 6 bp tags (20 tags = 400 potential F×R combinations) for each amplicon (fish). After amplification, we pooled amplicons approximately equimolarly and sequenced the pool (PGM). We then used the adjustable clustering method of Biedrzycka et al. (2017), implemented in the software amplisas (Sebastian et al., 2016; parameters as Phillips et al., 2018), to turn raw sequence data into individual genotypes. We followed this up by allocating MHC alleles to the supertype (ST) groups of Herdegen‐Radwan et al. (2020), based on physicochemical properties at positively selected sites reported by Lighten et al. (2017).

Copy number variation at the MHC means that the locus affinity (phasing) of alleles is rarely known, but an earlier cross‐breeding experiment (Phillips et al., 2018) allowed us to phase all of the alleles in this study's focal populations. That earlier experiment reported a single linkage block of 1–3 alleles in a de facto single locus, although only 1–2 alleles per block feature in the present study. We use the term “superhaplotype” to refer to haplotypes based on STs rather than alleles.

To provide a proxy for neutral genetic variation, all hosts in the experiment were genotyped at 15 microsatellite loci (Becher et al., 2002; Olendorf et al., 2004; Shen et al., 2007; Watanabe et al., 2003). Of these, eight were retained for the main analyses, with 4/15 dropped for <50% amplification success in at least one population and 3/15 dropped for significant departure from Hardy‐Weinberg equilibrium (Appendix S2). These loci are routinely used in guppy behavioural ecology and population genetics, including for comparisons against MHC variability (e.g., Herdegen‐Radwan et al., 2020; Lighten et al., 2017).

Tests of host and parasite populations

To assess how host population and parasite population affected the outcome of infection trials, without considering any explicitly genetic predictors, we first tested for biases in host death rate. For this, we used contingency table‐based analyses (χ2 and Fisher's exact tests), as we considered the death rate (5/114 fish, 4.4%) too low for logistic regression. We then tested for effects of host population and parasite population on the number of “worm‐days” experienced by fish that survived the experiment. Worm‐days were calculated as the area under a fish's infection trajectory graph (number of gyrodactylids against time), and are both ecologically pertinent and statistically tractable – more worm‐days can reasonably be considered a more intense infection, and the metric avoids needing to consider time series, temporal autocorrelation, zero inflation, or individual‐level random effects (Phillips et al., 2018). For fish that were reinfected, we retained the infection that reached the highest peak intensity (details in Appendix S3).

Worm‐days were analysed using general linear models (LMs; Gaussian errors) on log ‐transformed worm‐days, and multimodel inference implemented in the MuMIn package (Bartoń, 2016) of the statistical software “r” (R Development Core Team, 2016). We opted for LMs over the generalised linear models with negative binomial errors used by Phillips et al. (2018) because: (i) the interpretation of LMs is usually more intuitive (e.g., use of R 2 to quantify the proportion of explained variation); (ii) their post hoc options are more tractable and more widely known; and (iii) in the present study, LMs tended to produce residuals that slightly better reflected a normal distribution. In the Supporting Information, we show that no interpretive differences would have arisen had we used negative binomial models (Appendix S4). We used corrected Akaike Information Criterion (AICC) to rank models with all combinations of the following parameters: host population (factor, four levels); parasite population (factor, two levels); the interaction between host and parasite population; temperature (factor, three levels corresponding to date of infection; empirical temperature data and rationale for factor in Appendix S5); and fish size (standard length; continuous, z‐transformed). If the top‐ranked model was more than two AICC units clear of the second model, we deemed it the nominal best model and examined its coefficients for values significantly different from zero. If more than one model comprised the top two units of AICC, we used AICC‐weighted model averaging to estimate coefficients and their significance, implemented in MuMIn. We inspected LM assumptions for all models by adding five LM diagnostic statistics to the summary tables produced by MuMIn: skewness and kurtosis of residuals; Kappa and maximum variance inflation factor to assess multicollinearity; and the maximum Cook's distance value to check for influential data points.

As alternatives to using host population and parasite population, we tested two predictors based on host‐parasite allopatry/sympatry: one in which infections were considered “sympatric” if the host and parasite came from the same stream, and the other in which infections were considered ”sympatric” if the host and parasite came from the same island (i.e., Trinidad or Tobago). In exploratory analyses, we did not find a significant effect of parasite lineage within source location (Appendix S6).

As an indicator of possible relationships between population‐level genetic variability and worm‐days, we calculated the r ordered heterogeneity statistic (Rice & Gaines, 1994) for five population‐level genetic diversity metrics: phased and unphased diversity for MHC alleles and STs, and mean microsatellite diversity. Significant r values suggest an effect of a categorical predictor that is directional with respect to the ranking of categories by a third variable. Metrics were not derived from the present study's genotyped fish but from the much larger population genetics data set of Herdegen‐Radwan et al. (2020). We assess additional population variability metrics in Appendix S7.

Tests of individual‐level genetic predictors

Using all host populations, and restricting data to amplicons with ≥300 MHC sequence reads, we applied the same model ranking approach as above, with the same response variable and error distribution, to test for individual‐level effects of number of MHC alleles, number of MHC STs, and background (microsatellite) variability. Each genetic predictor was also tested for interactions with host population and parasite population. As our individual‐level multilocus microsatellite metric, we used 1 ‐ [homozygosity‐by‐loci] (henceforth “HL”). This weights loci by their expected heterozygosity (H E) when calculating multilocus heterozygosity, which may better capture background genetic variability when a microsatellite panel is small (Aparicio et al., 2006; in Appendix S8 we have repeated the analyses with alternative microsatellite metrics that apply different weightings, out of which HL produced the better fits). HL was calculated separately for each fish population, using allele frequencies from the experiment's genotyped fish. No population showed significant identity disequilibrium across the eight loci (g ≤ 0.07, p ≥ .40; Table S2.3; g tests performed in inbreedR (Stoffel et al., 2016), meaning there was not significant variance in individual hetero‐/homozygosity. Similarly, HL was not a significant predictor of MHC heterozygosity, although the relationship was positive (logistic regression: 0.37, SE = 0.28, z = 1.31, p = .19; note that we refer here to Mendelian heterozygosity of phased haplotypes, and not, as is often the case in the MHC literature, to the number of MHC variants that an individual carries).

A difference in effect size between MHC alleles and STs could be an artefact of aggregating a large number of alleles into a smaller number of STs. To assess this, we compared the observed coefficient for number of STs, taken from the best linear model to include the predictor, to a simulated distribution derived from repeating the model after reallocating alleles to STs at random (10,000 repeats; adapted from Lighten et al. (2017) and Herdegen‐Radwan et al. (2020)). We used the 335 alleles and 15 STs of Herdegen‐Radwan et al. (2020), and randomised ST membership size from a Dirichlet‐multinomial distribution (all α = 1; Appendix S9; R code included in data repository). To test whether the interpretation of HL was disproportionately influenced by any one microsatellite locus, we dropped each locus in turn, recalculated HL, and repeated the best model to include the predictor.

Tests of specific MHC STs

Finally, we tested for effects of specific supertypes on individual infection intensity. If a variant had at least three occurrences in more than one host population, we tested for an across‐population general effect, as well as the respective interaction. Population‐specific analyses were performed for all variants with 3+ instances in a population. If a variant had at least two occurrences for each gyrodactylid source within a host population, we tested the interaction. We did not analyse death rate with respect to individual MHC genotype, as only two dead fish gave MHC amplicons that met our genotyping quality criteria (details below).

Ethics statement

This experiment was conducted in accordance with Cardiff University's UK Home Office Licence PPL 303424. Tobago‐sourced wild fish were collected under Tobago House of Assembly Permit #004/2014. Permission to collect fish in Trinidad was granted by the Fisheries Division of the Ministry of Food Production and Fisheries, and fish were collected only from areas where guppies were previously identified as abundant.

RESULTS

We successfully infected 114 guppies with Gyrodactylus turnbulli, with 12–16 fish in each experimental block (four fish sources × two gyrodactylid sources; Figure 1). Overall mortality was low (five fish, 4.4%), and was entirely accounted for by fish from the Tobagonian populations of Scarborough Health Centre (“HC”; 2/32, 6.3%) and Roxborough (“Rox”; 3/27; 11.1%), with no fish dying from the Trinidadian populations Lopinot (“Lop”; 0/28) and Santa Cruz (“SC”; 0/27). Mortality was not significantly biased by either host population (bootstrap χ2: χ2 = 5.70, reps = 100 k, p = .12) or host island (Fisher's exact test: p = .06). Mortality by parasite source was evenly split (HC = 2/54, 3.7%; Lop = 3/60, 5.0%), and did not differ significantly across the two affected host populations (Fisher's exact test: p = 1).

FIGURE 1

Design of cross‐infection experiment, with four host (guppy) populations and two parasite (gyrodactylid) populations sourced from two different islands (Trinidad vs. Tobago). Sample sizes refer to numbers of infected hosts per treatment block, with each host receiving two gyrodactylid worms

The top‐ranked model of worm‐days among the 109 fish that survived the experiment, excluding genetic predictors, contained fish source, fish standard length, and an interaction between the two. Together, these explained 82% of variance (Tables 1, 2, 3). All pairwise differences between fish populations were significant in a post hoc test (t ≥ 4.65; p ≤ .001; ‘glht’ function of R package multcomp; Hothorn et al., 2008), with the rank order, in increasing infection intensity, of SC, Lop, HC, Rox (Figure 2a; Tables 2 and 3; infection trajectories in Appendix S11). Larger fish experienced significantly more worm‐days than smaller fish, though the slope differed significantly between populations. Although fish in SC were significantly smaller than the other three populations, multicollinearity was not problematic (Appendix S12). Of the two other models in the top two units of AICC (ΔAICC ≥ +1.34), both dropped the fish source × length interaction, retaining a significantly positive overall relationship between length and worm‐days (Tables 1 and 2). The second of these (ΔAICC = +1.95) added an interaction between fish and gyrodactylid source, but with a pattern not conforming to an obvious local adaptation scenario by either host or parasite (HC gyrodactylids produced heavier infections than Lop in SC and Rox fish, but lighter infections in Lop and HC fish; p ≥ .058; Table 2). We treat this latter interaction with extreme caution, however, as the model has very poor multicollinearity diagnostics (max. VIF = 19.0) and the interaction did not make the top set when the 12 fish without genotypes were excluded (see below and Appendix S14). There was poor support for a noninteracting effect of gyrodactylid source (ΔAICC ≥ +2.23), or for any effect of temperature (ΔAICC ≥ +4.35). Restricting local adaptation analyses to the two host sources from which parasites were collected (HC and Lop) changed the composition of the top set of models but did not produce a qualitative change in interpretation (Appendix S13).

TABLE 1

Selected models of the number of “worm‐days” experienced by guppies infected with Gyrodactylus turnbulli during our experiment, ranked by corrected Akaike Information Criterion (AICC)

Intercept	Fish source	Length	Gyro. source	Fish × length	Fish × gyro.	R 2	Res. skew.	Res. kurt.	Max. Cook	Max. VIF	Kappa	df	AICC	Weight
2.47	+	+	−	+	−	0.821	0.02	3.07	0.17	5.49	13.77	9	241.18	0.529
2.76	+	0.35	−	−	−	0.812	0.03	3.09	0.13	1.42	9.38	6	+1.34	0.271
2.59	+	0.35	+	−	+	0.820	0.01	3.05	0.11	18.98	20.67	10	+1.95	0.199
2.44	+	+	+	+	−	0.819	0.02	3.10	0.16	5.56	13.39	10	+2.23	‐
2.46	+	−	−	−	−	0.781	0.38	3.77	0.09	‐	5.07	5	+16.67	‐
4.36	−	0.61	−	−	−	0.129	0.44	2.40	0.06	‐	1.00	3	+164.93	‐
Sum of weights	1.00	1.00	0.20	0.53	0.20

Model ranking was conducted without genetic predictors and includes four fish source populations (see Table S13.1a for version restricting analysis to the two fish sources from which the experiment's gyros came). Selected models are: anything in the top two units of AICC (lowest AICC); the first model outside of this top set; the best model not to include length; and the best model not to include fish source. Models are general linear models (Gaussian error distribution) of log ‐transformed worm‐days. For continuous predictors, we give the regression slope when the parameter is present without an interaction, and “+” when it is part of an interaction; for categorical predictors and interactions, we indicate inclusion with “+”. Coefficients for fish length are for z‐transformed data. Importance is by Akaike weight and applies only to models in the top two units of AICC. “Res. skew.” and “Res. kurt.” are skewness and kurtosis of model residuals (Gaussian distribution has 0 and 3 respectively); “Max. Cook” is the maximum Cook's distance value, an indicator of over‐influential data points; “Max. VIF” (maximum variance inflation factor) and “Kappa” are multicollinearity diagnostics. Model‐averaged coefficients are in Table 2. Post hoc pairwise comparisons between all levels of fish source in the top‐ranked model are in Table 3

TABLE 2

Model‐averaged coefficients for models comprising the top two units of AICC model ranking of predictors of the number of worm‐days experienced by guppies infected with Gyrodactylus turnbulli during our experiment, in models without genetic predictors and including four fish source populations (Table 1)

Term	Slope	SE	z	p	Sum of weights
Intercept	2.57	0.22	11.73	<.001	–
Fish source
Lopinot	1.04	0.28	3.68	<.001	1
HealthCentre	2.20	0.27	8.28	<.001	1
Roxborough	3.92	0.29	13.46	<.001	1
Fish length (× fish source)
Fish length (no interaction)	0.35	0.08	4.27	<.001	0.471
Fish length (SantaCruz)	0.01	0.15	0.04	.969	0.529
× Lopinot	0.30	0.25	1.22	.225	0.529
× HealthCentre	0.53	0.20	2.70	.008	0.529
× Roxborough	0.49	0.24	2.02	.046	0.529
					Sum = 1
Gyro. Source × fish source
HealthCentre (SantaCruz fish)	0.39	0.27	1.43	.153	0.199
HealthCentre × Lopinot	−0.72	0.38	1.89	.058	0.199
HealthCentre × HealthCentre	−0.65	0.38	1.74	.082	0.199
HealthCentre × Roxborough	0.16	0.40	0.40	.686	0.199

Reference levels for fish and gyrodactylids are Santa Cruz and Lopinot respectively. Fish length (z‐transformed) was present in all models in the top set but interacted with fish source in 1/3 models, and we present this in the table as if length‐without‐interaction and length‐with‐interaction (the latter with Santa Cruz as the reference level) were two different predictors. Sum of weights (a.k.a. “importance”) is sum of Akaike weight. All p‐values are two‐tailed. Post hoc pairwise comparisons between all levels of fish source are in Table 3.

TABLE 3

Post hoc pairwise comparisons of number of worm‐days (log ‐transformed) between all levels of fish source in the top‐ranked model of Table 1

Term	Mean diff.	SE	t	p
Lopinot – SantaCruz	1.12	0.24	4.65	<.001
HealthCentre – SantaCruz	2.23	0.23	9.55	<.001
Roxborough – SantaCruz	4.08	0.23	17.39	<.001
HealthCentre – Lopinot	1.11	0.21	5.31	<.001
Roxborough – Lopinot	2.96	0.21	14.09	<.001
Roxborough – HealthCentre	1.85	0.20	9.19	<.001

FIGURE 2

Differences in worm‐days (log ) between (a) hosts (fish) from four populations, and (b) parasites (gyrodactylids) from two populations. Points are fitted means, errors bars are SEs for the means, and violins show partial residuals from the top‐ranked model to include the focal predictor. See Appendix S11 for average expected infection trajectories

All population‐level genetic diversity metrics had negative r values against population mean worm‐days, meaning populations with higher genetic diversity tended to have lower mean infection intensity (Table 4). The effect for microsatellite diversity was weaker than for STs (Table 4). Tests of additional genetic variability metrics are given in Appendix S7.

TABLE 4

Ordered heterogeneity testing by r (Rice & Gaines, 1994) of host population‐level genetic diversity metrics against population mean number of worm‐days

Metric	Marker class	MHC phasing	SantaCruz	Lopinot	HealthCentre	Roxborough	rSPc	p
ST div.	MHC ST	Unphased	0.84	0.83	0.72	0.14	−1.00	<.001
Msat div.(8)	Microsat.	‐	0.53	0.52	0.51	0.47	−1.00	<.001
S.hap. div.	MHC ST	Phased	0.84	0.84	0.72	0.14	−0.95	<.001
Haplo. div.	MHC allele	Phased	0.88	0.90	0.74	0.79	−0.60	.098
Allele div.	MHC allele	Unphased	0.89	0.91	0.74	0.79	−0.60	.098

r is calculated as the Spearman correlation coefficient (r) multiplied by 1 − (p‐value for the overall effect of host population). p‐values for r are derived from simulations (code in data repository). The host population columns (Santa Cruz‐Roxborough) are ordered to that mean worm‐days increases reading left to right (see also Tables 1 and 3). Additional metrics are tested in Appendix S7. All p‐values are two‐tailed.

Of the 109 fish that survived the experiment, 96 (88.1%) met our genotype quality criteria for inclusion in models of infection intensity. We recognise the theoretical possibility for the untyped fish to be a biased subset of genotypes, but the two supplementary analyses in Appendix S14 suggest that if such a bias exists, it is not problematic for our interpretation. The number of MHC STs was in all four models comprising the new top model set (Table 5). Its model‐averaged coefficient was significantly lower than zero (Tables 5 and 6), meaning individuals carrying more STs experienced infections of lower intensity. Microsatellite HL (our measure for overall genetic variation) was present in the third model (ΔAICC = +1.36; Table 5), with a nonsignificant negative coefficient (Table 6). Differences between populations in mean infection intensity were consistent with those of the nongenetic analysis (Tables 1 and 2; see also Appendix S14). Fish length retained its net positive relationship with worm‐days, but its interaction with fish source was no longer in the top‐ranked model (ΔAICC = +0.99; Table 5). Gyrodactylid source was not in the top set (ΔAICC = +2.31; Table 5), i.e. excluding it produced better fits. The first model to include number of MHC alleles had ΔAICC = +5.58.

TABLE 5

Models comprising the top two units of corrected Akaike Information Criterion (AICC) model ranking of predictors of the number of “worm‐days” experienced by guppies infected with Gyrodactylus turnbulli during our experiment, assessing models that were allowed to include an MHC and a microsatellite individual‐level genetic variability metric in addition to fish source (four‐level factor)

Intercept	Fish source	Msat HL	Length	N MHC als	N MHC STs	Gyro. source	Fish × length	Fish × N MHC STs	R2	Res. skew.	Res. kurt.	Max. Cook	Max. VIF	Kappa	df	AICC	Weight
2.83	+	–	0.42	–	−0.38	−	−	−	0.823	−0.28	2.92	0.12	1.73	8.63	7	208.31	0.399
2.58	+	–	+	–	−0.33	−	+	−	0.828	−0.28	2.80	0.21	5.66	12.78	10	+0.99	0.244
2.85	+	−0.07	0.43	–	−0.38	−	−	−	0.823	−0.26	2.93	0.11	1.79	8.85	8	+1.36	0.203
2.87	+	–	0.45	–	−0.63	−	−	+	0.827	−0.28	2.60	0.14	11.10	14.28	10	+1.89	0.155
2.81	+	–	0.42	–	−0.38	+	−	−	0.821	−0.28	2.91	0.11	1.75	8.24	8	+2.31	–
2.48	+	–	+	–	−	−	+	−	0.820	−0.23	2.76	0.25	5.35	12.96	9	+4.18	–
2.76	+	–	0.38	−0.19	−	−	−	−	0.812	−0.23	2.89	0.11	1.51	8.31	7	+5.58	–
Sum of weights	1.00	0.20	1.00	‐	1.00	‐	0.24	0.16

To provide additional context, we include: the first model outside the top set; the first model to include each genetic predictor not present in the top set; and the first model to include no genetic predictors. Models are general linear models (Gaussian errors) of log ‐transformed data. For continuous predictors, we give the regression slope when the parameter is present; for categorical predictors, we indicate presence with “+”. The coefficient for fish length is for z‐transformed data. Sum of weights (a.k.a. “importance”) is by Akaike weight and applies only to models in the top two units of AICC. Model‐averaged coefficients are in Table 6.

TABLE 6

Model‐averaged coefficients for models comprising the top two units of AICC model ranking of predictors of the number of worm‐days experienced by guppies infected with Gyrodactylus turnbulli during our experiment, in models that were allowed to include an MHC and a microsatellite individual‐level genetic variability metric in addition to fish source (Table 5)

Term	Slope	SE	z	p	Sum of weights
Intercept	2.78	0.20	13.56	<.001	–
Fish source
Lopinot	0.81	0.26	3.19	.001	1
HealthCentre	2.00	0.25	8.13	<.001	1
Roxborough	3.55	0.31	11.61	<.001	1
Microsats
Msat HL	–0.07	0.07	0.96	.335	0.203
MHC STs (× fish source)
N MHC STs (no interaction)	–0.37	0.15	2.49	.013	0.845
N MHC STs (SantaCruz)	–0.63	0.28	–2.22	.029	0.155
× Lopinot	0.19	0.35	0.55	.581	0.155
× HealthCentre	0.77	0.40	1.93	.057	0.155
× Roxborough	–0.11	0.47	–0.23	.822	0.155
					Sum = 1
Fish length (× fish source)
Fish length (no interaction)	0.43	0.09	4.97	<.001	0.756
Fish length (SantaCruz)	0.10	0.16	0.65	.520	0.244
× Lopinot	0.26	0.25	1.06	.291	0.244
× HealthCentre	0.44	0.20	2.21	.030	0.244
× Roxborough	0.49	0.25	1.95	.054	0.244
					Sum = 1

Models are general linear models (Gaussian error distribution) of log ‐transformed worm‐days. Fish source coefficients are given in reference to Santa Cruz. HL coefficient is for z‐transformed data. Fish length (z‐transformed) was present in all models in the top set but interacted with fish source in one quarter models, and we present this in the table as if length‐without‐interaction and length‐with‐interaction (the latter with Santa Cruz as the reference level) were two different predictors. The same approach was used for the interaction between fish source and number of MHC STs. Sum of weights (a.k.a. “importance”) is sum of Akaike weight. All p‐values are two‐tailed.

The stronger effect of number of STs relative to number of alleles was significantly greater than expected from randomised grouping of alleles (obs. coef. = –0.38, exp. = –0.13 ± 0.09 [SD], n. randomisations = 10 k, p = .003; details in Appendix S9), meaning it is unlikely to be a side‐effect of aggregating alleles into STs. No microsatellite locus exerted a disproportionate influence on HL in jackknife removal (Appendix S10).

Five supertypes were carried by 3+ fish in at least two populations (Tables 7 and 8) and were thus available for cross‐population testing for effects of specific STs on infection intensity. Three of these STs were in their respective top model set, of which two were in their top‐ranked model, of which one was present in all top‐set models and/or had a significant coefficient in its best model (Table 8; see also Appendix S15). Carriers of ST12 (SC and Lop) experienced significantly fewer worm‐days than noncarriers, with an effect size comparable to, and not confounded by, adding an extra supertype (Table S15.7).

TABLE 7

Per‐population counts of carriers of each MHC supertype (ST) among genotyped fish in our experiment. ST groupings are those of Herdegen‐Radwan et al. (2020), and STs 4, 8, and 9 were not detected in the present experiment's sample

MHC ST	Santa Cruz	Lopinot	Health Centre	Roxborough
ST01	2	5b	0	0
ST02	1	10b	9a	0
ST03	1	10a	0	0
ST05	0	0	18b	0
ST06	6b	2	0	0
ST07	13b	5b	5a	0
ST10	3b	11b	0	4d
ST11	1	1	0	0
ST12	9a	4a	0	0
ST13	2	0	0	0
ST14	3a	0	0	0
ST15	0	2	16b	20c,d

Instances of STs with 3+ carriers in a given population were tested for within‐population resistance/susceptibility effects: aST with an effect present in the top two units of its population's AICC‐ranked general linear models of worm‐days (see Table 9 for details of such instances); bST not present in its population's top model set; cST carried by every genotyped individual in a population; dPerfectly confounded with another genetic predictor. Tests of supertypes with 3+ carriers in 2+ populations are in Table 8.

TABLE 8

Descriptions of effects associated with single MHC supertypes (STs) in multiple populations (3+ carriers in 2+ populations) in AICC‐ranked general linear models of worm‐days

MHC ST	Populations	Res./susc.	First model	Sum of weights	Min. p	Interactions	See also
ST02	Lop, HC	Resistant	Top‐ranked model	0.390	.132	None	Table S15.4
ST07	SC, Lop, HC	Interaction	+0.03	0.590	.136	Some	Table S15.5
ST10	SC, Lop	–	+2.56	‐	–	–	Table S15.6
ST10	SC, Lop, Rox	–	–	‐	–	–	Appendix S15
ST12	SC, Lop	Resistant	Top‐ranked model	1.000	.023	None	Table S15.7
ST15	HC, Rox	–	–	–	–	–	Appendix S15

Res./susc. = whether carrying the ST is associated with reduced or increased infection intensities (“resistant” or “susceptible”) among the top two AICC units of a focal ST’s ranked models. Numerical values in “First model” are ΔAICC values relative to the top‐ranked model. Sum of weights, sum of Akaike weights of models containing the ST within the focal top model set (a.k.a. “importance”). Min. p, lowest p‐value for the ST’s effect among the top model set. “Interactions” indicates whether there are interactions between the ST and fish/gyro source. Tests of specific STs in population “Rox” are perfectly confounded with number of STs carried. See Table S15.2 for remarks on each test.

TABLE 9

Descriptions of single MHC supertypes (STs) with effects present in the top two units of a single population's AICC‐ranked general linear models of worm‐days

MHC ST	Populations	Res./susc.	First model	Sum of weights	Min. p	Interaction	See also
ST02	HC	Resistant	+1.61	0.174	.320	None	Table S15.11
ST03	Lop	Susceptible	+1.52	0.262	.118	None	Table S15.10
ST07	HC	Susceptible	+1.12	0.221	.234	None	Table S15.11
ST12	SC	Resistant	+1.92	0.089	.211	None	Table S15.9
ST12	Lop	Interaction	+1.26	0.293	.023	Some	Table S15.10
ST14	SC	Resistant	Top‐ranked model	0.487	.070	None	Table S15.9

An ST required 3+ carriers to be tested. “Res./susc.”, whether carrying the ST is associated with reduced or increased infection intensities (“resistant” or “susceptible”). Numerical values in “First model” are ΔAICC values relative to the top‐ranked model. Sum of weights, sum of Akaike weights of models containing the ST within the focal top model set (a.k.a. “importance”). Min. p, lowest p‐value for the ST’s effect among the top model set. “Interaction” indicates whether the ST interacts with gyrodactylid source. See Table 7 for other tested STs, and Table S15.3 for remarks on each test.

For population‐specific tests for effects of particular supertypes, there were 17 testable instances (3+ carriers in a population; 2–6 STs per population; Table 8). Of these, six were in their population's top model set (Table 9), of which only one was present in all models of a top set or had a significant coefficient. Lop hosts carrying ST12 experienced significantly fewer worm‐days, but only when infected with Lop (i.e., local) gyrodactylids. However, this result should be treated with caution, as the interaction is based on the minimum allowable sample size of two carriers per gyrodactylid source. Tests of specific STs in Rox hosts were perfectly confounded with number of STs – only two STs were detected in the genotyped fish, of which one (ST15) was carried by all individuals.

DISCUSSION

In our controlled parasite infection trials, we found strong differences between four host populations but no difference between two parasite populations. There were no significant host × parasite interactions, implying no significant local adaptation. Host populations with higher genetic variability (MHC supertypes [STs] and microsatellites) experienced gyrodactylid infections of lower intensity. There was a strong, consistent, negative relationship between individual MHC ST variability and infection intensity that could not be accounted for by microsatellite variability. This relationship was substantially stronger for MHC STs than for MHC alleles (nucleotide sequences). One ST showed a particularly strong effect, though the statistical power for such tests was constrained by low sample sizes.

The lack of evidence for local host or parasite adaptation contrasts with evidence for local adaptation to parasites in other fish systems (Bolnick & Stutz, 2017; El Nagar & MacColl, 2016), but is not out of keeping with the multitaxa meta‐analysis of Greischar and Koskella (2007) in which studies reporting no significant overall local adaptation slightly outnumber (56%) those that do. However, host population, but not parasite population, was a strong, significant predictor of the intensity of infection experienced by hosts. Comparable results of consistent differences between host populations have been reported in other fish species (e.g., Konijnendijk et al., 2013; Pérez‐Jvostov et al., 2015; Weber et al., 2017), as well as in laboratory guppies infected with the same parasite species (Smallbone et al., 2021). The precise reasons for such differences are usually unclear, though populations of stickleback (Gasterosteus aculeatus) have been shown to differ in gene expression profiles in response to Gyrodactylus spp. infection (Brunner et al., 2017; Robertson et al., 2017). In our study, infection intensity was associated with both population‐ and individual‐level genetic variability of hosts. Contrasts between the effects of MHC physicochemical STs, MHC alleles, and microsatellites indicate a functional role of the MHC in explaining these patterns, whilst also suggesting the MHC is not the whole story.

Higher host population genetic diversity was associated with lower population mean infection intensity. These associations were significant for MHC STs and microsatellites, but not for MHC alleles. Several previous studies have reported evidence for positive correlations between pathogen diversity and MHC polymorphism and/or positive selection in cyprinids (Šimková et al., 2006), birds (Minias et al., 2018), rodents (de Bellocq et al., 2008), and primates (Garamszegi & Nunn, 2011). Intraspecies, interpopulation examples include MHC polymorphism correlating positively with pathogen vector and ectoparasite abundance in the lizard Ctenophorus ornatus (see Radwan et al., 2014), and with length of time in which rabies has been present in racoon populations (Procyon lotor) (see Kyle et al., 2014). However, these examples should be compared to our study with caution, as all are observational/ecological rather than exposure‐controlled experiments, and all focus on MHC alleles without a parallel assessment of STs. The significant relationship for MHC STs, which agrees with another recent guppy‐gyrodactylid‐MHC study (Smallbone et al., 2021), suggests a history of stronger selection on functional MHC. Our observed relationship could also result from population demographic histories affecting diversity by drift, with which the significant alignment with microsatellite diversity would be consistent. Drift and selection, though, are not mutually exclusive (e.g., C. ornatus; see Radwan et al., 2014), and the much wider range of ST diversity relative to microsatellite diversity (0.14–0.84 vs. 0.47–0.53, respectively) is hard to reconcile with drift being dominant over selection on MHC functionality. Moreover, MHC allele diversity (range 0.74–0.91), which is known to be sensitive to drift (Herdegen‐Radwan et al., 2020; Lighten et al., 2017; McMullan & van Oosterhout, 2012; Radwan et al., 2010), was not significantly aligned with infection load.

Further evidence supporting the role of functional MHC diversity comes from our analyses of individual infection intensities. Stronger selection by parasites should not only maintain more variants in a population, but should also select for more variants expressed by an individual. The latter will be determined by both zygosity and the number of MHC genes in haplotypes (Bentkowski & Radwan, 2019; Minias et al., 2015), as more variants should widen the range of antigens that can be detected and responded to. Such effects have been reported in numerous other studies (Carrington et al., 1999; Oliver et al., 2009; Penn et al., 2002; Pierini & Lenz, 2018; Radwan et al., 2012), but others have found no such effect (Phillips et al., 2018), some have found the reverse (Ilmonen et al., 2007; Schwensow et al., 2007; Sepil et al., 2013), and, in species with high numbers of duplicated MHC loci, some have found intermediate numbers of MHC alleles to be associated with the lowest infection burdens (Kloch et al., 2010; Madsen & Ujvari, 2006; Wegner et al., 2003, 2003, 2008). In our experiment, individuals carrying more MHC STs experienced significantly fewer worm‐days: back‐transformed from the log scale (Table 6), one extra ST predicts a reduction in worm‐days over the study period by 31.4%. This effect was not confounded by differences in ST diversity between host populations – it applied within each population, and without a significant interaction. Moreover, it did not interact with parasite population, and, importantly, was independent of multilocus neutral (microsatellite) heterozygosity. In contrast, the effect of number of MHC alleles was relatively weak (17.2% worm‐day reduction, and nonsignificant), producing a poorer fit than a model with no genetic predictors. The ST effect is unlikely to be a consequence of aggregating alleles into STs, as the coefficient is significantly more negative than expected if alleles are clustered into STs at random. Collectively, this suggests that the relationship between the number of STs and infection intensity may be causal rather than correlational. If gyrodactylid infections select for individual‐level functional MHC variability in guppies, our analysis of individual‐level predictors supports a role for differences in past selection by parasites in causing differences between host populations in susceptibility to infection.

The controlled infection experiments of Smallbone et al. (2021), Phillips et al. (2018), and the present study highlight the need for better characterisation of the mechanistic role of MHC Class II in rapid primary immune responses to ectoparasites. Konczal et al. (2020) highlight the role of proinflammatory, MHC‐dependent Th17 pathway in response to gyrodactylid infection. Additionally, in vitro work on zebrafish skin tissue has suggested a role of antigen‐presenting skin cells in mediating production of proinflammatory cytokines (Lugo‐Villarino et al., 2010), and there is implied support for this from gene expression work on the salmonid response to ectoparasitic lice (e.g., Braden et al., 2015). Gyrodactylid infection is unlikely to provoke major upregulation of MHC Class II: it was not observed to do so in salmonids (Jørgensen et al., 2009) or guppies (Konczal et al., 2020), and its expression did not differ significantly between primary and secondary infection in goldfish (Zhou, Li, et al., 2018). However, Konczal et al. (2020) note its high constitutive expression, and its involvement in a gene co‐expression module that was itself correlated with gyrodactylid burden within 4–8 days of infection. The guppy‐gyrodactylid experiments to date cannot completely rule out a linkage effect driving the observed MHC patterns, but it would need to be close, physical linkage, which thus excludes MHC Class I in teleosts (Sato et al., 2000), in order to pass the breeding design of Phillips et al. (2018). Moreover, observing the strongest effects to be associated with physicochemical supertypes, which tend to be polyphyletic (e.g., Herdegen‐Radwan et al., 2020), argues for the MHC itself being responsible.

Interestingly, microsatellite variability itself also showed a weak negative effect in the top set of models with genetic predictors. Importantly, the addition of microsatellite variability neither weakened the effect of number of STs nor produced problematic multicollinearity. As with number of STs, the effect was not confounded by between‐population differences in variability. This suggests a general heterozygosity‐fitness correlation (HFC), of which there are numerous examples relating to infection susceptibility in both genome‐wide terms (e.g., Acevedo‐Whitehouse et al., 2003; Eastwood et al., 2017; Luikart et al., 2008) and with respect to other families of immune system genes (e.g., Hellgren et al., 2010; Lara et al., 2020; Levy et al., 2020). The relative weakness of this effect (an entirely heterozygous individual would experience 1.2% fewer worm‐days than an entirely homozygous individual) may be due to limited within‐population variance in genome‐wide heterozygosity (nonsignificant identity disequilibrium; David et al., 2007) coupled with a small number of loci. However, neither of these caveats justifies outright dismissal of our observed effect: although g correlates with HFC effect size (e.g., Miller & Coltman, 2014), HFCs can reach significance before g does (Szulkin et al., 2010), and weak microsatellite HFCs can hint at effects that become much stronger when large panels of neutral SNPs are available (e.g., Hoffman et al., 2014).

Santa Cruz and Lopinot hosts carrying ST12 experienced significantly fewer worm days than noncarriers. This effect was not confounded by number of STs carried by a host, suggesting ST12 may be a “resistant” ST. This should be treated with caution, however, as the number of Lopinot carriers is low (n = 4) and there is a suggestion that ST12 may interact with gyrodactylid source within this host population (Table S15.10). Indeed, the general weakness of evidence for effects of specific STs despite a strong effect of the number of MHC STs may be a result of our limited power to detect the former. Few MHC variants were shared among populations at high enough frequency to carry out meaningful analyses against the relatively small sample sizes (n = 12–16) of each host × parasite treatment block. Previous work has shown that the guppy‐gyrodactylid system is capable of producing such effects, but also that those effects are context‐dependent (Smallbone et al., 2021). If ST12 is taken as “resistant” in our study, this would underscore this context dependence, as it is not one of the STs identified as having a strong effect in Smallbone et al. (2021).

The lack of a significant effect of parasite source population, or of any interaction with host population, appears inconsistent with previous work on this study system showing that MHC variants that were novel to a given host population conferred a significant advantage in resisting local parasites (Phillips et al., 2018). That finding, which used replicated population crosses to control for non‐MHC genetic background, and which did not find significant interactions associated with those crosses, hinges on parasites being adapted to their local host's MHC composition. Host populations probably segregate at other gene families reported to be involved in fish immune responses to ectoparasites (Konczal et al., 2020; Lindenstrøm et al., 2004; Skugor et al., 2008; Zhou et al., 2018), and effects of these loci may mask any effect of MHC novelty in the present study. The results of a recent translocation study on sticklebacks imply such background effects (Bolnick & Stutz, 2017). While MHC alleles that were common in a local host population were more susceptible to local parasites (consistent with an advantage from introgressing MHC alleles), immigrant fish experienced higher parasite loads when MHC genotype was controlled for (Bolnick & Stutz, 2017). The lack of signal of local adaptation in our study may thus result from opposing signals of MHC and other genomic regions affecting resistance to gyrodactylids.

An alternative explanation for the lack of a significant effect of parasite source population in our study is the culturing of parasites (“gyro farming”). HC and Lop gyrodactylids were both cultured on a single, separate lineage of hosts for 18 days (3–4 farm cycles; potentially 9+ parasite generations) in order to obtain sufficient numbers for the experimental infections. This may have caused artificial selection, either on ability to infect novel hosts or to a specific set of host immune genotypes. However, we think this explanation is unlikely, as it would require the effective erasure of many generations of local adaptation in two large, natural populations of parasites, and it would require that this be achieved in a short time, with limited starting genetic variance (each culture lineage was founded from a single animal), on parasite‐naïve hosts.

Host mortality in the experiment was low and restricted to the two Tobago fish populations, and was not significantly biased by gyrodactylid origin. Low host mortality probably reflects our experiment's relatively benign conditions. In the wild, gyrodactylid infections expose guppies to a suite of additional selection pressures, including other infections, anything requiring efficient swimming (e.g., escape from predators and surviving flood events; van Oosterhout et al., 2007; Stephenson et al., 2016), and reduced reproductive opportunities (Kennedy et al., 1987).

Overall, our results suggest that pathogens may select for higher numbers of MHC supertypes at an individual level, and previous work implies that it can also promote MHC polymorphism within populations, independently of benefits derived from simply carrying more variants (Phillips et al., 2018). Both of these effects will lead to differences between populations in functional immunogenetic diversity. The wider genomic implications of such selection (e.g., interactions with other immunity genes, effects on neutral genetic variability, and the shaping of MHC phylogenetics) requires further investigation. While this is a considerable effort, our results highlight that our understanding of infection dynamics will remain incomplete unless we appreciate the differences in the history of selection imposed by pathogens. A lack of such understanding may limit our ability to predict consequences of emergent diseased threatening humans and wildlife (Altizer et al., 2003; Chabas et al., 2018; Ekroth et al., 2019; Penczykowski et al., 2016; Stephenson et al., 2017), and further research in this area should underpin the One Health approach (Daszak et al., 2000; van Oosterhout, 2021) in the coming decades.

CONFLICT OF INTEREST

The authors have no conflict of interest to declare.

AUTHOR CONTRIBUTIONS

Jacek Radwan, Cock van Oosterhout and Joanne Cable conceived the study, with input from Karl P. Phillips and Ryan S. Mohammed on experimental design. Jacek Radwan and Ryan S. Mohammed collected the fish, which Ryan S. Mohammed reared; Ryan S. Mohammed and Karl P. Phillips collected and reared parasites; Karl P. Phillips and Sebastian Chmielewski performed infection trials, with assistance from Karolina J. Przesmycka and Ryan S. Mohammed; Sebastian Chmielewski and Karolina J. Przesmycka assisted molecular work; Karl P. Phillips analysed the data; Jacek Radwan and Karl P. Phillips wrote the manuscript, with input from all coauthors.

Supinfo

Click here for additional data file.