Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Long-Range Gene Flow and the Effects of Climatic and Ecological Factors on Genetic Structuring in a Large, Solitary Carnivore: The Eurasian Lynx

Abstract

Due to their high mobility, large terrestrial predators are potentially capable of maintaining high connectivity, and therefore low genetic differentiation among populations. However, previous molecular studies have provided contradictory findings in relation to this. To elucidate patterns of genetic structure in large carnivores, we studied the genetic variability of the Eurasian lynx, Lynx lynx throughout north-eastern Europe using microsatellite, mitochondrial DNA control region and Y chromosome-linked markers. Using SAMOVA we found analogous patterns of genetic structure based on both mtDNA and microsatellites, which coincided with a relatively little evidence for male-biased dispersal. No polymorphism for the cytochrome b and ATP6 mtDNA genes and Y chromosome-linked markers were found. Lynx inhabiting a large area encompassing Finland, the Baltic countries and western Russia formed a single genetic unit, while some marginal populations were clearly divergent from others. The existence of a migration corridor was suggested to correspond with distribution of continuous forest cover. The lowest variability (in both markers) was found in lynx from Norway and Białowieża Primeval Forest (BPF), which coincided with a recent demographic bottleneck (Norway) or high habitat fragmentation (BPF). The Carpathian population, being monomorphic for the control region, showed relatively high microsatellite diversity, suggesting the effect of a past bottleneck (e.g. during Last Glacial Maximum) on its present genetic composition. Genetic structuring for the mtDNA control region was best explained by latitude and snow cover depth. Microsatellite structuring correlated with the lynx's main prey, especially the proportion of red deer (Cervus elaphus) in its diet. Eurasian lynx are capable of maintaining panmictic populations across eastern Europe unless they are severely limited by habitat continuity or a reduction in numbers. Different correlations of mtDNA and microsatellite population divergence patterns with climatic and ecological factors may suggest separate selective pressures acting on males and females in this solitary carnivore.

Introduction

Various studies on population genetics in large carnivores show diverse patterns of genetic structure. Some research emphasizes that, due to their high mobility, carnivore populations are maintained through long distance dispersal, which contributes to high gene flow across large spaces and low differentiation among populations [1][6]. There is also an increasing number of studies showing that gene flow does not occur equally in different directions, indicating the existence of various barriers for effective dispersal, including both physical obstacles and ecological conditions [7][16]. Genetic differentiation among carnivore populations has also been found to be related with the phylogeographic histories of species such as brown bear [17], [18] or wolf [19], [20].

Some of these discrepancies may result from using markers of different heritability, e.g. microsatellites and mitochondrial DNA (mtDNA) – the markers most frequently used in population genetic studies. For instance, a microsatellite-based study on the genetic structure of Scandinavian brown bears [21] found no support for the pattern of population subdivision previously based on mtDNA [22]. While mtDNA can be used to resolve taxonomy, questions on historical genetic variation and population structure, microsatellite markers are more suitable for inferring recent population history and contemporary gene flow [23].

From the research concerning the population genetic structure of the Eurasian lynx (Lynx lynx), two studies [24], [25] included both microsatellites and mtDNA. While Hellborg et al. [24] focused on understanding the causes of genetic differentiation within the Scandinavian and Baltic populations of lynx, the most recent study by Rueness et al. [25] covered almost the entire range of the species. This largest-scale study to date of spatial genetic patterns of the Eurasian lynx [25] has revealed a pronounced pattern of genetic structuring along the west-east axis of Eurasia, though the authors admitted the necessity of more fine-scaled sampling to account for the variability of ecological factors to explain the observed restrictions in gene flow. A previous study by Ratkiewicz et al. [26] based on the mtDNA control region, sampled the western edge of the species range (including the area between Norway and the Carpathian Mountains) and identified the existence of four genetic units that could have resulted from the marginal location of the populations studied, habitat fragmentation and/or demography. However, using a single genetic marker was not sufficient to assess the relative importance of these different factors. The importance of climatic factors on lynx genetic differentiation has been recently suggested by a review study [12]. Moreover, the Eurasian lynx was shown to shift its diet throughout its vast range with lynx in more northerly locations preying mainly on hares (Lepus timidus) (Russia, Finland) or domestic reindeer (Rangifer tarandus) (Norway) and the south-western populations feeding on Cervids [27]. This could potentially contribute to genetic divergence among these populations through local adaptations.

The aim of our study was to investigate the factors contributing to the patterns of genetic structuring in Eurasian lynx. To achieve this objective, we intensively sampled a large portion of the lynx's range (including sites at the periphery and within the core area of the species' distribution) and used three types of genetic markers: maternally inherited mtDNA, Y-chromosome markers and nuclear microsatellites. Our primary goal was to understand the relative role of factors affecting the genetic diversity in different populations. We predicted that: 1) lynx are capable of maintaining high rates of gene exchange providing that habitat is continuous; 2) genetic variability and distinctiveness of contemporary lynx populations in Eastern Europe are affected by both human activity (habitat fragmentation, population extermination) and past climatic events (Last Glacial Maximum); 3) genetic differentiation between lynx populations is correlated with climatic factors and/or 4) it is correlated with the main prey type of lynx. This knowledge will help in interpreting the importance of current and historical processes in shaping the distribution and genetic structure of the lynx population, as a model species of large, mobile carnivores.

Material and Methods

Ethics statement

The license for lynx live-trapping and blood sampling was obtained from the National Ethics Committee for Animal Experiments (no: DB/KKE/PL – 110/2001) and the Local Ethics Committee for Animal Experiments at the Medical University of Białystok (no: 52/2007). No animals were harmed during live-trapping and handling. An effort was made to minimize the time animals spent immobilized in the traps. All traps were equipped with radio- or GSM-alarm systems that allowed us to remotely control trapping and release the captured animals from the traps within 15 minutes to 1 hour.

Sampling

We analysed samples of Eurasian lynx from northern, central and eastern Europe bordered by Norway, the Carpathian Mountains and Russia (Kirov Republic) (Table 1). We arbitrarily assigned them to ten populations according to their geographic location and country of origin (Fig. 1). We are thus using the term “population” for simplicity, although we acknowledge that they may not represent populations in an ecological sense. One should note that Lithuanian lynxes were mostly sampled near to the Latvian boundary, while Belarusian lynxes were mostly sampled close to the Latvian and Russian borders (Fig. 1). Sampled lynxes span an area large enough to include a variety of habitat and climatic conditions. There is also a significant dietary diversification among the lynx populations over the sampled area (see S4 Table and references therein).

thumbnail
Figure 1. Distribution of Eurasian lynx sampling locations and mtDNA haplotypes.

Map showing distribution of sampling locations of the Eurasian lynx in north-eastern and central Europe. The color of each sampled individual denotes the haplotype of cr mtDNA and corresponds to the haplotype network in Fig. 2. Points are clustered into four groups as assigned by SAMOVA (both based on mtDNA and microsatellites) and shaded with different colors: Norway (yellow), BPF (green), Carpathians (blue) and remaining samples (pink). Names of arbitrarily assigned populations are given. Intensity of grey shading refers to the terrain ruggedness indicating mountainous areas. The background map was extracted from open access database available through USGS: http://srtm.usgs.gov/index.php). It is similar but not identical to the original image, and is therefore for representative purposes only.

https://doi.org/10.1371/journal.pone.0115160.g001

thumbnail
Table 1. Mitochondrial DNA- control region diversity indices for the Eurasian lynx samples studied.

https://doi.org/10.1371/journal.pone.0115160.t001

Altogether, 332 lynx individuals were sampled between 1992 and 2011. Sample sizes per population are given in Table 1. Samples were tissue collected from legally hunted individuals (Norway, Finland, Estonia, Latvia, Russia and Romania; n = 216), blood of live-trapped lynx (Białowieża Primeval Forest (BPF) in NE Poland; n = 18), hair samples (n = 54), animals found dead (Belarus, Lithuania, Poland (n = 25) and museum specimens (Poland, Slovakia, Ukraine; n = 19). Hairs were sampled non-invasively based on methods described by Schmidt and Kowalczyk [28]. Samples from northeastern Poland were assigned to two separate subpopulations (BPF and KARPF (Knyszyn, Augustów, Rominta, Piska Forests); Table 1) based on previous mtDNA results [26]. The Carpathian population consisted of samples from Poland, Slovakia, Romania and Ukraine which were pooled together due to the high geographic isolation of the Carpathian Mountains from the other populations in the distribution range and because only a single mtDNA haplotype shared by all individuals in this area (see Results).

DNA extraction, microsatellite amplification and error checking

Genomic DNA was extracted using the Genomic Mini Kit (A&A Biotechnology, Gdynia, Poland) according to the manufacturer's protocol and was stored at −20°C until used. The amplification of 13 microsatellite markers (twelve nuclear and one Y-linked) was carried out in three multiplex panels (Panel 1: Fca045, Fca090, Fca149, Fca391 and Fca559; Panel 2: Fca001, Fca008, Fca031, Fca043, F115; Panel 3: Fca077, Fca078 and Y-linked 278g21-4). The above-mentioned microsatellite loci were developed for domestic cats [29], [30] and successfully used in several lynx studies: [24], [25], [28], [31], [32]. The final reaction volume was 10 µl. Mastermix was created for 8-sample batches and included 36 µl Qiagen Multiplex Master Mix, 21 µl ultra-pure H2O (Qiagen, Hilden, Germany), 7 µl primer mix (containing 2pmol of each primer). Aliquots of 8 µl of this mix was then added to 2 µl of isolated DNA template. The thermocycling profile was as follows: 95°C for 15 min, followed by 30–45 cycles of denaturation at 94°C for 30 s, annealing at 57°C (Panel 1 and 3) or 56.5°C (Panel 2) for 90 s, extension at 72°C for 60 s and final extension at 60°C for 30 min. The PCR products were mixed with 10 µl ultragrade formamide and 0.2 µl GeneScan 500-LIZ size standard (Life Technologies, Carlsbad, CA, USA), denatured at 95°C for 5 min, immediately cooled and then separated using a four-capillary ABI 3130 Genetic Analyser (Life Technologies, Carlsbad, CA, USA). Allele sizes and genotypes were scored using GeneMapper 4.0 software (Life Technologies, Carlsbad, CA, USA). PCRs were repeated twice to test the consistency of the results for a subset (5%) of randomly chosen individuals. A total of 54 hair samples from Lithuania and Belarus were analyzed. Hair samples which gave low amplification rate and displaying contamination within hair traps were excluded. Microsatellite panels in the remaining hair samples (20) were amplified between two and four times until they gave reliable genotypes, e.g. heterozygotes at least three times out of four. Samples with inconsistent genotypes were discarded from the analysis. Possible genotyping errors due to stuttering, short allele dominance and null alleles were tested by MICRO-CHECKER 2.2.3. [33]. Multilocus genotypes of the hair samples were analysed in CERVUS 3.0.3 [34] to avoid using samples of the same individuals. Obtained genotypes are presented in S1 Appendix.

mtDNA and Y-chromosome analyses

A portion (613 bp) of the mtDNA control region (cr mtDNA) was amplified in 122 lynxes using the primer pair LGL283 and ISM015 [26], [35]. Additionally, 210 lynx cr mtDNA sequences from Ratkiewicz et al. [26] were used. Almost the entire cytochrome b gene, e.g. 1094 bp (mtDNA) was amplified for 12 individuals from Russia, Finland, Estonia, Latvia and Poland using lynx-specific primers: LcytbF: 5′ -CAC ATG GAA TTT AAC CAT GAC C -3′ and LcytbR: 5′-GAC TCT TCA TTT GAG GAG ACG. The entire 681 bp of the ATP6 gene (LATP_F: 5′-TCC AGA ACC TAA ATC CAC AAC C-3′ and LATP_R: 5′-GCA TGA GTT TGG TGG GTC ATT A-3′) was amplified in 8 lynxes from Norway, Finland, Latvia and Poland. All laboratory procedures and quality test controls for cr mtDNA followed [26]. We also amplified seven universal Y-linked markers (YCATS, Y-chromosome conserved anchor tagged sequences) according to [36], with the exception of marker DBY4, which had an annealing temperature of 48°C. We obtained clear PCR products for DBY4 (185 bp), DBY7 (275 bp), DBY8 (115 bp), SMCY7 (625 bp) and SMCY17 (125 bp) but no clear products were obtained for UBE1Y6 and UTY11. The five Y-chromosome markers that gave clear PCR products were sequenced in 43 males from Russia, Finland, Estonia, Latvia, Poland, Romania and Belarus according to [36].

Statistical analyses

For the autosomal microsatellite loci, FSTAT 2.9.3 [37] was used to estimate the number of alleles per locus (NA), allelic richness (AR), gene diversity (Hs) and inbreeding coefficient (FIS) in the populations. Differences among groups of samples in AR, Hs, expected heterozygosity (He), relatedness coefficient (r), and FST were tested using a permutation procedure (10000 iterations) in FSTAT. The same software was used to tests for linkage disequilibrium among 12 microsatellite loci within each population and for the total sample (13200 permutations). Hardy–Weinberg equilibrium probabilities were calculated using GENEPOP 4.0 [38]. Estimates of cr mtDNA haplotype diversity (h) and nucleotide diversity (π), number of segregating sites (s) and mean pairwise difference were calculated using ARLEQUIN v 3.1 [39]. Relationships among haplotypes were represented as a haplotype network obtained with the statistical parsimony method using the TCS v 1.21software [40]. Pairwise FST (for microsatellites) and ΦST (for mtDNA) values for the 10 populations were estimated using ARLEQUIN and their significance was tested using 1000 permutations corrected for multiple tests by Bonferroni correction. For comparative purposes, we have also calculated pairwise DEST [41] using SMOGD v. 1.2.5 [42] (results shown only as Supporting information (S3 Table)). Isolation by distance (IBD) patterns were determined by comparing pairwise FST/(1- FST) [43] for both nuclear and mtDNA markers to the logarithm of geographical distance (measured in a straight line between the central point of each population except the distances to Norway and Finland which were taken along an arbitral broken line omitting the Baltic Sea) using Isolation by Distance Web Service (IBDWS) (http://ibdws.sdsu.edu/~ibdws/ [44].

Principal coordinate analysis (PCA) was performed based on a sample-wise matrix of genetic distance for microsatellite and mtDNA data in GenAlEx v. 6.0 [45]. We then used PC1 scores with the geographic coordinates of lynx samples for the spatial analysis of genetic differentiation, using the kriging algorithm in the Surfer 10 software (http://www.goldensoftware.com/demo-downloads). Contour lines of the first axis of the PCA were interpolated and superimposed onto a geographic map of the study area. By interpolating the PC1 scores for each individual lynx throughout the study area, we tried to identify regions where genetic dissimilarity between individuals was considerably higher or lower than would be expected from the IBD effect alone. Such geographic areas may represent dispersal barriers or migration corridors, respectively.

Analysis of molecular variance (AMOVA; [46] using ARLEQUIN (with 10000 permutations) was performed to assess structuring within the data, where the sampling sites were grouped as a single population (using both, ΦST and FST). To explore patterns of genetic divergence in more detail, we applied the spatial AMOVA procedure using SAMOVA ver. 1.0 [47]. This allowed us to identify the grouping of sampling sites that maximized the ΦCT values (among group variance) and minimized ΦSC values (among populations within group variance). The significance of Φ-statistics was tested using 10000 permutations for K = 2 to K = 9 partitions of the sampling sites. The pRST - RST tests between identified groups were performed in SPAGEDI [48], using 20000 allele size permutations to identify if genetic differentiation was due to genetic drift or other evolutionary processes.

Contemporary migration rates (within the past one to three generations) between lynx populations were evaluated using a Bayesian approach implemented in BayesAss ver. 1.3 [49] according to the authors' recommendations. The analyses were performed with 3×106 iterations of which 106 were burn-in and a sampling frequency of 2000. We further applied a Bayesian clustering approach to infer the most likely number of genetically distinct groups of samples using the software STRUCTURE 2.3.4. [50] without prior information of the sampling locations. We assumed the admixture model with correlated allele frequencies, and specified burn-in of 106 iterations and 105 Markov Chain Monte Carlo (MCMC) replicates after the burn-in period. The program was run 5 times for each number of clusters (K) between 1 and 10 to verify the consistency of estimates across runs. To determine most likely number of genetic clusters we evaluated the probability logarithm of the data (lnP(D); [50]) and Delta K (ΔK; [51]) for each K averaged across runs, using the software Structure Harvester (Web version: v0.6.93, [52].

We applied Approximate Bayesian computation (using the program DIYABC 2.0.4; [53]) to estimate the relative likelihood of possible scenarios for the lynx population histories. The program uses reference tables (containing parameters based on known or estimated values) to establish scenarios from which simulated data sets could be compared to the observed values. As the program requires a priori defining populations, for the analyses we generally followed the subdivision of the lynx population suggested by the STUCTURE analysis into four subpopulations: 1) Norway, 2) BPF-KARPF, 3) Baltic and 4) Russia-Finland. However, we decided to distinguish a fifth sub-population, the Carpathian population, due to the analysis in SAMOVA. We tested nine possible scenarios. We set the effective population sizes (Ne) from 10 to 10000. Each competing scenario was given equal prior probability. We accepted the default mutation rate model prior distributions suggested in the software for microsatellites. For the mitochondrial dataset we assumed the mean mutation rate of 5.35×10−7 substitutions/site/year with SD  = 2.28×10−7 substitutions [54]. We set the number of simulated data sets at 500000. We compared the relative likelihoods of the scenarios by the logistic regression approach, with a 2% subset of the closest simulated data. We provisionally assessed the fit of the model to the data using principal components analysis (PCA).

We tested for a possible sex-bias in dispersal using a spatial autocorrelation analysis [55] in GenAlEx 6.5 [45] for each sex separately. The autocorrelation coefficient, pairwise r values, were divided into 20 distance classes (100 km each) with mean sample sizes in each distance class N = 254 (±182). A null distribution of r values for each distance class was obtained by permutation (N = 9999) and the confidence intervals (C.I.) about r were estimated by bootstrapping with replacement (N = 9999), and plotted in a correlogram. The extent of the detectable spatial genetic structure was approximated as the distance class at which r was no longer significant and the intercept crossed the x-axis.

To identify correlations between genetic distance between lynx populations and abiotic climatic, as well as biotic factors (lynx basic prey), we performed multivariate multiple regression analysis using DISTLM ver. 5 [56]. The mtDNA and the microsatellite population pairwise FST values were used to construct the response matrix and tested against the following predictor matrices: geographical distance, average temperature in January, number of days with snow cover (DSC), average snow cover depth (SCD) in January [57], and the winter index for the North Atlantic Oscillation (NAO) as an index of climate [58]. Marginal tests of each predictor were done, followed by conditional tests, where latitude was included as a covariable to the predictor variables. Then, sequential tests were done using a forward selection procedure to produce a combined model of genetic differentiation in the Eurasian lynx using DISTLM forward ver. 1.3 [59]. The P values were obtained from 10000 permutations. Abiotic climatic data were taken from the WorldClimate database (www.worldclimate.com), [60] and the Climate Atlas of Poland [60]. Data obtained from the literature on lynx diet is cited in S4 Table.

Results

Intrapopulation genetic variability

We obtained control region mtDNA sequences for 332 lynx samples. Sixteen 613 bp haplotypes of the cr mtDNA were observed, of which four (H14, H15, H16 and H18) were previously unreported (Fig. 2). All new haplotypes are provided in GenBank (accession nos.: H14 - KM000076, H15 - KM000077, H16 - KM000078, H18 - KM000080). The number of segregating sites (S) was 18 (17 transitions and one indel). The number of haplotypes per population ranged from 1 (Norway and the Carpathians) to 8 (Russia; Table 1). Cr mtDNA H1 was the most widespread haplotype and it was present in eight of the sampled locations (Figs. 1 and 2). It was fixed in Norway and very common in Estonia, Finland and Russia. On the other hand, H4 was fixed in the Carpathians and relatively common (0.231) in Russia and was also present at moderate to low frequencies in Estonia, Latvia, Belarus and KARPF in Poland (Fig. 1, see S1 Table for haplotype frequencies). It is worth noting that three haplotypes were unique to Russia and, surprisingly, three were exclusively present in the Lithuanian samples. Haplotype (h) and nucleotide (π) diversity in the overall sample was 0.78 and 0.388%, respectively (Table 1). The haplotype diversity (h) values ranged from zero (Norway and the Carpathians) to 0.82 (Belarus). The Norwegian and Carpathian populations showed no nucleotide diversity (π = 0.000%), whereas this parameter was highest in Latvia (0.468%).

thumbnail
Figure 2. mtDNA haplotype network of Eurasian lynx.

Haplotype network illustrating the relationship among 16 haplotypes of Eurasian lynx. Small black circles indicate missing haplotypes. Numbers denote the haplotypes. The size of the circles (except the haplotypes not found in this study) refers to the relative frequencies of a given haplotype in the whole sample. Colors of haplotypes correspond to Fig. 1.

https://doi.org/10.1371/journal.pone.0115160.g002

The cytochrome b gene (1094 bp) showed no polymorphism among 12 individuals from different populations (GenBank accession no.: KM000081). Likewise, no polymorphism was found among 681 bp of the entire ATP6 gene among 8 lynxes from different populations (GenBank accession no.: KM000083).

No sequence polymorphism was found in the five Y-linked chromosome markers (total: 1386 bp; GenBank accession no.: KM000084, KM000085), among 43 male lynx from seven sampling locations (Russia, Finland, Estonia, Latvia, Poland, Romania and Belarus) and the Y-linked microsatellite locus 278g21-4 possessed only a single allele (146 bp long) in all males studied.

The number of alleles for the microsatellite loci ranged from four to 18 and we did not find any significant (Bonferroni-corrected) linkage disequilibria between loci within any of the populations studied, as well as the frequencies of null alleles low (below 0.05, data not shown). Measures of genetic diversity for the microsatellite loci (A, AR and He) were lowest for the Norwegian and BPF populations, and largest for lynx from Russia (Table 2). The A, AR and He values for Norway and BPF were significantly lower (permutation test, p<0.05) than for the remaining lynx populations, possibly due to bottlenecks. Interestingly, the Carpathian population did not show any signs of reduced levels of polymorphism at the microsatellite loci; however it exhibited a positive and significant FIS value, as did the population from Belarus (Table 2).

thumbnail
Table 2. Microsatellite DNA diversity indices for the Eurasian lynx in Europe.

https://doi.org/10.1371/journal.pone.0115160.t002

Genetic differentiation among lynx sampling locations

Pairwise genetic differentiation values between sampling locations for the cr mtDNA ranged from 0.01 to 1.00 (ΦST) and the majority of comparisons were significant (Table 3). These results indicate that a high degree of genetic differentiation exists between lynx populations, especially those from Norway, BPF and the Carpathians. Nonsignificant values of ΦST were found between Estonia and Latvia, Finland and Estonia as well as Latvia and KARPF, whereas lynx populations in close geographic proximity in Lithuania and Latvia showed considerable and significant mtDNA differentiation (Table 3). For the cr mtDNA FST was not significantly different from zero between three population pairs only: Estonia and Russia, Latvia and Belarus, Belarus and KARPF (S2 Table). Similarly, various degrees of genetic differentiation were found at the microsatellite loci, with pairwise FST values ranging from 0.012 to 0.293 and the highest divergence between the populations of Norway, BPF and Carpathians with the other sampling locations (Table 3). All comparisons except one (the Latvia-Lithuania population pair) were significant. Comparative patterns of differentiation were obtained with DEST which ranged from 0.00 to 0.32 (S3 Table). The RST values (S3 Table) ranged from 0.004 to 0.234 and the majority of comparisons were significantly different from zero. The RST values were non-significant between the following population pairs: Finland - Lithuania, Estonia - KARPF, Estonia - Lithuania, Latvia - Lithuania, Lithuania - Belarus, Lithuania - KARPF, Lithuania - Carpathians, Lithuania - Russia, KARPF - BPF, KARPF - Russia.

thumbnail
Table 3. Pairwise differentiation between Eurasian lynx populations.

https://doi.org/10.1371/journal.pone.0115160.t003

Spatial autocorrelation among females and males showed a pattern of decreasing relatedness with increasing distance in the first two distance classes (0–200 km), with an x-intercept at 498 km and 586 km for females and males, respectively (Fig. 3). For females, the r values were positive and significant for the first four classes (up to 400 km), while for males the corresponding values were significant up to 500 km. Thus, the approximate scale of the positive genetic structure was between 400 and 500 km. The overall shape of the correlogram was similar for males and females (Fig. 3), though with slight indication of male-bias in dispersal distance. The relationship between r and distance was not significantly negative until 700 km for females and males (Fig. 3).

thumbnail
Figure 3. Results of the spatial autocorrelation analysis.

Correlograms of the average autocorrelation coefficient (r) for 20 distance classes of 100 km each for male (A) and female (B) Eurasian lynx. The dashed lines represent the 95% upper and lower bounds of the null distribution assuming no spatial structure. The error bars represent the 95% confidence intervals about r. Significant spatial structure is observed when r exceeds the null distribution and the error bars do not overlap zero.

https://doi.org/10.1371/journal.pone.0115160.g003

We did not find any statistically significant isolation by distance (IBD) pattern for mtDNA (r2 = 0.028, p = 0.28) nor for microsatellite loci (r2 = 0.076, p = 0.17) among our sampling locations (S1 Fig.). When we excluded lynx from Russia, the IBD pattern became significant for microsatellite loci (r2 = 0.394, p<0.05) but not for the cr mtDNA (r2 = 0.178, p = 0.08) (S1 Fig.).

Geographical structuring among lynx cr mtDNA haplotypes in the study area was highly supported by AMOVA results where all sampling sites were treated as a single group (ΦST = 0.570, p<0.001, FST = 0.464, p<0.001). SAMOVA was subsequently used to identify the subdivision that most likely explains the mtDNA structure observed in lynx. The data were best explained assuming four groups of lynx populations: (1) Norway, (2) BPF, (3) Carpathians and (4) the remaining 7 sampling sites (Fig. 1). This configuration maximized the among-group variation (44.24%) and minimized the variation among sampling locations within groups (13.67%); the variation within sampling locations was 42.10% (ΦCT = 0.442, p<0.05, ΦSC = 0.245, p<0.001, ΦST = 0.579, p<0.001).

For the microsatellite data, AMOVA supported the geographical structuring of the populations (FST = 0.10, p<0.001, RST = 0.09, p<0.001). The data were best explained in SAMOVA analysis assuming the same four groups of lynx populations: (1) Norway, (2) BPF, (3) Carpathians and (4) the remaining 7 sampling locations (Fig. 1). This configuration maximized the among-group variation (11.81%) and minimized the variation among sampling locations within groups (4.38%); the variation within sampling locations was 83.81% (FCT = 0.118, p<0.01, FSC = 0.050, p<0.001, FST = 0.162, p<0.001). RST did not differ significantly from pRST for any pairwise between-group comparisons (P>0.27, 1000 permutation tests).

Using the Bayesian clustering method based on multilocus lynx genotypes with STRUCTURE indicated that four genetic clusters were most likely using the ΔK method (S2 Fig.). However, with the probability logarithm of the data (lnP(D) method, although the likelihood values leveled off at K = 4, it reached a maximum at K = 8 suggesting further subdivision of the lynx population (S2 Fig.). About 75% of the individuals were assigned to genetic clusters using Q>0.80 as a threshold, with considerable genetic mixing between two of the different clusters. The genetic assignment plots produced in STRUCTURE clearly highlighted the genetic discontinuity among lynx populations. Assuming K = 4, the samples were assigned according to the following pattern: two sampling locations (the Norway and the North-Eastern Polish lynx, including the BPF and KARPF) were nearly uniformly composed of two different single genetic clusters, a third cluster was composed of lynx from Finland, Russia, Carpathians and to some extent from Belarus, whereas the fourth cluster included the remaining lynx (Estonia, Latvia, Lithuania and Belarus) (Fig. 4). Interestingly, irrespective of the number of clusters assumed, the Norwegian and the North-Eastern Polish lynx populations always formed clearly distinguishable genetic units.

thumbnail
Figure 4. Bayesian clustering of the Eurasian lynx population in north-eastern Europe.

Results of the STRUCTURE analysis assuming number of clusters from K = 3 to K = 10. Numbers correspond to classification of populations in Table 1. Note consistency of populations: 1 (Norway), 7 and 8 (KARPF and BPF, both NE Poland) being composed of two uniform genetic clusters across all values of K.

https://doi.org/10.1371/journal.pone.0115160.g004

Contemporary gene flow and possible scenarios of past lynx population history

Bayesian estimates of contemporary migration estimated in BayesAss were close to zero for Norway, BPF and the Carpathians. Some migration probably occurred from BPF into the KARPF population and to a certain extent from BPF to Lithuania. The lynx population from Belarus could act as a source population for the neighboring populations in the Baltic states. Regular recent migration was evident between Finland and Russia in both directions (Fig. 5).

thumbnail
Figure 5. Migration routes and rates of Eurasian lynx across north-eastern Europe.

Recent migration rates within the Eurasian lynx population in north-eastern and central Europe between arbitrarily assigned sampling populations, estimated using BayesAss. Directions and rates of migrations are shown with arrows and associated numbers. The numbers within the circles denote proportions of non-immigrants within the sampled populations (denoted with small digits which refer to the population numbers and names in Table 1). Grey shading represents the forest cover. It is prepared based on an open access GlobCOVER database (http://due.esrin.esa.int/globcover/) by extracting a range of data indicating forested and non forested areas (limited to values of 40–110). It is thus similar but not identical to the original image, and is therefore for representative purposes only. The lynx range (after Von Arx et al. [91]) is marked by continuous lines.

https://doi.org/10.1371/journal.pone.0115160.g005

The synthetic genetic map based on PC1 scores calculated from pairwise genetic distance between individuals, allowed us to identify potential dispersal barriers and migration corridors (Fig. 6). Microsatellite data set revealed areas of abrupt genetic change, indicating that possible barriers to lynx dispersal may exist in: (1) the northern part of the Scandinavian peninsula, (2) the gulf of Finland, (3) the highly fragmented habitats around NE Poland and (4) the deforested area north to the Carpathians (Fig. 6A). On the other hand, fairly homogeneous genetic landscape was visible between southern Finland, eastern Belarus and western Russia, suggesting a possible migration corridor in this part of the study area (Fig. 6A). The analogous migration corridor was also revealed based on the mtDNA data set, but in contrast to microsatellites, no barriers were detected in northern Scandinavia and north of the Carpathians (Fig. 6B). Relatively high genetic dissimilarity between sampled lynxes was visible on the mtDNA contour map near the gulf of Finland and in NE Poland (Fig. 6B).

thumbnail
Figure 6. Results of spatial analysis of genetic differentiation in the Eurasian lynx.

Maps showing contour lines of the first axis of the PCA performed on a sample-wise matrix of genetic distance for microsatellite (A) and mtDNA (B) data interpolated with the use of the kriging algorithm (Surfer12 software) and superimposed onto a geographic map of the study area. Black and white arrows show possible migration barriers or migratory corridors respectively. See Fig. 5 for other explanations.

https://doi.org/10.1371/journal.pone.0115160.g006

Using the approximate Bayesian computation with DIYABC we found highest support for two different (but very consistent through various approaches; see S2 Appendix), mutually exclusive scenarios of lynx population history depending on the genetic marker used. The analysis based on microsatellites suggested that the lynx from the Carpathian mountains could have been the source for all remaining populations (Fig. 7, Scenario 2). The first split, according to this scenario, has lead to the establishment of the North-Eastern Polish population of lynx (BPF and KARPF). During the second split a branch containing the Russian, Norwegian and Finnish lynx split off from the NE Polish branch. The latter population diverged also at a subsequent split giving rise to the Baltic population. During the last, most recent split, the Norwegian lynx separated from the Russian-Finnish branch.

thumbnail
Figure 7. Possible scenarios of Eurasian lynx population history.

Results of the approximate Bayesian computation conducted with the use of DIYABC 2.0.4 [53] to estimate the relative likelihood of alternative scenarios for the lynx's colonization history. A) Graphs in the upper panel illustrate the four final best supported scenarios from among nine scenarios proposed in the analyses (see S2 Appendix for details). The ten sampling locations analyzed in this study were grouped into five populations based on its distinctness suggested by the STRUCTURE and SAMOVA analyses: 1) Nor (Norway), 2) NE PL (BPF and KARPF), 3) Balt (Latvia, Lithuania, Estonia, Belarus), 4) Rus (Russia and Finland) and 5) CARP (Carpathians). Colors in the colonization scenarios indicate different (but unknown) population sizes (Ne). Zero means sampling time and t1-t4 mean relative times of past events of suggested population splitting. Graphs in the lower panel indicate the relative likelihoods of the four best scenarios compared by a logistic regression based on microsatellite (B) and mtDNA (C) data with a 2% subset of the closest simulated data. The best support is for the scenario 2 and 9 for microsatellite (B) and mtDNA (C) data, respectively. Note, although four scenarios were analyzed for mtDNA (C) the posterior probability of scenario 4 is not visible due to highly similar values with scenario 6.

https://doi.org/10.1371/journal.pone.0115160.g007

The mtDNA-based DIYABC analysis, however, showed that it is the Russian-Finnish population that was the source population from which the Norwegian lynx has split first (Fig. 7, Scenario 9). The Russian-Finnish lynx gave also rise to the Carpathian and Baltic populations at the second and third splits, respectively. During the most recent split the Baltic population has given rise to the NE Polish lynx.

Possible effect of ecological factors on lynx genetic structure

A test of the influence of geographical location on genetic differentiation among the lynx sampling locations showed that pairwise ΦST values (mtDNA) is correlated with latitude (46.5% of variation explained, P<0.001) but not with longitude. Similarly, the number of days with snow cover, snow cover depth, and the NAO index, as revealed by marginal tests, had significant effect on mtDNA differentiation among sampling locations (Table 4). When geographical coordinates were incorporated as covariables into the multiple regression analysis (conditional tests), only snow cover depth was correlated with the pairwise ΦST values between the lynx sampling locations. On the other hand, sequential tests showed that latitude and snow cover depth were the only significant factors (Table 4). The effect of snow (SCD and DSC taken as a set of variables) explained 59.3% of the variation (p<0.01) for the mtDNA data in marginal tests. No such relationship was found for the microsatellite data except in one analysis. There was a significant effect of the main prey and the share of red deer in the lynx diet explaining 82.5% (p<0.05; the same for sequential tests) and 26.6% (p<0.05) of the variation, respectively, in the marginal tests for microsatellites only (Table 4).

thumbnail
Table 4. Effects of abiotic climatic variables on genetic differentiation based on mtDNA and microsatellites among Eurasian lynx populations in Europe.

https://doi.org/10.1371/journal.pone.0115160.t004

Discussion

Population differentiation

The important finding of our research is the presence of similar patterns of genetic differentiation among populations based on two types of markers: mtDNA and microsatellites in the Eurasian lynx. The entire sample was consistently subdivided into four genetic units. SAMOVA also suggested that while lynx from a large area encompassing Finland, the Baltic countries, Belarus and Russia were all included in a single genetic unit, the three remaining sampling locations – Norway, the Carpathians and the Białowieża Primeval Forest (BPF) were separated into three clearly distinct genetic entities. Based on the different heritability of mtDNA and microsatellites, as well as previous research that suggested different genetic structure in the brown bear with use of the same types of markers [21], [22], we expected to find less consistency when considering both mtDNA and microsatellite variability in the Eurasian lynx.

Spatial autocorrelation analysis provided relatively little evidence for male-biased dispersal in the lynx. This is not fully concordant with the available telemetry data, suggesting that similar to other large felids [61][63], male lynx disperse more frequently and for longer distances than females [64][65]. This may mean lower effective dispersal in this species than expected from observational data. The existence of microsatellite genetic structure almost as strong as that found with mtDNA (which is inherited solely through females) indicates that there may be effective physical barriers to male dispersal and that female philopatry alone does not explain the observed pattern. Such barriers are likely to function in the form of habitat fragmentation. It is particularly visible on the synthetic genetic map in the case of the BPF population, which showed a remarkable genetic distinctness from the neighboring populations despite the short geographic distance between them. Indeed, the habitat seems most severely fragmented in this part of the lynx range [12], [66], and its effect was already reported in an earlier study based on six microsatellites [31]. Concordant genetic structure in the two types of markers was also recently found by Czarnomska et al. [14] in the wolf population inhabiting Poland and was explained by factors including natal-habitat biased dispersal and habitat disparities. The effect of habitat discontinuity on dispersal and therefore the population genetic structure in lynx may be more direct due to their stronger dependence on forest habitat [67], [68] as compared to that of wolves [69]. It is thus likely that effective dispersal may not be sufficient to counteract the effect of genetic drift in local, isolated populations of the Eurasian lynx. Indeed, no significant differences between RST and pRST values for microsatellite loci suggest genetic drift as the primary reason for observed genetic divergence among the lynx populations at microsatellite nuclear loci.

The strong genetic structure in the Eurasian lynx revealed by our analyses corresponds well with our previous study on mtDNA that showed large differentiation among the populations at the westernmost peripheries of its range [26]. However, sampling more locations in this study, particularly from the core of lynx range (Russia, Kirov), and applying microsatellite markers gave us more explanatory power with regards to plausible factors influencing the observed pattern of genetic variability in this carnivore. The isolation by distance (IBD) pattern appeared significant only for microsatellites for the data set without the Russian population. This suggests that the population differentiation at the westernmost part of the species' range may be dependent on the dispersal capabilities of lynx and IBD may thus explain only a portion of the observed pattern of population differentiation. With the inclusion of the Russian subsample the IBD was no longer significant. This suggests that there must be no particular genetic divergence between the core (Russia) and some peripheral parts (e.g. Finland, Belarus or Estonia) of lynx's range. Thus, the Eurasian lynx has the potential to maintain a panmictic population encompassing the vast territory of Eastern Europe, unless it encounters barriers due to habitat fragmentation. The spatial, individual-based, analysis of genetic distance showed that it is plausible that deforested areas at the western margins of the Eurasian lynx range act as effective barriers to dispersal between the BPF, Carpathians and the other studied populations. The high genetic gradients on the synthetic genetic maps along with low genetic diversity (for the microsatellites) of the Norwegian and the BPF samples suggests that they have both experienced genetic drift and are particularly restricted by barriers to dispersal. The barrier between northern Norway and Finland is most probably caused by intense hunting pressure in areas where reindeer are present [70].

Our study has also indicated that the genetic structure in the Eurasian lynx may be explained, at least in part, with climatic conditions, particularly the depth of snow cover, which was significantly correlated with mtDNA population differentiation, even when geographic coordinates were included as co-variables in the conditional tests. The lack of a barrier on the synthetic genetic map for mtDNA in northern Scandinavia throughout Finland and Russia coincides well with similar climatic conditions (e.g. high snow cover depth) in all northern lynx populations. Very similar results have recently been reported based on mtDNA for a smaller dataset of lynx samples from the western part of its range [12]. The mechanism suggested for this relationship involves the matrilineal inheritance of mtDNA as well as the fact that the climatic factors may also have affected the kill rate by female lynx [71]. Moreover, skull size in female lynx, but not males, was also found to be influenced by climatic factors [72]. Females in solitary cats are provisioning their offspring with food alone [73], so familiarity with local snow conditions may indeed influence their hunting efficiency and in effect the survival and development of their kittens.

The relationships between microsatellite population differentiation and climatic factors were less clear, with the only significant effect of snow conditions found in the sequential test. On the other hand, a significant effect on microsatellite differentiation in the marginal and sequential tests were found for prey, especially the proportion of red deer in the lynx's diet. Although drawing firm conclusions from the latter result can be somewhat speculative, red deer are killed more often by lynx in North-Eastern Poland (BPF and KARPF) than in other areas (S4 Table and citations therein), which coincides with the genetic distinctness of the lynx inhabiting this region.

The results of multivariate multiple regression, along with the IBD effect may suggest that despite generally concordant patterns of population subdivision, there might be different factors such as climate and prey type affecting mtDNA and the microsatellites genetic structure of lynx population, respectively. The effect of climate may have resulted in the matrilineal clusters of lynx being specifically adapted to hunting prey under local snow conditions. Moreover, as male lynx are more likely to kill larger prey [74], [75] one cannot exclude the possibility of lynx microsatellite genetic structure being additionally affected by male adaptations to killing red deer. The influence of prey specialization on predator population structure has already been well documented in the wolf [10], [11], and inferred for the Canada lynx (Lynx canadensis) [8].

Gene flow among lynx populations

This study revealed differing levels of gene flow within the range of the Eurasian lynx. Low genetic divergence and high rate of migrants identified between the eastern part of European Russia and some western lynx sampling locations suggested high gene flow along the east-westward direction, despite the distance between locations. Several independent analyses of our microsatellite data indicated possible migration corridors situated between Russia (Kirov and Arkhangelsk Oblast), Finland and Belarus. In contrast, little genetic exchange appears to have occurred between neighboring populations situated at the westernmost peripheries of the lynx range, e.g. Norway, the BPF and the Carpathians.

The possibility of recent gene flow between lynx inhabiting the far east of European Russia and the westernmost edge of the lynx range (approximately 1100–1300 km) may not be surprising considering the long distances large carnivores (including lynx) are capable of covering during dispersal and migrations [2], [63], [65], [76]. Nevertheless, our findings indicated, for the first time, the effectiveness of this mobility in a genetic sense. Previous genetic studies conducted at smaller scales were mainly focused on documenting how gene flow was limited in this carnivore, despite the high dispersal potential [24], [25], [26], [31].

Studies documenting long-distance, cross-continental movements in wildlife are rare due to difficulties in collecting representative material [2], [10], [15], [16]. The existence of a migration corridor in Eurasia was revealed recently for brown bear between Eastern European Russia and Finland by Keis et al. [15] and the authors suggested that the corridor they detected coincided with the southern border of the taiga biome. High gene flow across a large area of Eastern Europe was also found in the grey wolf [10], [16]. The feasibility of gene flow documented in large carnivores corresponds with the distribution of continuous forest cover, which spans Eastern Europe from Latvia, Estonia and Finland to Eastern Russia [66] and the habitat preferences of all three of these large carnivores [67], [68], [77].

The lynx population inhabiting eastern European Russia was characterized by high genetic variation, and therefore one may expect it to constitute a possible source population for these cats in Europe. However, we found present-day gene flow between Russian and Finnish lynx in both directions with migration rates slightly higher from Finland to Russia than in the opposite direction. The pattern of lynx population history revealed with DIYABC analysis of the mtDNA data supports the scenario of a Russian origin for the European lynx populations.

Genetic variability and contemporary and historical bottlenecks

The highest variability at both types of markers was found in sampling locations grouped into clusters situated within the core lynx range (Russia) or being directly connected to it through the continuity of the habitat present (Latvia, Estonia). Whereas high genetic variability in this cluster can be directly linked with the central location and habitat continuity [78], the reasons for the very low variability in each of the remaining three clusters (consisting of BPF, Norway and Carpathians) may have a unique reason in each case.

Among the three most divergent and least variable populations, the Norwegian and Carpathian lynx present different patterns of genetic variability. Whereas the lynx in both areas showed a complete lack of polymorphism in mtDNA, they had different levels of microsatellite variability. As the Norwegian lynx had one of the lowest reported values for allelic richness and heterozygosity, it is clearly concordant with the relatively recent bottleneck experienced by this population during the twentieth century [24], [79]. In contrast, the Carpathian lynx showed relatively high genetic variability in microsatellites, which indicates that no dramatic reduction of population numbers occurred recently. In this case the presence of only one haplotype of cr mtDNA may support the northern refugium hypothesis for the lynx with small effective population size during the last glacial maximum (LGM, [26], [80], [81]. This scenario is also in line with the statistical analyses in which the SAMOVA indicated that the Carpathian lynx formed a separate genetic group for both mtDNA and microsatellites, and the DIYABC analysis of microsatellite data, which suggested the Carpathian population to be the source for all the other lynx populations. The best supported scenario in this study, which implied the spread of the lynx from Carpathian mountains north and eastwards was also recently suggested by Rueness et al. [25]. The Bayesian approach in STRUCTURE suggested the Carpathian population to be composed of a single genetic cluster along with the Russian lynx. The haplotype H4 was found to be common in other sampling locations and frequent in the core population (Russia). Thus our results suggest that the present monomorphism in cr mtDNA in the Carpathian lynx could have been attributed to the bottleneck effect during the LGM. It is likely that the lynx population became extinct north of the Carpathians, and the few survivors reached these mountains and quickly established the population there. This might have deterred possible female migrants from settling in the existing social structure of the population [12] after the glacier retreated and resulted in the fixation of the single haplotype.

In the case of the BPF population, it is most likely that both its peripheral location and strong isolation from the core range were the main reasons for the low microsatellite and mitochondrial variability. A similar relationship between genetic variability and location relative to the core of the species range was reported for the Canada lynx [78] and the effect was explained with small population size and limited connections to other populations, which is typical in peripheries of a species' range [82]. The limited genetic variability observed in the BPF population might also be related to the short-term bottlenecks that most probably occurred during 19th and 20th centuries [83], causing genetic drift.

Conclusion

The Eurasian lynx is capable of maintaining a panmictic population across eastern Europe unless they are severely limited by habitat continuity or a reduction in numbers. The effect of climatic conditions (snow cover) on the genetic divergence between populations with respect to mtDNA, but not nuclear microsatellite loci, may suggest different selective pressures acting on males and females in solitary carnivores. However, further research is needed to identify such selective pressures.

Supporting Information

S1 Fig.

Results of isolation by distance (IBD) analysis of the Eurasian lynx populations. Two upper panels show IBD for microsatellite data for all ten putative populations (A) and excluding Kirov (Russia) (B). Two lower panels show IBD for mtDNA data for all populations (C) and excluding Kirov (Russia) (D). A significant relationship was only obtained for microsatellite data after excluding the Russian population.

https://doi.org/10.1371/journal.pone.0115160.s001

(PDF)

S2 Fig.

Number of genetic clusters of the Eurasian lynx population in north-eastern Europe. Evaluation of the most likely number of genetic clusters with the ΔK method (A) and the probability logarithm of the data (lnP(D) (B). ΔK clearly suggests that K = 4 is the most likely number of clusters, while lnP(D|K) shows leveling off from K = 4 and the highest likelihood at K = 8.

https://doi.org/10.1371/journal.pone.0115160.s002

(PDF)

S1 Table.

Haplotype frequencies of the mtDNA control region in ten populations of Eurasian lynx Lynx lynx.

https://doi.org/10.1371/journal.pone.0115160.s003

(DOC)

S2 Table.

Pairwise differentiation between Eurasian lynx populations. FST values based on 12 microsatellites (above diagonal) and FST values based on the mtDNA control region (below diagonal). Nonsignificant values are given in italics.

https://doi.org/10.1371/journal.pone.0115160.s004

(DOC)

S3 Table.

Pairwise differentiation between Eurasian lynx populations. DEST (above diagonal) and RST (below diagonal) values based on 12 microsatellites.

https://doi.org/10.1371/journal.pone.0115160.s005

(DOC)

S4 Table.

Habitat type, climate and lynx diet composition of the ten lynx populations studied. * References used for data on lynx diet: Norway [84]; Finland: [85]; Estonia and Latvia: [86]; Belarus: [87]; Poland (KARPF): [88]; Poland (BPF): [27]; Carpathians: [89]; Russia (Kirov): [90]. No data were available for Lithuania – for the computations we arbitrarily assumed the same prey composition as in Latvia due to very close geographical proximity. The data were recalculated from original sources considering the relative share of only key prey that occurred in all populations.

https://doi.org/10.1371/journal.pone.0115160.s006

(DOC)

S1 Appendix.

List of Eurasian lynx genotypes based on 12 microsatellite loci.

https://doi.org/10.1371/journal.pone.0115160.s007

(XLS)

S2 Appendix.

Procedures adopted in the Approximate Bayesian Computations (DIYABC 2.0.4 software) and the list of all scenarios tested.

https://doi.org/10.1371/journal.pone.0115160.s008

(DOC)

Acknowledgments

This study would not be possible without a large scale sampling conducted selflessly by numerous people, including: E. Bujko, K. Bozik, Z. Bosko, B. Dawidziuk, R. Fedaczyński, M. Kurzawa, S. Nowak, M. Shkvyria, R. Unici, J. Wa˛siński and T. Werkowski as well as numerous hunters who collected samples during the national programs of large carnivore control in Norway, Finland, Latvia, Estonia and Russia. We are very grateful to Astrid V. Stronen and Allan McDevitt for careful reading the manuscript, providing useful comments and correcting the English. Thanks to K. O′Mahony for his editing the English of the final version. M. Górny prepared the maps. We are deeply greatful to two anonymous reviewers for their very constructive comments, which greatly improved the first version of this paper.

Author Contributions

Conceived and designed the experiments: MR KS RK. Performed the experiments: MR MM APS VS JO PM LB IK HO KS. Analyzed the data: MR MM KS. Contributed reagents/materials/analysis tools: APS VS JO PM LB IK HO KS RK. Wrote the paper: KS MR MM.

References

  1. 1. Culver M, Johnson WE, Pecon-Slattery J, O'Brien S (2000) Genomic ancestry of the American puma (Puma concolor). J Hered 91:186–197.
  2. 2. Schwartz MK, Mills LS, McKelvey KS, Ruggiero LF, Allendorf FW (2002) DNA reveals high dispersal synchronizing the population dynamics of Canada lynx. Nature 415:520–522.
  3. 3. Kyle CJ, Strobeck C (2002) Connectivity of peripheral and core populations of North American wolverines. J Mammal 83:1141–1150.
  4. 4. Anderson CR, Lindzey FG, McDonald DB (2004) Genetic structure of cougar populations across the Wyoming Basin: metapopulation or megapopulation. J Mammal 85:1207–1214.
  5. 5. Castilho CS, Marins-Sá LG, Benedet RC, Freitas TO (2011) Landscape genetics of mountain lions (Puma concolor) in southern Brazil. Mamm Biol 76:476–483.
  6. 6. Row JR, Gomez C, Koen EL, Bowman J, Murray DL, et al. (2012) Dispersal promotes high gene flow among Canada lynx populations across mainland North America. Conserv Genet 13:1259–1268.
  7. 7. Cegelski CC, Waits LP, Anderson NJ (2003) Assessing population structure and gene flow in Montana wolverines (Gulo gulo) using assignment-based approaches. Mol Ecol 12:2907–2918.
  8. 8. Stenseth NC, Shabbar A, Chan KS, Boutin S, Rueness EK, et al. (2004) Snow conditions may create an invisible barrier for lynx. Proc Natl Acad Sci USA 101:10632–10634.
  9. 9. McRae BH, Beier P, Dewald LE, Huynh LY, Keim P (2005) Habitat barriers limit gene flow and illuminate historical events in a wide-ranging carnivore, the American puma. Mol Ecol 14:1965–77.
  10. 10. Pilot M, Jędrzejewski W, Branicki W, Sidorovich VE, Jędrzejewska B, et al. (2006) Ecological factors influence population genetic structure of European grey wolves. Mol Ecol 15:533–4553.
  11. 11. Musiani M, Leonard JA, Cluff HD, Gates CC, Mariani S, et al. (2007) Differentiation of tundra/taiga and boreal coniferous forest wolves: genetics, coat colour and association with migratory caribou. Mol Ecol 16:4149–4170.
  12. 12. Schmidt K, Ratkiewicz M, Konopiński M (2011) The importance of genetic variability and population differentiation in the Eurasian lynx Lynx lynx for conservation in the context of habitat and climate change. Mammal Rev 41:112–124.
  13. 13. Holbrook JD, DeYoung RW, Janecka JE, Tewes ME, Honeycutt RL, et al. (2012) Genetic diversity, population structure, and movements of mountain lions (Puma concolor) in Texas. J Mammal 93:989–1000.
  14. 14. Czarnomska S, Jędrzejewska B, Borowik T, Niedziałkowska M, Stronen A, et al. (2013) Concordant mitochondrial and microsatellite DNA structuring between Polish lowland and Carpathian Mountain wolves. Conserv Genet 14:573–588.
  15. 15. Keis M, Remm J, Ho SYW, Davison J, Tammeleht E, et al. (2013) Complete mitochondrial genomes and a novel spatial genetic method reveal cryptic phylogeographical structure and migration patterns among brown bears in north-western Eurasia. J Biogeogr 40:915–927.
  16. 16. Stronen AV, Jędrzejewska B, Pertoldi C, Demontis D, Randi E, et al. (2013) North-south differentiation and a region of high diversity in European wolves (Canis lupus). PLoS ONE 8:e76454.
  17. 17. Valdiosera CE, García N, Anderung C, Dalén L, Crégut-Bonnoure E, et al. (2007) Staying out in the cold: glacial refugia and mitochondrial DNA phylogeography in ancient European brown bears. Mol Ecol 16:5140–5148.
  18. 18. Saarma U, Ho SYW, Pybus OG, Kaljuste M, Tumanov IL, et al. (2007) Mitogenetic structure of brown bears (Ursus arctos L.) in northeastern Europe and a new time frame for the formation of European brown bear lineages. Mol Ecol 16:401–413.
  19. 19. Vila C, Amorim IR, Leonard JA, Posada D, Castroviejo J, et al. (1999) Mitochondrial DNA phylogeography and population history of the grey wolf Canis lupus. Mol Ecol 8:2089–2103.
  20. 20. Pilot M, Branicki W, Jedrzejewski W, Goszczyński J, Jedrzejewska B, et al. (2010) Phylogeographic history of grey wolves in Europe. BMC Evol Biol 10:104.
  21. 21. Waits L, Taberlet P, Swenson JE, Sandegren F, Franzen R (2000) Nuclear DNA microsatellite analysis of genetic diversity and gene flow in the Scandinavian brown bear (Ursus arctos). Mol Ecol 9:421–431.
  22. 22. Taberlet P, Swenson JE, Sandegren F, Bjarvall A (1995) Localization of a contact zone between two highly divergent mitochondrial DNA lineages of the brown bear (Ursus arctos) in Scandanavia. Conserv Biol 9:1255–1261.
  23. 23. Wan Q-H, Wu H, Fujihara T, Fang SG (2004) Which genetic marker for which conservation genetics issue? Electrophoresis 25:2165–2176.
  24. 24. Hellborg L, Walker CW, Rueness EK, Stacy JE, Kojola I, et al. (2002) Differentiation and levels of genetic variation in northern European lynx (Lynx lynx) populations revealed by microsatellites and mitochondrial DNA analysis. Conserv Genet 3:97–111.
  25. 25. Rueness EK, Naidenko S, Trosvik P, Stenseth NC (2014) Large-scale genetic structuring of a widely distributed carnivore - the Eurasian lynx (Lynx lynx). PLoS ONE 9(4):e93675.
  26. 26. Ratkiewicz M, Matosiuk M, Kowalczyk R, Konopiński MK, Okarma H, et al. (2012) High levels of population differentiation in Eurasian lynx at the edge of the species' western range in Europe revealed by mitochondrial DNA analyses. Anim Conserv 15:603–612.
  27. 27. Jędrzejewski W, Schmidt K, Miłkowski L, Jędrzejewska B, Okarma H (1993) Foraging by lynx and its role in ungulate mortality: the local (Białowieża Forest) and the Palearctic viewpoints. Acta Theriol 38:385–403.
  28. 28. Schmidt K, Kowalczyk R (2006) Using scent-marking stations to collect hair samples to monitor Eurasian lynx populations. Wildlife Soc B 34:462–466.
  29. 29. Menotti-Raymond M, David VA, Lyons LA, Schafer AA, Tomlin JF, et al. (1999) A genetic linkage map of microsatellites in the domestic cat (Felis catus). Genomics 57:9–23.
  30. 30. Luo S-J, Johnson WE, David VA, Menotti-Raymond M, Stanyon R, et al. (2007) Development of Y chromosome intraspecific polymorphic markers in the felidae. J Hered 98:400–413.
  31. 31. Schmidt K, Kowalczyk R, Ozolins J, Männil P, Fickel J (2009) Genetic structure of the Eurasian lynx population in north-eastern Poland and the Baltic states. Conserv Genet 10:497–501.
  32. 32. Davoli F, Schmidt K, Kowalczyk R, Randi E (2012) Hair snaring and molecular genetic identification for reconstructing the spatial structure of Eurasian lynx populations. Mamm Biol 78:118–126.
  33. 33. Van Oosterhout C, Hutchinson WF, Wills DPM, Shipley P (2004) MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Mol Ecol Notes 4:535–538.
  34. 34. Kalinowski ST, Taper ML, Marshall TC (2007) Revising how the computer program CERVUS accommodates genotyping error increases success in paternity assignment. Mol Ecol 16:1099–1106.
  35. 35. Hundertmark KJ, Shields GF, Udina IG, Bowyer RT, Danilkin AA, et al. (2002) Mitochondrial phylogeography of moose (Alces alces): late Pleistocene divergence and population expansion. Mol Phylogenet Evol 22:375–387.
  36. 36. Hellborg L, Ellegren H (2003) Y chromosome conserved anchored tagged sequences (YCATS) for the analysis of mammalian male specific DNA. Mol Ecol 12:283–291.
  37. 37. Goudet J (1995) FSTAT (Version 1.2). A Computer Program to Calculate F-Statistics. J Hered 86:485–486.
  38. 38. Rousset F (2008) Genepop'007: a complete reimplementation of the Genepop software for Windows and Linux. Mol Ecol Resour 8:103–106.
  39. 39. Schneider S, Roessli D, Excoffier L (2000) Arlequin: a software for population genetics data analysis. User manual ver 2.000. Geneva: Genetics and Biometry Lab, Dept. of Anthropology, University of Geneva. 111 p.
  40. 40. Clement M, Posada D, Crandall K (2000) TCS: a computer program to estimate gene genealogies. Mol Ecol 9:657–1660.
  41. 41. Jost L (2008) GST and its relatives do not measure differentiation. Mol Ecol 17:4015–4026.
  42. 42. Crawford NG (2010) SMOGD: software for the measurement of genetic diversity. Mol Ecol Resour 10:556–557.
  43. 43. Rousset F (1997) Genetic differentiation and estimation of gene flow from F-statistics under isolation by distance. Genetics 145:1219–1228.
  44. 44. Jensen JL, Bohonak AJ, Kelley ST (2005) Isolation by distance, web service. BMC Genetics 6:13.
  45. 45. Peakall R, Smouse PE (2012) GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research—an update. Bioinformatics 28:2537–2539.
  46. 46. Excoffier L, Smouse PE, Quattro JM (1992) Analysis of molecular variance inferred from metric distances among DNA haplotypes – application to human mitochondrial DNA restriction data. Genetics 131:479–491.
  47. 47. Dupanloup I, Schneider S, Excoffier L (2002) A simulated annealing approach to define the genetic structure of populations. Mol Ecol 11:2571–2581.
  48. 48. Hardy OJ, Vekemans X (2002) SPAGeDi: a versatile computer program to analyse spatial genetic structure at the individual or population levels. Mol Ecol Notes 2:618–620.
  49. 49. Wilson GA, Rannala B (2003) Bayesian inference of recent migration rates using multilocus genotypes. Genetics 163:1177–1191.
  50. 50. Pritchard JK, Stephens M, Donnelly P (2000) Inference of population structure using multilocus genotype data. Genetics 155:945–959.
  51. 51. Evanno G, Regnaut S, Goudet J (2005) Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol 14:2611–2620.
  52. 52. Earl DA, von Holdt BM (2012) STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv Genet Resour 4:359–361.
  53. 53. Cornuet J-M, Pudlo P, Veyssier J, Dehne-Garcia A, Gautier M, et al. (2014) DIYABC v2.0: a software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data. Bioinformatics 30:1187–1189.
  54. 54. Barnett R, Yamaguchi N, Shapiro B, Ho S, Barnes I, et al. (2014) Revealing the maternal demographic history of Panthera leo using ancient DNA and a spatially explicit genealogical analysis. BMC Evol Biol 14:70.
  55. 55. Smouse PE, Peakall R (1999) Spatial autocorrelation analysis of individual multiallele and multilocus genetic structure. Heredity 82:561–573.
  56. 56. Anderson MJ (2004) DISTLM v.5: a FORTRAN computer program to calculate a distance-based multivariate analysis for a linear model. New Zealand: Department of Statistics, University of Auckland.10 p.
  57. 57. Bednorz E (2004) Snow cover in eastern Europe in relation to temperature, precipitation and circulation. Int J Climatol 24:591–601.
  58. 58. Hurrell JW (1995) Decadal trends in the North Atlantic Oscillation and relationships to regional temperature and precipitation. Science 269:676–679.
  59. 59. Anderson MJ (2003) DISTLM forward: a FORTRAN computer program to calculate a distance-based multivariate analysis for a linear model using forward selection. New Zealand: Department of Statistics, University of Auckland. 10 p.
  60. 60. Lorenc H, Editor 2005) Atlas of climate in Poland. Warszawa: IMiGW. 116 p.
  61. 61. Waser PM, Jones WT (1983) Natal philopatry among solitary mammals. Q Rev Biol 58:355–390.
  62. 62. Smith JLD (1993) The Role of Dispersal in Structuring the Chitwan Tiger Population. Behaviour 124:165–195.
  63. 63. Biek R, Akamine N, Schwartz MK, Ruth TK, Murphy KM, et al. (2006) Genetic consequences of sex-biased dispersal in a solitary carnivore: Yellowstone cougars. Ecol Lett 2:312–315.
  64. 64. Schmidt K (1998) Maternal behaviour and juvenile dispersal in the Eurasian lynx. Acta Theriol 43:391–408.
  65. 65. Samelius G, Andrén H, Liberg O, Linnell JDC, Odden J, et al. (2012) Spatial and temporal variation in natal dispersal by Eurasian lynx in Scandinavia. J Zool 286:120–130.
  66. 66. Häme T, Stenberg P, Andersson K, Rauste Y, Kennedy P, et al. (2001) AVHRR-based forest proportion map of the Pan-European area. Remote Sens Environ 77:76–91.
  67. 67. Niedziałkowska M, Jędrzejewski W, Mysłajek RW, Nowak S, Jędrzejewska B, et al. (2006) Environmental correlates of Eurasian lynx occurrence in Poland – Large scale census and GIS mapping. Biol Conserv 133:63–69.
  68. 68. Podgórski T, Schmidt K, Kowalczyk R, Gulczyńska A (2008) Microhabitat selection by Eurasian lynx and its implications for species conservation. Acta Theriol 53:97–110.
  69. 69. Huck M, Jędrzejewski W, Borowik T, Jędrzejewska B, Nowak S, et al. (2011) Analyses of least cost paths for determining effects of habitat types on landscape permeability: wolves in Poland. Acta Theriol 56:91–101.
  70. 70. Kojola I, Aspi J, Hakala A, Heikkinen S, Ilmoni C, et al. (2006) Dispersal in an expanding wolf population in Finland. J Mammal 87:281–286.
  71. 71. Nilsen EB, Linnell JDC, Odden J, Andersen R (2009) Climate, season, and social status modulate the functional response of an efficient stalking predator: the Eurasian lynx. J Anim Ecol 78:741–751.
  72. 72. Yom-Tov Y, Kvam T, Wiig Ø (2011)Lynx body size in Norway is related to its main prey (roe deer) density, climate, and latitude. AMBIO 40:43–51.
  73. 73. Sunquist ME, Sunquist FC (2002) Wild Cats of the World. Chicago and London: The University of Chicago Press. 452 p.
  74. 74. Molinari-Jobin A, Molinari P, Breitenmoser-Würsten Ch, Breitenmoser U (2002) Significance of lynx Lynx lynx predation for roe deer Capreolus capreolus and chamois Rupicapra rupicapra mortality in the Swiss Jura Mountains. Wildlife Biol 8:109–115.
  75. 75. Okarma H, Jędrzejewski W, Schmidt K, Kowalczyk R, Jędrzejewska B (1997) Predation of Eurasian lynx on roe deer and red deer in Białowieża Primeval Forest, Poland. Acta Theriol 42:203–224.
  76. 76. Hénaux V, Powell LA, Hobson KA, Nielsen CK, LaRue MA (2011) Tracking large carnivore dispersal using isotopic clues in claws: an application to cougars across the Great Plains. Method Ecol Evol 2:489–499.
  77. 77. Fernández N, Selva N, Yuste C, Okarma H, Jakubiec Z (2012) Brown bears at the edge: Modeling habitat constrains at the periphery of the Carpathian population. Biol Conserv 153:134–142.
  78. 78. Schwartz MK, Mills LS, Ortega Y, Ruggiero LF, Allendorf FW (2003) Landscape location affects genetic variation of Canada lynx (Lynx canadensis). Mol Ecol 12:1807–1816.
  79. 79. Linnell JC, Broseth H, Odden J, Nilsen E (2010) Sustainably harvesting a large carnivore? Development of Eurasian lynx populations in Norway during 160 years of shifting policy. Environ Manage 45:1142–1154.
  80. 80. Gugolz D, Bernasconi MV, Breitenmoser-Würsten C, Wandeler P (2008) Historical DNA reveals the phylogenetic position of the extinct Alpine lynx. J Zool 275:201–208.
  81. 81. Sommer RS, Nadachowski A (2006) Glacial refugia of mammals in Europe evidence from fossil records. Mammal Rev 36:251–265.
  82. 82. Anderson GS, Danielson BJ (1997) The effect of landscape composition and physiognomy on metapopulation size: the role of corridors. Landscape Ecol 12:261–271.
  83. 83. Jędrzejewski W, Jędrzejewska B, Okarma H, Schmidt K, Bunevich AN, et al. 1996) Population dynamics (1869–1994), demography and home ranges of the lynx in Białowieża Primeval Forest (Poland and Belaruss). Ecography 19:122–138.
  84. 84. Sunde P, Kvam T (1997) Diet patterns of Eurasian lynx Lynx lynx: what causes sexually determined prey size segregation? Acta Theriol 42:189–201 Cited in S4 Table.
  85. 85. Pulliainen E, Lindgren E, Tunkkari PS (1995) Influence of food availability and reproductive status on the diet and body condition of the European lynxes in Finland. Acta Theriol 40:181–196 Cited in S4 Table..
  86. 86. Valdmann H, Andersone-Lilley Z, Koppa O, Ozolins J, Bagrade G (2005) Winter diets of wolf Canis lupus and lynx Lynx lynx in Estonia and Latvia. Acta Theriol 50:521–527 Cited in S4 Table.
  87. 87. Sidorovich VE (2006) Relationship between prey availability and population dynamics of the Eurasian lynx and its diet in northern Belarus. Acta Theriol 51:265–274 Cited in S4 Table.
  88. 88. Dojlida M (2006) Distribution, number and food of lynx in the Knyszyńska Forest. MSc thesis. Warszawa: Warszawa University of Life Sciences. [in Polish] Cited in S4 Table.
  89. 89. Hell P, Slamecka J (1996) Current status of the lynx (Lynx lynx) in Slovakia. Lynx in the Czech and Slovak Republics. Acta Sci Nat Brno 30:64–78 Cited in S4 Table.
  90. 90. Kozlovskiy IS (2003) The Vyatka-Kama Region. In: Matyushkin EN, Vaisfeld MA, editors. The Lynx. Regional features of ecology, use and protection. Moscow: Nauka. pp.168–217 Cited in S4 Table.
  91. 91. Von Arx M, Breitenmoser-Würsten C, Zimmermann F, Breitenmoser U (2004) Status and conservation of the Eurasian lynx (Lynx lynx) in 2001. Muri: KORA Bericht no 19:330.