- Research article
- Open Access
Alternated selection mechanisms maintain adaptive diversity in different demographic scenarios of a large carnivore
BMC Evolutionary Biologyvolume 19, Article number: 90 (2019)
Different population trajectories are expected to impact the signature of neutral and adaptive processes at multiple levels, challenging the assessment of the relative roles of different microevolutionary forces. Here, we integrate adaptive and neutral variability patterns to disentangle how adaptive diversity is driven under different demographic scenarios within the Iberian wolf (Canis lupus) range. We studied the persistent, the expanding and a small, isolated group within the Iberian wolf population, using 3 MHC class II genes (DRB1, DQA1, and DQB1), which diversity was compared with 39 microsatellite loci.
Both the persistent and the expanding groups show evidence of balancing selection, revealed by a significant departure from neutrality at MHC loci, significant higher observed and expected heterozygosity and lower differentiation at MHC than at neutral loci, and signs of positive selection. However, despite exhibiting a significantly higher genetic diversity than the isolated group, the persistent group did not show significant excess of MHC heterozygotes. The expanding group, while showing a similar level of genetic diversity than the persistent group, displays by contrast a significant excess of MHC heterozygotes, which is compatible with the heterozygote advantage mechanism. Results are not clear regarding the role of drift and selection in the isolated group due to the small size of this population. Although diversity indices of MHC loci correspond to neutral expectations in the isolated group, accelerated MHC divergence, revealed by a higher differentiation at MHC than neutral loci, may indicate diversifying selection.
Different selective pressures were observed in the three different demographic scenarios, which are possibly driven by different selection mechanisms to maintain adaptive diversity.
Genetic drift and natural selection are often seen as two counteracting evolutionary forces across populations . In large and stable populations, where the random effects of genetic drift are limited, balancing selection maintains high levels of adaptive variation affecting population fitness and evolutionary potential [2, 3]. By contrast, in small and isolated populations, or within the context of demographic events such as bottlenecks, genetic drift is predominant, masking possible effects of natural selection [4, 5]. Neutral molecular markers allow assessing interactions between demographic history and genetic diversity , but do not provide direct information on selective processes [7, 8]. Thus, evidence of selection and adaptive potential is best approached by contrasting neutral and functional molecular markers (e.g., [9, 10]).
In vertebrates, the major histocompatibility complex (MHC) is the prime candidate for pathogen resistance genes and contains some of the most polymorphic functional loci [11, 12]. The exceptional level of MHC polymorphism is believed to be driven by the antagonistic coevolution with pathogens and occurs through pathogen-mediated balancing selection [11, 12]. Beyond its clear significance in modulating pathogen resistance [2, 13, 14], MHC genes have been shown to influence other biological traits such as maternal-fetal interactions, kin recognition, life-time reproductive success and mate choice [15,16,17,18]. MHC is intimately linked with factors likely to affect individual fitness, population viability and evolutionary potential in changing environments. Thus, patterns of MHC diversity have been repeatedly used in a conservation context in populations of particular interest [19,20,21,22].
In Europe, abundance and distribution of large carnivores (bear, wolf and lynx) have been dramatically shaped by humans during the past few centuries. Although persecuted and driven to or close to extirpation in several countries, most of the remnant populations have stabilized or even expanded recently . One example is the Iberian wolf (Canis lupus), which erstwhile ranged over most of the Iberian Peninsula, but after decades of severe human persecution became confined to small and fragmented populations, representing about 1/5 of the former population range [23,24,25]. In the 1970s, the Iberian wolf population was estimated to be reduced to ca. 700 individuals, mostly restricted to the northwestern region [23, 26, 27]. After the 1970s, the population started to grow and expanded southwards and eastwards (Fig. 1) . Currently, the Iberian wolf is considered under a stable or positive demographic trend . It occurs in a core area where the species have always persisted in northwestern Iberia [26, 27, 29], and the adjacent re-colonized area (the expansion front), alongside two small and isolated populations [30,31,32]. Such a complex demographic history have translated into a cryptic population structure in the Iberian wolf at small spatial scales, with moderate level of genetic differentiation, including the differentiation of the re-colonized area .
Remarkable MHC class II variability has been shown in wolf populations [34,35,36,37]. However, results vary regarding the predominant role of selection over neutral forces. Some studies provided evidence that wolf MHC diversity is maintained by balancing selection, including studies on populations under demographic decline [34, 36], where parasite resistance was suggested as the possible driving force . In contrast, MHC diversity in the bottlenecked wolf population of Scandinavia was shown to be compatible with neutral evolution . Although MHC variability patterns of small and isolated populations may differ from that in large and outbred populations, the power of balancing selection acting on MHC can be outweighed by demographic events, such as bottlenecks and fragmentation with consequent genetic drift [38, 39]. Hence the signature of selection and/or population demographic effects is expected to vary across populations under different demographic histories, challenging the assessment of the relative role of different microevolutionary forces .
In this study, we integrate neutral and adaptive diversity patterns to disentangle how adaptive diversity is driven under three different demographic within the Iberian wolf range. Taking into account the genetic population structure of this population based on neutral markers , we used the permanent wolf range (i.e., persistent group) as the baseline for demographic stability and compared it to the population dominated by the re-colonization process after the 1970s bottleneck  (i.e., expanding group) and to the remaining small and isolated population in the South of Douro River in Portugal (i.e., isolated group). For the persistent group we expect to find natural selection acting to maintain adaptive variation at a stronger intensity than genetic drift; and therefore MHC diversity is expected to exceed neutral diversity. In this scenario, we can either detect balancing or diversifying selection. We may observe balancing selection when MHC exhibits lower population differentiation than neutral loci, due to enhanced effective migration rate and similar allele frequencies in MHC loci . Alternatively, we can detect diversifying selection when MHC exhibits accelerated divergence as a response to spatial or temporal heterogeneity in parasite abundance and diversity . Contrastingly, for expanding and isolated groups, genetic drift is expected to outweigh selection; and therefore similar diversity patterns between MHC and neutral markers are expected.
We found seven DRB1, four DQA1, and six DQB1 alleles across the 113 Iberian wolf samples (Additional file 1: Table S1), all previously reported in other wolf populations (Additional file 2: Table S2, Additional file 3: Table S3). Diversity at DRB1 and DQB1 was higher than diversity at DQA1 (Additional file 4: Table S4). Observed heterozygosity was higher than expected at all three loci considering the whole dataset, with significant departure from Hardy-Weinberg Equilibrium (Additional file 4: Table S4). Nucleotide and amino acid distances were higher for DRB1 and DQB1 than for DQA1, with amino acid distances higher than the nucleotide distances (Additional file 4: Table S4).
When clustering the three loci (DRB1/DQA1/DQB1), we observed seven haplotypes (Fig. 1, Table 1). Each DRB1 allele is found in haplotypic association with a different DQ pair, with exception of alleles DRB1*05401 and DRB1*05501, which share the same DQ pair, DQA1*00301/ DQB1*00401.
Five haplotypes are shared with other European wolf populations (Additional file 3: Table S3). We cannot confirm the presence of two haplotypes in other European wolf populations due to the absence of haplotype reconstruction in , where all alleles that compose these two haplotypes were reported. The three most frequent three-locus haplotypes are shared with Croatia, Italy and Finland (Fig. 2, Additional file 2: Table S2, Additional file 3 Table S3). Iberian haplotypes are distributed throughout the NJ tree (Fig. 2). Seven, six and four three-locus haplotypes were found in the persistent, expanding and isolated groups, respectively (Fig. 1, Table 1). The most frequent haplotype was different for each demographic group (Fig. 1, Additional file 2: Table S2). Sequence diversity as the number of segregating sites, nucleotide diversity, number of mutations and Watterson’s mutation parameter were similar among the three different demographic groups (Table 1).
All MHC loci revealed positive and significant departure from neutrality (Additional file 4: Table S4), indicating an excess of high-frequency segregating sites. Persistent and expanding groups revealed positive and significant departure from neutrality for Tajima’s D and Fu & Li D* tests (Table 1). The isolated group did not show significant departure from neutrality for Tajima’s D and Fu & Li D* tests (Table 1). The dN/dS ratio was significantly different from neutral expectations in all demographic groups both at all sites and at the peptide binding region for DRB1 loci, and only at all sites for DQB1 loci (Table 2). For DQA1, the dN/dS ratio was significantly different from neutral expectations at all sites in isolated group and at the peptide binding region in persistent and expanding group (Table 2).
Additionally, we found several MHC codons under positive selection in the three demographic groups using OmegaMap, a few having high posterior probability (PP > 95%, Additional file 5: Table S5). DRB1 and DQB1 loci had higher proportions of selected codons than DQA1 (Additional file 5: Table S5). We matched our alignment with the human orthologue  and observed that almost all selected codons with high posterior probability (DRB1: 6 in 7, DQA1: 1 in 2, DQB1: 4 in 4 codons) match the same peptide-binding sites described by these authors (Additional file 5: Table S5).
Genetic differentiation between demographic groups
We found a single model with support from the data to explain allelic richness and expected heterozygosity, with a deviance explained of more than 65% (Additional file 6: Table S6). Allelic richness is explained by a combination of demographic group and locus type, while expected heterozygosity is explained by demographic group, both including the effect of locus (Additional file 6: Table S6). For both the number of alleles and observed heterozygosity we found three models with support, presenting a deviance explained between 20 and 43% (Additional file 6: Table S6). Number of alleles is explained by demographic group or a combination of demographic group and locus type (Additional file 6: Table S6). Observed heterozygosity is explained by the null model or by locus type (Additional file 6: Table S6).
According to the selected models, allelic richness, expected heterozygosity and the number of alleles are significantly higher in the persistent than in the isolated group (Additional file 7: Table S7). Number of alleles is also significantly higher in the persistent than in the expanding group (Additional file 7: Table S7). Allelic richness and observed heterozygosity are significantly higher in MHC than in microsatellites.
Comparing MHC and neutral diversity
For both MHC and microsatellite loci, significant deviations from Hardy-Weinberg Equilibrium were observed for the persistent and expanding groups (Table 3). However this deviation is in opposite directions (Table 3). The persistent group showed significant higher expected than observed heterozygosity in both molecular markers, while the expanding group showed significant higher observed than expected heterozygosity in MHC (Table 3). Observed heterozygosity was significantly higher at MHC than at microsatellite loci for the persistent and expanding groups (p < 0.05, t-test, Additional file 8: Table S8). For the expanding group, the mean allelic richness and expected heterozygosity at MHC were also significantly higher than the ones at microsatellite loci (p < 0.05, t-test, Additional file 8: Table S8). For the isolated group, none of the diversity indices was significantly different between MHC and microsatellite loci (p > 0.05, t-test, Additional file 8: Table S8).
The first two components of PCA based on MHC loci did not differentiate among the three demographic groups (Fig. 3). Pairwise FST values using MHC data support the lack of genetic differentiation between persistent and expanding groups, while moderately significant differentiation is observed between all groups using Jost’s D differentiation (Table 4). Contrastingly, PCA based on microsatellite data showed clear genetic differentiation between all demographic groups, especially for the isolated group, while a partial overlap is observed between persistent and expanding groups (Fig. 3). All FST and Jost’s D pairwise comparisons between groups were significant for genetic differentiation using microsatellite loci (Table 4). Overall, pairwise genetic differentiation using FST was lower in MHC than in microsatellite loci, whereas for Jost’s D differentiation index, the isolated group exhibited higher values for MHC than for microsatellite loci (Table 4). The global co-inertia RV coefficient observed was high (0.98), indicating strong correlation between MHC and microsatellite matrices, but this coefficient is not different from what is expected by chance (p-value = 0.17).
This is the first report of MHC class II diversity in the Iberian wolf. Overall, Iberian wolves exhibit lower MHC diversity than their European counterparts, both in number of alleles and number of three-locus haplotypes (i.e., 6–7 alleles and 7 haplotypes in Iberia vs 6–13 alleles and 13–14 haplotypes in other European populations) [34,35,36]. The exception is the isolated Scandinavian wolf population, which went through a drastic bottleneck recovering from only three individuals . The reduced MHC diversity in the Iberian wolf has also been observed for other genomic regions [31, 43,44,45]. Long-term isolation and past bottlenecks in the Iberian population [27, 45] have certainly reduced neutral diversity and may also be at the origin of this depleted MHC diversity.
All MHC alleles have been previously identified in other canids [37, 46, 47], a pattern that is better explained by the trans-species polymorphism described in MHC canid phylogenies [37, 48] than by hybridization which is not common in the Iberian wolf . The observed 17 alleles were confined in just seven three-locus haplotypes of which five have previously been reported [34,35,36, 48]. The remaining two haplotypes were not unequivocal confirmed to be present in other wolf populations due to the absence of haplotype reconstruction in . The occurrence of conserved three-locus MHC haplotypes in Iberian wolves confirms the tight linkage of these gene loci, supporting a strong selective pressure maintaining haplotype combinations. The association observed between each DRB1 allele and a specific DQ pair supports the preferential association between alleles at these loci . A similar trend is observed in other European wolf population [34,35,36], though not as extreme as in Iberia, where seven DRB1 alleles yielded only seven haplotypes. However, similarly to other European populations, the seven haplotypes are widely distributed in the NJ tree, supporting maximal MHC diversity in the Iberian wolf.
Within Iberia, we found contrasting MHC diversity patterns according to different demographic trajectories. The persistent group exhibited significant higher genetic diversity than the isolated group, but did not show significant excess of MHC heterozygotes. However, we found evidence of balancing selection acting to maintain MHC adaptive variation in this group, compensating the effects of the severe population decline in the 19th and 20th centuries [50, 51]. This is revealed by a significant departure from neutrality at MHC loci, significant higher observed and expected heterozygosity and lower differentiation at MHC than neutral loci, and signs of positive selection with higher proportions of selected codons lying in peptide-binding regions, as inferred from humans. Similar results have been reported for Finnish wolves where, despite a strong genetic signal of population decline , the population still reveals a predominant role of balancing selection maintaining MHC diversity . Although selective advantage has been associated with a heterozygote advantage mechanism, MHC diversity can also be promoted by negative frequency-dependent selection . This model proposes that the emergence of a novel pathogen can promote the frequency of a rare or novel allele/haplotype [2, 3, 53]. The frequency shift of an adaptive allele/haplotype by such process can occur in very few generations . We can therefore not exclude that negative frequency-dependent selection may have contributed to maintain high allelic richness possibly compensating for the low number of heterozygotes observed.
The expanding group exhibited similar genetic diversity than the persistent group but show significant excess of MHC heterozygotes. Reduced genetic diversity is expected for marginal populations [10, 54,55,56], but increased heterozygosity at MHC supports a model of heterozygote advantage for this demographic group. The positive and significant result for neutrality statistics provides further evidence for balancing selection acting in this demographic group. The overdominance model postulates that heterozygotes have higher fitness than homozygotes due to the wider spectrum of MHC receptors able to induce parasite resistance . MHC heterozygotes have indeed been associated with resistance to several infections in wolves . Thus, the high MHC genetic diversity exhibited by the expanding group is expected to maintain its high adaptive potential.
The isolated group showed the lowest neutral and adaptive diversity and no significant excess of MHC heterozygotes, a likely result of the dramatic population decline that culminated in isolation in the early twenty-first century [26, 30]. The observed depletion of MHC diversity potentially increases population susceptibility to disease, and adds to current concern for its survival [12, 57]. Indeed, in a wolf serologic survey the isolated group showed no canine parvovirus antibodies, possibly related to high case-fatality rate . No significant departure from neutrality neither significant difference between MHC and neutral loci were observed, suggesting that MHC loci behave according to neutral expectations in the isolated group. Other examples are known for small and isolated populations where genetic drift outweighs the strength of selection [3, 38, 39]. However, one interesting result comparing the isolated group with the others may suggest selective forces acting in this group. The higher Jost’s D differentiation index observed for MHC relative to neutral loci may be explained by diversifying selection as consequence of spatial variation in pathogen-mediated responses  and/or due to stronger divergence due to short-term positive selection on beneficial alleles. Interestingly, the isolated wolf population shares several pathogens with wild and domestic species [59,60,61], but holds also new pathogens , which could be associated with the observed differentiation in MHC. Inconclusive results are probably due to the small sample size of this group and should be interpreted with caution .
By empirically testing adaptive diversity under three different demographic scenarios within the Iberian wolf range, we were able to show that these different demographic scenarios are associated with different selective pressures, which are possibly driven by different selection mechanisms to maintain adaptive diversity. Future studies are warranted to investigate the link between the Iberian wolf immunogenetic diversity and the pathogen community structure to better understand the mechanisms underlying adaptation in each demographic group.
Sampling and DNA extraction
A total of 113 wolf tissue samples were collected across the range of wolves in the Iberian Peninsula in the last two decades. Tissue samples were collected from several Portuguese and Spanish entities (please see Ethics section). To avoid close familiar relationships and the presence of hybrids between wolf and dog in our dataset we selected samples according criteria defined by . Samples were classified into one of three demographic groups: i) persistent group (n = 78), represented by samples from the northwest Iberian Peninsula where the wolf has been continuously present in the last decades, ii) expanding group (n = 25), represented by samples originally from the area where wolves were absent during the recent low of the population in the 1970s, representing a well-defined genetic cluster within Iberia (Castilla y León cluster) , and iii) isolated group (n = 10), represented by samples from the small and isolated wolf area south of the Douro river in Portugal (Fig. 1). Samples were stored in ethanol and DNA was further extracted using QIAGEN DNeasy Blood & Tissue Kit.
MHC amplification and sequencing
DLA class II genes DRB1, DQA1 and DQB1 were amplified and sequenced using intronic and locus specific primers described in . Polymerase chain reactions (PCR) were carried out in a total volume of 10 μl containing approximately 10 ng of DNA, 0.4 μM of each primer and 5 μl 2× MyTaq HS Mix (Bioline) (see Additional file 9: Table S9 for PCR conditions). PCR products were sequenced in both directions using BigDye® Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems) according to the manufacturer’s protocol and sequencing products were separated on a 3130xl Genetic Analyzer (Applied Biosystems).
MHC alleles and haplotype assignment
MHC sequences were phased in DnaSP v. 5.10  using PHASE  with the “recombination” model (−MR0) and 1000 iterations after 100 burn-ins. The threshold for accepting true alleles was set to 90%. This methodology has been shown to be accurate for MHC haplotype reconstruction , retrieving very low false positives . MHC alleles obtained for each locus were compared with available information in public datasets (NCBI: https://www.ncbi.nlm.nih.gov) for the genus Canis. As the three genes studied are chromosomally clustered in tandem, we analyzed each locus individually or as three-locus haplotypes. Three-locus haplotypes (DRB1/DQA1/DQB1) were reconstructed firstly in individuals homozygous for each locus, and then for heterozygous individuals basing inferences on haplotypes observed in homozygous state. Haplotype reconstruction was confirmed with available information for other canid populations.
MHC diversity was analysed for each locus and for the three-locus haplotypes. For both we estimated number of alleles (NA), mean number of effective alleles (Ne), mean observed (Ho) and expected (He) heterozygosities and inbreeding coefficient (FIS) using GENALEX 6.5 . Test of deviation from Hardy-Weinberg equilibrium was calculated using GENEPOP 4.6 . Allelic richness (AR) was estimated using FSTAT 220.127.116.11 , which performs rarefaction compensating sampling disparity, reducing all populations to a minimum sample size of n = 10. Sequence-level diversity values, such as number of haplotypes (h), number of segregating sites (S), nucleotide (π) and haplotype (HD) diversities, and Watterson’s mutation parameter (θw) were estimated either in each locus or each group (with the three-locus haplotypes) using DnaSP 5.10.
Molecular divergence analysis over all sequence pairs, also either in each locus or each group (with the three-locus haplotypes) was conducted using Mega 7  with a 1000 replicate bootstrap procedure. The best nucleotide and amino acid substitution models were chosen according to the Bayesian Information Criterion and used to compute nucleotide and amino acid evolutionary distances by the maximum-likelihood method through 1000 bootstrap replicates. The relationship between three-locus haplotypes of Iberian and other European wolf populations was reconstructed using the Neighbor Joining (NJ) method and p-distances in Mega 7. Three-locus haplotypes of other European wolf populations were retrieved from previous studies [34,35,36, 38].
Testing positive selection
To detect evidence of selection on MHC loci at the population-level, we estimated Tajima’s D  and Fu & Li D*  statistics using DnaSP 5.10. Under neutrality, D and D* are not different from 0, and significant negative deviations may indicate a recent selective sweep or population expansion while positive deviations indicate balancing selection or population contraction. Neutrality tests were computed for each MHC locus and for the three-locus haplotypes in the whole Iberian wolf range and in each demographic group. Deviations from neutrality were also assessed by a Z-test (non-synonymous substitution ≠ synonymous substitution; dN ≠ dS) implemented in MEGA 7, using Nei-Gojobori method with Jukes-Cantor correction for multiple substitutions and 1000 bootstrap replicates. Z-tests were computed for all sites, antigen-binding sites inferred from the human orthologue , and non-antigen binding sites for each demographic group. Additionally, we also used a “coalescent with recombination” Bayesian approximation implemented in OmegaMap  to find signatures of selection in each demographic group. OmegaMap allows both selection parameter (ω) and recombination rate (ρ) to vary along the sequence . Codons with excess of non-synonymous polymorphism compared to synonymous polymorphism are assumed under positive selection. We followed  to set distribution and parameter values of priors to represent neutrality. We performed two runs for each locus and each demographic group with 250,000 repeats and a burn-in of 20,000. Convergence of runs was check through summary function implemented in R 3.5.0 (R Core Team, 2018). Codons are considered as positively selected with posterior probabilities (PP) above 95%.
We selected a set of 39 microsatellite loci previously genotyped by  for all but three samples (two from the persistent and one from the isolated groups) sequenced for MHC to estimate neutral genetic diversity in each demographic group. This set included only loci with less than 25% missing data and with no significant departure from Hardy–Weinberg equilibrium after Bonferroni correction estimated using GENALEX 6.5. Genetic diversity, including number of alleles (NA), mean number of effective alleles (Ne), mean allelic richness (AR), mean observed (Ho) and expected (He) heterozygosities and inbreeding coefficient (FIS) were estimated using GENALEX 6.5 and FSTAT 18.104.22.168. For allelic richness all populations were reduced to a minimum sample size of n = 5. Test of deviation from Hardy-Weinberg equilibrium was calculated using GENEPOP 4.6.
Testing for genetic differentiation between demographic groups
To test for statistically significant differences of genetic diversity between the demographic groups, we performed linear mixed effects models (e.g.,[76, 77]) using lme4  and lmerTest  implemented in R platform (R Core Team, 2018). We used number of alleles, allelic richness, observed and expected heterozygosities as dependent variables, demographic group and locus type (MHC and microsatellites) as explanatory factors, and locus as random variable. Ten competing models were tested, including i) null, ii) demographic group-dependent, iii) locus type-dependent, iv) dependent on both demographic group and locus type, and v) dependent on both demographic group and locus type plus the association between these two factors. For each of these options we run two models accounting or not for locus random effect (for model syntaxes see Additional file 6: Table S6). We performed independent model testing for each diversity measure. We used Akaike Information Criterion corrected for small sample sizes, AICc, to determine the best model. Models with delta AICc (ΔAICc) lower than 2.0 were assumed as equally plausible. For each diversity measure we present estimates, standard errors and probability inferred from t-value, indicating the precision of the estimates.
Comparing neutral and functional diversity
Expectations that functional diversity exceeds neutral diversity were tested by a one-tailed t-test (p-value), following . The relationship between microsatellite and MHC loci was estimated using a co-inertia analysis (CoA) following . CoA is a multivariate method that estimates the covariance structure between two tables with the same observations (e.g., same individuals). For each genetic matrix (microsatellite and MHC loci) we performed a Principal Component Analysis (PCA) followed by factorial PCA using demographic groups as explanatory variable. For this analysis we used the adegenet and ade4 R packages [80, 81]. We then performed a CoA using the two most important PCA components. In a CoA plot, an arrow and a dot define each demographic group, where the length of the vector connecting arrow and dot is proportional to the divergence between the types of markers. The correlation coefficient RV was used to evaluate the strength of the relationship between the two matrices. Values of RV close to 1 indicate strong correlation between matrices. The significance of the RV coefficient was tested using 1000 permutations, where null hypothesis is that the two matrices are related by chance. Canonical weights for microsatellite and MHC alleles were used to disclose the contribution of each allele to the divergence of demographic groups.
MHC and neutral differentiation between groups was estimated through FST  and Jost’s D  using Arlequin 22.214.171.124  and the DEMEtics R package , respectively. FST and Jost’s D differentiation measures were both used because FST may be affected by highly variable markers .
Akaike Information Criterion corrected for small sample sizes
Mean allelic richness
- FIS :
number of haplotypes
Major histocompatibility complex
Number of alleles
Mean number of effective alleles
Principal Component Analysis
Polymerase chain reaction
Number of segregating sites
Watterson’s mutation parameter
Hartl DL, Clark AG. Principles of population genetics. Sunderland: Sinauer Associates; 1997.
Sommer S. The importance of immune gene variability (MHC) in evolutionary ecology and conservation. Front Zool. 2005;2:1–18.
Sutton JT, Nakagawa S, Robertson BC, Jamieson IG. Disentangling the roles of natural selection and genetic drift in shaping variation at MHC immunity genes. Mol Ecol. 2011;20:4408–20.
Ejsmond MJ, Radwan J. MHC diversity in bottlenecked populations: a simulation model. Conserv Genet. 2011;12:129–37.
Eimes JA, Bollmer JL, Whittingham LA, Johnson JA, van Oosterhout C, Dumm PO. Rapid loss of MHC class II variation in a bottlenecked population is explained by drift and loss of copy number variation. J Evol Biol. 2011;24:1847–56.
Holderegger R, Kamm U. F G. adaptive vs. neutral genetic diversity: implications for landscape genetics. Landsc Ecol. 2006;21:797–807.
Brumfield RT, Beerli P, Nickerson DA, Edwards SV. The utility of single nucleotide polymorphisms in inferences of population history. Trends Ecol Evol. 2003;18:249–56. https://doi.org/10.1016/S0169-5347(03)00018-1.
van Tienderen PH, de Haan AA, van der Linden CG, Vosman B. Biodiversity assessment using markers for ecologically important traits. Trends Ecol Evol. 2002;17:577–82. https://doi.org/10.1016/S0169-5347(02)02624-1.
Agudo R, Alcaide M, Rico C, Lemus JA, Blanco G, Hiraldo F, et al. Major histocompatibility complex variation in insular populations of the Egyptian vulture: inferences about the roles of genetic drift and selection. Mol Ecol. 2011;20:2329–40.
Rico Y, Ethier DM, Davy CM, Sayers J, Weir RD, Swanson BJ, et al. Spatial patterns of immunogenetic and neutral variation underscore the conservation value of small, isolated American badger populations. Evol Appl. 2016;9:1271–84.
Bernatchez L, Landry C. MHC studies in nonmodel vertebrates: what have we learned about natural selection in 15 years? J Evol Biol. 2003;16:363–77.
Ujvari B, Belov K. Major histocompatibility complex (MHC) markers in conservation biology. Int J Mol Sci. 2011;12:5168–86.
Hedrick PW, Lee RN, Buchanan C. Canine parvovirus enteritis, canine distemper, and major histocompatibility complex genetic variation in Mexican wolves. J Wildl Dis. 2003;39:909–13. https://doi.org/10.7589/0090-3558-39.4.909.
Schwensow N, Fietz J, Dausmann KH, Sommer S. Neutral versus adaptive genetic variation in parasite resistance: importance of major histocompatibility complex supertypes in a free-ranging primate. Heredity. 2007;99:265–77. https://doi.org/10.1038/sj.hdy.6800993.
Edwards SV, Hedrick PW. Evolution and ecology of MHC molecules: from genomics to sexual selection. Trends Ecol Evol. 1998;13:305–11. https://doi.org/10.1016/S0169-5347(98)01416-5.
Piertney SB, Oliver MK. The evolutionary ecology of the major histocompatibility complex. Heredity. 2005;96:7–21. https://doi.org/10.1038/sj.hdy.6800724.
Kalbe M, Eizaguirre C, Dankert I, Reusch TBH, Sommerfeld RD, Wegner KM, et al. Lifetime reproductive success is maximized with optimal major histocompatibility complex diversity. Proc R Soc B Biol Sci. 2009;276:925–34. https://doi.org/10.1098/rspb.2008.1466.
Sin YW, Annavi G, Newman C, Buesching C, Burke T, Macdonald DW, et al. MHC class II-assortative mate choice in European badgers (Meles meles). Mol Ecol. 2015;24:3138–50.
Aguilar A, Roemer G, Debenham S, Binns M, Garcelon D, Wayne RK. High MHC diversity maintained by balancing selection in an otherwise genetically monomorphic mammal. Proc Natl Acad Sci. 2004;101:3490–4. https://doi.org/10.1073/pnas.0306582101.
Marsden CD, Woodroffe R, Mills MGL, McNutt JW, Creel S, Groom R, et al. Spatial and temporal patterns of neutral and adaptive genetic variation in the endangered African wild dog (Lycaon pictus). Mol Ecol. 2012;21:1379–93. https://doi.org/10.1111/j.1365-294X.2012.05477.x.
Palomares F, Godoy JA, Lopez-Bao JV, Rodriguez A, Roques S, Casas-Marce M, et al. Possible extinction vortex for a population of Iberian lynx on the verge of extirpation. Conserv Biol. 2012;26:689–97.
Ishibashi Y, Oi T, Arimoto I, Fujii T, Mamiya K, Nishi N, et al. Loss of allelic diversity in the MHC class II DQB gene in western populations of the Japanese black bear Ursus thibetanus japonicus. Conserv Genet. 2017;18:247–60. https://doi.org/10.1007/s10592-016-0897-3.
Chapron G, Kaczensky P, Linnell JDC, von Arx M, Huber D, Andren H, et al. Recovery of large carnivores in Europe’s modern human-dominated landscapes. Science. 2014;346:1517–9.
Álvares FB, Blanco JC, Correia J, Cortés Y, Costa G, Llaneza L, et al. Wolf status and conservation in the Iberian Peninsula. In: Frontiers of wolf recovery. Colorado Springs: International Wolf Center; 2005. p. 66–7.
Núñez-Quirós P, Garcia-Lavandera R, Llaneza L. Analysis of historical wolf (Canis lupus) distributions in Galicia: 1850, 1960 and 2003. Ecologia. 2007;21:195–205.
Petrucci-Fonseca F. O lobo (Canis lupus signatus Cabrera, 1907) em Portugal. Problemática da sua conservação: Universidade de Lisboa; 1990.
Valverde JA. El lobo español. Montes. 1971;159:229–41.
Blanco JC, Cortés Y. Dispersal patterns, social structure and mortality of wolves living in agricultural habitats in Spain. J Zool. 2007;273:114–24. https://doi.org/10.1111/j.1469-7998.2007.00305.x.
Sommer R, Benecke N. Late-Pleistocene and early Holocene history of the canid fauna of Europe (Canidae). Mamm Biol. 2005;70:227–41. https://doi.org/10.1016/J.MAMBIO.2004.12.001.
Bencatel J, Ferreira CC, Barbosa AM, Rosalino LM, Álvares F. Research trends and geographical distribution of mammalian carnivores in Portugal (SW Europe). PLoS One. 2018;13:e0207866. https://doi.org/10.1371/journal.pone.0207866.
Hindrikson M, Remm J, Pilot M, Godinho R, Stronen AV, Baltrunaite L, et al. Wolf population genetics in Europe: a systematic review, meta-analysis and suggestions for conservation and management. Biol Rev Camb Philos Soc. 2017;92:1601–29.
López-Bao JV, Blanco JC, Rodríguez A, Godinho R, Sazatornil V, Alvares F, et al. Toothless wildlife protection laws. Biodivers Conserv. 2015;24:2105–8.
Silva P, López-Bao JV, Llaneza L, Álvares F, Lopes S. Blanco JC, et al. Sci Rep: Cryptic population structure reveals low dispersal in Iberian wolves; 2018.
Arbanasic H, Huber D, Kusak J, Gomercic T, Hrenovic J, Galov A. Extensive polymorphism and evidence of selection pressure on major histocompatibility complex DLA-DRB1, DQA1 and DQB1 class II genes in Croatian grey wolves. Tissue Antigens. 2013;81:19–27.
Galaverni M, Caniglia R, Fabbri E, Lapalombella S, Randi E. MHC variability in an isolated wolf population in Italy. J Hered. 2013;104:601–12.
Niskanen AK, Kennedy LJ, Ruokonen M, Kojola I, Lohi H, Isomursu M, et al. Balancing selection and heterozygote advantage in major histocompatibility complex loci of the bottlenecked Finnish wolf population. Mol Ecol. 2014;23:875–89.
Seddon JM, Ellegren H. MHC class II genes in European wolves: a comparison with dogs. Immunogenetics. 2002;54:490–500.
Seddon JM, Ellegren H. A temporal analysis shows major histocompatibility complex loci in the Scandinavian wolf population are consistent with neutral evolution. Proc R Soc B Biol Sci. 2004;271:2283–91. https://doi.org/10.1098/rspb.2004.2869.
Strand TM, Segelbacher G, Quintela M, Xiao L, Axelsson T, Höglund J. Can balancing selection on MHC loci counteract genetic drift in small fragmented populations of black grouse? Ecol Evol. 2012;2:341–53.
Alcaide M. On the relative roles of selection and genetic drift in shaping MHC variation. Mol Ecol. 2010;19:3842–4.
Loiseau C, Richard M, Garnier S, Chastel O, Julliard R, Zoorob R, et al. Diversifying selection on MHC class I in the house sparrow (Passer domesticus). Mol Ecol. 2009;18:1331–40.
Brown JH, Jardetzky T, Saper MA, Samraoui B, Bjorkman PJ, Wiley DC. A hypothetical model of the foreign antigen binding site of class II histocompatibility molecules. Nature. 1988;332:845–50.
Pilot M, Branicki W, Jedrzejewski W, Goszczyński J, Jedrzejewska B, Dykyy I, et al. Phylogeographic history of grey wolves in Europe. BMC Evol Biol. 2010;10:104. https://doi.org/10.1186/1471-2148-10-104.
Randi E. Genetics and conservation of wolves Canis lupus in Europe. Mamm Rev. 2011;41:99–111. https://doi.org/10.1111/j.1365-2907.2010.00176.x.
Montana L, Caniglia R, Galaverni M, Fabbri E, Ahmed A, Bolfíková BČ, et al. Combining phylogenetic and demographic inferences to assess the origin of the genetic diversity in an isolated wolf population. PLoS One. 2017;12:1–19. https://doi.org/10.1371/journal.pone.0176560.
Kennedy LJ, Barnes A, Short A, Brown JJ, Lester S, Seddon J, et al. Canine DLA diversity: 1. New alleles and haplotypes. Tissue Antigens. 2007;69(SUPPL. 1):272–88.
Kennedy LJ, Barnes A, Happ GM, Quinnell RJ, Bennett D, Angles JM, et al. Extensive interbreed, but minimal intrabreed, variation of DLA class II alleles and haplotypes in dogs. Tissue Antigens. 2002;59:194–204.
Kennedy LJ, Angles JM, Barnes A, Carmichael LE, Radford AD, Ollier WER, et al. DLA-DRB1, DQA1, and DQB1 alleles and haplotypes in north American gray wolves. J Hered. 2007;98:491–9.
Godinho R, Llaneza L, Blanco JC, Lopes S, Álvares F, García EJ, et al. Genetic evidence for multiple events of hybridization between wolves and domestic dogs in the Iberian Peninsula. Mol Ecol. 2011;20:5154–66.
Fernández JM, Ruiz de Azua N. Historical dynamics of a declining wolf population: persecution vs. prey reduction. Eur J Wildl Res. 2009;56:169–79.
Rico M, Torrente JP. Caza y rarifacación del lobo en España: investigación histórica y conclusiones biológicas. Galemys. 2000;12:163–79.
Aspi J, Roininen E, Ruokonen M, Kojola I, Vilà C. Genetic diversity, population structure, effective population size and demographic history of the Finnish wolf population. Mol Ecol. 2006;15:1561–76. https://doi.org/10.1111/j.1365-294X.2006.02877.x.
Eizaguirre C, Lenz TL, Kalbe M, Milinski M. Rapid and adaptive evolution of MHC genes under parasite selection in experimental vertebrate populations. Nat Commun. 2012;3:621–6. https://doi.org/10.1038/ncomms1632.
Ciborowski K, Jordan WC, Garcia De Leaniz C, Consuegra S. Temporal and spatial instability in neutral and adaptive (MHC) genetic variation in marginal salmon populations. Sci Rep. 2017;7 February:1–11. https://doi.org/10.1038/srep42416.
Eckert CG, Samis KE, Lougheed SC. Genetic variation across species’ geographical ranges: the central-marginal hypothesis and beyond. Mol Ecol. 2008;17:1170–88.
Volis S, Ormanbekova D, Shulgina I. Role of selection and gene flow in population differentiation at the edge vs. interior of the species range differing in climatic conditions. Mol Ecol. 2016;25:1449–64.
O’Brien SJ, Evermann JF. Interactive influence of infectious disease and genetic diversity in natural populations. Trends Ecol Evol. 1988;3:254–9.
Santos N, Almendra C, Tavares L. Serologic survey for canine distemper virus and canine parvovirus in free-ranging wild carnivores from Portugal. J Wildl Dis. 2009;45:221–6.
Figueiredo A, Oliveira L, Madeira de Carvalho L, Fonseca C, Torres RT. Parasite species of the endangered Iberian wolf (Canis lupus signatus) and a sympatric widespread carnivore. Int J Parasitol Parasites Wildl. 2016;5:164–7. https://doi.org/10.1016/j.ijppaw.2016.04.002.
Conceição-Neto N, Godinho R, Álvares F, Yinda CK, Deboutte W, Zeller M, et al. Viral gut metagenomics of sympatric wild and domestic canids, and monitoring of viruses: insights from an endangered wolf population. Ecol Evol. 2017;7:4135–46.
Müller A, Silva E, Santos N, Thompson G. Domestic dog origin of canine distemper virus in free-ranging wolves in Portugal as revealed by hemagglutinin gene characterization. J Wildl Dis. 2011;47:725–9.
Zhai W, Nielsen R, Slatkin M. An investigation of the statistical power of neutrality tests based on comparative and population genetic data. Mol Biol Evol. 2009;26:273–83. https://doi.org/10.1093/molbev/msn231.
Hall TA. BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucleic Acids Symposium Series. 1999;41:95–98. doi:citeulike-article-id:691774.
Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25:1451–2.
Stephens M, Donnelly P. A comparison of Bayesian methods for haplotype reconstruction from population genotype data. Am J Hum Genet. 2003;73:1162–9.
Bos DH, Gopurenko D, Williams RN, DeWoody JA. Inferring population history and demography using microsatellites, mitochondrial DNA, and major histocompatibility complex (MHC) genes. Evolution. 2008;62:1458–68. https://doi.org/10.1111/j.1558-5646.2008.00364.x.
Garrick RC, Sunnucks P, Dyer RJ. Nuclear gene phylogeography using PHASE: dealing with unresolved genotypes, lost alleles, and systematic bias in parameter estimation. BMC Evol Biol. 2010;10:118. https://doi.org/10.1186/1471-2148-10-118.
Peakall R, Smouse PE. GenAlEx 6.5: genetic analysis in excel. Population genetic software for teaching and research--an update. Bioinformatics. 2012;28:2537–9.
Rousset F. genepop’007: a complete re-implementation of the genepop software for windows and Linux. Mol Ecol Resour. 2008;8:103–6. https://doi.org/10.1111/j.1471-8286.2007.01931.x.
Goudet J. FSTAT (version 1.2): a computer program to calculate F-statistics. J Hered. 1995;86:485–6. https://doi.org/10.1093/oxfordjournals.jhered.a111627.
Kumar S, Stecher G, Tamura K. MEGA7: molecular evolutionary genetics analysis version 7.0. Mol Biol Evol. 2016;33:1870–4.
Tajima F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989;123:585–95.
Fu YX, Li WH. Statistical tests of neutrality of mutations. Genetics. 1993;133:693–709.
Wilson DJ, McVean G. Estimating diversifying selection and functional constraint in the presence of recombination. Genetics. 2006;172:1411–25. https://doi.org/10.1534/genetics.105.044917.
Wilson DJ, McVean G. Estimating diversifying selection and functional constraint in the presence of recombination. Genetics. 2006;172:1411–25.
Awadi A, Ben Slimen H, Smith S, Knauer F, Makni M, Suchentrunk F. Positive selection and climatic effects on MHC class II gene diversity in hares (Lepus capensis) from a steep ecological gradient. Sci Rep. 2018;8:1–14.
Demirbaş Y, Albayrak İ, Özkan Koca A, Stefanović M, Knauer F, Suchentrunk F. Spatial genetics of brown hares (Lepus europaeus Pallas, 1778) from Turkey: different gene pool architecture on either side of the Bosphorus? Mamm Biol. 2019;94:77–85.
Bates D, Mächler M, Bolker B, Walker S. Fitting linear mixed-effects models using lme4. J Stat Softw. 2014;67:1–48. https://doi.org/10.18637/jss.v067.i01.
Kuznetsova A, Brockhoff PB, Christensen RHB. lmerTest Package: Tests in Linear Mixed Effects Models. J Stat Softw. 2017;82:1–26. https://doi.org/10.18637/jss.v082.i13.
Dray S, Dufour A-B. The ade4 package: implementing the duality diagram for ecologists. J Stat Softw. 2007;22:1–20. https://doi.org/10.18637/jss.v022.i04.
Jombart T. Adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics. 2008;24:1403–5.
Slatkin M. A measure of population subdivision based on microsatellite allele frequencies. Genetics. 1995;139:457–62.
Jost L. GST and its relatives do not measure differentiation. Mol Ecol. 2008;17:4015–26. https://doi.org/10.1111/j.1365-294X.2008.03887.x.
Excoffier L, Lischer HEL. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and windows. Mol Ecol Resour. 2010;10:564–7.
Gerlach G, Jueterbock A, Kraemer P, Deppermenn J, Harmand P. Calculations of population differentiation based on GST and D: forget GST but not all of statistics! Mol Ecol. 2010;19:3845–52. https://doi.org/10.1111/j.1365-294X.2010.04784.x.
We acknowledge the Portuguese Institute for Nature Conservation and Forestry (SMLM/ICNF) for providing wolf tissue samples in Portugal. For collection of wolf samples in Spain, we thank Consejería de Medio Ambiente del Pricipado de Asturias, Consellería de Medio Ambiente de la Xunta de Galicia, Consejería de Medio Ambiente de la Junta de Castilla y León, Gobierno de Cantabria and Parque Nacional Picos de Europa. We are grateful to the Subject Editor, Dr. Helmut Schaschl, and to Dr. Franz Suchentrunk and another anonymous reviewer for valuable comments in a previous version of this manuscript. We thank Pedro Monterroso for helping implementing the linear mixed effects models. This is scientific paper number 23 of the Iberian Wolf Research Team (IWRT).
RGR was supported by a post-doctoral fellowship (project IF/00564/2012/CP0159/CT0003, funded by the Portuguese Foundation for Science and Technology, FCT), RG by the research contract IF/00564/2012 from FCT, and JVLB was supported by a Ramon & Cajal research contract (RYC-2015-18932) from the Spanish Ministry of Economy, Industry and Competitiveness. This work was partially supported by FCT (project PTDC/BIA-EVF/2460/2014 to RG) and by Norte Portugal Regional Operational Program (NORTE2020), under the PORTUGAL 2020 Partnership Agreement, through the European Regional Development Fund (ERDF) (NORTE-01-0145-FEDER-000007). These funding sources had no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.
Availability of data and materials
Microsatellite and MHC genotype data were uploaded to the Open Science Framework repository under the identifier https://doi.org/10.17605/OSF.IO/AKT5H.
Ethics approval and consent to participate
The Portuguese Institute for Nature Conservation and Forestry (SMLM/ICNF), Consejería de Medio Ambiente del Pricipado de Asturias, Consellería de Medio Ambiente de la Xunta de Galicia, and Consejería de Medio Ambiente de la Junta de Castilla y León, Gobierno de Cantabria were responsible for the collection of samples and for permits to use them in this study.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Table S1. DLA-DRB1, DLA-DQA1 and DLA-DQB1 alleles and their frequency (f) in the three assigned demographic groups and whole Iberian wolf range. (PDF 19 kb)
Table S2. Number of three-locus (DLA-DRB1/DQA1/DQB1) haplotypes (N) and their frequency (f) in the three assigned demographic groups and whole Iberian wolf range. (PDF 14 kb)
Table S3. Comparison of the DLA-DRB1, DQA1, DQB1 alleles and three-locus haplotypes (DRB1/DQA1/DQB1) in the Iberian (this study), Croatian , Italian , Finnish  wolf populations. Underlined: shared alleles and haplotypes between Iberia and at least one of the three other European wolf populations; *alleles and haplotypes found in other European wolf populations [37, 38] ¥alleles found in the North American wolf population [37, 48]; # haplotypes that are possibly found in other European wolf populations, three-locus haplotypes were not reconstruct but all alleles are present . NI: no DRB1 allele detected. (PDF 162 kb)
Table S4. Genetic diversity values for the three MHC loci in the Iberian wolf range. Sample size (n), mean number of alleles (Na), mean number of effective alleles (Ne), mean observed and expected heterozygosities (Ho and He), mean inbreeding coefficient (FIS), number of segregating sites (S), nucleotide diversity (π), number of mutations (η) Watterson’s mutation parameter (θW), neutrality tests of Tajima’s D and Fu &Li D*, and estimates of average nucleotide and aminoacid distances. Statistical significance: *P < 0.05, **P < 0.02. Standard error estimates shown in brackets. (PDF 139 kb)
Table S5. Positively selected codons in each Iberian wolf demographic group according to OmegaMap analysis and codons putatively belonging to the peptide-binding regions (PBR) suggested by  for DLA-DRB1, DLA-DQA1 and DLA-DQB1. Bold and underlined positions denote codons identified to be under positive selection with posterior probability > 0.95. (PDF 110 kb)
Table S6. Model selection information, including Akaike Information Criterion (AICc), delta AICc (ΔAICc) and weight of selected model for number of alleles (Na), allelic richness (AR), observed (Ho) and expected (He) heterozygosities. (PDF 128 kb)
Table S7. Estimates of combined diversity measures of each demographic group, standard errors (S.E.) and probability inferred from t-value, according to the selected model. (PDF 11 kb)
Table S8. Mean number of alleles (Na), mean number of effective alleles (Ne), mean allelic richness (AR), observed (Ho) and expected (He) heterozygosities, and mean inbreeding coefficient (FIS) at the microsatellite and MHC loci, in three demographic groups (persistent, expanding and isolated), including standard deviation (s.d). The probability of the values to differ between microsatellite and MHC loci was tested by a one-tailed t-test (p-value). Values with statistical significance are in bold (p < 0.05). (PDF 28 kb)
Table S9. Thermocycling profile for PCR amplification of the three MHC fragments used in this work. Annealing temperature (AT) was set initially at 64 °C (DRB1), 54 °C (DQA1) and 76 °C (DQB1). (PDF 97 kb)