Introduction

Red deer (Cervus elaphus) are among the most widespread ungulates in Europe and one of the most iconic game species. They have been heavily impacted by anthropogenic influences such as habitat fragmentation, translocations and selective hunting for centuries (Hartl et al. 2003). As such, red deer have been the target of many population and conservation genetic studies analyzing the genetic diversity and population structure in human-dominated landscapes (e.g., Kuehn et al. 2003; Pérez-Espona et al. 2008, 2009; Fickel et al. 2012 Frantz et al. 2017). The aims of these studies varied, and included the quantification of genetic diversity in isolated and sometimes inbred populations (e.g., Zachos et al. 2007), estimating the amount and genetic consequences of translocations (e.g., Pérez-Espona et al. 2009; Haanes et al. 2010), or characterizing the genetic impacts of postglacial recolonization (e.g., Krojerova-Prokesova et al. 2015).

In Schleswig–Holstein, Germany’s northernmost federal state, red deer are distributed across the north, southeast, and center of the state (Fig. 1). The local populations are managed in 12 administrative units. These units were not established on the basis of population structure, but rather were opportunistically located in areas with high red deer densities, mostly located around larger patches of forest (Meißner et al. 2008; Wotschikowsky 2010). Units located in close proximity to each other such as Barlohe (BAL), Iloo (ILO) and Schierenwald (SCW) are demarcated by spatial jurisdictions (e.g., municipalities or communities) or landmarks (e.g., rivers or roads), rather than by natural boundaries or population structure. Such administrative considerations are commonly included when delineating wildlife management units (Taylor and Dizon 1999). Yet, from a genetic standpoint, populations should only be treated as separate management units when their genetic divergence is high enough to suggest demographic independence, meaning that the rate of dispersal among populations must be low (Palsbøll et al. 2007). To emphasize that red deer management units in our study area (Schleswig–Holstein) are not based on population structure, we refer to them as administrative management units (AMUs) and distinguish them from units defined by genetic divergence (i.e., genetic management units; GMUs).

Fig. 1
figure 1

Map of Schleswig–Holstein (study area). Inset indicates location within Germany. The blue line in the center indicates the Kiel Canal. Broad dashed black lines represent major highways (Autobahn). Red deer management units are delineated with thin dashed black lines. Forested areas are indicated by dark green shading. Local deer management units of which samples were included are Northern Friesland (NFL), Elsdorf (ELD), Barlohe (BAL), Iloo (ILO), Schierenwald (SCW), Hasselbusch (HAB), Segeberger Heide (SEG), Duvenstedter Brook (DUV), Lauenburg West (LAW), East (LAE) and South (LAS) as well as Sachsenwald (SAW). The two reference areas Denmark (DK) and Mecklenburg-Western Pomerania (MWP) are delineated in red. Triangles represent larger cities throughout that area

Full size image

Historically, red deer within the AMUs in Schleswig–Holstein (SH) have been subject to various anthropogenic restrictions. For instance, until recently (i.e., 1980), red deer were only allowed to freely range in so called ‘designated red deer zones’ (Meißner et al. 2008; Wotschikowsky 2004, 2010). This policy was intended to prevent damages to crops and forests by red deer. Culling of all individuals outside these zones consequently limited gene flow between established populations (Ströhlein et al. 1993; Willems et al. 2016). Today, infrastructures such as fenced highways (Autobahn) or the Kiel Canal form potential barriers to gene flow across the entire state (Fig. 1). Additionally, estimated population sizes vary greatly among the AMUs (range 35–530; see Table 1) and many of them contain fewer than 100 individuals. The Hasselbusch population (HAB) was founded by dispersed individuals from the Segeberger Heide (SEG) in the second half of the nineteenth century (ca. 1870), but has been isolated from its source for decades. More recently, a fenced highway has prevented any potential migration between SEG and HAB (Meißner et al. 2008). A previous study found low genetic diversity as well as the first signs of inbreeding for the Hasselbusch AMU (Zachos et al. 2007). For example, multiple animals with brachygnathia inferior (shortened lower jaw), a condition linked to inbreeding depression, have been found in the HAB population (Zachos et al. 2007). Furthermore, there are influences of translocations: the Duvenstedt (DUV) population is not native but goes back to an enclosure population founded with red deer from Austria, Hungary and Poland which was released in the 1950s (Jessen 1988; Meißner et al. 2008). Within the last decade, red deer have dispersed from Denmark, established themselves south of the German border and are increasing in numbers (Reinecke et al. 2013). As a consequence, the latest red deer AMU established in Schleswig–Holstein was the Nordfriesland unit (NFL). In the neighboring state of Mecklenburg-Western Pomerania (MWP) located south-east of Schleswig–Holstein, red deer are more abundant and have been roaming the state with less restrictions while occupying a large area (Kinser et al. 2010). Therefore, an exchange of individuals from these populations could result in higher levels of genetic diversity in the three AMUs located in the Lauenburg area (LAW, LAE, LAS).

Table 1 Overview of sampled management units and reference areas including estimated population census size and number of samples (n) for each area
Full size table

Hunters and landowners participate in the management of red deer within the 12 AMUs in order to set different management goals such as hunting quotas (Wotschikowsky 2010). Therefore, managing the AMUs separately assumes that these units equate to GMUs, thus representing more or less disconnected (i.e., closed or genetically separated) populations that experience limited reproductive exchange of individuals with other populations (Moritz 1994). However, several recent studies have shown that if this implicit assumption is violated in wildlife management, actions in one management unit (MU) can substantially influence management effectiveness in neighboring units (Hemami et al. 2005; Robinson et al. 2008; Olea and Mateo-Tomás 2014; Stillfried et al. 2017). In such cases, management would need to be extended towards a larger spatial scale that includes multiple MUs and considers the degree of connectivity among them (e.g., Robinson et al. 2008; Wäber et al. 2013). Genetic approaches have been suggested for delineating more meaningful management units based on biological population entities (e.g., Moritz 1994; Palsbøll et al. 2007). Strong genetic sub-structuring or varying levels of genetic diversity among areas are still the metrics of choice commonly used to justify the separation of MUs (e.g., Wilting et al. 2015; Grosser et al. 2017; Gaillard et al. 2017). However, novel analytical tools now allow researchers to derive estimates of directed dispersal rates from genetic samples, which can provide important information on potential source-sink dynamics and gene flow (e.g., Draheim et al. 2016).

Overall, the history of red deer in SH and the different anthropogenic influences on the local populations raise the question of whether the current practice of managing each AMU as a separate, closed population is appropriate. In particular, it is questionable whether genetic diversity within AMUs is high enough in order to sustainably counteract genetic drift, thereby preventing a loss of genetic diversity and inbreeding. We expect some AMUs to be linked by dispersal and gene flow rates high enough to warrant management as a single unit. If this is the case, red deer AMUs in SH can be interpreted as a network of subpopulations where local populations are connected by gene flow of varying degrees (Pannell and Charlesworth 2000). If so, we should observe different levels of genetic exchange among AMUs and of genetic diversity within AMUs, with migration depending on connectivity among neighboring AMUs, and genetic diversity depending on a combination of connectivity and population size of neighboring AMUs.

To assess the genetic structure of red deer AMUs in Schleswig–Holstein, we make use of an extensive data set consisting of over 500 tissue samples collected over multiple years. We estimate measures of genetic diversity and test the hypothesis that diversity will vary between the AMUs in Schleswig–Holstein but still be relatively low compared to other populations throughout Europe (Zachos et al. 2016). For this, we also added samples from two reference areas located in the neighboring state of Denmark and the federal state Mecklenburg-Western Pomerania. By combining analyses of genetic differentiation and population structure with a new approach of genetically-derived estimates of relative migration rates (Sundqvist et al. 2016), we also delineate clusters of AMUs that are connected by gene flow and thus should be managed as one GMU. In order to further confirm the genetic structure of the AMUs, we correlate observed patterns of genetic diversity, differentiation and gene flow to available information on current population size and density at the local and regional scale.

Materials and methods

Study area

The study area extends over approximately 15,580 km2 and covers the entire mainland of the federal state of Schleswig–Holstein in Northern Germany, south of the border with Denmark (Fig. 1). The state comprises a mosaic of different types of land use, predominantly agriculture and pastures. Forested areas are scattered across the state but vary substantially in size and composition of tree species. Larger forest complexes such as the Segeberger Heide (SEG) form the core areas of the red deer distribution throughout the study area (Fig. 1). Often, red deer habitats are further characterized by mixtures of marshes, heathlands and moors. Administrative management units vary in size from 13,000 up to 48,000 ha (Reinecke et al. 2013). Distances between AMUs range from a few kilometers (< 5 km) up to 63 km between the NFL and ELD units. Available information suggests that local populations range in size from 30 to nearly 600 individuals within the AMUs (Table 1).

Schleswig–Holstein is not densely populated (182 people per km2; Statistisches Bundesamt 2018) compared to the German average (237 people per km2), with human settlements and villages scattered across the state. The landscape is fragmented by roads, major highways (Autobahn) and canals (e.g., the Kiel Canal), all of which form potential barriers to the movements of red deer (Pérez-Espona et al. 2008; Frantz et al. 2012).

Sampling

We obtained 279 genetic samples from red deer harvested during the hunting seasons of 2013 to 2015. In order to ensure a sufficient sample size across all 12 AMUs, we included 186 samples collected in previous studies (Zachos et al. 2007; Reinecke et al. 2013) during the years 2003 and 2004. Additionally, we used samples obtained from two reference areas for comparative purposes: (1) 34 samples from the Froslev forest located in Southern Denmark (DK) close to the German border, and (2) 46 samples from several forests within the federal state of Mecklenburg-Western Pomerania (MWP) neighboring Schleswig–Holstein in the Southeast (Fig. 1). This lead to a total sample size of 545 (149 female, 104 male, 292 with no sex ID) red deer individuals (overview on sampling periods and sample sizes provided in supplement S3). Since free ranging red deer can live up to over 12 years (e.g., Guinness et al. 1978) the gap between the two sampling periods corresponds to a maximum of only one deer generation.

All samples were re-genotyped for our marker set in order to be fully comparable. We only considered samples for which the spatially referenced location of harvest (e.g., the forest complex) was reported. Individuals from MWP originated from areas not directly neighboring our study area. Therefore, these samples were only included for comparative measures regarding genetic diversity whereas the DK samples were also used throughout the analyses on differentiation and gene flow.

DNA extraction and genotyping

DNA was extracted using the ‘all tissue DNA’ kit (Gen-Ial, Troisdorf, Germany) following the manufacturer’s instructions (final DNA-elution in 80 µl). DNA concentrations were measured spectrophotometrically using a NanoDrop1000 (PeqLab GmbH, Erlangen Germany).

To genotype each individual, we used a panel of 14 microsatellite loci (see supplement S1). One primer of each of the 14 pairs was 5′-labelled with a fluorescent dye (6-FAM or HEX). To save time and costs, primers were combined (after optimization) in multiplex mixes (CerMix1–CerMix4). CerMix1 contained primers for four loci (INRA6, C143, T40, and T115), CerMix2 combined three loci (C105, C180, and C229), CerMix3 combined four loci (T107, Haut14, ILSTS06, and BM757), and CerMix4 included three loci (CSSM14, FSBH, and BM1818). The genotyping reaction mixture (10 µl) consisted of 1 × buffer (Promega, Germany), 2 mM MgCl2, 1 µl multiplex primer mix [final concentrations per forward and reverse primers varied and were either 0.25 µM (INRA6, T115, T40, C180, C105, C229), 0.3 µM (T107, BM757), 0.5 µM (C143, Haut14), 1 µM (BM1818), 3.5 µM (CSSM14), 4 µM (FSHB), or 6 µM (ILSTS06)], 150 ng DNA, 0.25 U GoTaq polymerase (Promega, Germany) and 5.2 µl A.dest. (sterile). Cycling conditions were the same for all four multiplex mixes: 95 °C 5 min, 5x (95 °C 30 s, touchdown beginning at 63 °C, with a decrease of 2 °C per cycle down to 55 °C 90 s, 72 °C 30 s), 40x (95 °C 30 s, 55 °C 90 s, 72 °C 30 s), final extension at 60 °C for 30 min. Size of amplicons was determined by calibration using the GeneScan™ 500 ROX™ size standard. Separation of fragments was carried out on an A3130xl automated capillary sequencer using the software GeneMapper v.3.7 for allele scoring (all Applied Biosystems).

Genotyping error estimation

Microsatellite amplicons were screened for genotyping errors (large allele dropouts, stutter bands) and probability of null alleles being present using MICRO-CHECKER (version 2.23, Van Oosterhout et al. 2004). We tested all loci across all populations for consistent patterns of deviation from Hardy–Weinberg expectations (HWE) using GENEPOP (version 4.5.1; Rousset 2008). All pairs of loci were further checked for linkage disequilibrium within all sampling units applying the algorithms implemented in GENEPOP and ARLEQUIN (version 3.5; Excoffier et al. 2005) including Bonferroni correction for multiple comparisons (Rice 1989). Additionally, we calculated the number of identified alleles and estimated expected and observed heterozygosities as well as the polymorphic information content (PIC) for each marker using the adegenet R package (Jombart 2008). Monomorphic markers were excluded from further analyses.

Estimating genetic diversity

All statistical analyses were performed using the R environment (R Core Team 2017). We assessed the amount of genetic variation within each AMU by estimating expected and observed heterozygosities (HE, HO), allelic richness (AR) and the degree of heterozygote deficiency (FIS) in each management unit. Estimation of AR was based on rarefaction to correct for the smallest sample size (n = 12). Confidence intervals for AR and FIS metrics were obtained using bootstraps with 999 replications. All metrics were estimated applying the diveRsity package (Keenan et al. 2013). We estimated effective population sizes (NE) for all administrative management units using the NeEstimator v2 software (Do et al. 2014). NE values were based on the linkage disequilibrium method with bias correction developed by (Waples and Do 2008). The same critical thresholds (0.05, 0.02, 0.01) as in Zachos et al. (2016) were applied to correct for linkage of rare alleles with frequencies below these values. The NFL unit was excluded to avoid any potential bias in population size estimates due to low sample size below 15 individuals (Do et al. 2014).

Estimating genetic structure

We assessed genetic structure at the level of the AMUs based on pairwise FST values (Wright 1965) as well as the pairwise Jost's D metric (Jost 2008) using the strataG R package (Archer et al. 2017). While Jost’s D is more appropriate for quantifying genetic (allelic) differentiation of populations showing varying levels of genetic diversity, FST better reflects past demographic processes and fixation (Whitlock 2011; Jost et al. 2018). Significance of differences in pairwise comparisons was estimated with 9999 replications and subsequent Bonferroni correction.

To assess whether AMUs actually constituted genetically separate clusters, we applied a Bayesian clustering approach. Specifically, we used the program STRUCTURE (version 2.3.4, Pritchard et al. 2000) and tested for the presence of genotypic clusters (K), with the number of possible clusters ranging between K = 1 and K = 14, using an admixture model and correlated allele frequencies. After having checked for the likelihood to have converged, we estimated the probability for each K-value in five independent runs with 500,000 iterations as burn-in followed by 1,000,000 MCMC iterations. The optimal number of K was determined using log-likelihood plots and the ΔK method by Evanno et al. (2005) implemented in the STRUCTURE Harvester platform (Earl and von Holdt 2012). Individual likelihoods of cluster memberships (q) were averaged over the five runs using the CLUMPAK online program (Kopelman et al. 2015).

We used STRUCTURE in a hierarchical framework by re-running the clustering algorithm for each of the detected genetic clusters in the previous analysis (Coulon et al. 2008; Balkenhol et al. 2014). The procedure was repeated until the optimal number of inferred genetic clusters was equal to one (K = 1). By doing this, subtle structuring is more likely to be detected because the largest break in the dataset is reiteratively removed so that this strong signal does not blur a weaker signal at lower hierarchical levels (Janes et al. 2017). We performed the hierarchical STRUCTURE analysis with ‘sampling location’ (i.e., the AMU) as a prior (locprior; Hubisz et al. 2009). All AMUs from Schleswig–Holstein and the reference area from Denmark were included in this analysis, as these are the sampling areas among which gene flow can be substantial enough to form actual genetic clusters. MWP samples were from regions not directly in close vicinity to our study area and were hence excluded from the STRUCTURE analysis as well as the migration analysis.

Estimating directional migration rates

Relative, directional migration was estimated using the divMigrate method (Sundqvist et al. 2016) which is implemented in the diveRsity R package (Keenan et al. 2013). While other, more complex algorithms are available for estimating asymmetric migration rates (e.g., BayesAss, Rannala 2007; MIGRATE-N, Beerli 2004), we chose divMigrate (Sundqvist et al. 2016) because it can be calculated from standard measures of genetic differentiation and does not require multiple additional parameters to be estimated (Sundqvist et al. 2016). The method tests for significant directionalities in gene flow between pairs of populations based on asymmetric distributions of allele frequencies and generates an output with relative migration rates scaled to values between 0 and 1 (Sundqvist et al. 2016).

We chose the GST measure of genetic differentiation (Nei 1972) from the options provided by divMigrate since it is similar to the FST values applied above (Whitlock 2011). Again, the analysis was performed for all AMUs within the study area as well as the Danish reference population, which we included because of suspected ongoing migration from Denmark into Germany. Based on the pairwise migration rates, we calculated the mean immigration (I) and emigration (E) rates as well as their ratio (RI/E) for each AMU. RI/E > 1 would indicate that the rate of immigration in a population is higher than the emigration rate and vice versa for RI/E < 1. Finally, we note that the results of the divMigrate analysis do not necessarily represent actual migration but rather estimate the probability of the exchange of genes between two sampling locations (Marrotte et al. 2017; Bohling et al. 2019). Further, relative migration rates are estimated across all pairs of included populations and do not account for spatial context or distance between them.

Modeling of genetic patterns

In the next step, we used regression modeling to correlate genetic variation within and among the AMUs with available ecological and environmental information. Specifically, we tested whether genetic diversity, differentiation and migration rates could be explained by local population sizes (Si) or densities (Di) within each AMU i, or as a function of the cumulative sizes \(\left( {\sum\nolimits_{j = 1}^{n} {S_{j} } } \right)\) or cumulative densities \(\left( {\sum\nolimits_{j = 1}^{n} {D_{j} } } \right)\) of the three AMUs j (j = 1–3) closest to the focal AMU i. The first two indices, Si (number of individuals in AMU i) and Di (individuals per hectare in AMU i), assume that genetic patterns and migration are only influenced by local population characteristics (i.e., size or density). In contrast, the latter two indices essentially are metrics used to describe isolation of multiple, potentially connected populations, and assume that the existence of large or densely populated neighboring AMUs is important for explaining observed population structure (e.g., Balkenhol et al. 2013). The three closest AMUs were chosen to calculate the connectivity indices because this included, in all cases, all the directly neighboring management units that could potentially exchange dispersing individuals with the focal unit.

We used officially available population size estimates (Meißner et al. 2008; Ministerium für Energiewende, Landwirtschaft, Umwelt und ländliche Räume des Landes Schleswig–Holstein 2012; Reinecke et al. 2013) for each AMU to represent S, and estimated D by dividing population size by the area of potential red deer habitat in the AMU. Potential habitat for each AMU was based on official thematic landscape data (authoritative topographic cartographic information system, ATKIS) and included all patches of forest, heathland and moors within the range of each AMU (Reinecke et al. 2013).

We then modelled genetic diversity (AR), genetic differentiation (Jost’s D) and mean immigration (I) as well as emigration (E) rates as a function of the four different indices, as well as a null model (intercept-only). We chose AR as a measure of genetic diversity as it was corrected for varying sampling sizes across AMUs. Similarly, we chose Jost’s D as an estimate of genetic differentiation because it measures the fraction of allelic variation among populations and thus accounts for varying genetic diversities within AMUs (Jost et al. 2018). Finally, we chose immigration and emigration rates as measures of directional dispersal. To compare models, we used an information-theoretic approach based on Akaike’s Information Criterion corrected for small sample size (AICc; Akaike 1973; Burnham and Anderson 2002). The model with the lowest AICc value was deemed best, but models with ΔAICc ≤ 2 were considered equally plausible (Burnham and Anderson 2002).

Genetic drift and isolation by distance

Following Jordan and Snell (2008) we tested for the potential effect of drift in isolation assuming that historic drift as represented in low genetic variation in smaller populations caused higher levels of differentiation. Therefore, we expect to see a negative relationship between the mean pairwise FST values of each AMU with all other AMUs and their expected heterozygosities HE (i.e., AMUs with larger FST should show lower HE values than AMUs with smaller FST-values). We further correlated mean pairwise FST values with allelic richness (AR) as the predictor variable (Whiteley et al. 2010; Funk et al. 2016).

Finally, we tested for isolation-by-distance (IBD; Wright 1943) using a Mantel test between genetic distances (linearized FST, i.e., FST/1 − FST; and Jost's D values) and the natural log of the geographic distance among AMUs (Slatkin 1993). A significant IBD pattern in both FST and Jost’s D would indicate that gene flow occurs among AMUs but is spatially limited, which hints at subpopulations connected via dispersal (Hutchison and Templeton 1999; Aguillon et al. 2017).

Results

We excluded 65 samples from further analyses because of insufficient numbers of successfully sequenced loci (≤ 11 markers). Therefore, the final dataset consisted of 480 samples including 68 individuals from the two reference regions located in Mecklenburg-Western Pomerania (46 samples) and Denmark (22 samples; Table 1). Two (T40, C105) of the original 14 microsatellite markers were dropped as they had only two alleles and were near monomorphic in the vast majority of samples, with frequencies below 0.15 observed for one of the two alleles. The number of alleles of the remaining markers ranged between three and 14. Polymorphic information content ranged from 0.3 up to 0.86 with a mean PIC of 0.62 (SD 0.2) across all loci (more information on marker diversity is provided in the Supplement; File S2). None of the retained markers showed issues with null alleles or consistent deviations from HWE. We found no evidence for significant linkage for any of the compared pairs of loci across all sampling units. Private alleles were detected within samples from one reference area (MWP: three alleles) and from two management units (ILO and NFL one allele each).

Genetic diversity

We observed a mean expected heterozygosity of 0.59 (SD 0.04) and a mean allelic richness of 4.20 (SD 0.47) alleles with a minimum of 3.41 and a maximum of 5.12 alleles (based on 12 diploid individuals, see Table 1). The Hasselbusch administrative management unit (HAB) showed the lowest values regarding these two metrics. Samples from the two reference areas differed with regard to their genetic diversity with Denmark showing the lowest values of HE and AR (Table 1). The samples from Mecklenburg-Western Pomerania actually exhibited the highest estimates for all diversity metrics compared to DK and the AMUs from Schleswig–Holstein. We found no indications of significant heterozygote deficiency. With the exception of Barlohe (BAL) and Schierenwald (SCW), confidence intervals of all estimated FIS values were low and overlapped with zero (Table 1), conforming with expectations for random mating within AMUs.

Genetic structure

We observed a global fixation (FST) value of 0.09 and a global Jost’s D of 0.12 across all 12 AMUs of Northern Germany (p < 0.0001 for both values). Pairwise estimates of FST ranged between 0.01 up to 0.20 and Jost's D from 0.01 to 0.23 (Table 2). Overall, estimates of the two metrics agreed in most cases regarding the significant differentiation between the considered AMUs. However, not all AMUs were genetically differentiated. We were able to distinguish three groups of administrative management units that did not show significant structuring for both estimates. The first consisted of BAL, ILO and SCW, the second one included the three AMUs from the Lauenburg area (LAW, LAS and LAE), and the third group included NFL and DK, where the lowest level of differentiation was observed (FST: 0.01; Jost's D: < 0.01). In some pairwise comparisons, Jost’s D estimates differed from FST values, e.g., for the SCW-HAB value comparison, which indicated intermediate FST values, but very low differentiation based on Jost’s D (Table 2).

Table 2 Pairwise estimates of FST (above diagonal) and Jost's D values (below diagonal)
Full size table

We observed a complex, hierarchical genetic structure of three different levels based on the STRUCTURE analysis. Using the ΔK method (Evanno et al. 2005) the optimal number of genetic clusters K at the first level was three, essentially dividing the individuals into a northern (Cluster North), central (Cluster Center) and southern (Cluster South) group of origin (Fig. 2). The northern as well as the southern cluster was again split into another two subgroups whereas the central cluster comprised three different genetic groups at the second hierarchical level (Fig. 3). Finally, we only found additional substructures of K = 3 at the third level for one of the two southern clusters (South 2, Fig. 3; see also supplemental file S4). The majority of individuals were clearly assigned to the different clusters with high ancestry values (q) above 0.7 (Supplemental File S4).

Fig. 2
figure 2

First level of the hierarchical STRUCTURE analysis: the number of genetic clusters (K = 3) was determined with the Evanno method. Probabilities of group membership (Q-values) are presented for all individuals from the AMUs in Schleswig–Holstein and the reference area from Denmark (DK)

Full size image
Fig. 3
figure 3

Results of hierarchical STRUCTURE analysis. Upper part shows partitioning among clusters. The map presents the final results for all MUs showing the proportions of the most likely origin of the sampled individuals. The overall sample size of this analysis is 434 out of 480 individuals. Samples from MWP (n = 46) were excluded since they originated from regions not directly neighboring the study area

Full size image

Directional migration

Based on the divMigrate analysis, we observed variation in directionality and degree of gene flow among the AMUs and between some AMUs and the Danish reference area. Estimated rates of relative gene flow ranged from 0.04 up to 1 with an average of 0.15. A pairwise matrix with all directional estimates of gene flow is provided in the Supplement (File S5). We observed the highest rates of directional gene flow (> 0.2) between AMUs in the southeastern region (LAW, LAE, and LAS) as well as the central region (BAL, ILO, SCW; Fig. 4). The results further suggested that gene flow was more likely to occur from DK towards AMUs in the southern regions (e.g., the Lauenburg management units) than vice versa.

Fig. 4
figure 4

Direction and magnitude (indicated by arrow thickness) of estimated gene flow between management units (sources) based on the divMigrate analysis. Only results with relative migration rates above the average of 0.2 for the south-eastern region of Schleswig–Holstein without the DK reference area are shown here for illustrative purposes

Full size image

Differences with regard to directed migration rates were also detected by mean immigration and emigration rates (Table 3) with several AMUs either exhibiting similar rates of emigration (LAW, LAS, SAW) or very low values of overall gene flow (DUV, NFL). The HAB administrative deer management unit exhibited one of the lowest migration ratios (RI/E = 0.63), together with the reference area from Denmark (RI/E = 0.57). With an RI/E value > 1, local deer populations in five out of 13 AMUs potentially received more migrants than they produced (DUV, LAE, LAS, ILO, SEG), while the eight remaining ones contributed more migrants than they received (RI/E < 1:BAL, ELD, HAB, LAW, NFL, SAW, SCW and DK; Table 3).

Table 3 Mean immigration (I) and emigration (E) rates as well as their ratio (RI/E) estimated for all administrative deer management units in Schleswig–Holstein
Full size table

Influence of population size and neighboring population densities

In our regression analysis (Table 4), allelic richness was best explained by the neighboring population densities (adj. R2 = 0.57, p = 0.003); i.e., AR within AMUs increased with higher cumulative densities of red deer in the neighboring management units (Fig. 5). Mean emigration rates (E) were best explained by both neighboring population size and by the density of the neighboring management units (Table 4). Neither mean Jost’s D nor mean immigration rates (I) were explained by any of the indices.

Table 4 Meta-population study linking genetic metrics of diversity, differentiation and gene flow with estimates of meta-population structure
Full size table
Fig. 5
figure 5

Linear regression model showing the significant increase (adj. R2 = 0.57, p = 0.003) of allelic richness (AR) with higher cumulative densities of neighboring deer management units (deer per 100 hectares)

Full size image

Genetic drift and isolation by distance

Genetic differentiation based on mean pairwise FST values was negatively correlated with higher estimates of genetic diversity. This indicates that drift is potentially influencing genetic diversity and could drive divergence between AMUs in our study area. For example, we observed the highest r2 score of 0.73 (p < 0.001) between FST and expected heterozygosity (Fig. 6). Allelic richness also significantly decreased with higher values of mean FST (r2 = 0.58, p = 0.004). Further, we detected effects of spatially limited gene flow and significant isolation by distance. Results of the Mantel analyses indicated significant IBD among AMUs using both linearized FST (r = 0.42; p = 0.003) and Jost’s D (r = 0.28; p = 0.015), respectively.

Fig. 6
figure 6

Scatterplot showing the significant decrease of mean pairwise FST values and genetic diversity of administrative deer management units in Schleswig–Holstein based on allelic richness (AR) Results are based on a linear regression model (r2 = 0.58, p = 0.004)

Full size image

Discussion

We investigated the genetic structure and differentiation of administrative management units to determine whether the practice of managing the local red deer populations as separate populations is effective, or if future management should account for substructures and genetic exchange among them. This is particularly relevant when populations are low in abundance and experienced different types of restrictions in the past as in the presented study. We found that limited gene flow most likely caused by anthropogenic fragmentation and management goals (culling of individuals outside of designated deer areas) has led to genetic drift and decreased genetic diversity.

The low values of HO observed in our study are rare for red deer and usually only found in populations with long-term low effective population sizes such as the red deer from Sardinia or from Mesola in northern Italy (Hmwe et al. 2006). Both, HO and HE values, in the Hasselbusch AMU are among the lowest ever found in a population of this species (Zachos and Hartl 2011; Zachos et al. 2016). We note that comparisons of genetic diversity estimates across studies should be performed with caution (Reiner et al. 2019). Although our data set shares only four loci with the most comprehensive microsatellite study of red deer in Europe to date (Zachos et al. 2016), our study nevertheless shows that there is a clear trend towards low genetic diversity in red deer from Northern Germany.

A very similar pattern can be seen in NE values of all AMUs (Table 1). Although the NE values for several of the northern German deer AMUs are within the range of reported values from other European populations, many of the Schleswig–Holstein populations, again including HAB, are clearly at the lower end and below the effective population size threshold of 50 individuals, a value below which inbreeding depression is likely to occur (Frankham et al. 2010). However, observed FIS values were quite low (Table 1) with no clear signs of heterozygote deficiency and fixation. We assume that existing gene flow at short ranges seems to compensate for drift effects on genetic diversity in some cases. This assumption is supported by significant isolation by distance which indicates that drift and gene flow are in equilibrium at regional scales (Hutchison and Templeton 1999; Jordan and Snell 2008).

Ultimately, the lack of restrictions to dispersal, as experienced by red deer in our study area, should result in higher levels of genetic diversity. This was confirmed by relatively high values of genetic variability (HE, NE) in the reference population from MWP where red deer have not been restricted to declared red deer zones in the past and range throughout the state in higher abundance compared to Schleswig–Holstein (Kinser et al. 2010). Comparable values of genetic diversity were also confirmed for the three AMUs of the Lauenburg area (LAE, LAS, LAW) which could also be explained by the relatively large population size in Lauenburg in combination with gene flow from the east (Mecklenburg-Western Pomerania).

Furthermore, we observed a significant decrease in differentiation (FST) with higher levels of genetic diversity (HE, AR). This could be due to potential effects of historic drift in isolation and small populations which is expected to be the predominant cause for genetic differentiation (Jordan and Snell 2008; Whiteley et al. 2010; Funk et al. 2016). Only AMUs in close proximity did not exhibit any significant values of differentiation based on both FST and Jost’s D estimates (e.g., BAL, ILO and SCW or the Lauenburg populations). In a few instances we observed inconsistencies between the estimates of Jost’s D to FST (e.g., NFL and HAB and SCW not being significantly differentiated based on Jost’s D; Table 2). There is an ongoing debate on which measure is superior compared to others (Whitlock 2011; Leng and Zhang 2011; Verity and Nichols 2014). However, Meirmans and Hedrick (2011 suggest to use both measures complementary to each other as FST can describe demographic changes whereas Jost’s D infers actual differentiation.

Hierarchical structure and gene flow

We observed a hierarchical genetic structure comprising three main clusters: North, Center and South. The first cluster was located north of the Kiel Canal (an effective barrier to deer dispersal due to the steep embankments) and also included the Danish red deer. The small population of the NFL management unit was founded by red deer individuals dispersing from Denmark into northern Germany. The assignment of ELD to the northern cluster was surprising because its founders came from BAL crossing the Kiel Canal in the late 1960s (when the embankments were not yet in their present state; Meißner et al. 2008). However, low population size and genetic drift have apparently resulted in divergence from the central cluster, which is located just south of the canal. In addition, immigration from Denmark into the ELD population has been shown by means of genetic data (mtDNA sequences; Reinecke et al. 2013). Therefore, it seems in accordance with the population’s history to consider ELD as a separate northern sub-cluster.

The sub-structure of the central cluster can also be explained based on geography and historic background. The three sample sites just south of the Kiel Canal (BAL, ILO, SCW) have historically been separated from the ones further south because of limited dispersal of red deer outside the ‘designated red deer zones’ (Wotschikowsky 2010; Reinecke et al. 2013). The HAB population was founded by dispersed individuals from SEG but a fenced highway has prevented any potential migration between these two AMUs. Due to its low census and effective sizes, drift has been high in HAB, which is mirrored by its substantial differentiation from SEG today. The separate status of Duvenstedt (DUV) is understandable since the population is not native but was founded with red deer from other parts of Europe (Jessen 1988; Meißner et al. 2008). In the southeast, LAS is separated from LAW and LAE by a highway, but this is a relatively recent barrier (completion during the 1990s), obviously not yet reflected in the gene pools on either side. Interestingly, the SAW population comprises three different subclusters, two of which were only found there (South 2b and 2c). A red deer hunting enclosure located in that area perhaps suggests a similar historical development as with the DUV population. Individuals from other parts of Europe introduced into the private enclosure could potentially have escaped the fenced area and established themselves within the local population (cf. Frantz et al. 2017).

We observed high levels of gene flow between AMUs with low differentiation (Fig. 4), which were again the complex consisting of BAL, SCW, and ILO, as well as the local deer populations from the Lauenburg area in the south-eastern region of Schleswig–Holstein. Overall, diversity within and gene flow among AMUs was best explained by size and density of the surrounding local populations as our modeling analyses (Table 4) showed. Deer populations that were adjacent to larger or higher-density populations had higher rates of gene flow and higher levels of diversity.

Most of the AMUs with lower mean rates of immigration as compared to emigration rates (Table 3) are characterized by either small population sizes, low densities, or higher levels of isolation, suggesting that they probably received fewer genes from other populations in the past. Relative to the other populations they are therefore more likely to send out individuals (Bohling et al. 2019). Genetic similarity, for example due to historical reasons, will also lead to positive values of inferred migration. Since the HAB population was founded by migrants from SEG, migration values are not zero and reflect population history. This is in accordance with the Duvenstedt unit (DUV) showing no signs of migration to or from other populations because it was founded with non-native deer (Jessen 1988). Still, anecdotal reports of dispersing red deer further support the conclusion that there is some level of gene flow (Reinecke et al. 2013). Single individuals have been seen outside established population ranges, the Lauenburg red deer are known to be in contact with the neighboring populations in Mecklenburg-Western Pomerania to the east, and in 1986 and 1987, single stags migrated from Hasselbusch (HAB) to Barlohe (BAL) and from Duvenstedt (DUV) to Segeberg (Jessen 1988; Peters 2000; Zachos et al. 2007; Meißner et al. 2008). The latter is also supported by the results of the STRUCTURE analysis. Whether they successfully reproduced in SEG, however, is unknown.

Within the last decade red deer from Denmark have established themselves south of the German border and are increasing in numbers. We were able to detect first signs of genetic exchange between the NFL/DK population and the ELD management unit. This shows the high potential of the species to migrate throughout the state and establish new ranges.

Management Implications and Future Research

In summary, based on our analyses of genetic structure and gene flow, we were able to distinguish two major groups of AMUs which essentially represent single GMUs. First, in the central part of Schleswig–Holstein, the three AMUs of BAL, ILO and SCW form one genetically distinct cluster. Secondly, the same holds for the AMUs in the Lauenburg area (LAE, LAS, LAW) in the south-east of the state. This indicates a discrepancy between the current administrative delineation of management units and actual levels of genetic exchange among these areas (see also Fig. 4). Our results also show that observed genetic patterns (diversity and gene flow) in a local deer population are largely explained by the densities of populations in its close vicinity. Local management decisions that change local abundance could have genetic impacts not only on the local population but also on neighboring AMUs, especially if AMUs are interpreted as single GMUs when they are actually well connected to others. Therefore, future management of red deer populations in Schleswig–Holstein needs to incorporate parameters such as deer population sizes and habitat availability for neighboring administrative MUs. Data on dispersal or gene flow and population structure derived from genetic studies like ours should ideally be incorporated in current management practices or when new units for wildlife management are spatially delineated (Paetkau 1999; Lowe and Allendorf 2010).

Another important factor is temporal change in age- and sex-structure of the local populations (Langvatn and Loison 1999; Marjamäki et al. 2013). Recording these parameters could help to gain a better understanding of potential source-sink dynamics (Draheim et al. 2016). In particular, younger males are more likely to disperse at higher local densities (Loe et al. 2009). Therefore, future research should also focus on the proportions of young males in local populations and how density-dependent dispersal could potentially influence gene flow and the genetic differentiation of the subpopulations. For example, estimating dispersal among localities using capture-mark-recapture (Graves et al. 2014) or telemetry (Zeller et al. 2018) could be applied to assess the demographic effects of inter-population movements.

In particular, the exchange of individuals between isolated populations such as HAB needs to be enhanced in the near future to counteract the continuing loss of genetic diversity. HAB is not far away (approximately 10 km) from the larger GMU formed by SCW, ILO and BAL. Still, we observe high levels of differentiation and hardly any gene flow. The STRUCTURE analysis assigned one single individual sampled in HAB to the cluster of SCW, ILO and BAL (Fig. 2 and Supplement S4). Similar patterns can be observed for DUV and SEG which are also not far apart (ca. 15 km) but only two individuals sampled in DUV were assigned to the SEG cluster (Supplement S4). This leads to the conclusion that landscape characteristics between AMUs affect the genetic exchange among them and thus influence size and density of populations; we will need further analyses to identify landscape features that facilitate or impede natural dispersal among AMUs. Based on the results, migration corridors and locations for crossing-structures (e.g., green bridges) can then be identified to mitigate the effects of barriers and landscape resistance on the migratory movements of red deer (Chetkiewicz et al. 2006; Shirk et al. 2010).