Revealing the origin of wildcat reappearance after presumed long-term absence

Following severe population decline and local extinction due to massive habitat destruction and persecution, wildcats have recently reappeared in several parts of Germany’s low mountain region. It remains unknown how this reemergence occurred, specifically if local populations have been overlooked at low densities or if the species has successfully spread across the highly fragmented anthropogenic landscape. In the central German Rhön Mountains, for instance, wildcats were believed to be extinct during most of the twentieth century, however, the species was recently detected and subsequent genetic monitoring found the presence of a sizeable population. In this study, we used microsatellite and SNP genotypes from 146 wildcat individuals from 2008 to 2017 across a ~ 15,000 km2 area in the central German low mountain region to understand the population re-establishment of wildcats in the region. Bayesian clustering and subsequent analyses revealed that animals in the Rhön Mountains appear to be a mix from the two adjacent populations in the North and South of the area, suggesting a recent range expansion from two different directions. Both populations meet in the Rhön Biosphere Reserve, leading to an admixture of the northern, autochthonous, and the southern reintroduced wildcat population. While we cannot completely exclude the possibility of undetected population persistence, the high genetic homogeneity in the central German wildcat population and the lack of any signatures of past population decline in the Rhön favor a scenario of natural expansion. Our findings thus suggest that wildcats are well capable of rapid range expansion across richly structured landscape mosaics consisting of open land, settlements, and forest patches and document the potential of massive non-invasive genetic sampling when aiming to reconstruct the complex population and range dynamics of wildlife.


Introduction
Many populations of large-and medium-sized carnivores are currently re-expanding their ranges across the densely populated and anthropogenically modified landscapes of central Europe (Chapron et al. 2014). These altered landscapes exhibit high levels of habitat fragmentation, which is the disruption of continuous stretches of suitable habitat (Schadt et al. 2002). Habitat fragmentation can negatively affect the viability of populations and species. However, the extent to which species are affected varies significantly (Haddad et al. 2015). A species' reaction to cultural landscapes depends heavily on the size of suitable habitat, the distance between habitat patches, and the individual species' life history traits (Caruso et al. 2016). Carnivores, for instance, are at a heightened risk of being negatively impacted in fragmented landscapes due to their relatively large home ranges and low population densities (Crooks 2002;Noss et al. 1996;Woodroffe 1998). Therefore, investigating how carnivores successfully disperse and establish viable populations in altered landscapes is important to understand the potential for long-term existence of wildlife in human-dominated landscapes (Riley et al. 2003).
The European wildcat is a flagship species for nature conservation in Germany, as it primarily relies on highly structured natural deciduous and mixed forests (Driscoll and Nowell 2010). This reliance on connected forests has made it a pillar for the preservation of natural habitats. Recently, the wildcat has made a reemergence across the central German low mountain region (Steyer et al. 2016), providing a relevant example to understand dispersal dynamics in a highly fragmented landscape.
Historically, European wildcat populations were heavily depleted due to human persecution until the early twentieth century. In Germany, the species was restricted to few refugial areas within the low mountain region, such as the Pfälzerwald, Eifel, and the Harz Mountains (Eckert et al. 2010). After hunting pressure was reduced, distance between refugial areas, habitat fragmentation, and increasing anthropogenic development left local wildcat populations geographically isolated (Hille et al. 2000;Pierpaoli et al. 2003). Considering the possible negative effects of isolated populations, including inbreeding and genetic drift, efforts to connect these distinct wildcat populations have become a conservation priority (Klar et al. 2012;Mattucci et al. 2016). Given the concern for wildcat population connectivity, genetic monitoring relying on hair trapping was established across known wildcat territories (Klar et al. 2008;Simon and Hupe 2008;Steyer et al. 2013). The results of increased monitoring showed strong evidence for the recent recovery of wildcat populations in various areas Steyer et al. 2016). Most recently, the species had been documented in areas distant from the known source regions, such as in the Bayerischer Wald in the Southeast and the Lüneburger Heide in Northern Germany (unpublished data). These new detections of the species in areas outside long-standing refugia and reintroduction sites suggest that the wildcat is capable of rapid dispersal and establishment across anthropogenically modified landscapes (scenario 1). However, it remains unknown if the species may have persisted in small, overlooked populations that are now increasing locally (scenario 2). Revealing the origin of wildcat reemergence is of considerable importance, especially considering the first scenario, expansion from refugia and reintroduction areas, which would imply that effective dispersal and population establishment of mobile mammals such as wildcats may be less obstructed by habitat fragmentation and barriers within anthropogenic landscapes than previously assumed (Klar et al. 2008;Jerosch et al. 2018;Balkenhol and Waits 2009).
Here, we aim to disentangle the above-mentioned recolonization scenarios by investigating the fine-scale genetic population structure of wildcats in a region where the species was recently detected, namely the Rhön Biosphere Reserve (RBR). The RBR is located on the southeastern edge of published wildcat distribution in the central low mountain region ( Fig. 1)  ). The RBR is a biodiversity-rich landscape, dominated by meadows and open land, surrounded by multiple large highways (A7, A71, A4) (Jedicke 2013). The RBR is home to a variety of protected species, including the wildcat. Wildcats in the RBR were thought to be extinct during most of the twentieth century until first genetic detections of wildcats occurred between 2007 and 2009 (Birlenbach et al. 2009). Beginning in 2009, under the frame of various monitoring projects (Table S1), the RBR and surroundings were heavily searched for evidence of wildcat, resulting in the detection of a substantial presence of wildcat in this region (Reiners et al. 2014;Thein 2008). Scenario 1 would imply wildcats expanded from adjacent populations into the RBR. Given the regional history, there are two possible source populations; (i) the refugial population in the north and (ii) the reintroduced population in the adjacent Spessart Mountains (Fig.  1). The known wildcat distribution to the north of the RBR was documented since before 2009 as a known source of a stable wildcat population (Birlenbach et al. 2009). The wildcat reintroduction in the Spessart region was launched in 1984 under the frame of a long-term reintroduction project. Approximately 600 wildcats were released until 2008 (Büttner and Worel 1990). Additionally, in 2005, six wildcats were reintroduced in the Neuwirthauser Forest, which is at the southernmost part of the RBR. It remains unknown if this reintroduction was successful and contributed to the reemergence in the RBR or if the current population resembles dispersal of individuals from the north. Scenario 2 would indicate a local increase from a small, overlooked relict population within the Rhön Mountains, which would likely result in population sub-structuring, due to the loss of rare alleles in small isolated populations, indicating the presence that a relict population survived and is now expanding in the low mountain region (Excoffier et al. 2009).
The present study aims to disentangle the two proposed scenarios to shed light on the origin of wildcats in this low mountain region. For this, data was obtained from a largescale genetic monitoring program conducted over the past 10 years in Germany. From this monitoring program, we genotyped a subset of samples with SNP and microsatellite markers to reveal the genetic makeup in the RBR. We compared the genetic structure of wildcats in the RBR with both adjacent source and reintroduced populations (Steyer et al. 2016) to test the two scenarios explaining the emergence of the species in the RBR given above.

Study site
The study area comprised the Rhön Mountains, including the RBR and surroundings (Fig. 1). The reserve spans over 2433 km 2 , which was expanded in 2014 to include the southern Rhön with the NF. The Rhön is a low mountain range ranging from 250 to 950 m, with approximately 40% forested land.
Within the study area, the closest potential source populations were included; mainly the region surrounding the National Park Hainich in the north, and the Spessart Mountains, bordering the Rhön Mountains to the southwest. The Hainich Mountains range from 225 to 494 m and are the most expansive broad-leaved forest, the primary habitat of the wildcat, in Germany. The Spessart Mountains include one of Germany's most forested areas, peaking at 586 m. In total, our study consisted of approximately 15,000 km 2 of variable habitats.

Sample collection and laboratory methodologies
Multiple opportunistic and standardized monitoring projects collected samples from 2004 to 2017. Opportunistic wildcat roadkill samples were collected throughout the study period.
In addition, multiple monitoring projects began in 2008, most notably Rhoen Natur e.V., BUND Wildkatzensprung, and Biosphärenreservat Rhön projects, which resulted in intense monitoring within the RBR during this time (for a complete list, see supplementary Table S1). In this study, we combined microsatellite and SNP genotyping methods. We took 119 wildcat microsatellite genotypes from 2004 to 2013 from a German-wide study on wildcat population structure (Steyer et al. 2016) and added 49 additional genotypes from 2014 to 2017 to be further analyzed using a wildcat specific SNPtype™ marker panel (von Thaden et al. 2020). Microsatellite genotypes were considered for SNP genotyping if they had a minimum of 11 loci and showed amplification success rates of 80%. Additionally, this selection excluded potential hybrids as there are low rates of hybridization across  (Steyer et al. 2018). The samples were then selected to create an even distribution in location, time, and sample type. Of the selected genotypes, 64 samples originated from tissue of dead found wildcats, and the remaining 104 samples are hair samples from various monitoring projects.
We used the primers LF4 (Eckert et al. 2010) and H16498 (Kocher et al. 1989) to sequence a 110 bp fragment of the mitochondrial control region following the protocol in Steyer et al. (2016). All samples were analyzed using 14 microsatellite markers and a zink finger sex marker (Hartmann et al. 2013). For non-invasively collected samples, the multiple tube approach using three replicates per sample was applied (Hartmann et al. 2013). Individualization was carried out through a custom R script and duplicated individuals were removed before further analysis (see Steyer et al. 2016 for further details). The microsatellite genotypes utilized from the Steyer et al. (2016) study were compared with the samples collected between 2014 and 2017 to look for batch effects. As none was found, the sample sets were combined.
The SNP genotyping of wildcat samples (n = 168) was performed on a EP1 platform (Fluidigm Corp., USA) using microfluidic 96.96 Dynamic Arrays™. Detailed methods and wildcat-specific SNPtype™ Assays are presented in von Thaden et al. (2017) andvon Thaden et al. (2020). The SNPtype™ marker panel encompasses 84 loci selected for individual and population identification, 10 loci for hybrid detection, and 2 SRY-linked loci. All SNP experiments included four no template controls (NTC) per array and noninvasively collected hair samples were triplicated to detect potential errors. From 168 samples, seven samples did not show sufficient (> 70%) amplification success to be scored and used in downstream analysis. Subsequent individualization resulted in 146 individuals to be analyzed using the combined SNP and microsatellite data.

Data analyses
First, we tested the combined SNP and microsatellite data set for isolation by distance (IBD) effects. We performed an IBD analysis as implemented in GeneAlEx 6.5 (Peakall and Smouse 2012) to account for genetic differentiation solely based on geographical position of individuals. The program STRUCTURE v2.3.4 (Pritchard et al. 2000) was used to evaluate population genetic structure. After 100,000 steps of burnin, 200,000 MCMC steps were performed with admixture model and correlated allele frequencies using a range of K 1-10 with 10 iterations. We used Structure Harvester (Earl and vonHoldt 2012) with the Evanno method (Evanno et al. 2005) to determine the most likely number of population clusters. The replicates were consolidated with the software CLUMPP (Jakobsson and Rosenberg 2007) using the GREEDY algorithm.
Further review of the spatial structure was carried out with the ADEGENET package in R (v.3.4.2) using a spatial principal component analysis (sPCA) (Jombart 2008). No requirements of the data to meet Hardy-Weinberg expectations or linkage equilibrium are needed for this method. In addition to the genetic data, sPCA also uses spatial information and is particularly suitable for the analysis of weak genetic structures (Storfer et al. 2007). sPCA relies on Moran's I (Moran 1948(Moran , 1950 to identify spatial patterns within the genetic structure of the sampled individuals. The method distinguishes between global scores, which indicate gradients in allele frequencies, and local scores, indicating differences in neighboring samples. Standard genetic diversity indices including expected and observed heterozygosity, allelic richness, and population pairwise F ST values were calculated for the combined SNP and microsatellite genotypes using Arlequin version 3.5.2.2 (Excoffier and Lischer 2010).

Results
Of the 146 genotyped individuals, mtDNA haplotypes were successfully determined in 134 individuals and corresponded to SNG-HP-FS03/-04/-06/-22/-23 (Steyer et al. 2016). Haplotype SNG-HP-FS23 was found in eight individuals solely in the southern part of the study region (Fig. 2b). The other four haplotypes appeared in all parts of the study area.
Analysis of genetic structure based on combined SNP and microsatellite genotypes revealed no evidence of isolation by distance (IBD; Mantel test: r = 0.011, p = 0.48). Bayesian clustering implemented in STRUCTURE indicated K = 2 as the most likely number of clusters within the study area. Separation of these clusters could not be attributed to a geographical pattern, with individuals clearly showing representation from both clusters (Fig. 3).
Clustering with sPCA resulted in a significant global structure, indicating correlation between the genetic and geographic distances (p = 0.01). However, no significant local structure was found (p = 0.14). A plot of lagged scores from the first principal component suggested the global structure is linked to a north-south genetic border, with a transitional area within the RBR (Fig. 2a). Subsequent principal components show weaker genetic structure although the eigenvalues suggest the first principal component explains most of the genetic variance found in wildcat individuals from the study area (Fig. S1).
Standard measures of genetic diversity were carried out based on the two clusters defined by the sPCA (Table 1). When looking at two clusters, observed (0.41, 0.43) and expected (0.42, 0.43) heterozygosity values from the north and south, respectively, were highly similar. Allelic richness was also identical (2.45, 2.45) along with a comparable fixation index (0.047, 0.011). The RBR, which comprised much of the south cluster, showed no local genetic signature such as higher allelic richness or specific clustering in the region. We also carried out pairwise F ST analysis based on the two clusters, which were not significantly differentiated from each other ( Table 2). We also compared our dataset with a larger, microsatellite-only sample set of 439 individuals originating from a previous study (Steyer et al. 2016), which revealed some sub-structuring in the south, specifically between the reintroduced region and the RBR (Fig. S2).

Discussion
In recent decades, wildcat populations across low mountain regions in Germany showed signs of expansion, despite a considerably fragmented landscape (Hartmann et al. 2013;Würstlin et al. 2016). While the recolonization process has been well documented (Canters et al. 2005;Nussberger et al. 2018;Streif et al. 2017), few attempts have been made to distinguish between active range expansion and locally growing populations. We assessed regional population structure of wildcat based on available samples from various monitoring projects in the region (Table S1). We investigated population sub-structuring and fine-scale genetic diversity to shed light on two possible scenarios for the rapid appearance of wildcats within this low mountain region.
The number of wildcat individuals found within the study region between 2009 and 2014 ( Fig. 1) indicates the presence of a viable wildcat population within the RBR. In part, this high number of wildcat detections can be attributed to the intense monitoring activities by multiple concurrent projects (Table S1). However, it appears unlikely that the marked increase of wildcat evidence since 2009 can be explained solely by increased monitoring activities. To our knowledge, local experts have continuously looked for wildcat presence in the region, making it unlikely that a population remained undetected for decades in the Rhön (Franz Müller, pers. comm.). Roadkill monitoring has occurred since 2004 rather opportunistically, involving local authorities, conservationists, and hunters. In addition, the overall observed trend of recent wildcat expansion across various regions within Germany (Steyer Fig. 2 Map showing the spatial genetic structure as assessed by sPCA analysis and haplotype map. a Individual scores from the first principal component: large blue squares indicated a highly positive score, large red squares indicated highly negative scores. b Map of individuals within the study area with the haplotype 23 (red) and all other haplotypes (3, 4, 6, 22) The past reintroduction of wildcats in the Spessart allows for tracking possible wildcat expansion in this region through a unique mtDNA haplotype (Steyer et al. 2016). While haplotypes SNG-HP-FS03/-04/-06/-22 are the most common haplotypes in German wildcat populations (Steyer et al. 2016), SNG-HP-FS23 is restricted to the Spessart Mountains (Fig. 2b). While this gives clear indication that reintroduced wildcats have successfully established in the Spessart region, the distribution of this haplotype also documents a lack of substantial spread into the RBR. Our findings of the haplotype SNG-HP-FS23 being largely restricted to the Spessart reintroduction area suggests that the expansion into the Rhön area is mainly driven by male dispersal from the Southern reintroduction area and confirms the general observation of maledominated dispersal in carnivores (Støen et al. 2005).
Our results from population structure analysis indicate a highly admixed population that comprises the northern refugial and the reintroduced populations. The sPCA results indicate that wildcats in the northern area of the RBR are genetically indistinguishable from northern refugia, which derive from the central German wildcat population described in Steyer et al. (2016) (Fig. 2a). In addition, the data suggests that wildcats from the reintroduction are moving into the RBR. Interestingly, STRUCTURE results show a highly admixed population with no clear geographic boundaries (Fig. 3). This may be explained by overall weak substructuring (F ST < 0.05) (Hubisz et al. 2009;Stift et al. 2019). The F ST measures also showed no significant deviation between the northern and southern clusters defined by the sPCA. This hints at a highly admixed population throughout the study region.
We found no evidence of an overlooked population within the RBR, as no private alleles were discovered in the area and we did not find any genetic sub-structure separating the Rhön from adjacent regions. Thus, the most probable explanation for the observed genetic pattern is a recolonization from both (northern and southern) directions. In case of overall low substructure within the central German wildcat population, which has been found previously (Steyer et al. 2016) and is confirmed in this study, a scenario of local population size increase together with significant dispersal from adjacent areas appears feasible. Still, the lack of evidence for wildcat appearance prior to 2008 makes the first option more likely and both scenarios ultimately require substantial dispersal from adjacent regions, suggesting a significant permeability of the landscape for wildcats.
Our results imply that wildcats can disperse through human-dominated landscapes. Land use within the newly recolonized RBR comprises 41% forested land (Jedicke 2013), with the remaining being settlement, open meadows, or arable land. Our study suggests that wildcats can, within a few years or few generations, recolonize a region with a substantial proportion of open land. This reemergence implies the wildcat is not significantly isolated due to road infrastructure, and can establish home ranges in primarily agricultural habitats, which has been supported by recent studies showing similar results (Jerosch et al. 2017;Jerosch et al. 2018;Klar et al. 2009;Würstlin et al. 2016). The presence of green bridges, built in 2011, potentially facilitates the exchange between individuals on either side of the A7 highway, as has been shown in other regions (Pir et al. 2011). The main factor determining habitat choice in wildcats is thought to be distance to forest (Sarmento et al. 2006). Jerosch et al. (2018) highlight structural heterogeneity in open landscapes as a determining factor allowing wildcats to persist in human-dominated landscapes. Therefore, the presence of rich-structured mosaics of open land and forest patches might have provided suitable conditions for successful recolonization in the RBR.
Our findings confirm a recolonization process of the wildcat in the Rhön Mountains, which has important implications for current wildcat conservation strategies. The results mainly suggest that connecting suitable habitat through stepping stones of rich-structured patches, rather than continuous forest, may be sufficient for wildcat dispersal. Therefore, it is important to conserve landscapes of rich-structured mosaics of open land and forest patches, similar to the Rhön, as habitat fragmentation continues to occur. As a flagship species, efforts to create heterogeneous landscapes will help other forestdwelling species to migrate between patches.  Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.