Characterising the relationship between suitable habitat and gene flow for Myotis bechsteinii and Eptesicus serotinus in Britain

Habitat suitability models (HSM) have been used to understand the impacts of landscape-scale habitat connectivity and gene flow mostly by assuming a regular decrease in the cost of movement as habitat improves. Yet, habitat selection and gene flow are governed by different behavioural processes which may limit the reliability of this approach as individuals are likely to disperse through unsuitable habitat for breeding. The aim of this study was to identify the optimal relationship between gene flow and HSMs for two bat species (Myotis bechsteinii and Eptesicus serotinus) in Britain by testing a range of nonlinear negative exponential functions for the transformation of HSMs into resistance surfaces. We modelled habitat suitability using a hierarchical, multi-level approach that integrates models across three nested levels. Then, we measured the relationship between published genetics data of both species and six negative exponential transformations of the predicted outputs. The two most extreme transformations provided the best fit to genetic data for both M. bechsteinii (c = 32; R2 = 0.87) and E. serotinus (c = 16; R2 = 0.42). The negative linear transformations had the poorest fit. These results suggest that bats are able to disperse through areas of poor habitat for breeding, but will avoid the most unsuitable areas. We recommend comparing multiple transformations of HSMs at different resolutions to gain a more accurate representation of gene flow across heterogeneous landscapes and to inform cost-effective, targeted management.


Introduction
Land use change and intensification are major contributors towards the rapid and ongoing decline of biodiversity (Ceballos et al. 2017). Their negative effects can be amplified when coupled with climatic and other environmental change leading to population declines (De Chazal and Rounsevell 2009), the loss of genetic diversity, and ultimately species loss (Bailey 2007). Understanding how landscape composition and structure influence the genetic connectivity of a species is essential for wildlife conservation. This is particularly true for rare or elusive species for which we have limited knowledge, and species more at risk of genetic drift and local extinctions due to traits such as narrow niche requirements or low dispersal abilities (e.g. Razgour et al. 2011). Habitat Suitability Modelling (HSM) is increasingly used to fill these knowledge gaps Bellamy et al. 2013). This approach predicts the distribution of a species and provides information on the underlying drivers using occurrence records and environmental data layers (Guisan et al. 2017;Guisan and Zimmermann 2000). As well as understanding and predicting habitat suitability, effective conservation action requires information on how landscapes influence species movements and gene flow (Razgour et al. 2016). Although major barriers, such as mountain ranges or seas, are easy to identify and detect for some species (e.g. Garroway et al. 2011;McLeod et al. 2015;Moussy et al. 2015), the heterogeneous structure and composition of many landscapes will affect population processes in a more complex way. Therefore, detailed analytical methods are needed to help understand and predict these relationships. Landscape genetic approaches have been used to assess the amount of connectivity between habitats in a fragmented landscape (e.g. Clauzel et al. 2015;Ramirez-Reyes et al. 2016); characterise features affecting genetic variation; identify isolated populations; map and design potential wildlife corridors to facilitate dispersal and gene flow (Balkenhol et al. 2016;Razgour et al. 2016) and predict the impact of landscape changes on wildlife (Razgour et al. 2017).
Predicted habitat suitability often correlates with gene flow (e.g. Razgour et al. 2017). Most studies have used negative linear transformations (resistance decreases at a constant rate as suitability increases) of habitat suitability models into resistance surfaces (the degree to which movement of an organism is affected in a grid cell; Spear et al. 2010) to assess the relationship between gene flow and habitat suitability (Pullinger and Johnson 2010;Richards-Zawacki 2009). However, Beier et al. (2009) challenged this assumption as habitat selection and dispersal/gene flow are governed by different behavioural processes. The first is mainly driven by adult home ranges, whereas gene flow is primarily influenced by the dispersal or movement of individuals and will vary depending on the species' dispersal abilities (Balkenhol et al. 2016). It has been argued that dispersing individuals will be more tolerant to 'sub-optimal' habitat than individuals foraging within their home ranges. Therefore, the creation of resistance surfaces, where the species' ability to move decreases slightly as habitat deteriorates with the exception of lower suitability values where movement is almost impossible, is likely to better reflect dispersal and gene flow. Several studies applying negative exponential functions to transform HSMs into resistance surfaces ( Fig. 1) to test this hypothesis have proven successful for different species (Trainor et al. 2013;Keeley et al. 2016;Zeller et al. 2018). Yet, these are still not commonly applied when deriving HSMs into resistance surfaces.
Bats, being small and volant, are highly mobile and able to travel long distances within a short period of time in comparison to other terrestrial mammals (Ahlén et al. 2009;Bryja et al. 2009). Behaviours, such as autumnal swarming, where individual bats travel significant distances to aggregate at underground sites are also known to be linked with mating and facilitate gene flow (Kerth et al. 2003;Rivers et al. 2005). Such behaviours can be observed in highly specialist and sedentary species, such as Myotis bechsteinii (Dekeukeleire et al. 2016). Therefore, the transformation of HSMs into resistance surfaces by applying a negative exponential function, where only very poor habitat will act as a significant barrier for movement, as opposed to the more commonly negative linear transformation is likely to better predict landscape-scale genetic connectivity for most bat species.
The aim of this study is to identify the optimal relationship between gene flow and HSMs for two bat species exhibiting different reproductive and foraging strategies in Britain by testing a range of nonlinear negative exponential functions. We tested the functions on Myotis bechsteinii, a woodland specialist that performs autumnal swarming and Eptesicus serotinus, a more generalist species that can forage and roost in built-up areas but is not known to perform autumnal swarming. We hypothesised that higher transformations of HSMs would better reflect the relationship between suitable habitat and gene flow.

Habitat suitability models
We used the hierarchical, multi-level HSM framework developed by Bellamy et al. (2020) to model habitat suitability in Britain for M. bechsteinii, a woodland specialist that performs autumnal swarming, and E. serotinus, a more generalist species that is not known to occur at swarming sites. This method nests models across three levels (population range, home range, local habitat -See below) to limit predictor collinearity and enable context dependency. MaxEnt (Phillips et al. 2006), a presence-only HSM approach, was used to predict suitability at each level. All analyses were carried out using R (v. 3.4.2; R Core Team 2018) in R Studio (v.1.1.383; RStudio Team).
We used landscape variables reflecting environmental features and conditions likely to influence the distribution of temperate bats (See supplementary 1 for predictor variables; Boughey et al. 2011;Bellamy et al. 2013;Fuller et al. 2018;Rebelo et al. 2010). We performed an initial population range model focusing on climate at a 5 km resolution. This was followed by a summer home range (1 km resolution), and a local habitat (100 m resolution) model. The local habitat model was restricted geographically to southern Britain (i.e. were within both species' known current range), and employed land cover, landscape structure and land management variables.
We collected presence records between 2000 and 2016 during the months of April to September at a B 10 km resolution which were extracted from the National Biodiversity Network (NBN; https:// nbnatlas.org/) gateway, national and local bat monitoring schemes, and local record centres. For M. bechsteinii, we only used confirmed roosts and bats identified in hand, because the species' echolocation calls overlaps with that of other Myotis species making the use of acoustic data unreliable. The lack of information in the reporting of E. serotinus presence data throughout Britain meant that we also included confirmed acoustic records, as the species can be differentiated from others, along with roosts and bats identified in the hand. We filtered records to match the resolution of each model (i.e. 5 km, 1 km, and 100 m respectively) and to retain a single record per grid square. To reduce the amount of sampling bias towards heavily sampled areas (Graham et al. 2015), we used 10,000 pseudoabsence records of all bat species within the study extent from the same sources recorded over the same time period (Phillips et al. 2009). It should be noted that this approach does not entirely control sampling bias for M. bechsteinii as most bat records originate from acoustic data and M. bechsteinii can only be reliably identified in hand. The population range model was run first, and the output Habitat Suitability Index (HSI) values from each of the three levels were disaggregated to the resolution of the next lower level and incorporated as an explanatory variable. For each level, we removed highly correlated variables using the 'vifstep' stepwise function of the 'usdm' package (Naimi et al. 2014) and a conservative VIF threshold of three (Zuur et al. 2010). The package 'ENMeval' (Muscarella et al. 2014) was used to identify the optimal MaxEnt model settings. We tested combinations of feature types (L, linear; H, hinge; Q; quadratic; P, product) and disabled threshold features to reduce overfitting. We varied the regularisation multiplier in steps of 0.5, from 0.5 to 4. A final model was then created using the optimal settings to produce model predictions.

Landscape genetics
Myotis bechsteinii and Eptesicus serotinus genotypes were taken from previously published papers (Moussy et al. 2015;Wright et al. 2018), but excluded samples collected outside Britain. The M. bechsteinii dataset included genotypes from 260 individuals using 14 microsatellite loci collected at 8 sites. The E. serotinus dataset comprised 593 samples using 10 loci from 28 sites. We used the genetic landscape-shape analysis from the Alleles in Space software (Miller 2005) to identify patterns of genetic structure throughout the study areas (Moussy et al. 2015;Wright et al. 2018). The connectivity networks produced from this software, were based on Delaunay triangulations (Watson 1992;Brouns et al. 2003) and were used to build a three-dimensional map, where the Z-axis corresponded to genetic distance between populations.
To transform HSI values into resistance surfaces, we used a negative linear transformation and six transformation values on a negative exponential function ( Fig. 1) as reported by Trainor et al. (2013), Mateo-Sánchez et al. (2015) and Keeley et al. (2016): H = habitat suitability value, c = 1, 2, 4, 8, 16 and 32. Finally, because woodland connectivity is frequently used as an indicator of habitat connectivity for bats, we also created a woodland cover resistance surface where the highest resistance values (100) were associated with the lowest amount of woodland.
We used Circuitscape v4.0.5 (McRae and Nürnberger 2006) to calculate resistance distance matrices between study populations (focal nodes) using the ''pairwise'' modelling mode. To test the relationship between the transformed HSM variables and gene flow for M. bechsteinii and E. serotinus in Britain, we performed multiple regressions on distance matrices (MRDMs) using 10,000 permutations (Legendre et al. 1994;Goslee and Urban 2007). We performed MRDMs on the output matrices from Circuitscape, along with genetic and geographic distance matrices, as described by Santos (2017). Distance matrices were initially standardized from 0 to 1 using the 'scales' package (Wickham 2014). The effect of geographic distance was then removed by extracting the residuals of linear models of geographic distance against each resistance surface (Dyer et al. 2010). Univariate MRDMs of genetic distance were finally calculated for each resistance surface, in order to assess the correlation of different HSM transformations with gene flow. MRDMs were performed for all transformations for which the sea was assigned a resistance value of 70 and 120 to identify the optimal effect of the sea separating the Isle of Wight (an island on the coast of southern England) from mainland Britain on gene flow. For both species, suitable conditions predicted from the population range model were concentrated in southern Britain and were mainly influenced by temperature, precipitation and decreased woodland cover (Supplementary Information 2-8 for permutation importance and 3-9 for response curves). At the home range and local levels, disaggregated HSI values from preceding levels were the strongest explanatory variables with permutation importance's ranging from 76 to 86%. In the home range level model, high cover of both protected broadleaf and ancient woodland were associated with the presence of M. bechsteinii; whilst no clear patterns were identified with E. serotinus outside the importance of the preceding population range level (Supplementary Information 4-10 for permutation importance and 5-11 for response curves). Woodlands were also associated with M. bechsteinii at the local level as the species was primarily influenced by woodland patch shape and broadleaved woodland cover. Finally, E. serotinus habitat showed strong associations to areas close to buildings and hedgerows within a 100 m scale ( Fig. 2a and b; Supplementary Information 6-12 for permutation importance and 7-13 for response curves).

Landscape genetics
The genetic landscape-shape analysis revealed contrasting population structure between M. bechsteinii and E. serotinus. While only two M. bechsteinii populations showed signs of isolation, many E. serotinus populations were estimated to be isolated (Fig. 3).
For M. bechsteinii, the most extreme transformation of the local HSM output (c = 32) showed the highest correlation coefficient (R 2 = 0.87, p \ 0.001), while the second highest transformation had the   Table 2). The correlation coefficients did not all increase with increased transformation levels. In most models, coefficients were highest at c = 16 and, in the case of M. bechsteinii, these notably declined at c = 32 in the population and home range models ( Table 2). The models also revealed that woodland cover did not correlate with gene flow for both species (R 2 \ 0.05 and p [ 0.5 for both models).

Discussion
Maintaining gene flow between populations is crucial for securing a species' viability. However, understanding how small, mobile species, such as bats, move through the landscape in order to maintain gene flow is hard to measure directly. In this study, we assessed the strength of association between a range of nonlinear negative exponential functions derived from HSMs to describe the relationship between predicted suitable habitat and gene flow of two bat species (Myotis bechsteinii and Eptesicus serotinus) in Britain. We found that the two highest transformations (c = 32 and c = 16) of HSMs into resistance surfaces had the strongest relationship with gene flow for both species. Negative linear transformations of HSMs into resistance surfaces are frequently used for assessing the degree of permeability of different habitats to species, yet their sole use would have underestimated the ability for both species to maintain genetic connectivity in a fragmented landscape. Our results suggest that the design of corridors and other management efforts to maintain gene flow between populations of both species require a high negative exponential transformation of habitat suitability models. These findings also mirrored the few studies having previously assessed the use of such transformations on other species (Trainor et al. 2013;Mateo-Sánchez et al. 2015;Keeley et al. 2016;Zeller et al. 2018). These results also have important implications for bat conservation as they could provide new insight into the ability of populations to shift their distribution under climate change (Razgour et al. 2017).
The hierarchical framework has been found to improve the accuracy of HSMs by integrating predictors at their scales of effect (Bellamy et al. 2020). For landscape genetics, the framework also offers insight into the models that better represent gene flow, albeit we can only make assumptions on the reasons why certain models show higher correlations (e.g. higher resolution, higher AUC values, different environmental variables). Interestingly, the two highest transformations of the local level HSM model (highest resolution) were found to best explain gene flow. However, for M. bechsteinii, the correlations for the highest transformations for the population and home range models dropped significantly once the highest correlation was reached at transformation 16. This suggests that a wide range of transformations of different resolutions must be tested as correlations do not always increase along with HSM transformations.
The highest correlation level for M. bechsteinii was approximately twice as high as for E. serotinus. This could be explained by the higher performance of HSMs, which may be linked to the fact that the ecological niche of habitat specialists, such as M. bechsteinii, tends to be better predicted (Evangelista et al. 2008). Myotis bechsteinii was also associated with the highest transformation as opposed to the second highest for E. serotinus meaning that resistance values were higher throughout the landscape for E. serotinus. Some bat species, such as M. bechsteinii, are known to perform autumnal swarming which is linked with mating. While M. bechsteinii tend to forage over short distances in old growth woodlands (* 3 km, Schofield and Morris 2000), they will readily travel distances of up to 20 km to reach swarming sites (Dekeukeleire et al. 2016). In a fragmented landscape, as found in Britain, such movements would inevitably entail crossing some unsuitable habitat while the most unsuitable habitat will act as a barrier. The effective dispersal   Fig. 1 and the highest correlations appear in bold mechanism of M. bechsteinii has likely reduced the impact of habitat fragmentation. However, the isolation of certain populations raises questions on a possible dispersal threshold for the species. In other words, it is important to identify how much unsuitable habitat individuals are able to traverse before populations become genetically isolated.
As for E. serotinus, very little is known about their mating system and dispersal abilities. The species ranges over significantly larger areas (Catto et al. 1996) and will forage and roost in a broader range of habitats, including built-up areas, making it possibly less sensitive to human development than M. bechsteinii. However, it does not exhibit swarming behaviour like Myotis bats. Limited dispersal for mating could result in the greater genetic differentiation between populations observed in Britain as the species may be less likely to travel long distances and move through unsuitable habitat for mating. This could explain the stronger association with the second highest transformation (c = 16) as opposed to the highest for M. bechsteinii (c = 32). Differences between the dispersal abilities of both species can be observed on the Isle of Wight (UK) located 1-5 km from the south coast of England. Here, genetic structure between the Isle of Wight and Britain for E. serotinus has been identified (Moussy et al. 2015), while no apparent genetic structure was observed for M. bechsteinii (Wright et al. 2018). This reinforces the proposal that gene flow in M. bechsteinii is less restricted by the presence of unsuitable habitat than in the case of E. serotinus.
Whilst still being rarely performed the use of negative exponential transformations of HSMs tend to provide better insight into the dispersal mechanisms of wildlife (Trainor et al. 2013;Keeley et al. 2016;Zeller et al. 2018) and must continue to be encouraged. These may be particularly relevant for the design of wildlife corridors for species willing to travel through large areas of unsuitable habitat to establish new territories and breed. For example, terrestrial carnivores, such as Canis lupus or Puma concolor (Zeller et al. 2018), may disperse over hundreds of kilometres and will inevitably use unsuitable habitat. The added information from these models will help identify the most likely paths animals will use and where small improvements in habitat can be made to encourage the successful dispersal of a species as opposed to its establishment.

Conclusions
Our aim was to identify the optimal relationship between gene flow and HSMs for two bat species, M. bechsteinii and E. serotinus. For this we tested multiple negative exponential transformations of HSMs into resistance surfaces as opposed to the more commonly used single negative linear transformation. Our study highlights the importance of testing a multitude of transformations as the highest correlations for both species was identified from the two highest transformations (c = 16 and 32). These findings have important implications for bat conservation as the sole use of negative linear transformations would have underestimated the ability of both species to maintain gene flow between populations. Bats exhibit complex mating behaviour and may travel significant distances outside of high quality habitat for mating. It is therefore important that such transformations are applied on HSMs for landscape genetic studies as they may provide a more accurate representation of dispersal in a heterogeneous landscape.
Data availability Data used for this study was obtained from previously published resources (Moussy et al. 2015;Wright et al. 2018). Code and map outputs have been made available on figshare (https://doi.org/10.6084/m9.figshare.15036192).

Declarations
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.