Long term patterns of association between MHC and helminth burdens in the bank vole support Red Queen dynamics.

Major histocompatibility complex (MHC) genes encode proteins crucial for adaptive immunity of vertebrates. Negative frequency-dependent selection (NFDS), resulting from adaptation of parasites to common MHC types, has been hypothesized to maintain high, functionally relevant polymorphism of MHC, but demonstration of this relationship has remained elusive. In particular, differentiation of NFDS from fluctuating selection, resulting from changes in parasite communities in time and space (FS), has proved difficult in short-term studies. Here, we used temporal data, accumulated through long-term monitoring of helminths infecting bank voles (Myodes glareolus), to test specific predictions of NFDS on MHC class II. Data were collected in three, moderately genetically differentiated subpopulations in Poland, which were characterized by some stable spatiotemporal helminth communities but also events indicating introduction of new species and loss of others. We found a complex association between individual MHC diversity and species richness, where intermediate numbers of DRB supertypes correlated with lowest species richness, but the opposite was true for DQB supertypes-arguing against universal selection for immunogenetic optimality. We also showed that particular MHC supertypes explain a portion of the variance in prevalence and abundance of helminths, but this effect was subpopulation-specific, which is consistent with both NFDS and FS. Finally, in line with NFDS, we found that certain helminths that have recently colonized or spread in a given subpopulation, more frequently or intensely infected voles with MHC supertypes that have been common in the recent past. Overall, our results highlight complex spatial and temporal patterns of MHC-parasite associations, the latter being consistent with Red Queen coevolutionary dynamics.

INTRODUCTION

Genes of the major histocompatibility complex (MHC) of vertebrates code for transmembrane proteins crucial for molecular‐level self/nonself recognition and necessary for the initiation of adaptive, T cell‐mediated immunity. They are highly polymorphic, comprising dozens to hundreds of alleles which have evolved under selection favouring substitutions in amino‐acid positions involved in antigen binding (Kindt et al., 2007). MHC alleles are associated with resistance or susceptibility to disease (Matzaraki et al.,2017) and consequently MHC polymorphism (or its loss) has been implicated in determining the dynamics of epidemics (Alves et al., 2019), population viability (O’Brien et al., 1985; Radwan et al., 2010), as well as in driving sexual selection (Milinski, 2006) and speciation (Andreou et al., 2017; Eizaguirre et al.,2009). Despite its wide ecological and evolutionary importance, the evolutionary forces maintaining MHC polymorphism are still not fully understood (Radwan et al., 2020; Spurgin & Richardson, 2010).

Explanations focus on three, nonmutually exclusive mechanisms. The first one is heterozygote advantage (HA), which poses that having two distinct, codominantly expressed alleles per MHC locus will allow recognition of a broader spectrum of diverse parasitic species/genotypes (Doherty & Zinkernagel, 1975). The second mechanism that could maintain MHC polymorphism is negative frequency‐dependent selection (NFDS) resulting from the Red Queen dynamics of parasite and host adaptations and counteradaptations (Bodmer, 1972). Parasites are predicted to adapt to the most common MHC alleles in a given population, in a way that will allow them to escape recognition by the host immune system. Thus, having a rare MHC allele could be advantageous to its carrier, as parasites are less likely to be adapted to evading it. The third proposed mechanism rests on selection varying in time and space, which takes into account the fact that parasitic pressure on hosts can vary both spatially and temporally. The resulting fluctuating selection (FS) could promote polymorphism by favouring different MHC alleles in different geographic regions and at different times (Hedrick, 2002). The degree to which each of these mechanisms contributes to selection on MHC in natural populations is a subject of an ongoing debate (Radwan et al., 2020).

HA has been supported by experimental infection studies (Penn et al., 2002) (but see Ilmonen et al., 2007) and by studies examining associations with pathogen burdens observed in natural populations (Arora et al., 2020; Oliver et al., 2009; Osborne et al., 2015) but not universally so (for exceptions see for example: Dionne et al., 2009; Peng et al., 2021). Even when HA is observed, in natural populations it is in general difficult to discern the effect of heterozygosity alone, because rare (and possibly advantageous) alleles from the NFDS scenario are typically found in heterozygotes. Also, the predictions of the HA hypothesis are complicated by the fact that if a wider breadth of peptide‐binding is advantageous, natural selection should favour duplication and divergence of MHC loci. However, expansion of the MHC region might come at a certain cost, related to for example, tradeoffs associated with auto‐immunity. Hence, according to the optimality hypothesis, an optimal, rather than maximal individual MHC diversity may be favoured by selection (Nowak et al., 1992; Woelfing et al., 2009). Such nonlinear associations between MHC variant number and pathogen burdens have been observed in some studies (e.g., Wegner et al., 2003), but not in others (e.g., Radwan et al., 2012).

The role of spatiotemporal variation is less well documented. Parasite communities infecting vertebrate species are likely to differ both in time and space, with different parasites favouring different MHC alleles (Osborne et al., 2017). Gene flow between populations may then contribute to MHC population‐level diversity, at the same time preventing divergence between them (Bolnick & Stutz, 2017). NFDS has proven to be even more difficult to demonstrate, despite theoretical support (Borghans et al., 2004; Ejsmond & Radwan, 2015) and experimental endorsement of mechanistic prerequisites, for example, specific adaptation of pathogens to encountered MHC genotypes (Kubinak et al., 2012), selective advantage of rare or novel genotypes (Bolnick & Stutz, 2017; Phillips et al., 2018), sufficient fitness costs of carrying parasites that could drive change in allele frequencies (i.e., decreasing the susceptible and increasing the resistant ones) (Eizaguirre et al., 2012). Yet, compelling evidence for these processes leading to operation of NFDS in nature is scarce. NFDS has been mostly tested in snapshot studies of host and parasite populations, by looking at the association between infection status and allele frequency (Schad et al., 2005). However, the power of such tests is compromised, because adaptation requires time, resulting in time‐lags in host‐parasite coevolution (Dybdahl & Lively, 1998). As a consequence, a currently rare MHC allele that was common in the recent past can still be susceptible to infection and conversely, pathogens, particularly slowly‐evolving ones, might have not yet been adapted to a currently common allele if it has increased in frequency only recently. It is therefore not surprising that the results of such tests were often negative or conflicting (e.g., Jones et al., 2015; Kamath et al., 2014; Peng et al., 2021). Furthermore, associations of particular MHC alleles (but not their heterozygosity) with infection status are compatible with NFDS, but also with FS (Spurgin & Richardson, 2010) and both NFDS and FS can explain the pronounced temporal fluctuations in MHC allele frequency reported in some studies (Westerdahl et al., 2004).

Use of temporal data on MHC offers a more powerful test of the NFDS. For example, recently, Huang et al. (2022) found in a feral population of Soay sheep, monitored for over 20 years, that an initially rare haplotype conferred increased male total fitness and increased in frequency over consecutive seasons (although this pattern could also be consistent with directional selection). Temporal data on MHC‐infection associations also allow a test of time ×MHC type interactions, the presence of which would imply a dynamic nature of selection on MHC due to evolution of parasites, supporting the Red Queen process underlying NFDS. Furthermore, it allows testing of whether parasites are more successful in infecting host genotypes that were common in the recent past, thus accounting for a time required for parasites to adapt. Finally, temporal data occasionally reveal novel/colonizing parasites spreading in given host populations, as well as loss/extinction of species over time. Such events offer a strong prediction about frequency‐dependent parasite adaptation: host genotypes which were common at the time of colonization by a particular parasite species should be more prone to infection after the parasite has spread.

In this study, we investigated spatiotemporal patterns of associations between MHC and parasite infections using data accumulated through long‐term monitoring of helminth communities affecting bank voles (Myodes glareolus) in the Mazury Lake District Region in northeastern Poland (Behnke et al., 2001; Grzybek, Bajer, Bednarska, et al., 2015). Parasitological surveys were conducted in three sites moderately differentiated genetically (the microsatellite‐based F ST reported earlier was 0.066, Kloch et al., 2010). Previous snap‐shot analysis carried out within a single season revealed associations between MHC DRB variant and the oxyuroid nematode Aspiculuris tetraptera occurring in these populations (now referred to as A. tianjinensis, Behnke et al., 2015; Grzybek, Bajer, Bednarska, et al., 2015), as well as a nonlinear association of helminth species richness and the number of MHC DRB variants (Kloch et al., 2010). Extended monitoring has revealed that despite overall environmental similarity and some shared helminth species, these bank vole subpopulations have also maintained relatively stable differences in the composition of their helminth communities (Bajer et al., 2005; Behnke et al., 2001; Grzybek, Bajer, Bednarska, et al., 2015; Grzybek et al., 2015). This has allowed us to test the prediction of FS that different MHC alleles will be favoured in different populations due to their associations with resistance to spatially/temporally variable parasite species. Furthermore, temporal aspect of our data allowed us to test the prediction of NFDS that MHC variants that have been more common in the recent past will be associated with higher susceptibility to infection. Finally, during the monitoring period a new helminth species, Aonchotheca annulosa, appeared in our study population. Since its first record in 2002, it has risen steadily in prevalence and abundance in two of the studied sites (Grzybek, Bajer, Bednarska, et al., 2015). This allowed us to test the predictions of NFDS in a particularly powerful way, as the signal of adaptation to the most common MHC types should not be confounded by adaptation to this parasite in the past. As in the previous study, we focus here on MHC class II, because it presents peptide fragments from extracellular space, derived from for example, extracellular bacteria or macroparasites, such as helminths. In addition to the DRB genes studied by Kloch et al. (2010), we included DQB genes, which also code for classical MHC class II molecules, and in bank voles show stronger signatures of selection for amino‐acid substitutions compared to DRB (Migalska et al., 2019; Scherman et al., 2014). We explored the associations between bank vole MHC class II variation at both DQB and DRB and species richness, prevalence and abundance of helminths over 11 years of sampling. In addition to associations between MHC alleles and infection, Kloch et al. (2010) found that consistent with the optimality hypothesis, individuals with an intermediate MHC diversity over several duplicated DRB loci were infected with the fewest parasites. Therefore, in addition to testing predictions concerning spatial and temporal associations of particular alleles, we also test whether the optimality effect also holds across both DRB and DQB.

MATERIALS AND METHODS

Study species and study sites

The study species was the bank vole, one of the most widespread rodent species in European forests (Wilson & Reeder, 2005) utilized as a model species in numerous ecological and evolutionary studies (Larsen, 2016). The study was approved by the First Warsaw Local Ethics Committee for Animal Experimentation, which also has an overarching responsibility for fieldwork involving the trapping and culling of wild vertebrates for scientific purposes (licence numbers: 280/2003, 737/2007, 73/2010, 578/2014). The bank voles were sampled in mid‐August to mid‐September, over 11 years, at 3–4‐year intervals (i.e., in years 1999, 2002, 2006 and 2010). The methods used for trapping rodents, sampling and processing trapped animals have been thoroughly described (Behnke et al., 2001, 2008). Briefly, rodents were caught live in wooden traps. Trapping was carried out for 3–4 consecutive days at a time at each site. Trapped animals were brought to the University of Warsaw Field Station in Urwitałt for further examination. Three age classes were established (class 1— immature juvenile bank voles; class 2—mostly young adult bank voles; and class 3—breeding older animals), with use of principal components analysis (PCA) of a range of morphological measures, according to the method described in detail by Behnke et al. (2001) and Grzybek, Bajer, Bednarska, et al. (2015).

The three study sites, corresponding to three bank vole subpopulations, have already been described in detail by Behnke et al. (2001, 2008). They are located in northeastern Poland, in the Mazury Lake District. Site 1 is referred to as Urwitałt (N 53°48.153, EO 21° 39.784), Site 2 as Tałty (N 53°53.644, EO 21°33.049), and Site 3 as Pilchy (N 53°42.228, EO 21°48.499), after nearby settlements. Sites are located along NW‐SE transect (Tałty – Urwitalt ‐ Pilchy), spaced approximately 10 km apart, and separated by lakes, rivers, canals and pastures limiting vole migration (Figure S1).

Parasitological data

The methods used for identification and quantification of endoparasites, helminth community structure, prevalence and abundance of particular species, as well as comprehensive statistical analyses of the effects of intrinsic factors (e.g., sex) have been reported in Grzybek, Bajer, Bednarska, et al. (2015). Our current study focused on the associations between these helminth communities and MHC class II diversity. Except for the overall test of association between the individual number of DQB and DRB supertypes and helminth species richness (a total number of helminth species in a given individual), all the following analyses were conducted for the six most prevalent helminths over the course of these 11 years. This includes five nematode species: Heligmosomum mixtum, Heligmosomoides glareoli, Mastophorus muris, Aspiculuris tianjinensis (previously identified in bank voles as A. tetraptera), Aonchotheca annulosa and one cestode: Catenotaenia henttoneni. The animals from the three study sites were characterized by quite distinct combinations of infecting helminths. Detailed description and analyses of these communities is available in Grzybek, Bajer, Bednarska, et al. (2015), and a summary of spatiotemporal differences in parasite abundance and prevalence among sites (based on the subset of individuals analysed in the current study) is in Figure S2.

Bank vole MHC class II

Several studies to date have examined allelic diversity, the extent of gene duplication and the signatures of selective pressure acting on MHC class II genes of bank voles (Bryja et al., 2006; Migalska et al., 2019; Scherman et al., 2014). Overall, bank vole MHC class II is characterized by high allelic polymorphism, extensive duplication and intraspecific variation in the number of loci. While DP genes have not been identified in bank voles and DRA are close to monomorphic, DQA, DQB and DRB seem to be at least quadruplicated in some haplotypes (Axtner & Sommer, 2007; Bryja et al., 2006; Migalska et al., 2019; Scherman et al., 2014). Interestingly, unlike in many wild mammals, there are clusters of bank vole DRB genes that lack the hallmarks of classical, mammalian MHC class II genes, that is, positive selection in peptide binding sites and trans‐species polymorphism. In particular, Axtner and Sommer (2007) found two distinct, well‐defined gene clusters (“Cluster I” and “Cluster II”) that were characterized by negative (purifying) selection. Other DRB alleles (termed “Group III”) and all genes in the DQ group (both DQA and DQB) bear the characteristics of functional MHC class II molecules, with several codons shown to be under positive selection (inferred from high rates of nonsynonymous over synonymous substitutions), and widespread signals of trans‐species polymorphism (Axtner & Sommer, 2007; Bryja et al., 2006; Migalska et al., 2019; Scherman et al., 2014). In the current study, we chose to focus on DRB alleles from group III only and DQB alleles (see details in Supporting Information, Methods).

We used previously published primers for amplification of the second exons of DRB and DQB genes of bank voles, as these encode domains that contribute to the peptide binding groves of the MHC class II molecules and determine antigen‐binding specificity (Brown et al., 1993). Primers were slightly adjusted to account for recently acquired data on the diversity of expressed MHC class II molecules in this species (Supporting Information, Methods, Table S1). DNA from 1,022 bank vole tail clips, preserved in 20% DMSO in saturated saline, was extracted with MagJET genomic DNA kit (Thermo Fisher Scientific) and approximately 20–50 ng were used as template in subsequent PCR reactions. Details on PCR cycling and multiplexed library preparation for high‐throughput sequencing are in Supporting Information Methods. All samples were sequenced on MiSeq (Illumina) with MiSeq Reagent kit v2 (paired‐end, 2×150 cycles) in a single sequencing run, at the Institute of Environmental Sciences, Jagiellonian University, Kraków, Poland.

MHC genotyping was performed with an adjustable clustering method (Biedrzycka et al., 2017), using AmpliSAS software (Sebastian et al., 2016) for demultiplexing, clustering and filtering of reads. Details of the parameters used, repeatability and assessment of contamination are given in Supporting Information, Methods. After genotyping, DNA sequences were translated into protein sequences. Before downstream analyses, sequences with stop codons (pseudogenes) were removed, and variants coding for identical protein sequences were collapsed into one.

Next, protein variants were assigned to supertypes – groups clustering variants according to the physicochemical properties of peptide binding residues (PBRs) (Doytchinova & Flower, 2005). The MHC crystal structures have not been resolved in this species, and the exact residues involved in antigen binding are unknown, so we used positively selected sites (PSSs) as a proxy of PBRs. PSSs are identified as sites with an increased rate of nonsynonymous over synonymous substitutions, and in MHC, they are commonly assumed to be involved in antigen‐binding (Hughes & Nei, 1988). Although selection analyses have been performed previously for MHC class II genes of this species (Migalska et al., 2019; Scherman et al., 2014), we repeated them with all nucleotide MHC sequences identified in this study, combined with other bank vole DQB and DRB sequences available in Genbank (details in Supporting Information Methods). Similarly, in the supertyping procedure, all inferred amino acid MHC sequences from this study were used, in combination with other bank vole MHC sequences available in Genbank (sequence identifiers available in Supporting Information Data set), provided they covered PSSs. Altogether, we used 176 unique amino acid sequences for DQB, and 115 sequences for DRB. Supertyping was performed separately for DQB and DRB, according to the procedure described in Migalska et al. (2019). Details of supertyping are provided in Supporting Information Methods, and examples of amino acid sequences belonging to each supertype are shown in Figures S3 and S4.

Bank vole microsatellite data

Neutral, genome‐wide diversity was assessed with a set of microsatellite markers. We first preselected 18 previously described loci and validated their utility in our subpopulations by testing for departures from the Hardy‐Weinberg equilibrium, frequency of null alleles and linkage disequilibrium. Details of multiplexing, PCR conditions, genotyping and the validation procedure are in Supporting Information Methods, Tables S2 and S3. We then calculated individual multilocus heterozygosity (MLH) for 861 individuals based on 13 loci that remained after the validation procedure: Cg2C5, Cg3A8, Cg13F9, Cg13G2, Cg14E1, Cg15F7, Cg16A3, Cg16E4 (Rikalainen et al., 2008), Cg1F11, Cg4F9, Cg10A11, Cg12A7, Cg17C9 (Guivier et al., 2011, redesigned from Rikalainen et al., 2008). MLH in its “raw” form was calculated as a proportion of heterozygous loci. We also calculated MLH as standardized heterozygosity (SH) to correct for missing genotypes (Coltman et al., 1999) and as homozygosity‐by‐loci (HL) (Aparicio et al., 2006), which weights loci by their expected heterozygosity. Finally, we calculated pairwise F ST (population differentiation measure) between each subpopulation in a sampling year. Calculations were performed in Genepop v.4.6.9 (Raymond & Rousset, 1995; Rousset, 2008; dememorization = 10,000, batches = 100, iterations per batch = 5000).

Statistical analysis

In the first set of statistical models described in detail below, we examined correlations between individual‐level MHC diversity and helminth species richness. In the second set of models, we checked for specific associations between particular MHC class II supertypes and each individual helminth species in our study. In the third set of models, we investigated whether susceptibility to infection of MHC class II supertypes was related to their frequency in previous sampling seasons. In all these analyses, we focused on MHC supertypes, rather than allelic variants, because the numbers of alleles were too large (compared to sample size), to allow meaningful analyses using presence/absence of particular alleles as predictors. The use of supertypes preserves information on functional diversity, and is therefore a preferred way to reduce data set complexity in a biologically relevant fashion (instead of, for example, using just the most frequent alleles). Moreover, while in the first set of models we used all identified DQB and DRB supertypes, in the two latter sets of models we removed from the analysis the rarest supertypes and highly correlated supertypes (details of these exclusion are in Supporting Information Materials, Figure S5) ‐ leaving 14 DQB and seven DRB supertypes (i.e., DQ 01, DQ 04, DQ 05, DQ 07, DQ 08, DQ 09, DQ 10, DQ 11, DQ 14, DQ 15, DQ 16, DQ 17, DQ 19, DQ 22, DR 01, DR 08, DR 10, DR 12, DR 13, DR 19, DR 20). All statistical calculation were performed in R, v. 4.0.2 (R Core Team, 2020).

Individual MHC diversity and helminth species richness

We used Spearman's rank correlations to test whether standardized, microsatellite‐based heterozygosity (SH) correlated with either the number of MHC nucleotide variants or number of MHC supertypes. The overall association between the individual number of DQB and DRB supertypes and helminth species richness was estimated with generalized linear model (glm, stats package part of R, v. 4.0.2), using quasipoisson error distribution and a log link function. The model also included additional, “basic” explanatory variables: animal's sex, age, sampling site and year, as these factors had been shown previously to play a role in helminth infections of bank voles (Grzybek, Bajer, Bednarska, et al., 2015). We first tested if the model with MHC supertypes explained data better than the one with the basic variables alone. We did so by comparing via anova (stats, v. 4.0.2) a null model (with basic variables only) to the ones that contained either the sum of MHC class II supertypes (DQB +DRB), or that contained numbers of DQB and DRB supertypes as separate predictors. We included both linear and quadratic terms for MHC supertypes, to account for a possible “optimality” effect. Furthermore, we tested whether addition of individual, standardized heterozygosity (SH) as an explanatory variable improved the fit of the model.

Associations between particular MHC class II supertypes and prevalence/abundance of helminths

We explored potential associations between MHC class II supertypes and the prevalence or abundance of helminths in particular individuals using a multivariate method which simultaneously examines several helminth species as response variables – redundancy analysis (RDA). The advantage of RDA over the multiple glm approach is that it avoids problems associated with multiple testing, haunting analyses with many MHC variants as predictors. Here, we performed a partial RDA with MHC genotypes, sampling year and heterozygosity (SH) as explanatory variables, age and sex as conditioning variables (covariates), and helminth prevalence or abundance matrices as response variables. Partial RDA treats explanatory variables that are probably influential, but less pertinent to the aim of the study, as covariates. Their contribution to the variability in response variables is factored out, and then the ordination is performed on the residual variability. A previous study on the studied bank vole subpopulations revealed that the sampling site was the main determinant of the composition of helminth communities at the three sites (Grzybek, Bajer, Bednarska, et al., 2015). Consequently, ordinations of parasitological data set grouped data points primarily by sites (PCAs of parasite abundance data in Figure S6), and this aspect could dominate the ordination‐based analytical method. Thus, to focus our analysis on the potential effects of the MHC genotypes instead, we performed a separate RDA for each site. The analyses were performed with the rda function from the vegan package (v. 2.5–6) (Oksanen et al., 2019). We first constructed “full” models, with the helminth matrices comprising either prevalence or abundance of the helminth species present at each site as the response variables, and year of sampling (as factor), SH, DQB and DRB supertype presence‐absence data, and interactions between year and each supertype as possible explanatory variables. For analysis of helminth abundance, we first log‐transformed the data (as proposed by Anderson et al., 2006), and then standardized it by subtraction of a mean and division by the standard deviation. Presence of an interaction between MHC supertype and year would indicate a change in the effect that a given MHC supertype has over the years, suggesting for example, parasite adaptation to the particular supertype. Apart from the amount of constrained variance ‐ which is an equivalent to the R 2 of a multiple regression, ‐ we also report an adjusted R 2 (adjusted based on the number of predictors). Next, we assessed the significance of the full model using an ANOVA‐like permutation tests (anova.cca) for the joint effect of constraining (explanatory) variables. If the full model was significant overall or only marginally non‐significant (p < .1), we followed with a more refined selection of predictors, to identify the most important explanatory variables and to build “optimal” models. We used the forward model selection procedure (ordiR2step), which builds models forward (with scope = full model), by maximizing adjusted R2 at every step, until the adjusted R2 starts to decrease. The adjusted R 2 is based on simulations, therefore different runs of ordiR2step can give different results. The final reported values were averaged over 50 runs, and optimal models were built of predictors consistently kept by ordiR2step by the majority (>50%) of these runs.

MHC class II supertype frequency and infection status

Under NFDS parasites adapt to the most frequent MHC variants. This predicts that the more frequent a supertype was in a recent past, the higher the probability of hosts carrying that supertype being infected (i.e., higher prevalence of infection relative to hosts carrying other supertypes) and/or the higher the abundance of infection in the following years. An overview of the supertype frequencies in each site and their temporal fluctuations are shown in Figure S7. To investigate relationships between supertype frequencies in a preceding sampling year (i.e., year t‐1) and its infection status in a following sampling year (i.e., year t), we first calculated for each helminth species a supertype‐wise infection score (Biedrzycka et al., 2018). Infection score was the difference between the mean residual infection status of carriers and noncarriers of a given supertype in a given year. Residual infection status was the residual from a model accounting for sex and age only: that is, a binomial glm of prevalence or a quassi‐poisson glm of abundance of helminths, with sex and age as predictors. We did not include heterozygosity in these models, because it did not affect the parasite burdens in any previous test (see Results). Positive values of the infection score indicate that carriers were more often or more heavily infected in that year than noncarriers of a given supertype. Next, we checked if the infection score in year t was affected by the frequency of a given supertype in year t‐1 (e.g., infection score in year 2010 was affected by frequencies recorded at year 2006, etc.). The four years that elapsed, on average, between sampling years, probably constitute several generations of pathogens thus allowing time for adaptation. Although the exact longevity of the helminths in this study has not been reported, similar rodent helminths have generation times from two weeks (e.g., Heligmosomoides spp.) to several weeks (e.g., Mastophorus muris), and can live for several months (e.g., Heligmosomoides bakeri) (Bryant, 1973; Robinson et al., 1989). In the final step, for each helminth species, we constructed a set of linear models with infection score in year t as a response variable, and supertype frequency in year t‐1, sampling site, year of sampling and pairwise interactions between these terms. We further compared AICc values of all possible models (calculated with function dredge from MuMIn, v. 1.43.17, Bartoń, 2020), and chose the model with lowest AICc, provided that the ΔAICc between it and the next best model was ≥2. In situations in which several top models had similar AICc values (ΔAICc <2) we parsimoniously chose the simplest model among them. Since each helminth species was tested four times (infection score calculated based on prevalence and abundance and models for DQB and DRB supertype frequencies), we applied Bonferroni‐corrected p‐value of .0125.

RESULTS

MHC sequencing, genotyping and supertyping

During genotyping, the average number of processed reads per amplicon was 2557 (SD: 1069) for DQB, and 2663 (SD: 1065) for DRB. Out of 1022 individuals for which tail clips were available, we successfully determined MHC DQB genotypes for 1006 individuals, and DRB genotypes for 983 individuals. Information on animal sex, age, site (subpopulation) of origin and helminth abundance/prevalence was available for 972 individuals. After combining these data with MHC genotypes, we obtained complete information including DQB genotype for 910 individuals, DRB genotype for 902 individuals, and full MHC class II genotypes (both DQB and DRB) for 892 individuals (Supporting Information Data set). Among the fully genotyped individuals, there were 89 distinct nucleotide (nt) DQB variants (termed in the text as “alleles”, although we cannot assign them to loci, so they represent a collection of alleles derived from many loci). All these alleles were in an open reading frame and had no stop codons. They translated into 85 unique amino acid (aa) sequences, that belonged to 20 supertypes. The phylogenetic tree of alleles from this study is shown in Figure S8. There were also 87 distinct DRB nt alleles, but translation revealed that six had stop codons, and so were deemed potential pseudogenes and these were removed from further analysis. Also, at this stage, all the sequences that belonged to Cluster I (see Materials and methods and Supporting Information Methods) were removed. This left 64 distinct DRB nt alleles (“Group III” only) that translated into 45 distinct aa sequences, belonging to 11 supertypes. The phylogenetic tree of all functional alleles identified in this study, together with previously identified, expressed variants (Migalska et al., 2019) is shown in Figure S9; with designation of alleles belonging to Cluster I and II. The mean number of both nt and aa DQB variants per individual was 4.15 (SD: 1.21, range: 1–8), and this corresponded to 3.75 (1.10, 1–7) supertypes. The mean number of DRB variants was 4.34 (1.15, 1–8) for nucleotide sequences, 4.18 (1.08, 1–8) for amino acid sequences, and 3.01 (0.70, 1–5) for supertypes.

Not all alleles were present at all sites, and subpopulations differed in the numbers of “private” alleles (i.e., alleles not found elsewhere) – from only three DQB and four DRB such nt alleles found in Urwitałt, to 20 and 15 such alleles found in Pilchy (Table S4). Only between 34%–36% of nt alleles occurred in all three sites, and this number was only slightly raised when we looked at diversity at the amino acid level (~38%). In contrast, the majority of supertypes was shared between sites: 17 out of 20 DQB supertypes, and eight out of 11 DRB supertypes. Detailed statistics on allele sharing are in Table S4.

Microsatellite variation

Microsatellite‐based heterozygosity (calculated for individuals for which at least 9 out of 13 loci were successfully genotyped) was calculated for 853 individuals, for which complete information on hosts, their helminth infection status and MHC class II genotypes of both DQB and DRB were also available (Supporting Information Data set). Information about the sample size per site, year, sex and age class are provided in Table S5. All three measures of heterozygosity (MLH, SH, HL) were highly correlated (>0.98, p < .001), and so SH was used in all further analyses. Mean (across sampling years) pairwise F ST values between sites were moderate (0.060, range: 0.056–0.066), and higher than F ST values between consecutive sampling years within each site (mean for Pilchy: 0.011, Urwitałt: 0.013, Tałty: 0.005) (Table S6).

Individual heterozygosity, MHC diversity and helminth species richness

Individual standardized heterozygosity (SH) did not significantly correlate with either overall number of MHC class II nucleotide variants (rho < 0.001, p = .98) or sum of MHC class II supertypes (rho = 0.017, p = .63). There was also no correlation between heterozygosity and the numbers of DRB or DQB nucleotide variants/supertypes when these were considered separately (DRB: –0.009/0.026; DQB: 0.010/–0.005, all p > .4). SH did not significantly affect helminth species richness, and adding it as an explanatory variable to a basic model containing age, sex, site of origin (subpopulation) and year of sampling was not supported by comparison of nested models in ANOVA (F 1,843 = 0.308, p = .579). In contrast, the number of DQB and DRB supertypes did affect the number of helminth species found in individuals, but the relationship was complex. Adding numbers of DQB and DRB variants as separate predictors and incorporation of quadratic terms significantly improved the fit of the model ‐comparison to a model with basic variables (see Materials and Methods) with ANOVA: F 4,840 = 3.469, p = .008. Neither adding sums of DQB and DRB as one predictor, nor using only linear terms, significantly improved the fit compared to a model with basic variables only (F 2,842 = 1.980, p = .139 and F 2,842 = 0.975, p = .379, respectively). Quadratic terms were significant in both cases, but the numbers of DRB and DQB supertypes showed opposing effects on species richness. Individuals with an intermediate number of DRB supertypes had the lowest number of helminth species, while the opposite was true at intermediate numbers of DQB supertypes (Figure 1). Moreover, species richness significantly increased with age, was lower overall in 2010, and differed between sites (with animals from Tałty and Pilchy harbouring fewer species than those from Urwitałt). The estimated coefficients, t and p‐values are shown in Table S7.

FIGURE 1

Partial residual plot with a line representing the relationship between the number of (a) DQB and (b) DRB supertypes on individual helminth species richness estimated from a generalized linear model (see Materials and Methods). Results were visualized with package visreg (Breheny & Burchett, 2017); variables other than the focal (i.e., number of DQB or DRB supertypes) are held constant at the mean or at the most common factor

Associations between helminth communities and MHC class II supertypes

RDA analyses were used to test for a combined effect of MHC class II supertypes and covariates on prevalence and abundance of infection with multiple helminth species. The fraction of variance explained by the explanatory variables (i.e., constrained variance, the adjusted R 2 of optimal models) oscillated between 0.06–0.10 (Table 1), which was comparable to the portion of conditional variance that could be explained by model covariates (sex and age) – between 0.07 and 0.10. Heterozygosity (SH) was never included in the RDA’s optimal model, whereas sampling year was always included, and significantly affected helminth prevalence at each study site. In several cases RDA’s optimal models revealed that specific MHC class II supertypes also affected prevalence and/or abundance of helminths (Table 1), but no interaction between MHC supertype and year was retained in any of the optimal models. In Urwitałt, DRB supertypes 8 and 12 affected prevalence and/or abundance of helminths; in Tałty ‐ DQB supertype 19 affected prevalence; in Pilchy – DRB supertype 8 affected abundance. The nature of correlations between MHC class II supertypes and particular helminth species can be inferred from the angles between the vectors in triplots visualizing the results of the RDAs (Figure 2, optimal models): perpendicular angles mean no correlation, vectors pointing in the same direction – a positive correlation, and vectors pointing into opposing directions – a negative correlation. The most pronounced associations are visible for A. tianjinensis and A. annulosa. In Tałty, supertype DQB 19 is negatively associated with A. tianjinensis prevalence, whereas in Pilchy, a negative correlation is seen for the supertype DRB 8 and abundance of this species. Surprisingly, the same DRB 8 supertype is positively correlated with abundance of A. tianjinensis (and A. annulosa) in Urwitałt. At this site, supertype DRB 12 correlated negatively with both abundance and prevalence of A. tianjinensis and A. annulosa, having an opposing association to that of DRB 8.

TABLE 1

Summary of the redundancy analysis (RDA) results

Data set	Model	Parameter	Urwitałt (Site 1)	Tałty (Site 2)	Pilchy (Site 3)
Parasite abundance	Full model	Prop. constrained variance (R 2)	0.293	0.308	0.302
		Adjusted R2	0.070	0.028	0.080
		F value	1.312	1.099	1.356
		df	78 (213)	85 (184)	74 (207)
		Significance (Pr >F)	0.005**	0.200 n.s.	0.006**
	Optimal model	Prop. conditional variance	0.096	Not evaluated	0.076
		Prop. constrained variance (R 2)	0.107		0.104
		Adjusted R 2	0.094		0.093
		F value	7.672		8.822
		df	5 (286)		4 (277)
		Significance (Pr >F)	0.001***		0.001***
		Explanatory variables	year +DR_12 + DR_08		year +DR_08
Parasite presence‐absence	Full model	Prop. constrained variance (R 2)	0.276	0.312	0.314
		Adjusted R 2	0.036	0.044	0.101
		F value	1.148	1.161	1.468
		Df	78 (213)	85 (184)	74 (207)
		Significance (Pr >F)	0.050*	0.066.	0.001***
	Optimal model	Prop. conditional variance	0.066	0.102	0.089
		Prop. constrained variance (R 2)	0.078	0.0702	0.101
		Adjusted R 2	0.067	0.058	0.093
		F value	6.556	5.616	11.51
		df	4 (287)	4 (265)	3 (278)
		Significance (Pr >F)	0.001***	0.001***	0.001***
		Explanatory variables	year +DR_12	year +DQ_19	year
		Evaluated parasites	H. mixtum M. muris A. tianjinensis A. annulosa C. henttoneni	H. mixtum A. tianjinensis A. annulosa C. henttoneni	H. glareoli M. muris A. tianjinensis C. henttoneni

Full models, models with the parasitological matrices of either prevalence or abundance of helminth species present at each site as the response variables, and heterozygosity (SH), year of sampling, DQB and DRB supertype presence‐absence data, and interactions between year and each supertype as the explanatory variables; Optimal models, reduced models, built by adding explanatory variables to a null model, aiming at maximizing adjusted R 2 at every step, until the adjusted R 2 started to decrease; Prop. constrained variance (R 2), the amount of variance uniquely explained by the explanatory variables; Prop. conditional variance, the amount of variance explained by the covariates (sex and age). The remaining variance is unexplained by the model. The significance of the models was calculated with ANOVA‐like permutation tests for the joint effect of constraining (explanatory) variables. Df, degrees of freedom used (free remaining in brackets); Explanatory variables, explanatory variables included in the optimal model in forward model selection procedure. Significance codes: n.s. >0.1 >. >0.05 > * >0.01 > ** >0.001 ≥***

FIGURE 2

Triplots of the RDA results of optimal models that kept particular MHC class II supertypes among explanatory variables. Individual scores are plotted as circles (light red for individuals from Urwitałt, light yellow–Tałty, light blue–Pilchy, on‐line version of the article); in red (dashed lines) are denoted vectors of response variables (helminth abundance or prevalence); in blue (solid lines) are vectors representing the explanatory variables (MHC class II supertypes or sampling years). The angles between response and explanatory variables reflect their correlations (see main text for more details). RDA was plotted with symmetrical scaling (by square root of eigenvalues) of scores of individuals and their helminths. The amount of variance explained by the first two RDA axes is given in parenthesis

MHC class II supertype frequency and helminth infection status

In the majority of helminth‐MHC gene combinations, the supertype‐specific infection score in year t was best explained by an intercept only model (Table 2). Nonetheless, in five cases (out of 24), other terms were also included, but only three of these models had predictors with a significant (after Bonferroni correction) impact on the infection score (Table 2). These were: (1) prevalence‐based infection score of the locally new helminth, A. annulosa, and DQB supertype frequencies; (2) abundance‐based infection score of A. tianjinensis and DQB supertype frequencies; and (3) prevalence‐based infection score of H. glareoli and DRB supertype frequencies.

TABLE 2

Summary of the explanatory variables included in best parsimonious models explaining parasite infection scores

Infection score (year t)	MHC supertypes	Freqt−1	Site	Year	Freqt−1 ×site	Freqt−1 ×year	Year ×site	AICc	N. models with ΔAICc <2
A. annulosa PA	DQ	+	+	−	+	−	−	40.2	0
	DR	−	−	−	−	−	−	20.9	1
A. annulosa ABD	DQ	−	−	−	−	−	−	156.5	2
	DR	−	−	−	−	−	−	146.9	3
A. tianjinensis PA	DQ	−	−	−	−	−	−	130.0	1
	DR	+	−	−	−	−	−	72.8	1
A. tianjinensis ABD	DQ	+	+	−	+	−	−	444.6	0
	DR	−	−	−	−	−	−	310.9	1
C. henttoneni PA	DQ	−	−	−	−	−	−	72.4	3
	DR	−	−	−	−	−	−	40.2	1
C. henttoneni ABD	DQ	−	−	−	−	−	−	113.0	1
	DR	−	−	−	−	−	−	72.7	1
M. muris PA	DQ	−	−	−	−	−	−	84.0	2
	DR	−	−	−	−	−	−	52.1	2
M. muris ABD	DQ	−	−	−	−	−	−	163.3	3
	DR	−	−	−	−	−	−	57.5	3
H. mixtum PA	DQ	−	+	−	−	−	−	65.6	1
	DR	−	−	−	−	−	−	78.5	1
H. mixtum ABD	DQ	−	−	−	−	−	−	114.6	2
	DR	−	−	−	−	−	−	100.3	0
H. glareoli PA	DQ	−		−		−		55.8	1
	DR	+		+		+		31.9	0
H. glareoli ABD	DQ	−		−		−		68.7	0
	DR	−		−		−		43.9	0

The first column indicates which parasite's infection score at year t was the response variable and whether it was calculated based on presence‐absence data (PA) or abundance of infection (ABD) data; the second column indicates whether MHC supertype frequencies were calculated for DQB or DRB genes. Columns three to six represent possible explanatory variables; +/− indicates whether the stated term was included in the best model or not. H. glareoli is present only in one site (Pilchy), therefore neither site, nor its interaction terms could have been included in its models (thus empty fields in the table). The second‐to‐last column shows the AICc value of the chosen model. The last column shows how many models differed by less than 2 from the model with lowest AICc – if this values is >0, the simplest model among the “top” options was selected as the most conservative and parsimonious choice. Models with p‐values significant after Bonferroni correction are in bold.

In the first two cases the supertype frequency in a preceding sampling year (t‐1) significantly affected infection score in the following sampling year (t), but this effect was site‐specific (significant interaction term). In the first model, in Urwitałt, but not in Tałty, DQB supertype frequency in year t‐1 was positively associated with A. annulosa infection score in year t (Figure 3a, t = 2.811, p = .006, the whole model in Table S8). In the second model, DQB supertype frequency in year t‐1 was positively associated with A. tianjinensis infection score in year t in Tałty, but not in the two remaining sites (Figure 3b, t = 3.948, p < .001, Table S9).

FIGURE 3

The estimated effect of DQB supertype frequency in a preceding sampling year (t‐1) on (a) prevalence‐based infection score of A. annulosa in year t, (b) abundance‐based infection score of A. tianjinensis in year t. Positive values of the infection score indicate that carriers of a given supertype are more often or more heavily infected by the given parasite than noncarriers. Negative values indicate that carriers are less infected, compared to noncarriers, suggesting a protective effect. Results were visualized with package visreg; variables other than the focal (i.e., DQB supertype frequency) are held constant at the mean or at the most common factor, data points represent partial residuals and their number corresponds to the number of DQB supertypes analysed at each site times the number of years it was recorded. In (a) only two sites are shown, because A. annulosa is virtually absent from the third site (Pilchy), and data from this site were therefore not included in the model

The last case where explanatory variables were included in the best, parsimonious model concerns the prevalence‐based infection score of H. glareoli, a parasitic species present predominantly in one site, Pilchy (Grzybek, Bajer, Bednarska, et al., 2015). Here, the effect of DRB supertype frequency in year t‐1 on parasite infection score in the following sampling year changed over time (significant interaction terms for frequencyt‐1 × year; 2006: t = 3.370, p = .005, 2010: t = 3.792, p = .002, Table S10). Initially, more frequent supertypes in year 1999 were associated with lowest infection scores in 2002, but this changed sharply in the following years (Figure 4). In the two subsequent sampling seasons the relationship between frequency and infection score disappeared, and the infection scores became similar for both previously frequent and infrequent supertypes (Figure 4).

FIGURE 4

The estimated effect of DRB supertype frequency in a preceding sampling year (t‐1) on prevalence‐based infection score of H. glareoli in year t, shown for the subsequent sampling years. The results were visualized with package visreg, variables other than the focal (i.e., DRB supertype frequency) are held constant at the mean or at the most common factor, data points represent partial residuals and their number corresponds to the number of DRB supertypes analysed in each year. H.glareoli is only present in Pilchy, and so other sites were not taken into account

DISCUSSION

In the present study we used a data set spanning 11 years, reporting intestinal helminths infecting bank voles in three sites in the Mazury Lake District in northeastern Poland, to study spatiotemporal associations between complex helminth communities and MHC class II supertype diversity. Our analyses controlled statistically for age, sex and location, since each factor is known to influence helminth burden significantly, as earlier reported for this data set, and extensively discussed by Grzybek, Bajer, Bednarska, et al. (2015). We found that individual MHC diversity was associated with helminth species richness, but the patterns were inconsistent between genes. Furthermore, we showed that prevalence and abundance of helminths are significantly associated with MHC supertypes, adding to the existing evidence for such associations in natural populations (reviewed in Radwan et al., 2020; Spurgin & Richardson, 2010). More importantly, however, we found that such associations vary in space and time.

With respect to spatial variation, we found that in different subpopulations different supertypes (e.g., DRB 08 and DQB 19 in Pilchy and Tałty, respectively) were associated with prevalence/abundance of the same helminth species, A. tianjinensis. One possible explanation is that subpopulations of voles harbour different strains of A. tianjinensis, which might differ in the composition of molecules which are presented by MHC. Further support for this explanation is provided by the same supertype DRB 08 having a contrasting effect on A. tianjinensis infection in the Urwitałt and Pilchy subpopulations (Figure 2). Such a scenario is consistent with Red Queen dynamics whereby host‐parasite coevolution proceeds independently in different populations (rendering them “out of phase” in terms of MHC‐infection associations), but also with FS, if A. tianjinensis strains are stably associated with different sites. At present, lack of data on the genotypes of the helminths in this study does not allow us to differentiate between these hypotheses. However, some supports for FS stems from a weak negative association between H. glareoli, a helminth species found predominantly in the Pilchy subpopulation, and DRB 08 supertype. If this supertype is favoured by voles in Pilchy, but disfavoured in other populations, then such spatially fluctuating selection could contribute to its polymorphism, and moderate gene flow between populations (mean F ST = 0.06, a value close to that reported by Kloch et al., 2010) would be conducive to its maintenance (Lenormand, 2002). Indeed, DRB 08 was associated with increased infection risk with A. annulosa and A. tianjinensis in Urwitałt. At the same time, in Pilchy, it was clearly associated with lower risk of infection with A. tianjinensis (as mentioned above), and to a lesser degree, with C. henttoneni, which complicates interpretation of the role of local parasites in causing fluctuating selection on MHC. Overall, the spatial variation in MHC supertype – helminth associations that we detected is consistent with both NFDS and FS, but as highlighted by Spurgin and Richardson (2010), discriminating between these hypotheses is difficult with spatial data alone.

As for temporal patterns, FS is thought to result from temporal variation in parasite composition in time. In our data set, prevalence/abundance of some helminths did fluctuate, but, except for a newly‐emerged A. annulosa, roughly the same set of helminth species was observed at a given site across years (Grzybek, Bajer, Bednarska, et al., 2015; Figure S2). It thus appears that selection from different helminth species fluctuated at best moderately in time, at least on the timescale we investigated. As for NFDS, significant supertype × year associations could be expected when alleles lose resistance as a result of parasite adaptation. In our RDA models, such associations were not significant. However, the spread of escape mutations in parasites requires time, and consequently the signal of temporal resistance of particular MHC variants may be difficult to detect.

Our analyses of temporal patterns allowed a more direct test of NFDS: we predicted that as a result of adaptation to the most common MHC types, we would see a positive association between prevalence/abundance of helminths related to a supertype and its frequency in a previous sampling period. We found three significant associations, with the expected, positive association, seen in two of these three cases. Noteworthy, we did detect an association in the predicted direction for a helminth novel to all the Mazury subpopulations, A. annulosa, whereas for the two other helminths the association was either also positive, as predicted (A. tianjinensis in Tałty), or varied between sampling periods (H. glareoli in Pilchy). These results highlight that the patterns predicted by NFDS can be complicated by past history of selection and uncertainty as to the rate of parasite adaptation. As we argued in the introduction, novel parasites are not subject to the former constraint.

Since A. annulosa was first recorded among the monitored subpopulations in 2002, it remained virtually absent from Pilchy, but steadily rose in both prevalence and abundance in Urwitałt, and to a much lesser degree in Tałty (Grzybek, Bajer, Bednarska, et al., 2015; Figure S2). This more pronounced spread may be the reason why a significant association between DQB supertype frequency and its prevalence‐based infection score was observed only in Urwitałt (Figure 3a). Positive values of the infection score indicate that carriers of a given supertype were more likely infected than noncarriers, and so an increase in these values implies pathogen adaptation. Thus, the more frequent the occurrence of a supertype in a preceding sampling year, the more likely bank voles carrying that supertype were to be infected with A. annulosa in the following sampling year, a pattern predicted by NFDS.

In contrast to A. annulosa, A. tianjinensis was present in all the Mazury subpopulations from the very beginning of the sampling, although its prevalence and abundance differed markedly between the three sites (Grzybek, Bajer, Bednarska, et al., 2015; Figure S2). In 1999–2006 it infected about 60% of voles in Pilchy, whereas in the two remaining sites, this helminth was initially rare in 1999 (infected ≤10% individuals), but quickly rose in prevalence, peaking in 2006 (over 60% in Tałty, over 40% in Urwitałt). Changes in its abundance were less pronounced, except for a visible peak in Tałty in 2006 (three‐fold increase). It seems plausible that the fast spread of A. tianjinensis in Tałty was associated with adaptation to the most common MHC supertypes. Perhaps that is why this was the only site for which we detected a significant, positive association between recent DQB supertype frequencies and their abundance‐based infection scores (Figure 3b).

The temporally‐dynamic pattern of host‐parasite coevolution that makes capturing NFDS difficult in snapshot studies involving parasites that are well established in populations is illustrated by the association between DRB supertype frequencies and infection scores of H. glareoli, the helminth present almost exclusively in Pilchy. The significant year ×frequency association we detected appears to have resulted from the DRB supertypes that were more frequent in 1999 having a protective effect, with the pattern disappearing in the following years (Figure 4). The association found in 2002, which is opposite to the predictions of NFDS, might possibly reflect pre‐1999 coevolutionary history, when resistant MHC genotypes were increasing in frequency, possibly peaking around 1999–2002, and H. glareoli at Pilchy had not yet adapted sufficiently to evade them. The change in slope in following seasons might reflect progressing adaptation of the helminth to common MHC types, with the expectation that the relationship will become positive in the following season. Ultimate verification of this speculative scenario would require tracing the evasion genes in H. glareoli which are not yet known.

Both DRB and DQB occur in several copies in bank vole genomes, and the ranges of variants per individual we found are consistent with those reported in other studies (Galan et al., 2010; Migalska et al., 2019; Scherman et al., 2014) including a study that did not use high‐throughput sequencing but a mix of SSCP and Sanger sequencing (Axtner & Sommer, 2007). Coamplification of multiple loci precluded direct testing of HA, another important mechanism leading to the maintenance of MHC polymorphism (Doherty & Zinkernagel, 1975; Radwan et al., 2020). However, gene duplication and divergence, a common mechanism driving MHC evolution can achieve fixation of these favourable properties of heterozygotes, and thus should be favoured by selection in the long term (Bentkowski & Radwan, 2019). Yet, the limit to MHC expansion may result from trade‐offs with autoimmunity (Migalska et al., 2019; Nowak et al., 1992; Woelfing et al., 2009), leading to stabilizing selection on optimal MHC gene number. Previous research has provided some support for this optimality hypothesis: bank voles sampled in 2005 with an intermediate number of DRB alleles carried the lowest number of helminth species (Kloch et al., 2010). Here, we found a similar pattern when looking at the number of DRB supertypes alone, but we found an opposite effect for the numbers of DQB supertypes (Figure 2). No relevant functional differences between DR and DQ loci have been reported previously, making the observed relationship even more puzzling. Nonetheless, the contrasting shapes of association for the two gene families might explain why no association was seen when the sum of DQB and DRB supertypes was analysed, especially given very weak intra‐individual correlation between the numbers of DQB and DRB supertypes (0.16, p <.001). These results highlight that associations between MHC functional diversity and susceptibility to infection are not a rule, and particular MHC alleles/supertypes may have a higher impact on parasite resistance than their diversity (Kloch et al., 2013; Phillips et al., 2018). Information on the entire MHC class II haplotypes (not possible with current methodology) could also shed more light on the possible, underlying associations.

In conclusion, our findings on newly colonizing parasites, suggest adaptation of bank vole helminths to MHC consistent with the Red Queen coevolutionary dynamic, thus complementing still limited evidence for its action in nature. The positive association between supertype frequency and infections with the novel colonizing parasite A. annulosa (that we observed in a subpopulation in which this helminth had spread the most), provides strong support for NFDS. At the same time, our results for H. glareoli, which was already prevalent at the outset of our time series, show that the association between MHC supertype frequency and infection score can vary in time. This highlights difficulties in testing NFDS for established parasites, for which the phase of host parasite coevolution is unknown. The spatial patterns of MHC‐helminth associations that we found are consistent with the role of FS, although they could also result from NFDS. We note, however, how the patterns of MHC‐infection associations we found translate to selection pressures on MHC will depend on the effect of infection on fitness, which in our system remains unknown. Finally, our results indicate that the shape of selection for MHC diversity can vary not only among species (Radwan et al., 2020), but also between gene families within the same set of populations. This casts doubt on the dominating role of individual MHC diversity (or heterozygosity in the case of a single locus) as a universal mechanisms favouring MHC polymorphism.

AUTHOR CONTRIBUTIONS

Jacek Radwan and Magdalena Migalska designed the research; Magdalena Migalska and Karolina Przesmycka performed the laboratory procedures carried herein; Mohammed Alsarraf, Anna Bajer, Jerzy M. Behnke, Jolanta Behnke‐Borowczyk and Maciel Grzybek provided parasitological data and bank vole tissue samples; Magdalena Migalska, Karolina Przesmycka and Jacek Radwan analysed the data and wrote the article. All authors read, corrected and agreed the final version of the manuscript.

CONFLICT OF INTEREST

The authors declare no conflicts of interest.

BENEFIT‐SHARING STATEMENT

Benefits from this research accrue from the sharing of our data and results on public databases as described above.

Supplementary Material

Click here for additional data file.