Human‐modified canids in human‐modified landscapes: The evolutionary consequences of hybridization for grey wolves and free‐ranging domestic dogs

Abstract Introgressive hybridization between domestic animals and their wild relatives is an indirect form of human‐induced evolution, altering gene pools and phenotypic traits of wild and domestic populations. Although this process is well documented in many taxa, its evolutionary consequences are poorly understood. In this study, we assess introgression patterns in admixed populations of Eurasian wolves and free‐ranging domestic dogs (FRDs), identifying chromosomal regions with significantly overrepresented hybrid ancestry and assessing whether genes located within these regions show signatures of selection. Although the dog admixture proportion in West Eurasian wolves (2.7%) was greater than the wolf admixture proportion in FRDs (0.75%), the number and average length of chromosomal blocks showing significant overrepresentation of hybrid ancestry were smaller in wolves than FRDs. In wolves, 6% of genes located within these blocks showed signatures of positive selection compared to 23% in FRDs. We found that introgression from wolves may provide a considerable adaptive advantage to FRDs, counterbalancing some of the negative effects of domestication, which can include reduced genetic diversity and excessive tameness. In wolves, introgression from FRDs is mostly driven by drift, with a small number of positively selected genes associated with brain function and behaviour. The predominance of drift may be the consequence of small effective size of wolf populations, which reduces efficiency of selection for weakly advantageous or against weakly disadvantageous introgressed variants. Small wolf population sizes result largely from human‐induced habitat loss and hunting, thus linking introgression rates to anthropogenic processes. Our results imply that maintenance of large population sizes should be an important element of wolf management strategies aimed at reducing introgression rates of dog‐derived variants.


| INTRODUC TI ON
Domestication is a striking example of human-induced evolution.
Adaptation to new ecological niches created by humans as well as artificial selection pressures have resulted in extensive changes in domesticated organisms, including morphology, physiology, reproduction and trophic position (Milla et al., 2015;Solberg et al., 2020;Wilkins et al., 2014). The effects of these changes are not limited to domesticated organisms, but extend into the entire biosphere (Larson & Fuller, 2014). In particular, domesticated organisms have a profound effect on their wild relatives through resource competition, pathogen transfer and hybridization (Lescureux & Linnell, 2014;Turcotte et al., 2017). This effect is becoming increasingly pronounced with the ongoing population growth of domestic animals and plants (Foley et al., 2005;Gompper, 2014;Thornton, 2010).
Introgressive hybridization between domestic animals and their wild relatives has an important effect on gene pools and phenotypic traits of both groups. This process can be considered as a form of human-induced evolution, as it results from a combination of humaninduced processes that have taken place on different timescales, from ancient domestication processes to recent ecosystem changes affecting distribution, density and patterns of interbreeding between sympatric or parapatric taxa (Crispo et al., 2011;Grabenstein & Taylor, 2018). Although the process of hybridization between wild and domestic taxa is extensively documented (reviewed in McFarlane & Pemberton, 2019), its evolutionary consequences are poorly understood. In general, introgressive hybridization can have negative consequences such as the loss of unique adaptations or extinction via genetic swamping (Todesco et al., 2016). Alternatively, there may be positive consequences such as transmission of adaptive variation between species (e.g., Jones et al., 2018;Oziolor, 2019;Song et al., 2011), genetic rescue (Whiteley et al., 2015) and adaptive evolution (when hybridization-derived traits facilitate adaptation to novel environmental conditions; Hedrick, 2013).
Introgression from domesticated taxa is thought to have predominantly negative consequences for wild populations, given that genetic and phenotypic variation resulting from domestication may be deleterious for wild animals, and numerical preponderance of domestic animals over their wild relatives facilitates genetic swamping. In some instances, however, variants of immune system genes originating from domestic animals were shown to have positive fitness effects in their wild relatives. In North American grey wolves, a dog-derived variant of beta-defensin 103 immune system gene has increased in frequency and geographic range since its introduction through hybridization more than 1500 years ago, and currently shows signature of balancing selection in several populations (Anderson et al., 2009;Schweizer et al., 2018). Fitness measures such as lifespan and lifetime reproductive success show strong selective advantage for heterozygous individuals carrying the dog-derived variant (Coulson et al., 2011), which is likely related to an enhanced immune response (Schweizer et al., 2018). The Alpine ibex (Capra ibex ibex) has only two MHC DRB exon 2 alleles, one of which was shown to originate from domestic goats (Capra aegagrus hircus) (Grossen et al., 2014). The resulting improved immune response could have contributed to the successful reintroduction of this species following its near extinction (Grossen et al., 2014).
Although these examples show that genetic variation transferred from domestic animals may be adaptive in their wild relatives in some circumstances, it is unknown how anthropogenic introgression affects genome-wide adaptive variation in wild populations.
The evolutionary consequences of introgression from wild populations to their domesticated relatives are also largely unknown.
Cross-breeding populations of the grey wolf (Canis lupus) and the domestic dog (Canis lupus familiaris) provide an excellent model system for studying the consequences of anthropogenic introgression.
As the first domesticated species (Larson & Fuller, 2014), the dog has been living in human-modified habitats longer than any other domestic animal and acquired multiple adaptations to living as a human commensal. A notable example is the increased copy number of the alpha-2B-amylase (AMY2B) gene in comparison with wolves, which has facilitated digestion of starch-rich food (Axelsson et al., 2013) following the spread of prehistoric agriculture (Arendt et al., 2016).
The size of the dog population is positively correlated with the size of the human population and currently reaches one billion individuals worldwide, of which about 75% are free-ranging (Gompper, 2014).
In continental parts of the European Union, the number of dogs in rural areas was estimated at 18.4 million (Gompper, 2014), while the number of wolves in the same region at the same time was estimated at 12,000 (Chapron et al., 2014), resulting in a dog to wolf ratio of approximately 1500:1. This implies that most wolves are likely to encounter a dog during their lifetime. Such encounters can result in dogs being killed by wolves, or in some instances in mating and production of hybrid offspring (Lescureux & Linnell, 2014).
Wolf-dog hybrids are fertile and can reproduce with both wolves and dogs, resulting in introgression within both wolf and dog populations (e.g., Hindrikson et al., 2017;Kopaliani et al., 2014;Pilot et al., 2019). introgression from Tibetan grey wolves (Miao et al., 2017). The notion that wolf admixture can have beneficial effects for dogs resulted in intentional cross-breeding, which was practised by human societies across the world (Lescureux, 2018).
In contrast, introgression from domestic dogs into grey wolves is considered a conservation threat, as the accumulation of dogderived gene variants in wolves' gene pools can lead to a gradual loss of genetic distinctiveness and unique adaptive variation (Donfrancesco et al., 2019;Hindrikson et al., 2017). A general expectation exists that heritable traits originating from domestic dogs are maladaptive for wild wolves living in their natural ecosystems, since dogs have been subject to artificial selection and are adapted to the ecological niche of human commensal. However, for wolves living in human-dominated landscapes, some traits originating from domestic dogs may become advantageous. Currently, little is known about the effects of introgression of dog-derived variants into the gene pool of wolf populations, with the exception of adaptive introgression of a dog-derived CBD103 variant in North American wolves described above (Anderson et al., 2009;Schweizer et al., 2018).
The most likely sources of such introgression are free-ranging dogs (FRDs) inhabiting rural areas. Eurasian FRDs are not a product of admixture between breeds, but constitute a distinct and older genetic group Shannon et al., 2015). Although FRDs depend on anthropogenic food, they can survive without having any continued direct interaction with humans, they have no constraints in mate choice and therefore are not subject to ongoing artificial selection.
Accordingly, their genomes display lower levels of deleterious genetic variation than pure-bred dogs (Marsden et al., 2016). FRDs and purebred dogs show signatures of diversifying selection in genes related to reproduction, immunity and chemosensory perception, which may reflect adaptations of FRDs to independent survival (Pilot et al., 2016).
Some adaptive traits of FRDs may become beneficial to wild canids living in human-dominated landscapes, if acquired via hybridization.
For example, wolves living in habitats heavily transformed by humans are likely to be increasingly exposed to encounters with domestic dogs and dog-derived immune system gene variants could enhance immunity of wolves to dog-derived infectious diseases.
The aim of this study is to (a) assess introgression patterns in Eurasian populations of wolves and FRDs, (b) identify chromosomal regions with a significant deficiency or excess of introgressed ancestry in both canids and (c) assess whether genes placed within these regions are under selection, and are associated with phenotypic traits which differ between wolves and FRDs and for which introgression may result in increased fitness.

| Dataset
This study focused on Eurasian grey wolves and FRDs, known from earlier studies to have undergone introgressive hybridization (e.g., Fan et al., 2016;Hindrikson et al., 2017;Pilot et al., 2015). The dataset analysed in this study was obtained by merging genome-wide SNP genotypes from worldwide populations of grey wolves and FRDs as well as pure-bred domestic dogs from publicly available datasets Fitak et al., 2018;Frantz et al., 2016;Pilot et al., 2015Pilot et al., , 2019Stronen et al., 2015;Vaysse et al., 2011;Vernau et al., 2013). Although this study was focused on wolf-dog hybridization in Eurasia, we included data from North American grey wolves and purebred dogs in the analysed dataset to ensure sufficient representation of nonadmixed individuals. North American wolves show very limited genome-wide introgression from dogs Pilot et al., 2018) and most dog breeds of European origin show very limited or no wolf admixture vonHoldt et al., 2010).
All datasets included in the analyses were generated using the CanineHD Whole-Genome Genotyping BeadChip (Illumina).
The population structure analysis for the merged dataset using Admixture (Alexander et al., 2009) showed that individuals originating from different datasets but representing the same dog breeds or wolf populations cluster together, demonstrating that the merging was performed correctly. The CanineHD BeadChip produces genotypes at 167,989 autosomal SNP loci and 5660 X chromosome SNP loci. Some of the published datasets reported a reduced set of loci compared with the total number included in the CanineHD BeadChip, and/or did not report X chromosome loci. Therefore, the final merged and pruned dataset included data for 106,549 autosomal SNP loci. This dataset included 1526 individuals (506 wolves and 1020 dogs), each with a genotyping rate above 90% (Table S1).
With a few exceptions, sampling for this dataset does not cover exactly the same regions for wolves and FRDs ( Figure 1). However, Eurasian free-ranging domestic dogs originate from a recent geographic expansion, dated at about 15,000 years ago (Wang, Zhai, et al., 2016) and do not show strong population structure across Eurasia . Modern grey wolves were also shown to originate from an expansion of a single source population that began approximately 25,000 years ago (Loog et al., 2020) and display a limited population structure across Eurasia (e.g., Ersmark et al., 2016;Pilot, Dąbrowski, et al., 2014;Pilot et al., 2019), with the exception of Indian and Himalayan wolves (Sharma et al., 2004), that are not included in this research. Given that we have a broad geographic coverage for both wolves and dogs, our dataset provides a strong representation of genetic variability of both canids across Eurasia. Therefore, the inference of admixture between wolves and dogs should not be affected by the lack of exact geographic overlap between their sampling sites.
The Illumina Canine HD BeadChip was developed in collaboration with the LUPA Consortium and was based on a dataset that, in addition to dog breeds, included 15 grey wolves from Europe and North America . The wolves showed 118,256 segregating sites, of which 1471 sites were variable in wolves but not in any dog breed analysed. Because of the small number of wolves included in the reference panel, loci that show variation in wolves but not dogs are underrepresented in the chip, which may affect the admixture inference. However, the Canine HD BeadChip was previously used in studies on dog introgression into wolves (e.g., Galaverni et al., 2017) and wolf introgression into dogs (Caniglia et al., 2018).
This last study was focused on Czechoslovakian wolfdogs, and the inferred admixture proportions were consistent with the breed history. Moreover, the inference of admixture proportions in both wolf and dog populations based on the Canine HD BeadChip  is consistent with the inference from whole-genome data Freedman et al., 2014). Nevertheless, we note that the underrepresentation of genetic variation specific to wolves may result in the reduced power to detect wolf introgression into dogs compared with the reverse.

| Ancestry block analysis
Plink software (www.cog-genom ics.org/plink/ 1.9; Chang et al., 2015) was used to remove from the dataset SNPs with low variability (MAF < 0.01) and those with >10% missing data. For the purpose of the analysis in lAmP (see below), we also removed SNPs in strong linkage disequilibrium (LD; r 2 > 0.1). We assessed the presence of chromosomal ancestry blocks derived from hybridization in individual canids in two steps, initially using the software lAmP (Sankararaman et al., 2008) followed by elAi (Guan, 2014). lAmP allows ancestry block estimation without defining a priori ancestral nonadmixed populations. This unique feature was required for the analysis of our dataset, because most Eurasian wolf populations include individuals showing signature of past admixture with dogs Pilot et al., 2018), and some FRD populations and dog breeds show evidence of admixture with wolves Pilot et al., 2015). In lAmP, identification of ancestral populations was an integrated part of the admixture analysis.
In the lAmP analysis, we assumed a mixture proportion of 0.33:0.67, which was determined based on the frequency of wolves versus dogs in the dataset (506 vs. 1020 individuals). This ratio accommodated a conservative scenario where none of the individuals is admixed. We used a recombination rate of 1 × 10 −10 per base pair per generation and fraction of overlap between adjacent windows (offset) of 0.2. We assumed 10 generations since admixture, because power to detect recent admixture was diminished when more distant admixture events were assumed (Pilot et al., 2018).
Based on the lAmP results, we identified individuals that were not admixed or had very low levels of inferred admixture. These individuals were used to define reference genotype sets for wolves F I G U R E 1 Distribution of samples of Eurasian grey wolves (red circles) and free-ranging dogs (green circles) analysed in this study. Geographic locations of samples are precise except Mongolian wolves and free-ranging dogs from China and Portugal, which have approximate locations. The number of samples collected from the same locations is reflected by the circle size. The introgression pattern analysis was carried out for West Eurasian wolves (marked with the red frame) and all free-ranging dogs shown on the map that carried introgressed chromosomal fragments. Among East Asian wolves (black frame), only four admixed individuals were found, including an F1 hybrid. Grey wolf geographic range drawn according to Boitani et al. (2018) and Wang, Ma, et al. (2016) and dogs, which are required by the elAi analysis to estimate allele frequencies in nonadmixed populations. Thresholds to classify individuals as nonadmixed were set at 0.001 for dogs and 0.01 for wolves. The different thresholds account for different distributions of admixture proportions in wolves and dogs and were established so as to obtain at least 300 individuals in each nonadmixed reference sets. Due to the use of these thresholds, a small number of introgressed fragments present in the reference genotypes were considered as native and thus could not be identified as introgressed blocks in the individuals tested using elAi, potentially leading to the underestimation of admixture. The alternative of including in the reference genotype sets only individuals that were identified in lAmP as having no introgressed fragments would result in the reference genotype sets being smaller and not representative of all regional genetic variation (i.e., some wolf populations would not be represented and the dog reference set would consist mostly of pure-bred dogs). Using such reference sets could have resulted in the overestimation of admixture, and therefore we used the first, more conservative approach.
The resulting reference datasets included 304 nonadmixed dogs and 334 nonadmixed wolves, representing all primary geographic populations and breed groups studied. The remaining 717 dogs and 171 wolves (888 individuals), which were not included in the reference genotype sets, were further tested for admixture using elAi.
Although we used the set on nonadmixed individuals identified in lAmP as input for the elAi analysis, these analyses were otherwise independent. If lAmP incorrectly inferred nonexistent admixture, it was possible for elAi to identify all individuals as pure wolves and dogs. Accordingly, if lAmP underestimated admixture, it was possible for elAi to identify considerably higher admixture proportions compared with lAmP.
We used 20 expected maximization steps to estimate the parameters of the hidden Markov model implemented in elAi. elAi accounts for the possibility of continuous admixture throughout multiple generations, enabling a more realistic representation of wolf-dog admixture. Unlike lAmP, it does not require filtering for LD, and thus, ancestry could be inferred for all SNPs in the initial dataset (106,549 autosomal SNPs). elAi can detect ancestry blocks <1 cM and accounts for the presence of population substructure within each of the admixing entities. In our analysis, we assumed admixture between two main population clusters (wolves vs. dogs) during 100 generations, and the presence of 10 lower-layer clusters (advised to be five times the number of upper-level clusters).
Further, elAi does not require phasing or a recombination map, because it directly estimates cluster-switch rates between adjacent markers, which enables the direct inference of recombination rates at each locus from the data analysed (Guan, 2014). There is considerable variation in recombination rates between sexes and between individuals within species (e.g., Kong et al., 2010), and therefore, even a high-resolution recombination map will not necessarily be accurate if the analysed dataset is considerably different than the dataset used to construct the map. Hence, we selected elAi as the preferred software for the admixture analysis.

| Assessment of the effect of recombination rate on estimated admixture proportions
To assess whether the variation in local admixture proportions across the genome (inferred using elAi) is dependent on recombination rates, we tested for the correlation between these two variables, first for all SNP loci across 38 autosomal chromosomes and then for each chromosome separately. We used the average recombination rates for males and females provided by Campbell et al. (2016). We selected the recombination map from this study over other available high-density maps (Auton et al., 2013;Axelsson et al., 2012), because (a) this map was constructed using the Illumina CanineHD Whole-Genome Genotyping BeadChip, which was the same array as in our study, and (b) it is a pedigree-based map, which is expected to be more accurate than LD-based maps that can be affected by demographic patterns specific to the study populations. We used R 4.0.3 (R Core Team, 2020) to calculate Pearson's correlation coefficient and to fit a linear regression model to the data.

| Detection of chromosomal fragments with overrepresentation of introgressed variants
The elAi analysis included only 717 dogs and 171 wolves for which the admixture proportions estimated in lAmP were above the established thresholds. Pure-bred dogs were excluded from further analysis, as we were interested in natural introgression patterns only.
Both East Asian and North American wolves were also excluded due to small numbers of admixed individuals detected (four in each population). We excluded Mexican wolves as well, even though many individuals carried a small proportion of dog admixture, since this population is highly inbred (Fredrickson et al., 2007), and admixture may have resulted from a single event, possibly during captive breeding.
Analyses of introgression patterns focused on the remaining set of admixed individuals, which comprised 88 West Eurasian wolves and 201 Eurasian FRDs ( Figure 1). Based on elAi results, we calculated the mean admixture proportions within each autosomal chromosome and across autosomal chromosomes in both datasets. We identified chromosomal blocks with hybrid ancestry either greater or smaller than three standard deviations (SD) from the mean for each chromosome (see Figure S1), thus permitting identification of blocks with overrepresented or underrepresented hybrid ancestry, respectively. To reduce the false-positive rate, we only considered ancestry blocks that included at least 10 sequential SNPs.
Blocks with overrepresented or underrepresented hybrid ancestry were identified based on the SD at the level of individual chromosomes rather than the global SD across all autosomal loci, because chromosomes are natural genetic units with independent recombination. This approach is consistent with the previous step of the study, that is the detection of admixture and reconstruction of the distribution of introgressed blocks, which was done at the level of individual chromosomes. Accordingly, the detection of selection signatures was also based on patterns of extended haplotype homozygosity within individual chromosomes (see below), so our approaches for all the consecutive analyses were consistent. For comparative purposes, we also carried out the identification of overrepresented or underrepresented hybrid ancestry based on the global SD estimates, which is reported in the Supplementary Materials.
To estimate the false-positive and false-negative rate associated with the detection of blocks with overrepresented hybrid ancestry, we applied a random resampling approach. We randomly selected 38 chromosomal blocks by choosing a position within a selected autosomal chromosome (both determined randomly). A number of consecutive SNPs to be included within the block (counting from the selected position) was also chosen randomly from a range of 10 to 150. These 38 blocks were then assembled into a "randomised chromosome." This was done separately for wolves and FRDs. The "randomised chromosome" is representative of autosomal chromosomes from their population of origin and reflects the population's admixture proportions as well as the demographic processes and evolutionary forces that shaped the gene pool, but is free from chromosome-specific recombination patterns that could result from incomplete sampling of the parental population genetic diversity. We analysed introgression patterns in these "randomised chromosomes" in the same way as real chromosomes in order to identify chromosomal blocks with significantly overrepresented hybrid ancestry. We then estimated error rates by assessing whether (a) blocks with significantly overrepresented hybrid ancestry identified within a "randomised chromosome" were identified as having no overrepresented hybrid ancestry in the real chromosomes (falsenegative rate), and (b) blocks with no significantly overrepresented hybrid ancestry in the context of a "randomised chromosome" were identified as having overrepresented hybrid ancestry in the context of their real chromosomes (false-positive rate).

| Identification of loci under positive selection
The genotype data for each chromosome that contained ancestry blocks showing overrepresentation of introgressed alleles were phased using fastPHASE (Scheet & Stephens, 2006). For the phased data, we performed scans for signatures of selection across autosomal chromosomes using the Integrated Haplotype Homozygosity Score (iHS) statistics (Voight et al., 2006), as implemented in the R package rehh v. 3.1.2 (Gautier et al., 2017). The iHS statistics is derived from the extended haplotype homozygosity (EHH) statistics, which measures homozygosity decay in haplotypes carrying a specified "core" SNP at one end, with increasing haplotype length (Sabeti et al., 2002). An allele that rises rapidly in frequency due to selection will have high levels of haplotype homozygosity extending over a larger distance than expected under a neutral model (Sabeti et al., 2002). The Integrated Haplotype Homozygosity Score (iHS) is based on the integral of the observed EHH decay away from a specified core SNP until it reaches the value of 0.05 (Voight et al., 2006). We used the rehh package's functions scan_hh and ihh2ihs to calculate the iHS statistics and its two-sided p-value for each SNP.
Phasing and rehh analysis were conducted separately for wolf and FRD genotypes.
Introgression can change the distribution of allele frequencies and haplotype structure, which can potentially interfere with the detection of selection signatures using iHS (Booker et al., 2017).
However, while introgression can distort haplotype structure in any part of the genome, selection can only act on those parts of the genome that have a biological function (including protein-coding genes, long noncoding RNA and regulatory elements) and DNA regions in local proximity that may be affected by selective sweeps. Therefore, we assessed the error rate associated with the detection of adaptive introgression as the proportion of SNPs that are inferred to be under selection, but are located a sufficient distance from protein-coding genes and long noncoding RNA that they are unlikely to be affected by a selective sweep. We assumed this distance to be 100 kb, and functional information was based on the Ensemble CanFam3.1 dog genome annotation.
We considered the iHS results as significant if p < 0.05 and |iHS| > 2 and did not use a correction for multiple testing. The study that introduced the iHS statistics showed that |iHS| >2 is a powerful criterion to identify signals of selection and did not use or recommend corrections for multiple testing (Voight et al., 2006).
Here, we focused on signatures of selection within relatively short chromosomal blocks (161-4493 kb) showing overrepresentation of introgressed variants and did not attempt to identify loci showing signatures of selection across the entire genome. The signal of selection is not independent for individual SNPs, given that the iHS test is based on haplotype homozygosity. Therefore, loci showing signatures of selection are expected to be clustered (e.g., we found up to 14 significant SNPs within one gene). For these reasons, we report candidate genes identified using this test without correcting for multiple testing. However, in order to show how correction for multiple testing could affect our interpretation, we also report the candidate genes identified after applying the Bonferroni correction, based on the number of SNPs within each chromosome that were used to obtain the iHS statistics.

| Gene ontology enrichment analysis
Ensembl was used to identify coding genes located within chromosomal blocks displaying overrepresentation of hybridization-derived variants, based on CanFam3.1 dog genome assembly. We considered a gene as being located within a chromosomal block if the entire open reading frame or only a part of the reading frame was located within that block. We also generated an additional set of genes of interest by identifying human genes orthologous to canine genes found within those blocks using synteny analysis. The two gene sets were kept separate for further analyses. As each gene set was created separately for wolves and FRDs, we generated a total of four gene sets.
Each gene set was tested for Gene Ontology (GO) term enrichment using the web-based software g:Profiler (Raudvere et al.,

| 2439
PILOT eT aL. 2019). This analysis was carried out separately for wolves and FRDs. Genes identified based on the dog genome assembly were compared to the reference set of annotated genes from the dog genome. Human orthologues were compared to the reference set of annotated genes from the human genome assembly (GRCh38.p13).
A significance threshold of 0.05 with Benjamini-Hochberg correction was used, as well as the more conservative g:SCS (Set Counts and Sizes) false discovery rate correction method that accounts for multiple testing due to the overlap of functional terms (Reimand et al., 2007). We found differences in the frequency of individuals carrying introgressed ancestry blocks between geographic regions in both wolves and FRDs. In the case of wolves, the highest frequency of such individuals was found in Europe, lower in the Caucasus and Mongolia and the lowest in Yakutia. In the case of FRDs, high frequency of dogs carrying hybridization-derived variants was found in several geographically distinct regions, including Saudi Arabia, Mongolia and South-East Asia ( Figure S2).

| The effect of the recombination rate on admixture proportions
In the regression models fitted to the set of SNPs from across all autosomal chromosomes, we found no significant correlation be-  (Table S2).
We also assessed the correlation between the local admixture proportions in wolves versus dogs. If recombination had a significant effect on admixture proportions, the same effect should be expected in wolves and dogs, resulting in a positive correlation between their local admixture proportions. Instead, we found a weak, but significant negative Pearson's correlation (R = −0.042, 11,224 df, p = 7.384 × 10 −6 ) between the local admixture proportions in wolves versus dogs. The chromosome-level analysis did not reveal a consistent pattern: a significant negative correlation between local admixture proportions was found for 13 chromosomes, a significant positive correlation for 11 chromosomes and no significant correlation for the remaining 14 chromosomes (Table S2).

| Chromosomal blocks with significantly overrepresented hybridization-derived ancestry
Although average genome-wide admixture proportion was higher in West Eurasian wolves than Eurasian FRDs, the opposite pattern was observed with regard to the number and length of chromosomal blocks with significantly overrepresented hybridization-derived ancestry. This was only the case when outlier blocks were identified at the level of individual chromosomes; when the standard deviation was calculated across all autosomal chromosomes, a larger number of significantly overrepresented introgressed blocks was identified in wolves than in dogs (Table S3). This result was, however, biased by the presence of several regions in the dog genome with wolf admixture proportions 5-15 times higher than the global mean. As a result, the global standard deviation across all loci (0.23) was considerably higher than the mean standard deviation across 38 within-chromosome means (0.17). This prevented the detection of local outliers, with the exception of those few with extreme values that reached the global threshold. Therefore, we report the results of outlier block identification based on the global standard deviation in the Supplementary Materials (Table S3, Figures S3 and   S4), while the results reported below are based on the detection of outlier blocks within autosomal chromosomes. It should be noted that we did not analyse the admixture patterns within the X chromosome, because our dataset did not include X chromosome SNPs (see Section 2.1).
In West Eurasian wolves, we identified 16 blocks with overrep-  hybrid ancestry blocks) included 61 genes in wolves and 294 in FRDs.
Synteny analysis with the human genome identified larger numbers: 72 in wolves and 311 in FRDs. The chromosomal blocks discussed above do not include seven short blocks  for which the number of SNPs ranges between two and nine. These were removed from the analysis to reduce the number of false positives. Three of these blocks included one coding gene each, while the remaining blocks did not include any.
We also used the lAmP results to look at chromosomal blocks featuring underrepresentation of hybrid ancestry, using the analogous criterion as for overrepresentation of hybrid ancestry, that is frequency of hybridization-derived variants less than three SD below the mean. We did not find chromosomal blocks meeting this criterion in either wolves or FRDs, possibly because the threshold applied corresponded to strict "ancestry deserts," that is regions with no hybrid ancestry. When we applied a criterion of <0.1% of hybrid ancestry (as in Sankararaman et al., 2014), we identified two blocks with underrepresented dog ancestry in wolves and ten blocks with underrepresented wolf ancestry in FRDs (Table S4).
The error rate estimate based on the "randomised chromosomes" showed that in both wolves and FRDs, one out of the 38 randomly selected blocks was identified as a significant outlier, even though it was not significant in the original analysis ( Figure S5). This gives a false-positive rate of 3%. In FRDs, one chromosomal block that was identified as a significant outlier in the original analysis was not significant in the context of the "randomised chromosome," thus resulting in a false-negative rate of 3%. In wolves, we found no false negatives, but we found one chromosomal block that was identified as a significant outlier in both the original analysis and in the context of the "randomised chromosome" (true positive). The remaining 36 blocks were identified as true negatives in both wolves and FRDs.
The false-positive rate was therefore 3% and the average false-

| Functional characterization of the OHA genes
We carried out a GO enrichment analysis on the OHA gene sets for both wolf and FRD populations. Due to differences in annotation, the resulting set of overrepresented GO terms differed between the canine genes and their human orthologues (Table S5). The GO analysis for the OHA gene set in wolves, based on the canine assembly, showed overrepresentation of the molecular function "ATP-dependent peptidase activity." The analysis based on human homologues revealed overrepresentation of biological processes "amino acid neurotransmitter reuptake" and "glutamate reuptake." All these terms remained overrepresented with both the Benjamini-Hochberg correction and the more strict g:SCS correction.
The OHA gene set in FRDs contained a larger set of enriched GO terms compared to that observed in wolves (Table S5). Terms Other overrepresented GO terms in the OHA gene set for FRDs described more general biological processes involving PCDH genes: "cell-cell adhesion," "cell adhesion," "biological adhesion." Among the molecular functions overrepresented in this gene set were "cation F I G U R E 2 Distribution of dog ancestry in admixed West Eurasian wolves. x-axis shows SNP order along each autosomal chromosomes (without reflecting physical distances between SNP loci), and y-axis shows the proportion of dog admixture in wolves (with only admixed individuals considered). The solid horizontal line represents the mean dog admixture across autosomal chromosomes, and the dotted horizontal line represents the mean dog admixture within each chromosome. Chromosomal blocks with overrepresented dog ancestry are marked in red and are defined as having at least 10 sequential SNPs with the proportion of dog ancestry >3 SD above the mean, which was assessed at the level of individual chromosomes. Ancestry deserts are marked in orange  binding," "metal ion binding" and "calcium ion binding," which are common functions of PCDH genes and several other unlinked genes (Table S5).
The OHA gene set for FRDs also showed overrepresentation of the molecular function "olfactory receptor activity," which was a common function of 17 olfactory receptor (OR) genes identified in an analysis based on human homologues (with eight other OR genes excluded from the GO analysis by the software). These genes were located on three distinct chromosomes. The "olfactory receptor activity" was significantly enriched only when Benjamini- Hochberg correction was used and should therefore be considered with caution. In the GO analysis focused on canine genes, this term did not show significant overrepresentation (with 11 OR genes included), which is likely because the OR gene annotation is less complete in the dog genome assembly and is known to represent many pseudogenes.

F I G U R E 3 Distribution of wolf ancestry in admixed Eurasian free-ranging dogs (FRDs).
x-axis shows SNP order along each autosomal chromosomes (without reflecting physical distances between SNP loci), and y-axis shows the proportion of wolf admixture in FRDs (with only admixed individuals considered). The solid horizontal line represents the mean wolf admixture across autosomal chromosomes, and the dotted horizontal line represents the mean wolf admixture within each chromosome. Chromosomal blocks with overrepresented wolf ancestry are marked in red and are defined as having at least 10 sequential SNPs with the proportion of wolf ancestry >3 SD above the mean, which was assessed at the level of individual chromosomes. Ancestry deserts are marked in orange   (Table S6) and 12 (66.7%) within 100 kb from protein-coding genes. Two (11.1%) were more than 100 kb from any gene, which gives a false discovery rate of 11.1%, similar to the rate estimated in FRDs. We note that this error rate is a conservative estimate, given that functional annotation in the dog genome is less complete compared with the human genome Genes showing signatures of positive selection that are located within chromosomal blocks showing significant overrepresentation of hybridization-derived ancestry, are strong candidates for adaptive introgression and therefore will be referred to as "CAI genes" hereafter. In this study, we did not define the functional variants in the CAI genes and therefore we cannot infer the exact phenotypic effect of introgression of wolf-derived variants into FRDs, or in the reverse direction. However, the GO analysis and the data from published studies regarding the function of the CAI genes provide information on the general type of phenotypic traits that may be subject to adaptive introgression.
In grey wolves, we identified three positively selected SNPs within three genes from the set of 72 OHA genes (Table S6)  In the set of 72 CAI genes in FRDs, the majority of overrepresented GO terms were associated with calcium channel activity, with several genes (CACNA1C, RYR2, RYR3, TRPC4, TRDN, ATP7B, HTR2A) contributing to multiple terms ( Table 2). Most of these genes (CACNA1C, RYR3, TRPC4, HTR2A) are expressed in the brain and associated with neurodevelopment, behavioural traits and cognitive functions (Table 3). In addition, CACNA1C and HTR2Aencoded proteins are involved in viral infections, acting as receptors for influenza virus and JC polyomavirus, respectively (Assetta et al., 2013;Fujioka et al., 2018). RYR2 and TRDN proteins are functionally linked and constitute the main component of a calcium channel that supplies ions to the cardiac muscle enabling its contraction (Györke et al., 2004).
Another set of overrepresented terms was associated with the axoneme, which is the main cytoskeletal structural component of a cilium or flagellum (Ishikawa, 2017), with DNAH5, DNAH17, GAA-CCDC40 and LRGUK genes contributing to these terms (Tables 2   and 3). Formation of ciliary axonemal structures is necessary for regulating motility and beating of the cilia. Therefore, mutations in genes associated with axonemal architecture frequently lead to primary ciliary dyskinesia, characterized by recurrent infections of the respiratory tract and sperm immobility (Lee & Gleeson, 2011;Olbrich et al., 2002).
As explained in the Methods, candidate genes listed above were identified without correction for multiple testing. After applying a Bonferroni correction based on the number of SNPs tested for each chromosome, we found no loci under positive selection within chromosomal blocks with overrepresented hybrid ancestry in wolves and only two such loci in FRDs. These were located within the serotonin receptor subtype 2A (HTR2A) and the type 3 ryanodine receptor (RYR3) genes. HTR2A encodes one of the serotonin receptors, which plays a role in learning and memory (Harvey, 2003).
HTR2A is expressed widely in the brain and regulates levels of several hormones-oxytocin, ACTH, corticosterone, renin and prolactin (Van de Kar et al., 2001). RYR3 protein forms intracellular calcium channels in excitable tissues such as muscles and neurons (Santulli & Marks, 2015). It is expressed in a broad range of tissues, including the brain, where it affects synaptic plasticity (Balschun et al., 1999). The two strongest candidate genes for adaptive introgression in FRDs are thus well-described genes affecting neurobiological processes. The variation in introgression rates may be associated with human density and human footprint index (highest in Europe, lowest in Yakutia). Regions with high human densities may be associated with higher number of free-ranging dogs and thus higher encounter rate between wolves and dogs, which may facilitate admixture.

High frequency of wolves carrying signatures of dog admixture in
Europe could also result from strong and long-lasting hunting pressure, which resulted in the local extirpation of wolf populations in large parts of the continent (Dufresnes et al., 2018). Strong hunting pressure is thought to increase the frequency of hybridization due to the disruption to wolf pack structure (Moura et al., 2014;Rutledge et al., 2012) and thus persistent hunting within a given region over many generations may lead to recurrent hybridization and higher admixture proportions.
Hunting and anthropogenic habitat changes have also led to a reduction in wolf population size. In a small wolf population, a single back-cross event can provide a considerable contribution to the wolf gene pool. Since the size of the domestic dog population is linked to the human population size (Gompper, 2014), the imbalance between wolf and dog population sizes will continue to grow and is likely to lead to further increases in dog introgression into wolves, unless preventive measures are applied (see Donfrancesco et al., 2019;Salvatori et al., 2020).
Based on predictions from ecological studies on the effect of anthropogenic food on large carnivores, attraction to anthropogenic food sources may result in contemporary self-domestication TA B L E 2 Candidate genes under adaptive introgression in grey wolves ("wolves") and free-ranging domestic dogs ("FRDs") that contributed to enriched GO terms Note: N SNPs-the number of SNPs within the gene with signatures of positive selection inferred using rehh. For the complete list of candidate genes under adaptive introgression, see Table S5.
TA B L E 3 Functions of candidate genes under adaptive introgression in free-ranging domestic dogs ("FRDs") and grey wolves ("wolves") associated with the common groups of enriched GO terms (e.g., "calcium channel" denotes all GO terms related to calcium channel activity)

Population Gene GO terms Function
Wolves ABI1 Cell junction Plays an essential role in synapse formation (Proepper et al., 2007).

Wolves PLD1
Cell junction Plays a key role in neurotransmitter release and regulates dendrite morphogenesis (Zhu et al., 2012).
FRDs TRDN Calcium channel Physically links the RYR2 and CASQ2 proteins, enabling the regulation of RYR2 channel activity by CASQ2 (Györke et al., 2004).

RYR2
Calcium channel Ryanodine receptor, which forms intracellular calcium channels in excitable tissues such as muscles and neurons (Santulli & Marks, 2015). Primarily expressed in heart muscle, but affects also neurobiological processes; variation in this gene is responsible for sex differences in autism (Chen et al., 2017).

HTR2A
Calcium channel Encodes one of the serotonin receptors, which is a target of serotonergic psychedelic drugs and antipsychotic drugs (Moreno et al., 2011) that plays a role in learning and memory (Harvey, 2003) and is a receptor for JC polyomavirus (Assetta et al., 2013). Regulates levels of several hormones-oxytocin, ACTH, corticosterone, renin and prolactin (Van de Kar et al., 2001).

FRDs ATP7B
Calcium channel Copper transmembrane transporter (Braiterman et al., 2014) FRDs TRPC4 Calcium channel Forms a calcium-permeable cation channel, which plays a role in multiple processes including neurotransmitter release. Expressed in midbrain dopamine neurons and affects behavioural traits such as attention span and sociability (Illig et al., 2011;Rasmus et al., 2011). Involved in lung endothelial permeability (Tiruppathi et al., 2002).

FRDs CACNA1C
Calcium channel Encodes an L-type calcium channel Ca v 1.2, which is a critical mediator of brain development and experience-dependent brain plasticity. One of the most widely reproduced candidate genes for multiple neuropsychiatric disorders. Affects social behaviour and cognitive function (Kabir et al., 2017). Ca v 1.2 is a receptor for influenza virus (Fujioka et al., 2018).

RYR3
Calcium channel Ryanodine receptor, which forms intracellular calcium channels in excitable tissues such as muscles and neurons (Santulli & Marks, 2015). Expressed in a broad range of tissues, including the brain, where it affects synaptic plasticity (Balschun et al., 1999).
FRDs GAA and CCDC40 Axoneme assembly GAA and CCDC40 genes are located directly next to each other in mammalian genomes, with known insertions-deletions encompassing both genes (Amiñoso et al., 2013). CCDC40 regulates the assembly of the inner dynein arm and the dynein regulatory complexes and thus is essential for correct functioning of motile cilia (Becker-Heck et al., 2011).

FRDs DNAH17
Axoneme assembly Axonemal dynein-a cytoskeletal motor protein that moves along microtubules, driving the beat of eukaryotic cilia and flagella (King, 2018). A sperm-specific dynein associated with reduced sperm motility (Whitfield et al., 2019).
FRDs LRGUK Axoneme assembly Involved in the early stages of axoneme development and in multiple aspects of sperm assembly including sperm head shaping (Liu et al., 2015).

FRDs DNAH5
Axoneme assembly Axonemal dynein-a cytoskeletal motor protein that moves along microtubules, driving the beat of eukaryotic cilia and flagella (King, 2018). Is expressed in lungs and associated with respiratory ciliary disorders (Fliegauf et al., 2005).

| 2447
PILOT eT aL. (Newsome et al., 2017). Hybridization may considerably accelerate this process by enabling rapid acquisition of adaptations facilitating survival in human-modified landscapes (e.g., increased copy number of AMY2B gene facilitating starch digestion- Axelsson et al., 2013).
However, areas with frequent human-wolf interactions tend to have a high density of free-ranging domestic dogs, meaning that the eco-

| Factors affecting introgression patterns in free-ranging dogs
In Eurasian FRDs, the wolf admixture proportion was relatively low Ancient introgression has been reported both from wolves to dogs (Miao et al., 2017;Skoglund et al., 2015) and in the opposite direction (Bergström et al., 2020), but ancient dog-to-wolf introgression has been shown to be considerably more frequent (Bergström et al., 2020). At the early stages of dog domestication, the dog population sizes must have been small, and therefore, backcrossing events, even if rare, must have left a substantial signature in the dog gene pool. This may explain the presence of considerable admixture from a Pleistocene wolf lineage in Arctic and East Asian breeds (Skoglund et al., 2015). In the contemporary FRDs studied here, the average size of introgressed blocks is larger than in wolves, which may suggest that FRDs carry a larger proportion of introgressed blocks originating from recent admixture. However, a number of other factors could have affected the average ancestry block size, including selection on introgressed variants as well as demographic history. Linkage disequilibrium in FRDs is higher than in most wolf populations (e.g., Pilot et al., 2015), and this may affect the size of introgressed blocks.

| The effect of recombination on introgression patterns
In many admixed species and subspecies, local admixture proportions are positively correlated with recombination rates (e.g., Janoušek et al., 2015;Martin et al., 2019;Schumer et al., 2018). This is explained by the presence of variants that are deleterious in hybrids, which constitute barriers to introgression, in multiple loci across the genome (Martin et al., 2019). Recombination determines the extent of linkage between a barrier locus and surrounding neutral loci, and loci close to recombination hotspots are expected to have higher admixture proportions than those in regions of low recombination.
In the present study, we found no significant correlation between admixture proportions and recombination rates in the genome-wide analysis and found a significant positive correlation within only a small number of individual chromosomes. We also found few hybrid ancestry deserts (regions with <0.1% of hybrid ancestry), suggesting that the number of "barrier loci" between wolves and FRDs is small.
This may be the reason for the observed decoupling of introgression rates from recombination rates. Although wolves and dogs are ecologically distinct, their recent divergence, estimated at about 27-40 thousand years ago Skoglund et al., 2015;Wang, Zhai, et al., 2016), overlapping geographic ranges and continuous admixture, can all contribute to the low occurrence of barriers to introgression.
The observed differences in introgression patterns between wolves and dogs could potentially result from differences in the recombination rates between the two canids. To the best of our knowledge, genome-wide recombination maps are only available for domestic dogs, not for grey wolves, and we therefore could not test whether differences in recombination rates between wolves and dogs could affect the estimated size of admixed blocks. However, Muñoz-Fuentes et al. (2015) found no difference in the number and distribution of recombination breakpoints between wolves and dogs within 16 chromosomal regions containing genes associated with phenotypic traits that distinguish dogs from wolves (e.g., coat colour and length, leg length, "dewclaws"). They concluded that recombination patterns were not changed by strong directional selection associated with the domestication process. Therefore, it is unlikely that the shorter average length of admixed blocks in wolves relative to dogs resulted from differences in average recombination rates.
Recombination rate may also affect the selection inference (O'Reilly et al., 2008;Xiang-Yu et al., 2016). Positive selection was shown to reduce population-based estimates of recombination rate (O'Reilly et al., 2008) or, conversely, create false recombination hotspots (Reed & Tishkoff, 2006). Moreover, genome-wide scans to detect signatures of recent selection in humans identified loci located predominantly in regions of low recombination, implying a confounding effect of recombination rate on the power to detect selection (Reed & Tishkoff, 2006). The EHH-based methods reduce this effect by using genomic distances instead of physical distances, and by contrasting the EHH values of two alleles from each locus, thus removing the effect of local recombination rate (Sabeti et al., 2007;Voight et al., 2006;Xiang-Yu et al., 2016). Moreover, in our study recombination rates were not correlated with the local admixture proportions and therefore could not bias the inference of positive selection on introgressed variants.

| Chromosomal blocks with overrepresented introgressed variants
Although the genome-wide introgression rate from FRDs into West Eurasian wolves was over three times higher than that from wolves into Eurasian FRDs (2.7% vs. 0.75%), the number and average length of introgressed chromosomal blocks overrepresented in the gene pool were smaller in wolves than FRDs. Of the genes located within these introgressed blocks, the proportion showing signatures of adaptive introgression was four times smaller in wolves than FRDs (5.6% vs. 23%). This implies that introgression resulting from wolfdog hybridization yields proportionally larger adaptive advantage to FRDs than wolves. This may be due to the larger population size of FRDs compared to wolves, resulting in the increased efficiency of positive selection on weakly advantageous introgressed variants. An increased level of deleterious variation in dogs compared to wolves (Marsden et al., 2016) may also play a role.
In large FRD populations, genetic drift is weak and therefore gene variants with even a weak selective advantage may increase in frequency. Grey wolf population sizes are considerably smaller and adaptive variants are therefore more likely to be eliminated by drift. However, dogs have gone through multiple bottlenecks during their evolutionary history, including the bottleneck associated with the initial domestication process (Lindblad-Toh et al., 2005;Ostrander et al., 2017) as well as founder effects associated with expansion events across the world (e.g., Wang, Zhai, et al., 2016). As a result, the genetic variability of FRDs (measured by heterozygosity as well as the number and size of runs of homozygosity) is smaller than that of most wolf populations, although larger compared to that of pure-bred dogs Marsden et al., 2016). Therefore, FRDs may benefit from introgressive hybridization with wolves by increasing their adaptive variation. Moreover, domestic dogs have a higher genetic load than wolves due to the reduced ability of natural selection to remove weakly deleterious mutations during bottlenecks associated with domestication and breed formation (Marsden et al., 2016). Genetic load is the reduction in average fitness of genotypes in a study population compared to a reference. This reference may be an optimal theoretical genotype, the genotype of an individual with highest fitness in a population or an entire population that has a higher average fitness than the target population. In this case, the wolf population serves as a reference for the dog population. Purebred dogs have on average a two to three per cent greater genetic load than wolves, while the genetic load of FRDs is, as expected, intermediate between wolves and pure-bred dogs (Marsden et al., 2016). Wolf introgression may therefore result in a reduced genetic load in FRDs, providing an adaptive benefit.
It may be expected that introgression of deleterious variation from FRDs into wolf populations will result in purifying selection, thus creating chromosomal blocks where dog ancestry is underrepresented. We did not, however, detect any blocks with underrepresented dog ancestry in wolves when using a threshold of three SD below the mean, possibly because this threshold was too strict.
Using, instead, the criterion of dog ancestry <0.1% (previously used in Sankararaman et al. (2014) to identify Neanderthal ancestry deserts in humans), we identified only two dog ancestry deserts in wolves compared to ten wolf ancestry deserts in FRDs. Although the number of ancestry deserts identified depends on the threshold used to define them, the same threshold applied to wolves and FRDs enables a meaningful comparison of patterns observed in both taxa.
Most of the genetic load in dogs is mediated by weakly deleterious mutations that cannot be easily purged from gene pools (Marsden et al., 2016). It is therefore possible that mutations which are weakly deleterious in FRDs remain weakly deleterious when introgressed into wolves, and can thus be maintained in the wolf gene pool instead of being purged, as wolves have low effective population sizes Pilot, Greco, et al., 2014; see Table S7) and are thus affected by strong genetic drift.
Recessive mutations of large effect are known to occur in purebred dogs, but they are likely subject to strong purifying selection in FRD populations, and are thus unlikely to be transferred to wolves, as unsupervised hybridization between pure-bred dogs and wolves is rare. Our result suggests that only a small proportion of variation derived from FRDs has a strongly deleterious effect in wolves. The larger number of wolf ancestry deserts detected in FRDs may results from larger efficiency of selection against deleterious variants in FRDs, due to the larger population size compared to wolves.
Our finding that higher genome-wide admixture proportions are not necessarily associated with a larger number of blocks with overrepresented hybrid ancestry is consistent with the study by

| Genes within the chromosomal blocks with overrepresented introgressed variants in wolves
In the West Eurasian wolves, the set of 72 OHA genes (i.e., genes located within the chromosomal blocks characterized by overrepresented dog ancestry) was enriched for the GO terms "amino acid neurotransmitter reuptake" and "glutamate reuptake." Glutamate is the most abundant excitatory neurotransmitter in the vertebrate central nervous system and accounts for over 90% of synaptic connections in the brain (Platt, 2007). Glutamate receptors play a key role in the induction and maintenance of synaptic plasticity and are associated with learning and memory (Peng et al., 2011). They regulate genomic responses to dopamine stimulation in the neurons of the striatum, a part of the forebrain that coordinates motor functions as well as multiple cognitive functions such as action planning, motivation and reward perception (Wang et al., 2003).
Of 72 OHA genes we identified in West Eurasian wolves, only four (5.6%) showed evidence for adaptive introgression (Table 1), that is were characterized by overrepresentation of dog ancestry and signatures of positive selection. All four genes are involved in neurotransmission and neurodevelopment, and in humans, mutations in these genes are associated with disorders such as schizophrenia, Alzheimer's disease and congenital microcephaly (Table 3).
Overall, these results suggest that certain dog-derived behavioural or cognitive traits may be advantageous to wild wolf populations, highlighting the need for studies on behavioural consequences of wolf-dog hybridization based on direct observations.

| Genes within the chromosomal blocks with overrepresented introgressed variants in freeranging dogs
The set of 311 OHA genes from the FRDs included 35 OR genes located on three different chromosomes, and 13 protocadherin (PCDH) genes clustered on chromosome 2. The presence of multiple genes from the same functional groups is reflected in the results of the GO analysis. OR activity was a significantly enriched molecular function and "homophilic cell adhesion via plasma membrane adhesion molecules," a type of calcium-dependent cell adhesion specific to protocadherins, was a significantly enriched biological process.
The efficient function of both OR genes and PCDH genes is dependent on genetic variation within the gene families (Chen & Maniatis, 2013;Trimmer et al., 2019). Introgression from wolves is likely to increase that variability, thus alleviating the loss of genetic diversity in dogs resulting from the domestication bottleneck. Improvement in FRD olfactory abilities via introgression from wolves may facilitate detection of food sources, identification of unsuitable food, and detection of potential threats (e.g., humans or large predators).
Clustered PCDHs are involved in calcium-mediated transcriptional gene networks, are expressed primarily in the developing nervous system and play a key role in many neurodevelopmental processes, including axon guidance, creation of new synapses and dendritic self-avoidance (Garrett & Weiner, 2009;Lefebvre et al., 2012). The organization of the PCDH gene cluster enables the expression of multiple gene isoforms, facilitating the diversification of surface molecules in neuronal cells (Chen & Maniatis, 2013).
Therefore, clustered PCDHs are considered as "molecular barcodes for self-recognition by individual neurons in the vertebrate nervous system" (Chen & Maniatis, 2013). It is thus likely that adaptive benefits of wolf introgression in FRDs may result from both introduction of new adaptive variants and an increase in overall diversity of PCDH genes. PCDHs are under balancing selection in humans (Noonan et al., 2003) and this may be also the case in other vertebrates, which may explain why we found no signatures of positive selection in these genes.
Several studies regarding the genetic basis of animal domestication show that protocadherins displayed differential expression and allele frequency changes between domesticated and wild populations, as well as between populations that were experimentally selected for tame versus aggressive behaviour .
Comparative genomic analyses revealed diversifying selection on several PCDH genes between domestic and wild cats (Montague et al., 2014) as well as between tame and aggressive strains of silver foxes (Vulpes vulpes) originating from the Farm Fox Experiment . Genes from PCDHGA subfamily displayed differential expression in the brain between strains selected for either tame or aggressive behaviour in both silver foxes and rats (Rattus norvegicus) (Heyne et al., 2014;Wang et al., 2018). Tameness in animals is characterized by reduced fear of humans and lowered human-directed aggression and is, therefore, a necessary trait for wild animals living in close proximity of humans. Although FRDs gain various ecological benefits from living close to human populations, humans are also the main source of early-life mortality in FRDs (Paul et al., 2016). Extreme tameness associated with a complete lack of fear of humans may therefore lead to reduced fitness. Introgression from wolves could help ensure that tameness in FRD populations does not reach extreme levels where it could become maladaptive.
Interestingly, in wolves we also found overrepresentation of dog variants in a protocadherin gene. This gene, a nonclustered protocadherin PCDH15, plays a key role in the formation of sensory hair cells as well as the retina, and variants within this gene are associated with vision and hearing impairment in humans (Jacobson et al., 2008;Kazmierczak et al., 2007). Since we did not identify sequence level variants, we do not know the nature of the resulting functional change in canids. It will be interesting to determine such in subsequent studies, thus providing another example of how studies of canine genome evolution can inform disease studies in humans (see Ostrander, 2012).
Of the 311 OHA genes observed in FRDs, 72 (23%) showed signatures of adaptive introgression ( Table 1). The majority of resulting enriched GO terms were associated with calcium channel activity ( Table 2). Dysregulation of calcium-mediated transcriptional gene networks can disrupt the development of neuronal circuits (Kabir et al., 2017;Lohmann, 2009;Ramocki & Zoghbi, 2008). Multiple genes affecting calcium channel activity have been identified in earlier studies as candidate genes for neuropsychiatric disorders in humans and are shown to affect behavioural traits (Table 3 and references therein). Genes affecting behaviour were likely a selection target during the domestication process and introgression of wolfderived variants may introduce some wolf-like behavioural traits, which may be adaptive in free-living dog populations.
One of the candidate genes under adaptive introgression in FRDs, TRPC4, encoding a nonselective calcium-permeable cation channel, was found to be highly differentiated in Arctic sled dogs relative to other dogs and is therefore hypothesized to play a role in cold-climate adaptation, together with several other genes involved in calcium ion transport . Since our study includes dogs from different climate zones across Eurasia, the adaptive introgression observed on this geographic scale cannot be associated with climate adaptation. This suggests that TRPC4, as a gene with pleiotropic effects, could have contributed to both global and local adaptations in dogs.
Among the set of 72 genes showing adaptive introgression signatures, we also found overrepresentation of terms associated with the axoneme, the main structural component of a cilium. Genes associated with axoneme assembly are frequently involved in dyskinesia of primary cilia, which are specialized organelles responsible for calcium signalling within cells that regulate the hedgehog signalling pathways (Delling et al., 2013;Mukhopadhyay & Rohatgi, 2014).
Several genes involved in the functioning of primary cilia were previously shown to be under diversifying selection between FRDs and pure-bred dogs, suggesting their role in the independent survival of free-living populations (Pilot et al., 2016). This could be related to the role of primary cilia in reproduction (Lee & Gleeson, 2011;Olbrich et al., 2002), and/or their involvement in calcium signalling that affects neurodevelopment. Thus, introgression from wolves may result in the improvement of reproductive fitness of FRDs, or in the introduction of behavioural or morphological traits that may facilitate their independent survival.
The candidate genes discussed above were identified using iHS statistics (Voight et al., 2006) without correction for multiple testing. Two remaining candidate genes under adaptive introgression in FRDs after Bonferroni correction are HTR2A and RYR3, which are involved in calcium signalling, affect neurobiological processes and are associated with behavioural disorders in humans and other mammals (Matsuo et al., 2009;Serretti et al., 2007; Table 3). Therefore, even after applying a very conservative approach, we can maintain our conclusion that introgression from wolves has a particularly strong effect on behavioural traits in FRDs.

| CON CLUS I ON S AND IMPLI C ATI ON S FOR WOLF P OPUL ATION MANAG EMENT
Our results imply that introgressive hybridization with wolves is beneficial for FRD populations, given the evidence of adaptive introgression in multiple genes with important functions. In grey wolves, the introgression is driven mostly by drift, with a small number of positively selected genes being associated with brain function and behaviour. The predominance of drift may be a consequence of the small effective population sizes, resulting in the reduced efficiency of selection on weakly advantageous or against weakly disadvantageous introgressed variants. Variants associated with strong selective disadvantage are likely to be eliminated by natural and sexual selection so that they do not spread from F1 hybrids or recent- Finally, our results suggest that genes affecting neurobiological processes predominate among genes displaying higher than average introgression rates in both wolves and FRDs. Therefore, introgressive hybridization can affect behavioural or cognitive traits in both canids. Changes in behaviour of wild wolves can have important ecological consequences, and therefore, monitoring of wolf populations affected by hybridization should include behavioural observations.

ACK N OWLED G EM ENTS
We are grateful to Dr Kristien I. Brans, Dr Marty Kardos and two anonymous reviewers for their comments that allowed us to greatly improve the paper.

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

DATA AVA I L A B I L I T Y S TAT E M E N T
Data used in this study are available at the following repositories: