Neutral and Selective Processes Shape MHC Diversity in Roe Deer in Slovenia

Simple Summary Disease prevention and appropriate wildlife management are among the major challenges in wildlife conservation. In the present study, we made a first assessment of the variability of major histocompatibility complex (MHC) genes in roe deer in Slovenia and evaluated local population adaptation by comparing MHC variability with neutral microsatellites. We discovered three new MHC DRB exon 2 alleles in addition to seven previously described in the literature. Moreover, we found evidence of historical positive selection, as selection analysis indicated that approx. 10% of the encoded amino acids were subjected to episodic positive selection. This study provides the basis for further research on immunogenetic variation in roe deer and highlights opportunities to incorporate genetic data into science-based population management. Abstract Disease control and containment in free-ranging populations is one of the greatest challenges in wildlife management. Despite the importance of major histocompatibility complex (MHC) genes for immune response, an assessment of the diversity and occurrence of these genes is still rare in European roe deer, the most abundant and widespread large mammal in Europe. Therefore, we examined immunogenetic variation in roe deer in Slovenia to identify species adaptation by comparing the genetic diversity of the MHC genes with the data on neutral microsatellites. We found ten MHC DRB alleles, three of which are novel. Evidence for historical positive selection on the MHC was found using the maximum likelihood codon method. Patterns of MHC allelic distribution were not congruent with neutral population genetic findings. The lack of population genetic differentiation in MHC genes compared to existing structure in neutral markers suggests that MHC polymorphism was influenced primarily by balancing selection and, to a lesser extent, by neutral processes such as genetic drift, with no clear evidence of local adaptation. Selection analyses indicated that approx. 10% of amino acids encoded under episodic positive selection. This study represents one of the first steps towards establishing an immunogenetic map of roe deer populations across Europe, aiming to better support science-based management of this important game species.


Introduction
Genetic diversity is one of the most important mechanisms that enable a population's response and adaptation to environmental changes [1,2]. Assessment of genetic diversity (particularly of genes influencing immune response) is therefore urgently needed when developing effective conservation and management programs for wildlife species [3,4]. Since neutral genetic variation provides an incomplete picture of the evolutionary potential of populations (e.g., [5,6]), it is also important to monitor adaptive genetic diversity, which is defined as "genetic variation that produces an advantage in fitness" [7]. Polymorphism in the major histocompatibility complex (MHC), a family of highly variable genes, plays a key

Ethical Statement
All animals used in the study were legally harvested during regular hunting activities prescribed by the state of Slovenia in annual wildlife management plans. No animal was shot or otherwise killed for the purposes of this study.

Study Area and Sampling
Samples of 156 roe deer female yearlings were collected throughout Slovenia, i.e., in 58 hunting grounds from all 15 Slovene hunting management districts. Therefore, a full range of environmental factors and population traits of roe deer presence in the country was covered (Table S1).
All samples (either muscle tissue or uteri) were collected by hunters or wildlife researchers immediately after the harvest event within the hunting season (1 September-31 December) in the period 2013-2015, stored in alcohol solution, and frozen until transportation to the genetic laboratory at the University of Primorska.

DNA Extraction and MHC Genotyping by Sequencing
We extracted DNA from tissue samples using a peqGOLD Tissue DNA Kit (S-Line) and followed the manufacturer's instructions (VWR International, Leuven, Belgium). DNA extraction and quality control procedures are described in detail in Buzan et al. [44].
We amplified a 249-bp fragment of the exon 2 of the roe deer MHC DRB gene using primers LA31 (ATCCTCTCTCTGCAGCACATTTCC) and LA32 (TTCGCGTCACCTCGC-CGCTG), which were modified based on primers originally designed for cattle [45]. For the Ion Torrent PGM system, we designed extended primers as follows: forward primeradapter sequence (30 bp), 10-12 bp barcodes with a specific "GAT" linker to distinguish individuals and primer LA31; reverse primer-adapter sequence and primer LA32 without adapter sequence.
PCR amplification was performed in triplicates in 25 µL reaction mixtures of 50 ng genomic DNA, 0.5 mM of each primer, 0.5 µL dNTPs, 2.5 mM MgCl 2 , 0.5 µL HotStarTaq (Qiagen AllTaq core kit, Qiagen, Germany), 5.0 µL Q-solution, and 6.2 µL H 2 O. The PCR reaction started with an initial denaturation (2 min, 95 • C) followed by 35 cycles of 10 s denaturation at 95 • C, 30 s annealing at 60 • C, and 30 s extension at 72 • C. The final elongation was performed at 72 • C for 10 min. The amplicons from the triplicates were then pooled and purified with magnetic particles Agencourt ® AmPure ® (Agencourt Bioscience Corporation, A Beckman Coulter Company, Beverly, MA, USA). We measured concentrations of pooled and cleaned amplicons by Qubit 3.0 fluorometry using Qubit dsDNA HS (High Sensitivity) Assay Kit reagents (Invitrogen, Carlsbad, CA, USA). Samples were then normalized to 5 ng and combined into a single library. The library was again purified with Agencourt ® AmPure ® magnetic particles. The size and quality of the amplicons were determined using the Agilent DNA High Sensitivity Kit on a 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA), according to the manufacturer's recommendations. Before loading the library onto the chip, the library was normalized to 100 pM and sequenced using the Ion Torrent S5 on an Ion 530 chip (Thermo Fisher Scientific, Waltham, MA, USA).
For allele calling, we used the pipeline in the Amplicon Sequence Assignment (Am-pliSAS) web tool developed for high-throughput genotyping of duplicated polymorphic genes, such as MHC [46]. Filtering of raw data was performed with AmpliCLEAN by removing reads with a Phred quality score <20 and filtering of all reads <250 bp and >300 bp. AmpliSAS clusters true variants with their potential artefacts based on the platform-specific error rates. We used AmpliSAS's default parameters for Ion Torrent sequencing technology: a substitution error rate of 0.5% and an indel error rate of 1%. An accurate length was required to identify the dominant sequence within a cluster. We did not expect more than two DRB variants per individual, so we kept the "minimum dominant frequency" clustering threshold at 25%, based on previously published work on roe deer [25,34,35,47]. We discarded variants with a frequency <1% within an amplicon. True variants of the DRB exon 2 fragments were aligned and translated into protein sequences to check for evidence of pseudogenes, such as the presence of premature stop codons or indels. A maximum of 5000 reads per amplicon was used for genotyping to reduce computational load, and we repeated the analysis three times.

Sequence Diversity
The unique sequences were aligned, edited, and confirmed to be roe deer MHC DRB exon 2 alleles using MEGA X [48] by comparing them with alleles downloaded from GenBank (Table S2). DnaSP v.6.12 [49] was used to calculate the average number of nucleotide differences (k) and number of segregating (variable) sites (S). The average pairwise nucleotide distances (Kimura 2-parameter model; K2P) and Poisson-corrected amino acid distances were calculated in MEGA for the following: overall, Antigen Binding Sites (ABS), and non-ABS. The locations of the putative ABS and non-ABS were inferred from the human MHC II molecule structure [50].

Clusters/Populations Genetic Parameters and Molecular Variance
Values for nucleotide diversity (π), 4Nµ for autosomal genes of diploid organism theta (θ), Tajima's D based on the site frequency spectrum from DNA sequences, and allele number (A) were calculated using DNA SP. We visualized spatial patterns of genetic diversity across the landscape by interpolating and mapping the allelic richness (AR) for each population using the SPADs and R version 4.0.5 [51]. Observed and expected heterozygosity values (Ho and He) and pairwise F ST between groups were calculated and tested by Arlequin version 3.5 [52].
We tested conformance of the allele frequencies with Hardy-Weinberg expectations, both within each of the three clusters (defined by neutral markers; Figure S1; [44]), with groups/populations defined by geographical feathers, and overall, by using the complete enumeration algorithm of Louis and Dempster [53] as implemented in Genepop 3.5 [54]. In Genepop 3.5, tests addressing specific hypotheses of heterozygote excess and deficit were performed, and differences in allele and genotype frequencies between clusters/groups were assessed with Markov chain Monte Carlo approximations of Fisher exact tests [54].
MHC-based population structure in the study area was assessed by analysis of molecular variance (AMOVA; [55]) implemented in Arlequin with 10,000 permutations and by using program STRUCTURE 2.3.4 [56]. Mantel tests were performed in adegenet in R to examine the isolation-by-distance relationship between estimates of MHC F ST /(1 − F ST ) and the natural logarithm of the geographic distance.
We found no evidence for population genetic differentiation at MHC DRB exon 2 alleles ( Figure S2); for further analysis, we pooled individuals based on the location of their death into (i) ten predefined groups based on the different historical management of roe deer and geographical characteristics of Slovenia (for the purposes of the study, also called populations) and (ii) three genetically differentiated clusters revealed by neutral genetic markers (for details, see Table S1, Figure S1, and Buzan et al. [44]). This partitioning is important for science-based population management within the country as it provides additional insight into the adaptive immunity of roe deer within hunting management units.

Clustering of DRB Exon 2 Alleles
The evolutionary relationships between all DRB exon 2 alleles were analyzed by a median-joining network implemented with the software Network 4.6.0.0 (http://www. fluxus-engineering.com/sharenet.htm; accessed on 20 September 2021). All positions were equally weighted and the ε parameter was set to zero.
We also included published roe deer DRB exon 2 allele sequences from GeneBank, which were not found in Slovenia: Caca-DRB*0101, Caca-DRB*0202, Caca-DRB*0203, Caca-DRB*0204 (Table S2; [25,34]); in the analysis, they were represented as one individual. The MJ network was rooted using haplotypes of red deer (Cervus elaphus; GenBank number EU573258). We attempted to study functional MHC DRB diversity by clustering alleles into supertypes (i.e., groups of alleles assumed to have similar peptide-binding capacities). Therefore, we performed clustering based on amino acid polymorphism at the positively selected amino acid sites (PSSs). For the PSSs, we retained all codons with Bayes Empirical Bayes (BEB) posterior probability >95% in the M8 model ( Figure S3). The PSS of each allele was numerically characterized by a set of five physicochemical descriptors for each amino acid. To cluster the alleles into supertypes, we performed DAPC using the "adegenet" package. First, we defined clusters (k) using the find.clusters() function and kept all principal components (PCs). We compared different clustering solutions based on the Bayesian Information Criterion (BIC) values for increasing number of clusters. According to recommendation, the optimal number of clusters must be chosen as the number of clusters with the lowest BIC value, after which the BIC value decreases by a negligible amount. However, as the BIC values steadily decline, there was no point of stabilization or 'elbow' at which the value began to increase ( Figure S4). The likely reason for such result is that only six sites showed signs of positive selection, which did not provide enough resolution for reliable supertyping. We therefore decided that downstream analyses based on supertypes were not relevant and used a simplified approach that included alleles.

Detecting Signatures of Recombination and Selection in MHC DRB Exon 2 Alleles in Roe Deer
The codon-based Z-test of selection (implemented in MEGA X) was used to calculate the average nonsynonymous/synonymous substitution rate ratio (dN/dS = ω) across the entire sequence using the counting method of Nei and Gojobori [57], corrected for multiple substitutions [58].
Using a one-tailed Z-test with standard errors resulting from 10,000 bootstrap replicates, we calculated the rate of positive selection sites separately over the entire DRB exon 2 sequences and over the extracted antigen binding codons. We also used an approach with a more powerful maximum-likelihood method that allows the dN/dS ratio to vary among codon sites; this was implemented in the EasyCodeML program [59], which identifies codons affected by positive selection based on a Bayesian approach and was used in additional analysis for testing historical selection. The models implemented in this study were M2a and M8. Furthermore, we assessed the influence of positive selection on individual codons using single-likelihood ancestor counting (SLAC), fixed effects likelihood (FEL), mixed effects model of evolution (MEME), and fast unconstrained Bayesian approximation (FUBAR) methods performed in the Datamonkey 2.0 server (http://www.datamonkey. org/; accessed on 20 September 2021; see also [60]).

Determination of Genetic Variability among Groups and Clusters Based on Neutral Microsatellite Loci
We reran the analysis on 156 genotyped yearlings, previously used in Buzan et al. [44]. Population structure was determined using program STRUCTURE 2.3.4 [56,71]. Mean number of alleles (A), allelic richness (AR), observed (Ho) and expected (He) heterozygosity, deviation from Hardy-Weinberg equilibrium, pairwise F ST , AMOVA, and isolation-by-distance were estimated using the same software and protocols described by Buzan et al. [44].

Diversity of MHC DRB Alleles
We found 10 functional alleles for MHC DRB exon 2 coding for different amino acid sequences among tested individuals. No evidence of multiple locus amplification was found, confirming previous reports for roe deer [25,34]. A summary of the allele variants found in each individual is provided in Table S1. The three novel DRB exon 2 sequences identified in this study were named in accordance with previously established nomenclature as Caca-DRB*0402, Caca-DRB*, and Caca-DRB*0404, and were deposited in GenBank with accession numbers from OL355104 to OL355106. One of the new alleles, Caca-DRB*0402, has a deletion in codon 65 (del65), which has already been described in alleles Caca-DRB*0202, Caca-DRB*0301, Caca-DRB*0302, and Caca-DRB*0303 [25,34,35].
An analysis of nucleotide alignment revealed 22 segregating (variable) nucleotide sites distributed across 83 codons. The overall nucleotide evolutionary distance and the amino acid evolutionary distance were 5% and 9%, respectively. The average number of nucleotide differences among alleles was k = 9.833 (Table 1). Notes: k-average number of nucleotide differences; S-the number of segregating sites; ABS-antigen binding sites; non-ABS-nonantigen binding sites. * SD is given in parenthesis, except for S, for which the number of nonsynonymous sites is provided.

Clusters/Populations Genetic Parameters and Molecular Variance
We found no evidence for population genetic differentiation at MHC DRB exon 2 alleles ( Figure S2). Table 2 lists the values of diversity parameters for three genetic clusters of roe deer in Slovenia as recognized by neutral microsatellite diversity data (Table S1, Figure S1, [44]). The numbers of MHC alleles ranged from 8 in the southwestern cluster to 10 in the central cluster. Allelic richness was similar in all clusters (8.000-8.602). The expected heterozygosity was lower in the southwestern (0.414) than in other two clusters (0.524-0.593), but θ was higher in the first one (0.022) than in the central and the northeastern clusters (0.017). Tajima's D values were positive in all clusters and significant for overall data. The average pairwise F ST values between southwestern and central clusters were negative (equal to zero) and nonsignificant (−0.0029), while between southwestern/northeastern and central/northeastern clusters were positive but also nonsignificant (0.0089 and 0.0002, respectively). The average nucleotide diversity was π = 0.045, and was very similar in all three clusters ( Table 2). Notes: n-number of individuals; A-number of alleles; AR-allelic richness; Ho-observed heterozygosity; Heexpected heterozygosity; π-nucleotide diversity; θ-4 Nµ for autosomal genes of diploid organisms; Tajima's D (value in bold is significant; p < 0.05). Deviation from Hardy-Weinberg equilibrium was significant for overall and cluster data.
For the ten groups/populations (Table S3), the numbers of alleles in the populations ranged from 4 in C1 to 10 in C4 (Figures 1 and S5). Allelic richness was the highest in C2 (AR = 5.480) and the lowest in C1 (4.000). Observed heterozygosity in one of the northeastern (N2) population and the southwestern (S1) population were lower (He = 0.471 and  (Table S4), which means there is no genetic subdivision, except between populations N4 and C2, where the F ST value was positive and significantly different. of diploid organisms; Tajima's D (value in bold is significant; p < 0.05). Deviation from Hardy-Weinberg equilibrium was significant for overall and cluster data.
For the ten groups/populations (Table S3), the numbers of alleles in the populations ranged from 4 in C1 to 10 in C4 (Figure 1, Figure S5). Allelic richness was the highest in C2 (AR = 5.480) and the lowest in C1 (4.000). Observed heterozygosity in one of the northeastern (N2) population and the southwestern (S1) population were lower (He = 0.471 and 0.414, respectively) compared with other populations (0.672-0.849). Nucleotide diversity was the highest in C3 (π = 0.050) and the lowest in N3 population (0.039). Parameter θ varied from 0.022 to 0.030: it was lowest in C4 and C3 and highest in C1 and C2 populations. Tajima's D values were positive in most populations except C1 and N3, which indicates the possibility of a recent population increase after a bottleneck. The average pairwise FST values between most of the populations were negative and nonsignificant (Table  S4), which means there is no genetic subdivision, except between populations N4 and C2, where the FST value was positive and significantly different.  Hierarchical AMOVA revealed that most of the genetic variation in the population structure can be explained by within-population differences (98.9%), while only 0.7% of MHC variability is due to differences among populations within three clusters. However, this estimate was not statistically significant (Table S5). We did not detect significant isolation-by-distance among populations (Mantel test: r = 0.211, p = 0.406).

Recombination and Natural Selection in MHC DRB Alleles in Roe Deer
In recognized alleles, the overall number of codons was 83, of which 67 (81%) were non-ABS amino acids and 16 (19%) were ABS. Negative ω (dN/dS) for non-ABS indicates possible selective removal of deleterious alleles, which ultimately reduces variation in a population. The global estimates of ω, averaged across all codon sites using the codonbased Z-test of selection, showed the presence of positive selection at the DRB locus. The nonsynonymous mutation rate (dN = 0.05) exceeded the synonymous one (dS = 0.02) ( Table 3) identified up to eight codons predicted to be affected by positive selection. The selection models revealed different levels of selection pressure at the loci analyzed. The DRB locus showed signs of strong selection pressure with six positively (posterior probability >99%) selected codon sites (Table 4). Mean values of ω for individual codons at each locus are presented in Figure S3. Recombination analysis revealed two recombination events that resulted in two alleles originating from their combination process (Table 5). Table 3. Relative rates of nonsynonymous (dN) and synonymous (dS) substitutions (with standard errors) calculated in roe deer DRB exon 2 alleles (found in this study and previously published in [25,34]) for antigen binding sites (ABS) and non-ABS. Statistical significance (p-value) was tested using the one-tailed Z-test with standard errors resulting from 10,000 bootstrap replicates.

Main Parameters
Overall  Note: The codon sites (del65 is included in the numbering) inferred to be under selection with posterior probability >99% are listed in bold, and sites with posterior probability >95% are in standard font. ω values that represent synonymous vs. nonsynonymous substitutions (dN/dS) for particular codons are given in Figure S3.

Evolutionary Relationships among DRB Exon 2 Alleles in Roe Deer
The median-joining (MJ) network among 14 alleles (eleven previously described [25,34,35] and three novel ones) showed that the number of mutational steps between adjacent alleles ranged from one (several cases) to eight (between Caca-DRB*0401 and Caca-DRB*0403) ( Figure 2). The maximum number of mutational steps across the network between Caca-DRB*0101 and Caca-DRB*0404 was 42, with one hypothetical allele between Caca-DRB*0202 and Caca-DRB*0203/Caca-DRB*0304. The MJ network grouped seven alleles with a deletion of codon 65 close together in one group and the other five without deletion of this codon in a separate group.

Genetic Diversity and Genetic Structure Based on Neutral Loci
Reanalysis of genotyped data on 156 yearling females confirmed previously published results on 241 roe deer individuals, including adults [44]. The best population assignment model resulting from the structure analysis (using K after [56]) divided individuals into three groups (K = 3), which matched fairly well with their geographical origin ( Figure S1, Table S1) as well as with the population structure previously found using a more heterogenous sample set [44].

Evolutionary Relationships among DRB Exon 2 Alleles in Roe Deer
The median-joining (MJ) network among 14 alleles (eleven previously described [25,34,35] and three novel ones) showed that the number of mutational steps between adjacent alleles ranged from one (several cases) to eight (between Caca-DRB*0401 and Caca-DRB*0403) (Figure 2). The maximum number of mutational steps across the network between Caca-DRB*0101 and Caca-DRB*0404 was 42, with one hypothetical allele between Caca-DRB*0202 and Caca-DRB*0203/Caca-DRB*0304. The MJ network grouped seven alleles with a deletion of codon 65 close together in one group and the other five without deletion of this codon in a separate group. Figure 2. Evolutionary relationships among DRB exon 2 alleles in European roe deer represented by a median-joining network. Alleles are indicated by circles whose size is proportional to the number of roe deer individuals. Number of mutations separating nodes is represented by slashes that cross network branches. A small black circle indicates a hypothetical allele. Alleles marked in blue have been found only in Slovenia (our study), red-marked alleles are known in roe deer from other countries/studies and were also found in Slovenia, and green-marked alleles were not found in our study (represented in the network only by one individual per allele). Nonsynonymous changes are red-marked and indicated by amino acid change (protein variants) coded from A to Q [50]. Amino acid position 86 (nucleotides position 256/257) was proved to be under positive selection by four of our tests.

Genetic Diversity and Genetic Structure Based on Neutral Loci
Reanalysis of genotyped data on 156 yearling females confirmed previously published results on 241 roe deer individuals, including adults [44]. The best population as- Figure 2. Evolutionary relationships among DRB exon 2 alleles in European roe deer represented by a median-joining network. Alleles are indicated by circles whose size is proportional to the number of roe deer individuals. Number of mutations separating nodes is represented by slashes that cross network branches. A small black circle indicates a hypothetical allele. Alleles marked in blue have been found only in Slovenia (our study), red-marked alleles are known in roe deer from other countries/studies and were also found in Slovenia, and green-marked alleles were not found in our study (represented in the network only by one individual per allele). Nonsynonymous changes are red-marked and indicated by amino acid change (protein variants) coded from A to Q [50]. Amino acid position 86 (nucleotides position 256/257) was proved to be under positive selection by four of our tests.
The average number of alleles per population ranged from 3.3 to 6.0. Due to smaller sample size, allelic richness across populations (3.58-3.88) was lowest as found previously [44], with the highest value in Slovenske gorice/Prekmurje (N4) and the lowest in the coastal region (S1). A similar pattern was also found for Ho (0.58-0.67) and He (0.60-0.66), with Dinaric Mountains population (C3) having the lowest Ho value (Table S6). The three genetic clusters defined by STRUCTURE had very similar genetic diversity parameters (Ho varied between 0.61 and 0.66 and AR between 6.02 and 6.12, respectively; Table S7). Pairwise F ST values between groups/populations ranged from 0.004 to 0.069, and more than half of them were significantly different from zero (52%; Table S8). The highest F ST value was observed between population S1 (coastal Slovenia) and one of the central (C2) and northern (N4) populations, i.e., similar to those found by Buzan et al. [44]. Between genetic clusters, the F ST values were not significant (Table S9).
On the other hand, AMOVA did not strongly support the three-group structuring revealed by the STRUCTURE, as the among-group variance was low and not significant. Nevertheless, it supported a significant association of individual's genotype with its geographical position, as differences among populations within groups were significant (Table S10; [44]).
Microsatellite genetic distances between individuals were positively correlated with the geographical distances between them (Mantel test: t = 10.03, p < 0.001, R 2 = 0.008). Although this finding confirms to some extent that the geographical distance among individuals has some effect on the spatial genetic structure of roe deer in Slovenia, its influence is rather weak, which was already shown by Buzan et al. [44].

Discussion
Our study confirmed that roe deer in the study area (i.e., the contact zone between the Alps and the Dinaric Mountains in central Europe, where data on the MHC gene polymorphism of this species have been completely lacking) have maintained a high level of MHC diversity compared with those in northern Europe [34] despite reductions in abundance and effective population size in the past [73]. This may be due to strong historical demographic fluctuations consistently attributed to founder events [74] and has also been previously described in roe deer from France [25,75]. This coincides with a high variation in biochemical markers in roe deer from central Europe [76], which Mikko et al. [34] already hypothesized should lead to more MHC alleles being found there.
Roe deer in Slovenia as well as other roe deer populations studied [25,34,75] have lower MHC diversity compared with other wild ungulates (Alpine chamois (Rupicapra rupicapra) [77]; Pyrenean chamois (Rupicapra pyrenaica) [78]; red deer [26,79]; wild boar (Sus scrofa) [80]). This is probably an influence of Pleistocene glaciation [34,81], when gene pool depletion was caused by inbreeding combined with genetic drift [82]. However, the history of the species in the study area probably also plays an important role, since, during the 20th century, roe deer from the central part of Slovenia increased its range both in the Sub-Mediterranean and Karst regions in the south, as well in the open agricultural areas of the Sub-Pannonian region in the northeast [44]. Therefore, gene flow has been primarily limited within our study area (all of Slovenia), which is also supported by the pronouncedly philopatric spatial behavior of the species (e.g., [37,38]).
The most common allele found was Caca-DRB*0301 with a deletion of codon 65, which was rare in northern Europe [34]. Interestingly, four other alleles with a deletion of codon 65 were found in Slovenian roe deer: three have already been found in France [25], but Caca-DRB*0404 is a new one. Mikko et al. [34] suggested that these alleles evolved from a common ancestor of the Bovidae and Cervidae (approx. 20 million years ago). The presence of del65 alleles in these two clearly distinct families suggests the possibility of a crossspecies mode of inheritance of a del65 allelic lineage that arose before the split of bovids and cervids [35]. However, it is extremely difficult to determine with certainty whether shared sequence motifs between unrelated MHC alleles are due to common ancestry or convergent evolution [83]. Nevertheless, the del65 mutation is interesting from a functional point of view because it could affect antigen binding of exon 2 of the DRB gene [35,84]; thus, the mutation appears to be functionally acceptable or even favored.
The evolutionary distances between nucleotides and amino acids were much higher in ABS than in non-ABS (17% vs. 2%; 43% vs. 3%), indicating a direct effect on binding properties for antigenic peptides [85]. Mutations in the ABS domain are usually associated with differential disease resistance [86]. The observed DRB exon 2 alleles in roe deer showed similar nucleotide and amino acid divergence at the ABS and non-ABS positions as previously described by Mikko et al. [34] and Quéméré et al. [25].

Genetic Diversity and Population Structure
The lack of spatial genetic differentiation in the MHC DRB is interpreted as a consequence of balancing selection where genetic structuring at the MHC loci is expected to be low because MHC polymorphism is maintained over the long-term across populations, even in the case of restricted gene flow [23,87]. In our study, this is supported by STRUC-TURE and AMOVA analyses ( Figure S2, Table S5), a positive and significant Tajima's D, and lower nucleotide diversity compared with haplotype diversity (Table S3).
Probably, MHC genetic diversity was affected by balancing selection and to a lesser extent by neutral processes, but we did not found evidence of local adaptation, which has already been suggested for roe deer by Quéméré et al. [25]. Our conclusion is supported by several facts that showed weaker patterns of genetic structuring for MHC compared with microsatellite loci (for the latter, see also [44]): (i) Cluster analyses revealed no structure in the MHC ( Figure S2, Table S5), whereas three evident genetic clusters were found using microsatellites ( Figure S1). (ii) Only one of the pairwise F ST comparisons was significant for MHC (between C2 and N4 populations), while F ST values were much higher for microsatellites and were significant in 22 cases, i.e., 52% (Table S4, Table S8). (iii) We found no evidence of isolation by geographic distance for MHC (p = 0.41), while a strong and significant pattern was confirmed for neutral microsatellite loci (p < 0.001; see also [44]).
When roe deer individuals were grouped into three clusters defined by microsatellite data (Table S3), the number of MHC DRB alleles and observed heterozygosity did not differ significantly among those clusters ( Table 2). In contrast, the number of MHC DRB haplotypes varied significantly among ten groups/populations predefined on the basis of differences in historical management and geographic characteristics of Slovenia ( Figure S5). In all studied populations, we found moderate MHC diversity (He = 0.286-0.833) relative to the observed heterozygosity in other European roe deer populations, which ranged from low in northern Europe [34] to relatively high in France [25]. We also found the same number of alleles (10) and polymorphic sites (22), but two more nonsynonymous sites, as in the three populations from France [25], whereas only four alleles were found in northern European populations, i.e., in Norway and Sweden [34].
Heterozygosity and allelic richness were maintained in all populations (Table 2), despite demographic changes, which were proven with the neutral genetic diversity data (Table S6 and [44]). Isolation and significant effects of genetic drift may contribute to lower allelic richness in the Dinaric Mountains. This may be due to environmental conditions unsuitable for roe deer in this region, i.e., the predominance of dense old forests with a closed canopy (mainly Abieti-Fagetum), which cover tens of thousands of hectares along this mountain chain [72]. The low allelic richness in one of the southern populations (S1) and the sub-Pannonian population (N2) (Figure 1) may also be related to the history of the species in the 19th century, when roe deer in southern Slovenia experienced a genetic bottleneck, while the sub-Pannonian populations had a smaller effective size and were less likely to require individuals to disperse over long distances and/or roam over agricultural land, which is also indicated by neutral genetic diversity data.
As most wildlife populations have undergone complex historical demographic changes and range shifts in recent decades in response to human-mediated land-use changes, it is necessary to clarify their consequences on adaptive genetic diversity [88]. Due to past bottlenecks and recent founder effects, large shifts in allele frequencies can occur on the colonization front without selection. The low MHC diversity observed in the most southern population could also be due to declining parasite-mediated selection [89]. Lower parasite pressure during the colonization process in the dry and warm sub-Mediterranean climate may have caused roe deer to shift energetic resources from host immunity to reproduction [90]; however, we believe this was not the case in Slovenia, as the reproductive potential of the species in the coastal area is lower than in other parts of the country, particularly in the northeast [44,91]. Thus, genetic drift and migrations could be the dominant contemporary forces shaping MHC variation in roe deer (Table 4). Positive selection most likely affects only a few sites over time [92,93]. We detected sites under positive selection and an excess of nonsynonymous substitutions in amino acid residues, suggesting that balancing selection played a role in maintaining this MHC sequence polymorphism. This result is consistent with most of the empirical and theoretical work on MHC in small populations or those experiencing bottlenecks [94,95].

Selection and Recombination
Both maximum-likelihood codon-based selection models support historical selection at the MHC DRB exon 2. Importantly, all codons detected under positive selection were at ABS sites [96], which are involved in pathogen binding recognition [97]. Recombination was detected at one site near an ABS codon. We found a signature of positive selection by the overall ratio ω = 1.40, which is similar to that found in other wild ruminants in Europe (Alpine and Pyrenean chamois [78]; red deer [26,79]). Although our results indicate that positive selection likely played an important role (the amino acid in position 86 was proved to be under positive selection by four of our tests; Figure S3), this does not mean that selection is still acting on current populations. The excess of nonsynonymous mutations may take many generations/cohorts to disappear after the selection process stopped [98]. Short-lived balanced polymorphisms favor rapid adaptive response of populations to changing selection pressures through small adjustments in the level of gene expression [99] and occur in high frequency in the genome, especially in alternative combinations of alleles involved in pathogen resistance [100].
The multigene/multipathogen association approach developed by Quéméré et al. [75] suggested that, in roe deer, innate immunity plays a key role in pathogen-mediated directional selection rather than heterozygote advantage in MHC genes. However, they also stressed the limitation of their work because of insufficient statistical power to test the role of rare variants and the lack of a long time series [101]. In line with this, our study also revealed a strong need for continuing research into MHC genes, their diversity, and consequences for the resilience and fitness of individuals. Indeed, by employing the most important game species (e.g., roe deer or other wild ungulates) as models it may be possible to overcome the issue of a large sample set covering longer time series.

Conclusions
Our results suggest that MHC DRB exon 2 diversity in roe deer was influenced by balancing selection and to a lesser extent by neutral processes; however, we did not find evidence of local adaptation. Overall, this study contributes to a greater understanding of how roe deer (populations) respond to selection pressures across their range. However, further research would be needed to investigate the relationships between immune genetic diversity and pathogen resistance, particularly with respect to habitat variation in biotic and abiotic factors that likely support different pathogen communities and, thus, different pathogen-mediated selection systems. Analysis of MHC variation provides a good framework for studying the local adaptations and genetic health of wildlife populations but provides only a partial understanding of how species adapt to exposure to diseases and/or other stressors [75]. Therefore, studies of other genes associated with innate immunity are needed to provide a more comprehensive understanding of the genetic adaptive potential of species.
Further studies are also needed to investigate the association between the genetic diversity found here and parasite tolerance/resistance; body condition; and ultimately, individual fitness, with the aim of elucidating the mechanisms of balancing selection such as heterozygote advantage and spatiotemporal fluctuating selection [102]. Such research may have important implications for wildlife management and epidemiology [29] but may also help to reduce risks for human and domestic animal health from a One Health perspective.
Supplementary Materials: Supplementary materials are available online at https://www.mdpi.com/ article/10.3390/ani12060723/s1: Figure S1. Genetic structure (from the STRUCTURE) of Slovene roe deer from 10 sampling areas/populations; Figure S2. Results from the STRUCTURE analysis: a) The highest natural log likelihood from STRUCTURE HARVESTER suggests that the number of genetic groups (K) is most likely 1 when considering MHC DRB exon 2 alleles in roe deer female yearlings across Slovenia. b) Analysis for different number of clusters (K = 2-5) of Slovene roe deer, based on MHC; Figure S3. Distribution of positively selected sites on the exon 2 DRB gene of European roe deer estimated by EasyCodeML (model 2); Figure S4. Choice of the number of clusters (K) for discriminant analysis of principle components (DAPC) for MHC supertypes in roe deer shows no point of stabilization or 'elbow' at which the value begins to increase; Figure S5. Plot of the frequency of MHC DRB exon 2 alleles in ten studied geographic groups (i.e., "populations") of European roe deer; Table S1. Basic data on European roe deer females included in the study; Table S2. MHC DRB exon 2 alleles of European roe deer obtained from GenBank; Table S3. MHC DRB exon 2 genetic diversity in roe deer across its distribution range in Slovenia; Table S4. Pairwise FST values and p values (after Bonferroni corrections) between roe deer geographic groups/populations in Slovenia based on MHC DRB exon 2; Table S5. Parameters of analysis of molecular variance of MHC DRB exon 2 data in ten roe deer groups/populations (overall and within three K-clusters) in Slovenia; Table S6. Genetic diversity in Slovenian roe deer populations based on 11 microsatellite loci; Table S7. Genetic diversity in three Slovenian roe deer genetic clusters revealed by STRUCTURE based on 11 microsatellite loci; Table S8. Pairwise FST values and p values (after Bonferroni corrections) between roe deer geographic groups/populations in Slovenia based on microsatellite data; Table S9. Pairwise FST values between roe deer clusters in Slovenia based on microsatellite data; Table S10. Parameters of analysis of molecular variance of microsatellite data in ten roe deer groups/populations (overall and within three K-clusters) in Slovenia. (Reference [103] is cited in the supplementary materials).   Table S1.