Introduction

An expected consequence of global change is the modification of pathogen communities and host–pathogen interactions (Kutz et al. 2005; Altizer et al. 2013; Dobson et al. 2015), especially in ecosystems where climate has been limiting the expansion of pathogens (Callaghan et al. 2004), e.g., in Arctic environments (Kutz et al. 2005, 2009). Anthropogenic activities could also lead to changes in the pathogen regime through transport and introduction of new pathogens, or a switch of pathogens from domestic to wild host species (Peeler and Feist 2011). It is crucial to understand how such alterations will impact individual fitness in wild species to better predict potential changes in population dynamics, and to adapt conservation and management policies.

For hosts, immune functions are costly and often generate trade-offs between energy allocation to immunity and other functions, e.g., maintenance, growth, or reproduction (Soler et al. 2003; Moreno-Rueda 2011; Souchay et al. 2013; Griesser et al. 2017). In many cases, pathogen-infected individuals are in poorer condition (Moreno-Rueda 2011), have lower reproductive success (Griesser et al. 2017), and lower survival (Souchay et al. 2013) in comparison with their healthy counterparts. Because infections can impact fitness, immune responsiveness is probably under strong stabilizing selection that tends to minimize both the costs associated with immune response and pathogen-induced damages (Lochmiller and Deerenberg 2000; Moret 2003). Conversely, genes responsible for pathogen recognition usually are under balancing selection as organisms need to cope with various and often-changing immune challenges (Bernatchez and Landry 2003; Sommer 2005).

In vertebrates, genes of the major histocompatibility complex (MHC) play a key role in pathogen recognition (Janeway et al. 2001). They code for cell-surface glycoproteins that bind to and present antigens to T cells, thus allowing a specific defense against recognized pathogens (Klein 1986). MHC genes have high nucleotide diversity, meaning that several sites are polymorphic, especially those from the antigen-binding site (ABS). These regions also present a large number of alleles, are frequently duplicated, and may present individual copy number variations (Hughes and Yeager 1998; Sommer 2005; Piertney and Oliver 2006). This diversity allows for the recognition of various pathogens, and it is assumed that individual genetic diversity at MHC is positively correlated to immunocompetence and, as a result, to individual performance (Hughes and Yeager 1998; Piertney and Oliver 2006). Based on this assumption, several non-mutually exclusive mechanisms have been suggested to be involved in the maintenance of MHC polymorphism. These mechanisms are related to pathogen-driven balancing selection (i.e., heterozygote advantage, divergent allele advantage, negative frequency-dependent selection, and fluctuating selection) or to sexual selection (e.g., MHC-disassortative mate choice; Bernatchez and Landry 2003; Sommer 2005; Milinski 2006). Under the heterozygote-advantage hypothesis, heterozygote individuals at MHC genes recognize and can elicit a specific immune response against a larger array of pathogens in comparison with their homozygote counterparts (Doherty and Zinkernagel 1975). It has also been suggested that alleles that present high divergent physicochemical properties at their ABS have a lower overlap in the pathogen arrays they recognize (Lenz et al. 2013; Pierini and Lenz 2018). As a result, the advantage conferred by heterozygosity would be even greater in individuals carrying divergent alleles, and selection would act to maintain these highly divergent alleles (the divergent allele-advantage hypothesis; Wakeland et al. 1990; Pierini and Lenz 2018). Alternatively, under the negative frequency-dependent selection (or rare allele-advantage) hypothesis (Clarke and Kirby 1966), pathogens are selected to escape the defense of the most common MHC alleles in their hosts, therefore favoring the less common alleles. Due to their low frequencies in populations, rare alleles are generally found in heterozygote genotypes, giving a selective advantage to heterozygotes. In this scenario, the selective advantage of rare alleles diminishes as they increase in frequency (Clarke and Kirby 1966; Takahata and Nei 1990). The spatial and temporal heterogeneity of the pathogen regime can also prevent the fixation of alleles, as alleles are differentially selected in space or time (fluctuating selection; Hill 1991). Finally, it has been shown that MHC-disassortative mate choice (based on odorant cues associated with MHC genes; Penn 2002) could also act to maintain MHC diversity, as individuals choose to mate with partners that will optimize or maximize their progeny’s MHC diversity (Reusch et al. 2001; Baratti et al. 2012; Huchard et al. 2013). It is hypothesized that MHC-disassortative mate choice is either a response to selective pressures exerted by pathogens (Reusch et al. 2001) or a way to avoid inbreeding (Huchard et al. 2013).

Most of the mechanisms described above would result in an advantage for individuals with high MHC diversity and consequently, several studies have investigated the link between MHC diversity and individual performance (or heterozygosity–fitness correlations (HFC); Bernatchez and Landry 2003; Sommer 2005). One challenge in doing so is that MHC diversity could be linked to another locus genotype (e.g., Aguilar et al. 2004) or to genome-wide diversity (e.g., Rico et al. 2015). As a result, associations between genetic diversity at MHC and individual performance could result from (1) a functional advantage of certain MHC genotypes (the “direct effect hypothesis”; Mitton 1997; David 1998), (2) an effect of another functional locus affecting fitness that is in linkage disequilibrium with MHC genes (the “local effect hypothesis”; Hansson and Westerberg 2008), or (3) inbreeding or outbreeding depression affecting fitness and genome-wide diversity (the “general effect hypothesis”; David et al. 1995; Hansson and Westerberg 2002). Discriminating among these three hypotheses can be challenging, but one might assume that local effects are less likely to occur for functional loci like MHC genes, which are themselves directly submitted to selection (Bernatchez and Landry 2003; Sommer 2005). In addition, some studies have ruled out the general effects by comparing MHC and genome-wide diversity and their respective relationship with fitness (e.g., Osborne et al. 2015). Due to their diversity and role in immunity, MHC genes can be used to study the direct effects of genetic diversity on individual fitness. This is especially of interest in the case of declining populations, as pathogens can have a significant impact on population dynamics (Albon et al. 2002).

Across the world, most caribou and reindeer (Rangifer tarandus) populations are declining (Vors and Boyce 2009). Climate change and anthropogenic development were raised as potential causes of decline (Vors and Boyce 2009), but only a few studies have investigated the role of pathogens and parasites on individual performance of caribou and reindeer individuals (Stien et al. 2002; Côté et al. 2005; Ballesteros et al. 2012; Pachkowski et al. 2013) and populations (Albon et al. 2002). Understanding how caribou-immune genes are linked to individual fitness could disentangle how changes in pathogen communities due to climate change (Kutz et al. 2005) might affect caribou population dynamics.

The objectives of our study were twofold: (1) to study correlations between MHC diversity and performance in migratory caribou, and (2) to test a new approach to disentangle general and direct effects of MHC diversity. We assessed the effect of genetic diversity of the MHC class II DRB exon 2 locus (hereafter MHC–DRB) on individual fitness components (body mass and survival) in two declining herds of migratory caribou (R. tarandus). This locus corresponds to the antigen-binding region of the most polymorphic gene of the MHC in caribou (Kennedy et al. 2010; Taylor et al. 2012), and it is the most polymorphic immune gene among ungulates (e.g., Quéméré et al. 2015; Portanier et al. 2019), making it potentially related to the performance of individuals. We hypothesized that MHC diversity would increase immunocompetence and positively affect individual performance. We proposed a comprehensive approach to discriminate between direct and general effects at MHC genes, comparing the relative influence of five MHC diversity indices on individual performance. These indices cover a gradient from allelic genetic diversity, estimated from DNA sequences, to functional diversity, inferred from the translated peptides. We predicted that these five indices would correlate differently with performance if general or direct effects were involved. Indeed, as inbreeding produces a correlation of heterozygosity across all loci (Balloux et al. 2004), we expected allelic diversity indices to correlate with genome-wide diversity if inbreeding was occurring (general effects). As a result, a correlation between allelic diversity indices and performance would give support to the general effect hypothesis if other indices did not correlate with performance. We did not expect to observe such pattern, as we had previously shown that the caribou populations studied here do not suffer from inbreeding depression (Gagnon et al. 2019). Instead, we predicted that direct effects were more likely to be involved in a potential correlation between MHC diversity and performance. As such, we expected that the functional indices of diversity, evidencing the level of functional divergence between alleles, would more strongly correlate with performance than the allelic diversity indices. A schematization of the research framework is presented in Fig. 1.

Fig. 1: Schematization of the research framework.
figure 1

Different indices of MHC genetic diversity were obtained from individuals of two herds of migratory caribou in northern Quebec and Labrador (Canada), and were independently correlated to fitness components (survival and body mass).

Full size image

Materials and methods

Study area and data collection

We sampled individuals from the two herds of migratory caribou in Northern Quebec and Labrador (Canada): Rivière-George herd (RG) and Rivière-aux-Feuilles herd (RAF). These two populations use large annual ranges (Fig. 2b) that are subdivided into smaller seasonal ranges used successively during the year. Individuals migrate annually between their calving ground, located in the northern part of their annual range, and their winter range, located in the southern part (Le Corre et al. 2017). In doing so, they encounter various environmental conditions, including different vegetation assemblages and varying intensities of human land use (Plante et al. 2018). Although the entire herd may gather at the same place at the same time, there are various wintering areas used by different animals (Le Corre et al. 2020). There is relatively low fidelity to the wintering area, and animals frequently switch among areas (Le Corre et al. 2020) and between herds (Boulet et al. 2007). RG and RAF herds are only slightly genetically differentiated (FST: 0.002 based on 6384 SNPs), and present no intrapopulation substructure (Gagnon et al. 2019). This prevents bias linked to the presence of a substructure, either spatial or genetic, in the herds we sampled. Both herds experienced large fluctuations in population size with strong declines in recent years (about 99% and 70% declines for RG and RAF herd, respectively; see Fig. 2a for details on population dynamics). Between 1996 and 2016, 152 caribou from RG (59 males and 93 females) and 252 from RAF herds (92 males and 160 females) were captured with a net gun fired from a helicopter following the guidelines from the Canadian Council on Animal Care and sampled for DNA (biopsies, hair, and/or blood samples). All samples were frozen, except for some biopsies that were conserved in ethanol (>70%). Blood samples were conserved with EDTA in vials. At capture, we estimated caribou’s age from the wear pattern of its incisor teeth (Hamlin et al. 2000). We assessed annual survival from capture to 2017 for 254 adults (≥2 years; RG n = 69, RAF n = 185) who were fitted with tracking collars using different satellite devices (Argos, Globalstar, or Iridium). Mortalities were signaled by collars after 12–24 h of inactivity and confirmed in the field when possible based on visual cues, when retrieving the collars. Among all marked caribou, 253 animals (1–12.75 years, RG n = 82 and RAF n = 171) were also weighed at capture to the nearest 0.1 kg using a hanging scale. Annual survival (from capture to 2017) and body mass (collected in January–March, June, October, or November) were used as indices of individual performance and included in models as response variables. Some of the sampled individuals (n = 114) were not included in the mass and survival models because they were in one of the following situations: (a) they were captured as calves, and their performance was too dependent on maternal effects (Skogland 1990; Veiberg et al. 2017) to be included in our models (n = 109), (b) they were not weighed at capture and were eventually hunted, so they could not be included in mass or survival models (n = 4), and (c) their age was unknown (n = 1). These individuals were nevertheless genotyped and included in the analysis described below under “Assessing population structure and selection.”

Fig. 2: Population sizes and distribution of Rivière-aux-Feuilles (RAF) and Rivière-George (RG) migratory caribou herds, northern Quebec and Labrador (Canada).
figure 2

Population sizes (a) were estimated based on post-calving photographic aerial surveys (Ministère des Forêts, de la Faune et des Parcs du Québec). Their confidence intervals were computed using a number of parameters recorded during surveys, and are inherent to the method (Des Meules and Brassard 1964; Le Hénaff 1976; Couturier et al. 1996, 2004; Ministère de Faune, des Forêts et des Parcs unpublished). We delineated ranges (b) using 100% minimum convex polygons based on ARGOS locations collected on 54 females for RG and 60 females for RAF in 2010.

Full size image

DNA extractions

We isolated genomic DNA from muscle (n = 361), blood (n = 12), and hair (n = 31) samples using the DNeasy Blood & Tissue Kit (Qiagen, Inc., Valencia, CA, USA) following the manufacturer’s protocol, with minor modifications for hair and blood samples (described in Gagnon et al. 2019). Genomic DNA concentrations were quantified using a Qubit 2.0 fluorometer (Life Technologies), and all samples standardized at 200 ng/µl before library preparation.

Illumina MiSeq library preparation and sequencing

We used Illumina MiSeq technology to sequence a 250-bp fragment of exon 2 of MHC–DRB gene (the complete exon has 267 bp in caribou; Kennedy et al. 2010) for 404 caribou. We amplified this locus with the LA31 and LA32 primer pair designed by Sigurdardóttir et al. (1991) for bovids, but successfully used in caribou from Alaska (Kennedy et al. 2010) and for Peary caribou (R. tarandus pearyi; Taylor et al. 2012). Kennedy et al. (2010) showed that this pair of primers amplifies only one copy of MHC–DRB in caribou. To multiplex individuals in the same final library, we modified forward and reverse primers by adding to the 5′-end of the sequence a unique barcode of 8 bp. We also added two to four random nucleotides (NN, NNN, or NNNN) to the barcodes to circumvent the low sequence diversity issue in amplicon sequencing with Illumina technologies (Mitra et al. 2015). We designed a total of 24 forward barcodes and 40 reverse barcodes, resulting in 960 unique combinations. Two or three independent PCR were performed for each individual, using different barcode combinations. We used 130 barcode combinations as negative controls consisting of PCR mix without DNA template to estimate tag switching and sequence–sample contamination. Each PCR reaction was performed in a total volume of 20 µl containing 16 µl of AmpliTaq GoldTM 360 Master Mix (ThermoFisher Scientific), 2 µl of a mix of two tagged primers (5 µM each), and 2 µl of template DNA (200 ng/µl). The thermocycling conditions were as follows: 10 min at 95 °C, 35 cycles of 30 s at 95 °C, 30 s at 55 °C, 45 s at 72 °C, 7 min at 72 °C, and hold at 4 °C. All reaction products were pooled together and purified using Promega Wizzard® SV Gel and PCR Clean-Up System following the protocol provided with the kit. The final library was prepared using the TruSeq Nano DNA HT Sample Prep Kit from Illumina using the MetaFast protocol (Fasteris SA, Geneva, Switzerland). It was then sequenced in a full run with a 250-bp paired-end MiSeq technology (Illumina).

MHC–DRB genotyping

Raw reads were first processed with the package OBITools (Boyer et al. 2016) to align forward and reverse reads, assign sequences to their individual sample, and dereplicate MHC fragments (see SI for complete details on those steps). We used R (v. 3.4.0, R Core Team 2017) to obtain genotypes for all individual replicates. Sequence variants shorter than 240 bp and those that were not associated with at least 20 reads in one sample were discarded. We discarded variants that were covered by less than 17% of the reads at the individual level as they were considered to be sequencing artifacts (see SI for details about how this threshold was defined). Finally, we compared replicated genotypes for mismatches.

MHC sequence analysis

We used MEGA7 (v. 7.0.18, Kumar et al. 2016) with the “ClustalW” algorithm to align sequences. Individual heterozygosity was coded as a binary variable describing whether MHC–DRB was heterozygous (coded as 1) or homozygous (coded as 0). Then, we calculated nucleotide diversity as the average number of substitutions between the two alleles (haplotypes) carried by an individual, i.e., the number of substitutions divided by the total sequence length (250 bp). Heterozygosity and nucleotide diversity were computed with DNAsp (v. 5.10.1, Librado and Rozas 2009). Similarly, we estimated amino acid (a.a.) diversity as the mean number of a.a. substitutions between the peptide sequences produced by the two alleles carried by an individual, i.e., the number of a.a. substitutions divided by the total length of the peptide (82 a.a.). We computed the number of a.a. substitutions using MEGA7. Furthermore, functional heterozygosity was coded as a binary variable describing whether an individual carried one (0) or two (1) different MHC–DRB supertypes.

Supertypes are functional groups of alleles that share similar physicochemical properties at their ABS (Sandberg et al. 1998; Doytchinova and Flower 2005; Buczek et al. 2016). To define these groups, we identified positively selected sites (PSS) in all sequences available on GeneBank for R. tarandus MHC–DRB (n = 39, including ours, after removing duplicates; accession numbers are presented in Table S1). We used the MEME (Murrell et al. 2012), FUBAR (Murrell et al. 2013), and FEL (Pond and Frost 2005a) methods implemented on the Datamonkey server (Pond and Frost 2005b; Delport et al. 2010; Weaver et al. 2018), with default parameters. For MEME and FEL methods, we used α = 0.1, and for FUBAR, we used a posterior probability threshold of 0.9. We also used PAML to compare two pairs of site models: M1a–M2a and M7–M8 using a likelihood-ratio test of positive selection with α = 0.05 (v. 4.8, Yang 2007). As the likelihood-ratio test revealed positive selection, the Bayes Empirical Bayes procedure was then used to identify sites under positive selection (Yang 2007). We retained all PSS that were identified by at least two methods (n = 16; Table S1). We then used the functions “find.clusters” and “dapc” of the R package adegenet (v 2.0.1, Jombart 2008; Jombart and Ahmed 2011) to cluster alleles based on the physicochemical properties (Sandberg et al. 1998) of the PSS following the procedure described by Buczek et al. (2016). We used the Bayesian Information Criterion to identify the optimal number of groups (K) above which increasing K did not lead to considerable improvement of the model fit (SI Fig. S4). We refer to the combination of supertypes carried by an individual as its functional genotype.

Finally, for each individual genotype, we computed the point accepted mutation (PAM) distance, a corrected distance between two nucleotide sequences that considers the likelihood of each amino acid substitution in the primary structure of the protein to be accepted by natural selection, given that amino acids have varying physicochemical properties (Dayhoff and Eck 1968; Dayhoff et al. 1978; Kosiol and Goldman 2005). This distance is considered to be the most “functional” index of MHC diversity along the gradient of genetic diversity depicted in Fig. 1. PAM distance was estimated with the function “Rprotdist” implemented in the R package RPhylip (Felsenstein 2013; Revell and Chamberlain 2013).

Relationship between MHC and genome-wide diversity

To verify our assumption that global genomic diversity would be more strongly correlated with allelic diversity indices compared with functional diversity indices, we compared genome-wide estimates of genetic diversity between homozygotes and heterozygotes for MHC alleles or supertypes. The standardized multilocus heterozygosity (sMLH) was derived from >22k SNPs (Gagnon et al. 2019) for most (n = 350) of the individuals comprised in this study. We used linear models with sMLH as a response variable and MHC–DRB heterozygosity or functional heterozygosity, herd, and the interaction between the two as fixed-effect variables. We conducted post hoc comparisons using the “lsmeans” function in the emmeans R package (Lenth 2019), with p value adjusted for multiple comparisons with the “Tukey” option. Because of the large differences in sample size between the RAF and RG herds, we equalized sample sizes resampling 1000 times without replacement of the same number of individuals in the two herds (n = 122) to test the correlation between genome-wide and MHC diversity. p values were estimated as the proportion of significant effects on 1000 samples.

Heterozygosity–fitness correlations

Body mass and survival were used as response variables in different sets of models testing the effect of MHC diversity on performance. We defined one set of models for each diversity index (heterozygosity, nucleotide diversity, amino acid diversity, functional heterozygosity, and PAM distance, Fig. 1), resulting in five sets of models, to allow for each index to be selected in the best model of its set. Each set included (1) a null model (with only random effects), (2) a baseline model (with only relevant, nongenetic variables, described below), (3) a simple effect model (with the baseline variables and the simple effect of the genetic diversity index), (4) a model with an interaction between the genetic diversity index and age in addition to the baseline variables, and (5) a model with an interaction between the genetic diversity index and sex in addition to the baseline variables. Nucleotide distance, amino acid diversity, and PAM distance were scaled between 0 and 1 to improve convergence. We defined an additional set of models to test for the effect of the functional genotype, by replacing the genetic diversity indices in the models described above by a categorical variable representing functional genotype. For all models, we ensured that the assumptions of normality and structure of residuals were respected, and that collinearity was not an issue (Zuur et al. 2010). We used Akaike’s information criterion for small-sample sizes (AICc) to rank models, and calculated ΔAICc and AICc weights (ωi) with the R package AICcmodavg (Mazerolle 2017). We selected the model with the lowest AICc value. We computed conditional and marginal R2 (cR2 and mR2 respectively) for the selected models using the “r.squaredGLMM” function implemented in the MuMIn R package (Bartoń 2019).

For the body mass models, we used linear mixed-effect models and the “lmer” function implemented in the package lme4 (Bates et al. 2015) of R (v. 3.4.0, R Core Team 2017). We log-transformed body mass to increase homogeneity of variances. We determined that the baseline variables for body mass models should include age, √age, sex, herd, and month of capture as fixed-effect variables because these factors are known to influence body mass in caribou (Parker 1981; Couturier et al. 2009a, b). As body mass does not increase linearly with age and forms a plateau after a few years, we included the √age term in the models to accurately model this relationship. We also included year of capture and individual identity as random factors to account for repeated measures and annual variation.

For the survival models, we fitted generalized mixed-effect models with a binomial-distributed error (logit link function), using the “glmer” function implemented in the package lme4 (Bates et al. 2015). Annual survival was binary coded as 1 if an animal was alive during a given year and as 0 the first complete year the animal was dead. The baseline variables for survival models included age, age2, sex, and herd as fixed-effect variables (Loison et al. 1999; Couturier et al. 2010), and individual identity and year as random-effect variables. The inclusion of age and age2 accounted for the nonlinear decrease in survival probability with age. We centered the age variable on its mean in order to avoid high collinearity between age and age2. Although body mass is often a predictor of survival, we did not include it as a fixed-effect variable in our survival models as individuals were not weighed at the same time (month and year) and at the same age, and because the reproductive status was unknown and could have affected mass (Chan-McLeod et al. 2000).

Assessing population structure and selection

We tested for temporal changes in allelic and supertypic frequencies over time, indicating potential fluctuating selection. We first defined cohort as a group of ≥5 individuals born the same year (Ortego et al. 2011). Consequently, not all years were included in the analysis: 13 cohorts for RAF and 12 cohorts for RG herd. Then, G tests implemented in the software GenoDive (v. 2.Ob27, Meirmans and van Tienderen 2004) were used to test for temporal differentiation between cohorts within each herd in terms of allelic or supertypic frequencies. We adjusted p values for multiple comparisons using the Benjamini and Yekutieli (2001) correction. We also computed observed and expected values of heterozygosity and functional heterozygosity, as well as Fis by cohort in GenoDive to assess departure from Hardy–Weinberg equilibrium (HWE). All tests mentioned here were based on 1000 permutations. The analysis described in this section were conducted using all captured individuals for which we had estimated the age on the field, including some calves that were not included in previous models testing the effect of MHC diversity on mass and survival.

Results

MHC–DRB genotyping

Only 0.067% of all raw reads were associated with negative control combinations of barcodes, although negative controls corresponded to 13.5% of all samples in our libraries. No control sample presented a sufficient number of reads to be retained after filtration, hence suggesting that tag switching and sequence–sample contamination were not an issue. We successfully genotyped 391 individuals out of 404 originally sequenced. All individual genotypes were supported by at least two concordant replicates with a mean number of reads of 8915 ± 3072 reads (SD) per replicate. Only one pair of replicates that had passed the filters was not concordant and was eliminated.

We identified 20 allelic sequences (or haplotypes, GenBank Accession numbers: MH171262–MH171281) out of which 2 and 4 were private alleles to RG and RAF herds, respectively (Tables 1 and S2; SI Fig. S5). MHC allele frequencies did not significantly differ between herds (SI Fig. S5). Nine alleles had already been identified in Grant’s Caribou (R. tarandus grantii) in Alaska (Kennedy et al. 2010), in Peary Caribou in Canadian Arctic (Taylor et al. 2012), and in other Reindeer populations (Mikko et al. 1999, Wei and Happ (unpublished, GenBank), Djakovic (unpublished, GenBank)). Eight alleles were variants (1–2-bp differences) or subtypes (3–4-bp differences) of previously identified alleles, and three were novel, unidentified alleles (SI Table S2). Among these alleles, 38 (15.2%) nucleotidic sites were polymorphic, and all alleles coded for different peptide chains, with 25 (30.1%) variable amino acids. We clustered the alleles in four supertypes that comprised 7, 5, 5, and 3 alleles each (SI Table S2). MHC supertype frequencies did not differ between herds (SI Fig. S6).

Table 1 Summary statistics of genetic diversity in Rivière-George (RG) and Rivière-aux-Feuilles (RAF) herds for both sexes (Females = F, Males = M).
Full size table

Relationship between MHC and genome-wide diversity

We did not find significant differences of genome-wide diversity (sMLH) between MHC–DRB homozygotes and heterozygotes (p value = 0.98), or between the functional homozygotes and heterozygotes (p value = 0.94). No other significant effects were found in the models.

Effect of MHC diversity on body mass

In all model sets predicting body mass (genetic diversity indices and functional genotype), the selected model was the baseline model (mR2 = 0.92; cR2 = 0.95; SI Tables S3 and S4), i.e., the one including age, √age, sex, herd, and month of capture as fixed-effect variables. Next-ranked models (SI Tables S3 and S4) presented no significant effect of genetic variables. In the baseline model, age, √age, sex, herd, and month of capture (June) had a significant effect on body mass. For each increase of 1 year in age, log(mass) decreased by −0.19 to −0.27, and for each increase of one unity of √age, log(mass) increased by 1.08–1.36. Roughly, this translated into a steep increase in weight for the first few years of life and a stabilization of the weight after 4 years around 95 kg. Males were heavier than females, and individuals from RG herd were heavier than RAF individuals. Individuals were generally heavier in Fall months compared with Winter or (early) Summer months. All estimates of the mass baseline model are listed in SI Table S5.

Effect of MHC diversity on survival

The baseline model (mR2 = 0.06, cR2 = 0.24) was selected for all sets of candidate models predicting survival (genetic diversity indices and functional genotypes), except for the one that included PAM distance (SI Tables S6 and S7). All next-ranked models of all sets did not present any significant effect of the genetic variables, except for one significant interaction between age and the functional genotype composed of supertypes 1 and 3 in the functional genotype × age model (ΔAICc = 13.65; ωi = 0). The selected model in the PAM distance set was the one including simple effects of PAM distance in addition to the baseline variables (ωi = 0.47; mR2 = 0.08; cR2 = 0.25; SI Table S6). In this model, only herd and PAM distance had a significant effect on survival probability (Fig. 3). Individuals from the RG herd had a probability of survival 0.41 times lower than that of those from the RAF herd (95% CI odd ratio (OR) = [0.25, 0.67], Fig. 3). In addition, individuals with maximal PAM distance in our dataset had a probability of survival 0.38 times lower than that of individuals with minimal PAM distance (95% CI OR = [0.18, 0.88], Fig. 3). Based on the predictions of the model, individuals with minimal PAM distance had a survival probability on average 15 ± 11% [SE] superior to that of individuals with maximal PAM distance (prediction for individuals of all ages, sexes, and herds, Fig. 4). Nongenetic variables followed the same trend in the baseline model, with similar statistical significance. All estimates of the survival baseline model and the survival PAM distance model are listed in Tables 2 and S8.

Fig. 3: Odds ratios and 95% confidence interval (CI) of the fixed-effect variables of the PAM distance model of migratory caribou survival (n = 215 individuals).
figure 3

Odds ratio 95% confidence interval that does not cross 1 is considered statistically significant. Levels of reference are Female for sex and Rivière-aux-Feuilles for herd. PAM distance model: Survival ~ age + age2 + sex + herd + functional heterozygosity + (1|year) + (1|identity).

Full size image
Fig. 4: Annual survival probability of migratory caribou according to PAM distance at MHC-DRB.
figure 4

The four pannels show annual survival probability of females (a, b) and males (c, d) from the Rivière-aux-Feuilles (a, c) and Rivière-George (b, d) herds according to their age and PAM distance. Survival probabilities of individuals aged 3–11 years were comprised between those of individuals aged 2 and 12 years. In the best model (ωi=0.47), PAM distance (estimate±SE: −0.95±0.40) and herd (estimate±SE: −0.89±0.25) had significant effects. Age, age2, and sex did not have a significant effect.

Full size image
Table 2 Estimates of the baseline model of survival (n = 215 individuals).
Full size table

MHC population structure and selection over time

We did not detect departure from HWE for the MHC alleles and supertypes in any of the cohorts (all p values > 0.05; SI Table S9). Multiple comparisons between cohorts of each herd revealed that there was no significant temporal change of allelic or supertypic frequencies over ~ 20 years of sampling (G test, α = 0.05; all p values equal to 1 after Benjamini and Yekutieli (2001) correction).

Discussion

We investigated the effect of immune gene variability (MHC–DRB) on performance-related traits in two declining herds of migratory caribou, and showed that various indices of genetic diversity correlated differently to indicators of performance. Body mass was not associated with MHC diversity, but we found a negative effect of PAM distance on both male and female survival. According to our model, the survival probability of individuals with the maximal PAM distance in the dataset was considerably lower than that of individuals with the minimal PAM distance (mean difference = 15 ± 11% [SE]).

Effects of MHC diversity on body mass

The lack of correlation between body mass and MHC–DRB diversity could be due to multiple factors, aside from the absence of a real effect. The link between MHC genotype and body mass is indirect, and the negative effects of infections or immune responses on body condition may be buffered by changes in resource allocation (Derting and Compton 2003) or increased food intake (Moret and Schmid-Hempel 2000; Lukasch et al. 2017a, b). In addition, ungulates usually present large variation in body mass over seasons (Parker et al. 1993; Chan-McLeod et al. 2000), and although we accounted for month of capture in the analyses, the correlation design we used was less powerful than sampling individuals of the same sex–age class at the same time. In our study, the reproductive status of individuals was unknown, which could have also influenced body mass, especially for females (Chan-McLeod et al. 2000). Such seasonal fluctuations, in conjunction with variance in the reproductive status and possible behavioral or physiological responses, may have obscured the potential effect of MHC diversity on body mass. Last but not least, like most phenotypic traits, body mass is a polygenic trait (Kemper et al. 2012; Fuller et al. 2020), meaning the effect of one single gene (let alone that of one single exon) is expected to be relatively low.

Effects of MHC diversity on survival

Our most striking finding is the negative association between PAM distance and annual adult survival probability (or negative HFC). Because natural selection acts on phenotypes and not on genotypes (Endler 1986), PAM distance is more likely to show an association with fitness-related traits as it is more closely related to the phenotype than heterozygosity or nucleotide diversity, for instance. As PAM distance was the only index that predicted survival in our models, we conclude that direct effects caused the observed relationship.

The negative HFC we observed here is surprising as higher diversity at immunity-related genes is usually associated with higher performance (e.g., de Assunção-Franco et al. 2012; Osborne et al. 2015; Portanier et al. 2019), especially when general effects can be ruled out. Some studies have shown, however, that lower MHC diversity can be beneficial to individuals. For instance, Ilmonen et al. (2007) found that mice (Mus musculus domesticus) that were heterozygous at MHC genes responded poorly to Salmonella (Salmonella enterica) infections: they suggested that it could be due to a gene-dose effect that conferred an advantage to homozygotes or to depleted T-cell repertory, and reduced immunocompetence in individuals with maximum MHC diversity. Major differences between Ilmonen et al. (2007) experimental study design and our correlative study design make it unlikely that these mechanisms were involved in the relationship observed in our study. Indeed, it is likely that the gene-dose effect suspected by Ilmonen et al. (2007) would be observable only when hosts are infected with a single pathogen species, whereas in our study, individuals were most likely coinfected by multiple pathogen species (Pollock et al. 2009; Steele et al. 2013; Simard et al. 2016). Furthermore, the heterozygote mice studied by Ilmonen et al. (2007) had maximal level of heterozygosity at all class I and II MHC genes, making it plausible that these individuals suffered from reduced T-cell repertory (Nowak et al. 1992). Although we have not genotyped the whole caribou MHC genes, it is unlikely that wild individuals would present such high levels of MHC diversity, especially considering that caribou has very infrequent copy number variation in MHC genes (no known copy number variation in DRA, DRB, and DQB and potential copy number variation in DQA; Kennedy et al. 2010).

Outbreeding depression is another mechanism that could result in MHC homozygote advantage (Bos et al. 2009), if general effects were involved. Indeed, outbreeding can break up good combinations of alleles (e.g., adapted to local conditions) and result in reduced performance in outbred individuals (Bonneaud et al. 2006; Freeland et al. 2011). Here, the fact that the allelic diversity indices of MHC diversity (i.e., heterozygosity and nucleotide diversity) were not correlated with either genome-wide diversity or performance suggests that general effects were not involved. In addition, we found no evidence of such process affecting the populations studied here in previous work, based on the absence of correlation between genome-wide diversity and performance (Gagnon et al. 2019).

Last, the apparent homozygote advantage can result from single alleles associated with better performance (Charbonnel et al. 2010; Lukasch et al. 2017a, b) or a combination of alleles associated with worst performance (McClelland et al. 2003; Acevedo-Whitehouse et al. 2018) that increase mean homozygote relative performance. Here, no functional genotype showed any association with higher or lower performance, suggesting that the negative association between MHC diversity and survival was not mediated by one supertype. Nevertheless, Lukasch et al. (2017a, b) showed that highlighting the effect of a single MHC allele in a multiple-infection context can be complex, as each allele can have positive and negative effects in response to different immunological challenges.

The mechanism linking caribou MHC–DRB genotype and survival remains unclear, although the results suggest that it most likely originates from direct effects. Assuming that individuals with high functional diversity have higher immunocompetence than individuals with low functional diversity, we can hypothesize that they experience higher costs of immunity. Acevedo-Whitehouse et al. (2018) discussed this hypothesis, arguing that if higher MHC diversity translates into better antigen recognition and higher antibody production, then it should also increase the risks of immunopathology (i.e., damages caused by the host immunity to the host tissues). In this context, caribou with high MHC functional diversity would not only benefit from the recognition of a broader array of pathogen (divergent allele-advantage hypothesis) but also suffer the costs associated with more frequent or more intense immune responses. This idea is in line with observations made in European badgers (Meles meles; Sin et al. 2016), suggesting that total investment in immune functions (acquired and innate) was positively correlated with MHC diversity. Furthermore, it was suggested that relative investment in innate and acquired immune functions might be mediated by infection prevalence in a way that balances the costs and benefits of maintaining these functions (Moret 2003). Possible changes in pathogen communities due to climate change (Kutz et al. 2005; Altizer et al. 2013; Dobson et al. 2015) could have shifted the optimal investment in immune functions to the detriment of individuals with higher functional diversity at MHC.

We recognize that our sample sizes for each cohort were probably too small to detect selection, rare alleles (Allendorf and Luikart 2007), and departure from Hardy–Weinberg. This weakness could explain why we detected no departure from the HWE even though we would have expected the levels of genetic diversity to be lower than under neutrality, given the disadvantage of individuals with high PAM distance. In addition, because not all individuals were captured at their first year of life, selection may have acted before our sampling, making it more difficult to detect afterward. The identification of multiple PSS at MHC–DRB and the maintenance of allelic diversity during the 20 years of sampling (~5 generations), despite the strong demographic decline, suggest, however, that MHC–DRB diversity was maintained by balancing selection. As survival is only one component of fitness, trade-offs between life-history traits such as survival and reproduction, could conceal a more complex relationship between MHC diversity and fitness. For instance, Lo et al. (2015) showed that reproductive females of the Taiwan field mouse (Apodemus semotus) were more infected with parasites than nonreproductive ones, and that females treated with an anthelmintic drug had reduced survival. The authors suggested that treated females were in better condition and may have invested more in their current reproduction at the expense of their own survival and future reproduction. It is not impossible that such immunity-mediated trade-off between survival and reproduction could occur in RG and RAF herds and explain that positive selection occurs, despite lower survival in individuals presenting high MHC functional diversity. Until further studies investigate the relationship between immune gene diversity, survival probability, and reproductive success in migratory caribou, these results should be treated with caution.

Our results suggest that MHC diversity could have a sufficiently large impact on adult survival to also affect population dynamics. Adult survival is usually very stable in large ungulates, and variation in this parameter can have profound impacts on population trends (Gaillard et al. 2000). The predicted effect of PAM distance on survival probability (15 ± 11% [SE]) could be large enough to make the difference between a stable and a fast-declining population (Gaillard et al. 2000). It could support the conclusions of Albon et al. (2002) who found that pathogens could play an important role in population dynamics of caribou. To confirm this hypothesis, further studies should investigate the link between MHC diversity and more direct proxies of immunocompetence (e.g., performance at multiple immune-challenge assays) in caribou.

Overall, our research highlights the importance of using various indices of genetic diversity when investigating the effects of genetic diversity on performance, as this may lead to different outcomes. Our results also open new perspectives for future research on the genetic interaction between host survival, immunity, and pathogens, and on the mechanisms that maintain genetic variability at MHC genes in natural declining populations.