Conservation genetics of a wide-ranged temperate snake: same species, different locations, and different behaviour

Even though reptiles are threatened worldwide, few studies address their conservation, especially snakes. The goal of our study was to measure the genetic structure of a widely distributed temperate reptile, the smooth snake Coronella austriaca using eight microsatellite markers in two different areas at the core (Alsace, north-eastern France) and at the edge (Wallonia, southern Belgium) of its range. We sampled 506 individuals in 38 localities (respectively 10 and 28). Analysis of genetic structure conducted with a clustering method detected three clusters in Alsace, one group gathering all populations but two. In Wallonia, differentiation was observed on both sides of the Meuse river and in the Southern Ardenne region (southernmost sampling sites). Spatial autocorrelation analysis showed that individuals share parental relationship up to a distance of 2.8 km in Alsace and up to 10 km in Wallonia. Isolation by distance was detected in Wallonia but the distance explained a very limited part of the differentiation (r = 0.033), whereas no isolation-by-distance pattern was detected in Alsace. Even though genetic differentiation between populations separated by large rivers, highways, or crop elds was detected, dispersal between populations seem currently sucient to avoid any kind of genetic drift in both regions. These results are strongly contrasting with a previous study in England, suggesting sharp local variation of genetic structuring and diversication between location within the same species, probably related to the position in the distribution area and different densities.


Introduction
In the scope of the current biodiversity crisis, most studies on animal conservation focus on endangered and highly threatened species (Foden et Reading and Jofré 2020). Thus, it is crucial to understand the evolutionary mechanisms involved in such adaptations. According to the IUCN Red List, one of the major threats to terrestrial wildlife is the destruction of habitats (Pimm et al. 2014; IUCN 2020), which induces fragmentation and loss of connexion between populations. Small vertebrates such as reptiles are particularly sensitive to such disturbances, and it has been shown that some populations of common species are particularly exposed to fragmentation, and thus can become severely threatened (Driscoll 2004;Guiller 2009).
As amphibians, non-avian reptiles face declines at a global scale (Gibbons et al. 2000). For example, in Europe, about 20% of species are listed in the threatened categories of the IUCN (Cox and Temple 2009).
All the more, several species of snakes are facing decline at a global scale, among which some are found in Europe (Reading et al. 2010). Fragmentation, alteration, and destruction of habitats are highlighted as the main threats to reptile populations in Europe (Corbett 1989;Cox and Temple 2009;Reading et al. 2010). These processes lead to isolation of populations and rupture in gene ow (Frankham et al. 2002;Dixo et al. 2009), and are often associated with a reduction of population size, genetic drift, and loss of genetic diversity (Frankham et al. 2002). Therefore, adaptation capacities of populations faced with biotic or abiotic events such as diseases, unusual climatic events, or pollution, become reduced and thus can speed up the decline process at a local level (Frankham et al. 2002). The cumulative effects of these factors can lead to extinctions of populations and, in the case of species with a rather limited range, to the extinction of species. In vertebrates, such rapid extinctions have been observed, for example within amphibians (Pounds and Crump 1994;Stuart et al. 2008;Collins and Crump 2009). Few studies have been carried out in order to assess the conservation status and the genetic structure of populations of snake species (Mullin and Seigel 2009). Still, snakes can be indicators of the level of fragmentation in a given landscape (Guiller and Legentilhomme 2006). Moreover, a limited number of studies were conducted on the same species but in different locations, in order to determine if species can react differently in different habitats or if their dispersal behaviour is xed (Lane and Shine 2011). Consequently, evaluating the genetic structure within different regions in the same species could provide valuable information on the biology of the species, but also on how species behave in different environments.
Snakes are particularly sensitive to habitat fragmentation and direct destruction (Bonnet et al. 1999). Additionally, reptiles, and particularly snakes, are not often used as model species in conservation biology (Bonnet et al. 2002;Mullin and Seigel 2009). Therefore, we decided to examine a cryptic and rather poorly studied species given its wide range, the smooth snake, Coronella austriaca. Its distribution is rather well de ned (Sillero et al. 2014), and a recent study suggest that the widely distributed taxon currently assigned to Coronella austriaca throughout Europe, western Asia and the Middle East might in fact be a In this context, our study aimed at evaluating through genetic structure the magnitude of fragmentation on Coronella austriaca in two mainland environments with similar population densities but in different climates and locations within the core of the range of the species (Alsace, north-eastern France), and at the northern edge of the range (Wallonia, southern Belgium). The main questions we tackled were to assess the genetic relations of different population patches distributed throughout Alsace and Wallonia using microsatellite markers and evaluate if, as in England, the genetic isolation is clear between proximate populations. We hypothesised that as the smooth snake seems to be a poor disperser (Völkl and Käsewieter 2003;Pernetta et al. 2011;Dick and Mebert 2017), genetic structure should be marked even between geographically close populations. Also, we expect that even if the habitats and climate are somewhat different, the general pattern of genetic differentiation should be quite similar between regions.

Study area
The study area comprised of two densely inhabited (roughly 220 inhabitants per km 2 ) regions within the range of the Western 1 Coronella austriaca clade described by Jablonski et al. (2019): Alsace and Wallonia (Fig. 1). The sampling area in Alsace (north-eastern France) was located between 47°44'N and 48°3'N and 7°E and 7°27'E, at elevations ranging from 200 to 500 m a.s.l (Fig. 1A). The sampling area in Wallonia (southern part of Belgium) covered a territory of about 10,000 km 2 . The latitude was comprised between 49°32' and 50°38'N, and the longitude between 4°26' and 5°56'E ( Fig. 1B). Altitudes of the studied habitats were comprised between 70 and 400 m a.s.l. Both regions contain several environmental elements that could induce habitat fragmentation, such as networks of motorways, roads, railways, wide surfaces of crop elds and vineyards (Alsace), and large surface of coniferous forest culture (Wallonia).

Sampling design
We km, and the two most distant ones were separated by 85 km. Snakes were found either by sight or with the use of arti cial shelters (70cm x 70cm dark rubber plates in Alsace and of various materials in Wallonia). In order to avoid double sampling, a photograph of the pileus (dorsal side of the head) of each specimen captured enabled individual recognition on the eld (Sauer 1994(Sauer , 1997 Forward dyed primers were used in order to analyse them with an automatic sequencer (AB3130xl Applied Biosystem). Allele lengths were then read with the software PEAK SCANNER v.  (Goudet 1995). Genetic comparison between sampling sites were conducted with an ANOVA (with and without locus or population as factor) in R v3.6.3 (R Development Core Team 2016). We also performed a hierarchical structural analysis of gene diversity (AMOVA) to assess the molecular variance among sampling sites, among individuals and within individuals with GenAlEx, again for each region separately. We evaluated population subdivision for both region with a Bayesian clustering approach implemented in the software GENELAND v4.0.3 (Guillot et al. 2005). This method can be used to infer the number of genetic clusters (K) from the individual genotype distributed dataset in a spatial framework. We rst performed ve independent MCMC runs with K ranging from 1 to 10 for the Alsace region and from 1 to 30 for Wallonia, with the following parameters: 500 000 MCMC iterations, 5 000 thinning, maximum rate of Poisson process xed at 100, uncertainty attached to spatial coordinates xed to 0.2 km. Then, we ran the MCMC model 100 times with the same parameters, ve times rst to determine the best K value and the sixth simulation was conducted with the best K value only as suggested by Guillot et al. (2005). From the last simulations, we selected the 10 runs with the highest mean logarithm value of posterior probability, and calculated the posterior probability of population membership for each pixel of the spatial domain for each of these 10 runs, using a burnin of 10%. The number of pixels was set to 100 for the X axis and 350 for the Y axis for the Alsace region and respectively 350 and 350 for Wallonia, in order to avoid having two sampling sites in the same pixel. Finally, we computed posterior probability of population membership for each pixel of the spatial domain and the modal population of each individual. We then ran again standard population genetic analysis based on the number of populations inferred with GENELAND by calculating pairwise F ST and F IS with FSTAT, and tested the signi cance of the inferred structure by performing an AMOVA with GenAlEX.
For each region, we tested isolation of sampling sites by distance (IBD) with a Mantel test (Mantel 1967) by confronting corrected genetic differentiation [F ST /(1-F ST )] with the log values of the geographic distances between each sampling site (Rousset 1997). This test was implemented in R with the mantel.rtest function from the ade4 package (Dray and Dufour 2007) and 10 000 repetitions. We combined the results of the IBD for both regions and also from the UK with the data included in the article of Pernetta et al. (2011). In addition, we performed a spatial autocorrelation analysis separately for both regions using SPAGeDI v.1.3 (Hardy and Vekemans 2002) in order to determine correlation between geographic distances and genetic relatedness measured by Moran's I-statistic (Moran 1950; Sokal and Wartemberg 1983). We assigned geographic coordinates for each locality. Distance classes were chosen in order to provide similar numbers of pairwise comparisons for each class, separately for each region.

Null alleles and genotypic disequilibrium
We detected the presence of a probable null allele in Alsace for Ca30, due to an excess of homozygosity. Therefore, Ca30 was discarded from our following analyses that were consequently conducted with eight microsatellites markers (Ca16, Ca19, Ca40, Ca43, Ca45, Ca612, Ca63 and Ca66).
We detected an excess of homozygosity in Wallonia in Ca40 and Ca45, suggesting the presence of null alleles for these two markers. Moreover, genetic disequilibrium was signi cant between Ca43 and Ca612, as well as between Ca43 and Ca16. We consequently used only six microsatellite markers (Ca16, Ca19, Ca30, Ca612, Ca63 and Ca66) for this dataset.

Genetic variation and diversity
The number of alleles varied from 2 (Ca63) to 17 (Ca40) in Alsace (Table 1) and from 3 (Ca63) and 19 (Ca34) in Wallonia (Table 2). Globally, the allelic richness seems a bit higher in Alsace, even if four markers were different between both dataset (Ca30, Ca40, Ca43 and Ca45). The expected heterozygosity was also slightly higher in the Alsace region, even if the observed heterozygosity is similar. However, A R , H O or H E were not signi cantly different between regions when considering the ve similar loci in both populations, even when considering locus or populations as cofactors (for all P > 0.189).

Isolation by distance and spatial autocorrelation
Mantel's correlation test did not reveal any effect of IBD between the studied sites (Mantel test: r = -0.023; P = 0.86; Fig. 2) in the Alsace region, but the distance was negatively correlated with the distance at the margin of statistical signi cance when tested within the rst group obtained with GENELAND (Mantel test: r = 0.079; p = 0.080). In Wallonia, IBD was detected at the whole scale (Mantel test: r = 0.032; P = 0.0003; Fig. 2), but within the groups identi ed with GENEALEX, a signi cant signal of IBD was only present in the group5 (Mantel test: group4 r = -0.007; P = 0.730; goup5: r = 0.389; P = 0.032; group6: r = 0.121; P = 0.173).
The spatial autocorrelation between Euclidian distance and relatedness (measured as Moran's index) was signi cant for the distance classes between 0 and 2.8 km for the Alsace region. This indicated that smooth snakes in Alsace are more related to each other within a distance of 2.8 km. For Wallonia, this distance is even bigger, as a signi cant autocorrelation was detected up to a distance of 10 km.

Genetic differentiation in Alsace and Wallonia
Globally, the genetic differentiation between populations in both Alsace and Wallonia is limited or follow natural, geographical or historical isolation (River Meuse in Wallonia, or the southern part of the Ardennes). Indeed, as no or only weak isolation by distance was detected (opposite to the results in UK; Pernetta et al., 2011; see below), it is likely that other processes are in uencing the genetic differentiation between populations, such as some landscape elements (e.g., rivers) or historical events (e.g., climate uctuations) that we could not detect with this study.
In some speci c cases, populations are more differentiated than expected based on the genetic pattern observed in both regions. In Alsace, the south-westernmost group (cluster 3, Fig. 1A) seems to be isolated from the others, according to GENELAND and presenting a signi cant pairwise F ST with all other sampling. Also, it has a negative F IS value (-0.048), which, though not signi cant, indicates a propensity to outbreeding. This result could be caused by a reduction of the size of this population, due to a recent isolation event, and thus conduction to some genetic drift. Indeed, this sampling site is located west and north of two main highways. These elements might represent recent physical barriers to gene ow.
Similarly, other sub-populations included in cluster 1 separated by major highways or by large areas of crop elds also show high and signi cant F ST values (e.g., highways between pop4 and pop8; crop elds between pop4 and pop5; see Table 2). Therefore, it would then be possible that populations of C. austriaca located at the southwest of Alsace are isolated from the rest due to the fragmentation caused by highways. It has been shown that average sized and small species of snakes tend to avoid crossing roads (Andrews and Gibbons 2005) or are killed when trying to cross roads (Bonnet et al. 1999). Our results tend to suggest that it should be the case for C. austriaca, as highways could indeed constitute a strong barrier to gene ow if no underpasses are found along large sections of such roads.
In Wallonia, a signi cant but weak (r = 0.033) signal of isolation by distance was detected. Indeed, it could be related to the global structure detected with GENELAND, with the occurrence of four clusters, three of them representing well-separated regions [southern region (cluster 6); central region (cluster 4) and the edge of the central region (cluster 5)]. Even if some signal of IBD could be detected within clusters 5 and 6, it was signi cant only for cluster 5. Moreover, the grouping realised with GENELAND is not based on distances as sampling sites of clusters 4 and 5 are sometimes very close. This splitting is probably more related to historical reasons. For instance, the group6 gathered all the populations from the southern part of the Ardennes, where the species is not very common as habitats are cold and mainly composed of forests. We can hypothesise the differentiation between clusters 4 and 5 is resulting from the occurrence of the Meuse River. Indeed, all but one (sampling site 28) sampling sites from the group5 is on the north-western part of the Meuse River. Within group4, all sampling site except sampling site 3 are on the shore or south-eastern part of this river. The Meuse River is the largest river in Belgium; it probably acted as a barrier to the movement of C. austriaca for several centuries. Cluster 7, that gathered individuals from the single sampling site 21, does not present particular geographic barriers with other populations of the group5 that could explain its genetic differentiation. Local monitoring in the sampling site 21 highlighted a strong increase of individuals during the last years, with the lack of smooth snake 30 years ago (Graitson et al. 2012). We can consequently hypothesise that this population undergone a strong founder effect with the colonisation by a very limited number of individuals only 2-3 generations ago, which could explain the signi cant F IS value. Contrary to what is observed in Alsace, the difference between the groups does not seem to be explained by a barrier effect induced by motorways, which are nevertheless present in the sampled area, but more by geographical elements. We believe that the hilly terrain in Wallonia offers more possibilities to cross the motorways through underpasses.
Such a low genetic differentiation was unexpected, as strong differentiation was detected in several species of snakes with similar ecology requirements, even within putative interconnected habitat: for example, ecologically interconnected populations of Nerodia sipedon, an aquatic colubrid from North America, which has a similar home range as C. austriaca (between 1 and 4 ha), showed a marked genetic differentiation, maybe resulting from a high degree of philopatry (Prosser et al. 1999 over an even larger distance (maximum distance between sampled sites ≈ 125 km) and with the occurrence of at least four genetic groups. It is to note that the F ST values between the three studies did not result from the exact same set of genetic markers (2/8 were similar for the three studies; 3/8 between Alsace or Wallonia and England), but the genetic diversity, number of alleles and allelic richness are similar between all loci, suggesting that the use of different markers would have only a limited impact on the comparison. Also, the effect of isolation by distance was signi cant in Southern England (r=0.511, p<0.05), whereas no effect was detected in Alsace and only a weak signi cant signal in Wallonia ( g. 2).
We expected to nd a stronger effect of isolation in both Alsace and Wallonia due to the larger distances between populations if a similar genetic pattern as in England has been detected, which was not the case. This observed discrepancy obtained at a different scale should lead to further studies at the same spatial scale and with the same set of microsatellites in order to avoid artefacts due to large variation in distance between populations. However, the comparison of Pernetta et al. (2011) and this project clearly suggests that, within a species, genetic structure can strongly vary between habitats or regions. Such Therefore, it would be interesting to investigate the genetic structure of populations of Coronella austriaca in other parts of its distribution limits and in similar habitats (i.e., lowlands), in Scandinavia or in Western France for example, to detect if the dispersal behaviour varies in function of the positions within the distribution of the species or more due to local geographical elements.
Moreover, the studied populations in both regions are still large enough and rather widespread to avoid a strong genetic drift, as shown by the similar level of genetic diversity and limited F ST values (Table 1 and 2), contrary to the populations studied in England (Pernetta et al. 2011