Combining phylogenetic and demographic inferences to assess the origin of the genetic diversity in an isolated wolf population

The survival of isolated small populations is threatened by both demographic and genetic factors. Large carnivores declined for centuries in most of Europe due to habitat changes, overhunting of their natural prey and direct persecution. However, the current rewilding trends are driving many carnivore populations to expand again, possibly reverting the erosion of their genetic diversity. In this study we reassessed the extent and origin of the genetic variation of the Italian wolf population, which is expanding after centuries of decline and isolation. We genotyped wolves from Italy and other nine populations at four mtDNA regions (control-region, ATP6, COIII and ND4) and 39 autosomal microsatellites. Results of phylogenetic analyses and assignment procedures confirmed in the Italian wolves a second private mtDNA haplotype, which belongs to a haplogroup distributed mostly in southern Europe. Coalescent analyses showed that the unique mtDNA haplotypes in the Italian wolves likely originated during the late Pleistocene. ABC simulations concordantly showed that the extant wolf populations in Italy and in south-western Europe started to be isolated and declined right after the last glacial maximum. Thus, the standing genetic variation in the Italian wolves principally results from the historical isolation south of the Alps.


Introduction
Human activities have deeply shaped the structure of landscapes in continental Europe for millennia, often reducing the extension of pristine ecosystems to small isolated fragments [1]. Nowadays habitat fragmentation is a major threat to the survival of natural animal populations [2]. Carnivores and ungulates in particular need wide extensions of suitable habitat and are among the most sensitive species to habitat fragmentation. Ungulates have been intensively hunted for meat while large carnivores were persecuted as pest predators, and both declined for centuries simultaneously with the fragmentation and changes of their natural habitats [3]. However, the ongoing rewilding wave in continental Europe demonstrates that declining population trends can be interrupted and even reversed: endangered top predators and ungulates are now locally abundant, and food chains in forest ecosystems are partially restored [4]. Therefore, it is interesting to understand the consequences of past declines and fragmentation in the perspective of population expansion and forthcoming rejoining.
The genetic consequences of protracted population declines are theoretically well known but not easily predictable in empirical case-studies [5,6]. According to the strength of bottlenecks and time of isolation, genetic diversity (number of alleles, heterozygosity) is progressively lost by random drift. Loss of selectively neutral genetic diversity is correlated to the effective size of the population (Ne), which is usually much lower than the observed census population size (Nc) [7]. Estimating Nc in natural populations is not simple, and estimating Ne is even more difficult due to the variability of its several determinants: reproductive success, sex-ratio, departures from random mating and variations of these factors from one generation to the next. Various simulation methods, including the Approximate Bayesian Computation (ABC) [8], have been implemented to evaluate the role of demographic parameters in determining the dynamics of genetic diversity in complex historical scenarios [9]. However, if Ne is not too small, the outcomes of genetic drift can be contrasted by balancing or frequencydependent selection and by hitchhiking on functional gene complexes [10]. Although the fitness consequences of small Ne and low standing genetic variation are controversial, the assumption that adaptive potential and evolvability are positively correlated to heterozygosity is a precautionary principle widely accepted in conservation genetics [11].
In this study, we reassessed the amount and origin of genetic variation in a wolf (Canis lupus) population that remained isolated and declined in peninsular Italy for centuries [12,13]. Wolves disappeared from the Alps in the early 1900s and were strongly reduced in the peninsular regions until the early 1970s, when less than 100 individuals survived in two isolated subpopulations in remote mountain areas of the central-southern Apennines [14]. Both legal protection (granted to the wolf in the 1970s) and deep socio-ecological changes (industrialization, urbanization, and the abandonment of marginal agricultural lands in mountains and hills) boosted spectacular re-expansions of forests, wild ungulates and consequently wolves [4]. Permanent wolf packs rapidly established along the entire Apennines chain, reaching the western Alps in the early 1990s, then colonizing new areas in lowlands and also in southern Italy. Currently, the Italian wolf population consists of c. 1500 (± 300) individuals and at least 320 documented packs [15].
Regardless of the strength of the ongoing wolf expansion trend, the amount of standing genetic variability should have been largely determined by the duration of the last bottleneck, because there has not been enough time for mutations to reconstitute the lost variation [13]. However, the effective isolation period of the Italian wolf population south of the Alps is controversial. Genetic and genomic data suggest that it could have been effectively isolated for 9500-19 000 years (as estimated using microsatellite data, depending on inferred Ne values) [13], 3200-5600 years (as estimated from genome-wide SNP data) [16], or 2800-7000 years (from whole genome sequences, depending on the mutation rate applied [17]), much longer than the c. 100 years estimated from the observed isolation south of the Alps. Genetic analyses indicate that the autosomal genetic diversity in the Italian wolves is c. 30% lower than in other wolf populations in Europe [16]. The mtDNA genetic diversity is limited to a single controlregion (CR) haplotype named W14 that is widespread uniquely in the Italian wolves [18], although a recent study suggests that a second rare haplotype named W16 is also naturally present in the population [19]. Monitoring programs and molecular analyses revealed an increasing occurrence of wolf hybridization with free-ranging domestic dogs in sectors of the wolf range in Peninsular Italy [20,21]. The uncertain length of isolation and the occurrence of hybridization make it difficult to understand if the observed genetic diversity is explained by the last bottleneck or by admixtures.
In this study we genotyped wolves sampled from the Italian and nine other populations in Europe at four mtDNA regions (the control-region CR, and three coding genes: ATP6, COIII and ND4) and 39 autosomal microsatellites [20]. Aiming to identify possibly-admixed genotypes and clarify the origin of the observed genetic diversity, we included also 69 village dogs and 74 known wolf x dog hybrids sampled from areas within the wolf distribution range in Italy and Estonia. Specifically, we aimed to: 1) confirm the attribution of the rare mtDNA CR haplotype W16 to the Italian wolf population [19], increasing the sample size of wolves and hybrid canids from different European population, and 2) test if the standing genetic variation in the Italian wolves is mainly determined by a historical isolation south of the Alps dating back to the end of the last Pleistocene glaciation, or to the most recent anthropogenic bottleneck about one century ago. In order to validate the most likely hypotheses, we performed Bayesian assignment analyses to investigate the partition of autosomal genetic diversity among wolf populations, then we applied Bayesian phylogenetic procedures to estimate the divergence times among wolf and dog mtDNA sequences and ran ABC simulations to reconstruct past demographic scenarios and infer splitting dates of extant wolf populations in Europe.

Ethics statement
No animals were sacrificed for the only purposes of this study. In Italy, all samples of founddead wolves were collected by specialised technician personnel on behalf of the Italian Ministry of Environment (MATTM) and Italian Institute for Environmental Protection and Research (ISPRA). Wolf stool samples from Czech Republic and Slovakia were collected by Friends of the Earth organisation (FoE CZ), which is monitoring the wolf population in the Carpathian Mountains. In Slovakia, FoE CZ has permission to collect non-invasive samples of wolves, issued by Regional Office Trenčín, Department of Environment, No. OU-TN-OSZP1-2014/49/3475. Carpathian wolf tissue samples were legally culled during the open hunting season in Slovakia within a quota set by the local authorities, in conformity with regulation No 344/2009 Coll. The wolves were shot during individual patrols or collective hunts. The use of poisoned bait or leg-hold traps is strictly forbidden according to hunting law. All hunters had permission for hunting, and we confirmed that the culls were reported before quota fulfilment. Croatian and Slovenian samples were obtained either from animals killed in traffic accidents or culled during regular hunting management according to quotas defined by the Croatian Commission for monitoring large carnivore populations and approved by the Croatian Ministry for Environmental and Nature Protection. Iberian wolf samples from Portugal were from found-dead individuals collected by specialised technician personnel on behalf of the Portuguese Institute for Nature Conservation (ICNF). Wolf samples from Spain were obtained from animals that were road-killed or culled during regular hunting management according to quotas defined by Principado de Asturias authorities. Samples from Estonia, Latvia, Finland, Poland, Greece, and Bulgaria were also collected from animals found dead or legally harvested by hunters for purposes other than this project. No ethics permit was required since the sample collection involved dead animals. All samples were collected by specialised technician personnel. No ethic permit is also required to collect stool samples in these countries. Dog blood samples were obtained by veterinaries with the assistance of the owners and all the possible efforts to minimize animal suffering. The owners of the dogs gave permission for their animals to be used in this study. Salivary samples were obtained through buccal swabs by specialised technicians.

Sampling
We genotyped 190 wild-living wolves, 69 village dogs and 74 known wolf x dog hybrids from Italy and several other countries in Europe (Table 1). Moreover, we analysed other five unrelated wild-living wolves, sampled in different areas of the central-northern Apennines, that had the rare mtDNA haplotype W16 [19,20] (Table 1). Wolf samples were collected from 1990 to 2015 in Italy (sample size n = 34), Spain and Portugal (Iberian Peninsula; n = 20), Slovenia (20), Croatia (20), Greece (15), Bulgaria (17), Czech Republic and Slovakia (20), Poland (16), Estonia (10), Latvia (10), and Finland (9). The 69 village dogs were sampled in Italy from northern and central Apennine areas, had size and shapes similar to shepherd dogs, did not belong to certified breeds, and were selected independently of genotypic information. The known wolf x dog hybrids were sampled in Italy (68) and Estonia (6), and were previously identified by genetic or morphological analyses [20,22]. All wolves had the typical wolf coat colour pattern and did not show any apparent morphological sign of hybridization. We stored tissue and blood samples in 10 volumes of 95% ethanol or five volumes of a Tris/SDS buffer [23] at -20˚C, respectively. We extracted DNA samples using the QIAGEN DNeasy tissue extraction kit (Qiagen Inc, Hilden, Germany) in a robotic liquid handling system MULTIP-ROBE II EX (Perkin-Elmer, Weiterstadt, Germany). All the DNA samples used in previous studies were reanalyzed. Negative extraction controls (no DNA in the test tubes) were used to check for laboratory contaminations.

Microsatellite genotyping
We selected a panel of 39 canine autosomal microsatellites (seven tetranucleotides and 32 dinucleotides) that were used in some of the most recent studies on wolf population genetics and hybridization in Europe ( [24] and references therein). These microsatellites mapped on 26 different chromosomes (S1 Table) and were not in linkage disequilibrium in the studied populations. The panel includes 15 markers from the Finnzymes Canine Genotypes™ Panel 1.1 multiplex kit (Finnzymes, Thermo Fisher Scientific, Vantaa, Finland). One of them, the Amelogenin, was used to sex the individuals. The microsatellites were amplified in eight PCR multiplexes using the Qiagen Multiplex PCR Kit (Qiagen, GmbH-Hilden, Germany). Negative (no DNA) and positive (samples with known genotypes) PCR controls were used to check for laboratory contaminations. To confirm allele calls, all samples were independently analysed twice to check for the occurrence of allelic dropout and false alleles, which were never observed. The amplicons were analysed in an ABI 3130XL automated sequencer (Applied Biosystems; Foster City, California, USA) and allele sizes were estimated using the software GENEMAPPER 4.0. Details on the selected markers, primers and PCR profiles are available in the S1A Appendix.

Microsatellite variability, genotype clustering and assignment testing
Observed and effective allele numbers (Ao and Ae), observed and expected heterozygosity (Ho and He), F-statistics [25], and tests for departures from Hardy-Weinberg (HWE) and linkage equilibria (LE) were computed in GENALEx 6.01 [26] and GENETIX 4.05 [27]. The individual multilocus genotypes were clustered and assigned to their most likely population of origin using: 1) a principal coordinate analysis (PCoA in GENALEx); 2) a discriminant principal component analysis (DAPC in ADEGENET) [28]; and 3) a Bayesian clustering model (minimizing departures from HWE and LE in the genetic clusters) implemented in STRUCTURE 2.3.4 [29,30]. We used STRUCTURE to infer the optimal genetic partition of the sampled groups, assuming K from 1 to 15 with four independent runs for each K with 500 000 Monte Carlo Markov Chains (MCMC) steps and discarding the first 50 000 steps as burn-in, using the admixture and independent allele frequency models, and no prior information (option usepopinfo not activated). We used CLUMPAK (http://clumpak.tau.ac.it) to identify the highest rate of increase in the posterior probability LnP(D) of the clusters between each consecutive K [31] and to aggregate the individual membership probability (qi) from the four MCMC replicates [32].
All the amplifications were performed in 10 μL total reactions containing 20-40 ng/μL DNA, 1X PCR buffer with 2.5 mM Mg 2+ , 0.3 μM of primer mix (forward and reverse) and 0.25 units of Taq Polymerase (5 PRIME Inc., Gaithersburg, USA). Amplifications were performed with an initial DNA denaturation step at 94˚C for 2 minutes, followed by 45 cycles of denaturation at 94˚C for 15 seconds, annealing at 55˚C for 15 seconds, extension at 72˚C for 30 seconds, and final extension at 72˚C for 10 minutes. Amplicons were purified using ExoSAP-IT (Affimetrix, Inc., Cleveland, Ohio, USA) and sequenced in both directions in an ABI automated DNA sequencer 3130XL (Applied Biosystems). Sequences were visually corrected in SEQSCAPE 2.5 and aligned in GENEIOUS 7.1 (Biomatters Ltd., Auckland, New Zealand). GENEIOUS was also used to fix alignment ambiguities, mainly caused by indels in the mtDNA CR. The four mtDNA regions were concatenated in a multi-fragment alignment of 2164 bp. Taking into account the presence of indels, identical haplotypes were collapsed using DNASP 5.10.01 [35], that was also used to estimate haplotype (H) and nucleotide (π) diversity [36] for each of the four regions and for the concatenated sequences.

Phylogenetic analyses
In addition to the new 210 mtDNA sequences (S2 Table), we downloaded from GenBank the homologous sequences of 18 extant wolves, four ancient canids and 322 dogs (S3 Table). We then aligned these sequences to construct Neighbor-Joining (NJ) [37], Maximum Likelihood (ML with heuristic search) [38] and Bayesian (BT) [39] phylogenetic trees that were rooted using as an outgroup a coyote sequence (Canis latrans, GenBank access number DQ480509). NJ and ML phylogenetic analyses were done in PAUP Ã 4.0 beta [40]; the BT was computed in MRBAYES 3.2 [41]. For the NJ and ML analyses, we selected the best evolutionary model using the modeltest option [42] and the Akaike Information Criterion [43]. We obtained internode supports by 1000 bootstrap replicates [44] in NJ trees, and by 100 bootstrap replicates in ML trees, using the faststep search in PAUP. We identified the best-fit evolutionary model for Bayesian analyses using PARTITIONFINDER [45]. MRBAYES 3.2 was run for 2x10 6 generations, with a sampling frequency of 100 generations, and with one cold and three heated MCMC (temperature = 0.45; first 10% of the trees excluded to ensure convergence) [46]. To check for convergence of parameter for Bayesian analyses, we used Tracer 1.6 (http://tree.bio.ed.ac.uk/ software/tracer). We estimated in GENEIOUS the frequency of invariable (I) sites, the parameters of the gamma (γ) distributions and the transition-transversion (Ti/Tv) ratios for the NJ, ML, and BT analyses.

Estimates of mtDNA haplotype divergence times
We used BEAST 2.4.2 [47] to estimate the coalescent time to the most recent common ancestor (TMRCA) of the main mtDNA clades identified in the phylogenetic trees of the concatenated sequences, including all the wolf and dog haplotypes. We applied a strict molecular clock with the coalescent extended bayesian skyline model [48], and using as priors the ages of the four ancient canid sequences (S3 Table; S1B Appendix). Model parameters and trees were sampled every 10 000 over a total of 100 000 000 iterations in two independent MCMC chains. The first 10% iterations were discarded as burn-in. We used TRACER to check for MCMC convergence. When the two independent runs converged on the posterior distributions and reached stationarity, we combined the sampled trees into a single tree file with LOGCOMBINER 2.

ABC simulations
We used microsatellite data to run Approximate Bayesian Computation simulations (ABC) [49] implemented in the software DIYABC 2.1.0 [50] to model plausible demographic scenarios and estimate divergence times (in generations) among wolf populations sampled from European countries and corresponding to the clusters identified by STRUCTURE. We selected three wolf population samples for modelling full ABC simulation scenarios: WIT (pop1), WIBP (pop2) and WDIN (pop3), excluding any sample with possible traces of dog admixture (S2 Table). Samples from WBALK, WCARP and WBALT were not used in the simulations because these populations are still in connection one another and also with unsampled wolf populations in eastern Europe [51,52].
According to Pilot et al. [16] and Fan et al. [53], southern European populations diverged very closely in time, and their effective sizes steadily decreased in the last tens of thousands of years. Therefore, we tested four demographic scenarios (S5 Fig), assuming that the three populations split simultaneously (scenarios 1 and 2) or sequentially (scenarios 3 and 4) and that the three populations passed through a bottleneck (scenarios 2 and 4) or not (scenarios 1 and 3).
We ran 6 x 10 6 simulations for each scenario using uniform prior distributions of the effective population size and time parameters with default mutation settings. We selected the following summary statistics for all the microsatellites: a) one sample: mean number of alleles, mean genetic diversity, mean size variance; b) two samples: mean number of alleles, mean genetic diversity, Fst, shared allele distance (S4 Table).
Scenarios were compared by estimating posterior probabilities with the logistic regression method in DIYABC using 1% of the simulated datasets. For the best models, posterior distributions of the parameters were estimated with a logit-transformed linear regression on the 1% simulated datasets closest to the observed data. Scenario confidence was evaluated by comparing observed and simulated summary statistics. Finally, the goodness-of-fit of the posterior parameters for the best performing scenarios was tested via the model checking option with default settings, and significance was assessed after Bonferroni correction for multiple testing [54].

Microsatellite variability and cluster analysis
All the 39 microsatellites were polymorphic in the sampled groups, confirming previously published results [20]. The PCoA plotting of the multilocus genotypes showed that the Italian wolves and dogs are sharply distinct from one another and from all the other wolf populations (Fig 1A). The DAPC plot confirmed this pattern of population clustering (not shown). The first discriminant function in DAPC further showed that wolves from the Iberian Peninsula clustered in a separate sector from the other regions in Europe (Fig 1B). Results of the multivariate analyses are supported by STRUCTURE results, where the best clustering were obtained at K = 3 and at K = 7, wherein the optimal subdivision for the European wolf populations was observed (Fig 2; S1 Fig). Dogs and wolves cluster separately at K = 2, and the Italian wolves were the first to cluster separately from the other wolves at K = 3, followed by the Iberian wolves at K = 4, then by wolves from Slovenia and Croatia at K = 5, from Greece and Bulgaria at K = 6, from Carpathian region, Estonia, Latvia and Finland at K = 7 (S2 Fig). We assume that these seven clusters represent the main genetic subdivision among the sampled populations, and hereafter we will refer to the following wolf populations: Italian Peninsula (WIT), Iberian Peninsula (Spain and Portugal; WIB), Dinaric regions (Slovenia and Croatia; WDIN), Balkanic regions (Greece and Bulgaria; WBALK), Carpathians (Czech Republic; WCARP), and Baltic Countries (Estonia and Latvia plus Finland; WBALT) ( Table 1). WIT and WIB did not show signatures of admixed ancestry, which, however, was apparent in some individuals sampled in other clusters (Fig 2). The five wolves carrying the W16 CR haplotype were totally assigned to the WIT cluster with q wit > 0.99 and 90% CI = 0.95-1.00 (samples W943, W1223, H1122, W1816 and W1906 listed at the end of S2 Table) and were confirmed to belong to the WIT population.
The numbers of observed and effective microsatellite alleles were the lowest in WIT (Ao = 3.9; Ae = 2.  (Table 2).   Table). The Italian wolves (WIT) showed the lowest number of haplotypes (N), and the lowest haplotype and nucleotide diversity (H and π) at each of the four mtDNA regions ( Table 2). Nucleotide diversity was lower than 1.4% in all the sampled populations, suggesting recent origins of their mtDNA diversity. The concatenated mtDNA sequences yielded two distinct haplotypes in the Italian wolves, five in WIB, three in WDIN, 10 in WBALK and six in WBALT ( Table 2). The two Italian wolf haplotypes, WH14 and WH19, differed by a single nucleotide substitution at position 15 629 of the complete Canis lupus mitochondrial genome (GenBank access no AB499825). Haplotype WH14 includes the shorter mtDNA CR haplotype W14 repeatedly described in the Italian wolf population [13,18,20,52]. The concatenated haplotype WH19 includes the shorter mtDNA CR haplotype W16 already described by Boggiano et al. [55] and Randi et al. [20] (S2 and S5 Tables) in five of our analyzed samples. We did not detect WH19 in any other individual (S6 Table). All the 40 sequenced wolf x dog hybrids sampled in Italy showed the same wolf haplotype WH14, suggesting a preferential wolf maternal-biased hybridization [22]. The BT, NJ and ML phylogenetic trees of the concatenated mtDNA sequences showed very similar topologies, although node supports and relationships among the less divergent haplotypes differed (S3A, S3B and S4 Figs). The two closely related Italian wolf haplotypes WH14 and WH19 belong to the same basal clade named A1 in the phylogenetic trees. Clade A1 is strongly supported in the NJ (bootstrap = 84.7), ML (bootstrap = 70) and BT (posterior probability = 1.0) trees, and includes other five haplotypes (named WH15, WH17, WH18, WH20, and WH21) that were sampled in wolves from Slovenia, Greece, Bulgaria and Poland. All these haplotype share the same ATP6 (A3) and COIII (C3) haplotypes, but WH14, WH17 and WH19 are the only haplotypes that carry the ND4 haplotype N5, whereas the other four haplotypes carry the haplotype N4 (S5 Table). None of these haplotypes were found in any other extant wolf population analysed so far. Clade A1 is the sister clade of the strongly supported clade A2 that includes five haplotypes found only in ancient dog breeds (Swedish and Norwegian elkhounds; S3 Table) and the haplotype S14.5k that was identified in an ancient wolf sample from Switzerland dated at c. 14 500 years ago [56]. Clade A is a sister clade to all the other modern wolf and dog clades in the BT, suggesting that it includes ancient mtDNAs currently distributed mainly in southern European wolf populations. Clades A1-A2 are closely related to, but not nested within, haplogroup A2 as defined by Pilot et al. [57]. With a few exceptions, the other clades are composed exclusively of dog or wolf haplotypes [58,59].

Bayesian estimates of the mtDNA TMRCA
The Bayesian tree generated by BEAST (Fig 3) was similar to the other trees. All the posterior probabilities of the main internodes were > 0.81, except for node 2 (P = 0.54) and node 6 (P = 0.53). Consistently with Thalmann et al. [56] the ancient wolf-like specimens from Belgium (B30k) and Russia (R18k) were the most basal haplotypes in the tree (node 1 and 2. The extant wolf haplotypes split into the two clades A and B 45 400 years ago (node 3). Haplotypes within clade B coalesced 41 100 years ago (node 4), whereas clade A1 and A2 coalesced 28 300 years ago (node 5). Clade A2 includes the historical sequence S14.5k [56] and the dog haplotypes Dog1, that diverged 26 000 years ago, whereas the haplotypes included in clade A1 coalesced 6800 years ago (node 7).

Results of the ABC simulations
ABC simulations provided the best support for scenario 2 (simultaneous population splitting with bottlenecks) that clearly better performed than the other three (S6 Fig). The best scenario showed non-significant P-values for all the posterior parameters after Bonferroni correction (S7 Table). Under this scenario the median values of the divergence time showed that the three wolf populations have been genetically isolated for the last 6830 generations (5% quantile (q050) = 3240 generations-95% quantile (q950) = 9600 generations) ( Table 3). Assuming a wolf generation time of three years [60,61], the TMRCA of these populations is 20 490 years ago, while their bottlenecks were estimated around 15 030 years ago (S7 Fig; Table 3). The current effective populations sizes (N1, N2 and N3; S5 Fig) were much lower than the corresponding effective populations sizes before the demographic decline (Table 4, Fig 4). In particular, the bottleneck led the Italian wolf population to decline by c. 1.9 times while the Iberian and Dinaric populations c. 4.4-3.0 times, respectively.

Discussion
During the last few centuries large carnivores and ungulates declined in south-western Europe as consequence of deforestation, overhunting and direct persecution [4,62]. In peninsular Italy wolves survived in isolation after protracted range contraction and strong demographic decline. In the first half of the 1970s, during the most recent bottleneck, only c. 100 individuals remained in isolated mountain areas in central-southern Italy [14]. Several sources of genetic data concordantly highlighted the consequences of the population bottleneck. Microsatellite markers and genome-wide SNP screenings showed that the Italian wolves have c. 33%-42% less autosomal genetic variability than all the other wolf populations in Europe [13,16,52]. An exception is the isolated Iberian wolf population, which also shows low genetic variability [16,63]. Furthermore, uniparental genetic diversity is also very low in the Italian wolves, and until recently, only one mtDNA CR and two Y-chromosome haplotypes were described [20,64]. Inferences from those estimates of genetic variability indicated that wolves might have been isolated in peninsular Italy for thousands of generations, thus ruling out a predominant effect of the most recent anthropogenic bottleneck [13,53]. However, the time of the effective genetic isolation south of the Alps is controversial, with estimates ranging from c. 3000 to 19 000 years depending of the molecular markers (microsatellites, SNPs or entire genomes), or on the assumptions about the effective population size and mutation rates [13,17,53]. The  genetic consequences of the population bottlenecks and the origin of the standing genetic variation in the Italian wolves are still not satisfactorily understood. The results of our study confirm that the mtDNA CR haplotype WH19 is not originated via introgressive hybridization with dogs, but instead it was likely already present in the Italian propulation before the 20 th century bottleneck [19]. In Bayesian admixture analyses, the autosomal multilocus genotypes of the five samples with haplotype WH19 were entirely assigned to the Italian wolf population (S2 Table). Moreover, we did not find this haplotype in any other wolf population worldwide, nor in dogs. Haplotype WH19 is identical to WH14, except for one SNP at the control region. Closely related to these two haplotypes is haplotype WH17, found in two wolves from Greece. These three haplotypes belong to the same monophyletic mtDNA clade A1 (Fig 3), which includes four other related haplotypes identified in wolves sampled mainly in southern Europe and in the Balkans: Peninsular Italy (haplotypes WH14 and WH19), Greece (WH15, WH17 and WH18), Bulgaria (WH18), Slovenia (WH20), and Poland (WH21), consistent with Pilot et al. [57]. This phylogeographic pattern could represent ancient relic lineages, as indicated by the basal position of clade A in the phylogenetic trees (Fig 3 and S4 Fig). The origin of clade A1 seems extremely recent, dating back to 6800 years ago (95% HDP = 2600-12 400 years), when it coalesced with its sister clade A2. This latter includes the haplotype S14.5k, which was identified in an ancient wolf sample from Switzerland dated at c. 14 500 years ago [56], and five haplotypes found only in ancient Scandinavian dog breeds (S3 Table; Fig 3). Clade B can also be divided in two sub-clades, with clade B1 including the haplotypes found in clade 1 from Pilot et al. [57], and clade B2 containing haplotypes that are intermediate between clades 1 and 2 from Pilot et al. [57], consistent with Thalmann et al. [56]. Moreover, our clades A and B coalesced 45 400 years ago (95% HDP = 35 100-57 200 years), coherently with Thalmann et al. [56] and Koblmüller et al. [65], suggesting that the entire mtDNA diversity in extant wolf populations has been generated during the last glaciations (the Würm glaciation in Europe), early Holocene. From a geographical point of view, the haplotypes of the highly diverse haplogroup B are widespread in Europe, Asia and North America (confirming [57]). In contrast, haplotypes in clade A have more limited distributions centered in southern-central European countries. This pattern, together with archaeological, morphological and genetic findings, suggest that extant Old and New World wolf populations expanded during or right after the last glacial period, balancing the local extinction of ancient wolf ecomorphs [53,56,[65][66][67][68]. However, it is not known if the turnover of wolf populations has been the consequence of a generalized megafaunal extinction wave due to climate changes or to hunting pressure and/or competition with expanding modern human populations [69][70][71]. Historical recolonization patterns might have been blurred by the extreme wolf potential to disperse or by recent anthropogenic eradication of many local populations [72]. Samples and genetic data of wolves from Eastern Europe and Asia are still scanty, preventing the reconstruction of the presumably complex population dynamics in those areas. For instances, while Italy is the single area where only haplogroup A is found, the poorly sampled sectors of Carpathian, Balkan and Caucasian regions show the presence of haplotypes from both A and B haplogroups, suggesting that these are (or were) areas of ancient population admixture (this study and [52,57]). Furthermore, the extinction of ancient mtDNA lineages and specialized wolf ecomorphs has been documented in North America, but it is not well described in Europe [67,73,74]. The lack of information on the extent of genetic and phenotypic variation in Paleolithic wolf populations in Eurasia makes the identification of wolf populations and the areas of origin of domesticated dogs difficult [75]. A better comprehension of these phenomena will help to reconstruct the wolf phylogeography in Eurasia that, at present, is still incomplete.
The phylogeographic patterns emerging from the mtDNA data are completed by multivariate and Bayesian analyses of the multilocus autosomal wolf genotypes. Results indicate that the main geographical populations in central and southern Europe are genetically distinct. The Italian and the Iberian wolves were the first ones to cluster separately in both multivariate and Bayesian clustering, confirming Lucchini et al. [13], Stronen et al. [52], and Pilot et al. [16]. Wolves sampled from populations in the eastern countries, and particularily from the Balkan and Carpathian regions show signatures of admixture, supporting previous mtDNA and autosomal SNP findings [16,52]. Recent genomic analyses [17] suggest that the most plausible demographic scenario should assume closely sequential splitting wolf population bottlenecks. ABC simulations indicate that southern European wolf populations (Iberian, Italian and Dinaric) might have split simultaneously at 20 490 and bottlenecked at 15 030 years ago (Table 3; Fig 4). Simulations also showed that all the current southern European wolf effective population sizes are significantly lower than sizes before bottlenecks with the Italian population that declined by c. 1.9 and the other populations c. 3.0-4.4 times, coherently with inferences from genomic data [53].

Conclusions
Despite the uncertainty that is usually associated to TMRCA estimates computed from limited number of markers (as in this study), or from limited number of samples (as in the genomic studies published so far [17,53], results concordantly show that currently fragmented wolf populations in south-western Europe were effectively isolated right after the last glacial maximum (c. 20 000-16 000 years ago). During this process, the Italian wolf population underwent a sharp demographic decline, which strongly reduced its standing genetic diversity. Thus, the result of the ancient isolation compromised the genetic pool of the Italian wolves long before the most recent anthropogenic bottleneck. However, these findings do not rule out the possibility that a recent demographic decline impoverished even more the genetic variation of this population.