Past, present, future: tracking and simulating genetic differentiation over time in a closed metapopulation system

Genetic differentiation plays an essential role in the assessment of metapopulation systems of conservation concern. Migration rates affect the degree of genetic differentiation between subpopulations, with increasing genetic differentiation leading to increasing extinction risk. Analyses of genetic differentiation repeated over time together with projections into the future are therefore important to inform conservation. We investigated genetic differentiation in a closed metapopulation system of an obligate forest grouse, the Western capercaillie Tetrao urogallus, by comparing microsatellite population structure between a historic and a recent time period. We found an increase in genetic differentiation over a period of approximately 15 years. Making use of forward simulations accounting for population dynamics and genetics from both time periods, we explored future genetic differentiation by implementing scenarios of differing migration rates. Using migration rates derived from the recent dataset, simulations predicted further increase of genetic differentiation by 2050. We then examined effects of two realistic yet hypothetical migration scenarios on genetic differentiation. While isolation of a subpopulation led to overall increased genetic differentiation, the re-establishment of connectivity between two subpopulations maintained genetic differentiation at recent levels. Our results emphasize the importance of maintaining connectivity between subpopulations in order to prevent further genetic differentiation and loss of genetic variation. The simulation set-up we developed is highly adaptable and will aid researchers and conservationists alike in anticipating consequences of conservation strategies for metapopulation systems.


Introduction
Many wildlife species of high conservation concern face threats from habitat loss and habitat degradation (Lowe et al. 2005;Brook et al. 2008), subsequently resulting in population fragmentation. As populations become increasingly fragmented, they may enter a metapopulation system of geographically distinct subpopulations. Without sufficient migration, effects of genetic drift will result in increasing genetic differentiation between subpopulations, leading to reduced genetic diversity and increased extinction risk of small subpopulations (Frankham et al. 2010). Establishing or enhancing sufficient gene flow within such metapopulation systems, however, can delay or even reverse increasing differentiation (Lowe and Allendorf 2010), and therefore presents itself as a primary conservation strategy (Holderegger et al. 2019).
Decisions on the conservation of threatened metapopulations are ideally informed by predictions based on empirical estimates of migration rates and genetic differentiation between subpopulations (Bennet 2003;Kettunen et al. 2007). Conservation genetic analyses are essential to gain the basic insights needed for such decisions. Repeatedly performing genetic analyses over several time periods recognizes trends, which exceed one-time analyses (Schwartz et al. 2007). In addition to comparisons of past and present Florian Kunz and Annette Kohnen contributed equally to this work. genetic data, which facilitate valuable conclusions regarding observed genetic differentiation and its drivers, future projections might be particularly informative to support conservation decisions.
Parameterized with recent field data, such future projections allow for different scenarios to be designed and simulated forward in time. However, designing meaningful scenarios and parameterizing simulations is far from trivial. Populations of conservation concern are often reduced in size, display sex and stage biases, or suffer from reduced gene flow to neighbouring populations (i.e. increasing isolation) (Frankham et al. 2010). Therefore, most assumptions commonly used for simulations, like panmictic mating within infinitely large populations, usually do not hold (Hoban 2014). It is thus crucial to design simulations based on realistic life history parameters and demographics, while at the same time balance model complexity with the risk of underfitting (Hoban 2014). A simulation engine should allow for the specification and variation of these parameters when simulations are seeded with field data of species of conservation concern (Hoban et al. 2012;Hoban 2014). Simulations based on realistic scenarios ultimately enable researchers and conservationists to pinpoint priority areas for conservation actions and evaluate their potential effectiveness. Therefore, simulations of genetic data together with population models represent excellent tools to investigate genetic differentiation between populations and subpopulations, particularly within closed metapopulation systems.
A well suited model system to study and simulate effects of different migration scenarios on genetic differentiation is the Western capercaillie (Tetrao urogallus) in the Black Forest (Germany). Due to its specific range and habitat requirements, the capercaillie is considered as an umbrella species for forests rich in biodiversity (Suter et al. 2002;Pakkala et al. 2003), serving as a flagship species for conservation management (Mollet et al. 2008;Suchant and Braunisch 2008). Worldwide, capercaillie occurs over a large range (Coppes et al. 2015). However, in Western-and Central Europe populations are mainly restricted to mountain ranges and many of them are declining or became extinct in the past (Klaus et al. 1989;Storch 2001;Coppes et al. 2015). Here, we focus on the capercaillie population in the Black Forest, which underwent a severe decline in the past decade (Coppes et al. 2019). This population is well suited for studying genetic differentiation as: (1) it represents a metapopulation system consisting of four subpopulations; (2) the metapopulation system can be considered as a closed system, as distances to the next adjacent populations exceed the species' dispersal capacities; (3) it has been under investigation for over two decades, enabling the assessment of genetic differentiation over time.
We tracked genetic differentiation of capercaillie over time within the Black Forest metapopulation system, using two datasets from different time periods (a historic dataset, sampled from 1999 to 2004; and a recent dataset, sampled from 2013 to 2017). We then ran individual-based forward simulations informed by these field data, implementing different scenarios. Thereby, we aimed to predict genetic differentiation as a function of migration patterns and to model effects of isolation or of the re-establishment of connectivity on future genetic differentiation. We therefore included two from today's perspective realistic yet hypothetical scenarios, one simulating the extinction of one subpopulation, the other one simulating re-establishment of connectivity within the metapopulation. This allowed us to derive conservation actions and to highlight the strength of the simulation engine as a tool to analyse metapopulation systems.

Study area and focal species
This study was performed in the Black Forest, a lower mountain range (up to 1500 m a.s.l.) in south-western Germany (Fig. 1). The area harbors a population of capercaillie, a large forest dwelling grouse, inhabiting open to semi-open conifer-dominated forests at elevations above 800 m a.s.l. (Klaus et al. 1989;Storch 2001;Graf et al. 2009;Zohmann et al. 2014). Capercaillie numbers and range in the Black Forest have been declining over the past decades. While in 1900 the number of lekking males was estimated to 3800 individuals (Coppes et al. 2019), it declined to around 1300 individuals by mid-century (Roth et al 1990) and only 570 individuals in 1971 based on the first census across the Black Forest (Roth 1974). Since 1983, yearly censuses indicated a further decrease to a historic minimum of 167 lekking males in 2018 (Coppes et al. 2019). Since 1989, the range of capercaillie in the Black Forest decreased from 607 km 2 in 1993 to 344 km 2 in 2018 (Coppes et al. 2019). The main cause for this decline is assumed to be habitat deterioration (Kämmerle et al. 2020), but also increasing predation pressure (Kämmerle et al. 2017), increasing human disturbance ) and climate change (Braunisch et al. 2013) are considered as driving factors. Capercaillie mainly disappeared from small and isolated patches of its range (Kämmerle et al. 2017), while remaining occurrences are increasingly fragmented from each other (Coppes et al. 2019). Currently, the capercaillie population of the Black Forest is split into four, geographically separated subpopulations (i.e. North, Central, East, South, see Fig. 1), which are delineated by topography and landscape characteristics and by median dispersal distances (5-10 km, Storch and Segelbacher 2000). A previous study, including data from 1999 to 2004, indicated effects of barriers for gene flow between the northern and the southern part of the Black Forest area while overall genetic differentiation was found to be weak (Segelbacher et al. 2008). As dispersal between the Black Forest and its nearest neighbouring populations (i.e., Vosges mountains in France, Jura mountains in Switzerland) is highly unlikely due to large distances and unsuitable landscape features (i.e. intensive agricultural land, settlements), the four subpopulations within the Black Forest should be considered as a closed island population (Segelbacher et al. 2003).

Sampling and genotyping of historic and recent data
Datasets from two time periods were included in our analysis: faecal and feather samples collected from 1999 to 2004 (from here on referred to as historic dataset) and from 2013 to 2017 (from here on referred to as recent dataset). For the historic dataset, sampling and laboratory procedures are described in Segelbacher et al. (2008). Molted feathers of 213 individuals were collected non-invasively from 1999 to 1 3 2004 and were genotyped at ten microsatellites. In order to correspond to the recent dataset, we excluded one individual that showed missing values at more than three loci, resulting in 212 individuals used in the simulations runs (Table 1).
The recent dataset comprised 1278 faecal samples, which were collected non-invasively from December to April in the years 2013 to 2017. Sampling was conducted up to five days after the last snowfall to ensure high quality of DNA. Samples were genotyped at eleven pre-selected microsatellites (for details see the online appendix). Following the multipletubes approach (Navidi et al. 1992;Taberlet et al. 1996), a consensus genotype was accepted when at least two out of three replicates resulted in the same alleles for heterozygote loci and three out of three replicates for homozygote loci. In case of ambiguous results, additional three replicates were done and samples that still remained ambiguous after these six replicates were discarded afterwards. To check for contamination, negative controls were included in each extraction and PCR batch. Extractions and amplification were performed in separate rooms with regularly sterilized equipment. Samples that failed to amplify at more than three loci were excluded from further analyses.

Genetic diversity, population structure and recent migration rates
To check whether the used microsatellites were appropriate for distinguishing individual genotypes, the probability of identity (P ID ) and the probability of identity of siblings (P IDsib ) (Waits et al. 2001) were calculated using GenAlEx 6.503 Smouse 2006, 2012). To compare genetic diversity between datasets (sharing four loci), summary statistics were calculated per locus using GenAlEx 6.503 and allelic richness was calculated per locus and subpopulation using PopGenReport 3.0.4 (Adamack and Gruber 2014).
Both datasets were checked for deviations from Hardy-Weinberg equilibrium and for linkage disequilibrium using Genepop 4.2 on the web (Raymond and Rousset 1995;Rousset 2008) with default Markov Chain parameters (1000 step dememorisation, 100 batches, 1000 iterations). We corrected for multiple testing (sensu Narum 2006) by using the false discovery rate approaches by Benjamini and Hochberg (1995) and Benjamini and Yekutieli (2001), calculated in R 3.6.0 (R Core Team 2019). Micro-Checker 2.2.3 (Van Oosterhout et al. 2004) was used to test for large allele dropout, presence of null alleles, and stuttering.
Summary statistics of genetic diversity were calculated using GenAlEx 6.503. Allelic richness and private allelic richness were calculated using HPrare 1.1 (Kalinowski 2005).
Genetic population structure of each dataset was inferred using multiple approaches. We performed a Bayesian clustering approach on each dataset, implemented in Structure 2.3.4 (Pritchard et al. 2000), to estimate the number of clusters K and assign individuals accordingly, based on their genotypes. Each analysis was run with 200,000 burn-in and 500,000 MCMC iterations, with 25 iterations for each K from 1 to 7. Information about sampling location (population origin) was implemented using the locprior option (Hubisz et al. 2009). As sample sizes were uneven between subpopulations, an alternative ancestry prior (separate alpha for each subpopulation, and initial alpha of 0.25) was used (Wang 2017). The most likely number of clusters was estimated considering both delta K (Evanno et al. 2005) and mean log likelihood LnP(K) (Wang 2017), calculated by Structure Harvester Web 0.6.94 (Earl and VonHoldt 2012). Post-processing and visualization of runs were done using Clumpak (Kopelman et al. 2015). Population structure was further analysed using the multivariate discriminant analysis of principal components (DAPC, Jombart et al. 2010) implemented in R package adegenet 2.1.1 (Jombart 2008;Jombart and Ahmed 2011) in R 3.6.0. This multivariate method first transforms the data using a principal component analysis, followed by a discriminant analysis on the retained components. The resulting discriminant functions visualize the highest between-group variation, while at the same time minimizing within-group variation (Jombart et al. 2010). DAPCs were run including a priori information (population origin) and based on k-means clustering implemented in adegenet. To address pairwise differentiation between populations, we calculated mean pairwise G'' ST (Meirmans and Hedrick 2011) as fixation index and mean pairwise Jost's D est (Jost 2008) as differentiation index with the R package mmod 1.3.3 (Winter 2012). We calculated 95% confidence intervals based on 1000 bootstrapped iterations, and corrected for bias following the method implemented in diveRsity 1.9.9 (Keenan et al. 2013). An AMOVA was calculated using Arlequin 3.5.2.2 with 10,000 permutations. We used BayesAss 3.0.4 (Wilson and Rannala 2003) to estimate recent migration rates and 95% credible intervals. We conducted ten independent repeats of 50 × 10 6 iterations (including 5 × 10 6 iterations burn-in) with a sampling frequency of 2000, each initiated with a different random seed for each dataset. In order to keep the acceptance rates for proposed changes between 40 and 60%, delta values were adjusted to Δm = 0.1, Δa = 0.17 and Δf = 0.17. Convergence of chains was confirmed using Tracer 1.7.1 (Rambaut et al. 2018) and by checking for concordance between repeats. We used the Bayesian deviance (Meirmans 2014) to search for the best fitting model (selecting the one with the lowest Bayesian deviance) (Faubet et al. 2007).

Forward simulation of migration scenarios
We performed forward simulations, using individual-based models implemented in the R package rmetasim 3.1.7 (Strand 2002). rmetasim provides a flexible environment that incorporates demographic parameters, migration rates and genetic parameters in the modelling procedure. It offers the possibility to seed simulations using existing genotype data. Stepwise forward simulations can then be run, tracking the fate and genotype of individuals per step, which in our case refers to one year (Strand 2002). After running the simulations, we extracted the genotypes of all individuals in each population and calculated observed and expected heterozygosity using adegenet and pairwise G'' ST and Jost's D est using mmod 1.3.3.
We built a stage-based transition model for capercaillie including three stages: juveniles, adult females and adult males. Although some studies indicate slight differences in juvenile survival rates between sexes (Wegge 1980;Klaus et al. 1989;Hörnfeldt et al. 2001), we lumped sexes into a single stage to reduce model complexity. Demographic parameters of capercaillie are multifactorial and differ between years and habitats (Klaus et al. 1989;Grimm and Storch 2000;Kangas and Kurki 2000;Åhlen et al. 2013;Jahren et al. 2016;Augustine et al. 2020). We therefore based the parameterisation on estimates reported from a nearby comparable population (in the Bavarian Alps, Grimm and Storch 2000). Survival rates were set to 0.73 and 0.85 for adult females and males respectively, and 0.36 for juveniles. Reproduction rate was set to 1.5, based on mean clutch size (7), clutch survival (0.65), hatching success (0.95) and chick survival rate (0.342, averaged over both sexes, Grimm and Storch 2000).
Density dependence is assumed to negatively affect survival and reproduction in capercaillie (Kangas and Kurki 2000;Sachot et al. 2006). We thus considered effects of density dependence by reducing female and juvenile survival rates and reproduction rate at carrying capacity. Carrying capacities of subpopulations were estimated by dividing patch sizes by the average size of a female's home ranges (following Sachot et al. 2006), which resulted in 100 to 400 individuals for each subpopulations (North: 400, Central: 100, East: 100 and South; 200). However, the resulting estimates should be seen as rough proxies, as the extent of available habitats can vary as a consequence of changing forest management practices. The simulations were initially seeded using the genotypes, sexes, and subpopulation assignment of the historic and recent datasets. We used all loci represented in the datasets and assumed a general mutation rate of 0.00045 (Whittaker et al. 2003), following the stepwise mutation model.
We ran two simulations, each implementing several scenarios (Table 2). We initialized the first simulation (historic to recent, abbreviated as HR) with the historic dataset and then ran the simulation for 15 years. We implemented three scenarios and compared their results with the genetic differentiation from the recent dataset: The first scenario did not comprise migration between the subpopulations (HR_1: no migration). The second scenario contained migration rates that were derived from empirical analyses of recent migration using BayesAss (HR_2: evidence-based estimates). In the third scenario, we predefined identical migration rates between all subpopulations corresponding to the highest empirically observed migration rate (HR_3: ideal migration). This allowed us to explore a theoretical minimum and maximum differentiation between subpopulations and to evaluate the simulations' set-up (i.e. by comparing scenario HR_2 to the recent dataset). We then initialized the second simulation (recent to future, abbreviated as RF) with the recent dataset, and ran the simulation for 35 years into the future. Thereby, we implemented the same three scenarios as for the HR simulation and two further scenarios derived from recent conservation considerations. This allowed us to contrast potential effects of different migration scenarios on the population differentiation in approximately 2050.
For scenario RF_1, we had to increase the female survival rate slightly to 0.77, so that the simulation would not result in subpopulation numbers too low to be analysed-especially for the recently small subpopulation Central, which was initialized with only 31 individuals (Table 1). We implemented scenario RF_4 to explore effects of isolation/ extinction of the subpopulation East. This subpopulation has experienced the largest decline in the past decades and showed the lowest number of lekking males in 2018 (Coppes et al. 2019). In contrast, we implemented scenario RF_5 to explore effects of an increased migration from and to subpopulation North. According to our analyses of historic and recent population structure, gene flow from and to North markedly decreased in recent times. Therefore, scenario RF_5 explores the potential conservation action of re-establishing gene flow from and to North.
As rmetasim is a stochastic simulation program, we replicated each scenario of both simulations 1000 times (confirmed by an a priori power analyses using the R package pwr 1.2.2, Champely 2018). We analysed each replicate by calculating observed and expected heterozygosity using adegenet 2.1.1, and the pairwise fixation index G'' ST as well as the pairwise differentiation index Jost's D est using mmod 1.3.3. We calculated the mean and the 95% confidence intervals (function: qnorm) over all 1000 replicates per scenario (Hoban et al. 2012).

Microsatellite data and observed genetic diversity
After correcting the historic dataset for multiple testing, no deviations from Hardy-Weinberg equilibrium and no evidence for linkage disequilibrium were found within subpopulations. Micro-Checker 2.2.3 indicated the presence of null alleles at locus BG5. Repeating further analyses with omission of this locus did not change the conclusion of results, thus we show results based on data including this locus.
Within the recent dataset, the genotyping success rate was 88%. Locus sTUT2 showed high amounts of missing values (> 30% in all four subpopulations) and was therefore discarded from further analyses. In addition, locus sTUT1 deviated from Hardy-Weinberg equilibrium, Micro-Checker 2.2.3 indicated the presence of null alleles. Therefore, sTUT1 was also discarded from further analyses, reducing the number of microsatellite loci to nine in the recent dataset. We found no evidence for linkage disequilibrium.
The remaining loci in both datasets were powerful enough to detect individuals (historic dataset: P ID = 5.6 × 10 −7 , P IDsib = 2.1 × 10 −3 ; recent dataset: P ID = 1.7 × 10 −8 , P IDsib = 7.1 × 10 −4 ). The loci TUT3, TUT4, BG15 and BG18 (shared in both datasets) displayed highly comparable variability, as shown by their allelic richness per subpopulation. In general, the recent dataset was more variable than the historic dataset due to a high allelic richness of locus sTUD5 compared to a lower allelic richness of the loci BG6 and TUT10 (Table 1 in the online appendix). Genetic diversity was generally higher in the recent dataset compared to the historic dataset (Table 1). However, this result might be affected by the overall higher genetic variability of the microsatellites used in the recent dataset, rather than biological reasons. There was no difference between subpopulations in both datasets. F IS values were not significant for all but one subpopulation (North in the recent dataset), indicating random mating within subpopulations.

Population structure and genetic differentiation from past to present
Results of Structure 2.3.4 without a priori information on sample origin did not show any differentiation of individuals into subpopulations for both datasets (not shown). Using the locprior models for the historic dataset, the most likely number of clusters K ranged from 1 to 3 (with LnP(K) indicating one cluster while ΔK indicated three clusters). Irrespective of the number of clusters, Structure could not distinguish between subpopulations, as all individuals shared comparable proportions of clusters ( Fig. 1 on the online appendix).
Within the recent dataset using locprior informed models, the most likely number of clusters K was not clear either. K ranged from 5 to 6 (with LnP(K) indicating five clusters, while ΔK indicated six clusters). Individual assignments however displayed a specific pattern of differentiation between the subpopulations (Fig. 2). While North appeared differentiated as a whole, some individuals of Central appeared to be similar to individuals of East. This indicated migration from East to Central. South displayed an admixed pattern, with individuals appearing similar to individuals of East as well.
The DAPCs showed an increase in differentiation between subpopulations from the historic dataset to the recent dataset, especially for the subpopulations North and Central (Fig. 2). Explained variance retained by the PCA principal components were 96% for the historic dataset and 90% for the recent dataset. The proportion of reassignment, which is the ability of the DAPC to reassign individuals into their original clusters, increased from 61 (historic dataset) to 75% (recent dataset), indicating an increased strength of the genetic signal of differentiation.

Fig. 2 DAPC and
Structure results for the historic and recent datasets. DAPCs are shown for the first vs. the second discriminant function (x-axis and y-axis, respectively), Structure results are shown for the most probable number of clusters K inferred by Evanno's ΔK (Evanno et al. 2005) and LnP(K) (Wang 2017) (K = 3 for the historic dataset and K = 5 for the recent dataset) with the models including location information. Further Structure plots for K = 3-7 for both datasets are shown in the online appendix, Fig. 4 The comparison of the mean pairwise fixation index G″ ST (Meirmans and Hedrick 2011) and differentiation index Jost's D est (Jost 2008) revealed an increase at all pairings from the historic dataset to the recent dataset (Fig. 3a). Both indices were highly positively correlated (hence results for Jost's D est are presented in the online appendix, Fig. 4a).
AMOVAs for the historic dataset and the recent dataset both yielded high amounts of variance within individuals, with only 1.69% of variance (p < 0.001) for the historic dataset and 2.88% (p < 0.001) for the recent dataset due to differences between subpopulations (Table 2 in the online appendix).
Log likelihood was comparable between BayesAss runs, as was Bayesian deviance. Proportions of nonmigrants per subpopulation did not approach either 66 or 100%, altogether indicating a good fit of the Bayes-Ass model (Meirmans 2014). Significant recent migration (Table 3) was highest between East and South (both directions, m BayesAss = 0.11/0.13) and from East to Central (m BayesAss = 0.13). While uni-directional significant migration rates were found from South and Central to North and from Central to South, they manifested themselves only in low rates (ranging from m BayesAss = 0.07 to 0.04).

Simulations of genetic differentiation
The scenarios for both simulations resulted in stable population numbers, and therefore proved to be useful for subsequent analyses. The three scenarios used for simulation HR resulted in different pairwise indices of genetic fixation and differentiation, although partially overlapping confidence intervals indicate variability (Fig. 3b). Scenario HR_1 (no migration) and scenario HR_3 (ideal migration) resulted in highest and lowest pairwise indices, respectively. Indices resulting from scenario HR_2 (evidence-based estimates) were comparable to the indices calculated from the recent dataset. Therefore, our simulation parameters (including the evidence-based estimates for migration rates) proved to be appropriate to generate plausible results of forward simulated population structure.
Scenario RF_1 (no migration) and scenario RF_3 (ideal migration) resulted in the lowest and highest indices of genetic fixation and differentiation (Fig. 3c). Compared to recent pairwise indices of fixation and differentiation (Fig. 3b), indices approximately doubled over the simulated 35 years. Scenario RF_2 (evidence-based estimates) predicted increasing genetic fixation and differentiation when current migration rates are maintained. Scenario RF_4 (isolation of East) predicted higher pairwise indices of differentiation compared to scenario RF_2. The simulation yielded population numbers for East ranging from less than 10 individuals to 0. This was in line with the scenario assumption of isolation and subsequent extinction of subpopulation East. However, pairwise indices of fixation and differentiation based on such low numbers are not informative and are thus not displayed in detail. Scenario RF_5 (reestablishment of connectivity to North) resulted in overall reduced population differentiation, with some pairwise indices of fixation and differentiation being comparable to recent levels.
Observed and expected heterozygosity per subpopulation indicated an increasing trend with increasing migration rates ( Table 3 in the online appendix). Within simulation HR, scenarios HR_1 and HR_2 resulted in slightly reduced heterozygosity whereas heterozygosity of scenario HR_3 remained unchanged, compared to the heterozygosity of the historic dataset. Within simulation RF, scenarios RF_1 and RF_2 again resulted in slightly reduced heterozygosity whereas heterozygosity of scenario RF_3 remained unchanged, compared to the heterozygosity of the recent dataset. Within scenario RF_4, heterozygosity slightly decreased in the remaining subpopulations, while within scenario RF_5, heterozygosity in the subpopulation North increased, compared to the heterozygosity of scenario RF_2. However, observed differences in heterozygosity were generally small.

Discussion
Comparing two datasets, sampled in two different time periods, we were able to trace changes in genetic differentiation between subpopulations within a closed metapopulation system. Combining this with individual-based forward simulations allowed us to explore future genetic differentiation as a consequence of different migration scenarios. Our simulations have shown that genetic differentiation within a metapopulation system is highly dependent on gene flow. With about 15 years between the two sampling periods (roughly 2000 to 2015), we found already increased genetic differentiation between the subpopulations. Simulations for another 35 years (to 2050) revealed a further increase in genetic differentiation if migration patterns stay the same.

Tracking genetic differentiation from past to present
Our results revealed an increase in population structure and genetic differentiation in Black Forest capercaillie within Fig. 3 Comparison of pairwise G'' ST and bias corrected confidence intervals over datasets and simulation scenarios. a Contrasts the genetic differentiation of the historic dataset, the recent dataset and the scenario RF_2 (recent to future), which is based on recent migration rates. b Shows the scenarios of simulation HR (historic to recent) and contrasts those to the recent dataset. c Shows the scenarios of simulation RF (recent to future) ◂ a short period of approximately 15 years (Fig. 2). For the historic dataset, we found no signs of population structure or increased differentiation between specific subpopulations, which is in line with previous studies (Segelbacher et al. 2008). While the slight differentiation of the subpopulation North corresponded to our expectations, high rates of historical migration could have been counterbalancing differentiation effects (Segelbacher et al. 2008), ultimately preserving a nearly panmictic structure (Lowe and Allendorf, 2010). In the recent dataset, however, subpopulation North appeared differentiated from the three other subpopulations. Subpopulation North is geographically separated from the other subpopulations by a large valley with low mountains and poor habitat suitability for capercaillie (Braunisch and Suchant 2007;Coppes et al. 2019). Additionally, increasing habitat deterioration (Kämmerle et al. 2020) and human disturbance ) might act as driving factors of the observed increase in subpopulation differentiation. The subpopulations South and Central appeared still wellconnected with high bi-directional migrations rates, while migration from East to Central was one-directional. These results indicate that the differentiation between North and the other subpopulations has increased due to limited migration between the subpopulations. Conservationists should thus treat these first signs of differentiation (Segelbacher et al. 2008) as early warning signs of declining functional connectivity and plan management strategies to increase migration rates or to reduce factors hindering migration of the target species.
Our HR simulation showed that the increased differentiation can be attributed at least partially to the ongoing segregation of subpopulations and subsequent loss of gene flow. Comparing scenarios of evidence-based migration rates with simulated optimal conditions (i.e. HR_2 vs. HR_3, Fig. 3b) highlights the effect of migration rates on genetic differentiation. The long-term negative population trend of capercaillie in the Black forest, paired with ongoing habitat contractions (Coppes et al. 2019), resulted in small and fragmented subpopulations. However, within a metapopulation system, migration between subpopulations is essential to compensate for small population sizes and increased extinction risks (Frankham et al. 2010). While single large populations can still harbour low levels of genetic diversity when being isolated (Rutkowski et al. 2017), well-connected metapopulation systems can preserve higher levels of genetic variability even when subpopulation sizes are low (Alstad 2001;Allendorf et al. 2012). Accordingly, we found a slight trend of increasing levels of heterozygosity with increasing migration rates. Yet, the scenario with optimal migration rates (HR_3) could only maintain heterozygosity on historic levels, while the scenario using evidence-based migration rates (HR_2) resulted in decreased heterozygosity compared to historic levels.
When assessing population structure and genetic differentiation, potential time lags have to be considered. The manifestation of signals in genetic patterns is highly dependent on a species' dispersal capability (Landguth et al. 2010). Given median dispersal distances of capercaillie of about 5 to 10 km (Storch and Segelbacher, 2000), the detected increase in differentiation might not be related to specific fragmentation events in the past 15 years, but rather represents cumulative effects within a much longer time period.

Simulating genetic differentiation from present to future
Forecasting genetic differentiation using evidence-based migration rates, we observed a further increase by approximately 2050 (Fig. 3a). However, it is likely that our forecasts based on these rates might even underestimate genetic differentiation. The migration rates were appropriate in simulating the historic dataset to recent times, yet the population decline has accelerated over the entire study period (Coppes et al. 2019), thereby reducing gene flow to a larger extent. Additionally, our simulations featured stable population numbers (cf. Grimm and Storch 2000), further raising serious conservation concerns. Keeping population numbers stable, which is per se a challenging conservation target, evidently does neither prevent further genetic differentiation sufficiently nor preserve the metapopulation's genetic variability.
By implementing different migration rates between the subpopulations, we explored two realistic scenarios. Scenario RF_4 addressed effects of a potential extinction of one subpopulation (East) on the genetic differentiation between the remaining subpopulations. Considering the low population size of the subpopulation East (Coppes et al. 2019), stochastic effects could have a high impact on the subpopulation, and potentially lead to its extinction. Although initiated with recent migration rates for the remaining pairs of subpopulations, the simulation resulted in levels of genetic differentiation that are comparable to the scenario with no migration (RF_1). The subpopulation East therefore seems to act as an important core area for gene flow, connecting at least the southern subpopulations. The genetic differentiation between the subpopulations Central and North increased as well in this scenario (compared to RF_2), although the subpopulation East is not directly connecting these two subpopulations. However, with the simulated extinction of the subpopulation East, subpopulation Central became relatively isolated as well, probably driving the observed increase in genetic differentiation. This scenario clearly indicates that changes in one subpopulation might have far-reaching effects on all other subpopulations within a metapopulation system. Addressing a potential conservation action, we simulated effects of an improved connectivity from and to the subpopulation North on the entire metapopulation system within scenario RF_5. Based on existing data of habitat suitability (Braunisch and Suchant 2007) and functional connectivity for capercaillie in the Black Forest (Braunisch et al. 2010) as well as on our own analyses, the re-establishment of corridors to the subpopulation North and particularly corridors between North and Central appeared to be important for preventing further isolation. Our simulations showed a general positive effect of this management action in terms of overall decreased genetic differentiation compared to RF_2. While all pairwise comparisons of subpopulations including the subpopulation North were particularly affected, also the genetic differentiation between the remaining subpopulations in the southern regions of the Black forest was reduced. Additionally, heterozygosity was higher especially for the subpopulations North and Central, compared to scenario RF_4. This again supports the idea, that differentiation between all pairs of subpopulations in an entire metapopulation system might be affected by changes in one single subpopulation.
Increasing migration rates between subpopulations may only lead to maintaining rather than reducing genetic differentiation in the future. Landguth et al. (2010) found a strong nonlinear relationship between losing a historic barrier's signal in the genetic pattern and the species' dispersal capabilities. In highly mobile species, the signal will be lost within several generations. By contrast, species with limited dispersal capabilities were found to accomplish the same within tens or hundreds of generations. This might also be true for capercaillie with comparatively low median dispersal distances of 5 to 10 km (Storch and Segelbacher 2000). Accordingly, dispersal of two male individuals from the subpopulation North to the subpopulation East and of one female individual from the subpopulation East to the subpopulation Central has been found, thereby representing rare cases of migration between subpopulations (with n = 3 out of 1278 samples).

Conclusions for conservation
In the light of the drastic decline in population size and habitat area of capercaillie in the Black Forest (Coppes et al. 2019), conservation strategies must be planned and implemented urgently. While we found genetic differentiation to be driven by migration rates between subpopulations, preservation and increase of population size are important prerequisites to enable migration. Yet our simulations clearly showed that genetic differentiation can increase severely even in the case of stable population numbers. We therefore emphasize the need for conservation strategies that aim at re-establishing functional connectivity between subpopulations. Related actions might include improving habitat suitability within corridors or creating stepping stones. Our simulations have shown that increased migration rates between subpopulations can counteract genetic differentiation even in metapopulation systems with realistic demographic parameters and low population numbers. Re-establishment of functional connectivity could help maintaining genetic differentiation at its current level. Further, we strongly recommend the continuation of periodic genetic monitoring which will allow for adjusting and improving the predictions of future genetic differentiation along with predicting consequences of conservation actions.
Individual-based simulations, as used in the present study, can be valuable tools for both scientists and practitioners (Epperson et al. 2010;Hoban 2014;Holderegger et al. 2019). The flexibility of simulations allows for a variety of scientific questions to be examined and for conservation actions to be designed and tested (Hoban et al. 2012). We showed how simulations can support analyses of genetic differentiation in metapopulations systems and how realistic scenarios can provide valuable insights for species conservation.
Frank Zachos and two anonymous reviewers for valuable comments on the manuscript.
Funding Open access funding provided by University of Natural Resources and Life Sciences Vienna (BOKU). This study is partly funded by the Ministry of the Environment, Climate Protection and the Energy Sector Baden-Württemberg and the Ministry for Rural Affairs and Consumer Protection Baden-Württemberg and co-funded by Elektrizitätswerk Mittelbaden AG & Co. KG, EnBW-Energie Baden-Württemberg AG, ENERCON, German Wind Energy Association, Ökostromgruppe Freiburg, the Swedish Environmental Protection Agency (Vindval) and Windkraft Schonach GmbH. The funding organisations had no influence on the manuscript, study design, methods or the interpretation of the results.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Informed consent All authors approved this version of the manuscript for publication.
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://creat iveco mmons .org/licen ses/by/4.0/.