High genetic diversity of immunity genes in an expanding population of a highly mobile carnivore, the grey wolf Canis lupus, in Central Europe

The aim of the study was to assess the effect of changes in population size and migration on variation in functional immunity genes in the previously bottlenecked population of the grey wolf, Canis lupus.


| INTRODUC TI ON
Genetic diversity of endangered species has long been an interest of conservation biologists, as it reflects the adaptive potential of a species (Hoelzel et al., 2019). Events such as demographic bottlenecks or limited connectivity between local populations deplete genetic variation, decreasing viability of endangered species (Radwan et al., 2010). The population growth and expansion also affect genetic diversity, often in a counter-intuitive manner. Generally, an expansion may be summarized as multiple founder events resulting in lowered genetic variation in newly established populations and increased genetic variation among populations (Excoffier et al., 2009).
However, the evolutionary consequences of these processes vary considerably depending on the context: for example, a decrease in genetic diversity in edge populations may hinder adaptation to new environmental conditions (Agashe et al., 2011;Bridle & Vines, 2007), whereas increased genetic drift may either accelerate or impede adaptation due to alleles surfing on the expansion front (Burton & Travis, 2008;Klopfstein et al., 2006).
Traditionally, neutral markers such as microsatellites have been used for characterizing genetic conditions of endangered species, but analyses based on such markers do not reflect the footprint that adaptation to local conditions leaves on the genetic diversity. The selection may maintain or even increase genetic diversity counteracting drift and demographic processes (Holderegger et al., 2006) but on the other hand, it may decrease the diversity through fixation of adaptive alleles in local populations. Studies focused on coding regions, encompassing genes important for adaptation, are crucial for recognizing the role of microevolutionary changes providing a theoretical background for the management of populations of endangered species (Mable, 2019). Furthermore, genetic diversity within functional loci may affect ecological processes such as species productivity, population resilience or pathogen resistance (Hughes et al., 2008). In particular, the pathogen resistance seems crucial for species of conservation concern (McCallum & Dobson, 2002).
The genes of the major histocompatibility complex (MHC) have long been used in conservation studies as a marker of functional diversity within a species reflecting its potential to cope with infections (e.g. Castro-Prieto et al., 2011). The MHC genes give also insight into the evolutionary potential of expanding invasive species (Biedrzycka et al., 2020). The importance of genetic diversity in immunity genes is supported by significant associations between certain alleles within immunity genes and susceptibility to infections that have been reported in many species, including wolves (Niskanen et al., 2014). However, MHC accounts for only a fraction of the variation among individuals in pathogen resistance (Jepson et al., 1997) and further studies, encompassing other immunity components, are needed to understand the role of functional variation for threatened species.
Demographic processes, in particular changes in population size, affect allele frequencies and population structure. The impact of demographic bottlenecks on population genetic diversity depends on the history of the declining species. Factors, such as location of the population (insular versus mainland), bottleneck characteristic (its duration and intensity) and prior history of bottlenecks, all are strong determinants of genetic diversity remaining in the threatened species (e.g. Sutton et al., 2015;Taylor & Jamieson, 2008). Strong, longlasting bottlenecks result in an excess of low-frequency variants, while recent, weak and short bottlenecks produce an excess of intermediate frequency variants, which may mimic balancing selection (Cutter, 2019). If a bottleneck is followed by an expansion, genetic variation becomes partitioned between local populations through genetic drift and adaptation to local conditions. Populations colonizing new areas are expected to exhibit lower genetic diversity and to differ in the frequency of genetic variants from the source population because the migrants usually are not a random sample from the source population, and selection for certain traits in successful migrants is likely (reviewed in Excoffier et al., 2009). As the population expansion proceeds, the connectivity between local populations increases, as do the population size what results in increased overall genetic diversity (Dlugosch & Parker, 2008). Nonetheless, our understanding of the evolutionary processes operating in populations that experienced a severe size reduction followed by quick recovery and range expansion is limited, yet it is crucial for proper evaluation and further development of conservation strategies.
In this paper, we characterized variation at seven genes coding components of immune response in a recovering population of the grey wolf Canis lupus in Central Europe. The innate immunity genes were Toll-like receptors (TLRs), and adaptive immunity loci were parts of the major histocompatibility complex, called DLA in dogs.
Wolves are apex predators, playing a crucial role in forest ecosystems in a top-down cascade (Jędrzejewska & Jędrzejewski, 1998). Their interference with humans led to their persecution and eradication from most countries in the 19th and early 20th century. Thanks to the conservation efforts, their numbers have been steadily increasing, and at the end of 20th century, they are recolonizing areas from which the species had been exterminated (Chapron et al., 2014). In Central Europe (Poland, Slovakia, and Czechia), wolves constitute three genetically distinct populations Szewczyk et al., 2019): (a) Baltic, east to the Vistula river; (b) recently established Central European (west to the Vistula river); and (c) Carpathian. The Baltic and Central European populations inhabit lowlands, the latter undergoing rapid expansion colonizing Western Poland (Nowak & Mysłajek, 2016), western Czech Republic (Hulva et al., 2018), Germany (Reinhardt et al., 2019) and areas further west (Andersen et al., 2015;Lelieveld et al., 2016). Recent studies have shown that the Central European population originates from animals inhabiting the western edges of the Baltic population (Szewczyk et al., 2019). Baltic population had never been eradicated, yet it experienced a bottleneck lasting at least 20 years (Flousek et al., 2014;Wolsan et al., 1992), and between 1960 and 1980 its size dropped below 50 individuals (Okarma, 1989(Okarma, , 1993Sumiński, 1975;Wolsan et al., 1992). In the highlands, the Carpathian population is also spreading westwards, and past bottlenecks that happened in the 20th century are visible in contemporary genetic diversity (Hulva et al., 2018). To some degree, the Carpathian population is isolated from both lowland populations, and despite lack of physical barriers dispersals into lowlands are rare Hulva et al., 2018;Szewczyk et al., 2019).
Our aim was to assess the effects of migrations and changes in population size on adaptive genetic variation in previously bottlenecked, large, highly mobile carnivore. We hypothesized the recently established populations, as well as those living close to the continuous species range, to have lower genetic diversity in functional loci compared with the source or long-established populations due to the drift overwhelming selection. Alternatively, we expected the diversity in the functional loci to be maintained by selection and adaptation to local conditions. We also hypothesized the connectivity between local populations to affect the genetic diversity.

| Sample collection
The samples were collected through non-invasive approach (scats, urine, hair, or oestrus blood) during long-term wolf tracking study from 2009 to 2016, focused on areas of known wolf presence Nowak & Mysłajek, 2016;Szewczyk et al., 2019). We also collected tissues from wolves killed in traffic accidents, poached (illegally shot or snared) or found dead due to natural causes. The study also involves 17 samples from wolves hunted legally in Slovakia. All samples used in the current study constituted a subset of 130 samples described by Szewczyk et al., (2019) who studied neutral variance. For the present work, we selected samples of best DNA quality, regardless of their population F I G U R E 1 Location of the samples analysed in the current paper. Colours depict populations recognized by Linnell et al. (2008) and Szewczyk et al. (2019), basing on ecological data and neutral genetic markers. Blue-Baltic population, red-Central European, yellow-Carpathian. Samples genotyped in MHC are marked with diamonds, and samples genotyped in TLR are circles of origin, in order to maximize the number of samples genotyped in the immunity genes. In particular, rather long TLR amplicons require not degraded templates, thus the number of genotyped samples varied between loci (Table S1). Distribution of sampling locations is shown in Figure 1, and further information about samples is provided in the Supplementary Materials, Table S2.

| Amplification of neutral and immunity loci
From faeces and urine, DNA was extracted using QIAamp DNA Stool Mini Kit, and to extract DNA from FTA cards, QIAamp DNA Investigator Kit was used. Ethanol-preserved muscle tissues from roadkills or poached individuals were first washed in water, dried with a paper tissue and processed using Qiagen DNeasy Tissue Kit.

Primers sets and PCR conditions are given in the Supplementary
Materials, Table S3. Amplified fragments were sized using the ABI3730xl sequencer and analysed using Peak Scanner 1.0 software.
MHC alleles were reconstructed using the PHASE algorithm implemented in DNASp v. 5.0 (Librado & Rozas, 2009;Stephens et al., 2001) with default settings: 100 burn-in iterations, 100 iterations and 10 thinning intervals. Using blast, we compared the obtained alleles with the Immuno Polymorphism Database (IPD-MHC, Maccari et al., 2020). DLA haplotypes (DRB1-DQA1-DRB1) were reconstructed in individuals where all three DLA loci were genotyped, following method described by Kennedy et al., (2007): first, we identified haplotypes in homozygous animals, and then, those haplotypes were subtracted from other combinations present in our data.
Four Toll-like receptors (TLR1, TLR2, TLR6, TLR9) were amplified using primers and protocols designed for the dog (House et al., 2009;Mercier et al., 2014), as descibed in Table S3. Amplified fragments (2200-3000bp) were pooled for each individual in similar quantities and purified twice using CleanUp kit (Aabiot). The library was prepared using Nextera XT DNA according to the manufacturer's protocol and sequenced in a single run on Illumina MiSeq using MiSeq Reagent Kit v3 (300 cycles) producing 2x150 bp reads. Resulting pair-end reads were mapped using bwa-mem as described by Kloch et al., (2018), and variants were called in two-step procedure using FreeBayes v1.1.0 (Garrison & Marth, 2012). SNPs that could not be phased based on their physical position on a read were phased computationally using PHASE. Details of the mapping and variant calling are given in Supplementary Materials. Overall, we obtained 6 469 400 reads of an average coverage >50,000× (Table S4). As a reference, we used TLR sequences annotated in the dog genome assembly CanFam3.1 (ENSCAFG00000024010, ENSCAFG00000008351, ENSCAFG 00,000,016,172, ENSCAFG00000023201).
Reading frames for MHC-DLA and TLR loci, and intron/exon structure in TLRs were resolved based on alignment with dog sequences. In DLA, the location of antigen-binding sites (ABS) was determined using an alignment to functionally annotated codons in human HLA (Reche & Reinherz, 2003). In TLRs, using SMART (Simple Modular Architecture Research Tool, Letunic and Bork, 2018) we identified domains, in particular leucine-rich repeats (LRRs) which form functionally important parts of the protein.

| Genetic variation and signatures of selection acting on immunity genes
Basic statistics summarizing polymorphism in the immunity genes, such as the number of variable sites (S), allele diversity (Ad), average number of nucleotide differences (k) and nucleotide diversity per site (π) were calculated using DnaSP. Observed and expected heterozygosity, and F IS values were calculated in Genetix 4.05 (Belkhir et al., 2004). Allelic richness with rarefaction was calculated in R using library pegas (Paradis, 2010). Within DLA class II complex, the statistics were calculated separately for each locus, as they may experience different evolutionary forces shaping their diversity.
We used several tests to infer selection acting on the immunity genes in wolves. First, we analysed departures from neutrality.
Classic neutrality tests, such as Tajima's D (Tajima, 1989), are sensitive to changes in population size. As analysed data included samples from an expanding population, such tests could not be applied for inferring selection. Instead, we used compound neutrality tests proposed by Zeng et al., (2006), Zeng et al., (2007) implemented in software DH (available from the website http://zeng-labve rsus.group. shef.ac.uk) which are robust against demographic assumptions. We used dog sequence as an outgroup, and we performed calculations with default 10 000 coalescent simulations.
Second, we applied codon-based tests which are capable of detecting sites under selection within the sequence based on dN/dS ratio (Kosakovsky Pond & Frost, 2005). For the computations, we used DataMonkey server (Delport et al., 2010). First, we tested for possible recombination events that might have affected the diversity of the studied loci using Genetic Algorithm Recombination Detection (GARD, Kosakovsky Pond et al., 2006), and the identified recombination effects were accounted for in the models. Next, we ran three phylogenetically controlled models of selection: (a) MEME, the Mixed Effects Model of Evolution, which allows for the detection of episodic positive selection (Murrell et al., 2012) Finally, we performed an F ST outlier analysis (Beaumont & Nichols, 1996) in Bayescan (Foll & Gaggiotti, 2008). This software uses a logistic regression model to partition FST coefficients into a population-specific component (beta) and a locus-specific compo-

| Population structure
In the analysis of the population structure, we focused on the management units (MU, henceforth called "populations", delimited by Linnell et al., (2008) and then applied in the assessment of the wolf conservation status in Europe for the IUCN Red Data List of Endangered Species (Boitani, 2018). Recently, Szewczyk et al., (2019) confirmed such an approach using mitochondrial and microsatellite DNA markers and a large set of samples.
First, we analysed the similarities between alleles in each immunity locus using minimal spanning networks implemented in PopArt (Leigh, JW, Bryant D (2015). PopART: Full-feature software for haplotype network construction. Methods Ecol Evol 6(9):1110-1116.). We also constructed haplotype networks for three-locus DLA haplotypes using only individual genotypes in three DLA loci. Differentiation between populations was analysed using AMOVA and pairwise FST between populations calculated in Arlequin v. 3.1 (Excoffier et al., 2005) and discriminant analysis of principal components, DAPC implemented in R package adegenet (Jombard et al., 2010). Next, we used DAPC to explore population structure with no prior assumption about the number of groups. The algorithm groups observations minimizing variation within groups and maximizing between groups. The most probable number of clusters is identified as a minimum in a plot of the Bayesian information criterion (BIC) score of a model versus. number of groups (Jombart et al., 2010). Arlequin and adegenet analyses were also applied on microsatellite loci to provide background measures of the neutral variation.
To assess whether genetic diversity in lowland areas recolonized recently differs from the Baltic population that acts as a source of migrants (Szewczyk et al., 2019), we computed observed heterozygosity and allelic richness separately for areas recolonized in different time (see Figure 1). This part of the analysis included the year when a given area was recolonized (e.g. Nowak & Mysłajek, 2016;, not the year when a sample was collected. Whenever sample size from areas colonized in a given year was low, they were grouped with samples collected in areas recolonized in the following or previous year. Next, we used DAPC to visualise differences in genetic variation between those areas. Our data set contained relatively large number of young adult wolves killed in road traffic collisions. Szewczyk et al., (2019) noticed that dispersers could be overrepresented among this group, probably due to the fact that dispersing wolves often use sub-optimal habitats characterized by high road density .
Thus, we run additional analysis to check for discrepancies between geographic location where the animal was found and genetic assignment of this individual to populations sensu Szewczyk et al., (2019).
Membership analysis was run in STRUCTURE 2.3.4 with the same parameters as described by Szewczyk et al., (2019).

| Genetic polymorphisms and signatures of selection in immunity genes
Overall, we found no evidence for decreased polymorphism of the immunity genes in wolves from Central Europe, despite the Note: Only coding parts of the sequenced genes were considered.
Abbreviations: Adiv, allele diversity; k, average number of nucleotide differences between sequences; n, number of alleles; n, number of individuals; S, number of segregating sites; π, nucleotide diversity per site. Only coding parts of the sequences were analysed.
TA B L E 1 Diversity indices of the immunity loci in the grey wolf in Europe confirmed bottleneck between 1960 and 1980. All DLA alleles were previously reported from wolves (Table S5), and we recognized 11 three-loci DLA haplotypes (Table S6).
Within DLA class II complex, DLA-DQA1 was less polymorphic than DLA-DQB1 and DLA-DRB1 in terms of all measures of genetic diversity (Table 1), and MHC in general had considerably higher number of segregating sites and nucleotide diversity than TLR.
Surprisingly, among TLRs, the least polymorphic locus was TLR2 that had the lowest number of segregating sites, the fewer alleles and the lowest allele diversity (Table 1).
At the population level, the observed heterozygosity in DLA was high, exceeding 0.6 in all populations and all loci. In TLR, the observed heterozygosity was lower than expected in all loci and all populations, yet the low F IS values did not show a considerable deficiency of heterozygotes (Table 2). Generally, populations with higher allelic richness had more private alleles, yet in TLRs, the Carpathian population had more private alleles than the lowland populations ( Table 2).
The compound neutrality tests (Zeng et al., 2007) that are robust against demography, recombination and background selection did not give significant results in either test, although in DLA-DQA1 and DLA-DQB1 the Tajima's D > 2 suggested deficiency of rare alleles (Table S7). We

| Population structure in the immunity loci
The allele networks in all immunity loci consisted of a few frequent alleles shared between populations and several less frequent ones usually unique for a single population ( Figure S1). Such a star-like topology suggests that the frequent alleles located in the centre of the network originated earlier in time, eventually mutating into newer, less frequent alleles. In DLA, the distances between alleles and haplotypes were much higher than in TLR, which reflects the high polymorphism within this gene family.
Generally, the topology of the haplotype networks did not correspond to the geographic locations from where the samples TA B L E 2 Variation at the population level in seven immunity loci compared with neutral variation at 13 microsatellite loci in the grey wolf in Europe  originated. In all loci, most alleles were shared between populations, although in some TLR loci clusters grouping alleles present in a single population were visible. For instance, a cluster of three alleles unique for the Carpathian population was visible in TLR1.
In all genes, AMOVA indicated that the majority of the observed variance (>90%) was at the intra-population level, and the differences between populations accounted for a small fraction (1%-5%) of the observed variation (Table S8). In neutral microsatellite markers, 11% of total variance was explained by the variation among populations, and the overall Fst was much higher than in any immunity locus. Next, we used DAPC to describe population structure best reflecting the genetic variation in the studied loci, with no prior assumption about the number of groups. In both gene families, the model fit indicated by BIC decreased with an increasing number of groups suggesting no structure. In DLA, the BIC curve stabilized after reaching 9 clusters, yet the BIC score reached negative values after 1 cluster. In TLRs, six clusters seemed to best describe the variance in data, but detected clusters did not correspond to their actual geographic origin, as shown in Figure 3b.  Figure 4, the Carpathian-lowland admixture was low, but relatively high proportion of individuals from the Central European population clustered with the Baltic group. These were usually young adult road-killed wolves not related to adjacent family groups which cannot be considered effective migrants, as they died before managing to establish own territories and breed. To make sure that their presence did not affect estimates of the genetic diversity in immunity loci in the Central European population, we re-calculated the analysis for the lowland populations removing the migrants. As shown in Table S10, neither the allelic richness nor the number of private alleles changed in DLA loci after the migrants were removed.
In TLR, the number of private alleles was slightly lower in Central

| Expansion of the lowland population and the genetic variation in immunity genes
In this part of the analysis, we focused on the temporal changes in the variation at immunity genes in the expanding lowland population. First, we compared heterozygosity and allelic richness in areas colonized recently (i.e. after the year 2000) with those parameters from areas with permanent wolf population. Generally, heterozygosity and allelic richness did not differ in MHC loci but decreased in TLR1 and TLR6 in recently colonized areas ( Figure 5).
Using DAPC, we found no differences between areas colonized in different years ( Figure 6). The clades overlapped and their position in the graph did not suggest any directional changes that would be visible if the recently colonized areas suffered from allele loss or founder effect.

| Population structure and variation at immunity genes
In our study, we analysed genetic diversity at seven immunity loci in an expanding population of the grey wolf in Central Europe. Based on neutral markers, wolves inhabiting this area have been assigned to three genetic groups (populations: Baltic, Central European and Carpathian, Szewczyk et al., 2019) corresponding to management units delineated earlier based on geographic criteria (Linnell et al., 2008).
If population structure in immunity genes and neutral microsatellites is similar, this suggests that demographic processes overweight selection. We find no support for this hypothesis.
Despite clear population structure in neutral markers, described by Szewczyk et al., (2019) and confirmed in the current paper, the population structure in immunity loci was weak as indicated by several Thus, further studies are needed to verify whether differences in microbial communities between lowlands and highlands may drive TLR diversity in wolves.
The recently emerged Central European population is smaller than the source Baltic population, and it has been established less F I G U R E 6 Genetic differences between areas colonized in different years in lowland populations (Baltic and Central European) visualized using DAPC algorithm, run separately for each family of the immunity genes than 20 years ago what-considering the wolf generation timegives little time for new functional alleles to emerge and spread in the population. It still can be viewed as a subsample from the Baltic population and thus we expected to find lower genetic diversity in immunity loci in the Central European population, as not all alleles present in the Baltic could be carried west with the migrants. This hypothesis was not supported in the current study. It is important to note, however, that we compared genetic diversity in functional loci between Central European population and westernmost (Polish) part of the Baltic population. Other studies suggested that neutral genetic diversity was higher in areas further east, in Lithuania and Belarus (Szewczyk et al., 2019). Thus, the diversity indices of the Baltic population presented here might be slightly underestimated.
In the present paper, we used a priori definition of populations based on geographic location augmented with studies of microsatellites (Hulva et al., 2018;Szewczyk et al., 2019). The studied populations show some degree of interpopulation admixture-for instance, individuals of lowland origin had been observed in Carpathians (Hulva et al., 2018;Szewczyk et al., 2019). The analysis presented in Figure 4 confirmed that the Carpathian-lowland admixture was low, but there were some individuals from the Central European population clustered with the Baltic group. The analysis of genetic diversity after excluding the potential migrants confirmed high heterozygosity in the Central European population. The migrants may introduce new alleles to a population leading to a decreased population structure. However, this is possible only if demographic events outweigh selection, and we found no such evidence in our study, where studied loci differed in levels of polymorphism.
The measures of genetic diversity in Toll-like receptors were comparable to studies from other mammals, with allele diversity exceeding 0.6 (e.g. Kloch et al., 2018;Quéméré et al., 2015).
Heterozygosity in TLR loci was high compared with genome-wide estimates of the diversity which for Eastern European wolves is 0.21-0.24 based on 61K SNP set (Pilot et al., 2014). The only exception in the current paper was TLR2, but this result can be explained by a lower number of samples genotyped in this locus. TLR variation has not yet been characterized in wolves in the context of population genetics. In a paper reporting SNP diversity in canid TLRs, Cuscó et al., (2014) included 335 dogs of various breeds and 100 wolves from Italy and Russia. Only 6.25% of the non-synonymous variants were wolf-specific which suggests trans-species polymorphism.
Allele frequencies differed between wolves and dogs, and in wolves, the frequency of potentially damaging or deleterious variants was much lower, what may indicate that selective pressure associated with pathogens acts stronger on free-living wolves compared with domesticated dog breeds.
Genetic variability in DLA class II complex has been reported in several wolf populations in Europe. We found from 6 to 9 alleles per locus which are comparable to studies from areas maintaining high wolf population size (Niskanen et al., 2014) and higher than in southern Europe where wolf populations have been isolated for a long time (Galaverni et al., 2013;Rocha et al., 2019). Frequent DLA-DQA1 and DLA-DQB1 alleles reported in our study were also abundant in Italy, Spain, Finland and Russia (Galaverni et al., 2013;Niskanen et al., 2014;Rocha et al., 2019). Interestingly, DLA-DRB1 alleles frequent in Italy (Galaverni et al., 2013) and Iberian Peninsula (Rocha et al. 2019) were either rare or absent in our study area which suggests possible northern origin of those variants. This is not surprising, as the Central European population originated from the Baltic population inhabiting north-eastern Poland and Baltic countries.
Previous studies of DLA variation in wolves in Europe encompassed isolated populations in Iberia, and Italy (Arbanasić et al., 2012, Galaverni et al., 2013Rocha et al., 2019) except for Finland and Karelia (Niskanen et al., 2014). On the other hand, scenario described in the current paper was similar to the study from Finland (Niskanen et al., 2014), where local wolves are supplied by migrants from the Baltic population, which is connected to the large Asian wolf population.
Contrary to our predictions, we did not find evidence for diminished genetic diversity in the recently established Central European population, nor in the western edge of the Baltic population, despite lowered heterozygosity in TLR1 and TLR6 in areas colonized after 2015. High genetic diversity in the westernmost part of the Baltic population may be explained by its good connectivity with areas located further east, where wolves are widespread (Sieber et al., 2015;Stronen et al., 2013) and are able to disperse over long distances (Byrne et al., 2018) due to the less dense human infrastructure than in Western Europe (Ibsch et al., 2016). In contrast, the area inhabited by the Central European population constituting the westernmost edge of continuous wolf range in Northern European Plains (Chapron et al., 2014) consists mainly of the human-dominated landscapes Reinhardt et al., 2019) that are supposed to limit dispersal (Tucker et al., 2018). However, ecological studies revealed that in Western and Central Europe, the recolonization was initialized by individuals which dispersed over highly fragmented landscapes with numerous anthropogenic barriers (Huck et al., 2010(Huck et al., , 2011, and established family groups in areas distant from the source populations (Nowak & Mysłajek, 2016;Reinhardt et al., 2019).
Recent studies showed that dispersal tendency is associated with certain genetically based traits such as boldness (Cote et al., 2010) which may facilitate the settlement in human-altered habitats.
Moreover, the survival of the pups in Central European population is higher than in both Baltic and Carpathian populations (Jędrzejewska et al., 1996;Nowak & Mysłajek, 2016;Nowak et al., 2008). This suggests the high quality of migrants and pre-saturation dispersal (Lidicker, 1975). Our findings reporting their high genetic diversity in functional loci involved in pathogen resistance support this claim.

| Selection acting on immunity genes in wolves
The evidence from multiple species shows that generally the TLRs are under purifying selection but specific TLR codons show evidence of positive selection (Alcaide & Edwards, 2011;Areal et al., 2011;Fornůsková et al., 2013;Grueber et al., 2013;Quéméré et al., 2015).
In wolves from central Asia, purifying selection dominated TLR evolution, although several positively selected codons were detected in TLR1, 2, 6 and 9 (Liu et al., 2017). Relatively high TLR variation was shown in coyotes and critically endangered red wolf. In both species, signatures of purifying selection were found in TLR6, and alleles in this locus were associated with Ehrlichia exposure (Brzeski, 2015).
Similarly, we found that most codons are affected by negative per- MHC genes usually exhibit high levels of trans-species polymorphism (TSP) as a result of long-term persistence of allelic lineages (Klein et al., 2007). In TLRs, the TSP is rare, although the genes from pattern recognition receptor (PRR) family, despite being conserved, are the candidates for TSP as considerable non-synonymous polymorphism with predicted functional significance has been recently documented both on interspecific and intraspecific levels (Těšický & Vinkler, 2015). In particular, the same TLR alleles were detected in wolves and in dogs, but in each species they differed in frequencies (Cuscó et al., 2014). In Europe, wolves may cross with dogs, yet the levels of introgression are generally low (<7%, Dufresnes et al., 2019). The hybridization usually occurs locally, where access to potential wolf mates is limited and the number of free ranging dogs is high (Fabbri et al., 2014;Hurford et al., 2006;Wabakken et al., 2001).
Thus, we do not expect that possible hybridization affects significantly the variation within immunity genes in wolves. Nonetheless, some level of adaptive introgression cannot be excluded, because it is not possible to distinguish alleles of trans-specific and hybrid origin from those present only in a given species. Even though wolfdog hybrids may be identified using microsatellites, and we did so in the current study making sure no hybrids were genotyped in immunity loci, the distinction is only possible for F1. Generally, identification of back-crosses and hybrids is not reliable beyond second generation, even using large SNP panels (Godinho et al., 2015;Pilot et al., 2018). Thus, the role of introgression in maintaining diversity in immunity genes in wolves remains unclear.

| Significance of studying loci under selection in species of conservation concern
Contrary to expectations, in the present study we did not find evidence for a selection acting on neither immunity loci using compound neutrality tests except for significant Tajima's D > 2 in DLA-DQA and DLA-DQB. However, we identified codons under selection, both positive and negative. Significant values of Tajima's D and Gu and Li D* have been reported in Finland (Niskanen et al., 2014) and in the Iberian Peninsula (Rocha et al., 2019). These findings should be treated with caution, as the significant values may indicate both, selection and changes in population size. As we a priori knew that the size of studied population is not constant, we used compound neutrality tests that are robust against demography, therefore suitable for an expanding population (Zeng et al., 2007), yet none of them was significant. A number of recent studies provide evidence suggesting that drift is the main force shaping diversity at TLR loci in recently bottlenecked or expanding populations (Grueber et al., 2013;Hartmann et al., 2014;Quéméré et al., 2015).
Testing for selection may be difficult under different demographic scenarios, in particular in populations colonizing new areas (White & Perkins, 2012). Unfortunately, the impact of changes in population size on genetic variation mimics the effect of selection (Cutter, 2019). In particular, the population expansion results in increased Ne and skewed distribution of site frequency spectrum resembling the effects of the positive or negative selection (Cutter, 2019). The balancing selection, which is expected to act on immunity genes, gives an opposite effect promoting intermediate allele frequencies. Thus, the effects of demography and selection may cancel resulting in no significant departures from neutrality. In the current study, we found signal of selection in most loci at a single codon level, similarly to Niskanen et al., 2014. The ratio of non-synonymous to synonymous substitutions used in codonbased tests reflects past rather than contemporary selection (Garrigan & Hedrick, 2003). Most codon-based tests used in the current study detected pervasive, that is past selection, especially in TLR genes. Thus, it seems easier to trace and detect footprints of historical selection in populations that undergone different demographic scenarios, while testing for contemporary selection in expanding species is challenging (White & Perkins, 2012).
We were also not able to detect any outlier loci based on F ST statistic. Similarly as in the case of previous selection tests, the detection of outliers may be affected by non-equilibrium demographic histories. Our results do not discount the importance of selection acting at functional loci such as TLRs or DLA, but rather highlight the effects of demography on allelic variation. We also show that neutrality tests do not seem to be a suitable tool for inferring selection, and codon-based tests, augmented with other methods such as comparing haplotype networks with neutral spatial population structure, should be used to determine the variability of functional loci in species of conservation concern with high dispersion abilities. would like to thank the Reviewers, in particular Reviewer 3, whose work greatly improved the manuscript.

CO N FLI C T O F I NTE R E S T
Authors declare no conflict of interest.

PEER R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13360.