Ecological significance of dorsal polymorphism in a Batesian mimic snake.

Batesian mimicry is the process in which harmless species adopt the appearance of a dangerous, aposematic species. In some prey species, both Batesian mimetic and non-Batesian morphs coexist, presupposing that both morphs have to be evolutionarily advantageous. The viperine snake, Natrix maura, exhibits a zigzag dorsal pattern and antipredatory behavior that mimics European vipers. This snake also has a striped dorsal pattern that coexists with the zigzag pattern. We have examined whether individuals belonging to different geographically structured clades were more likely to exhibit a certain dorsal pattern, and whether the zigzag pattern has a protective function by exposing artificial snakes to predation in natural environments, in addition to comparing antipredatory behavior between zigzag and striped snakes also in natural environments. Our results indicate that the striped pattern was not geographically structured, but habitat-dependent. Aerial predators less frequently attacked zigzag plasticine models than striped or unpatterned models. We detected a shift in antipredator behavior between the 2 morphs, as Batesian mimicking N. maura responded to an approaching potential predator by remaining immobile or fleeing at shorter distances than did striped ones. We conclude that Batesian mimics maintain the cryptic and aposematic value by resembling vipers, whereas in open habitats the non-Batesian mimic has altered its antipredator behavior to maintain its fitness.

In the predator–prey context, prey have developed multiple strategies to survive (Ruxton et al. 2004). Among others, aposematism consists in the development of warning signals designed to dissuade potential predators (Poulton 1890; Rowe and Guilford 2000). Interestingly, some species have evolved to mimic aposematic species to gain some functional advantage. In Batesian mimicry (Bates 1862), harmless species adopt the appearance of a dangerous, aposematic species, to protect themselves from predation, although they are really inoffensive (Maynard Smith 1998). Batesian mimicry appears to be a widespread taxonomic phenomenon (Ruxton et al. 2004). Edmunds (1974) provides many examples of mimetic associations in butterflies, some of them being Batesian mimics, and many spider species mimic ants (Mclver and Stonedahl 1993). Among vertebrates, there are examples among fish (McCosker 1977), amphibians (Kuchta 2005), and birds (Rowe et al. 1986), although Batesian mimicry is rare in mammals (but see Pough 1988). In this study, we address an example of Batesian mimicry in snakes.

Snakes are highly polymorphic in dorsal patterning and an excellent model group to address the evolutionary drivers of dorsal pattern design (Cox and Davis Rabosky 2013), with camouflage, aposematism, and thermoregulation being the major evolutionary drivers of dorsal pattern diversification (Allen et al. 2013). Moreover, many inoffensive snake species resort to Batesian mimicry by reproducing the aposematic dorsal pattern of dangerous species (Brodie and Brodie 2004). A classic example of aposematism occurs in European vipers of the genus Vipera, which present a distinctive set of morphological and behavioral characteristics, such as a short and stout body, head triangulation, coiling, snorting, and zigzag dorsal pattern, which their potential predators relate to their high intensity of defences (Wüster et al. 2004; Niskanen and Mappes 2005; Valkonen et al. 2011a). The zigzag dorsal pattern is likely a multifunctional protective coloration; it provides protection against predation because of its aposematic role, although this pattern in certain settings would also favor crypsis (Cott 1940; Andrén and Nilson 1981; Shine and Madsen 1994; Santos et al. 2014).

The elaborate and consistent set of morphological and behavioral traits of most European viper species has been copied by the inoffensive viperine snake Natrix maura (Linnaeus 1758) (Rollinat 1934) as defensive mimicry. This water snake is distributed by the Western Palaearctic and mimics morphological and behavioral traits of European vipers with which it coincides geographically to some degree, for example, Vipera aspis (Linnaeus 1758), V. latastei Boscá 1878, and V. seoanei Lataste 1879 (Aubret and Mangin 2014; Santos et al. 2014). For example, N. maura exhibits a zigzag dorsal pattern (hereafter zigzag) formed by a wide and dark dorsal band, zigzag shaped, from the neck to the tail, which contrasts with lighter background body colors. Moreover, N. maura displays head triangulation, a spring-like body posture (coiling), snorting, body inflation, and striking (Rollinat 1934).

Notably, N. maura also shows dorsal pattern variability, since some specimens have a striped pattern of 2 dorso-lateral light bands (longitudinally striped sensuvon Helversen et al. 2013), or bilineata dorsal pattern (hereafter bilineata) (Schätti 1982). Based on mitochondrial markers, these bilineata individuals are genetically indistinguishable from the zigzag individuals of the same population (Carranza S, personal communication). Intraspecific variability in snakes with one of the patterns being striped is quite common in snakes, especially in Colubridae (Wolf and Werner 1994), and can be geographically structured or sympatric (Kark et al. 1997). The striped pattern tends to enhance a visual escape strategy, as snakes can confuse predators by not providing body references in tracking its movement (Pough 1976; Allen et al. 2013); this dazzle effect, the motion dazzle camouflage, has recently been supported experimentally (Hogan et al. 2016). Accordingly, compared with non-striped snakes, striped ones usually show a high frequency of broken tails as a consequence of failed predation events (Pleguezuelos et al. 2010).

In the absence of temporal variation in selective forces, the maintenance of polymorphism (variability within a geographical location) and of polyphenism (variability between geographical locations) presupposes that the morphs should thrive under different environments (Wolf and Werner 1994; Shine and Harlow 1998). For N. maura, the coexistence (syntopy) of bilineata and zigzag snakes also suggests the use of different antipredator strategies (i.e., fast-moving vs. static snakes; Allen et al. 2013). Wolf and Werner (1994) suggested that striped snakes tend to thrive in open habitats, while those with alternative dorsal patterns (i.e., blotched) prosper in structurally more complex habitats (see also Pough 1976). Unfortunately, studies on the ecological drivers and the adaptive value of intraspecific dorsal variability in snakes under an experimental or comparative framework have been scarcely addressed (but see Andrén and Nilson 1981; Shine and Harlow 1998; Cox and Davis Rabosky 2013).

The objective of this study is to uncover the evolutionary drivers of the dorsal pattern dimorphism in N. maura. While the studies on the function of dorsal patterns in model species (e.g., true vipers) are common (see review in Allen et al. 2013), far fewer studies have analyzed the same topic in mimetic species (but see Greene and McDiarmid 1981; Cox and Davis Rabosky 2013). Similarly, few studies have addressed antipredator responses against predators other than humans (Gregory 2016). We specifically address the following questions:

(1) Is the occurrence of zigzag and bilineata dorsal patterns in N. maura geographically structured? Schätti (1982) showed that the bilineata pattern appeared fundamentally in coastal areas, increasingly to the south and being completely absent in some populations. Later, Guicking et al. (2002) studied the phylogeography of the species, providing a genetic framework for the analysis of the spatial distribution of this trait. These authors detected 3 geographically structured lineages from the Pliocene, 2 in Africa, and one in Europe. If the occurrence of the bilineata pattern in N. maura was not geographically structured and occurred in all the lineages, we hypothesize that the dorsal dimorphism would have no relationship with the phylogeography of the species.

(2) Is the relative frequency of zigzag and bilineata patterns associated with any landscape appearance? In some snakes, strong selection promotes variability across spatial scales (Cox and Davis Rabosky 2013). Despite the latitudinal cline of the bilineata frequency cited by Schätti (1982), geographically close populations exhibit opposite frequencies of the zigzag and bilineata patterns (Duguy and Saint-Girons 1993). Based on the field experience of the authors, the striped morphotype appears to be more common in wetlands and other open landscapes (see also Duguy and Saint-Girons 1993; Santos 2009), but there are no studies on the association between habitat type and the frequency of the 2 dorsal pattern across the range of N. maura.

(3) Have bilineata snakes lost the Batesian-mimic role by no longer resembling a viper? The zigzag pattern in N. maura is presumed to have an aposematic, even cryptic function, due to its viper mimesis, while the bilineata dorsal pattern in snakes does not (Brodie 1992; Wolf and Werner 1994; Allen et al. 2013). Thus, in the absence of behavioral modulation (i.e., by snake plasticine models) and maintaining similar background, we hypothesize that the non-Batesian (bilineata) individuals will be more frequently attacked by predators than will the Batesian (zigzag) ones.

(4) Do bilineata and zigzag N. maura display different antipredatory behavior when confronted with predators? In snakes, dorsal pattern and antipredator behavior interact (Brodie 1992), and dorsal pattern diversification is related mainly to behavior, rather that habitat choice (Allen et al. 2013). For this reason, the maintenance of striped individuals in a population is expected when this pattern can confer some antipredatory and other evolutionarily advantages (Brodie 1992; Wolf and Werner 1994). Thus, we hypothesize that, linked to the loss of the viper mimesis in striped N. maura, individuals would display different antipredatory strategies with regard to their dorsal pattern: whereas Batesian mimic (zigzag) individuals are expected to mirror viper strategy, that is, by staying immobile, non-aposematic (bilineata) individuals are expected to try to escape, as occurs in other striped snakes (Jackson et al. 1976; Brodie 1992; Wolf and Werner 1994; Pleguezuelos et al. 2010; Allen et al. 2013).

Materials and Methods

Geography and spatial ecology of the variability

The spatial distribution of the variability of N. maura was examined by the analysis of bilineata and zigzag dorsal pattern frequencies in 917 individuals from 12 populations (Table 1), covering all the phylogenetic clades (Guicking et al. 2002) and most of the range of the species (Figure 1A). We examined vouchers from scientific collections (nine localities) and field records of the authors and colleagues from specific sites (3 localities) corresponding to either of the 2 contrasting natural landscapes where the species thrives (Santos 2009), that is, open and not-open habitats, with sample size >15 (Table 1). Open habitats were marshlands and peat bogs where vegetation, depending on site, was composed of broadleaf cattail Typha latifolia, common reed Phragmites australis, bulrush Scirpus holoschoenus, great fen-sedge Cladium mariscus, perennial saltwort Sarcocornia fruticosa, and sedges (Carex sp.), reaching a height of 0.4–1.5 m, with grassland and reeds dominating in the banks. Closed habitats were rivers and streams with banks dominated by stones and thick vegetation, and, depending on sites, made up of willows (Salix sp.), silver poplar Populus alba, elm Ulmus minor, ash Fraxinus angustifolia, honeysuckle Lonicera arborea, and elmleaf blackberry Rubus ulmifolius, with a height of 1.0–15.0 m (Figure 1B,C).

Table 1.

Localities considered in this study for the analysis of the dorsal pattern in the viperine snake, Natrix maura, within the Western Mediterranean

Locality	Latitude	Longitude	Habitat	Clade	Distance to viper (km)	N	ZZ	BL	Sampling
Algeria	36°46′ N	8° 13′ E	Open	Algerian	30	16	10	6	EBD
Rif Mountains	35°07′ N	5°17′ W	Not-open	Moroccan	20	41	40	1	Private collection S. Fahd
Baetic Mountains	37°12′ N	3°42′ W	Not-open	Iberian	15	77	76	1	DBAG
Burgos province	42°45′ N	3°48′ W	Not-open	Iberian	2	21	21	0	EBD
Cádiz province	36°35′ N	6°14′ W	Open	Iberian	26	28	14	14	Field sampling, S. Busack
Cazorla Mountains	37°54′ N	2°56′ W	Not-open	Iberian	1	28	28	0	EBD
Delta Ebro	40°42′ N	0°45′ E	Open	Iberian	35	318	183	135	Field sampling
Doñana	36°59′ N	6°26′ W	Open	Iberian	5	113	69	44	EBD
Río Frío River	37°09′ N	4°12′ W	Not-open	Iberian	3	147	147	0	DBAG
Matarranya River	41°07′ N	0°12′ E	Not-open	Iberian	10	50	50	0	Field sampling
Charca Suárez	36°43′ N	3°12′ W	Open	Iberian	18	37	21	16	Field sampling/DBAG
Puertollano pools	38°39′ N	4°06′ W	Open	Iberian	15	41	38	3	Field sampling

Notes: We include the geographic coordinates for the barycentre of the localities, the habitat (2 possibilities, open when marshlands or peat bogs, and closed when others), the phylogenetic clade, the distance to the closest viper population (Vipera latastei), the sample size (N), the frequency of the bilineata (BL) and zigzag (ZZ) dorsal pattern, and the provenance of the sample. Abbreviations are as follows: EBD, Estación Biológica de Doñana, Sevilla, Spain; DBAG, departamento de biología animal, Granada University, Spain; CRBA, Centre de Recursos de Biodiversitat Animal, Barcelona University, Spain.

Figure 1.

Geographic distribution of snakes belonging zigzag and bilineata patterns in the 12 sites studied along the Natrix maura range (A). On the top, examples of freshwater closed (B) and open (C) habitats where N. maura inhabits. On the bottom, living examples of snakes with the 2 dorsal patterns (D, E; photos by Raúl León). Pie charts depict the frequency of bilineata (red) and zigzag patterns (gray). Locality names in red refer to open (wetland) habitats, and in black to closed (river) habitats.

For each individual examined, body size (snout–vent length, SVL, to the nearest millimeter), sex (according to Feriche et al. 1993), and dorsal pattern (2 possibilities, zigzag and bilineata, Figure 1D,E) were recorded. Although in Batesian mimicry systems mimic and model do not need to be strictly sympatric (Pfennig and Mullen 2010), it is also true that if the number of models decreases, the benefits gained by mimics also decrease (Valkonen and Mappes 2014). The closest viper species to all populations considered in this study is V. latastei, both in North Africa and the Iberian Peninsula, and both species shared the visually oriented and most specialized predator (snake-eater) in the study area, Circaetus gallicus; thus, V. latastei is the model species mimicked by N. maura. Vipera latastei is a rather shy species with abundance hard to measure. However, to include the effect of the model species in the Batesian system, we have considered in our statistical analysis the distance of the closest viper population to each N. maura population, based on our field experience, and by the tool for measuring distances in Google Earth (Table 1). Although snake pigmentation fades somewhat with preservation in liquid, all museum specimens were easily classified as zigzag or bilineata. Given that snake vouchers shrink due to fixative and the preservative process (Barry 2011), we estimated SVL for museum vouchers without fresh measurements by applying the regression of fresh against preserved SVL for this species (Santos et al. 2011).

Predation experiment with plasticine models

To test whether bilineata snakes lost their Batesian mimicry role compared with zigzag snakes, we performed a field-based predation experiment with plasticine models (see similar experiments in Wüster et al. 2004; Niskanen and Mappes 2005; Valkonen et al. 2011b). We made 12 bilineata and 12 zigzag plasticine models of snakes plus 12 patternless models as control (thereafter control). Models were made by melting non-toxic plasticine (JOVER®) in a silicone mould of a N. maura voucher from the study area (355 mm SVL, the average body size for this population; N = 290). The background body color of the models was made by melting brown, green, and white plasticine colors (ratio 10:5:1, respectively). The pattern in zigzag models was painted with a permanent marker (EDDING 500®), 5 mm width, while the 2 mm wide dorso-lateral stripes of the bilineata model were made by melting white, green, and yellow plasticine colors (ratio 5:2:1, respectively). Control models exhibited only the ground color (Figure 2).

Figure 2.

Plasticine models used for the experimental study of predation upon the viperine snake, Natrix maura. From top to bottom, zigzag, striped (bilineata) and uncolored (control) model.

Experiments with plasticine models were conducted in 2 localities from the south-eastern Iberian Peninsula, where both dorsal patterns coexist, Charca Suárez, a 28 ha marshland, 200 m from the sea shore (3°32′W, 36°43′N, 1 m altitude), and Padul, a 85 ha peat bog, inland and 30 km from the other site (3°36′W, 37°00′N; 725 m altitude). The landscape of both study sites was formed by wetland vegetation, such as T. latifolia, P. australis, S. holoschoenus, the introduced giant cane Arundo donax, and C. mariscus in Charca Suárez. Based on the list of N. maura predators (Santos et al. 2011), during the study period, nine potential avian predators were reported in Charca Suárez (Ardea cinerea, A. purpurea, Ardeola ralloides, Bubulcus ibis, C. gallicus, Egretta garzetta, Larus michahellis, L. ridibundus, Porphyrio porphyrio; from the data of the authors), and seven potential avian predators in Padul (A. cinerea, A. purpurea, B. ibis, C. gallicus, Circus aeruginosus, Egretta garzetta, Nycticorax nycticorax; Pérez-Contreras J, personal communication). Attacks were recorded as made by mammals when tooth marks were observed, and recorded as made by birds when bill marks were detected (Niskanen and Mappes 2005). As predators, birds are primarily visually oriented while mammals are primarily odour oriented. For this reason, marks caused by birds and mammals were recorded and analyzed separately. To avoid insect parasitism and damage, models were sprayed with insect repellent (Autan®; Valkonen et al. 2011b), and to avoid model displacement by predators, they were pegged to the soil by an iron wire that ran the full length of the model.

Field experiments were performed in 2013 in Padul, and 2013, 2015, and 2016 in Charca Suárez, during May and June, a period with high diurnal activity of N. maura (Santos and Llorente 2001). We conducted 12 trials (4 in Padul and 8 in Charca Suárez) to check predation on plasticine models. At each trial, 27–36 models (9–12 sets of models with one model of each dorsal pattern per set) were placed on a grassland background on the banks of the wetlands. Sets were placed 20–50 m apart, and models were spaced roughly 6 m apart, in random order. Models were placed in the early morning (08.00–09.00), and experimental trials lasted from 2 to 5 days, with no differences in the duration of trials between Padul and Charca Suárez (Mann–Whitney U-test; Z = 0.17, P = 0.9). Models were checked on a daily basis in late afternoon (19.00–20.00) during the trials. More than one attack on a single model during a trial was considered a single-attack event to reduce pseudoreplication and bias caused by potential multiple attacks from an individual predator on a single plasticine model (Valkonen et al. 2011b). Thus, the sampling units were the models in each trial. A total of 104 models were examined in Padul (4 trials × 27 models minus 4 models stolen by humans), and 276 in Charca Suárez (6 trials × 36 models, and 2 trials with 30 models).

Antipredatory behavior experiments

As the previous predation experiments with plasticine models suggested that bilineata snakes lost the advantages of the mimesis with vipers (see results below), we also examined differences in antipredatory behavior between bilineata and zigzag individuals. For this, we conducted a field experiment in Charca Suárez, a study area harboring both types of dorsal patterned N. maura. Surveys were conducted by a single researcher along a 1,500-m transect in the wetland banks, during May and June 2013, 2015, and 2016 (90 h of surveys). Snakes were searched when basking very close to the water (known as the K thermoregulatory strategy; Hailey and Davies 1987). When a snake was spotted, it was observed through binoculars (8 × 30), its body size was estimated (SVL, to the nearest cm, after training), and its dorsal pattern was classified as zigzag or bilineata. The sex of the individuals was disregarded because of the lack of sexual differences in the frequency of both dorsal pattern (Duguy and Saint-Girons 1993) and in the burst speed in lab trials on defence behavior (Hailey and Davies 1986). Afterward, the researcher approached slowly and directly at a constant speed of 0.8 m s−1, until the snake began to flee from the researcher. The fleeing distance was measured by a metric tape to the nearest 10 cm.

Body temperature can affect reptile behavior, particularly when deciding whether to move or remain still when confronted by a predator (Hertz et al. 1982). Unfortunately, some basking individuals (all bilineata) fled to deep water when the researcher approached, precluding capturing them to record body temperature. Thus, we recorded air temperature 1 m above ground (Hibok 14 thermometer, to the nearest 0.1°C) as a surrogate of snake body temperature. Although body temperature of basking individuals of this species is usually above air temperature, there is a good correlation between the 2 temperatures (Hailey and Davies 1987), particularly during spring (Santos 2009). Snakes escaping also precluded marking for individual identification; however, we assumed that the large size of the study area minimized pseudoreplication of the data by recording the escape distance of the same individual several times.

Statistical analyses

Two analyses were conducted to identify which factors could explain the occurrence of each dorsal pattern of the viperine snake across its range: 1) with the entire data set (all localities, including museum specimens and live ones), we performed a generalized linear model (GLZ) with a binomial distribution of the dependent variable and a logit function, using dorsal pattern (zigzag and bilineata) as a dependent variable, and sex and body size of individuals (SVL) as factors. This analysis was aimed at examining the effects of factors like sex and body size on the dorsal pattern occurrence regardless of locality. 2) For each population, we used the percentage of snakes with a bilineata pattern as a dependent variable (percentages were arc-sin transformed), and we performed a general linear model (GLM) with the habitat type (open and closed habitats), the distance (km) to the closest V. latastei population, and the latitude, as factors. This second analysis assessed how landscape type, model organism proximity (viper), and geographic location affected the dorsal pattern frequency of the bilineata pattern among the populations. To meet the assumption for a linear model, residuals of the dependent variable on the grouping factors were examined.

Log-linear analyses were run to test differences in predation (depredated and non-depredated) among the 3 plasticine models (zigzag, bilineata, and control) and the 2 study sites (Charca Suárez and Padul peat bog). We tested the association among these 3 categorical variables (predation, model type, and locality) in a multidimensional contingency table. Log-linear analysis uses a likelihood ratio chi-square statistic. The algorithm used generates several models to test interactions among all variables and selects the least complex model that best accounts for the variance in the observed frequencies. We were specifically interested in testing the association of Predation× Locality and Predation×Dorsal pattern. The results were interpreted by checking the odds-ratio scores in expected values of partial and marginal association tests. The odds ratio computes the likelihood ratio statistic of the model containing a particular term or not (significant when P < 0.05).

A GLZ with a Poisson distribution was used to compare fleeing distance of snakes of both dorsal patterns when confronted with a potential predator. Fleeing distance is a continuous variable; however, due to inaccuracies of some distance measurements we considered fleeing distances in 0.5 m intervals, using a Poisson distribution to perform the model. GLZs were performed with air temperature (surrogate of body temperature) and body size as covariates, due to their potential effect on escape speed of this species. Body condition (BC) can be considered a surrogate of snake’s locomotor capacity (Jayne and Bennett 1990). Thus, we hypothesized that differences in BC between zigzag and bilineata could be linked to variation in escape distance to predators, that is, lower BC for snakes that flee at longer distances, and higher BC for snakes that allow a closer approach of predator. The residuals of the regression analysis between log body mass and log SVL were used as a measure of BC. We were unable to capture most of the individuals during our trials to check differences in escape behavior between zigzag and bilineata snakes in the face of an approaching predator. For this reason, BC was compared between zigzag and bilineata snakes from the same locality, the Ebro Delta, where body mass and length were measured for a large sample size. BC was calculated separately for males and females (Santos and Llorente 2004). In the “Results” section, means are followed by ±1 SD.

Results

Of the 917 snakes examined, 687 were zigzag and 213 bilineata (76% and 24%, respectively). However, there was interpopulation variation in the frequency of the 2 dorsal patterns across the N. maura range (Figure 1A). The distribution of both dorsal patterns was puzzling; some sites were devoid of bilineata individuals, and bilineata frequencies differed markedly between pairs of nearby sites. The bilineata dorsal pattern was present in regions with the 3 phylogenetic groups (Figure 1A and Table 1).

The GLZ demonstrated that neither SVL (Wald Statistic = 0.27, P = 0.6) nor sex (Wald Statistic = 0.09, P = 0.8) explained the occurrence of the bilineata pattern on snakes. We repeated this analysis with only animals examined at sites where at least one bilineata snake appeared, and the results were the same, that is, there were no sexual or ontogenetic effects in dorsal pattern occurrence. The GLM indicated that the percentage of individuals with bilineata pattern at each of the 12 populations was related to the habitat type (F = 14.45, P = 0.005); bilineata snakes were frequent only in open habitats, whereas this pattern was almost absent from closed habitats (Figure 1A). A nonparametric test using the percentage of occurrence in each population confirmed significant differences between the 2 major habitat types in the frequency of the 2 dorsal patterns (Mann–Witney U-test, Z = 2.86, P = 0.004). By contrast, we found no significant spatial effect, that is, geographic latitude of each population (F = 0.20, P = 0.7), nor any relation with the distance to the closest V. latastei population, that is, the model species for the mimicry (F = 0.77, P = 0.4).

From the 380 predation trials made with plasticine models (276 in Charca Suárez and 104 in Padul), 100 attacks by avian predators were recorded (26.3% of the models). In total, bilineata was the most frequently attacked model (n = 45; 35% of the bilineata models examined), the control the second most attacked (n = 32; 25% of the control models examined), and the zigzag pattern was the least attacked (n = 23; 19% of the zigzag models examined). The log linear analysis showed a significant interaction between Locality and Attacks, and between Model and Attacks (Table 2). The marginal frequency tables showed that the occurrence of attacks was higher in Padul than Charca Suárez, and for bilineata than for the other 2 models. When mammal attacks were examined, we found no differences in the frequency of attacks made on the 3 plasticine models (χ2 = 1.77, P = 0.18).

Table 2.

Partial and marginal association tests of the log linear analysis for the attack condition (attacked and not attacked by birds), locality (Padul and Charca Suarez) and dorsal pattern (bilineata, zigzag, and control) of plasticine models of Natrix maura in 2 localities of south eastern Iberian peninsula where both natural morphotypes for dorsal pattern (bilineata, zigzag) are represented

	df	Partial X2	P	Marginal X2	P
Site (S)	1	79.40	<0.0001	79.40	<0.0001
Dorsal pattern (DP)	2	0.08	0.96	0.08	0.96
Attacks (A)	1	87.28	<0.0001	87.28	<0.0001
S * DP	2	0.45	0.80	0.22	0.90
S * A	1	35.88	<0.0001	35.65	<0.0001
DP * A	2	9.11	0.01	8.89	0.01

Note: Sample size in the results section.

Fleeing distances were recorded for eight bilineata (average flight distance 3.0 ± 1.1 m), and 17 zigzag N. maura individuals (average flight distance 0.9 ± 0.7 m). The GLZ showed that bilineata individuals fled at the farther distances from the approaching potential predator than did zigzag individuals, irrespective of body size or air temperature (Table 3). It bears noting that 6 out 17 zigzag snakes recorded did not flee at all, and exhibited aposematic behavior such as head triangulation, coiling, snorting, and false attacks, only after they were gently captured by hand. The remaining 11 individuals that fled to a very close distance from the researcher did not exhibit aposematic behavior.

Table 3.

General linear models of differences in flight distance between bilineata and zigzag Natrix maura phenotypes using snake body size (SVL, snout–vent length, in millimeter) and air temperature (°C) as covariates in the population of Charca Suárez, south-eastern Iberian peninsula

	Estimate	Standard error	Wald Stat	P
SVL	−0.000874	0.002	0.25682	0.61
Air temperature	−0.006432	0.007	0.79421	0.37
Dorsal pattern	0.601609	0.164	13.32934	<0.0001

No differences were found in BC between zigzag and bilineata individuals either in males (Student’s t = 0.18, df = 67, P = 0.9) or in females (Student’s t = 0.10, df = 53, P = 0.9). This suggests that body shape of zigzag and bilineata individuals did not vary.

Discussion

Dorsal pattern and geography

In the interpretation of animal coloration, tests of geographical variation are critical (Cox and Davis Rabosky 2013; Santos et al. 2014). Our results on the geographic distribution of the 2 dorsal patterns of the viperine snake were puzzling, although in contrast to Schätti (1982), we found that bilineata individuals were present in many inland regions. Notably, bilineata individuals also appeared in populations representing the 3 phylogenetic lineages of the species, north-western Africa east of the Muluya River, north-western Africa west of the Muluya River, and the European clade (Guicking et al. 2002, 2008). This suggests the character to be African and ancient, as the species evolved in that continent and clades separated 3.5–4.0 Myr (Guicking et al. 2002)—although it may also have evolved repeatedly (by convergence) in response to similar ecological demands in the different clades, as is common in snakes (Allen et al. 2013).

Dorsal pattern and habitat

The polymorphism for the dorsal pattern in N. maura seems to be linked to the habitat. Regardless of sex and body size, bilineata individuals were almost exclusively in open habitats, whereas zigzag individuals were found in all types of habitats (see also Wolf and Werner 1994; Kark et al. 1997; Cox and Davis Rabosky 2013), although the most frequent in closed habitats. This result reveals why previous studies failed in fitting a geographic pattern for the presence of both dorsal patterns in N. maura (Schätti 1982; Duguy and Saint-Girons 1993); mimicry is not restricted to one sex or to ontogenetic stage, and variation occurred across habitat, rather than geography. In their review of the frequency of the striped dorsal pattern in snakes, Wolf and Werner (1994) concluded that snakes bearing this phenotype inhabited habitats that were less structurally complex than those bearing blocked or cross-banded dorsal patterns (see also Pough 1976).

Camouflage is frequently selected in snakes (Allen et al. 2013), with different habitats favoring different morphs (Kark et al. 1997; Shine and Harlow 1998). We speculate that the zigzag dorsal pattern would be beneficial in concealing individuals in more closed habitats, where a chaos of images is created by lights and shadows (Cuthill et al. 2017; Endler and Mappes 2017). While not proven, this interpretation has been posed several times for this dorsal pattern in vipers (Shine and Madsen 1994; Niskanen and Mappes 2005). In contrast, striped snakes would be favored in more open habitats, where the combination of lights and shadows would create straighter forms on the background (Wolf and Werner 1994). As N. maura inhabits 2 highly contrasting aquatic habitats, wetlands and rivers (Santos 2009), we propose here that heterogeneity in the vegetation structure surrounding water bodies would be the driving factor for the emergence of dorsal pattern dimorphism in this species.

Dorsal pattern and predator attacks

Our experiment with plasticine models within the same locality and habitat confirmed that predators responded differentially to the snake dorsal pattern. On a grassland background, control and bilineata models were more frequently attacked by aerial, visually oriented predators (birds), than were zigzag models. This was the same result found with viper plasticine models (Valkonen et al. 2011b), confirming that snakes with the zigzag dorsal pattern better avoided attacks than did snakes with other dorsal patterns (Wüster et al. 2004; Niskanen and Mappes 2005), at least when prey are static (Lindell and Forsman 1996). Our result supports the idea that the camouflage provided by the zigzag pattern or that predator wariness to the model organisms (i.e., vipers) can drive a mimetic system in snakes (Greene and McDiarmid 1981; Cox and Davis Rabosky 2013).

In the Padul peat bog, plasticine replicas were more frequently attacked than those in Charca Suárez. Although the potential predator communities were similar, we were unable to adequately measure predator density at either site, and therefore we have no suitable explanation for this difference.

Dorsal pattern and snake behavior

One of the benefits of the high and constant body temperature that N. maura maintains during its K thermoregulation is to retain a potential for instant high muscular performance (Hailey and Davies 1986). Our study demonstrated that, when a potential predator approached, bilineata individuals fled from longer distances than did zigzag individuals. Some of the latter even remained static until they were captured by the researcher. Thus, viperine snakes with different dorsal patterns varied their antipredatory strategies, suggesting adaptations to increase their fitness. Bilineata individuals, lacking the camouflage or the aposematic effect of the zigzag pattern, begin to flee when a potential predator approaches. The dorsal stripe design is common in snakes, especially in the fastest species (Brattstrom 1955). Snakes bearing this dorsal pattern rely on the confounding effect on predators of the dorsal stripes when the snake moves (Pough 1976; Allen et al. 2013), making the formation of a search image more difficult (Hogan et al. 2016). That is, the perception of speed (von Helversen et al. 2013) and direction (Jackson et al. 1976) is disrupted, and the location of the focal part of the body, the head (Langkilde et al. 2004), is more difficult (Brodie 1992; Hu et al. 2009). This interpretation fits the link between the bilineata dorsal pattern and open habitats found in this study. In open habitats such as wetlands, the vegetation has less structural complexity and less irregularly shaped patterns of shadowing. In these micro-scenarios, we argue that the zigzag dorsal patterned viperine snakes would begin to leave evolutionary room for striped snakes, which adopt longer escape distances from approaching predators (Allen et al. 2013).

Crypsis or aposematism?

Experiments with plasticine models did not resolve the question that differences in predation rates between morphotypes were the consequence of the cryptic or aposematic role of the zigzag pattern, or even whether both roles are acting. Several studies have demonstrated the aposematic value of the zigzag pattern (Wüster et al. 2004; Niskanen and Mappes 2005; Valkonen et al. 2011b). However, simultaneous adaptive benefits of this pattern cannot be ruled out, for example, a distance-dependent function for crypsis (Valkonen et al. 2011b). We suggest 3 complementary reasons to explain the cryptic role of the zigzag pattern in the viperine snake system: 1) the frequency of zigzag snakes is independent of the distance to the aposematic viper model, taking into account that the aposematic value of a design depends on the presence of the model (Pfennig et al. 2001). 2) The maintenance of the dorsal polymorphism in N. maura is habitat-dependent, suggesting some type of concealment value for the zigzag pattern in more structured habitats. This habitat dependence is supported by the study of Santos et al. (2014), who verified that the zigzag design in V. latastei varied geographically in relation to the lithology and land uses. 3) The observed differences in escape behavior and approach distance between zigzag and bilineata viperine snakes. Selection from visually guided predators acts on the interaction between dorsal pattern and behavior (Pough 1976; Lindell and Forsman 1996). The fact that most zigzag viperine snakes remain still before the predator (the researcher) begins any prey-catching behavior, and that only after handling displayed indisputable aposematic behaviors, such as head triangulation, snorting, and (false) attacks (authors, personal observation), suggests that the snakes firstly rely on the cryptic role of the zigzag dorsal pattern. If this camouflage strategy fails, the aforementioned signaling mechanisms trigger, and the aposematic role of the zigzag dorsal pattern appears. Indeed, we conclude that the zigzag dorsal pattern in the viperine snake (as well as in vipers, the model organism) would have a double function, concealing individuals against the background and deterring predators. We also experimentally confirm the Batesian mimesis of the viperine snake to the European vipers as model organisms.

We acknowledge that our results on attacks are based on plasticine models, and those on fleeing distance of live individuals are derived from a rather small sample size, which makes it advisable to discuss these issues under alternative hypotheses. The Batesian mimesis of the viperine snakes on vipers do not limit to the zigzag dorsal pattern, and includes other behavioral aposematic traits such as coiling, striking, snorting, and head triangulation (Rollinat 1934; Aubret and Mangin 2014). However, as we used the same cast (that was not coiled neither showed head triangulation) for the 3 types of models (Figure 2), we deduce that the higher frequency of bird attacks to bilineata and control models was exclusively the consequence of the differences in dorsal design among models. The use of static behavior, typically balling up, as defence in N. maura is more frequent in small individuals (Hailey and Davies 1986). However, we rule out the confounding effect of body size in zigzag and bilineata dorsal patterned individuals when a potential predator approaches, as we failed to find body-size differences between the individuals of different dorsal patterns recorded for this analysis. The same null hypothesis was confirmed for the effect of sex, as sexual dichromatism is often absent in snakes (Shine and Madsen 1994). Lastly, in their laboratory experiment, Hailey and Davies (1986) also found that starved snakes (for 4–5 weeks) were more likely to use static defence than post-absorptive ones. We were unable to consider this factor in our trials, as most bilineata individuals fled to the water, so that the weight or feeding status of these specimens could not be checked. However, we minimized this possibility, as the study site harbored a good population of birds and N. maura, likely for its abundant trophic resources (unpublished data of the authors). Moreover, during field surveys over 3 years, we never found an emaciated individual of this species in the study area (n = 60).

Ecological significance of dorsal polymorphism

Assuming the ecological advantages of the zigzag pattern in N. maura, what are the evolutionary drivers of the maintenance of dorsal polymorphism in this species? In this study, we have ruled out the possibility that a historical (phylogeographic) factor could give rise to dorsal variability. Moreover, the result for the distance to the closest viper population also suggests no effect of model-species proximity on the frequency of the viperine snake dorsal pattern. By contrast, we have demonstrated that the coexistence of both designs is habitat dependent, suggesting that snake populations generate individuals of both patterns, and these would increase their fitness (survival) depending on the habitat structure (Pyron and Burbrink 2009). According to Levins (1968), phenotypic polymorphism increases fitness when a species occurs in heterogeneous environments. The landscape heterogeneity of the Mediterranean region, in some places composed of a patchy distribution of freshwater environments surrounded by open or forested habitats, would have been the scenario to allow coexistence of alternative N. maura dorsal patterns adapted to different antipredatory strategies.

Author contributions

X.S. and J.M.P. conceived the idea. X.S. analyzed the data, and all authors contributed to collecting the data and writing the manuscript.