Research Article
Print
Research Article
Genetic identification, morphology and distribution of Natrix helvetica subspecies in southern and western Switzerland (Reptilia, Squamata, Serpentes)
expand article infoAndreas Schild, Hannes Baur§, Stefan T. Hertwig§, Uwe Fritz|, Sylvain Ursenbacher#¤
‡ University of Bern, Bern, Switzerland
§ Natural History Museum Bern, Bern, Switzerland
| Museum of Zoology (Museum für Tierkunde), Dresden, Germany
¶ University of Basel, Basel, Switzerland
# University of Neuchâtel, Neuchâtel, Switzerland
¤ Balaton Limnological Research Institute, Tihany, Hungary
Open Access

Abstract

Most of Switzerland is inhabited by the nominotypical subspecies of the barred grass snake (Natrix helvetica helvetica), which is characterized by mitochondrial DNA lineage E. Only in the northeast of the country, the common grass snake (N. natrix) occurs and hybridizes with N. h. helvetica in a narrow contact zone. However, we discovered that in southern and western Switzerland barred grass snakes representing another mtDNA lineage (lineage C) are widely distributed. Lineage C is typical for Alpine populations of the southern subspecies N. h. sicula. Our microsatellite analyses of the Swiss samples revealed differences between the two subspecies and also a substructure with two clusters in each subspecies. Furthermore, we discovered a contact and hybrid zone of N. h. helvetica and N. h. sicula along the northern shore of Lake Geneva and also confirm that interbreeding with alien common grass snakes (N. n. moreotica, mtDNA lineage 7) occurs there. This finding is of concern for nature conservation and measures should be taken to prevent further genetic pollution. Using morphometrics, we found no differences between the two subspecies of N. helvetica, while N. natrix was slightly distinct from N. helvetica.

Key words

Microsatellites, mitochondrial DNA, morphometrics, nuclear DNA, taxonomy

Introduction

Grass snakes constitute a complex of three species which were regarded as conspecific for many decades (Kabisch 1999; Pokrant et al. 2016; Kindler et al. 2017). Natrix astreptophora (Seoane, 1884) occurs in the North African Maghreb region, the Iberian Peninsula and adjacent France. Natrix helvetica (Lacepède, 1789) lives in Western Europe north of the Pyrenees and in Britain and Italy. In Central Europe N. helvetica ranges eastward approximately to the Rhine region. Natrix natrix (Linnaeus, 1758) occupies the largest distribution range of the three species, from east of the Rhine region across Fennoscandia, the Balkan Peninsula and large parts of the Near East to Lake Baikal in Central Asia. Traditionally, many grass snake subspecies have been recognized, based on often somewhat fuzzy morphological features such as body proportions, coloration and size (see reviews in Kabisch 1999; Kindler et al. 2013; Fritz and Schmidtler 2020). Most of these subspecies are currently no longer recognized (Kindler et al. 2018a; Schultze et al. 2020; Asztalos et al. 2021b). According to the current view (Asztalos et al. 2021b), two subspecies of N. natrix occur in Central Europe, the nominotypical subspecies N. n. natrix and N. n. vulgaris Laurenti, 1768. However, the two subspecies hybridize in large parts of southern and southeastern Central Europe, while the populations in northeastern Switzerland represent N. n. vulgaris with introgressed mitochondria of N. n. natrix (Asztalos et al. 2021a, 2021b). For N. helvetica, four subspecies are currently distinguished. Natrix helvetica helvetica is distributed north of the Alps, while N. h. sicula (Cuvier, 1829) occurs south of the Alps, i.e., on the Italian Peninsula and on Sicily. Two further subspecies live in Corsica (N. h. corsa [Hecht, 1930]) and Sardinia (N. h. cetti Gené, 1839; Schultze et al. 2020). Natrix helvetica sicula crossed the Alps at least twice and occurs in the Inn River drainage of Austria and southernmost Bavaria (Glaw et al. 2019; Asztalos et al. 2021a). In addition, in Switzerland one record is known from beyond the Simplon Pass in the canton Valais (Kindler and Fritz 2018). Natrix helvetica sicula harbours several deeply divergent mtDNA lineages, reflecting ancient divergence processes on the Italian Peninsula and Sicily that began approximately 6.8 million years ago (Kindler et al. 2018a; Schultze et al. 2020). The populations of N. h. sicula relevant for the present study possess mtDNA lineage C (Schultze et al. 2020).

There are many regions across the distribution range of grass snakes where non-native individuals were introduced (France: Asztalos et al. 2020; Germany and Great Britain: Kindler et al. 2017; Italy: Schultze et al. 2020; Netherlands: van Riemsdijk et al. 2020; Asztalos et al. 2021c). Within the scope of a study in Switzerland, Dubey et al. (2017) genetically uncovered alien grass snakes (N. natrix) allegedly originating from Rijeka, western Croatia, which escaped in the 1970s from an outdoor reptile park in Lausanne (46.5427°, 6.6423°). However, the mtDNA lineage detected there (lineage 7 of Kindler et al. 2017, corresponding to N. n. moreotica [Bedriaga, 1882]) does not occur in the putative source region, but distinctly further southeast in the Balkans, western Anatolia and Cyprus (Asztalos et al. 2021b). Additionally, Dubey et al. (2017) found grass snakes yielding mtDNA lineage C at the same site in Lausanne. In contrast, Chèvre (2015) and Kindler and Fritz (2018) detected this lineage only in the cantons of Valais and Ticino, raising the question of whether these records also refer to non-native snakes or whether clade C, i.e., N. h. sicula, might have a wider distribution in Switzerland.

Dubey et al. (2017) assumed that grass snakes with lineage C also escaped from the reptile park and established in the region. However, it is also possible that N. h. sicula had passed the Alps in the Holocene and naturally occurs in western Switzerland. As grass snakes of clades C (N. h. sicula) and 7 (N. n. moreotica) have potentially hybridized at Lausanne with the native N. h. helvetica (clade E), it is of utmost importance for nature conservation to find out whether clade C is also native in western Switzerland or has been introduced.

The aim of the present study is to determine the natural distribution of N. h. sicula in Switzerland. To do so, we collected DNA samples from wild snakes and museum specimens to determine their mtDNA lineage. In addition, we used nuclear DNA markers (microsatellites) for Bayesian cluster analyses to estimate the amount of admixture between snakes corresponding to different mitochondrial lineages. In accordance with the concept of integrative taxonomy, we also used various characters of pholidosis and colour pattern as well as morphometric measurements to identify possible external morphological differences in our sample.

Materials and methods

Study area and sampling

The immediate study area was limited to southern and western Switzerland. It comprised the cantons of Vaud, Valais and Ticino (Fig. 1) as lineage C has been detected there before (database of info fauna – karch; https://www.infofauna.ch/). The sampling scheme was habitat-specific, i.e., we focused on preferred habitats of grass snakes, like ponds, rivers, lakes and spots with previous sightings. Snakes were caught by hand, measured on site, and released directly afterwards at the exact location of capture. DNA samples were taken with buccal swabs (dried and stored at -20 °C) or scale clipping (one to four ventral scales stored in 70% alcohol). Following Thorpe (1979) and Chèvre (2015), individuals with a total length of less than 50 cm were classified as juveniles. To reduce handling time, all necessary body parts of adult snakes were photographed for later morphological analyses. Furthermore, museum specimens from Vaud, Valais and Ticino were studied; a few additional individuals from the adjacent cantons were also included. From museum specimens, liver or muscle tissue was sampled. A list of samples is provided in Suppl. material 1: table S1.

Figure 1. 

Distribution of mitochondrial lineages of grass snakes in southern and western Switzerland (a) and microsatellite clusters according to our STRUCTURE analyses (b). Symbols and colours of mitochondrial lineages correspond to Kindler et al. (2013). Some symbols are slightly shifted from the sampling location, enhancing readability. The four samples in the bottom right corner originate from northern Italy beyond the map sector. Borders within Switzerland denote cantons. Abbreviations for cantons mentioned in the text: BE – Bern, GE – Geneva, SZ – Schwyz, TI – Ticino, UR – Uri, VD – Vaud, VS – Valais. Lake Geneva, mentioned in the text, is located in the southwesternmost part of Switzerland (cantons of Geneva and Vaud) and adjacent France. The maps were created by modifying a Wikimedia map (https://commons.wikimedia.org/wiki/File:Reliefkarte_Schweiz2.png).

DNA extraction and purification

Samples were incubated with ATL buffer and proteinase K (Qiagen) in a heat block for 16–20 h at 56 °C. Scales were previously placed in water for 24 h to remove alcohol. After digestion, liquid from swab tips was extracted using a centrifuge. The DNA was purified following the protocol “Purification of total DNA from Animal Tissues (Spin-Column Protocol)” of the DNeasy Blood and Tissue Kit (Qiagen) using a Qiagen robot.

mtDNA sequencing

To determine the mtDNA lineage, the cytochrome b gene (cyt b) and the NADH dehydrogenase subunit 4 gene with adjacent regions coding for tRNAs (ND4) were used, as in previous studies on grass snakes (e.g., Kindler et al. 2013, 2017; Chèvre 2015; Pokrant et al. 2016; Glaw et al. 2019; Schultze et al. 2020; Asztalos et al. 2021a, 2021b). Amplification followed Kindler et al. (2013) and Chèvre (2015). PCR products were sequenced by LGC Genomics GmbH (Berlin, Germany). Sequences were processed with CodonCode Aligner (https://www.codoncode.com) and compared with GenBank sequences to determine the mitochondrial lineage.

Microsatellite analyses

The same thirteen microsatellite loci were used as in Chèvre (2015) and Kindler et al. (2017) and genotyped according to their protocols. PCR products were analysed on an ABI 3130xl Genetic Analyser (Applied Biosystems) at the Zoological Institute of the University of Basel. The Microsatellite Plugin in GENEIOUS PRIME 2020.0.3 (https://www.geneious.com) was used to visualize peaks and determine allele lengths.

For inferring the nuclear genomic identity of the 73 successfully processed samples, the Bayesian clustering approach based on the Monte Carlo Markov chain (MCMC) algorithm implemented in the software STRUCTURE ver. 2.3.4 (Pritchard et al. 2000; Pritchard and Wen 2002) was used. STRUCTURE assumes unlinked microsatellite loci at linkage equilibrium and divides the dataset into partitions (K) optimized for the presence of Hardy-Weinberg equilibrium. After a burn-in of 100,000 generations, MCMCs were run for 200,000 iterations, ten times per K between one and ten. The optimal number of K was determined using both the ΔK method (Evanno et al. 2005) in STRUCTURE HARVESTER software (Earl and vonHoldt 2012) and the L(K) approach of Pritchard and Wen (2002). The best STRUCTURE run (highest likelihood) with the optimal K was used both to determine genotypic identity and to assess admixture. Snakes with an assignment ≥ 80% to a specific cluster were treated as pure. FST values between clusters were calculated with FSTAT ver. 2.9.3 (Goudet 1995) using only genotypically non-admixed individuals.

Morphological analyses

Only snakes exceeding 50 cm in total length and with a genotypic cluster assignment ≥ 80% were used for morphological examinations. Following Chèvre (2015), morphological variables with strong geographic variation were selected (Suppl. material 1: table S2) and analysed along with landmark data to examine for possible morphological differences between microsatellite clusters. The dataset was enriched with additional data from photographs provided by M. Chèvre (21 genotyped N. h. helvetica, mtDNA lineage E; 20 genotyped N. n. vulgaris, mtDNA lineage 3 from northeastern Switzerland), so that 31 males, 38 females and three sex-undetermined grass snakes were available for morphology.

Geometric morphometrics of landmark data

To obtain landmark data, standardized pictures were used showing the right and dorsal sides of the head of adult snakes. Fixing the focal length and manual focus of the camera ensured that the scale of the pictures was identical, which was double-checked using a ruler in the pictures. Photographs were taken twice to calculate mean landmark coordinates, which reduces potential inaccuracies due to slight shifts in photographing and placing landmarks. M. Chèvre provided only a single photograph per snake, for which landmark coordinates were produced twice to account for imprecise landmark placing. Mostly easily identifiable junctions of scales were chosen as landmarks to facilitate the workflow (Suppl. material 1: fig. S1).

Some landmarks were removed for analysis because the sample size was too small for 27×2 coordinate variables. Landmark 8 was removed because the temporal scale was sometimes divided and/or small, so it did not reach the 7th supralabial scale. Landmarks 12, 14, 16 and 18 were excluded as they all have other landmarks in close proximity. Lastly, landmarks 21, 24 and 27 were removed because they are located at the edge and might already be influenced by the curvature of the head. Therefore, only nineteen landmarks (1–7, 9–11, 13, 15, 17, 19, 20, 22, 23, 25, 26) were finally used (Suppl. material 1: figs S3, S4).

Landmarks were placed in the software TPSDIG2 ver. 2.30 (Rohlf 2017) and its coordinates were saved in tps files created by the software TPSUTIL32 ver. 1.74 (Rohlf 2013). In the statistical software R ver. 3.4.1 (R Core Team 2017), the function estimate.missing () from the package GEOMORPH (Adams et al. 2019) was used to interpolate missing landmarks with the thin-plate spline method. Then, mean values were calculated for each landmark per specimen and side.

All analyses of the mean landmark coordinates were performed in MORPHOJ ver. 1.07a (Klingenberg 2011), similar to the procedure described in Sidlauskas et al. (2011). First, a least-square Procrustes Fit was calculated. Allometric correction was then performed using a linear regression with the log centroid size as explanatory variable and the Procrustes coordinates as the response variable. The regression included a permutation test with 10,000 rounds and pooled regression within clusters. To examine shape changes, the regression residuals were included in a Canonical Variate Analysis (CVA) with a permutation test for pairwise distances of 10,000 iterations. Scatterplots of CV scores were checked for clustering of groups and wireframe graphs were used to visualize shape changes. The starting and target shape of wireframe graphs were placed next to each other, as suggested by Klingenberg (2013), to be able to objectively examine shape changes. Besides visually plotting the shape differences, the CVA also runs both Mahalanobis and Procrustes permutation tests (10,000 iterations) to check for the significance of shape differences.

Analysis of distance measurements

Distance measurements (SVL, TL, HL and HW; Suppl. material 1: table S2) were taken in the field and analysed using Multivariate Ratio Analysis (Baur and Leuenberger 2011; Baur 2024). The ‘shape PCA’ enables the examination of shape changes depending on size. The ‘PCA ratio spectrum’ then allows the interpretation of principal components (PCs) in terms of ratios and shows the most discriminating ratio with respect to a particular shape PC. This approach has been used in several studies to find morphological differences among taxa (László et al. 2013; Baur et al. 2014; Huber and Baur 2016; Gebiola et al. 2017; Waser et al. 2017). The package MICE (van Buuren and Groothuis-Oudshoorn 2011) was used to replace missing values through multiple imputation within variable groups after excluding individuals with > 25% missing data.

Analysis of scale counts and colour markings

Scale counts (VS, SCS, PTS and GS; Suppl. material 1: table S2) taken in the field or from photographs were analysed using a standard PCA on the correlation matrix of the data. Colour markings (LBS, LBL, LBW, NMS, NMW, NMUC, NMLC and NMPS; Suppl. material 1: table S2) were quantified with the number of coloured scales as the unit and analysed using a standard PCA on the correlation matrix of the data.

Linear discriminant analysis

Distance measurements, scale counts, measures for colour markings and three additional variables (RelRedPos, TPOS and BW; Suppl. material 1: table S2) were compared in a linear discriminant analysis to check whether microsatellite clusters are distinguishable. The function lda () from the package MASS (Venables and Ripley 2002) was used with equal priors and CV=TRUE. This was repeated for samples of N. helvetica only to test whether the discrimination could be improved.

Results

Distribution of mtDNA lineages

The distribution of mtDNA lineages in southern Switzerland is shown in Fig. 1a (for details of each sample, see Suppl. material 1: table S1). The cantons of Ticino and Valais are solely inhabited by Natrix helvetica with lineage C, typical for Alpine N. h. sicula. In contrast, the snakes in the cantons north of the Alps harbour mainly lineage E, typical for N. h. helvetica. However, north of the Alps a few lineage C individuals were recorded as well (cantons of Bern, Schwyz and Vaud). The northern shore of Lake Geneva seems to be a broad contact zone of lineages C and E. In addition, in this region, in Lausanne, some N. n. moreotica (mtDNA lineage 7) were caught, as already described by Dubey et al. (2017).

Microsatellite clusters

According to Pritchard and Wen (2002), the optimal number of clusters K has the highest likelihood L(K) value, which is here K=5 (Suppl. material 1: fig. S2, top left). In contrast, the ΔK method (Evanno et al. 2005) revealed K=2 as the best solution but also inferred a second pronounced peak for K=5 (Suppl. material 1: fig. S2, bottom right). The ΔK method reliably identifies the uppermost hierarchical level of genetic partitioning (Evanno et al. 2005); in our case, this corresponds to samples representing the two species of grass snake, N. helvetica and N. natrix. For inferring genotypic partitions (clusters) within N. helvetica and N. natrix, either subsets corresponding to each species can be examined separately in STRUCTURE or the STRUCTURE runs using the highest L(K) value can be inspected. Considering the highest L(K) value and the second peak of the ΔK approach, we present here the results for STRUCTURE runs using K=5. These five clusters correspond in our dataset to one cluster for N. natrix (i.e., genotypes of N. n. moreotica from Lausanne plus N. n. vulgaris from northeastern Switzerland) and four clusters within N. helvetica (for details, see Suppl. material 1: table S1).

Fig. 1b shows the geographic distribution of the five microsatellite clusters for western and southern Switzerland. One cluster (grey in Fig. 1b) represents the alien N. n. moreotica with mtDNA lineage 7 from Lausanne. Natrix helvetica sicula (mtDNA lineage C) is divided into two clusters, one in Valais (C-VS, red in Fig. 1b) and another one in Ticino (C-TI, green in Fig. 1b). Translated into the hydrographic net, the green cluster is found in river valleys connected to the great pre-Alpine Italian lakes and the Po drainage, and the red cluster is confined to the eastern (i.e., Alpine) Rhone drainage, from eastern Lake Geneva upstream. Natrix helvetica helvetica (mtDNA lineage E) is also divided into two clusters, one in the Swiss Plateau (Mittelland; E-ML, blue in Fig. 1b) and another cluster along Lake Geneva (E-GE, yellow in Fig. 1b). Admixture between clusters is evident in particular, but not only, in geographic contact zones, mainly the Lake Geneva region and the adjacent Rhone valley near Montreux. Most of these admixed snakes are N. helvetica (admixture among the four respective microsatellite clusters; 16 individuals from the cantons Geneva, Bern, Vaud and Ticino; Suppl. material 1: table S1). Another N. helvetica from the canton Bern (Langnau im Emmental) shows mito-nuclear discordance. This snake harbours a mitochondrial haplotype of lineage C combined with a microsatellite assignment to cluster E-ML. However, there are two further snakes with mito-nuclear discordance resulting from hybridization of N. helvetica with the alien N. n. moreotica in the region of Lausanne. These two snakes are having predominantly N. natrix genotypes combined with mitochondrial haplotypes of N. helvetica (lineage C; Suppl. material 1: table S1).

Genetic differentiation (FST) values are similar between clusters, except for the slightly lower values for E-GE/E-ML and for E-ML/natrix, whereas the highest value was observed between C-VS/E-GE (Suppl. material 1: table S3).

Morphology

For morphological and landmark analyses, snakes representing the clusters E-ML and E-GE were merged in one cluster E, while C-VS and C-TI were kept separate. C-VS and C-TI are geographically divided by mountainous regions difficult to cross for grass snakes, whereas E-ML and E-GE are in contact and admixing. Additionally, the number of samples for E-GE is very limited.

Mahalanobis and Procrustes permutation tests revealed significant morphological differences between all clusters (Table 1). However, Fig. 2 shows that there is always some overlap among all clusters, except for the lateral landmarks, where N. natrix was distinct. Also, there are no obvious shape differences visible in the wireframe graphs (Suppl. material 1: figs S3, S4).

Table 1.

Statistical test results of Canonical Variate Analysis (CVA) of lateral and dorsal landmarks (Fig. 2) using 10,000 permutations for each test. E, C-VS and C-TI represent Natrix helvetica clusters derived from microsatellite analyses (Fig. 1); natrix refers to samples from northeastern Switzerland (N. n. vulgaris, mtDNA lineage 3).

Lateral landmarks
Mahalanobis distances among groups Procrustes distances among groups
natrix E C-VS natrix E C-VS
E 4.8305 E 0.0308
C-VS 8.2776 5.6166 C-VS 0.0380 0.0285
C-TI 6.5247 4.5882 5.2595 C-TI 0.0431 0.0324 0.0237
P values from permutation tests P values from permutation tests
natrix E C-VS natrix E C-VS
E < 0.001 E < 0.001
C-VS < 0.001 < 0.001 C-VS < 0.001 < 0.001
C-TI < 0.001 < 0.001 < 0.001 C-TI < 0.001 < 0.001 0.0783
Dorsal landmarks
Mahalanobis distances among groups Procrustes distances among groups
natrix E C-VS natrix E C-VS
E 2.7561 E 0.0232
C-VS 3.5488 3.7053 C-VS 0.0236 0.0268
C-TI 4.4217 4.1624 2.9620 C-TI 0.0375 0.0389 0.0264
P values from permutation tests P values from permutation tests
natrix E C-VS natrix E C-VS
E < 0.001 E < 0.001
C-VS < 0.001 < 0.001 C-VS < 0.001 < 0.001
C-TI < 0.001 < 0.001 < 0.001 C-TI < 0.001 < 0.05 < 0.001
Figure 2. 

Canonical Variate Analysis (CVA) of dorsal (a) and lateral (b) landmark coordinates for grass snakes assigned to the microsatellite clusters of Natrix helvetica (E, C-VS, C-TI) and N. natrix (N. n. vulgaris, mtDNA lineage 3, from northeastern Switzerland). Only individuals with a genotypic cluster assignment ≥ 80% are included; circles represent 95% confidence ellipses.

Shape PCA and standard PCAs show no differentiation of mtDNA lineages E, C-VS and C-TI of N. helvetica (Fig. 3). Only N. natrix is slightly distinct because of smaller and differently shaped markings. Similar results are obtained for the linear discriminant analysis. Natrix natrix also represents here the most distinct cluster with a percentage of 86.7% of correctly allocated individuals. For all other clusters, the percentages are below 70% (Suppl. material 1: table S4). Mean values and standard deviations of each morphological trait are summarized for the two N. helvetica subspecies in Suppl. material 1: table S5.

Figure 3. 

Shape PCA of distance measurements (a), standard PCA of scale counts (b) and standard PCA of colour marking measurements (c) for grass snakes assigned to the microsatellite clusters of Natrix helvetica (C-TI, C-VS, E) and N. natrix (N. n. vulgaris, mtDNA lineage 3, from northeastern Switzerland). Only individuals with a genotypic cluster assignment ≥ 80% are included.

Discussion

After the discovery of two putatively alien mtDNA lineages of grass snake around Lausanne (lineages C and 7 of Kindler et al. 2013), concerns were raised about genetic pollution of native populations of Natrix helvetica (Dubey et al. 2017). It was clear that the common grass snakes harbouring mtDNA lineage 7, now identified with the subspecies N. natrix moreotica (see Asztalos et al. 2021b), are allochthonous and originate from the former Vivarium at Lausanne. However, for the other mitochondrial lineage (C), typical for Alpine representatives of N. helvetica sicula (Schultze et al. 2020; Asztalos et al. 2021a), there remained the possibility of a natural occurrence. Dubey et al. (2017) also assumed that grass snakes with clade C had escaped from the Vivarium in the 1970s, together with the ancestors of what is now called N. n. moreotica (in Dubey et al. 2017, the inappropriate name N. n. persa is still used; see Asztalos et al. 2021b). However, the Alpine lineage of N. n. sicula is widely distributed in northern Italy, including the Alps, and its range extends northwards across the Alps to Tyrol and southernmost Bavaria (Glaw et al. 2019; Schultze et al. 2020; Asztalos et al. 2021a; Neumann et al. 2024). Until recently, only two records for N. n. sicula from Switzerland were published, one from Ticino and the other from beyond the Simplon Pass (Niedergesteln, Valais; Kindler and Fritz 2018; Schultze et al. 2020). Yet, these records and the wide Alpine distribution of N. h. sicula suggest that this subspecies could not only occur naturally in Switzerland, but that it has a much wider distribution than previously thought. This is supported by our present investigation.

Our study shows that N. h. sicula is widely distributed in the cantons of Ticino and Valais and ranges to the canton of Vaud, along Lake Geneva, where it hybridizes with N. h. helvetica, as evidenced by microsatellite genotypes. Therefore, N. h. sicula should no longer be considered as non-native around Lausanne, as supposed by Dubey et al. (2017) for the records of mtDNA lineage C there.

The two microsatellite clusters of N. h. helvetica and N. h. sicula correspond to local population structure. It is possible that our cluster E-GE (Fig. 1b) matches the southern cluster revealed by Asztalos et al. (2020), who also identified two distinct clusters for N. h. helvetica. Their southern cluster from France could reach upstream in the Rhone Valley to Lake Geneva, whereas our cluster E-ML could match the northern cluster of Asztalos et al. (2020). However, the sample from Geneva in Asztalos et al. (2020) was also assigned to the northern group and the differentiation of the cluster E-GE could instead also be related to a local effect (isolation by distance) and the tendency of STRUCTURE to cluster samples with similar ancestry. In contrast, until now, no substructure has been described for Alpine N. h. sicula, but the distribution pattern of our two clusters C-VS and C-TI makes sense biogeographically. Each cluster matches another drainage system (Alpine Rhone Valley vs. Po drainage).

North of Lake Geneva, we were not only able to detect the natural contact and hybridization zone of the two subspecies of N. helvetica. Our microsatellite data (Fig. 1b; Suppl. material 1: table S1) also confirm genetic pollution of N. helvetica from alien N. n. moreotica. Nature conservation should take action to prevent wider introgression of alien genes in this region. Non-native snakes and their hybrids should be removed before they further compromise the genetic identity of the native populations of N. helvetica, which is currently classified as “endangered” in the most recent Swiss Red List of Reptiles (BAFU & info fauna 2023).

The distribution and hybrid zone of N. h. helvetica and N. h. sicula in Switzerland can be explained in a biogeographical context. Kindler et al. (2018b) suggested that the nominotypical subspecies survived the last glaciation in southern France, from where it expanded its range to more northern regions during the Holocene warming. Natrix helvetica sicula, on the other hand, was inferred to have survived the last glaciation in a distinct ‘microrefugium’ in northeastern Italy (Kindler et al. 2013; Schultze et al. 2020), but it was impossible to endure in the Alpine Rhone Valley, which was completely covered by ice during the Late Glacial Maximum (Seguinot et al. 2018). As a consequence of the Holocene warming, N. h. sicula crossed the main Alpine chain not only in the east to Tyrol and southern Bavaria (Glaw et al. 2019; Asztalos et al. 2021a), but also in the west, as we know now. There, it established a secondary contact and hybrid zone with N. h. helvetica along Lake Geneva. This scenario implies transalpine dispersal at altitudes of about 2000–2200 m a.s.l., i.e., at altitudes that were only very recently colonized by the species (database of info fauna – karch, which hosts more than 25,000 records for grass snakes in Switzerland and is the national reference centre for Swiss amphibians and reptiles). Thus, it can be speculated that N. h. sicula crossed the western Alps during a Holocene period, which was at least as warm as present, probably via the Simplon Pass (approx. 2000 m a.s.l.).

In this context, two other northern records of grass snakes with lineage C are difficult to interpret (Langnau im Emmental, canton Bern, and Muotathal, canton Schwyz, directly at the border to canton Uri). For the snake from Muotathal no genotype is available, but the snake from Langnau im Emmental is, according to its microsatellite genotype, a pure representative of the cluster E-ML. Also, two other snakes from the canton Bern have an admixed genotype (Fig. 1; Suppl. material 1: table S1), suggesting that some N. h. sicula left indeed their genetic footprint there. It is well known that grass snakes are spread with building material, etc. or have been voluntarily moved (e.g., Ahnelt et al. 2021; Asztalos et al. 2021c), so these northern records are not necessarily evidence for wide-reaching natural dispersal and introgression across the Alps.

In contrast to genetic data, our morphological analyses revealed only a weak differentiation among the studied grass snakes. Only the two species N. helvetica and N. natrix could be morphologically discriminated with some confidence, while the used traits were not helpful in discriminating the two subspecies of individual genetic clusters within N. helvetica. This does not contradict the validity of the involved two subspecies because morphological traits that can be distinguished by humans are neither necessarily biologically relevant nor a prerequisite for taxonomic distinctness (compare, for instance, Kindler and Fritz 2018; Dufresnes et al. 2023, 2024). Also, we could have missed some relevant traits in coloration and pattern that were only recently highlighted (Fritz et al. 2023). According to these authors, some individuals of N. h. sicula show a “spotted” colour pattern that never occurs in the nominotypical subspecies. It cannot be excluded that further traits exist that help to identify the different subspecies even in the field, but it remains a challenge to disentangle this situation.

Acknowledgements

We thank M. Chèvre for advice, data, photos and his lab protocols and V. Zwahlen for assistance with processing microsatellites. Special thanks go to all researchers and institutions that provided samples: M. Chèvre, S. Dubey, E. Gallice, T. Gil, G. Meier, B. Meister, D. Muri, the Musée de la nature Sion (S. Gerber), the Museo cantonale di storia naturale Lugano (N. Zambelli), and the Naturhistorisches Museum Bern (R. Hagmann). We thank all authorities of the cantons of Vaud, Valais and Ticino for the necessary permits. The figures were enhanced by M. H. Fischer. An anonymous reviewer, P. Mikulíček and the editor R. Jadin provided helpful comments on the manuscript.

Additional information

Conflict of interest

The authors have declared that no competing interests exist.

Ethical statement

Captures in the wild were conducted according to the local authorities’ requirements (veterinary authorisation number GR_2013_15kü).

Funding

This study was funded by the Swiss Bundesamt für Umwelt (BAFU).

Author contributions

Conceptualization: SU. Data curation: AS, SU. Formal analysis: AS, HB. Funding acquisition: SU. Investigation: AS. Methodology: SU. Project administration: SU. Resources: HB, STH. Supervision: SU. Visualization: AS, UF. Writing – original draft: AS, UF. Writing – review and editing: AS, HB, STH, SU, UF.

Author ORCIDs

Hannes Baur https://orcid.org/0000-0003-1360-3487

Stefan T. Hertwig https://orcid.org/0000-0003-1617-8180

Uwe Fritz https://orcid.org/0000-0002-6740-7214

Sylvain Ursenbacher https://orcid.org/0000-0001-5093-6403

Data availability

All data that support the findings of this study are available in the main text or the Suppl. material. No new haplotypes were identified in the present study, which is why no sequences were uploaded to the European Nucleotide Archive (ENA).

References

  • Ahnelt H, Romanova T, Klinge A, Böhme W, Fritz U, Asztalos M (2021) The common grass snake (Natrix natrix) on Sylt: Human-mediated colonization of a North Sea island. Salamandra 57: 285–290.
  • Asztalos M, Ayaz D, Bayrakcı Y, Afsar M, Tok CV, Kindler C, Jablonski D, Fritz U (2021b) It takes two to tango – Phylogeography, taxonomy and hybridization in grass snakes and dice snakes (Serpentes: Natricidae: Natrix natrix, N. tessellata). Vertebrate Zoology 71: 813–834. https://doi.org/10.3897/vz.71.e76453
  • Asztalos M, Glaw F, Franzen M, Kindler C, Fritz U (2021a) Transalpine dispersal: Italian barred grass snakes in southernmost Bavaria – This far but no further! Journal of Zoological Systematics and Evolutionary Research 59: 1136–1148. https://doi.org/10.1111/jzs.12471
  • Asztalos M, Schultze N, Ihlow F, Geniez P, Berroneau M, Delmas C, Guiller G, Legentilhomme J, Kindler C, Fritz U (2020) How often do they do it? An in-depth analysis of the hybrid zone of two grass snake species (Natrix astreptophora and Natrix helvetica). Biological Journal of the Linnean Society 131: 756–773. https://doi.org/10.1093/biolinnean/blaa152
  • Asztalos M, Wielstra B, Struijk RPJH, Ayaz D, Fritz U (2021c) Aliens in the Netherlands: Local genetic pollution of barred grass snakes (Squamata: Serpentes: Natricidae). Salamandra 57: 173–179.
  • BAFU & info fauna (2023) Rote Liste der Reptilien. Gefährdete Arten der Schweiz. Bundesamt für Umwelt BAFU, Bern, 32 pp.
  • Baur H (2024) Morphometric methods for species recognition. In: Heraty JM, Woolley JB (Eds) Chalcidoidea of the World. CABI Press, Oxford, in press.
  • Baur H, Kranz-Baltensperger Y, Cruaud A, Rasplus J-Y, Timokhov AV, Gokhman VE (2014) Morphometric analysis and taxonomic revision of Anisopteromalus ruschka (Hymenoptera: Chalcidoidea: Pteromalidae) – An integrative approach. Systematic Entomology 39: 691–709. https://doi.org/10.1111/syen.12081
  • Chèvre M (2015) The contact zone of the grass snake (Natrix natrix) in Switzerland. Master thesis, Université de Neuchâtel, Neuchâtel, 97 pp.
  • Dubey S, Ursenbacher S, Schuerch J, Golay J, Aubert P, Dufresnes C (2017) A glitch in the Natrix: Cryptic presence of alien grass snakes in Switzerland. Herpetology Notes 10: 205–208.
  • Dufresnes C, Ambu J, Galán P, Sequeira F, Viesca L, Choda M, Álvaret D, Alar B, Suchan T, Künzel S, Martínez-Solano I, Vences M, Nicieza A (2024) Delimiting phylogeographic diversity in the genomic era: Application to an Iberian endemic frog. Zoological Journal of the Linnean Society: zlad170. https://doi.org/10.1093/zoolinnean/zlad170
  • Dufresnes C, Poyarkov N, Jablonski D (2023) Acknowledging more biodiversity without more species. Proceedings of the National Academy of Sciences of the United States of America 120: e2302424120. https://doi.org/10.1073/pnas.2302424120
  • Earl DA, vonHoldt BM (2012) STRUCTURE HARVESTER: A website and program for visualizing structure output and implementing the Evanno method. Conservation Genetics Resources 4: 359–361. https://doi.org/10.1007/s12686-011-9548-7
  • Fritz U, Grismer LL, Asztalos M (2023) Hybrid zones of Natrix helvetica and N. natrix: Phenotype data from iNaturalist and genetics reveal concordant clines and the value of species-diagnostic morphological traits. Vertebrate Zoology 73: 383–395. https://doi.org/10.3897/vz.73.e103319
  • Fritz U, Schmidtler JF (2020) The Fifth Labour of Heracles: Cleaning the Linnean stable of names for grass snakes (Natrix astreptophora, N. helvetica, N. natrix sensu stricto). Vertebrate Zoology 70: 621–665. https://doi.org/10.26049/vz70-4-2020-07
  • Gebiola M, Monti MM, Johnson RC, Woolley JB, Hunter MS, Giorgini M, Pedata PA (2017) A revision of the Encarsia pergandiella species complex (Hymenoptera: Aphelinidae) shows cryptic diversity in parasitoids of whitefly pests. Systematic Entomology 42: 31–59. https://doi.org/10.1111/syen.12187
  • Glaw F, Franzen M, Oefele M, Hansbauer G, Kindler C (2019) Genetischer Erstnachweis, Verbreitung und südalpine Herkunft der Barrenringelnatter (Natrix helvetica spp.) in Bayern. Zeitschrift für Feldherpetologie 26: 1–20.
  • Huber C, Baur H (2016) Nebria (Patrobonebria) incognita n. sp. and Nebria (P.) hiekeiana n. sp., two new species from the Western Himalaya, with remarks on Nebria (P.) desgodinsi (Coleoptera, Carabidae, Nebriinae). Entomologische Blätter und Coleoptera 112: 203–214.
  • Kabisch K (1999) Natrix natrix (Linnaeus, 1758) – Ringelnatter. In: Böhme W (Ed.) Handbuch der Reptilien und Amphibien Europas. Band 3/IIA: Schlangen II. Aula-Verlag, Wiebelsheim, 513–580.
  • Kindler C, Böhme W, Corti C, Gvoždík V, Jablonski D, Jandzik D, Metallinou M, Široký P, Fritz U (2013) Mitochondrial phylogeography, contact zones and taxonomy of grass snakes (Natrix natrix, N. megalocephala). Zoologica Scripta 42: 458–472. https://doi.org/10.1111/zsc.12018
  • Kindler C, Chèvre M, Ursenbacher S, Böhme W, Hille A, Jablonski D, Vamberger M, Fritz U (2017) Hybridization patterns in two contact zones of grass snakes reveal a new Central European snake species. Scientific Reports 7: 7378. https://doi.org/10.1038/s41598-017-07847-9
  • Kindler C, de Pous P, Carranza S, Beddek M, Geniez P, Fritz U (2018a) Phylogeography of the Ibero-Maghrebian red-eyed grass snake (Natrix astreptophora). Organisms, Diversity & Evolution 18: 143–150. https://doi.org/10.1007/s13127-017-0354-2
  • Kindler C, Fritz U (2018) Phylogeography and taxonomy of the barred grass snake (Natrix helvetica), with a discussion of the subspecies category in zoology. Vertebrate Zoology 68: 269–281. https://doi.org/10.3897/vz.68.e31615
  • Kindler C, Graciá E, Fritz U (2018b) Extra-Mediterranean glacial refuges in barred and common grass snakes (Natrix helvetica, N. natrix). Scientific Reports 8: 1821. https://doi.org/10.1038/s41598-018-20218-2
  • László Z, Baur H, Tótmérész B (2013) Multivariate ratio analysis reveals Trigonoderus pedicellaris Thomson (Hymenoptera, Chalcidoidea, Pteromalidae) as a valid species. Systematic Entomology 38: 753–762. https://doi.org/10.1111/syen.12026
  • Neumann A, Asztalos M, Fritz U, Glaw F (2024) A spotlight on the hybrid zone of grass snakes (Natrix helvetica sicula and Natrix natrix) in southern Bavaria – the Prien Valley. Salamandra 60: 17–28.
  • Pokrant F, Kindler C, Ivanov M, Cheylan M, Geniez P, Böhme W, Fritz U (2016) Integrative taxonomy provides evidence for the species status of the Ibero-Maghrebian grass snake Natrix astreptophora. Biological Journal of the Linnean Society 118: 873–888. https://doi.org/10.1111/bij.12782
  • Pritchard JK, Wen W (2002) Documentation for STRUCTURE Software: Version 2. University of Chicago Press, Chicago, IL, 29 pp.
  • R Core Team (2017) R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna. https://www.R-project.org/
  • Rohlf FJ (2013) TPSUtility program, version 1.56. Department of Ecology and Evolution, State University of New York, Stony Brook, NY.
  • Rohlf FJ (2017) TPSDig2, version 2.30. Department of Ecology and Evolution, State University of New York, Stony Brook, NY.
  • Schultze N, Spitzweg C, Corti C, Delaugerre M, Di Nicola MR, Geniez P, Lapini L, Liuzzi C, Lunghi E, Novarini N, Picariello O, Razzetti E, Sperone E, Stellati L, Vignoli L, Asztalos M, Kindler C, Fritz U (2020) Mitochondrial ghost lineages blur phylogeography and taxonomy of Natrix helvetica and N. natrix in Italy and Corsica. Zoologica Scripta 49: 395–411. https://doi.org/10.1111/zsc.12417
  • Sidlauskas BL, Mol JH, Vari RP (2011) Dealing with allometry in linear and geometric morphometrics: A taxonomic case study in the Leporinus cylindriformis group Characiformes: Anostomidae) with description of a new species from Suriname. Zoological Journal of the Linnean Society 162: 103–130. https://doi.org/10.1111/j.1096-3642.2010.00677.x
  • Thorpe RS (1979) Multivariate analysis of the population systematics of the ringed snake, Natrix natrix (L). Proceedings of the Royal Society of Edinburgh 78B: 1–62. https://doi.org/10.1017/S026972700001294X
  • van Riemsdijk I, Struijk RPJH, Pel E, Janssen IAW, Wielstra B (2020) Hybridisation complicates the conservation of Natrix snakes in the Netherlands. Salamandra 56: 78–82.
  • Waser LE, Schweizer M, Haas A, Hertwig ST (2017) From a lost world: An integrative phylogenetic analysis of Ansonia Stoliczka, 1870 (Lissamphibia: Anura: Bufonidae), with the description of a new species. Organisms Diversity & Evolution 17: 287–303. https://doi.org/10.1007/s13127-016-0294-2

Supplementary material

Supplementary material 1 

Supplementary information

Andreas Schild, Hannes Baur, Stefan T. Hertwig, Uwe Fritz, Sylvain Ursenbacher

Data type: docx

Explanation note: fig. S1. Landmark positions. fig. S2. STRUCTURE HARVESTER results indicating the optimal number of microsatellite clusters (K). fig. S3. Wireframe graph showing shape changes along the canonical variate 1 (CV1) for dorsal landmarks. fig. S4. Wireframe graph showing the shape changes along the canonical variate 1 (CV1) for lateral landmarks. table S2. Morphological variables measured. Illustrations from Chèvre (2015). table S3. Pairwise FST values between microsatellite clusters of Natrix helvetica (E-GE, E-ML, C-VS and C-TI) and N. natrix. table S4. Morphological assignment of specimens (horizontal rows) to microsatellite clusters of Natrix natrix and N. helvetica (E, C-VS, C-TI, vertical columns) based on all morphological variables using a Linear Discriminant Analysis (LDA). table S5. Mean and standard deviation (SD) of different morphological traits for Swiss barred grass snakes (Natrix helvetica).

This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (3.29 MB)
login to comment