Genetic structure of recently fragmented suburban populations of European stag beetle

Abstract Habitat loss and fragmentation due to urbanization can negatively affect metapopulation persistence when gene flow among populations is reduced and population sizes decrease. Inference of patterns and processes of population connectivity derived from spatial genetic analysis has proven invaluable for conservation and management. However, a more complete account of population dynamics may be obtained by combining spatial and temporal sampling. We, therefore, performed a genetic study on European stag beetle (Lucanus cervus L.) populations in a suburban context using samples collected in three locations and during the period 2002–2016. The sampling area has seen recent landscape changes which resulted in population declines. Through the use of a suite of F ST, clustering analysis, individual assignment, and relatedness analysis, we assessed fine scale spatiotemporal genetic variation within and among habitat patches using 283 individuals successfully genotyped at 17 microsatellites. Our findings suggested the three locations to hold demographically independent populations, at least over time scales of relevance to conservation, though with higher levels of gene flow in the past. Contrary to expectation from tagging studies, dispersal appeared to be mainly female‐biased. Although the life cycle of stag beetle suggests its generations to be discrete, no clear temporal structure was identified, which could be attributed to the varying duration of larval development. Since population bottlenecks were detected and estimates of effective number of breeders were low, conservation actions are eminent which should include the establishment of suitable dead wood for oviposition on both local and regional scales to increase (re)colonization success and connectivity among current populations.


| INTRODUC TI ON
Land use changes such as urbanization have significantly driven habitat loss and population fragmentation. Such fragmentation can reduce gene flow among populations and decrease population sizes. This can in turn increase the rate at which genetic variation is lost by genetic drift and the negative effects of inbreeding depression on survival and reproduction, ultimately increasing the risk of population extinction (Allendorf et al., 2013). In contrast, population viability may be enhanced by connectivity to surrounding populations. Inference of patterns and processes of population connectivity derived from spatial genetic analysis has proven invaluable for conservation and management (van Strien et al., 2014).
However, a more complete account of population dynamics may be obtained by combining spatial and, often neglected, temporal sampling. Firstly, where genetic differences are low, temporal stability of the differentiation adds confidence that this reflects a true signal (Waples, 1998). Second, the use of temporal samples is one of the best ways to estimate the crucial demographic factor effective population size. Finally, analysis of temporal samples permits recognition of population turnover (i.e., when local populations become extinct, but the habitat patch is recolonized). Discriminating between situations whereby genetic differences reflect permanent population isolation or metapopulation dynamics (i.e., dynamic balance between migration, extinction, and recolonization) is important for judicious conservation, particularly for species with low colonizing ability (Holyoak & Ray, 1999).
The European stag beetle, Lucanus cervus L. (Coleoptera: Lucanidae; Figure 1), is a saproxylic beetle attributed with limited vagility. Although it is distributed widely across Europe, a decline is presumed in many countries and regions, and populations have gone extinct or are threatened at the northern edge of its range (Harvey et al., 2011). Accordingly, the species has been designated "near threatened" in the European Red List (Cálix et al., 2018) and is protected by the Habitats Directive of the European Union since 1992 and by the Bern Convention since 1979 (Luce, 1996). Habitat loss and fragmentation have been identified as the main threats for this species (Della Rocca et al., 2017). Moreover, habitat continuity was found to be of major importance to maintain its populations (Thomaes, 2009;Thomaes et al., 2018).
Emergence can be observed from May to June, and adults are active for a few weeks up to three months (Harvey et al., 2011). Although the life cycle of stag beetle suggests its generations to be discrete, larval development, especially the third instar stage, can vary in duration, even locally due to the microclimate and habitat quality.
So far, evidence showed that dispersal distances in the species seem to range from a few hundred meters up to a maximum of five kilometers (Rink & Sinsch, 2007). As female stag beetles usually stay close to their site of emergence, males are considered the main dispersers (Rink & Sinsch, 2007;Tini et al., 2017Tini et al., , 2018. Female dispersal distances become even smaller when suitable oviposition sites are abundant and nearby, according to Tini et al. (2018). However, in suburban areas and under less suitable habitat conditions, dispersal distances can increase (Rink & Sinsch, 2007;Thomaes et al., 2018).
Sex-biased dispersal may be scale dependent, since long-distance and short-distance dispersal are likely to be caused by different reasons. While inbreeding is likely avoided by dispersal on a local scale (e.g., Lebigre et al., 2010), dispersal over a longer distance could be induced by the need to colonize new sites or by high local competition or predation risk (Lawson Handley & Perrin, 2007).
Dispersal behavior of stag beetle has to date been studied primarily using radio-telemetry and mark-and-recapture methods (e.g., Chiari et al., 2014;Hawes, 2008;Rink & Sinsch, 2007;Thomaes et al., 2018;Tini et al., 2018). These methods, however, are limited due to low sample sizes, short time windows, and difficulties to register long dispersal events. Furthermore, they do not monitor effective dispersal (i.e., dispersal followed by reproduction). Genetic markers have proven useful for the monitoring of "real-time" dispersal and gene flow over space and time in other systems but have not been applied to many saproxylic species, especially within a fine scale habitat fragmentation context. The objective of this study was to assess population genetic structure and sex-specific dispersal among stag beetle within a region that has undergone recent habitat fragmentation linked to urbanization. Until the 60s, stag beetle was quite common in the surroundings of Brussels, capital of Belgium (Thomaes et al., 2008(Thomaes et al., , 2010. As in other regions of NW Europe (Rink & Sinsch, 2006), the relict populations now occur in more suburban habitats such as gardens, parks, tree-lined lanes, afforested slopes, and edges of large woodlands, such as Sonian Forest (Thomaes et al., 2010). Due to the species' limited dispersal capacity, we expected to find evidence of restricted gene flow. Analysis of temporal samples was used to investigate if spatial structure was temporally stable, indicating demographically independent populations. Alternatively, genetic change over time could point to demographic instability which may F I G U R E 1 Male European stag beetle at Jagersveld, Watermaal-Bosvoorde (Belgium) include population turnover events. We tested for the occurrence of recent genetic bottlenecks and reduced effective population sizes that could be considered coincident with census population declines.
Based on the results, implications for conservation were discussed.

| Study area and sampling
The suburban municipalities Beersel, Watermaal-Bosvoorde, Overijse, and Sint-Genesius-Rode hold the remaining populations of Lucanus cervus south of the city of Brussels (Belgium) of a once more continuous population from Halle to Leuven (35 km) (Thomaes et al., 2010). This area has been well investigated and most if not all remaining populations are well known (Thomaes, 2009;Thomaes et al., 2010, R. Cammaerts, unpublished data). Jagersveld situated in Watermaal-Bosvoorde was visited twice per day (by RC), and other sites were visited at least once per week (by AT and a network of volunteers) between the end of May and mid-September. During  Table 1). For Watermaal-Bosvoorde, we chose samples taken three years apart starting from 2002 until 2008 (Table 1), which largely concurs with the generation time of stag beetle in the area (Rink & Sinsch, 2008;Smit & Hendriks, 2005). For Overijse, we selected samples taken one to two years apart (Table 1). This variation in number of years between sampling enabled us to test if there is also a temporal genetic structure. Only 10 samples could be collected in F I G U R E 2 Sampling locations of European stag beetle. Top: The inset shows a regional map of Belgium with the study area indicated with a blue rectangle; the red dots are the sampling locations. Middle: number of females, males, and samples of unknown sex depicted as pie charts. Bottom: number of samples per sampling year depicted as pie charts. The size of the pie charts is relative to the total. Resolution Layers of 2015 with a 20m resolution (https://land.coper nicus. eu/pan-europ ean/high-resol ution -layers) Sint-Genesius-Rode. A single sample was also collected at the edge of the Sonian Forest near Hoeilaart in 2015. Two small populations in Beersel were not included in this study, located 350 m north and 2 km west of the Sint-Genesius-Rode population. Also, a few sites that lie in close vicinity of Jagersveld (Watermaal-Bosvoorde) have not been sampled. When determinable, the gender of each sampled stag beetle was recorded. The samples were either stored in a freezer (−20°C) or preserved in absolute ethanol and stored in a refrigerator at c. 4°C. Some of the samples were also analyzed for the phylogeographic study of Cox et al. (2019).

| DNA extraction and genotyping
The tissue used for DNA extraction depended on availability, but was mostly muscle tissue from legs.  Lcerv-3, Lcerv-4, Lcerv-6, Lcerv-7, Lcerv-8, Lcerv-9, Lcerv-16, Lcerv-17, Lcerv-20, Lcerv-21, Lcerv-25, Lcerv-28, Lcerv-29, Lcerv-30, Lcerv-31, and Lcerv-36, described by McKeown et al. (2018). We included the primer sets in four multiplex PCRs and one simplex PCR as described by Cox et al. (2019) and performed the PCR reaction under the same conditions. Genotyping analysis was achieved using an ABI 3,500 Genetic Analyzer (Applied Biosystems) and with the GENEMAPPER v.4.0 software package. To test for reproducibility, 7% of the samples were blindly replicated two to five times within and across well plates. Samples with fewer than twelve scored loci were discarded. To investigate possible deviations from Hardy-Weinberg equilibrium (HWE), we used the test available in the program GENEPOP 4.6 (Rousset, 2008). GENEPOP was also used to assess the presence of null alleles with the maximum likelihood method following the expectation maximization (EM) algorithm of Dempster et al. (1977), and a second software program, ML-NULLFREQ (Kalinowski & Taper, 2006). GENEPOP was further used to test for linkage disequilibrium (LD) for each pair of loci. We conducted these tests for each population (location), for each population and year of sampling, and after excluding potential migrants (see further). We implemented a correction for multiple testing with TA B L E 1 List of sampling locations with year of sampling and associated population genetic statistics. When samples were discarded, the adjusted number of samples was given between parentheses the false discovery rate method (FDR) (Benjamini & Hochberg, 1995) with a nominal level of 5%. Before further analysis was conducted, replicate genotypes were detected because some samples may have been collected from the same individual due to the nature of the samples (prey rests). Only one genotype among replicates was kept in the final dataset.

| Spatiotemporal genetic structure and dispersal
We calculated the following estimates of genetic diversity using the R package hierfstat 0.04-22 (Goudet, 2005) for each population and for each year separately as well when sample size was at least seven: rarified allelic richness (AR), observed (H O ) and expected heterozygosity (H E ), and the inbreeding coefficient (F IS ). Number of private alleles was estimated using R package poppr 2.8.5 (Kamvar et al., 2014(Kamvar et al., , 2015. Calculations were conducted using R statistical language (R Core Team, 2019). Mean relatedness coefficients of Queller and Goodnight (1989) were also calculated for each spatial group or population within and among sampling years using COANCESTRY 1.0.1.9 (Wang, 2011). The significance of the comparison of relatedness within a sampling year with relatedness among a set of sampling years was assessed by drawing dyads ad random 1,000 times while keeping the properties of the original groups (with e.g., dyads be- Average relatedness and difference in relatedness among groups is then calculated each time, delivering a distribution with which the observed difference is compared. Genetic differences between samples were also quantified using the R package diveRsity 1.9.90 (Keenan et al., 2013) with F ST (Weir & Cockerham, 1984) and D est (Jost, 2008). Unbiased confidence intervals of 95% were calculated from 1,000 bootstraps. They were calculated among populations but also among sampling years of the same location, when sample size was at least seven.
Population structure was investigated using the Bayesian program BAPS v. 6.0 which allows the inclusion of individual geographic coordinates as a prior (Cheng et al., 2013;Corander et al., 2008). The program was run ten times for each value of K = 1-15. An admixture analysis (Corander & Marttinen, 2006) was performed using 100 iterations, a minimum of three individuals per population, 200 reference individuals for each population, and 20 iterations of reference individuals. We analyzed individual-based genetic structure with a principal component analysis (PCA) summarizing the allele frequency data. This was followed by a spatial PCA (sPCA), a useful method to find small scale spatial structure . Like PCA, this approach is independent from Hardy-Weinberg assumptions or linkage disequilibrium. This multivariate method uses allele frequencies and accounts for their genetic variability and spatial autocorrelation, calculated with Moran's I (Moran, 1948(Moran, , 1950 on the basis of a connection network. We used a distance based connection network, the neighborhood by distance graph, with a minimum distance d1 of 0 and a maximum distance d2 of 15 km to create a closed network. Moreover, the sPCA analysis was performed solely on the samples of Overijse where multiple sites could be sampled. This permitted a local spatial genetic analysis. In this case, the maximum distance of the network was 1 km. Global and local spatial structures (i.e., positive and negative spatial autocorrelation, respectively) were tested with 999 permutations and with randomly distributed allele frequencies as the null hypothesis. Analyses were performed with the R package adegenet 2.1.1 (Jombart, 2008). The spatial pattern of genetic variation was also investigated using spatial autocorrelation analyses on a local scale. We assessed the genetic similarity between pairs of individuals at different distance classes jointly for individ-  (Barton, 2020).
To characterize dispersal between sites assignment tests were performed in GENECLASS2 (Piry et al., 2004) to detect first generation migrants. We only included the data from 2008 and 2009 to be able to include all three populations and to stay as much as possible within one generation. Assignment tests were performed under the criteria L home /L max , the ratio of the likelihood of drawing that individual's genotype from the population in which it was sampled and the maximum assignment likelihood calculated for the individual considering all populations . This criterion holds when all subpopulations are sampled. Because there are two small subpopulations near Sint-Genesius-Rode not included in this study, the L home criterion could be better suited. The Bayesian method of Rannala and Mountain (1997) together with the resampling algorithm of Paetkau et al. (2004) was used to simulating 10,000 individuals (α = 0.01). We further used sibship analyzed with COLONY 2.0.6.5 (Jones & Wang, 2010;Wang, 2004;Wang & Santure, 2009) to discover potential dispersal events across years through the iden-

| Sex-biased dispersal
To assess if the dispersal patterns are different for males and females, we limited the data to those individuals of known sex (see Table 3 for numbers of males and females per population and year).
Again, we analyzed spatial autocorrelation for individuals sampled during 2007, 2008, and 2009 in Overijse with SPAGEDI as mentioned above, but now for males and females separately. Secondly, we used the following statistical descriptors as implemented in R package hierfstat: F IS , F ST , mean corrected assignment index (mAIc) and variance of corrected assignment index (vAIc) (Goudet et al., 2002). Due to mixture of resident and immigrant individuals in the sample resulting in a Wahlund effect, a larger heterozygosity deficit and higher F IS values are expected for the more dispersing sex (Goudet et al., 2002). F ST values on the other hand should be lower for the dispersing sex. The assignment index (AI) refers to the probability that a multilocus genotype or individual occurs within a sampling locality (Favre et al., 1997;Paetkau et al., 1995). The correction of the AI is performed by subtracting population means after log-transformation to remove population effect that may arise from different levels of genetic diversity in each population. This creates a distribution centered at zero (Favre et al., 1997), with positive values of AIc for resident individuals and negative values for potential dispersers (Goudet et al., 2002). The dispersing sex will have lower mAIc values and higher vAIc values than the philopatric sex.
Significance for the comparisons of all four statistics was assessed using 1,000 permutations. We performed the tests on all populations and on each population separately, except for Sint-Genesius-Rode due to the small sample size, using samples from all available sampling years, which could lead to using potentially temporally different populations. Also, this set of samples does not comply with the assumption of sampling individuals after dispersal events and before reproduction. We therefore repeated the analysis on partitioned data sets comprising subsets of years. Finally, we calculated relatedness among individuals and the difference in mean relatedness of each sex within locations using COANCESTRY as described above. It is expected that the group with the smallest relatedness coefficient is the dispersing sex.

| Effective population size and bottlenecks
There are several methods to calculate effective population size (N e ), each with their own assumptions. Because these assumptions are often violated or because they cannot all be tested for, it is generally recommended to use multiple methods to estimate N e . In this study, we used two single sample estimators and two temporal methods.
The sibship assignment method in COLONY, a single sample method, was used with the same settings as described earlier. In addition, the linkage disequilibrium method LDNE with a bias correction (Waples & Do, 2008) implemented in NeEstimator 2.01 (Do et al., 2014) was used, assuming random mating. A minimum allele frequency of 0.05 was chosen and jackknife-based corrected 95% confidence intervals were calculated. As temporal methods, we used the method developed by Jorde and Ryman (2007)

| Genotyping
Only one out of 149 samples from Overijse did not provide a genetic profile of good quality and was discarded (

| Spatiotemporal genetic structure and dispersal
Gene diversity was the lowest in Sint-Genesius-Rode and the highest in Overijse, although AR was quite similar in Sint-Genesius-Rode and Watermaal-Bosvoorde (Table 1) (Table 2).
The Bayesian, spatial clustering using BAPS and based on 16 microsatellites resulted in three groups, largely dividing the individuals among the sampling locations (Figure 4). Only a few individuals were assigned to another cluster than their sampling location, but there was hardly a sign of admixture, that is, apparently misassigned individuals were migrants rather than hybrids. No clustering according to sampling years was detected. When loci Lcerv-16 and Lcerv-30 were excluded, the optimal number of clusters became two as the separation of Sint-Genesius-Rode from Overijse resolved. However, differentiation between these sites was evident in F ST and sPCA based results.
The PCA results showed Watermaal-Bosvoorde to be distinct from the other two locations but with overlap, whereas Overijse and Sint-Genesius-Rode appeared to be very similar ( Figure 5). When subsets by location and year were analyzed as offspring with individuals collected in earlier years as potential parents, the same full-sib pairs were identified with a probability > .95 and none of the parental genotypes were assigned as mother or father of the investigated offspring.

| Sex-biased dispersal
Spatial autocorrelation analysis for samples collected during [2007][2008][2009] in Overijse showed no significant relationship between Moran's I and geographic distance for the female stag beetles (Figure 7). For males, spatial autocorrelation was significant in the third distance class, from 500 to 1,000 m. The regression slope on the logarithm of spatial distance was only significant for the males (b = −0.005, R 2 = 0.0002, P obs < exp =.0098), not for the females (b = −0.037, and mAIc was significantly lower for the females from Overijse in 2016 (p = .032) and for the females when analysing all samples together (p = .0404).
The relatedness among females (r = 0.062) was significantly lower than among males (r = 0.082) in Watermaal-Bosvoorde when all sampling years were analyzed together (p < .02). The difference in relatedness in Overijse for males and females was slightly higher

| Effective population size and bottlenecks
The single sample estimates of N e are given in

TA B L E 3
Results of the randomization tests for sex-biased dispersal the maximum likelihood approach, MLNe, to estimate N e over different time periods, the values for Overijse varied from 89 to 527, but always with the maximum possible N e as upper limit of the confidence intervals (Table 6). Gilbert and Whitlock (2015) found that while MLNe provides accurate estimates of N e , the coverage probability of the confidence intervals were inaccurate. When Overijse was assumed to be isolated, no estimate for N e below the assumed maximum N e was obtained. Using Watermaal-Bosvoorde as source population usually resulted in much lower estimates than when Sint-Genesius-Rode acted as source population. Estimates for Watermaal-Bosvoorde were stable among migration scenarios and were comparable to those obtained by the other methods, although N e estimated with temporal methods concerns a different time period than when using single sample methods.
A significant signal for a recent genetic bottleneck was found under the IAM model for all populations and years, except for Overijse in 2015 (Table 7). The Wilcoxon signed rank test performed using the TPM model was only significant for Watermaal-Bosvoorde. The allele frequency distribution showed a mode shift distortion in 2016 for the samples of Overijse which suggests the occurrence of a recent bottleneck.

| D ISCUSS I ON
There have been numerous studies of the effects of habitat fragmentation on genetic diversity in animals; however, this has been rarely studied among insects. In general, insects represent a group that has received scant study of fine scale spatiotemporal population processes (but see Bretman et al., 2011;López-Uribe et al., 2015;Melosik et al., 2020;Oleksa et al., 2015;Schauer et al., 2018;Zytynska et al., 2018). This study assessed fine scale spatiotemporal genetic variation within and among habitat patches that have been recently fragmented by urbanization using a suite of F ST , individual assignment, and kinship analyses. The results revealed numerically small to moderate but statistically significant genetic differentiation among patches. A key feature was that this differentiation was evident over multiple temporal samples and therefore indicates "biologically meaningful" restricted gene flow (Waples, 1998). In addition to restricted gene flow, assignment and kinship analyses confirm that contemporary dispersal between patches is highly restricted. Firstly, clustering analyses delineated groups that aligned with patch membership. Secondly, assignment analyses revealed low rates of migration. Collectively these results indicate that the recently fragmented areas now represent demographically independent populations from one another, at least over decadal time scales of relevance to conservation. Bottleneck and N e estimators suggest recent genetic declines and local populations to comprise small numbers of individuals that are successfully breeding. There were also indications of further spatial structuring within sites and social cohesion.

| Population genetic structure
Population genetic studies of saproxylic insects have typically reported weak genetic structuring over larger geographical ranges than studied here (Komonen & Muller, 2018 (Thomaes et al., 2010).
There was no clear temporal structuring for either of two patches Overijse and Watermaal-Bosvoorde. Genetic differentiation between years was largely nonsignificant. Likewise, Snegin (2014) found genetic similarity among groups of individuals sampled in con- Analysis results on subsets with offspring and potential parental genotypes suggested this full-sib relationship was not mistaken for a parent-offspring relationship. Although an increase of the larval stage by six years is highly unlikely, this finding supports a close relationship among sampling years. Collectively this indicates that the two patches for which sampling was sufficient to test temporal patterns, are largely self-sustaining and not prone to extinction-recolonization processes.

| Fine scale structure
Understanding the role of geographical distance is important for defining management units. When we explored the local spatial structure in Overijse, no clear IBD pattern was present, while sPCA results showed a subtle division of northern and southern samples (all years together). This could have resulted from local, nonclinical landscape features, as the alluvial valley of the river IJse is unsuitable for oviposition (Thomaes et al., 2008). Spatial autocorrelation analyses showed only significantly more related individuals within the sampling location (distance class 0-150 m) and again in the distance class of 500 m to 1 km. As eight of the eleven full-sib pairs we identified were collected in the same year and location, this supported the social cohesion present within the smallest distance class. The 500 m to 1 km distance class coincides with the majority of maximum dispersal distances recorded for the species, but usually with higher distances for males than females (Rink & Sinsch, 2007;Thomaes, 2009

| Sex-biased dispersal
According to previous studies, female stag beetles exhibited smaller dispersal distances than males (Rink & Sinsch, 2007;Tini et al., 2018). Rink and Sinsch (2011) (Thomaes et al., 2018). While female stag beetles may stay close to their site of emergence when suitable oviposition sites are abundant and nearby (Tini et al., 2018), they can traverse larger distances when such sites are rather scarce or already occupied by other larvae (Harvey, 2007). Long-distance dispersal across generations was also suggested by Cox et al. (2019) on the basis of the maternally inherited mtDNA. A possible explanation for the discrepancy between the genetic and radio-telemetric outcomes is that the telemetric studies did not compensate for the longer female activity period (Harvey et al., 2011)(adult activity period is 2-3 weeks for males and up to 2 months for females in the Belgian populations (Thomaes et al., 2010)). In telemetric studies the initial period of adult activity after emergence is followed, as battery lifetime is limited to 10-15 days.

| Genetic bottleneck
In this study we employed heterozygosity excess and mode shift tests to identify population declines occurring over recent time scales. Specifically, for heterozygosity excess it is estimated that bottleneck footprints will relate to events occurring with <4N e generations, while for mode shift this time frame is even shorter (Cornuet & Luikart, 1996). Heterozygosity excess tests revealed significant bottleneck signatures for both IAM and TPM for all Watermaal-Bosvoorde tests. Significant heterozygosity results were also reported under the IAM in a number of cases for Overijse (2008,2009,2016), with the 2016 sample also reporting a significant mode shift. In general the biggest criticisms of the bottleneck detection tests employed here, and other such tests, have been their high Type II error with genetic studies often failing to identify bottlenecks despite compelling demographic data (Busch et al., 2007;Le Page et al., 2000;Mardulyn et al., 2008;Peery et al., 2012;Queney et al., 2000;Steinfartz et al., 2007).
Furthermore, it is important to note that the detection of a recent bottleneck may be compromised in situations where ancestral diversity was already low, as suggested for European stag beetle (Cox et al., 2019). Therefore, significant results reported here strongly support the populations have both experienced or are experiencing recent genetic bottlenecks.
Small population size is an important factor in loss of genetic diversity, exacerbated by genetic bottlenecks, which may increase the risk of inbreeding (Allendorf et al., 2013). Here we generated estimates of contemporary effective (or breeding) population sizes using two temporal methods and two single samples estimates (LD and sibship). Although stag beetle was believed to be a species with discrete generations, our results seemed to suggest otherwise. Stag beetles might be considered more as a semelparous species with variable age at maturity, described as the Salmon-model by Waples (2005). In this case, single sample estimates are of N b (effective number of breeders) and are affected by N b in prior years when sampling adults. Estimates of N b derived by the temporal methods, on the other hand, apply to a number of years prior to both sampling years (Waples, 2005).
Estimates of N e or rather of N b were generally small, particularly for Sint-Genesius-Rode and Watermaal-Bosvoorde. In Watermaal-Bosvoorde, N b estimates were similar across methodologies. Monitoring of stag beetle in Watermaal-Bosvoorde started in 1962 (R. Cammaerts, unpublished data) and showed a continuous population present. Initially Watermaal-Bosvoorde held only one stag beetle population, in a forested slope a few 100 meters north of Jagersveld. It was assumed Jagersveld and other nearby sites were colonized by stag beetles from that particular slope. Jagersveld became eligible, man-made habitat for stag beetle when oak logs for fencing were incorporated in a school garden. The clear signs of a reduction in population size could therefore be attributed, at least partly, to a founder event instead of a bottleneck. The lowest N b estimate was found for Sint-Genesius-Rode. Despite the low sample size, which obliged us to use only the single sample estimators, the estimate of N b is in agreement with the lower gene diversity found in the population.
Overijse holds the highest gene diversity and allelic richness of the three populations, though slightly lower levels than those found in surrounding areas of West Europe (Cox et al., 2019).

| Implications for conservation
The data indicate that the three sampled areas are reciprocally isolated in terms of dispersal and gene flow. Furthermore, populations may comprise limited numbers of successfully breeding individuals and exhibit signatures of recent demographic declines. As such, these populations must be considered vulnerable and managed separately. The reported isolation means that rescue events may be limited and/or unpredictable. Small, recently isolated populations are at risk of reduced viability owing to demographic and genetic (inbreeding) effects which can lead to such populations entering an "extinction vortex" (Gilpin & Soulé, 1986). In this case, if one population crashes the available habitat patch may not be readily recolonized, which in turn could have a domino effect on the persistence of other populations within a metapopulation system. Conservation efforts should focus on increasing habitats (dead wood with the necessary fungi) within and between populations. The lack of such habitats within the suburban landscape here may be a driver of the slight female bias in dispersal as they search for oviposition sites.

ACK N OWLED G M ENTS
We Furthermore, we thank Prof. Christian Vandermotten for his input.

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
Microsatellite genotypes were deposited in Dryad: https://doi. org/10.5061/dryad.573n5 tb5h.