Long-distance Eurasian lynx dispersal – a prospect for connecting native and reintroduced populations in Central Europe

Dispersal is a key process for the maintenance of intraspecific genetic diversity by ensuring gene flow within and between populations. Despite the ongoing expansion of large carnivores in Europe, lynx populations remain fragmented, isolated, and threatened by inbreeding and loss of genetic diversity. In the course of large carnivore monitoring in the Czech Republic, several biological samples of Eurasian lynx were collected outside the permanent occurrence of this species. Using microsatellite genotyping we identified these as four dispersing lynx males and applied multiple methods (Bayesian clustering in STRUCTURE, Principal Component Analysis (PCA), frequency-based method in GENECLASS2, and machine-learning framework in assignPOP) to assign them to possible source populations. For this we used genotypes from five European lynx populations: the Bohemian-Bavarian-Austrian (N = 36), Carpathian (N = 43), Scandinavian (N = 20), Baltic (N = 15), and Harz (N = 23) population. All four dispersers were successfully assigned to different source populations within Europe and each was recorded at a distance of more than 98 km from the edge of the distribution of the source population identified. Such movements are among the longest described for lynx in Central Europe to this point. The findings indicate the ability of lynx males to disperse in human-dominated landscape thus facilitation of these movements via creation and/or protection of potential migratory corridors together with protection of dispersing individuals should be of high importance in conservation of this iconic predator in Central Europe.


Introduction
Dispersal has been recognised as a key process in the dynamics and evolution of natural populations. Through the redistribution of individuals, dispersal is the main factor leading to gene flow within and between populations and thus maintenance of genetic diversity (Bullock et al. 2002;Clobert et al. 2012). However, in human-dominated landscape, populations are often segregated to isolated patches surrounded by an unsuitable matrix, where movement between the patches can be challenging (Fahrig 2003;Seidler et al. 2015;Crooks et al. 2017).
Species with low densities but extensive home-ranges, such as large carnivores, are substantially vulnerable to anthropogenic changes in the landscape. A human-dominated landscape may pose a serious obstacle in particular for long-distance dispersers, reducing the frequency of such dispersal movements. That said, these long-distance dispersal events may often remain undetected, as they are difficult to record, particularly in elusive species (Trakhtenbrot et al. 2005).
Despite the ongoing expansion in numbers and distribution of large carnivores in Europe (Chapron et al. 2014), only a few long-distance movements have been recorded (Bartoń et al. 2019). Grey wolves Canis lupus, were recorded to disperse up to 1092 km (Wabakken et al. 2007) and brown bears Ursus arctos, up to 467 km (Støen et al. 2006). Compared to bear or wolf, the Eurasian lynx Lynx lynx, has been described as a conservative disperser (Zimmermann et al. 2005), although long-distance movements have been recorded, with maximum up to 428 km in Scandinavia (Samelius et al. 2012) and up to 200 km in the Alps (Groff et al. 2009). While the landscapes of the Alps and Scandinavia are comparably remote, very little is known about the ability of lynx to disperse within human-altered landscape (Schmidt 1998;Zimmermann et al. 2005Zimmermann et al. , 2007Samelius et al. 2012;Herrero et al. 2020). Generally, as many other solitary felids (Janečka et al. 2007;Gour et al. 2013;Fattebert et al. 2015;Wultsch et al. 2016), lynx seem to follow a male-biased dispersal pattern (Samelius et al. 2012;Krojerová-Prokešová et al. 2019), with males dispersing far more widely than females from a given source population, although this has not been unequivocally proven in all areas (Zimmermann et al. 2005Herrero et al. 2020).
In Central Europe native populations of lynx occur only in the Western Carpathians and in north-eastern Poland, on the edge of the distribution range of the Baltic population. Besides these native populations, several populations have been established in central, western and southern Europe based on reintroduction of lynx mainly of West-Carpathian origin (Dinaric, Bohemian-Bavarian-Austrian (BBA), Jura, Alpine, and Vosges-Palatinian populations; Kaczensky et al. 2013;Fig. 1). The population in the Harz Mountains in Germany and that in the Kampinoski National Park in central Poland were both based on the release of lynx of captive origin. Although reintroduction in Central Europe has already started in the 1970s (Kaczensky et al. 2013), established populations are mostly isolated and of small size (fewer than 200 or even less than 100 individuals; Chapron et al. 2014). Further, in each case, the number of released animals was low and in some cases the animals were even closely related. This has inevitably resulted in low levels of genetic diversity and an increased risk of inbreeding depression (Breitenmoser-Würsten and Obexer-Ruff 2003;Sindičić et al. 2013;Bull et al. 2016;Mueller et al. 2020). However, it has been shown that even single dispersing individuals may considerably enhance the genepool and population viability in small populations (Frankham et al. 2002;Bull et al. 2016;Åkesson et al. 2016).
In this study we thus focused on the detection of potential long-distance dispersers of Eurasian lynx within a Central European region, namely the Czech Republic. The country is occupied by the native Carpathian population in the north-east and the reintroduced BBA population in the south-west. The majority of the Czech landscape is highly urbanised with large open areas of agricultural land, dissected by busy motorways. The forested mountain areas, characteristically associated with species such as lynx, are located mainly around the borders. During extensive monitoring of large carnivore populations, biological samples of Eurasian lynx were collected outside the known areas of permanent occurrence of the lynx in the Czech Republic. Knowledge about the origin of these lynx individuals is important for understanding movement patterns within such a fragmented landscape towards reviewing the future conservation strategy for lynx in the entire Central European region as well as for planning for improvements in landscape permeability using mitigation measures like green bridges or underpasses/overpasses. The aim of this study was thus to assign the lynx samples of unknown origin to their possible source populations. We genotyped individuals recorded outside the main distribution range and lynx from potential source populations. We then applied a variety of assignment methods to determine the origin of the focal lynx individuals. Based on our findings, we discuss the conservation implications and potential for long-term natural connectivity of native and reintroduced lynx populations in a single Central European metapopulation.

Sampling and genotyping
Evidence for the presence of lynx individuals of unknown origin ('dispersers') were found outside the permanent distribution range of the lynx in the Czech Republic ( Fig. 1; Table 1). Five non-invasive scat samples (two in Jizera Mts, two in the Moravian Karst, one in Ore Mts) were collected during extensive monitoring of large carnivores, especially wolf. One blood sample was collected during immobilisation of lynx in the Moravian Karst and one tissue sample was taken from the carcass found after traffic collision on the D1 motorway. These samples were genotyped as described below.
Tissue, blood, and scat samples were stored in 96% ethanol. DNA from tissue and blood samples was extracted using the Genomic DNA Mini kit Tissue (Geneaid Biotech Ltd., New Taipei City, Taiwan). The QIAamp DNA Stool Mini Europe is delineated according to Chapron et al. (2014) and for the Czech Republic updated according to Kutal et al. (2017). Dashed arrows show the possible course of the dispersal based on the results of genetic assignments. As a background we used the forest type cover layer from © Copernicus Land Monitoring Service 2018, European Environment Agency (EEA), European Union Kit (Qiagen GmbH, Hilden, Germany) was used to extract DNA from faecal samples. Dedicated laboratories for DNA extraction and PCR setup were used and we followed strict rules and procedures to prevent contamination. In total, 15 microsatellite loci and the amelogenin marker for sex determination were used. PCR amplification of non-invasive samples was repeated in independent multiple parallels according to Taberlet et al. (1996) and Adams and Waits (2007) with requirements of a minimum of three positive amplifications for homozygotes and two for heterozygotes. We followed the same protocols for isolation, genotyping, PCR conditions and fragment analysis as are given in Krojerová-Prokešová et al. (2019).
To be able to assign the origin of possible dispersers, a reference dataset of 137 unique individual genotypes from five European populations ( Fig. 1; Table S1; Table S2) was used. The reference samples were genotyped on the same microsatellite set of loci following the same protocols as  (77) were taken from individuals found dead, poached, or legally hunted. Non-invasive samples (60; 42 scat, 17 hair, 1 blood samples) were collected during regular monitoring in the Bohemian Forest, the Carpathians, and the Harz Mountains (Table S1).

Analysis of population structure and individual assignment
Differentiation between comparative populations, necessary for correct assignment of dispersers, was quantified by pairwise F ST based on Weir and Cockerham (1984) in the hierfstat package (Goudet 2005) in R software (R Core Team 2013). Confidence intervals of F ST values (97.5% CI) were estimated using 999 bootstrap replicates in the same R package.
To infer the structure of populations and origin of dispersers, we used the Bayesian clustering analysis implemented in STRU CTU RE v2.3.4 (Pritchard et al. 2000). The program was run with 10 independent runs for each value of K from 1 to 10, with 1 000 000 MCMC iterations and initial burn-in of 100 000. An admixture ancestry model with a correlated allele frequency model, without using sampling locations as prior information, was applied. The most optimal number of genetically distinct clusters was estimated using the ΔK method of Evanno et al. (2005) and evaluated also with the aid of estimators accounting for uneven sampling and hierarchical structure (Puechmaille 2016) in online application StructureSelector (Li and Liu 2018). The same application, integrating Clumpak program (Kopelman et al. 2015), was used to generate graphical representations for specific K.
As an alternative approach for investigation of population divergence and assignment of dispersers to potential source populations we performed a Principal Component Analysis (PCA) using package adegenet v.2.0.0 (Jombart 2008) in R software.
Further, solely for individual assignment we applied the frequency-based method of Paetkau et al. (1995) within the program GENECLASS2 (Piry et al. 2004). Genotypes of all individuals from reference populations were used within one file while genotypes of dispersers were inputted separately. Probability of assignment was performed by simulating 100 000 individuals with Monte Carlo resampling method ) and setting the type I errors to 0.05 (Piry et al. 2004).
Finally, to investigate the origin of the dispersers we applied a machine-learning approach in the R package assignPOP v1.1.4 (Chen et al. 2018). The program assign-POP was developed to overcome issues associated with nonindependence and imbalance of datasets (Chen et al. 2018). The assignPOP approach included data evaluation, where all individuals from reference populations were randomly divided into training and test sets and the assignment accuracies were estimated via Monte-Carlo cross-validation based on the following parameters: proportion of individuals used in training set: 0.5, 0.7, and 0.9; proportion of loci used in training set: 0.1, 0.25, 0.5, and 1; loci sample method: F ST ; iterations: 1000; and model: linear discriminant analysis (LDA). Simulations performed best when all loci and individuals were used (Fig. S1). Therefore, for assignment test we used all individuals from reference populations and all microsatellite loci with chosen model LDA.

Microsatellite genotyping
DNA isolation and subsequent genotyping was successful for all seven samples collected outside the 'permanent' distribution of lynx in the Czech Republic. Within these samples we identified four lynx males (Table 1; Fig. 1).

Structure of populations
The pairwise F ST values confirmed differentiation among all possible source populations sampled (Table 2) with the lowest difference apparent between the BBA and Carpathian populations (F ST = 0.13) and highest value between the BBA and Harz populations (F ST = 0.40).
The Bayesian clustering analysis detected the most probable number of clusters at K = 2 and then at K = 5 (Fig. S2a). The best model for K = 5 was supported also by Puechmaille approach (Fig. S2b). Structure of two clusters was in concordance with delineation of two lynx subspecies: the Northern lynx (Lynx lynx lynx) and the Carpathian lynx (Lynx lynx carpathicus). Further clustering at K = 5 corresponded to presumed structure of the five populations sampled (Fig. 2).
PCA analysis likewise clearly distinguished the five possible source populations (Fig. 3), although the differentiation between the Carpathian and the BBA populations (itself founded by lynx from the Carpathians) became clearer when PCA was repeated using only the samples from these two populations (Fig. S3).

Assignment of dispersers
Lynx L3 and L4 were unequivocally assigned to the same probable population of origin by all four methodological approaches -L3 to the Carpathian population (confirmed also by parentage analysis, Appendix 1) and L4 to the Harz population. L1 was classified as of BBA origin except for GENECLASS2 approach, where it had the highest probability to be of Carpathian origin (0.84); an origin from the BBA population however also received relatively high support (0.52). L2 was identified as Baltic lynx by majority  Table 3). All individuals were sampled at a distance of more than 98 km (98-456 km) from the border of the distribution range of the identified source population, indicating long-distance dispersal (Table 3).

Population origin of dispersers
Except for two cases, all four genetic approaches provided consistent assignment for all dispersers. In the first case, GENECLASS2 software was unable to distinguish between BBA and Carpathian origin of lynx L1. In the second case, program assignPOP assigned individual L2 both to Baltic and Harz population with almost the same probability. Differentiation between these two pairs of populations is in any case the most problematic (Table 2) as the BBA population originated from the Carpathian population and the Harz population was based on the release of captive lynx of mixed origin (von Arx et al. 2009), potentially including also ancestors from the Baltic region.
According to the majority of approaches, male L1, found dead on the D1 motorway in Central Bohemia, originated from the BBA population (103 km straight line distance from the BBA distribution range) from where it is presumed he was unsuccessfully trying to disperse. D1 motorway (the Fig. 4 Probabilities of individual assignment test in assignPOP for all four dispersers (L1-L4) Table 3 Comparison of four approaches for the assignment of dispersers to their putative population of origin *Estimated minimum distance is the distance between sampling site and the border of the distribution range of the presumed source population measured in ArcGIS 9.3.1 as a straight-line distance to the centroid of the nearest EEA square with continuous lynx occurrence. In the case of L2 a distance of 337 km corresponds to the population in Kampinoski NP, 456 km to the Baltic population and 283 km to the Harz population. For L3, a distance of 122 km corresponds to the distance to the nearest EEA square of its maternal home range (Appendix 1) main motorway of the Czech Republic) is known to represent an important migration barrier for large mammals (Anděl et al. 2010). During the seventeenth and eighteenth century lynx were believed to be exterminated in Jizera Mts and in adjacent Krkonoše Mts, but single individuals were spotted there again in 2002 and 2006, respectively (Flousek et al. 2014). Similarly, on the other side of the border with Poland, single lynx have been recorded between 2000 and 2018 (Mysłajek et al. 2019). Due to the mountain wildlife corridor with suitable habitat from the Western Carpathians to this area, it was believed that these lynx are dispersers from the Carpathians (Kratochvíl and Vala 1968). Similarly, lynx L2 sampled in Jizera Mts in the northern part of the Czech Republic was expected to be of Carpathian origin. Surprisingly, our genetic data do not support this presumption as lynx L2 was classified as a Baltic lynx. The edge of the native Baltic population in the north-eastern Poland is 456 km far away from the sampling site of the L2 and to reach the Jizera Mts, the individual would have had to cross highly urbanised central Poland with disrupted migration corridors and many barriers. As an alternative, it is perhaps possible that the smaller population occurring in Kampinoski National Park in central Poland (337 km straight line distance), could have been the possible source of this lynx. This population was founded between 1993 and 2000 (Böer et al. 2000), when 31 individuals from zoos and wildlife parks in Germany, Sweden and Finland were released there. Unfortunately, more detailed information about their ancestry is missing (Böer et al. 2000), thus we can only speculate if lynx of the Baltic origin may have been released there. An origin of L2 from within the Kampinoski National Park may, however, be supported by observations of the high dispersal ability of the original lynx released in that area. Soon after their release some of the lynx dispersed between 10-50 km to adjacent nature reserves (Böer et al. 2000), single individuals even up to 160 km (Reklewski 2006).
There is also a possibility that lynx L2 could come from the Harz population, as was indicated in assignPOP, but this is not the most parsimonious explanation and 3/4 of used approaches do not support this hypothesis. The problem with the correct assignment of L2 in assignPOP may also be influenced by the reference samples used for the Baltic population. We used 15 samples from Latvia, but recent genetic analyses revealed slight differentiation of the Baltic population in Poland to that in the central part of the distribution range (Ratkiewicz et al. 2014;Lucena-Perez et al. 2020).
Lynx L3 in the Moravian Karst was firstly recorded in September 2016. After that the lynx was captured there, fitted with a GPS collar and radio-tracked from June 2017 till June 2018 (Duľa and Krofel 2020). During autumn 2018, after losing his GPS collar, he disappeared from all monitoring sites established in the whole area, where he was previously regularly observed using camera traps (Duľa, pers. comm.). Our results showed that lynx L3 was of Carpathian origin (minimum dispersal distance 98 km), who, after dispersal, had settled in an area halfway between the Carpathian and BBA populations, suggesting this area might constitute a possible stepping-stone population between these two isolated populations. The results of the genetic assignment were further supported by parentage analysis (Appendix 1), which identified the parent pair of L3 in the Moravian-Silesian Beskids in the Western Carpathians. The pair of lynx identified as the parents is one of a few pairs who gave birth to the majority of juveniles in this area during 2011) and the centre of its maternal home range was in the straight line distance of 122 km from the site in the Moravian Karst where L3 had been recorded.
There are just a few documented records of lynx since 2017 in Ore Mountains, where the non-invasive sample from L4 was collected. It was initially expected the lynx would have come from the BBA population along the forested borderland between the Czech Republic and Germany, but genetic assignment suggested its origin in the Harz Mts in Germany (169 km distant). This movement shows the potential for gene flow between the BBA and the Harz population, both isolated populations with low (Bull et al. 2016), or declining levels of genetic diversity (Mueller et al. 2020). Admixture between different subspecies could pose a risk because it can lead to low or in an extreme case zero fitness of progeny through loss of local adaptations (Lynch 1991). On the other hand, fitness can be enhanced due to new genetic combinations and increased adaptive potential (i.e., heterosis), especially in inbred populations (Anderson 1949;Abbott et al. 2013;Frankham 2015). In any case, the effect of potential future crossbreeding of the expanding Harz population and neighbouring BBA population, should be closely monitored, since of all the populations included in this study these are the most genetically distinct from one another (F ST = 0.4).

Factors affecting dispersal
Assuming the origins of the four lynx individuals sampled in the Czech Republic are accurately identified, the distances travelled from the population of presumed origin (98-456 km) are among the longest described, so far, for lynx dispersal events in Central Europe. For Jura and Alpine populations, the mean dispersal length of subadult lynx is reported at 26-63 km (max. 129 km; Zimmermann et al. 2005). Long-distance dispersals (up to 124 km) are known from the central part of the Baltic population (Latvia, Estonia). The longest dispersal distances, with mean 47-148 km (up to 428 km), have been recorded in Scandinavia, where prey density is lower and home ranges are 1 3 larger in comparison to other European populations (Samelius et al. 2012).
All four dispersers were males, which is consistent with the general expectation of male-biased dispersal in polygynous mammals (Greenwood 1980;Dobson 1982) and has previously been reported in the Eurasian lynx (Samelius et al. 2012;Krojerová-Prokešová et al. 2019) as well as in other felids, e.g., leopard (Fattebert et al. 2015), jaguar (Wultsch et al. 2016), bobcat (Janečka et al. 2007), and tiger (Gour et al. 2013). Female lynx seem to disperse less frequently and for shorter distances than males, even though male-biased dispersal was not confirmed in some areas (Zimmermann et al. 2005Herrero et al. 2020). Except for Northern Europe (Fennoscandia, Baltic states), where dispersal distances are generally longer (Samelius et al. 2012;Bagrade et al. 2016;Herrero et al. 2020), there seem to be no published records about female dispersal longer than 100 km (see more details given in the supplementary material of the review of Bartoń et al. (2019)). The reluctance of female lynx to disperse over longer distances, particularly in fragmented landscapes, thus reduces the species ability to colonise new areas (Port et al. 2020).
Reported expansion of the distributional range of large carnivores in Europe relates primarily to distribution of wolves and bears (Chapron et al. 2014); by comparison, lynx populations are mostly stable or decreasing and only three populations are slowly increasing in numbers and range (Alpine, Jura and Harz population; Large Carnivore Initiative for Europe 2020). Except for the recent sporadic occurrence of a lynx in the military area Libavá (Kutal et al. 2017), lynx L3 represents the only well-documented expansion of lynx westwards from the Carpathians despite the presence of suitable habitats in neighbouring areas, e.g. mountains along the Czech-Polish border -Jeseníky Mts, the Orlické Mts, Krkonoše Mts and Jizera Mts. Human intervention (traffic mortality, poaching) at the population periphery probably plays a considerable role in limiting the West-Carpathian lynx expansion (Krojerová-Prokešová et al. 2019).

Conservation implications
To facilitate natural movement between lynx populations and to maintain the viability of these populations, there is an urgent need to ensure landscape permeability via creation and/or protection of potential migratory corridors together with protection of dispersing individuals. Establishment of transboundary conservation strategies including appropriate population monitoring focusing on genetic diversity, inbreeding status and demography of the populations will help to form and protect exchange of individuals between the various European lynx populations. Source populations should be large enough and in favourable status to supply enough dispersers to facilitate gene flow; it has been shown, that if there is a contact with a larger population and there are suitable vacant territories and prey availability, the population can successfully expand and retain gene flow even between the areas with potential barriers to movement (Chapron et al. 2014;Bagrade et al. 2016).
The human-mediated translocation of animals between the populations and founding of new small or mediumsized populations in suitable areas in Central Europe, which may act as stepping-stone populations, is another option in reconnecting isolated lynx populations in different areas into one viable metapopulation. Moreover, as suggested by Port et al. (2020), consideration should be given to the translocation of a few females into areas accessible by male dispersers, or regularly visited by them, as a starting point for the development of new steppingstone populations, as males will be attracted into the area and they will have tendency to stay there. However, reintroduction of lynx populations into new areas should be well-planned (selection of suitable founders with regard to their population origin and genetic status, habitat quality and landscape connectivity in the surroundings) and used very cautiously due to a risk associated with capture and release of animals, high mortality of lynx exploring new areas after release and the conflicts arising with the acceptance of lynx by local key stakeholders (Červený et al. 2019), whose perception highly affects the success of all these conservation efforts.
Author contributions BG wrote the first draft of the manuscript and performed statistical analyses. BG and JKP performed laboratory analyses and designed the study. All authors collected samples and edited/ approved previous versions of the manuscript.
Funding This article is based upon work from COST Action G-Bike (CA18134), supported by COST (European Cooperation in Science and Technology). The study was financially supported by INTER-EXCELLENCE -INTER-COST (LTC20021), Interreg V-A SR-CR (304021D016), Interreg Central Europe (CE1001; 3Lynx) and by Institutional Research Plan (RVO: 68081766).