Phylogeography supports lineage divergence for an endemic rattlesnake (Crotalus ravus) of the Neotropical montane forest in the Trans-Mexican Volcanic Belt

Abstract

The formation of the Trans-Mexican Volcanic Belt (TMVB) and Pleistocene climatic fluctuations have been shown to influence the diversification of lineages and species distributed throughout central Mexico. In some taxa, however, evidence of lineage diversification is not easily recognized, as often is the case in reptiles. Here we present a phylogeographic study on a Mexican endemic rattlesnake species (Crotalus ravus), with the aim of understanding how distinct lineages are distributed across the TMVB. Genetic (mtDNA) and genomic (ddRADseq) data were generated from samples across the species’ range to evaluate phylogeographic structure, estimate phylogenetic relationships and divergence times, and perform environmental niche modeling (ENM). Both datasets recover strong phylogeographic structuring of two distinct lineages on an east-west axis, with an estimated Pleistocene divergence (~1.47 Myr). The ENM suggest that the distribution of the two lineages experienced expansion and reduction events throughout recent evolutionary time. We attribute the diversification of C. ravus lineages to geological events associated with the formation of the TMVB, as well as Quaternary climate changes, both of which have been previously recognized in co-distributed taxa in the TMVB. This work emphasizes the existence of cryptic diversification processes in a morphologically conserved species distributed in a region of complex climatic and orogenic heterogeneity.

endemism, mtDNA, phylogepography, snake, SNPs, TMVB

INTRODUCTION

Rattlesnakes are distributed across the American continent with about 55 species (Uetz et al., 2020), and Mexico has the greatest diversity (~45 spp.) (Campbell & Lamar, 2004; Uetz et al., 2020). Mountain highland species are the most diverse (~30 spp.; Campbell & Lamar, 2004; Uetz et al., 2020) and can be considered the centre of diversity for rattlesnakes (Bryson et al., 2011c; Blair et al., 2016). Traits such as philopatry, low vagility, conserved morphology and microhabitat selection have promoted studies at the inter- and intraspecific levels in venomous snakes (Bryson et al., 2011a, c; Blair & Sánchez-Ramírez, 2016; Blair et al., 2018), with few studies at the infraspecific level (Sunny et al., 2018).

Past geological and topographic events such as volcanism and mountain formation have largely influenced allopatric speciation processes in reptiles and particularly rattlesnakes (Bryson et al., 2011b, c, 2012a, b; Bryson & Riddle, 2012; Blair et al., 2018). Climate fluctuations in the Trans-Mexican Volcanic Belt (TMVB) region have had a direct influence on geographical distribution patterns across taxa over time (Mastretta-Yanes et al., 2015). These environmental fluctuations likely exerted a major influence on organisms with narrow thermal requirements such as reptiles, which are sensitive to climate change (Huey, 1982; Russildi et al., 2016). During the Pleistocene, strong oscillations between glacial and interglacial periods directly influenced the montane forests of central Mexico (Ramírez-Barahona & Eguiarte, 2013). Glacial periods caused the downward movement of montane forests, allowing formerly isolated populations of lowland temperate-affinity taxa to reconnect by dispersion (McDonald, 1993; Marshall & Liebherr, 2000; Flores-Villela & Goyenechea, 2001), while interglacial periods interrupted the genetic flow between these populations (Anducho-Reyes et al., 2008; Mastretta-Yanes et al., 2015), in some cases promoting parapatric speciation processes in some taxa (Mastretta-Yanes et al., 2015). Two main models for the influence of climate on tropical mountain species were proposed regarding the existence of pleistocenic refugia during the Last Glacial Maximum (LGM). The dry forest model (Haffer, 1969) states that decreasing temperatures during the glacial periods favoured the expansion of montane lineages, but the aridity of the lowlands that characterized these periods of cooling limited them to areas of intermediate elevation in the mountains, where temperature and humidity were stable (Haffer, 1969; Burnham & Graham, 1999; Ramírez-Barahona & Eguiarte, 2013). In contrast, the moist forest model proposes that the decrease in temperature was not necessarily accompanied by aridity from lack of rainfall, and that in the tropical highlands this factor had little or no effect on species distributions during the climate fluctuations of the Pleistocene (Colinvaux et al., 2000; Bush & Silman, 2004; Ramírez-Barahona & Eguiarte, 2013). Bryson et al. (2011b) found that climate oscillations during the Pleistocene were the main driver of species divergence in the Crotalus intermedius mountain rattlesnake complex, allowing for the dispersal of these snakes throughout the mountains in cold periods (i.e. dry forest model). Bryson et al. (2011a) and Leaché et al. (2013) concluded that the genetic structure of gopher snakes (Pituophis) and fence lizards (Sceloporus bicanthalis) is partially explained by past climate fluctuations, which together with uplift in the TMVB formation, promoted population diversification and isolation.

Crotalus ravus (Cope, 1865) is an endemic mountain rattlesnake distributed across the TMVB with a distinctive arrangement of head scales and with extensive geographic colour variation (Campbell & Armstrong, 1979). Crotalus ravus inhabits pine-oak forests and grasslands at elevations of 1500 to 3000 m.a.s.l. (Campbell & Lamar, 2004). This rattlesnake belongs to the group of mountain pygmy rattlesnakes known as the ‘Crotalus triseriatus’ group (Bryson et al., 2011c; Alencar et al., 2016; Blair et al., 2016) and has recently been established as a sister species to Crotalus brunneus (Blair et al., 2018), which is distributed in the mountains of the state of Oaxaca and was considered its subspecies until recently.

We used mitochondrial (ND4 and 12S), nuclear DNA (SNPs) and ecological niche modelling to test whether the genetic and ecological divergence of the C. ravus complex was shaped by past volcanic activities and climate fluctuations during the Quaternary in the TMVB. We hypothesize that geological and volcanic activity and Quaternary climate change in the TMVB promoted diversification in two ways: either by (1) generating new phylogeographical barriers leading to allopatric speciation; or (2) past climate fluctuations creating new habitats and facilitating dispersal from adjacent areas and the subsequent divergence of taxa.

MATERIAL AND METHODS

DNA extraction and mitochondrial sequencing and processing

We collected a total of 37 tissue samples of the mountain pygmy rattlesnake C. ravus across its range (Table 1; Fig. 1). Genomic DNA was extracted following the DNeasy Blood and Tissue Kit protocol from Qiagen (Qiagen Valencia, CA). We sequenced 1033 base pair fragments from two mitochondrial loci NADH dehydrogenase 4 and flanking tRNAs (ND4) and the 12S ribosomal subunits (12S) using the primers proposed in Arevalo et al. (1994). The protocol for DNA extraction, amplification and purification of each gene considered, is described in Appendix S1 in the Supporting Information.

Table 1.

List of individuals of C. ravus considered for this study. The bold letters represent the sequences downloaded from GenBank to complete the representativeness of the sampling

ID	Country/state	Locality	Latitude	Longitude	Elevation (m)	Voucher	GenBank number (ND4)	GenBank number (12S)	Population
MEX_Joq_West_6	Mexico, Estado de México	Joquicingo de León Guzmán	19.0421	-99.5317	2619	AYCB 39	MN527363	MN527037	Western
GUE_Tax_West_7	Mexico, Guerrero	Taxco	18.5609	-99.6110	2040	MZFC 25112	MN527364	MN527038	Western
MEX_Joq_West_8	Mexico, Estado de México	Joquicingo de León Guzmán	19.0392	-99.5319	2653	AYCB 40	MN527365	MN527039	Western
CMX_Mag_West_9	Mexico, Mexico City	Magdalena Contreras	19.2879	-99.2670	2711	MZFC 20902	MN527366	MN527040	Eastern
MEX_Joq_West_10	Mexico, Estado de México	Joquicingo de León Guzmán	19.0392	-99.5319	2660	AYCB 41	MN527367	MN527041	Western
PUE_Tep_East_11	Mexico, Puebla	Tepanco de López	18.5504	-97.5601	1816	ROM 47070	HQ257813	HQ257567	Eastern
PUE_Zac_East_12	Mexico, Puebla	Zacatepec	19.2636	-97.5326	2372	ROMFC228	HQ257881	HQ257637	Eastern
PUE_Toc_East_13	Mexico, Puebla	Tochimilco	18.8887	-98.5724	2051	ROM 47048	HQ257824	HQ257578	Eastern
VER_Lim_East_14	Mexico, Veracruz	El Limón Totalco	19.5002	-97.3448	2351	ROM 47001	HQ257843	HQ257597	Eastern
TLA_Hua_East_15	Mexico, Tlaxcala	Huamantla	19.3139	-97.9233	2511	ROM 47066	HQ257838	HQ257592	Eastern
TLA_Hua_East_16	Mexico, Tlaxcala	Huamantla	19.4402	-97.8825	2598	MZFC 26315	MN527369	MN527042	Eastern
TLA_Hua_East_17	Mexico, Tlaxcala	Huamantla	19.4402	-97.8825	2601	MZFC 26318	MN527370	MN527043	Eastern
PUE_Chi_East_18	Mexico, Puebla	Chignahuapan	19.7448	-98.0314	2708	MZFC 4797	MN527371	MN527044	Eastern
TLA_Atl_East_19	Mexico, Tlaxcala	Atlangatepec	19.5296	-98.2060	2514	MZFC 24360	MN527372	MN527045	Eastern
PUE_Zap_East_20	Mexico, Puebla	Zapotitlán de Salinas	18.3273	-97.4753	1476	Cravus_E8	MN527373	N/A	Eastern
CMX_Xoc_West_21	Mexico, Mexico City	Xochimilco	19.2402	-99.1282	2350	MZFC 25244	MN527374	MN527046	Eastern
PUE_Ori_East_22	Mexico, Puebla	Puebla	19.0244	-98.1638	2178	JAFF 2235	MN527375	N/A	Eastern
PUE_Ori_East_23	Mexico, Puebla	Puebla	19.0154	-98.1263	2282	Cravus5	MN527376	N/A	Eastern
CMX_Hui_West_28	Mexico, Mexico City	Huixquilucan de Degollado	19.3813	-99.3578	2777	MZFC 14287	MN527377	MN527048	Eastern
MEX_Ame_East_29	Mexico, Estado de México	Amecameca	19.1223	-98.7667	2482	MZFC 6286	MN527378	MN527049	Eastern
PUE_Zap_East_30	Mexico, Puebla	Zapotitlán de Salinas	18.3510	-97.4435	1599	MZFC 16469	MN527379	MN527050	Eastern
PUE_Zap_East_31	Mexico, Puebla	Zapotitlán de Salinas	18.3320	-97.5134	1580	JAC 22511	MN527380	MN527051	Eastern
TLA_Pan_East_35	Mexico, Tlaxcala	Huiloapan	19.3750	-98.2509	2331	ADM 2	MN527384	MN527055	Eastern
PUE_Tla_East_36	Mexico, Puebla	Tlatlauquitepec	19.8965	-97.5194	2181	MZFC 28471	MN527385	MN527056	Eastern
VER_Teo_East_37	Mexico, Veracruz	El Limón Totalco	19.4814	-97.3507	2408	MZFC 30316	MN527386	MN527057	Eastern
HID_Tep_East_38	Mexico, Hidalgo	Tepeapulco	19.7805	-98.5528	2548	MZFC 34494	MN527387	MN527058	Eastern
MEX_Mex_West_39	Mexico, Estado de México	San José El Totoc	18.9791	-99.3655	2499	MZFC 34495	MN527388	MN527059	Western
TLA_Zit_East_40	Mexico, Tlaxcala	Zitlaltepec	19.2283	-97.9129	2588	MZFC 34496	MN527089	MN527060	Eastern
CMX_Cec_West_41	Mexico, Mexico City	Santa Cecilia Tepetlapa	19.2142	-99.1045	2495	MZFC 34497	MN527390	MN527061	Eastern
MEX_Joq_West_42	Mexico, Estado de México	Joquicingo de León Guzmán	19.0476	-99.5453	2754	MZFC 34498	MN527391	MN527062	Western
PUE_Tec_East_43	Mexico, Puebla	Santa María Zotoltepec	19.6914	-97.9141	2441	MZFC 34499	MN527392	MN527063	Eastern
PUE_Zau_East_45	Mexico, Puebla	Santiago Zautla	19.7110	-97.6200	2474	MZFC 34501	MN527394	MN527065	Eastern
HID_Tep_East_46	Mexico, Hidalgo	Tepeapulco	19.7805	-98.5528	2548	N/A	MN527395	MN527066	Eastern
MEX_Mex_West_47	Mexico, Estado de México	San José El Totoc	18.9791	-99.3655	2499	MZFC 34514	MN527396	MN527067	Western
PUE_Zot_East_48	Mexico, Puebla	Santa María Zotoltepec	19.6757	-97.8670	2268	MZFC 34500	MN527397	MN527068	Eastern
PUE_Que_East_49	Mexico, Puebla	Quecholac	18.9633	-97.6261	2308	MZFC 34502	MN527398	MN527069	Eastern
PUE_Ixt_East_50	Mexico, Puebla	Vista Hermosa	19.7091	-97.8216	2853	MZFC 34515	MN527399	MN527070	Eastern
PUE_Rit_East_51	Mexico, Puebla	Santa Rita Tlahuapan	19.3568	-98.5936	2670	MZFC 34503	MN527400	MN527071	Eastern
MEX_Joq_West_52	Mexico, Estado de Mexico	Joquicingo de León Guzmán	19.0475	-99.5479	2763	MZFC 34504	MN527401	MN527072	Western
MEX_Joq_West_53	Mexico, Estado de México	Joquicingo de León Guzmán	19.0475	-99.5479	2759	MZFC 34505	MN527402	MN527073	Western
MEX_Zac_ West _55	Mexico, Estado de México	Zacamulpa	19.3560	-99.3372	2684	ROM 47050	HQ257823	HQ257577	Eastern
MEX_Joq_West_56	Mexico, Estado de México	Joquicingo de León Guzmán	19.0376	-99.4992	2602	MZFC 25110	MN527403	MN527074	Western
PUE_Nop_East_61	Mexico, Puebla	Santa María Ixtiyucan	19.1712	-97.7761	2358	CHJ 1401	MN527404	MN527075	Eastern
PUE_Nop_East_62	Mexico, Puebla	Santa María Ixtiyucan	19.1712	-97.7761	2361	CHJ 1394	MN527405	MN527076	Eastern

Figure 1.

Map showing the geographic distribution of samples of C. ravus included in this study. The star represents Mexico City. Haplotype network obtained from mitochondrial DNA data. The eastern population is in blue and the western population is in grey. The numbers on the network are each individual’s ID.

Forward and reverse sequences were assembled by eye in SEQUENCHER v.5.4.6. (http://genecodes.com/ 2016). We downloaded twelve sequences of C. ravus available on GenBank, all of which were aligned and edited manually using PhyDE v.0.9971 (http://phyde.de/download.html 2010) (Table 1). For all mtDNA analyses we used the concatenated matrix, except for genetic diversity.

Genetic diversity and haplotype networks

We obtained the haplotypes (H), genetic (h) and nucleotide (π) diversity for each locus and for the concatenated matrix in DNASP v.5.10 (Librado & Rozas, 2009). We estimated the population differentiation (FST) and the hierarchical percentages of molecular variation using an AMOVA in ARLEQUIN v.3.5 (Excoffier & Lischer, 2010), between the eastern and western populations of C. ravus based on results of the analyses described below. Genealogical relationships among haplotypes were calculated using statistical parsimony implemented in TCS v.2.1 (Clement et al., 2000). We generated one network for the concatenated matrix. A 95% connectivity limit was set, and gaps were treated as missing data.

Demographic history and SPLITSTREE

We calculated Tajima’s D (Tajima, 1989) and Fu’s Fs (Fu, 1997) for the population expansion test for two sets of C. ravus populations (east-west, see Results).

A Bayesian skyline plot (BSP; Drummond et al., 2005) of changes in effective population size (N_e) over time was generated with BEAST v.2.4.2 (Bouckaert et al., 2014), using only the ND4 locus because its substitution rate is known. Analysis specification is presented in the Supporting Information (Appendix S2).

We performed a phylogenetic network analysis using the Neighbour-net method implemented in SPLITSTREE v.4.13 (Hudson & Bryant, 2006). A distance-based method for the construction of phylogenetic networks, Neighbour-net is based on the neighbour-joining algorithm (a network or graph allows for multiple paths). The Neighbour-net graph was calculated based on the ND4 and 12S markers using a Neighbour-net algorithm with uncorrected distances.

Phylogenetic relationships and divergence times

We estimated the phylogenetic relationships among mtDNA sequences of C. ravus and the C. triseriatus group. Sequences of C. brunneus (N = 5), Crotalus exiguus (N = 4), C. triseriatus (N = 6), Crotalus aquilus (N = 6), Crotalus lepidus (N = 2), Crotalus lepidus klauberi (N = 2), Sistrurus miliarius (N = 1) and Sistrurus catenatus (N = 1) were downloaded from GenBank and used as outgroups. We used JMODELTEST v.2.1.10 to find the best-fit model of substitution for mtDNA according to AIC criteria with 95% confidence credibility (Darriba et al., 2012). We used MRBAYES v.3.2.7 (Huelsenbeck & Ronquist, 2001) for over 50 million generations with a sample frequency of 1000 steps, and a 25% burn-in. Convergence of the four chains was evaluated in TRACER v.1.7.1 (Rambaut et al., 2018). The results were visualized in FIGTREE v.1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/).

We constructed a calibrated phylogenetic tree under a coalescent approach in BEAST v.2.4.2 (Bouckaert et al., 2014) using the concatenated mtDNA matrix and all outgroup sequences. Analysis specifications can be found in the Supporting Information (Appendix S3).

RADSeq library generation and computational analysis

Loci were obtained from the genome of 27 individuals of C. ravus through the ddRADSeq process (Peterson et al., 2012). More details of the library preparation and bioinformatic analysis are given in the Supporting Information (Appendix S4).

Stacks Pipeline and raw reads analysis

We processed the ddRADseq data largely following the pipeline described in Streicher et al. (2014) with a few modifications described below. We used STACKS v.1.8.4 (with default parameters except where noted) to process the raw data, first using process_radtags to remove low-quality reads (removing reads whose average quality was 90% or less). Second, we used ustacks to generate loci (stacks) for each individual, requiring a minimum of five reads per stack. We then built a catalogue of loci using cstacks, allowing up to four mismatches per locus between individuals. We matched sequences from individuals to the catalogue using sstacks and then produced data sets with either 50% or 30% missing data for downstream analyses using the program populations.

Genetic diversity of the SNPs data

We estimated nucleotide (π) diversity within the eastern and western populations sets of C. ravus individuals in ARLEQUIN v.3.5 (Excoffier & Lischer, 2010). An AMOVA was run to see the percent variation in the east-west configuration of populations and population differentiation was calculated (FST). We applied two configurations to this analysis, first for individuals with shared lineages (see results below) in the eastern population and then those attached in the west.

Genetic structure and phylogenetic network using SNPs

The SNPs of the 27 individuals were genetically analysed. We ran STRUCTURE v.2.3.4 (Pritchard et al., 2000) under the admixture model with correlated allele frequencies. Ten independent chains were run for each K, from K = 1 to K = 10. Burn-in length was 500 000 and the number of MCMC replications after the burn-in was 1 000 000. The most likely number of populations was determined estimating the DeltaK (∆K) and the log likelihood of K, ln P(K) = L(K) between successive K values (Evanno et al., 2005).

Additionally, we used the conStruct R package (Bradburd et al., 2018) to visualize the individuals clustering under a spatial model taking into account the combination of discrete and continuous genetic variation. The specifications of the analysis followed the protocol of Marshall et al. (2021), only the value of the K model was modified, which was from 1 to 3 following Structure results.

Finally, we used SPLITSTREE v.4.16.1 (Hudson & Bryant, 2006) to generate an unrooted phylogenetic network to reveal all possible relationships among individuals. The Neighbour-net graph was calculated based on 899 SNPs using a Neighbour-net algorithm with uncorrected distances.

Ecological niche modelling and niche comparison

We generated niche models to evaluate current and past climatic suitability for C. ravus in Mexico using five bioclimatic variables from WorldClim (bio4, bio9, bio15, bio18 and bio19) after eliminating most collinear variables using the variance inflation factor (VIF: Marquaridt, 1970). We ran eight model algorithms with the sdm R package (Naimi & Araújo, 2016): MaxEnt, MARS, GBM, RF, CART, SVM, GLM and GAM. We calibrated these models using 5-fold cross-validation and ten bootstrapping replications for a total of 50 repetitions per algorithm with a 70–30 data-partitioning scheme. We randomly created a set of 800 pseudo-absences (100 times the number of presences) to validate the models, and finally, we generated an ensemble model weighting each algorithm by maximizing True Skill Statistic (TSS). We transferred models to 11 past climate change scenarios available in the PaleoClim database (Brown et al., 2018), encompassing the following time horizons: Holocene: Meghalayan (4.2–0.3 kya), Northgrippian (8.3–4.2 kya), Greenlandian (11.7–8.3 kya), and Pleistocene: Younger_Dryas_Stadial (12.9–11.7 kya), Bølling-Allerød (14.7–12.9 kya), Heinrich_Stadial (17.0–14.7 kya), the Last Glacial Maximum (LGM ~21 kya) and the Inter-Glacial (LIG ~121 kya). We summed the suitability values for each clade across the region and tested suitability similarity across geography between both linages using Schoener’s D. We also used a Principal Components Analysis (PCA) to transform the environmental space of the 19 bioclimatic variables into a two-dimensional space (Strubbe et al., 2015) and tested niche similarity in ecological space by using Schoener’s D with a randomization test (100 replicates) in the ecospat R package (Di Cola et al., 2017).

RESULTS

Genetic diversity and haplotype networks

Forty-four mtDNA sequences of C. ravus were obtained, resulting in a matrix with a length of 1033 bp (Table 1). The ND4 locus was 711 bp long and the 12S locus was 372 bp long. The entire known distribution range of the species for individuals was represented by the data of the two mitochondrial markers (only three individuals had missing information re. the 12S locus; Fig. 1).

A total of 38 haplotypes were recovered in a single network using the concatenated mtDNA matrix (Fig. 1). The haplotype network recognized two main haplotype groups, one with individuals distributed in the east (13 haplotypes) and another with individuals distributed in the west (25 haplotypes). The two haplogroups were separated by seven mutational steps (from H7 to H26). Of all the haplotypes, four were shared by two or more individuals from the same population of C. ravus (eastern: H4, H12, H13; western: H31), while only the H21 haplotype was shared by an individual from the eastern population (46) and one from the western population (42); however, the latter seemed to be more closely related to the eastern population haplotypes. Finally, one haplotype collected in the east (17) was grouped with haplotypes found in the west (H36) (Fig. 1). Based on the results of the haplotype network with the mtDNA and those of the STRUCTURE analysis with the SNPs (see below), we report our results based on the two lineages (western and eastern).

Genetic diversity (h) and nucleotide diversity (π) based on the concatenated matrix mtDNA were high for the western populations (h = 0.90; π = 0.004) and low for eastern populations (h = 0.37; π = 0.003), and were moderate for both lineages analysed together (h = 0.66; π = 0.009), respectively (Table 2). The AMOVA analysis grouping in an east-west configuration (both groups) indicated that the highest percentage of genetic variation was explained between populations (61.63%), while within populations only 38.37% was explained (Table 3). The FST value was high and significant (0.62, P < 0.01) indicating genetic differentiation between eastern and western lineages. Values of Fu’s Fs and Tajima’s D for the eastern lineage of C. ravus were significant and negative (-3.98, P < 0.02; -2.12, P < 0.01, respectively), suggesting a recent demographic expansion, while for the western lineage values for both estimators were negative but not significant (Fu’s Fs = -0.19, P > 0.05; Tajima’s D = -1.03, P > 0.05) (Table 2). Bayesian skyline plots suggested a recent demographic expansion in the eastern lineage of C. ravus. Demographic expansion of the eastern lineage occurred ~3000 years BP with a narrow 95% confidence value (Fig. 2). The Neighbour-net tree recognized the eastern and western lineages, as did the parsimony network (see Supporting Information, Appendix S5). In fact, individual Hua_Este_17 collected in the eastern population was clustered into the western populations, while individual Joq_Oeste_42 collected in the western populations was clustered into the eastern populations.

Table 2.

Demographic indices of concatenated mtDNA (ND4 + 12S). Indices with an asterisk (*) indicate significance values less than 0.002

Population	N	Tajima’s D	Fu’s F	h	π
East	29	-2.11706*	-3.97994*	0.37 ± -0.01	0.003 ± -0.00
West	15	-1.03336	-0.19370	0.90 ± -0.00	0.004 ± -0.00
C. ravus	44	-0.087	-1.127	0.66 ± -0.00	0.009 ± -0.00

Table 3.

Results of analysis of molecular variation (AMOVA) of the populations of C. ravus from mtDNA (ND4 + 12S). Statistically significant values (< 0.05) are represented by an asterisk (*).

Source of variation	d.f.	Sum of squares	Variation (%)	Fixation index
Among populations	1	41.600	61.63
Within populations	43	54.000	38.37	FST = 0.62*
Total	44	95.600	100.00

Figure 2.

Bayesian skyline plot for the eastern C. ravus population. The y-axis is the effective population size and the x-axis is time (Mya).

Phylogenetic relationships, divergence times and demographic analysis

The phylogenetic relationships indicate C. ravus is a monophyletic group (PP = 1.0) (Fig. 3). Within the C. ravus clade, eastern individuals were recovered as a monophyletic group with a high support value (PP = 1.0), while western individuals appear as a paraphyletic group. The two lineages of C. ravus together appear as a sister clade to the Oaxaca rattlesnake C. brunneus according to Blair et al. (2018). Our time-calibrated phylogeny of the mtDNA matrix using BEAST revealed a Pliocene-Pleistocene history between C. ravus and C. brunneus, with a divergence time of ~2.46 Myr (4.2–0.7 Myr, 95% HPD). The two east-west lineages were recovered, though monophyly was forced, and in both groups divergence occurred ~1.47 Myr (2.7–0.4 Myr, 95% HPD) in the Pleistocene (Fig. 4).

Figure 3.

Bayesian tree obtained from the concatenated matrix of the ND4 and 12S loci. The blue bar indicates the population of the eastern lineage and the grey is the western lineage within the clade of C. ravus. The (*) above each branch denotes a 1.0 posterior probability value.

Figure 4.

Divergence times obtained from the coalescence analysis. The shaded bars represent the 95% confidence intervals for each calculated divergence age. The arrow indicates the divergence event between the eastern and western populations of C. ravus.

Genetic diversity from the SNPs data

We obtained two matrices with different numbers of missing data, one with 30% and the other with 50%, and we analysed 899 and 3194 SNPs, respectively. Although both matrices produced similar results, we only show the results for the matrix of 899 SNPs. The AMOVA revealed a genetic structure with a greater within-population percent variation (72.46%), while for mitochondrial DNA the fixation index was high (FST = 0.28, P < 0.05), indicating high genetic differentiation between eastern and western lineages (Table 4). Nucleotide diversity was higher in the eastern lineage than the western (Table 5). Similar results were found for the SNPs data set, even with 50% of the data missing (Supporting Information, Appendix S6).

Table 4.

Results of analysis of molecular variation (AMOVA) of the populations of C. ravus from SNPs data with 30% of missing data. Statistically significant values (> 0.05) are represented by an asterisk (*)

Source of variation	d.f.	Sum of squares	Variation (%)	Fixation index
Among populations	1	397.153	27.54
Within populations	52	1953.403	72.46	FST = 0.28*
Total	53	2350.556	100.00

Table 5.

Nucleotide diversity from SNPs data with 30% missing data. Statistically significant values (< 0.05) are represented by an asterisk (*)

Population	N	No. of loci	No. of loci used	No. of polymorphic loci	π
East	17	899	432	306	*0.172 ± 0.084
West	10	899	455	192	*0.153 ± 0.077

Genetic structure and phylogenetic network using SNPs

The results of STRUCTURE were used to infer population groups (K values) from our SNPs data set for C. ravus with 30% of the data missing. Structure analysis recognized two main groups, where the greatest increase in the DeltaK value was K = 2. These two recognized populations correspond to the mtDNA lineages, eastern and western populations throughout the TMVB. Two individuals were mixed (21 and 28, western population), specifically in the surrounding mountains close to Mexico City (Fig. 5B).

Figure 5.

Genetic structure obtained from the SNPs data. A, population structure analysis with the optimal K = 2 model estimated with conStruct. Each sample is represented by a pie-chart with colours indicating ancestral proportions. B, STRUCTURE analysis, the y-axis represents the ID of each individual and the x-axis is the percentage of lineage present in each snake sampled.

The results of the conStruct analysis show that the variation between the two C. ravus lineages could not be explained by a pattern of isolation by distance (Fig. 5A; Supporting Information, Appendix S7). The K = 2 model was suitable for the input SNP data as in the STRUCTURE analysis, where the two lineages present a west-east configuration with some individuals presenting admixture of both lineages in their genetic composition.

The Neighbour-net tree obtained with the SNPs recovered the same two lineages as the Structure results. Individuals 21 and 28 were slightly separated from the eastern lineage (Fig. 6).

Figure 6.

Network analysis inferred from the program SPLITSTREE with the SNPs data. Terminals not grouped in any population represent individuals with mixed lineages.

Ecological niche modelling and niche comparison

The niche models yielded variable performance metrics and those with the highest predictive accuracy were maxent, brt, rf and glm (AUC values > 0.9; TSS values > 0.7) and were used to generate an ensemble model for the current potential distribution (Fig. 7A). The transfer model to Holocene and Pleistocene climate scenarios suggests that suitable habitat for C. ravus was more expanded and continuous at the Younger Dryas Stadial (12.9–11.7 kyr), Bølling-Allerød (14.7–12.9 kyr), Heinrich Stadial (17.0–14.7 kyr) and LGM (21 kyr) than those at the LIG (140–120 kyr), the Holocene: Meghalayan (4.2–0.3 kyr), Northgrippian (8.3–4.2 kyr), Greenlandian (11.7–8.3 kyr) and Present conditions (Fig. 7B). The suitable habitat was fragmented during the LIG and Meghalayan, when climates appeared to be warmer than during the cold periods of the Pleistocene. Eastern and western lineages showed low values of Schoener’s D index (D = 0.175; Supporting Information, Appendix S8) indicating that niche conservatism likely due to high climatic similarity across the region (niche randomization test with P > 0.05).

Figure 7.

A, potential current distribution for C. ravus. B, paleoclimatic reconstructions of for C. ravus (see main text for full details).

DISCUSSION

Diversity and genetic differentiation within C. ravus lineages

The genetic and phylogeographic analyses using mtDNA and SNPs for C. ravus revealed the presence of two highly structured lineages in an east-west configuration on the eastern sector of the TMVB. This pattern is shared with other reptiles (Bryson et al., 2011b; Bryson & Riddle, 2012), plants (Mastretta-Yanes et al., 2018) and birds (Rodríguez-Gómez et al., 2015). The mtDNA showed strong genetic structure (FST = 0.62) between lineages, and high genetic diversity (h = 0.90) and low nucleotide diversity (π = 0.004) for the western lineage compared to the eastern lineage, which had low genetic and nucleotide diversity (h = 0.37; π = 0.003). High genetic diversity and low nucleotide diversity could suggest rapid population growth from an ancestral population with a small effective size (Avise, 2000), as supported for the western lineage in the star-like network topology suggesting recent colonization of the western TMVB. Conversely, the eastern lineage had low genetic and nucleotide diversity, suggesting that populations may have experienced a prolonged or severe demographic bottleneck in recent times (Avise, 2000). This contrasts our BSP results and Fu’s Fs and Tajima’s D; however, the demographic growth time for the eastern population is very recent (3000 years BP), which could indicate that the population (group) has been recovering recently. The latter could be explained by the high volcanic activity that the region has undergone from the Pleistocene to the Present, along with the geographical distribution of C. ravus in a highly active volcanic region, coupled with climate changes that have affected the mountain forests on the TMVB (Demastes et al., 2002).

Although habitat transformation by volcanic activity is high in the region of Valley of Mexico and can limit gene flow and promote population decline (Kim et al., 1998; Sunny et al., 2018), the populations of the western lineage seem to have a high to moderate degree of genetic diversity. This is similar to the findings in dusky rattlesnakes (C. triseriatus) in the State of Mexico (Sunny et al., 2018), where there is high genetic diversity regardless of the human impact in the area. We observed that C. ravus tends to stay within disturbed areas despite human presence and potential risks (A. Cisneros-Bernal pers. obs.). This is particularly true if these areas have patches of mountain pasture vegetation or are surrounded by pine-oak forests, known habitats where this species probably waits for prey that feeds on crops, such as small rodents (Sánchez-Herrera, 1980).

The results of the AMOVA of the mtDNA and SNPs differ regarding where the greatest variation is found. The mtDNA suggests that it is between lineages (populations) (61.63%), while the SNPs suggest that it is within lineages (72.46%). Despite the above, the two markers recover the structure of two lineages with an east and west configuration, which indicates that the evolutionary histories of the two lineages of C. ravus are independent.

Thus, both genetic markers support two C. ravus lineages: an eastern lineage and a western lineage. The formation of these lineages likely is the result of the past complex topographic dynamics of the TMVB. In particular, the formation of stratovolcanoes—a complex fault system in the Basin of Mexico (eastern sector; revised in Arce et al., 2019)—in the middle and last formation periods of the TMVB (Ferrari, 2012; Arce et al., 2019) could have created conditions that promoted the isolation of plant and animal lineages (Mastretta-Yanes et al., 2015), while the glacial periods of the Pleistocene promoted the descent of montane forests, thus connecting them to previously isolated lineage populations, enabling gene sharing in relatively recent times (Mastretta-Yanes et al., 2018). In both Structure and Construct analysis (21 and 28 samples, respectively) and the haplotype network (H21), mixed individuals may be the result of incomplete divergence due to the retention of ancestral polymorphism or a recent divergence between the two lineages, or incipient gene flow (Tang et al., 2012). SPLITSTREE network analyses also recovered the same pattern.

Phylogenetic relationships and divergence times

Bayesian Inference (BI) analysis supports C. ravus as a monophyletic group (Fig. 3). Crotalus brunneus from the Oaxaca population is a sister group to C. ravus, which together are closely related to C. exiguus from the Guerrero population. Crotalus brunneus and C. exiguus had been considered subspecies of C. ravus until recently; Blair et al. (2018) used Ultraconserved Elements (UCEs) to propose elevating them to species status. The phylogenetic relationships of C. ravus have always been controversial, since it is phylogenetically within the genus Crotalus (Bryson et al., 2014; Alencar et al., 2016; Blair & Sánchez-Ramírez, 2016; Blair et al., 2018; Zaher et al., 2019), even though it shares morphological characteristics (nine large symmetrical scales on the head) with the genus Sistrurus, which is considered plesiomorphic (McCranie, 1988; Murphy, 2002).

The inclusion of more species from the C. triseriatus group (Bryson et al., 2011c) in Bayesian analysis allowed us test the existence of cryptic diversity within C. ravus. In other species, such as C. triseriatus, unexpected diversity was found, resulting in a complex of species (Bryson et al., 2011c, 2014). For C. ravus, neither of the two genetic lineages (eastern, western) were associated with another species, so the two lineages are phylogenetically related to each other; however, non-reciprocal monophyly was supported (Fig. 3), perhaps due to the fact that the divergence between the two lineages is recent and there has not been enough time for a complete divergence (Tang et al., 2012).

The estimated divergence times place the origin of the two C. ravus lineages in the Pleistocene (~1.47 Myr). This period was distinguished by glacial-interglacial cycles that allowed the upward and downward shifts of pine-oak forests, mountain grasslands and cloud forests distributed in the TMVB (Vázquez-Selem & Heine, 2004; Ramírez-Barahona & Eguiarte, 2013; Mastretta-Yanes et al., 2015). High volcanic activity was characteristic of the last period of formation of this mountain range in central-eastern Mexico (i.e. 1.5 Myr for TMVB; Ferrari, 2012; Mastretta-Yanes et al., 2015). Based on our estimated divergence time for C. ravus and C. brunneus (~2.46 Myr; Fig. 4), we suggest that the eastern and western genetic lineages of C. ravus were already being influenced by the confluence of the Pleistocene climate changes and the complex geological dynamics of the TMVB, as suggested for other Crotalus snakes (Blair et al., 2018; Bryson et al., 2011). Furthermore, other animal taxa with an eastern geographical distribution on the TMVB, such as salamanders (Parra-Olea et al., 2012) and the volcano rabbit (Osuna et al., 2020), had similar diversification times around 1.5 and 1.4 Myr, respectively.

Demographic and climate dynamics in the Pleistocene

The shifting distribution of C. ravus likely responded to recurrent past climate oscillations and it is probable that the observed anthropogenic climate change had impacts on the snake survival in this region (Böhm et al., 2013). According to our models, C. ravus experienced reductions in its potential distribution on the TMVB in the LIG, while in the LGM and Early Holocene there was a population expansion in its potential distribution (Fig. 7B). It is possible that the less severe Pleistocene climate of Mesoamerica (i.e. different to that of North America) and the long-term stability of suitable habitat conditions allowed C. ravus to undergo a LGM expansion and remain stable until recently (3000 years BP), when it appears to have undergone an expansion. Other taxa exhibited a similar pattern in the mountains of Mexico (Ornelas et al., 2016; Mastretta-Yanes et al., 2018).

We found that suitable habitat for C. ravus was fragmented and reduced during the LIG, unlike during the LGM (Fig. 7B). This suggests that C. ravus likely experienced population and geographic expansion when climate conditions were cooler than present. The dry refugia model generates areas with cooler conditions during the LGM that shifted species migrations (expansion) to the lowlands but the prevalence of arid conditions led to species displacements and contractions into glacial refugia; subsequently, population expansions and long-distance recolonizations took place with the increasingly warm, humid conditions in the Holocene (Ramírez-Barahona & Eguiarte, 2013). This is consistent with our results for the eastern lineage where genetic diversity was low (h = 0.37) and the neutrality tests suggest demographic expansion, while the western lineage seems stable, as observed for other taxa (Ornelas et al., 2016). However, the high niche similarity in these lineages likely is due to the climatic stability of the cloud and pine-oak forests during the Pleistocene (Mastretta-Yanes et al., 2015; Ornelas et al., 2016). If this similarity in the niche preference of the two lineages continues, it could have negative implications due to current climatic changes (Esparza-Estrada et al., 2022).

The results of the Bayesian skyline plots, Fu’s Fs and Tajima’s D tests indicate a pattern of recent demographic expansion in the western lineage of C. ravus species (Fig. 2). The proximity of the western lineage’s populations to the urban area of Mexico City has limited their expansion; changes in land use in areas near large cities has created a mosaic of rural areas surrounded by growing urban areas (Wigle, 2010). The eastern lineage has a larger distribution in the states of Puebla, Tlaxcala and Veracruz than the western lineage does. Conserved areas of pine-oak forests and grasslands are more abundant there (INEGI, 2017) and in surrounding rural areas where, according to our observations, C. ravus can tolerate these habitats (mainly agricultural and livestock areas). Thus, accelerated urbanization is another factor that could also be affecting the genetic structure and demographic expansion of the C. ravus lineages.

If urbanization continues expanding to peripheral areas (Liu et al., 2014, 2016), this may substantially increase habitat loss and fragmentation (Scolozzi & Geneletti, 2012; Güneralp & Seto, 2013), likely leaving the extant lineages of C. ravus completely isolated. Further study of the historical demography is necessary to establish the potential impacts from future urbanization in the Basin of Mexico in both C. ravus lineages because this region is overpopulated and it is projected to expand in the future (INEGI, 2017).

CONCLUSION

Multiple phylogeographic studies documented cryptic diversity in the TMVB (Demastes et al., 2002; Bryson et al., 2011c; Mastretta-Yanes et al., 2015, 2018; Blair et al., 2018; Osuna et al., 2020). Our results with C. ravus are congruent with volcanism and other geological events, plus Quaternary climate change, as stronger triggers of species diversification in the TMVB. The estimated divergence in C. ravus lineages is temporally and geographically congruent with the early formation of the TMVB, suggesting that the development of this mountain range played a critical role in the diversification of this lineage. More genetic sampling in each locality will help us better understand which processes gave rise to the pattern observed in our data.

SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this article at the publisher’s web-site.

Appendix S1. Protocol of DNA extraction, amplification and purification of each gene in this study.

Appendix S2. Specification for the Bayesian skyline plot analysis.

Appendix S3. Analysis specifications of the calibrated phylogenetic tree.

Appendix S4. RADSeq library preparation and bioinformatic analysis specifications.

Appendix S5. Neighbour-net tree obtained from mitochondrial DNA data. The eastern haplogroup is in blue and the western haplogroup is in grey. Individuals represented by red dots were grouped into populations other than those where they were collected.

Appendix S6. Demographic indices and analysis of molecular variation (AMOVA) obtained with the SNPs data with 50% of the data missing. Indices with an asterisk (*) indicate significance values of less than 0.05.

Appendix S7. Bar plot of the lineage’s composition of each individual of C. ravus obtained in construct with a K = 2 model. The y-axis refers to the percentage of mixing and the x-axis refers to the individuals.

Appendix S8. Randomization test for niche comparisons between clades. The red line indicates the observed Schoener’s D value (0.175).

ACKNOWLEDGEMENTS

We thank the two external reviewers for their comments and observations on the manuscript. A.Y.C.-B. would like to thank the members of the Fujita Laboratory at University of Texas at Arlington, Texas, for their help in the sequencing work, especially T.J. Firneno, and also F. Ramírez-Corona of the Taller de Sistemática y Biogeografía, Facultad de Ciencias- Universidad Nacional Autónoma de México, and M.E. Muñiz Díaz de León of the Taller de Biología de Plantas, Universidad Nacional Autónoma de México. A.Y.C.-B. also thanks all members of the Herpetology Laboratory at Facultad de Ciencias-Unversidad Nacional Autónoma de México and the Evolutionary Ecology of Amphibians at Facultad de Estudios Superiores, Iztacala and Reptiles Laboratory at Facultad de Estudios Superiores Iztacala for providing help and samples. Finally, A.Y.C.-B. would like to thank the Posgrado en Ciencias Biológicas at UNAM, as well as CONACYT (project no. 669772) and PAPIIT (IN216218) for funding. J.A.V. is grateful for financial support from the Programa Universitario de Investigación en Cambio Climático (PINCC-UNAM) through a project grant. Collection permits were issued by SEMARNAT (FAUT-0015) to O.F.-V. We have no conflicts of interest to declare.

DATA AVAILABILITY

The mitochondrial genetic data used in this work are available in GenBank. The ID numbers of each sequence in GenBank are referred to in Table 1. ddRADSeq data (SNPs) are available at the following link on the open access platform Open Science Framework (OSF): http://osf.io/3upg/.