Located in the south-western part of Slovakia, the study area lies in the Danubian Plain (Podunajská rovina), part of the Danubian Lowland, belonging to the Pannonian biogeographical region. The landscape was formed by the tributaries of the Danube (Little Danube and others) and by the southern Váh, Nitra and Žitava Rivers. Soil properties and climatic conditions make this region ideal for agriculture. In the area, wheat, sugar beets, sweet corn, vegetables, fruits and tobacco are grown. In the late 19^th century, the study area was located in the centre of wet meadows interlaced by meanders of Žitava River, as shown by the Third Military Survey of the Habsburg Empire (1869–1887) (2018). Over the past 130 years, extensive flood protection measures and construction of the channel network fragmented the area of natural meanders of the Žitava River, which changed the direction and hydrodynamics of their flow and led to the aridification of the surrounding ecosystems. By 1926, the original wet meadows had dried up and a network of meanders had formed into a compact wetland habitat crossed by two channels. The consequences of these interventions were still visible in 1964, but the size of wetland habitats had become further reduced by 1990 (Topographic Maps of Czechoslovakia 2018). Today, remnants of the original meander system exist only as a few patches, where the Pannonian root vole still survives. The extension, shape and connectivity of recent wetland habitats scattered across the agricultural landscape are determined by water levels and temporary floods. In the summer and winter of 2010, the extreme precipitation raised water levels and caused the flooding of the Danube (Slovak Hydrometeorological Institute 2011). The last highest water level rise of the river was recorded in spring 2013 (Jakubcová et al. 2014; Matoková and Smrtík 2014), which also influenced our study area. The localities included in this study are situated in the area interlaced by three channels between Patince, Chotín and Marcelová Villages (Fig. 1).

Figure 1.

The distribution of the sampled localities of root vole subpopulations in the Danubian Lowland. The sampling sites found in the habitat fragments are indicated by letter and number codes. The insert shows the location of the study area in Slovakia.

The research was carried out in different types of habitat fragments like waterlogged areas overgrown by Carex spp. and Phragmites spp. at the edges of channels, oxbow lakes and remnants of former tributaries intersecting large areas of agricultural lands. Animals were captured with live traps baited with apple and cereals. One line transect of 25 traps with 10 m intervals was established for five consecutive nights in each trapping site. Trapping was conducted three times a year (spring, summer, autumn), from 2014 to 2017. Traps were inspected once a day. Upon capture, each vole was investigated for body weight, age and reproductive status before release. At the first capture, the tail tip of each vole was clipped for genetic analyses. The clipped tail tips were put immediately in 96% ethanol and preserved at -85 °C in the laboratory until DNA extraction. Trapping and sampling methods were realised in agreement with the rules of State Nature Protection of Slovak Republic “Species and habitats monitoring of European importance within the Habitats Directive and the Birds Directive” project. Due to the rarity and endangered status of the Pannonian root vole, sample sizes were low in some trapping transects and, therefore, a total of 69 tissue samples from nine sites were used in molecular analyses.

Five characteristics related to wetland habitat fragments were measured or calculated in six time periods between 2004 and 2019 (January 2004, March 2011, April 2014, March 2017, August 2017 and March 2019): 1) number of all habitat fragments presumably suitable for the Pannonian root vole; 2) number of temporarily suitable habitat fragments appearing occasionally, depending on water levels; 3) total area of all habitat fragments measured in hectares; 4) average size of habitat fragments (ha); and 5) overall connectivity of our region of interest. Suitable habitat fragments were delimited based on the subspecies’ known habitat preference for humid, densely vegetated areas (see Introduction), which markedly differed from the vegetation of agricultural parcels in the study area. The size of particular habitat fragments (ha) was calculated in QGIS software 3.4.12-Madeira (QGIS Development Team 2019). To calculate the connectivity of the whole area (total size: 2873.8 ha), the equivalent connectivity (EC) of PC (probability of connectivity) probabilistic index was used (Saura et al. 2011a, b), where Euclidian edge-to-edge distances with 700 metres’ maximum dispersal distance threshold parameter (5% probability for a species to exceed 700 m) were set as a connection between habitat fragments. This distance was set as the average dispersal distance observed on the root vole by Steen (1994). All connectivity indices were calculated in Conefor 2.6 (Saura and Torné 2009). The maps of habitat fragments were created as digitalised topographic maps of Google Earth in map scale 1:5000. Individuals from the same locality were treated as a putative subpopulation (hereafter referred to as subpopulation).

Due to different trapping efforts at each study location, the size of each subpopulation was evaluated as the relative abundance (rA) of individuals captured at the location per 100 trap-nights (C/100TN), based on data without recaptures. The number of captured specimens (N) was recalculated into the transformed rA index assuming a random (Poisson) distribution (rA = (–ln (1-N/100) 100) of small mammals to remove the saturation effect caused by single traps (Caughley 1977). Relative abundance values were averaged over trapping occasions.

We tested the correlation between average rA of Pannonian root vole and average habitat fragment size using Pearson correlation analysis, where both data were log-transformed due to non-normal distribution.

DNA extractions were performed using commercial Isolate II Genomic DNA Kit (Bioline) according to the manufacturer’s protocol with the following modifications: during pre-lysis, samples were incubated overnight at room temperature; after adding preheated Elution Buffer G (70 °C), the elution step covered the incubation of samples at room temperature for 30 minutes and after that at 70 °C for 5 minutes before elution.

The thirteen microsatellite loci included in the analyses were developed for Microtus arvalis: Mar003, Mar016, Mar049, Mar063, Mar076 (Walser and Heckel 2008) and for Alexandromys oeconomus: Moe1, Moe2, Moe3, Moe4, Moe5, Moe6, Moe7 and Moe8 (Van de Zande et al. 2000). Microsatellite markers were grouped into four multiplex sets (Set1A: Moe1, Moe2, Moe3; Set1B: Moe4, Moe5, Moe6; Set2: Moe7, Moe8; Set3: Mar003, Mar016, Mar049, Mar063, Mar076).

Multiplex PCR reactions were performed in 12 μl volumes containing 2 μl (~80 ng) of DNA and a volume of 10 μl of the following mixture: 3.9 μl of KAPA2G Fast Multiplex Mix (KAPA Biosystems), 0.8 μl of BSA, 0.5 μl (concentration of 10 pm/ng) of each primer and RNase-free water to fill the volume to 10 μl. Amplification of DNA was carried out using peqSTAR 96X Universal thermal cycler (Peqlab).

To amplify microsatellites in Set 1A, Set 1B and Set 2 the PCR reaction consisted of the initial step at 94 °C for 7 minutes, 30 cycles including: denaturation at 94 °C for 1 minute, annealing at 55 °C in case of Set 2 and 60 °C in case of Set 1A and Set 1B for 2 minutes and extension at 72 °C for 90 seconds, followed by a final step at 72 °C for 10 minutes. Microsatellites in Set 3 were amplified in a PCR reaction including the following steps: 7 minutes at 95 °C, 35 cycles of 30 seconds at 94 °C, 1 minute at 57 °C and 90 seconds at 72 °C and after the cycles a final step of 10 minutes at 72 °C. To prepare the genotyping procedure 1 μl PCR product of each sample was mixed with 12 μl formamide and 0.3 μl GeneScan 500-LIZ size standard (Applied Biosystems). After a denaturation step of 5 minutes at 95 °C, a cooling step was implemented. Genotyping was carried out using ABI PRISM 310 Genetic Analyser (Applied Biosystems) and microsatellite genotypes were examined using GeneMapper software v.4.0 (Life Technologies).

We successfully genotyped 69 individual samples and the amplification success varied amongst markers (94.2–100%). The presence of null alleles may cause significant heterozygote deficit and deviation from the HWE. We therefore estimated the proportion of null alleles (NA) at each locus in each subpopulation using the programme FREENA (Chapuis and Estoup 2007). All genotypic distributions were in accordance with HWE expectations for all loci and, as we did not detect null alleles based on FREENA, consequently, all loci were included in the analyses.

Bayesian clustering of microsatellite genotypes was performed using STRUCTURE v.2.3.2 (Pritchard et al. 2000). Ten independent runs were performed for each value of K ranging from one to ten under a model assuming admixture and correlated allele frequencies (Falush et al. 2003). Each run comprised a burn-in period of 100,000 replications followed by a run length of 1,000,000 Markov Chain Monte Carlo (MCMC) iterations. The results of replicated runs for each value of K from one to ten were combined using STRUCTURE HARVESTER v.0.6.94 (Earl and vonHoldt 2012)and the optimal value of K was assessed by the inspection of log-likelihood values and according to the ΔK method developed by Evanno et al. (2005). Twenty independent runs were conducted with K fixed at the estimated optimal number of clusters where a burn-in of 100,000 and 1,000,000 MCMC iterations were used. The outputs of replicated runs were combined using the Greedy algorithm in CLUMPP v.1.1.2 (Jakobsson and Rosenberg 2007) and DISTRUCT v.1.1 (Rosenberg 2004) was used to visualise cluster assignments.

The mean number of alleles (A), observed (H_O) and expected (H_E) heterozygosity (Nei 1978) and inbreeding coefficients (F_IS) were calculated using GENETIX v.4.05.2 (Belkhir et al. 2004). The allelic richness (AR) was estimated using the rarefaction procedure for the lowest sample size (n = 10) in the programme FSTAT v.2.9.4 (Goudet 2003). Deviation of subpopulations from Hardy-Weinberg equilibrium (HWE) was calculated with exact tests assessing heterozygote deficiency and excess in GENEPOP v.4.2 (Rousset 2008). The basic level of significance was set to P = 0.05 and for multiple comparisons, we applied a Bonferroni procedure compensating for the risk of an inflating type 1 error.

The programme FREENA was used to estimate global F_ST, by performing 10,000 permutations. In addition, a Monte Carlo test of likelihood ratio G-statistic (Goudet et al. 1996) was performed using package hierfstat v.0.04-22 (Goudet and Jombart 2015) in R v.3.5.0 (R Core Team 2018) to test the presence of genetic structuring. Estimations of pairwise F_ST were implemented in FSTAT v.2.9.4 according to Weir and Cockerham (1984). Significant differences of F_ST estimators from zero were tested using 100,000 permutations without the assumption of Hardy-Weinberg equilibrium and significance level was set to P = 0.05 and adjusted using the Bonferroni correction. In addition, isolation by distance was tested by assessing the correlation between the geographical distance matrix (given in kilometres) and pairwise F_ST/(1-F_ST) estimates using Mantel’s Test with 10,000 permutations in the programme ISOLDE in GENEPOP v.4.2. The GENETIX v.4.05.2 package was used to investigate genetic relationships amongst all genotyped individuals of the nine subpopulations by factorial correspondence analysis (FCA), a method that identifies the linear combination of variables (allele frequencies at different loci) that captures the most variation between observations (individuals or populations) and visualises genetic relationships in 2D space determined by FCA axes. Analysis of Molecular Variance (AMOVA) (Excoffier et al. 1992) was implemented in GENALEX v.6.5 (Peakall and Smouse 2006, 2012) to test the proportion of genetic variance amongst individuals and subpopulations. Statistical significance of the variance components was assessed with 999 permutations.

Between 2004 and 2019, we identified 26 permanent habitat fragments as suitable habitats for root voles, in 14 of which their presence were confirmed. We have also identified several temporarily suitable fragments, the number of which varied seasonally and annually. All measured characteristics of potential habitats (number of fragments and temporary fragments, total area of all habitat fragments, mean fragment size and connectivity of the whole area) changed over the six time periods of study (see Table 1). In general, connectivity had higher values in periods when the number and size of habitat fragments were higher. In March 2004 and March 2011, the number of temporary fragments increased, filling the matrix of agricultural land between permanent fragments with new suitable habitat (see Fig. 2). The highest values of parameters were observed in 2011, except for the average fragment size. After this period, the total area of habitat fragments, the number of temporary fragments and the connectivity decreased and permanent fragments were thus more isolated. Values of parameters remained relatively similar in the following years. Partial changes in habitat fragment characteristics and connectivity are detectable in seasonal comparison between March and August 2017, with higher values in spring. In a detailed view of the study sites in Fig. 2, it can be seen that, while in August 2017, sites SK3, SK4, SK5 and SK6 were relatively small and isolated from each other, in March 2004, 2011 and 2017, the same sites were larger and only a channel interrupted their direct connection.

Figure 2.

Change in the number of suitable habitat fragments for the Pannonian root vole in the north-western part of the study area over time. The green box in the insert shows the boundaries of the larger maps.

The average occupancy of fragmented habitats by the Pannonian root vole varied spatially (Fig. 3). While some fragments had a high relative population size, others had very low. The analysis did not confirm a positive correaltion between average rA of Pannonian root vole and average size of fragmented habitats (r = 0.15, N = 9, P = 0.74).

Figure 3.

The relative abundance (rA) of Pannonian root vole subpopulations in sampled habitat fragments. Relative abundance values were averaged over trapping occasions and average rA was visualised on the map as graduated symbol size.

Genetic structuring inferred from STRUCTURE analysis is presented in Fig. 4 for K values ranging from 2 to 5. The ΔK method indicated that the optimal number of clusters was 5. However, this method can evaluate ΔK only for K > 1 and the inspection of log-likelihood scores revealed that K = 1 had a relatively high likelihood, similar to other values of K. Generally, assignment probabilities were roughly symmetric, which indicated that none of the clustering analyses captured the real population structre. Notably, when K was fixed at 5, SK1 individuals sampled in 2014 (individuals 3, 6, 7 and 8) and an additional individual from SK5 (individual 26 sampled in 2015) had relatively high assignment probabilities to the same cluster and the same for individuals 57 and 60 sampled in SK9 in 2017. However, not all individuals sampled in 2014 and 2015 were assigned to this cluster (indicated with grey in Fig. 4).

Figure 4.

Genetic structure of the sampled root vole subpopulations in the Danubian Lowland. The graph is based on STRUCTURE runs when K was fixed at 2–5. Each individual is represented by a line proportionally divided into colour segments corresponding to its membership in certain clusters. Black lines separate the individuals from different habitat fragments.

Clustering analysis did not give a strong evidence of structuring; therefore, measures of genetic diversity were calculated for all subpopulations separately and for all samples pooled together. Genetic diversity and HWE were not calculated for locations SK3, SK4 and SK7 due to their small sample size (N = 3).

The number of alleles per locus in subpopulations ranged from 2 (locus Moe4) to 16 (locus Moe7), and the mean number of alleles per locus (A) ranged from 4.31 to 5.62 (Table 2). Allelic richness (AR) across subpopulations ranged from 4.22 to 5.46, being the highest in the SK1 subpopulation and lowest in the SK8 subpopulation. H_O (range between 0.562 and 0.669) was the highest in subpopulations SK8 and SK9, while the lowest value was shown for SK5. In contrast, H_E (range between 0.624 and 0.657) showed another pattern, with the highest value for SK1 and the lowest value observed for SK8. Two subpopulations, SK8 and SK9 had heterozygosity excess, based on observed and expected heterozygosity values. Subpopulations SK1, SK2, SK5 and SK6 showed significant heterozygote deficiency on the basis of F_IS (significant positive values). Subpopulations SK1, SK2 and SK5 deviated from HWE according to heterozygote deficiency exact tests in GENEPOP (P = 0.05) after Bonferroni correction (Table 2). Heterozygote excess exact tests were not significant in any of the cases (data not shown). Mean number of alleles and expected heterozygosity indicated higher genetic diversity and the exact test showed a significant deficiency of heterozygotes when all samples were pooled together.

Subpopulations SK3, SK4 and SK7 were not included in F_ST analyses because of their small sample size (N = 3). The global F_ST for six subpopulation samples was 0.025 (95% CI: 0.01–0.041). The overall G-test was significant (P < 0.001), indicating genetic structuring amongst locations. Pairwise F_ST values were relatively low, although variable. The highest F_ST values were observed for SK8 and SK9. Pairwise genetic differentiation was not significant in most of the comparisons, except in cases where one subpopulation of the pair was always SK6, SK8 or SK9 (Table 3). Significant F_ST was observed between SK1 and SK8, but when the individuals of SK1 sampled in 2014 were excluded from pairwise F_ST analyses, this difference became non-significant (F_ST = 0.048). Microsatellite-based genetic distances did not correlate with geographical distances between subpopulations (a = 0.0254, b ~ 0, P = 0.6), thus spatial separation (IBD) was not confirmed.

The FCA plot, based on individual genotypes, clearly separated SK1 along the first factorial axis (explaining 20.1% of variation) from all other subpopulations. The second axis (explaining 17.9% of variation) mainly separated the individuals from SK1 and SK9, while individuals from SK6 showed only a weak segregation (Fig. 5A). Subopulation SK8 was separated along the third axis (explaining 14.6% of variation, Fig. 5B).

Figure 5.

Two-dimensional plots of FCA performed for nine subpopulations showing the 1^st and 2^nd (A) and the 1^st and 3^rd (B) axes. The proportion of explained variance is written in parentheses on each axis.

In the Analysis of Molecular Variance, significant genetic variation was attributed to the differences between subpopulations (4.4%, P < 0.001) and most of the variability occured within subpopulations (95.6%).

The number, size and shape of habitat fragments in the studied region varies in time as the result of exogenous factors (precipitation, surface water levels, agricultural activities). These dynamic changes have an effect on fragment connectivity, suggesting that connectivity was positively influenced by the number of habitat fragments. As we have also noted, permanent habitat fragments, relatively distant from each other at one time, can change size and shape and become neighbouring habitats at another time. In addition, the temporary fragments can play the role of stepping stones during vole movements. Thus, despite the constant presence of habitat fragments and channel-side vegetation, fragment connectivity can vary seasonally and yearly, as can change the possibility of individuals’ replacement between the studied subpopulations. In Norway, root voles increased dispersal distance as a response to fragmentation, but it was less affected by connectivity (Bjørnstad et al. 1998). Detailed data on the dispersal ability of the Pannonian root vole are still missing, although they have been shown to be able to cross barriers when the quality of original habitats deteriorated and they may follow stepping stones when searching for new habitats (Horváth and Herczeg 2013). Moreover, Kratochvíl and Rosický (1955) observed that Pannonian root voles in sedge (Carex sp.) survived on hummocks when the water level increased (50–70 cm). When the water level decreased, voles spread across the lower parts of these hummocks, while some individuals dispersed to new suitable fragments. In view of this, a possible explanation for the observed genetic patterns could be that the animals found shelter on hummocks during floods in 2010 and probably 2013, thus opening the way for admixture. Later, when water levels fell, presumably admixed individuals resettled lower areas.

Levels of genetic differentiation between the habitat fragments varied, but were mostly non-significant and we found no support for isolation by distance between subpopulations. Bayesian clustering in STRUCTURE did not reveal pronounced genetic structuring, indicated by approximately equal assignment probabilities to different clusters in all cases of K from 2 to 5. This result suggested that dispersal effectively prevents marked genetic subdivision between studied habitat fragments, which can be additionally confirmed by the lack of isolation by distance between fragments. Given the small geographical scale and landscape pattern of the study area, we would expect gene flow between localities to maintain very low or no differentiation between subpopulations. In a study conducted in the Netherlands, Mauritzen et al. (1999) concluded that ditches are likely to favour linear movements of root voles and may enhance connectivity in an agricultural landscape. Most habitat fragments sampled in our study were directly connected by reed- and sedge-lined channels that are likely to serve as dispersal corridors and permanent or temporary habitat fragments in the study area also likely to facilitate dispersal (Ambros et al. 2016). Although the longest geographical distance along the channels separating our sampling sites from their nearest neighbours was nearly 6 km, as for SK1 and SK2, the presence of root vole was confirmed in other permanent habitat fragments between these two sampling sites (which could not be included in our study due to the small number of samples) and most of the habitat fragments sampled were relatively close to each other. Movement distances of male root voles can exceed several hundred metres in a short time or more than two kilometres within a few days (Steen 1994; Andreassen et al. 1996).

Consistent with the changing possibility of individuals’ replacement between the studied fragments, AMOVA results also showed a low, but significant, percentage of variability between subpopulations. In addition, signs of genetic differentiation were detected between subpopulations SK8, SK9 and SK6, based on significant pairwise F_ST values and the FCA analysis confirmed the separation of these samples. Results may therefore indicate that dispersal is not unhindered between all subpopulations and root vole individuals in the network of studied habitat fragments may not be viewed as a panmictic population.

In a detailed study, Hulejová Sládkovičová et al. (2018) already noted that, despite the high genetic variability of the subspecies, local populations of Pannonian root vole may show signs of genetic depletion. Although our analyses were limited by small sample sizes due to the rarity of the Pannonian root vole, our results indicate lower genetic diversity in subpopulations (in terms of allelic richness, mean number of alleles per locus and expected heterozygosity) than what was measured in other populations from south-western Slovakia (Hulejová Sládkovičová et al. 2018). Similarly, reduced genetic diversity (based on the mean number of alleles per locus and the expected heterozygosity) was observed compared to the populations of another isolated Western European root vole subspecies, namely A. o. arenicola (Van de Zande et al. 2000). Expected heterozygosity and mean number of alleles in the total studied population indicated decreased variability compared to other Slovakian populations (Hulejová Sládkovičová et al. 2018). Low genetic diversity of the studied root vole population putatively reflects a pronounced geographical isolation of the surviving population and the decrease in effective population size that has likely occurred. In addition, flood events may have reduced the size of the vole population (Kratochvíl and Rosický 1955), which may have created a bottleneck, although its influence on the genetic diversity of the studied population cannot be clearly confirmed on the basis of our results alone. The reduced genetic diversity is in accordance with the view of Hulejová Sládkovičová et al. (2018) who suggested that the Danube River and associated branch network could facilitate vole dispersal, but local populations embedded in the agricultural landscape and further from the Danube are likely to remain isolated. Moreover, based on the evaluation of the distribution of Pannonian root vole, our study area is close to the eastern boundary of the subspecies’ range (Gubányi et al. 2009; Ambros et al. 2016) and populations occurring here are likely isolated by the distance from other populations to the west.

Subpopulations SK8 and SK9 tend to have lower levels of allelic richness, which is consistent with the possibly lower probability of dispersal through the agricultural land matrix compared to other sites. However, we did not find significant deviation from the Hardy-Weinberg equilibrium in these subpopulations. We observed deviations from the Hardy-Weinberg equilibrium in SK1, SK2 and SK5 and significant positive F_IS coefficients in the same subpopulations, together with SK6, which may result from the social structure of root voles. Matriline-based groups in root vole populations (Tast 1966) give the possibility of sampling kin individuals, which may result in departures from the Hardy-Weinberg equilibrium (Aars et al. 2006). On the other hand, Stewart et al. (1999) argue that genetic composition (and departures from the Hardy-Weinberg equilibrium) of water vole (Arvicola amphibius) colonies is under the influence of yearly fluctuation and suggested that genetic drift is a plausible cause resulting in that pattern. They interpreted high levels of genetic diversity as a sign of gene flow between colonies. Similarly, despite the relatively large number of analysed individuals, Pilot et al. (2010) observed yearly changes in departures from the Hardy-Weinberg equilibrium in a root vole population in Poland; nevertheless, a high level of genetic diversity was maintained over the years. In addition, they found no genetic differentiation between years regardless of changes in density, kin structure and deviations from HWE, indicating that genetic composition is stable over time in large, demographically-stable populations. In contrast, studies of water vole metapopulations revealed temporal genetic differentiation, but genetic variability remained high as a result of gene flow (Stewart et al. 1999; Aars et al. 2006).

Regarding individuals 3, 6, 7 and 8 (sampled in 2014 and 2015) in SK1, highlighted by clustering analyses, their high assignment probability to a separate cluster may reflect their distinct origin. Habitat fragment SK1 is the closest to the Danube River amongst the studied fragments and the floods in 2010 or 2013 potentially facilitated dispersal from further areas and the aforementioned individuals might be immigrants or their descendants. Other individuals that had relatively high assignment probabilities to the same cluster were captured in SK5 (individual 26 from 2015) and SK9 (individual 57 and 60). These can be found at a few kilometres distance from SK1, but given the small spatial scale, it is not unlikely that these specimens may be the offspring of dispersing individuals. Alternatively, it is also possible that genetic drift over time changed the genetic composition of subpopulations in the fragmented landscape; hence, some (but not all) samples collected in 2014 and 2015 were highlighted by STRUCTURE clustering. However, we are not able to declare which possibility is more plausible without genotyping individuals from other areas and without temporal analysis of samples.

The genetic diversity of the local subpopulations in the study area is relatively low; although the studied subpopulations probably have connections with each other due to the effect of extensive floods and the network of fragments and channels in the agricultural landscape, their reduced genetic variability is detectable compared to the pooled genetic variability of other populations of Pannonian root vole occurring closer to the more uninterrupted marshlands in Szigetköz, Hanság and Neusiedlersee Regions (Hulejová Sládkovičová et al. 2018). Based on the pattern of remnant habitat fragments, it is likely that root vole subpopulations in our study area form a metapopulation system and the observed genetic patterns do not contradict this possibility. However, studies on extinction and recolonisation patterns and migration rates would be essential to support this view (Gaggiotti 2004; van der Merwe et al. 2016).

Landscape changes and habitat destruction resulted in the fragmented distribution of root vole habitats in the study area and fluctuating surface water levels induce considerable changes in habitat size, quality and connectivity to this day. Only one fragment (SK2) in our study area is protected as a Special Protection Area. However, for the long-term persistence of root vole populations, it would be critical to ensure legal protection of habitats. The importance of protected core areas has been demonstrated for water vole metapopulations (MacPherson and Bright 2011) and mainland or source population demography has been shown to shape habitat use of root voles (Glorvigen et al. 2013). Based on the low differentiation, it seems unlikely that the studied vole subpopulations formed a mainland-island type metapopulation system (Stewart et al. 1999), but the signs of admixture highlight the vital role of habitat fragments and gene flow in the maintenance of genetic diversity. On the other hand, Van de Zande et al. (2000) proposed that weak differentiation may be an indicator of progressive isolation of local populations and accordingly, the prevention of fragmentation should be considered in the conservation management of A. o. arenicola. It has been suggested for other vole species and subspecies that persistent gene flow between habitat fragments is crucial for the maintenance of genetic variability (Telfer et al. 2003). Therefore, it is necessary to protect not only habitat fragments with current occurrences of the Pannonian root vole, but also other potentially-suitable habitat fragments, to form a habitat network and to ensure the possibility of gene flow (Neuwald 2010). Water levels in the study area and, as we demonstrated, total size of permanent habitat fragments and number of temporary habitats are in dynamic change, which further emphasises the contribution of these fragments to the survival of the subspecies in the agricultural landscape (Kalmár and Riezing 2017).

It is expected that the overall genetic diversity of the subpopulations will decrease as a result of their small size and isolated location in the agricultural matrix. This implies that the restoration of habitats and corridors is indispensable for the long-term preservation of diversity, as has been stressed earlier (Thissen et al. 2015a; Hulejová Sládkovičová et al. 2018). As a habitat specialist, the Pannonian root vole is sensitive to habitat quality changes, which are likely to be the key factors determining habitat use (Glorvigen et al. 2013; Horváth and Herczeg 2013). Thus, it is important to preserve optimal water level conditions for the root vole, for example, by controlled construction of infrastructure that may cause aridification and the disappearance of the subspecies (List et al. 2010; Kalivodová et al. 2018). For this reason, it is necessary to take the needs of the subspecies into consideration when infrastructural projects are planned and managers should be involved in landscape and land-use planning. Finally, long-term monitoring is recommended to provide further information about population and subpopulation parameters (Hayes et al. 2017), genetic diversity and structure (Proença-Ferreira et al. 2019) and on the changes of habitats (Martensen et al. 2017) to elaborate an effective conservation management action plan that ensures the long-term survival of the Pannonian root vole.