Occasional long-distance dispersal may not prevent inbreeding in a threatened butterfly
BMC Ecology and Evolution volume 21, Article number: 224 (2021)
To set up successful conservation measures, detailed knowledge on the dispersal and colonization capacities of the focal species and connectivity between populations is of high relevance. We developed species-specific nuclear microsatellite molecular markers for the grayling (Hipparchia semele), a butterfly endemic to Europe and of growing conservation concern in North-West Europe, and report on its population genetics, in a fragmented, anthropogenic landscape in Belgium. Our study included samples from 23 different locations nested in two regions and additional historical samples from two locations. We assessed contemporary, long-distance dispersal based on genetic assignment tests and investigated the effect of habitat loss and fragmentation on the population genetic structure and genetic variation using data of nine microsatellite loci.
Detected dispersal events covered remarkably long distances, which were up to ten times larger than previously reported colonisation distances, with the longest movement recorded in this study even exceeding 100 km. However, observed frequencies of long-distance dispersal were low. Our results point to the consequences of the strong population decline of the last decades, with evidence of inbreeding for several of the recently sampled populations and low estimates of effective population sizes (Ne) (ranging from 20 to 54 individuals).
Our study shows low frequencies of long-distance dispersal, which is unable to prevent inbreeding in most of the local populations. We discuss the significance for species conservation including future translocation events and discuss appropriate conservation strategies to maintain viable grayling (meta) populations in highly fragmented, anthropogenic landscapes.
In recent years, terrestrial insect declines have been reported worldwide, with habitat loss and fragmentation as main drivers of population decline in anthropogenic landscapes . Habitat loss and fragmentation result in increased distances between local populations and between such populations and vacant habitat patches. In such spatially heterogeneous environments, movement through the matrix typically comes with higher costs than movements within habitat patches . A longer search time for resource patches in the landscape matrix will increase the cost of dispersal (e.g. mortality risks, predation risks, energetic reserve exhaustion, and deferred costs), which may further reduce the physiological conditions of immigrants; all these factors will contribute to reduced fitness of immigrants after they dispersed through the hostile matrix . Habitat fragmentation can thus hinder or even prevent dispersal, and thus gene flow between populations. Furthermore, it can prevent natural colonisation when distances between populations or to potentially suitable habitat patches become greater than the dispersal ability of the species . Subsequently, severely reduced gene flow and the effects of genetic drift can cause a loss of genetic diversity and increase genetic differentiation between historically connected populations, independent of local adaptation. Complete isolation of populations, among other factors such as population size, may eventually lead to an increase in inbreeding and extinction probability .
The grayling (Hipparchia semele) is a butterfly species endemic to Europe that is highly affected by habitat loss and fragmentation, in particular in NW Europe. It is listed as ‘Threatened’ in Flanders (northern Belgium)  and as ‘Vulnerable’ both in the UK  and the Netherlands . In most NW European countries, the grayling mainly occurs in coastal dunes and inland heathlands which are fragmented and restricted to typically small and isolated remnants [9, 10] embedded in an anthropogenic and agricultural landscape matrix. In Belgium, for example, coastal dune areas decreased by 36% during the twentieth century mainly due to building for tourism . Similarly, heathlands strongly decreased in northern Belgium; only 11% of the heathland surface was left in 1995 compared to 1965 and the decline is still ongoing . Heathlands are negatively impacted by eutrophication through aerial nitrogen deposition, and the abandonment of traditional land use, such as sheep grazing and burning [9, 13]. European heathlands and coastal dunes with herbaceous vegetation ('grey dunes') are of high conservation priority and they are protected under the EU Habitats Directive . The grayling can be considered as an umbrella species since it requires rather large patches of open vegetation of different succession stages, which offers habitat opportunities for several other invertebrates of European heathlands and coastal dunes .
For the conservation of species in fragmented biotopes that need additional measures relative to standard vegetation management, species action plans should address the coordination of conservation and management actions on a landscape-scale, encompassing a network of suitable habitats and metapopulations . Conservation efforts should thus not only concentrate on maximizing the habitat quality of currently occupied patches, but also of vacant potentially suitable habitat which could be (re-)colonized  and on the improvement of connectivity between local populations. This helps species to move around the landscape and increase the rate of colonization and gene flow . Clearly defined spatial conservation units are a handy tool for guiding conservation actions on a local scale. For the delimitation of such conservation units, detailed knowledge on the dispersal and colonization capacities of the focal species are required.
In demographic studies, dispersal is often estimated with direct methods such as mark-release-recapture (MRR) studies, the dynamics of patch colonization and patch extinctions or data on range expansions. Other mobility indices of insect species, such as mobility assessment (e.g. very sedentary versus highly mobile) and relative flight speed, are usually based on expert judgements . These methods are often constrained by manpower and sampling area . This makes rare long-distance movements difficult to detect within such studies, which may lead to the underestimation of the dispersal ability of species (e.g. ). Conversely, genetic methods provide indirect estimates of realized dispersal. Assignment tests based on molecular markers, among other indirect genetic methods, provide an alternative and powerful approach for getting more insight in the functional connectivity and gene flow between populations and long-distance movements .
While previous studies on the grayling focussed on habitat use and mobility [15, 23, 24], metapopulation dynamics [15, 25], disturbance tolerance  and local adaptation to climate variables , this study addressed the species’ potential for gene flow and dispersal ability based on neutral molecular markers. In this study, our aim was to (1) estimate contemporary, long-distance dispersal of the grayling in northern Belgium based on genetic assignment tests and (2) investigate how habitat loss and fragmentation affect the population genetic structure and patterns of genetic variation. We hereto developed species-specific polymorphic nuclear microsatellite markers, which will be useful tools for other population genetic studies of the grayling. Finally, we recommend appropriate conservation, translocation and reintroduction strategies to restore viable populations of the grayling in northern Belgium, a region characterized by highly fragmented landscapes under strong anthropogenic pressure.
Within Europe, the grayling Hipparchia semele (Linneaus, 1758) has a widespread distribution, but it is absent in the northernmost parts and the highest European mountains . See Additional file 1: S1 for a map of the distribution of the grayling in North-West Europe. In Belgium, its occurrence is largely restricted to northern Belgium (Flanders), as most of the grayling populations in southern Belgium (Wallonia) have gone extinct . Due to heathland fragmentation and the loss of heathland and grey dune dynamics, the species has suffered a decline in its distribution in northern Belgium of 47% since 1991 , 55% in the Netherlands since 1950  and 62% in the UK since 1970 . The grayling is a long-living butterfly species; adult individuals have a lifespan of up to 40 or even 60 days . The species has one generation per year. The eggs hatch around August, and larvae grow in four instars from August to the following June.
Current populations of the grayling were located by mapping observations made by skilled volunteers using an online data portal (http://waarnemingen.be). For our genetic study, we sampled the locations in northern Belgium where a minimum of ten individuals were observed in one year in the time period 2015–2019. Sampling occurred during the summer of 2020. We assumed that each sampling location contained a single deme; a local population consisting of closely related individuals that form a distinct gene pool, and where mating is random. We sampled 23 different locations nested in two different regions: the Belgian coast in north-west Belgium (five locations) and the inland Campine region located in the north-east of Belgium (18 locations) (Fig. 1, Table 1) (Additional file 1: S2). On two other inland locations that met our criteria (KSV, GRS), were no individuals found. The sampling areas ranged from 4197m2 to 487 435m2 with a mean of 76 865m2. The maximum distance between samples within a location ranged from 89 to 3132 m with a mean of 668 m. The distance to the nearest population ranged from 644 m to 46 730 m with a mean of 6800 m. Although we sampled all known recent populations in Belgium, we cannot exclude that we have missed a few recently colonized small populations within the regions.
Sampling, DNA extraction and amplification
Per sampling location, we took wing-clips (2–3 mm2 tissue) of 30 grayling butterflies, except for the location Maatheide (MAH) and Teutelberg (TEB), where only 2 and 25 individuals were detected, respectively. This resulted in a total of 657 non-destructive wing-clip samples. Wing-clip sampling does not reduce flight ability nor does it affect survival rate (e.g.), although there may be some short-term behavioural changes associated with capture and manipulating individuals . Additionally, we included 24 and 20 grayling individuals collected in 2001 at the locations Mechelse Heide (MEH) and Teutelberg (TEB), respectively. They were stored at −80 °C. These historical samples allowed a comparison of genetic diversity at two time points of the latter populations. The wing-clips were preserved in cryotubes filled with 96% ethanol until analysis.
Before extraction of genomic DNA (Additional file 1: S3), each wing fragment was air-dried for one hour to vaporise all the 96% ethanol. Genomic DNA was extracted from each wing fragment, which was homogenized in a solution of 100 µL 6% Chelex InstaGene Matrix solution (Biorad) and 5 µL proteinase K (> 600 mAU/ml, Qiagen). Microsatellite development and genotyping services were carried out by AllGenetics & Biology SL (www.allgenetics.eu) (Additional file 1: S4). For each sample, 20 polymorphic nuclear microsatellites were amplified. PCR products were run on an ABI 3500 analyser with the GeneScan-600 LIZ size standard (Applied Biosystems) and analysed using Geneious Prime 2019.3.2 (https://www.geneious.com). Details on microsatellites and PCR conditions are given in Additional file 1: S5. Assessment of the genotyping error was performed by including duplicate DNA-extractions for 24 wing-clip samples and, in addition, two to five independent duplicated PCR-amplifications of 21 wing-clip samples.
In total, 701 butterfly wing clips were genotyped. Five microsatellite loci were discarded from the analyses due to stutter peaks and subsequent difficulty in allele scoring. Hence, we worked with 15 microsatellite loci. Next, we examined the assumptions of Hardy–Weinberg equilibrium and of no linkage disequilibrium (LD). We used GENEPOP v4.3  to test for significance of deviations from Hardy–Weinberg equilibrium at individual loci, corrected for multiple testing by the Bonferroni-method . Next we tested for LD for each pair of loci in each population (using the Markov chain method and default parameter settings) and estimated the frequency of null alleles (r) by using the Dempster method  (Additional file 1: S6). A high amount of null alleles is common in butterfly studies (e.g. ), which we also found in this study. To check the influence of loci with deviation from Hard-Weinberg proportions and null alleles on further analyses on the robustness of the results, we followed Waples . We therefore compared results for overall (G’ST, DEST) and population-specific genetic diversity (FIS, Ho, He, Ar), with and without each of these loci. Six loci were discarded during these steps and the two samples of location Maatheide (MAH) were discarded due to missing data at more than 3 loci. The following data analyses were performed based on nine microsatellite loci and 641 unique samples from 24 sampling locations (two and 22 locations sampled in 2001 and 2020, respectively). A total of 63 alleles, with an average of 3.6 alleles per locus, were observed and the dataset had an overall proportion of missing data of 3.7%. For additional information, see Additional file 1: S7.
Genetic assignment methods of dispersal inference make use of individual genotypes and population allele frequencies to estimate where individuals were born or not . GENECLASS2  uses individual assignment tests to determine the putative first-generation dispersers. We used the Frequency based method (Additional file 1: S8) and the Bayesian method as assignment criterion, of which the latter is assumed to be the most efficient method . We did not sample populations in the neighboring countries (the Netherlands, Germany and France) and we cannot exclude that we have missed a few recently colonized small locations in Flanders. We therefore selected the L_home criterion (with L_home being the likelihood of drawing the individual genotype from the population where the individual has been sampled), which is recommended when it is known that not all populations have been sampled . This criterion is however less powerful than the L_home/L_max criterion (the ratio of L_home to the maximum likelihood observed for drawing the individual genotype from any population, including where the individual was sampled) . We ran 1000 Monte Carlo resampling algorithms. Covered distances of long-distance movements were measured as the distance between the location where the individual was assigned to (L_home) and the location where the individual was sampled.
For each sampled location, we calculated the average number of individuals genotyped per locus (N), average number of observed alleles per locus (A), effective number of alleles (Ae), average observed heterozygosity (Ho), Nei’s overall expected heterozygosity (He) and the mean inbreeding coefficient (FIS) using GenAlEx v6.501  and the R package DiveRsity . We additionally calculated FIS corrected for the potential presence of null alleles with the program INEst 2.2 . Next, we used the rarefaction method implemented in HP-Rare v1.0  to calculate the average number of alleles per locus (Ar) and number of private alleles (Ap) adjusted for sample size, based on a minimum sample size of 13 individuals. To detect significant differences in genetic diversity indices, ANOVA analysis were conducted in R . We tested for recent bottlenecks using the program BOTTLENECK v1.2.02  which evaluates deviations of He from the values expected at mutation–drift equilibrium (He > He eq) (Additional file 1: S9). We estimated the contemporary, local effective number of individuals (Ne) for each sampling location by using the linkage disequilibrium method (LDNe) implemented in the program NeEstimator v2  and the sibship assignment method in the program Colony2 [49, 50] (Additional file 1: S10). The correlation of the area of the site and the distance to the nearest population to Ne, FIS, Ar and He were tested with univariate linear regressions in R .
We used the R package DiveRsity to calculate global, regional and pairwise genetic differentiation between populations. We calculated G’ST  and additionally DEST  based on 1000 bootstraps. Both G’ST and DEST are zero (or slightly negative) when there is no differentiation between populations and one at complete differentiation. Differences in the distribution of genotypes between all population pairs were assessed through calculation of 95% confidence intervals (C.I.) of the G’ST values with a bias-corrected bootstrapping method (1000 bootstraps). Using the ISOLDE program in GENEPOP v4.3, we tested for an overall and regional pattern of genetic heterogeneity over geographic distances with the isolation-by-distance-test (IBD) . We used linearized FST-estimates (FST/(1–FST)) for pairs of populations and the Euclidian geographic distance in kilometres between sampling locations. We set the number of permutations for Mantel test to 100 000 and the rest of the parameters on default parameter settings. Next, we used GenAlEx v6.501 to estimate levels of hierarchical structuring within populations, among populations and among regions by analysis of molecular variance (AMOVA) with populations nested within two regions (coastal and inland heathland) (999 permutations). Bayesian analyses of population structure was performed with the programs STRUCTURE v2.3.4  and BAPS v6  (Additional file 1: S11). Lastly, we further investigated the genetic structure among populations using a Principal Component Analysis (PCA) in the program GenAlEx v6.501.
The population assignment test using the Bayesian method identified seven putative first-generation long-distance dispersers (p ≤ 0.01) (Table 2). Putative dispersal events mainly occurred between inland populations (range: 13–69 km, mean: 33 km). Remarkably, one putative disperser originated in a coastal population (ZWI) and moved to an inland population (TLM), covering a distance of 142 km. This putative long-distance dispersal event was also detected by the STRUCTURE analysis (Fig. 2) where the latter putative disperser was assigned to the coastal population cluster. An extra quality and scoring check of the genetic profiles of this putative long-distance disperser confirmed this result. We detected 16 dispersers with a permissive p-value of 0.05 (Table 2). Among these, an even longer dispersal distance was detected. A putative disperser originated in an inland population (BAL) and moved to a coastal population (WEH), covering a distance of 188 km. While the latter method provides the same power as the L_home/L_max criterium with a p-value of 0.01, it also increases the error rate . We further only consider migrants with p ≤ 0.01. Results of the Frequency based method are given in Additional file 1: S8.
We found higher levels of genetic diversity (in terms of allelic richness) in the inland populations than in the coastal populations (Ar: 3.26 and 2.83 resp., ANOVA test p < 0.05). Genetic diversity ranged from 2.8 to 4.3 alleles/locus and values of expected heterozygosity (He) ranged from 0.35 to 0.46. Corrected allelic richness values (Ar) ranged from 2.6 to 3.7 and the corrected number of private alleles (Ap) ranged from 0.0 to 0.1 (Table 1). Genetic diversity indices are shown in Additional file 1: S12. Evidence of inbreeding calculated with the R package DiveRsity was found in four of the five (80%) coastal populations (mean FIS: 0.161, range: 0.0859–0.277) and in twelve of the seventeen (70%) recent inland populations (mean FIS: 0.139, range: 0–0.259) (Table 1). The historical populations sampled in 2001 showed no evidence of inbreeding. After correcting FIS-values for the potential presence of null alleles, we still found evidence of inbreeding in two of the five (40%) coastal populations (mean corrected FIS: 0.086, range: 0.039–0.131) and in five of the seventeen (30%) recent inland populations (mean corrected FIS: 0.067, range: 0.016–0.137) (Additional file 1: S13). Only three coastal populations (SLD, TYZ and ZWI) showed weak evidence of genetic bottlenecks (under IAM p-values < 0.05, ZWI under TPM p-value < 0.05) (Additional file 1: S9). Estimation of Ne based on the sibship assignment method produced reliable results, as the two replicate runs yielded very similar values. Ne was estimated between 20 (ZWI) and 54 (ZBN) individuals using the Random Mating model (Table 1). We found a significant positive correlation between Ne and the area of a sample site (r = 0.64, p = 0.002). We further detected a positive trend between Ar and the area of a sample site (r = 0.40, p = 0.068) and a negative trend between Ar and the distance to the nearest population (r = − 0.40, p = 0.068). The other tested correlations were not significant.
The global genetic population differentiation was low (global G’ST = 0.059 (95% C.I.: 0.046–0.073) and DEST = 0.014 (95% C.I.: 0.009–0.020)). Regionally, genetic differentiation was even lower between the inland populations (regional G’ST = 0.039 (95% C.I.: 0.025–0.054) and DEST = 0.007 (95% C.I.: 0.002–0.013)) and slightly higher between the coastal populations (regional G’ST = 0.060 (95% C.I.: 0.032–0.090) and DEST = 0.017 (95% C.I.: 0.008–0.029)). Both global and regional differentiation were significant (G’ST: p < 0.001, global DEST: p < 0.001, regional DEST: p < 0.05). Overall, pairwise population G’ST-values ranged from − 0.012 to 0.160 (Additional file 1: S14). Regionally, pairwise population G’ST-values between inland populations ranged from − 0.012 to 0.114 and between coastal populations from 0.000 to 0.128. There was significant genetic differentiation (p < 0.05) between 97 of the 276 population pairs. IBD-analysis showed that genetic variability was only globally structured as indicated by the significant linear relationship (p = 0.00021, R2 = 0.28). Within the inland region, we found no significant isolation by distance (p = 0.185, R2 = 0.027). The IBD-analyses of the coastal region showed a positive trend (p = 0.087, R2 = 0.55). So there is weak evidence for isolation-by-distance in this region (Additional file 1: S15). The AMOVA indicated that the genetic variance was the highest within populations (95%), but only explained by region and population by 3% and 2% respectively (Additional file 1: S16). The Bayesian cluster analysis performed with BAPS (Additional file 1: S11) and STRUCTURE showed the highest probability for three clusters (K = 3) (Fig. 2). Individuals of the coastal populations were clearly clustered together. Since the inland populations showed more genetic admixture, STRUCTURE analyses showed no clear differentiation between the two inland clusters. Remarkably, the inland MEH population sampled in 2001 was assigned in the same cluster as the coastal populations, while the samples collected in MEH in 2020 clustered with neighbouring inland locations. Similarly, in the PCA populations were clustered according to region, separated by the first axis, with exception of the samples of the inland population MEH collected in 2001; the latter clustered together with the coastal populations (Additional file 1: S17).
We estimated long-distance dispersal based on microsatellite markers of the grayling Hipparchia semele inhabiting highly fragmented heathland or dune vegetation patches imbedded in an anthropogenic landscape in Belgium. Remarkably, maximum dispersal distance exceeded 100 km. Nevertheless, observed frequencies of long-distance dispersal events were low. The low genetic diversity reflected the recent population decline of the species in northern Belgium. While we found no strong evidence of severe bottlenecks, we detected inbreeding in several populations and we found low estimates for effective population sizes. Our results suggest that former studies likely have underestimated the dispersal distances that the grayling can cover. However, occasional long-distance dispersal events may be insufficient to preserve populations in this highly fragmented landscape in the long-term.
High dispersal ability
The grayling is thought to be a fairly mobile butterfly species. Previous mark-release-recapture (MRR) studies observed covered distances of 1.2 km , 1.5 km  and 4 km  and colonisation distances up to 10–15 km . We found seven putative first-generation dispersers (7 out of 599; 1.2%), of which six migrated between the inland heathland populations and one between a coastal and an inland population. The dispersal distances between the inland populations ranged from 13 to 69 km with a mean of 33 km, while the distance covered between the coastal population and the inland population was 142 km.
The dispersal distances observed in this study are substantially longer than the ones reported in previous MRR studies. In this study, we covered the entire distribution area of the grayling in northern Belgium (covering 12 625 km2). In contrast, previous MRR studies in the region [15, 31, 57] focussed on dispersal at a smaller spatial scale within one Belgian nature reserve. While the results of MRR studies can give more insight in local dispersal patterns (meters to several kilometres), the results based on genetic markers provide more insight in the dispersal ability of species at the regional spatial (landscape) scale (up to several hundreds of kilometres). Other studies found similar differences between dispersal distances obtained by observational methods such as MRR and those obtained from genetic marker studies (e.g. ). For example, for the sedentary butterfly Phengaris (Maculinea) alcon, Vanden Broeck et al.  found evidence of a dispersal distance of 2.9 km. In the same region, however, the maximum recorded movement of P. alcon from MRR studies was 0.5 km and the maximum distance recorded from spontaneous colonization data so far was 1.7 km . To get a complete overview of the dispersal behaviour of a species, both direct (Mark-Release-Recapture) and indirect methods (genetic analyses) should ideally be implemented.
Rare long-distance dispersal events
Although we provide evidence for dispersal ability of the grayling over unexpected long distances, the observed frequencies of long-distance dispersal events were low (1.2%). This observed low frequency of long-distance dispersal can, at least partly, be explained by the low genetic differentiation between populations, which limits the power of the assignment tests to detect dispersal events . Furthermore, since we possibly did not sample all populations from which migrants may have originated (nearby populations in neighbouring countries (Fig. 1) and recently colonised populations within the region), we tested for the likelihood that an individual grayling was foreign-born, i.e. born in another population than where it was sampled . The latter also results in a lower power to detect migrants compared to a sample including all potential source populations. Therefore, the real number of dispersers is likely to be underestimated in this study.
However, the majority of the detected long-distance dispersal occurred between inland populations, while no dispersal between the coastal populations was detected. Even with less strict detection methods, the vast majority of detected dispersal events occurred between the inland populations, with sparse dispersal events within the coastal region and between regions (Table 2; Additional file 1: S8). In contrast to the inland region, we found weak evidence of isolation-by-distance in the coastal region, supporting the observed lack of dispersal events between sampled locations within the coastal region. Increasing long-distance dispersal events would decrease the isolation and would likely lower the strength of inbreeding among the coastal populations. For the inland region, the sampled locations consist of a metapopulation with interaction by gene flow between the subpopulations. Here, sampled locations are less isolated compared to the sampled locations of the coastal region, despite longer distances between inland populations. This indicates a higher connectivity and lower extinction risk for the grayling in the inland region compared to the coastal region. The differences in dispersal events observed between the coastal and inland region could be explained by differences in the landscape matrix (for detailed maps, see Additional file 1: S18). Low long-distance dispersal events in the Belgian coastal region can likely be attributed to the high levels of urbanization. Buildings, parking lots and roads separate small patches of grey dune habitat by many kilometres  and may act as barriers limiting gene flow. A study on the effect of urbanization on the species richness of different taxa found that butterflies, and the most mobile and specialist species in particular, were strongly negatively affected by urbanization . However, our findings are in contrast with the results of a population genetic study of the endemic butterfly Atrytonopsis sp. on the highly urbanized coast of North Carolina (USA). Based on Amplified Fragment Length Polymorphic (AFLP) molecular markers, Leidner and Haddad  found that not coastal urbanization but natural barriers (ocean and forests) limited butterfly dispersal between coastal populations. It must be noticed that, in the latter study, conclusions on actual migration rates were based on the population genetic structure obtained by the program STRUCTURE  and on Isolation-by-distance (IBD) analyses. These methods are however considered unsuitable for making reliable deductions about actual dispersal events. Generally, STRUCTURE yields poor individual assignments to source populations and it does not reliably recover the actual population structure when sampling is unbalanced . Furthermore, patterns of IBD can arise from limited dispersal but also from historical demographical processes such as founder effects and historical gene flow rather than contemporary dispersal rates [52, 63]. Due to these concerns, we used genetic assignment tests instead. Genetic assignment tests are based on the distribution of observed allele frequencies to draw inference about where individuals were or were not born, allowing direct, real-time estimates of dispersal . These assignment tests are useful for the identification of immigrant individuals in the current generation and thus the potential for contemporary gene flow, which was the interest of this study.
In contrast to the coastal region, the grayling inland populations are principally surrounded by a matrix of woodland, heathland and meadows. Woodlands are often considered to be dispersal barriers for butterflies of open habitats (e.g. heathlands, grasslands), but Nowicki et al.  showed that long-distance butterfly dispersal is not always prevented by forests. Moreover, towards the southern and eastern part of the grayling’s distribution range, woodland and trees are a part of its habitat . The inland populations may therefore be located in a more favourable matrix compared to the coastal populations. We can however only speculate about possible barriers in the different regions. A study on landscape genetic analysis, which takes spatial structures into account, should provide more robust evidence for potential barriers for the grayling.
Fragmentation effects on genetic diversity
When habitat fragmentation hampers gene flow, populations can suffer from loss of genetic diversity reflected by low effective population sizes, resulting in inbreeding and extinction risks. We did find evidence of inbreeding in several populations and low estimates for effective population sizes (range: 20–54 individuals). Inbreeding increases levels of homozygosity and exposes deleterious recessive alleles, which can result in inbreeding depression and reduced population fitness [5, 66, 67]. These negative consequences on the population level can have implications for metapopulation dynamics . As a loss of populations reduces the total number of migrants, increased population extinctions caused by inbreeding would decrease the colonization probability of unoccupied habitat patches . An entire metapopulation, and especially a small metapopulation, might thus suffer from an increased extinction risk caused by inbreeding effects [68, 69].
For further investigating the potential effects of fragmentation on genetic diversity, we included samples collected two decades ago (dated 2001) from two inland populations (MEH and TEB) and compared the genetic diversity with recent samples of the corresponding locations (dated 2020). Levels of genetic diversity, in terms of effective population size and allelic richness, were similar between the recent and historical samples for both locations. However, for the location MEH, there was no evidence of inbreeding in the historical sample, while we found weak evidence of inbreeding in the current population. This may indicate that the levels of diversity did not change over the past two decades, but that a consistent low average population size (we estimated low effective population sizes and the census population sizes have been decreasing over decades ) caused an increase in inbreeding over time. Additionally, a remarkable result is the clustering of the 2001 population of MEH with the coastal populations. This might indicate that the two regions were more genetically similar in the past likely because of a higher connectivity between the regions.
Previous studies on the two heathland butterfly species Phengaris alcon  and Plebejus argus (unpublished data) occured in the same inland Campine region of this study. The results of these studies and our study on H. semele were obtained using the same methods. P. alcon is listed as ‘Critically endangered’, while P. argus is listed as ‘Endangered’ in Flanders . In terms of effective population size and genetic diversity, we found similar results to both previously studied species. The inbreeding coefficients in our study were similar to the ones found in the P. argus populations. The study on P. argus also found a low number of dispersers within the inland region and no significant Isolation-by-distance.
Our inbreeding values should however be interpreted with caution. While we took into account the presence of null alleles during the selection of the loci for analyses, there is still evidence for null alleles in a small percentage of loci x population combinations. The presence of null alleles causes an increase in homozygosity, which means that levels of inbreeding could be inflated (e.g. ). Additionally, substructure within populations can also lead to a deficiency of observed heterozygotes (the Wahlund effect ), and an increased coefficient of inbreeding [38, 67]. Though we cannot exclude a reduction in heterozygosity caused by subpopulation structure in some of the sampled locations, individuals were generally caught in fairly small sampling areas (Additional file 1: S2). H. semele is also a mobile species that can easily move distances of 100 m and more within a couple of days , covering a large area within a sampling location.
Implications for conservation
Our study corroborates the results of demographic studies that show strong population declines of the grayling principally caused by habitat loss and fragmentation, particularly in North-West Europe (e.g. [6, 7, 25, 30, 71]). Climate change may also have a negative impact on reproductive success , as extreme droughts have already occurred in Belgium during the last few years (e.g. 2018 and 2019). Conservation actions are needed to prevent a further decline and local population extinctions of the grayling in northern Belgium.
Based on the habitat characteristics and the genetic differentiation and mobility, we can assign grayling populations in northern Belgium to two different Functional Conservation Units (FCUs) (cf. ); the coastal conservation unit (FCU 1) and the inland heathland conservation unit (FCU 2). Although these conservation units, separated by at least 70 km, are potentially connected by gene flow and thus not completely isolated from each other, we consider them as different spatial entities in which specific management and restoration measures should be implemented. Estimates of the effective population sizes were higher in larger sample locations. Additionally, the allelic richness of the populations showed a trend of being higher in larger areas and in populations with shorter nearest neighbour distances. Within each conservation unit, management actions should therefore focus on enlarging actual and restore potential habitat areas.
For the coastal conservation unit (FCU 1), as current inter-population connection is likely extremely low, we suggest that the emphasis should be on creating new habitat between existing populations in order to improve population network connectedness. Furthermore, translocation actions such as reinforcements should be considered, as spontaneous exchange of individuals between populations is likely hindered by anthropogenic barriers to gene flow along the Belgian coast and because of the observed inbreeding. The inland conservation unit (FCU 2) consists of a larger metapopulation, in terms of the number of demes and occupied territory, and connectivity between sub-populations is likely higher and less hindered by anthropogenic barriers compared to the coastal populations. However, also for FCU 2, restored and/or new heathland habitat have likely a low chance to be colonized spontaneously within a short time frame. Over the past years, there has been a strong decline of the grayling in Northern Belgium , which is most noticeable in the central and western part of the inland region. Relying on spontaneous colonisation and gene flow for the recovery of populations would therefore be too a considerable risk. Therefore, we recognize targeted re-introductions and translocations as valuable management actions to improve the connection between the populations in the western and eastern part of the inland region.
Translocations should be conducted according to IUCN guidelines . To avoid a high risk of outbreeding depression , founder populations should consist of individuals from the same conservation unit. Founder individuals should be chosen from non-inbred populations with high genetic diversity, as inbreeding and reductions in genetic diversity are likely to lower fitness and thus the persistence of the founder population . For the coastal conservation unit (FCU 1), Schipgatduinen (SGD) could act as a source sub-population for translocations as it is one of the largest suitable areas (the risk of damaging the local populations is low) and it is a population with a low inbreeding value. However, levels of genetic diversity in SGD are low and therefore combining different sub-populations as a source could be considered as an alternative. For the inland conservation unit (FCU 2), the sub-population Zwarte Beek Noord (ZBN) appears the best source population for translocations. It has the highest estimated level of genetic diversity and the highest estimated effective population size of all the sub-populations. To assess translocation success, an appropriate monitoring scheme needs to be implemented.
This genetic study on the grayling showed that dispersal distances covered by grayling butterflies were much larger than previously assumed. Long-distance movements are however scarce. Estimated effective population sizes were low and several populations showed evidence for inbreeding, which reflect the consequences of strong population declines during the last decades. Although genetic differentiation was found to be low, both within as between regions, there was a significant differentiation between the coastal and inland region. Therefore, we suggest to treat the two regions as as two separate conservation units. Urgent measures (nature management, translocations, reinforcement or supplementation) are needed for the sustainable conservation of the species in northern Belgium.
Availability of data and materials
The dataset supporting the conclusions of this article is available in the Dryad repository, https://0-doi-org.brum.beds.ac.uk/10.5061/dryad.5tb2rbp4h.
Analysis of molecular variance
Functional conservation unit
- MRR studies:
Principal components analysis
Cardoso P, Barton PS, Birkhofer K, Chichorro F, Deacon C, Fartmann T, et al. Scientists’ warning to humanity on insect extinctions. Biol Conserv. 2020;242:108426.
Dennis RLH. Butterfly habitats, broad-scale biotope affiliations, and structural exploitation of vegetation at finer scales: the matrix revisited. Ecol Entomol. 2004;29:744–52.
Baguette M, Van Dyck H. Landscape connectivity and animal behavior: functional grain as a key determinant for dispersal. Landscape Ecol. 2007;22(8):1117–29.
Hanski I. Habitat connectivity, habitat continuity, and metapopulations in dynamic landscapes. Oikos. 1999;87:209–19.
Frankham R, Ballou JD, Briscoe DA. Introduction to conservation genetics. Cambridge: Cambridge University Press; 2010.
Maes D, Vanreusel W, Jacobs I, Berwaerts K, Van Dyck H. Applying IUCN Red List criteria at a small regional level: a test case with butterflies in Flanders (north Belgium). Biol Cons. 2012;145(1):258–66.
Fox R, Warren MS, Brereton TM, Roy DB, Robinson A. A new red list of british butterflies. Insect Conserv Divers. 2011;4(3):159–72.
van Swaay CAM. Basisrapport Rode Lijst Dagvlinders 2019 volgens Nederlandse en IUCN-criteria. Wageningen: De Vlinderstichting; 2019.
Webb NR. The traditional managment of European heathlands. J Appl Ecol. 1998;35:987–90.
Exeler N, Kratochwil A, Hochkirch A. Restoration of riverine inland sand dunes: implications for the conservation of wild bees (Apoidea). J Appl Ecol. 2009;46:1097–105.
Provoost S, Bonte D. Levende duinen: een overzicht van de biodiversiteit aan de Vlaamse kust. Brussel: Instituut voor Natuurbehoud; 2004.
Van Landuyt W, Vanhecke L, Hoste I, Hendrickx F, Bauwens D. Changes in the distribution area of vascular plants in Flanders (northern Belgium): eutrophication as a major driving force. Biodivers Conserv. 2008;17(12):3045–60.
Rose RJ, Webb NR, Clarke RT, Traynor CH. Changes on the heathlands on Dorset, England, between 1987 and 1996. Biol Cons. 2000;93:117–25.
European Commission. Interpretation manual of European Union Habitats-EUR27. Brussels: European Commission DG Environment; 2007.
Maes D, Ghesquiere A, Logie M, Bonte D. Habitat use and mobility of two threatened coastal dune insects: implications for conservation. J Insect Conserv. 2006;10(2):105–15.
Bourn NAD, Bulman CR. Landscape scale conservation, theory into practice. In: Kuhn E, Feldmann R, Thomas JA, Settele J, editors. Studies on the ecology and conservation of butterflies in Europe general concepts and case studies 1. Dorset: Butterfly Conservation; 2005. p. 111–2.
Maes D, Vanreusel W, Talloen W, Dyck HV. Functional conservation units for the endangered Alcon Blue butterfly Maculinea alcon in Belgium (Lepidoptera: Lycaenidae). Biol Cons. 2004;120(2):229–41.
Ellis S, Wainwright D, Berney F, Bulman C, Bourn N. Landscape-scale conservation in practice: lessons from northern England. UK J Insect Conserv. 2010;15(1–2):69–81.
Stevens VM, Turlure C, Baguette M. A meta-analysis of dispersal in butterflies. Biol Rev Camb Philos Soc. 2010;85(3):625–42.
Schneider C. The influence of spatial scale on quantifying insect dispersal: an analysis of butterfly data. Ecological Entomology. 2003;28:252–6.
Vanden Broeck A, Maes D, Kelager A, Wynhoff I, WallisDeVries MF, Nash DR, et al. Gene flow and effective population sizes of the butterfly Maculinea alcon in a highly fragmented, anthropogenic landscape. Biol Cons. 2017;209:89–97.
Kim KS, Sappington TW. Population genetics strategies to characterize long-distance dispersal of insects. Journal of Asia-Pacific Entomology. 2013;16(1):87–97.
Schirmel J, Fartmann T. Coastal heathland succession influences butterfly community composition and threatens endangered butterfly species. J Insect Conserv. 2014;18(1):111–20.
Tropek R, Cizek O, Kadlec T, Klecka J. Habitat use of Hipparchia Semele (Lepidoptera) in its artificial stronghold: necessity of the resource-based habitat view in restoration of disturbed sites. Pol J Ecol. 2017;65(3):385–99.
van Strien A, van Swaay C, Kéry M. Metapopulation dynamics in the butterfly Hipparchia semele changed decades before occupancy declined in The Netherlands. Ecol Appl. 2011;21(7):2510–20.
Bonte D, Maes D. Trampling affects the distribution of specialised coastal dune arthropods. Basic Appl Ecol. 2008;9:726–34.
Middlebrook I, Hardy PB, Botham MS, Dennis RLH. The importance of unique populations for conservation: the case of the great orme’s head grayling butterfly Hipparchia semele (Linnaeus, 1758) (Lepidoptera: Satyrinae). J Insect Conserv. 2019;23(2):381–91.
Kudrna O, Harpke A, Lux K, Pennerstorfer J, Schweiger O, Settele J, et al. Distribution atlas of butterlfies in Europe. Halle: Gfs; 2011.
Fichefet V, Barbier Y, Baugnée J-Y, Dufrêne M, Goffart P, Maes D, et al. Papillons de jour de Wallonie (1985–2007). Region Wallonne: Gembloux; 2008. p. 320.
Maes D, Fajgenblat M, Herremans M, Piesschaert F, Vantieghem P, Jacobs I, et al. IUCN Rode Lijst van de dagvlinders in Vlaanderen 2021. Rapporten van het Instituut voor Natuur- en Bosonderzoek Instituut voor Natuur- en Bosonderzoek, Brussel. 2021.
Vanreusel W, Cortens J, Van Dyck H. Herstel van dagvlinderpopulaties in en om het Nationaal Park Hoge Kempen. Universiteit Antwerpen (UIA-UA)—in opdracht van afdeling Natuur van het Ministerie van de Vlaamse Gemeenschap, Wilrijk. 2002.
Kingsolver JG. Experimental analyses of wing size, flight, and survival in the western white butterfly. Evolution. 1999;53(5):1479–90.
Hamm CA, Aggarwal D, Landis DA. Evaluating the impact of non-lethal DNA sampling on two butterflies, Vanessa cardui and Satyrodes eurydice. J Insect Conserv. 2009;14(1):11–8.
Rousset F. genepop’007: a complete re-implementation of the genepop software for Windows and Linux. Mol Ecol Resour. 2008;8(1):103–6.
Hochberg Y. A sharper Bonferroni procedure for multiple tests of significance. Biometrika. 1988;75(4):800–2.
Dempster AP, Laird NM, Rubin DB. Maximum likelihood from incomplete data via the EM algorithm. J Roy Stat Soc. 1977;39(1):1–38.
Keyghobadi N, Roland J, Strobeck C. Genetic differentiation and gene flow among populations of the alpine butterfly, Parnassius smintheus, vary with landscape connectivity. Mol Ecol. 2005;14(7):1897–909.
Waples RS. Testing for Hardy-Weinberg proportions: have we lost the plot? J Hered. 2015;106(1):1–19.
Paetkau D, Slade R, Burden M, Estoup A. Genetic assignment methods for the direct, real-time estimation of migration rate: a simulation-based exploration of accuracy and power. Mol Ecol. 2004;13:55–65.
Piry S, Alapetite A, Cornuet JM, Paetkau D, Baudouin L, Estoup A. GENECLASS2: a software for genetic assignment and first-generation migrant detection. J Hered. 2004;95(6):536–9.
Cornuet JM, Piry S, Luikart G, Estoup A, Solignac M. New methods employing multilocus genotypes to select or exclude populations as origins of individuals. Genetics. 1999;153:1989–2000.
Peakall R, Smouse PE. GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research—an update. Bioinformatics. 2012;28(19):2537–9.
Keenan K, McGinnity P, Cross TF, Crozier WW, Prodöhl PA, O’Hara RB. diveRsity: AnRpackage for the estimation and exploration of population genetics parameters and their associated errors. Methods Ecol Evol. 2013;4(8):782–8.
Chybicki IJ, Burczyk J. Simultaneous estimation of null alleles and inbreeding coefficients. J Hered. 2009;100(1):106–13.
Kalinowski ST. Do polymorphic loci require large sample sizes to estimate genetic distances? Heredity (Edinb). 2005;94(1):33–6.
R Core Team. R: A language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2020.
Cornuet JM, Luikart G. Description and power analysis of two tests for detecting recent population bottlencks from allele frequency data. Genetics. 1996;144:2001–14.
Do C, Waples RS, Peel D, Macbeth GM, Tillett BJ, Ovenden JR. NeEstimator v2: re-implementation of software for the estimation of contemporary effective population size (Ne ) from genetic data. Mol Ecol Resour. 2014;14(1):209–14.
Jones OR, Wang J. COLONY: a program for parentage and sibship inference from multilocus genotype data. Mol Ecol Resour. 2010;10(3):551–5.
Wang J. A new method for estimating effective population sizes from a single sample of multilocus genotypes. Mol Ecol. 2009;18(10):2148–64.
Hedrick PW. A standardized genetic differentiation measure. Evolution. 2005;59(8):1633–8.
Jost L. G(ST) and its relatives do not measure differentiation. Mol Ecol. 2008;17(18):4015–26.
Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155:945–59.
Corander J, Marttinen P. Bayesian identification of admixture events using multilocus molecular markers. Mol Ecol. 2006;15(10):2833–43.
Tinbergen N. The courtship of the grayling Eumenis (=Satyrus) semele (L.) (1942). London: George Allen and Unwin Ltd; 1972.
Dennis RLH, Sparks TH, Shreeve TG. Geographical factors influencing the probability of Hipparchia semele (L.) (Lepidoptera: Satyrinae) occurring on British and Irish off-shore islands. Glob Ecol and Biogeogr Lett. 1998;7:205–14.
Segers N. Mobility and habitat use of the butterfly Hipparchia semele (Lepidoptera, Satyrinae) in the National Park Hoge Kempen (Belgium) [Master Thesis]. Antwerpen: Universiteit Antwerpen; 2012.
Harper GL, Maclean N, Goulson D. Microsatellite markers to assess the influence of population size, isolation and demographic change on the genetic structure of the UK butterfly Polyommatus bellargus. Mol Ecol. 2003;12(12):3349–57.
Provoost S, Ampe C, Bonte D, Cosyns E, Hoffmann M. Ecology, management and monitoring of dune grasslands in Flanders, Belgium. Littoral 2002 The Changing Coast. 2002:11–22.
Concepción ED, Moretti M, Altermatt F, Nobis MP, Obrist MK. Impacts of urbanisation on biodiversity: the role of species mobility, degree of specialisation and spatial scale. Oikos. 2015;124(12):1571–82.
Leidner AK, Haddad NM. Natural, not urban, barriers define population structure for a coastal endemic butterfly. Conserv Genet. 2010;11(6):2311–20.
Puechmaille SJ. The program structure does not reliably recover the correct population structure when sampling is uneven: subsampling and new estimators alleviate the problem. Mol Ecol Resour. 2016;16(3):608–27.
Meirmans PG. The trouble with isolation by distance. Mol Ecol. 2012;21:2839–46.
Nowicki P, Vrabec V, Binzenhöfer B, Feil J, Zakšek B, Hovestadt T, et al. Butterfly dispersal in inhospitable matrix: rare, risky, but long-distance. Landscape Ecol. 2014;29(3):401–12.
Maes D, Vanreusel W, Van Dyck H. Dagvlinders in Vlaanderen: nieuwe kennis voor betere actie! Tielt: Lannoo nv; 2013.
Drees C, De Vries H, Härdtle W, Matern A, Persigehl M, Assmann T. Genetic erosion in a stenotopic heathland ground beetle (Coleoptera: Carabidae): a matter of habitat size? Conserv Genet. 2009;12(1):105–17.
Zilko JP, Harley D, Hansen B, Pavlova A, Sunnucks P. Accounting for cryptic population substructure enhances detection of inbreeding depression with genomic inbreeding coefficients: an example from a critically endangered marsupial. Mol Ecol. 2020;29(16):2978–93.
Nieminen M, Singer MC, Fortelius W, Schöps K, Hanski I. Experimental confirmation that inbreeding depression increases extinction risk in butterfly populations. Am Nat. 2001;157:237–44.
Nonaka E, Siren J, Somervuo P, Ruokolainen L, Ovaskainen O, Hanski I. Scaling up the effects of inbreeding depression from individuals to metapopulations. J Anim Ecol. 2019;88(8):1202–14.
Wahlund S. Zusammensetzung von population und korrelationserscheinung vom stand-punkt der vererbungslehre aus betrachtet. Hereditas. 1928;11:65–106.
Bos F. Heivlinder Hipparchia semele. Natuur van Nederland. 2006;7:235–7.
Salgado AL, DiLeo MF, Saastamoinen M, Rasmann S. Narrow oviposition preference of an insect herbivore risks survival under conditions of severe drought. Funct Ecol. 2020;34(7):1358–69.
IUCN/SSC. IUCN guidelines for reintroductions and other conservation translocations, version 1.0. IUCN Species Survival Commission, Gland. 2013.
Lynch M. The genetic interpretation of inbreeding depression and outbreeding depression. Evolution. 1991;45(3):622–9.
The authors thank Nico De Regge, Irina De Landtsheer, Koen Devos, Floris Van Den Broeck, Axel Neukermans, Nicolas Vanermen, Dries Bonte, Joost Sturtewagen, Johan Debuck, Arjen Breevaert, Luc De Bruyn, Ann Milbau, Ive Van Krunkelsven, Stern Van Krunkelsven, Marcel Van Waerebeke, Jos Gorissen, , Michel Huysmans, Oscar Maes for sampling assistance and Johan Lamaire, Koen Marechal, Guy Vileyn, Dirk Raes, Jef De Winter, Harry Thys, Willy Pardon, Ghis Palmans, Ernesto Zvar, Michel Broeckmans, Corina Cools, Manu Vermeulen, Jan Appermont and Patrick Reynders for permission during fieldwork. We would also like to thank Chris van Swaay from the Dutch Butterfly Conservation for providing the distribution data of the grayling in North-West Europe. We further really appreciate the constructive comments of the reviewers of BMC Ecology & Evolution. Lastly, we would like to thank the Agency for Nature and Forests (ANB) for financing the study.
This study was funded by the Agency of Nature and Forest (ANB) of the Flemish Government. None of the funding sources were involved in the study design, data collection, analysis, writing or submission of this paper.
Ethics approval and consent to participate
For the execution of the fieldwork, we obtained a permission for deviation of The Flemish Government Decree of 15 May 2009 on species protection and management (‘het Soortenbesluit’). Permissions for entering sampling locations were obtained through local nature managers.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
S1. Distribution of Hipparchia semele in North-West Europe. S2. Samples and sampling sites of Hipparchia semele. S3. DNA extraction. S4. Microsatellite development. S5. Microsatellite loci and PCR-conditions. S6. Estimation of the frequency of null alleles with the Dempster method. S7. Hardy-Weinberg equilibrium and Linkage Disequilibrium testing. S8. Population assignment test: Frequency based method. S9. Genetic bottlenecks. S10. Estimation of effective population size (Ne). S11. STRUCTURE and BAPS analysis. S12. Genetic diversity statistics. S13. The effect of null alleles on the FIS-values. S14. Pairwise G’ST values. S15. Isolation-by-distance analyses. S16. Levels of hierarchical structuring within populations, among populations and among regions estimated by analyses of molecular variance (AMOVA). S17. Plots of the Principal Component Analysis (PCoA). S18. Detailed maps of the land-use in northern Belgium.
About this article
Cite this article
De Ro, A., Vanden Broeck, A., Verschaeve, L. et al. Occasional long-distance dispersal may not prevent inbreeding in a threatened butterfly. BMC Ecol Evo 21, 224 (2021). https://0-doi-org.brum.beds.ac.uk/10.1186/s12862-021-01953-z