The genetic structure and connectivity in two sympatric rodent species with different life histories are similarly affected by land use disturbances

The negative impact of habitat fragmentation due to human activities may be different in different species that co-exist in the same area, with consequences on the development of environmental protection plans. Here we aim at understanding the effects produced by different natural and anthropic landscape features on gene flow patterns in two sympatric species with different specializations, one generalist and one specialist, sampled in the same locations. We collected and genotyped 194 wood mice (generalist species) and 199 bank voles (specialist species) from 15 woodlands in a fragmented landscape characterized by different potential barriers to dispersal. Genetic variation and structure were analyzed in the two species, respectively. Effective migration surfaces, isolation-by-resistance (IBR) analysis, and regression with randomization were used to investigate isolation-by-distance (IBD) and the relative importance of land cover elements on gene flow. We observed similar patterns of heterozygosity and IBD for both species, but the bank vole showed higher genetic differences among geographic areas. The IBR analysis suggests that (i) connectivity is reduced in both species by urban areas but more strongly in the specialist bank vole; (ii) cultivated areas act as dispersal corridors in both species; (iii) woodlands appear to be an important factor in increasing connectivity in the bank vole, and less so in the wood mouse. The difference in dispersal abilities between a generalist and specialist species was reflected in the difference in genetic structure, despite extensive habitat changes due to human activities. The negative effects of fragmentation due to the process of urbanization were, at least partially, mitigated by another human product, i.e., cultivated terrains subdivided by hedgerows, and this was true for both species.


Introduction
Habitat loss and fragmentation have negative impacts on populations, and are considered as one of the main causes of biodiversity loss and therefore a major issue in conservation biology (Fischer and Lindenmayer 2007;Wilson et al. 2016;Wu 2013). In particular, anthropogenic habitat fragmentation has modified the distribution and population sizes in many different organisms (Crooks et al. 2017;Haddad et al. 2015) with local and/or global reduction of genetic diversity and connectivity (DiBattista 2008;Leigh et al. 2019). Monitoring the genetic consequences of human activities that increase habitat fragmentation is therefore important to develop appropriate conservation and management strategies (Hoban et al. 2020).
The major consequence of habitat loss and fragmentation is to create discontinuities (i.e., patchiness) in the distribution of critical resources (e.g., food, habitat, water) or in environmental conditions (e.g., microclimate) (Segelbacher et al. 2010). Such discontinuities reduce connectivity among populations (Kindlmann and Burel 2008), threatening their long-term viability due to genetic (e.g., reduced evolutionary potential and inbreeding depression) and demographic factors (e.g., demographic stochasticity) (Balkenhol and Waits 2009). Habitat fragmentation may also have different short-term consequences in different species, for example, by reducing the suitable habitats or increasing the predation success, but these effects poorly predict long-term responses (Evans et al. 2017). Gene flow among subpopulations is necessary to alleviate the adverse genetic consequences of population fragmentation, reducing genetic drift and maintaining local genetic variation (Frankham et al. 2017). From a conservation perspective, inferring the functional connectivity of populations across landscapes becomes crucial (Segelbacher et al. 2010). Identifying the areas where gene flow is either facilitated or prevented, and the landscape factors responsible for that, is a high priority (Sork 2016;Trombulak and Baldwin 2010).
One interesting opportunity to investigate the causes and the genetic consequences of fragmentation is represented by sympatric species with partially overlapping ecological niches (Homola et al. 2019;Varudkar and Ramakrishnan 2015). Different species, in fact, may respond very differently to the same landscape matrix (Frantz et al. 2012;Gangadharan et al. 2016;Nevill et al. 2019;Robertson et al. 2018;Thatte et al. 2020). They may also react differently to the fragmentation of their previously continuous habitat, and these differences may be reflected in the geographic distribution of their genetic variation. In this work, we investigate the effects of habitat fragmentation present in agricultural landscapes in Central Italy on the genetic structure of two sympatric rodent species, the wood mouse (Apodemus sylvaticus) and the bank vole (Myodes glareolus).
The wood mouse is a generalist species known to inhabit a wide range of habitats including forests, hedgerows and agricultural fields (Marsh and Harris 2000;Montgomery and Dowie 1993;Mortelliti et al. 2010). In contrast, the bank vole is a "forest specialist", i.e., it is more strictly associated with forest habitats, from mature stands to recently coppiced woodlands (Capizzi and Luiselli 1996;Ecke et al. 2002). In general, specialist species tend to be more affected by habitat fragmentation than generalist species. In fact, highly dispersed resources are harder to reach by the specialist species (Braschler and Baur 2005;Nupp and Swihart 2001;Youngentob et al. 2012), which also suffer from competitive exclusion by the generalist species (Sozio and Mortelliti 2016). Accordingly, the specialist bank vole seems to prefer sites with high connectivity (Mortelliti et al. 2009;Sozio and Mortelliti 2016), and the generalist wood mouse can also be found in highly fragmented habitats, being able to move across cultivated fields (Sozio and Mortelliti 2016;Sozio et al. 2013). We currently do not know whether these differences directly correspond to a stronger genetic structure in the bank vole compared to the wood mouse, and if (and how) different natural or anthropogenic habitat features have different relative impacts on gene flow. Our study addresses these questions following three steps: (i) neutral genetic markers will be used to estimate the genetic diversity and the population structure separately in each species; (ii) patterns of increased and reduced gene flow will be analyzed and compared between species; and (iii) species-specific landscape features with the largest influence on genetic variation patterns will be identified.
Woodland patches, consisting of mixed deciduous forest dominated by downy and turkey oaks (Quercus pubescens and Quercus cerris), were embedded in an agricultural matrix (mainly wheat fields) crossed by a network of hedgerows providing structural connectivity to habitat patches. The S2 highway and a railway bisect the study area, potentially acting as barriers to wildlife movements (Grilo et al. 2009). Urban areas are also present and represent approximately 5% of the total area. Twelve trapping sessions were conducted over a 2 year period, with bi-monthly trapping from April 2011 to February 2013. During each session, grids were trapped for three consecutive nights. A total of 199 bank voles and 194 wood mice samples were obtained. Sample sizes for each of the 15 different woodland patches is reported in

Genotyping
Genomic DNA was extracted from the mouse ear lobe samples using the NucleoSpin® Tissue kit (Macherey-Nagel,

Genetic diversity
Descriptive statistics of nuclear genetic diversity were estimated separately for each population (woodland patch) in each species. The mean number of alleles, and the observed and expected heterozygosities, were estimated using Genalex 6.5 (Peakall and Smouse 2006), and the same program was used to test for deviation from Hardy-Weinberg equilibrium. Allelic richness (AR) was calculated using the rarefaction procedure in the Fstat 2.9.4 software (Goudet 2001). Arlequin 3.5.2.2 (Excoffier and Lischer 2010) was used to test for linkage disequilibrium between each pair of loci for each sampling population following a likelihood-ratio statistic, whose null distribution was obtained by a permutation procedure. We applied sequential Bonferroni corrections to account for multiple comparisons (Rice 1989 (Weir and Cockerham 1984) was estimated for each pair of sampling population with Arlequin. Statistical significance of the F ST values was tested using 10000 permutations, and P values were multiplied by the total number of comparisons following the conservative Bonferroni approach for multiple testing.

Genetic structure
A Bayesian clustering method was used to identify the number of genetic groups (STRU CTU RE v2.3.4; Pritchard et al. 2000). A burn-in length of 50,000 iterations and a run length of 100,000 iterations were used in an admixture model with correlated allele frequencies among populations testing each K value between 1 and 15. Each K value was run 10 times. The optimal K value was determined using the ∆K method (Evanno et al. 2005) implemented in the online tool STRU CTU RE Harvester (Earl and VonHoldt 2012). CLUMPP (Jakobsson and Rosenberg 2007) was then applied to average the multiple runs given by STRU CTU RE and to verify correct label switching. To display the results, the output from CLUMPP was visualized with DISTRUCT (Rosenberg 2004).

Visualizing deviation from Isolation by Distance
Genetic diversity between populations often exhibits patterns consistent with Isolation-by-Distance (IBD) (Wright 1943), where populations far apart in the geographic space receive less gene flow than neighbouring ones. Given the ubiquity of this phenomenon (Kuchta and Tan 2005;Sharbel et al. 2000) it is interesting to see locations where this does not hold true, as they might represent barriers or zones of high contact. Global deviation from IBD can be identified, for example, studying the decrease of similarity or autocorrelation with geographic distance. However, specific deviations in some areas, but not in others, cannot be easily investigated and visualized by standard methods. One recent solution to this problem comes from the use of Estimated Effective Migration Surfaces (EEMS; Petkova et al. 2015). EEMS employs individual based migration rates in order to visualize zones with higher or lower migration with respect to the overall rate. These areas represent locations in which the pattern of gene flow predicted by IBD is facilitated or hindered. The region under study was first divided in a grid of demes and the individuals were assigned to the deme closest to their sampling location. The matrix of effective migration rates was then computed by EEMS based on the stepping-stone model (Kimura and Weiss 1964) and on resistance distances (McRae 2006). We used the EEMS script for microsatellites analysis runeems_sats available from Github (https:// github. com/ dipet kov/ eems) to construct EEMS surfaces for the bank vole and the wood mouse. Considering that the number of demes simulated during the grid construction phase can influence the scale of the deviation from the overall migration rate, we averaged three runs with 50, 100, 200, 300 and 400 demes to produce the final EEMS surface. Each single run consisted of 200,000 burn in steps followed by 1,000,000 MCMC iterations sampled every 1 3 10000 steps. We plotted the averaged EEMS and checked for MCMC convergence using the rEEMSplots package in R.

Isolation by resistance
Understanding the effect of environmental components on the genetic makeup of natural populations is the goal of landscape genetics, which integrates population genetics, landscape ecology and spatial statistics (Holderegger and Wagner 2008;Manel et al. 2003;Storfer et al. 2007). One of the techniques more commonly used in landscape genetics to identify discontinuities in gene flow and determine the relative resistance to movement imposed by different landscape elements is Isolation-by-Resistance (IBR; McRae 2006). IBR offers a conceptual model in which landscape resistance is the analogue of electrical resistance, and the movements of individuals and flow of genes are analogues of electrical current (Amos et al. 2012). It greatly extends the ability to model multiple complementary paths of connectivity, while being sufficiently computationally efficient to allow its use over large landscapes at relatively fine resolution (McRae and Beier 2007;McRae et al. 2008). In order to analyze the effect of specific landscape components on gene flow, we tested for the presence of IBR. We first constructed a raster grid encompassing all our study area reclassifying the land cover based on features that were a priori most likely to affect gene flow in both the bank vole and the wood mouse: woodland, urban areas, cultivated terrain and hedgerows ( Fig. 1). We obtained a raster file encompassing the entire studied area with a resolution of 10 m that was previously classified based on major land cover types (CORINE Land Cover; Büttner et al. 2004). Aerial photographs were used to add the locations of hedgerows with the same resolution (Sozio and Mortelliti 2016). We also included in our raster grid the major roads intersecting our study area from Open-StreetMap (http:// www. opens treet map. org) and the railways tracks from the DIVA-GIS database (http:// www. diva-gis. org/ gdata). In order to determine the relative importance of land cover elements in hindering or facilitating gene flow, we modified the created raster under two different scenarios. The first set (resistance set) was aimed at determining the resistance caused by a specific land cover feature with respect to the others. We assigned a varying maximum resistance (RE max ) to a target component, keeping the other landscape features to a uniform minimum resistance (RE min = 1). The second set of grids (permeability set) was built to establish the possible role of a specific landscape feature in facilitating the connection between different populations. Due to the challenges in determining the effective resistance posed by each environmental feature to gene flow, we employed a optimization procedure encompassing a wide range of resistance values (Ortego et al. 2015;Roffler et al. 2016). We assigned a minimum resistance value to a target landscape component and a varying RE max to all remaining features. For both set of grids we employed eight maximum resistance values (RE max = 5, 10, 50, 100, 500, 1000, 5000 and 10000) obtaining a total of 96 different surfaces. We computed pairwise resistance distances between populations for both the bank vole and the wood mouse using the different sets of grids. Distances were obtained considering the eight-neighbour cell connection scheme in CIRCUITSCAPE v4.0 (McRae, Shah, and Mohapatra, 2013) with the sampled woodland patches as focal regions. We also computed an IBD scenario considering a homogeneous resistance surface (all RE = 1) (Castillo et al. 2014;Ortego et al. 2015). We then compared the resistance and the F ST matrices using multiple matrix regression with randomization (MMRR) (Wang 2013). We considered a consistent trend of increase of R 2 , at least for 4 levels of resistance values, with P values always smaller than 0.05 in single tests, as evidence of the impact of the corresponding land cover feature. For each landscape variable, the most supported model was identified as the one corresponding to the highest supported R 2 value. In case of plateau, we preferred the model corresponding to the onset of the plateau (Castillo et al. 2014). Statistical significance of the coefficients was determined using 9999 permutations with the MMRR function (Wang 2013). All statistical analyses were conducted in R.

Genetic diversity
All loci were polymorphic in both species. The average expected heterozygosity was very similar in the two different sets of markers typed in the two species (0.74 in the bank vole and 0.72 in the wood mouse), and the number of alleles varied between 2 and 16 in the wood mouse and between 3 and 11 in the bank vole, respectively. All the genetic variation statistics are reported in Table 1.
No systematic deviation from linkage equilibrium was observed between loci for any population in both species, and none of the pairwise tests was significant after Bonferroni correction. Some loci showed evidence of the presence of null alleles, but only in some populations. We analyzed the effect of these alleles by comparing matrices of pairwise F ST values computed from the complete data set with values corrected for null alleles as estimated by FreeNA. Multilocus global F ST values had identical values when calculated with and without correcting for null alleles in both species (wood mouse: F ST = 0.03, with or without correction; bank vole: F ST = 0.08, with or without correction), with identical or very similar confidence intervals in the two analyses (0.01-0.05 in wood mouse, with and without correction, 1 3 0.07-0.09 and 0.06-0.08 in bank vole, with and without correction, respectively). Multilocus pairwise F ST values with and without correction were also highly correlated (wood mouse: r = 0.99; p = 0.001; bank vole: r = 0.99; p = 0.001; Mantel test). We decided therefore to use the complete data set for all downstream analyses. Pairwise F ST values in the wood mouse were significant after sequential Bonferroni correction only in 7 out of 105 comparisons, all involving the PRV population (with F ST values never larger than 0.08) ( Fig. 2; Table S2). On the contrary, the bank vole showed a finer geographic structure. Approximately half of the F ST values were significant, with the highest divergence values observed in comparisons including PRV, and, as reported above, the average F ST was much higher than that estimated in the wood mouse ( Fig. 2; Table S2).

Genetic structure
The most likely partition implied three genetic groups (K = 3) in both species. Here we present individual assignment plots for K equal to 2, 3, and 4 ( Fig. 3a, b) to better visualize different aspects of the genetic structure, and we also report the geographic distribution of the most supported clades in both species (Fig. 3c). In the wood mouse (Fig. 3a), the isolation of PRV already suggested by the pairwise F ST matrix was supported at different values of K. With the most supported K = 3, or with K = 4, a large fraction of individuals and populations (with the exception of PRV) showed a mixed ancestry. In the bank vole (Fig. 3b), populations appeared more internally homogeneous, with three distinct genetic groups prevailing in the northern areas (ALB, BRN, FDT, FRR and GST), in the western areas (API, IUG, MCD, PRV and YAH), and in a single eastern population (CRC), respectively, and the other populations having a more mixed and less geographically localized genetic composition. The results are robust with increasing burn-in and replicates.

Visualizing deviation from IBD
The spatial visualization of the geographic areas with higher or lower gene flow compared to IBD expectations was similar in the two species (Fig. 4). The main pattern consisted of a central area of reduced gene flow, centered around PRV, extended only in the bank vole towards the southern and the eastern borders of the region. These branches of reduced migration clearly produced the higher genetic structure observed in the bank vole when compared to the wood mouse, with the latter having a much higher connectivity in most of the areas we considered.

Isolation by resistance (IBR)
Both the wood mouse and the bank vole populations showed in the IBR a significant pattern of IBD (Tables S3-S6). However, we also found consistently higher association between pairwise F ST and resistance distance in some models including land cover features ( Fig. 5; Tables S3-S6). Highest R 2 values in the best models were always smaller in the wood mouse, reflecting probably the fact that gene flow in relation to resistance (a, b) and permeability (c, d) distance matrices plotted against resistance values for different landscape features. Circles with black outline showed significant P-values to produce a small additional fit to the model. The second set of resistance values (permeability), designed to highlight factors facilitating connectivity between different populations, supported the role of cultivated areas in both species, with a 50% increase of R 2 compared to IBD at RE = 1/100 (see Fig. 5c,d; Tables S5-S6). In addition, woodlands seem also to play a role especially in the bank vole, and areas comprising and surrounding major roads increased the fit of the model only in the wood mouse. A summary of the major results obtained in the two species is in Table 2.

Discussion
In this study we investigated the relationships between human-related changes in habitat amount and configuration (i.e., habitat structure), habitat use and genetic structure. We applied an identical sampling scheme within the same fragmented area to two sympatric rodent species, the wood mouse and the bank vole. Our main result is that the same factors appear to limit (urban areas) and favour (woodland and cultivated areas) gene flow in both species despite their difference in ecological niche. The genetic structure in the same region before human intervention is not known, but it is interesting to note that the joint and compensatory effects of urban and cultivated areas did not modify the relative levels of genetic structure expected in the two species: the generalist wood mouse has a population structure much more genetically connected than the forest-specialized bank vole. More specifically, gene flow favored by cultivated areas likely increases the genetic exchanges in the wood mouse even above the high level expected in natural conditions, now limited only by urban areas. In the bank vole, on the other hand, cultivated areas possibly act compensating the genetic fragmentation due to the loss of the preferred niche (woodland) and the increase of urban areas. Any further replacement of woodlands with urban and not cultivated areas will clearly reduce this compensatory effect. Overall, we conclude that the difference between wood mouse and bank vole is still reflected in the difference between their current genetic structure, but if woodlands will be further replaced by urban settlements and not cultivated areas this difference will likely increases.

Genetic diversity
Habitat fragmentation did not produce a detectable loss of genetic variation in the two species. Levels of diversity in different populations are comparable to those reported for other rodent species (Czarnomska et al. 2018;Dominguez et al. 2021;Gerlach and Musolf 2000;Martin Cerezo et al. 2020). When the global genetic divergence between populations was analyzed, the wood mouse showed much weaker population structure than the bank vole. This pattern is expected considering that, at a short geographic scale (distances < 30 km), genetic structure is commonly found only in rodents with a specialized ecological niche (Bani et al. 2017;Fasanella, et al. 2013;Kozakiewicz et al. 2009;Łopucki et al. 2022).
With the exclusion of the population sampled in PRV (see below), the wood mouse appears rather homogenous at this geographic scale, indicating that gene flow was not prevented by the human-induced fragmentation of their natural habitat. This result reflects the enormous capacity of adaptation and mobility in this species, which can be found in all types of forests and even in cultivated fields during certain periods of the year (Harris and Woollard 1990;Szacki et al. 1993;Tew 1994). Conversely, populations of the bank vole sampled in the same patches showed the presence of a significant genetic differentiation with a lower degree of genetic admixture and higher F ST values. Similar studies on bank vole confirmed that there is a significant reduction of gene flow already at geographical distance of about 8 km (Stacy et al. 1997), and that environmental features, such as seasonal temperature variations, can contribute in a decisive way in increasing the genetic structure of this species (Gȩbczyński and Ratkiewicz 1998).
Fragmentation of the natural environment as a consequence of human activity (e.g., urbanization, agriculture, pastures) can reduce the genetic diversity of populations by decreasing population size and gene-flow, while increasing genetic drift. Generalist species can be considered as "urban adapters" (Blair 2001), since they can occupy novel environments/ecosystems given their wider niche breath. In white-footed mice (Peromyscus leucopus), urbanization increased population differentiation while suffering minimal loss in genetic diversity (Munshi-South and Kharchenko 2010) and keeping IBR between patches of urban vegetation (Munshi-South 2012). Our results also show a similar pattern in which both species retained their genetic diversity and showed population structure. While we cannot assess the presence of historical population differentiation prior to environment fragmentation, we show that different habitat disturbance similarly affected sympatric species with different ecological specializations.

Spatial patterns of gene flow
Isolation-by-distance was significant for both species indicating that geographic distance is an important factor for isolating the various populations. An additional shared feature appears to be the isolation of PRV in all analyses, supporting the hypothesis that individuals of both species have difficulty reaching this area. This result may be related to the fact that urban areas are highly diffused around PRV, and the IBR analysis suggested that they act as a barrier for both species.
For the bank vole, we found additional areas of enhanced or reduced gene flow, in comparison with the IBD pattern in the background. Specifically, three main areas showed higher gene flow than expected, corresponding to western, eastern and northern patches. Barriers separating them are composed of a mix of different environmental features but the IBR modelling suggests that urban areas play a major role, as observed for other rodent species (Mapelli et al. 2020;Munshi-South 2012).
Interestingly, our results show that the predicted difference in genetic structure (less structured generalist vs. more structured specialist) was maintained despite habitat alterations. We hypothesize a sort of compensatory effect, similarly acting in both species, between a decrease of gene flow in urban areas, while cultivated terrains increased population connectivity. Compensatory dynamics in disturbed habitats have been observed before in minks (Neovison vison) in Scotland and in the western slimy salamander (Plethodon albagula) in the USA (Oliver et al. 2016;Peterman et al. 2014). For salamanders, a compensatory dispersal behavior was observed through suboptimal habitats (i.e., high resistance surfaces) (Peterman et al. 2014). For our species, we observed higher gene-flow than predicted in some patches, which may be due to a compensatory behavior (increased gene-flow) through the landscape so previous patterns of population connectivity are maintained. This pattern, combined with compensatory landscape modifications (reduced connectivity in urban areas and increased in cultivated terrains) may have contributed to the maintenance of previous relative differences in the two investigated species. Further genetic monitoring of these populations can elucidate if this pattern will be maintained in time.
Our analysis highlighted the importance of woodland patches acting as a corridor between different sampled locations in both studied species, with particular emphasis for the bank vole. Urban areas were shown to be the major barriers experienced by gene-flow in the bank vole and wood mouse, albeit to a different degree. Moreover, railways and roads (never wider than 10 m in this area) cannot be considered as barriers to the dispersal of these species, consistently with previous studies Redeker et al. 2006). Indeed, roads appear as a factor that favors gene flow in the wood mouse. This may be because, for this species, the size of the roads present in the study area should not be considered as a barrier and/or that roads, in the environmental matrix, were included in (or surrounded by) a suitable ground. Similarly, cultivated fields do not limit dispersal, but may even play a role as corridors (Johnson and Munshi-South 2017;Miles et al. 2019;Tattersall et al. 2001). The only anthropogenic factor that seems to negatively affect the dispersal pattern is the presence of urban areas. Clearly, if woodlands will be further reduced by urbanization, genetic drift due to fragmentation could increase in our two studied species.
Spatial scale effects may be of particular importance in landscape genetic studies, especially on species having dispersal abilities limited. Different spatial scale could have revealed other environmental factors, but our sampling was spatially designed to test only a few factors related to anthropic activities.

Conclusions and implications for conservation
Overall, our results show that despite extensive habitat changes due to human activities, levels of genetic variation are quite high in both species, and their difference in the dispersal abilities is still reflected in the difference of genetic structure. The wood mouse, a generalist species with high dispersal ability, shows in fact higher genetic connectivity than the bank vole, which is a less mobile species closely linked to woodland areas. Nevertheless, we found also that cultivated fields and urban areas modifies the natural dispersion patterns in both species, probably in a way that will, in the future, increase the difference between their genetic structure. Our study supports the view that patterns of gene flow can be similarly affected, even in sympatric species with different ecological characteristics, by the same changes of land use. Locally, this implies that future monitoring efforts should prioritize the bank vole, the species with the highest genetic structure where genetic fragmentation is more likely to increase due to urbanization.