Historical biogeography of the genus Rhadinaea (Squamata: Dipsadinae)

Abstract Multiple geological and climatic events have created geographical or ecological barriers associated with speciation events, playing a role in biological diversification in North and Central America. Here, we evaluate the influence of the Neogene and Quaternary geological events, as well as the climatic changes in the diversification of the colubrid snake genus Rhadinaea using molecular dating and ancestral area reconstruction. A multilocus sequence dataset was generated for 37 individuals of Rhadinaea from most of the biogeographical provinces where the genus is distributed, representing 19 of the 21 currently recognized species, and two undescribed species. Our analyses show that the majority of the Rhadinaea species nest in two main clades, herein identified as “Eastern” and “Southern”. These clades probably diverged from each other in the early Miocene, and their divergence was followed by 11 divergences during the middle to late Miocene, three divergences during the Pliocene, and six divergences in the Pleistocene. The ancestral distribution of Rhadinaea was reconstructed across the Sierra Madre del Sur. Our phylogenetic analyses do not support the monophyly of Rhadinaea. The Miocene and Pliocene geomorphology, perhaps in conjunction with climate change, appears to have triggered the diversification of the genus, while the climatic changes during the Miocene probably induced the diversification of Rhadinaea in the Sierra Madre del Sur. Our analysis suggests that the uplifting of the Trans‐Mexican Volcanic Belt and Chiapan–Guatemalan highlands in this same period resulted in northward and southward colonization events. This was followed by more recent, independent colonization events in the Pliocene and Pleistocene involving the Balsas Basin, Chihuahuan Desert, Pacific Coast, Sierra Madre Occidental, Sierra Madre Oriental, Sierra Madre del Sur, Trans‐Mexican Volcanic Belt, and Veracruz provinces, probably driven by the climatic fluctuations of the time.


| INTRODUC TI ON
Inferring the evolutionary history of groups in a particular region is the first step in elucidating the processes by which the region's fauna originated (Colston et al., 2013). Multiple biogeographical and phylogeographical studies of groups with broad distributions through North and Central America are illustrative in this context (e.g., Burbrink et al., 2008;Fritz et al., 2012;Hofmann & Townsend, 2017;Phillips et al., 2015). Within this area, the Mexican Transition Zone (MTZ) is a complex, distinguishable area where Neotropical and Nearctic biotas overlap, spanning the region from the southwestern deserts of the United States and northern Mexico to the dry and humid forests of the Nicaraguan lowlands (Morrone, 2010).
Across the North and Central American regions, multiple geological and climatic events resulting from a complex orogeny and paleoclimatic conditions acted as geographical or ecological barriers associated with the speciation and diversification of many taxa (e.g., Daza et al., 2010;Ferrusquía-Villafranca & González-Guzmán, 2005;Vanzolini, 1970). Some of the events that are considered of major importance or that have received the most attention are as follows: (a) The Mississippi River Basin (MR) (Burbrink et al., 2000), which was involved in the divergence of many marine and terrestrial taxa during the Pleistocene (Soltis et al., 2006) (Ferrari et al., 2012;Gómez-Tuena et al., 2007) in four major orogenic events (Ferrari et al., 2012) that undoubtedly affected both the timing and tempo of the biota diversification (Bryson et al., 2012a(Bryson et al., , 2012b; (d) the faulting and marine introgressions across the Isthmus of Tehuantepec (IT) in southeastern Mexico around 3 Ma (Mulcahy et al., 2006), a region which is a narrow lowland area that has been identified as a biogeographical barrier for many upland taxa (Castoe et al., 2009); (e) the Nicaraguan depression (ND), an area that presented different states of terrestrial conformation during the Neogene (2.5-23 Ma) (Funk et al., 2009) and probably presented a lowland biogeographical barrier to some taxa (Daza et al., 2010); (f) the Panama Isthmus in southern Central America, another narrow area that was completely conformed during the Pliocene (3.5 Ma), which has separated numerous taxa between Central and South America (Mendoza et al., 2019); and (g) the climatic fluctuations during the Pleistocene (0.01-2.5 Ma) (Vanzolini, 1970) that conditioned the diversification of a variety of taxa across the American continent through the repeated expansion and contraction of coniferous forests, leading isolated populations of forest-adapted taxa to speciation (Haffer, 1969(Haffer, , 1997.
These events, in addition to other physiographic conditions, such as river drainages within the major sierras, basins, and faults, are considered to act as biogeographical barriers (Bryson, Murphy, et al., 2011;Bryson et al., 2007;Daza et al., 2010;León-Paniagua et al., 2007), yet the effectiveness of these barriers in isolating lineages throughout the past several million years remains to be clarified (Bryson et al., 2012a(Bryson et al., , 2012b. The colubrid snakes of the genus Rhadinaea are slender, diurnal, medium-to small-sized snakes that are characterized by longitudinal dark stripes along the dorsal scales (Figure 1), a small subpreocular scale inserted between the corners of two supralabial scales at the anteroventral edge of the orbit, the same number of longitudinal dorsal scale rows throughout the body, and maxillary teeth without grooves posterior to the diastema (Myers, 1974(Myers, , 2011Palacios-Aguilar & García-Vázquez, 2020). Currently, 21 species of Rhadinaea are recognized and arranged into five morphological groups Myers, 2011;Palacios-Aguilar & García-Vázquez, 2020). The genus is present in North and Central America from southeastern United States to Panama, with discontinuities across the Chihuahuan Desert, southern Guatemala, El Salvador, Honduras, and central Nicaragua (Table 1) (Myers, 1974(Myers, , 2011. The most speciose and widely distributed of these groups is the decorata group, represented by 12 species mainly distributed over the Mexican Sierras (García-Vázquez, 2012;García-Vázquez et al., 2009;Luría-Manzano et al., 2014;Pérez-Higareda et al., 2002;Sánchez-García et al., 2019;Torres-Carvajal et al., 2019); the taeniata group is endemic to Mexico and is composed of three species distributed in central Mexico (Canseco-Márquez & Gutiérrez-Mayén, 2010;García-Sotelo et al., 2018;Myers, 1974); the flavilata group is comprised of two species with allopatric distributions in North America (Auth et al., 1999;Lares et al., 2013;Walley, 1999); and the calligaster and vermiculaticeps groups are restricted to Central America (Myers, 1974). In addition to these 21 species, the existence of two undescribed species has been suggested based on their morphology and previous analyses (pers. obs.).   Bailey, 1939;Cadle, 1984;Myers, 1974;Zaher et al., 2019). Additionally, Rhadinophanes is a monotypic genus that is considered the sister group of Rhadinaea and Coniophanes (Cadle, 1984;Myers & Campbell, 1981), a relationship that has not been tested using molecular data. To explore their monophyly and evolutionary history, Palacios-Aguilar and García-Vázquez (2020) suggest that it is necessary to include more representatives of the genera using more molecular markers.
To date, the five groups composing Rhadinaea are sorted out based on their color patterns, which are the most informative character to distinguish them, even among species within Rhadinaea (Myers, 1974(Myers, , 2011. However, the reciprocal monophyly of the groups is still to be assessed . On the other hand, according to Myers (1974), the origin of Rhadinaea is related to the dispersal of an ancestor related to Rhadinella (former godmani group), in which geographical isolation and the subsequent evolution of terminal populations occurred, probably after unfavorable climatic changes or flooding that created barriers of lowland regions, such as the IT and the ND (Myers, 1974); yet, other geographical barriers corresponding to the distribution of the genus, such as the formation of the TVB or climatic events, are not discussed.
In this work, we describe the phylogenetic relationships of Rhadinaea to evaluate the role of major orogenic events and Pleistocene climatic fluctuations on lineage diversification. Samples of the five recognized groups and two undescribed species of Rhadinaea were included (R. cf. marcellae and R. cf. taeniata). Two mitochondrial and two nuclear loci were sequenced. We inferred the phylogenetic relationships and a time-calibrated tree from these data. Finally, ancestral ranges were reconstructed at each divergence event. The resulting patterns of diversification are discussed in the context of mountain formations and climatic change.

| Taxon sampling and laboratory methods
In this study, we cover 37 samples of Rhadinaea, including most of the currently recognized species of the genus ( Figure 2; Appendix S1),

F I G U R E 2
Sampling localities for the genetic samples used in this study (see Appendix S1). Black lines indicate political boundaries; red lines indicate biogeographic regions (Escalante et al., 2013;Löwenberg-Neto, 2014;Morrone et al., 2017). Alleghanian Region ( with the exception of R. sargenti and R. vermiculaticeps, two species that are rarely found within biological collections. Samples from two undescribed species from the TVB were also included. To increase geographic representativeness for species with a wide distribution range, a sample per biotic region was included, following the biogeographical regionalization of North America (Escalante et al., 2013), Mexico , and Neotropical America (Löwenberg-Neto, 2014).
Partial sequences of the mitochondrial gene coding for the cytochrome b (cytb), NADH dehydrogenase subunit 4 (ND4); complete sequences of the noncoding tRNA-His and tRNA-Ser; partial sequences of the noncoding tRNA-Leu; partial sequences of the nuclear genes coding for oocyte maturation factor (cmos) and dynein axonemal heavy chain 3 (DNAH3) were obtained for all the 37 individuals of Rhadinaea and five outgroup individuals. Loci were selected as they previously showed to be informative at different levels of divergence between snakes (Bryson, García-Vázquez, et al., 2011;Bryson, Murphy, et al., 2011;Lawson et al., 2005;Myers et al., 2017).
Primer sequences for cmos were given by Lawson et al. (2005) and Saint et al. (1998); for cytb by Burbrink et al. (2000), de Queiroz et al. (2002), and Slowinski and Lawson (2002); for DNAH3 by Townsend et al. (2008); and for ND4 by Arévalo et al. (1994), Forstner et al. (1995. Also, additional internal primers were designed in this study. See Appendix S2 for primer sequences, technical details on DNA sequencing, and sequence edition.

| Phylogenetic inference
With regard to testing incongruences among mitochondrial and nuclear loci, we performed Bayesian inference (BI) and maximumlikelihood (ML) separate analyses for each nuclear and mitochondrial locus. The best-fitting substitution models and partition schemes were selected jointly using the Bayesian information criterion in the software PartitionFinder 2 (Lanfear et al., 2016). BI analysis was conducted using MrBayes 3.2.7a (Ronquist et al., 2012) with four Monte Carlo Markov chains (MCMC), sampling every 5,000 generations for 100 million generations. Output parameters were visualized using Tracer 1.7.1 (Rambaut et al., 2018) to ascertain stationarity and convergence. The first 25% of generations were discarded as burn-in to obtain a majority rule consensus tree using the command sumt. ML analysis was conducted using raxmlGUI (Silvestro & Michalak, 2012) under the GTRGAMMA model (Stamatakis, 2006) with 1,000 nonparametric bootstrap replicates to assess nodal support.
Additionally, all nuclear and mitochondrial datasets were combined into one dataset after independent analyses established the congruence between topologies and levels of support between the IB and the ML analyses, using the same partitions and models of evolution suggested by PartitionFinder 2 (Lanfear et al., 2016), and the concatenated matrix was analyzed under BI and ML approaches following the above specifications and settings in each program. Nodes were considered strongly supported if their Bayesian posterior probability (pp) was ≥0.95 and their bootstrap (bs) value was ≥70% (Huelsenbeck & Rannala, 2004).

| Divergence times
Divergence dates and phylogeny were estimated simultaneously using a relaxed Bayesian molecular clock framework implemented in

| Ancestral area reconstruction
Ancestral ranges at each divergence event were reconstructed using the Bayesian binary Monte Carlo analysis (BBM) and dispersalextinction-cladogenesis (DEC) implemented in RASP 2.0 (Yu et al., 2011). This program can determine the probability of an ancestral range at a node by averaging a posterior set of trees and thereby accounting phylogenetic uncertainty (Bryson et al., 2013). A total of 16,000 post-burn-in trees and the condensed maximum clade credibility tree obtained in the BEAST analysis were loaded from the divergence time analysis into RASP, removing the outgroups and the Rhadinella clade using the "remove selected groups" tool that is implemented in the same program. According to the distribution of each species (Table 1) for the nodes in the phylogeny were estimated, whereas the analysis was conducted for 10 million generations by sampling each 1,000 using ten chains, and the first 25% of the generations were discarded as burn-in for BBM. For DEC analysis, we set the dispersal rate all equal among the defined areas, performing 100 replicas.

| Phylogenetic inference
Individual genealogies obtained using IB and ML did not share the same topology, as these genealogies show numerous nonsupported nodes and polytomies; however, we did not observe discordances between ML and IB reconstructions for each gene in shared supported nodes, where all of these were congruent (Appendix S4).
Due to this situation, we assumed that the phylogenetic information that each sequence dataset provided could be combined to obtain a more accurate phylogeny, whereas congruent data can be combined to yield phylogenies that do not represent organismal history accurately (Cunningham, 1997;Hipp et al., 2004).

| Divergence times
Our multilocus analysis produced a reconstruction for Rhadinaea with moderate resolution and node support (75% of nodes with pp > 0.95). In the calibrated tree, the same clades as phylogenetic analyses were recovered, as well as the relationships between the clades. The dated phylogeny suggests that the diversification of  pulveriventris (14.2 Ma), followed by five splits among these species.
Our estimates placed the remaining divergences within the Eastern and southern clades during the Pliocene and Pleistocene (Figure 4).

| Ancestral area reconstruction
The BBM and DEC analyses using the dated multilocus phylogeny obtained in BEAST were conducted excluding R. calligaster (see Discussion). For the DEC analysis although some ambiguity and possible alternative resolutions exist, the results were consistent with BBM analysis (Appendix S5). Because BMM shows highest resolution at nodes, we considered it most likely for the hypotheses here.
Also, the estimation of ancestral area marginal probabilities taking into account phylogenetic uncertainty (Bayesian-like) has been suggested to reduce uncertainty in the biogeographical reconstruction (Nylander et al., 2008). However, we comment about the major difference of DEC respect to BBM.
The RASP analyses showed that the diversification of Rhadinaea

| Phylogenetic relationships within Rhadinaea
Our molecular-based phylogeny generally does not support the traditional taxonomy of Rhadinaea according to morphology (Myers, 1974(Myers, , 2011. Our analyses place Rhadinaea calligaster as the sister group of Rhadinella (Figures 3 and 4). This relationship is not bers (Myers, 1974), which makes it difficult to accurately diagnose the members of this genus given the lack of exclusive diagnostic characters and the existence of variable morphological traits, such as the form of the neck collar, pigmentation, and longitudinal patterns of the dorsal scales (Myers, 1974;Sánchez-García et al., 2019).
This situation has resulted in the further modification of the scientific understanding of the composition of the genus (Myers, 2011;Palacios-Aguilar & García-Vázquez, 2020). About this, Myers (1974) described the relationships of this monospecific group as obscure and mentioned that morphological features like a bilobated hemipenis, the absence of a subpreocular scale, and a pale bar from the eye to the corner of the mouth indicate an ancestry from the godmani group (now Rhadinella). This relationship is plausible, given our results. Therefore, this leads us to infer that Rhadinaea is not monophyletic, in contrast with the results obtained by Palacios-Aguilar and García-Vázquez (2020). Whereas Rhadinaea and Rhadinella are recovered as reciprocally monophyletic in their study, here we only retrieved the monophyly of Rhadinella, which is yet to be verified including a larger number of taxa.
The classification of the species of Rhadinaea into groups defined by Myers (1974Myers ( , 2011 shows inconsistencies that involve the decorata, flavilata, and taeniata groups. Our phylogenetic hypothesis presented the decorata and taeniata groups as polyphyletic and the flavilata group as paraphyletic (Figures 3 and 4), indicating that the relationships among these species need a further revision due to the low nodal support in the clades containing these species and that they probably do not represent natural groups. Following these observations, we also suggest that the phylogenetic position of R. calligaster needs to be revisited in order to confirm the close relationship with Rhadinella observed in our analyses, taking into account the evidence that suggests that this group may represent another Rhadinella species.
We recovered the majority of the species with more than one sample as monophyletic with strong support values, except for Rhadinaea marcellae (Figures 3 and 4). The species that did not show a monophyletic pattern are R. montana and R. quinquelineata, and appeared to be paraphyletic groups (Figures 3 and 4). This pattern might be the result of an unclear delimitation of these species, as previously suggested by Dixon et al. (2011), with respect to the previously proposed morphological series consisting of R. gaigeae, R. montana, and R. quinquelineata (Dixon et al., 2011;Myers, 1974), species that show a very similar morphology and are codistributed over the SMOr (Canseco-Márquez et al., 2004;Medina-Romero et al., 2016;Myers, 1974). Because of this situation, we decided to include each haplotype of these species in divergence time estimation and ancestral area reconstruction analyses to explore their relationships and diversification. Furthermore, we suggest that additional studies of these taxa including better sampling, additional loci, and explicit testing of alternative species hypotheses using coalescent methods for species delimitation (e.g., Fujita et al., 2012;Yang & Rannala, 2010) are needed to properly address their status and kinship.
Concerning the closest relatives of the genus, it has been shown that Rhadinaea is closely related to Amastridium, Coniophanes, Rhadinella, and Tantalophis (Daza et al., 2010;Palacios-Aguilar & García-Vázquez, 2020;Pyron et al., 2013). Pyron et al. (2013) found a strong relationship between Rhadinaea and Coniophanes and between both genera and Rhadinella, according to Palacios-Aguilar and García-Vázquez (2020). In our study, these last two relationships were recovered in the analyses, even though the relationship of with a low support value in the BI analysis. Due to the complex morphology of this species and uncertainty in its phylogenetic position (Myers & Campbell, 1981), we cannot assess whether this species is more closely related to Rhadinella or Rhadinaea calligaster given the low support value in our analyses (Figure 3), in which we note that Rhadinophanes monticola is not closely related to the main clades of Rhadinaea we identified herein (Figures 3 and 4). This outcome could be clarified, considering a broader perspective, including other related colubrid species, as in other studies (e.g., Daza et al., 2010;Pyron et al., 2013).

| Historical biogeography
Based on our results of BBM and DEC, it appears that colubrid snakes of the genus Rhadinaea have had a relatively long history in North and Central America. The SMS is an extensive mountain system that has been present in the Mexican territory since its formation due to the Laramide geologic activity (DeMets & Stein, 1990)  has experienced changes during similar times (Ornelas et al., 2010(Ornelas et al., , 2013. This mosaic-like landscape has been associated with centers of diversification along elevational gradients and is believed to be closely related to some divergences between some vertebrate taxa  Myers, 1974), and the heterogeneity of the SMS (Bryson et al., 2017;Luna-Vega et al., 1999;Santiago-Alvarado et al., 2016).
In this same period, two dispersal events from the SMS to other areas took place. First, we note inside the Eastern clade a northward dispersal of the MRCA of R. flavilata, R. gaigeae, R. montana, and R. quinquelineata (Figure 4), which became widespread during the middle Miocene. This event corresponds temporally with the first episode of volcanic and orogenic activity that originated the TVB in the central portion of Mexico (Bryson et al., 2012a(Bryson et al., , 2012bFerrari et al., 2012;Ferrusquía-Villafranca & González-Guzmán, 2005;Gómez-Tuena et al., 2007). This mountain system is known for its influence on the diversification of various montane taxa, creating new montane habitats (Bryson et al., 2012b) and probably allowing further colonization of more mesic-adapted lineages (e.g., García-Vázquez, Nieto-Montes de Oca, et al., 2018;Milstead, 1960) such as R. flavilata, in coordination with the low temperatures posterior to the MCO (Zachos, 2001). Secondly, we observe a dispersal event during the Miocene inside the southern clade toward the south by the MRCA of R. decorata and R. pulveriventris ( Figure 5). This divergence is generally consistent with the formation of the Chiapan-Guatemalan highlands of northern Central America, which formed during two different time intervals (Campbell, 1999). The uplift of the extensive northern Central American plateau occurred during the late Miocene to early Pliocene, from approximately 10-3.8 Ma (Rogers et al., 2002). Additionally, the formation in the late Pliocene of a younger chain of volcanoes along the western portion of the Central American plateau (Williams, 1960) had a significant impact on the local biota, both through extinction and the resulting climatic change, creating cloud forest conditions on the windward (south) slopes and rain shadow conditions in the interior valleys (Campbell, 1999). Moreover, an abrupt turnover from xeric, subhumid vegetation to humid forests occurred in south-central Chiapas and extended along the coast to south-central Guatemala (Campbell & Vannini, 1988 (Castoe et al., 2009), where significant biogeographic barriers in Central America, such as the Motagua-Polotchic fault (Marshall, 2007) and the Nicaraguan Depression (Marshall, 2007;Rogers et al., 2002), played an essential role in the diversification of several taxa (Campbell, 1999;Devitt, 2006;Parra-Olea et al., 2004;Perdices et al., 2005;Savage, 1982), probably this explain the vicariance event found by DEC analysis.
Along with these critical orogenic events, the Miocene climate change appears to have played an important role, sparking evolutionary radiations in some successful modern lineages, including colubrid snakes, and segregating the species along latitudinal and altitudinal environmental gradients (Greene, 1997;Van Devender & Spaulding, 1979). These conditions probably have had an influence in other colonization events during this time within the Eastern and southern clades, such as the colonization of the SMOr by the MRCA of R. gaigeae, R. montana, and R. quinquelineata ( Figure 5), as well as a posterior divergence in the SMOr, perhaps produced by changing ecosystems associated with the wetter climate during this period through the Northern Mexican highlands (Bryson et al., 2013(Bryson et al., , 2017Retallack, 2001) as seen in other reptile groups (e.g., Gerrhonotus [ García-Vázquez, Nieto-Montes de Oca, et al., 2018] Figure 5).
This last hypothetical species is only known from the westernmost portion of the TVB in Sierra San Juan, a too complex area with a high number of endemisms (Escalante & Llorente, 1985;Miranda & Luna-Vega, 2006).
During the Pliocene, the last episodes of formation of the TVB  (Paillard, 2017;Vanzolini, 1970). These cool intervals are believed to have caused the expansion of some pine, pine-oak, and humid forests to lower elevations due to temperature fluctuations, provoking an extension of the Mexican montane flora to lower elevations of at least 1,000 m (Jaramillo- Correa et al., 2009;McDonald, 1993;Ornelas et al., 2013;Sosdian & Rosenthal, 2009). Therefore, these orogenic and climatic factors in more recent conjunction are considered to be the processes that have had a more significant impact in a high number of montane taxa (as discussed in Bryson et al., 2012b), and throughout these Pliocene and Pleistocene periods, most extant species of Rhadinaea among the Eastern and southern clades originated and colonized other regions.
Within the MTZ, the contact between the SMOr, SMS, TVB, and VER regions is characterized by a very complex biotic interchange, presenting a significant amount of shared floristic and faunistic elements (Espinosa et al., 2004;Marshall & Liebherr, 2000), some of which showed a recent dispersal to adjacent provinces during the Pliocene and Pleistocene (Cavender-Bares et al., 2011).

| CON CLUS IONS
Biogeographical studies seek to explain the distribution of spe- Smith and an anonymous reviewer for commenting of manuscript.
All contributors observed state, federal, and international laws concerning the collection and transport of live or preserved specimens.
All specimens were collected under a collecting permit provided by the Secretaría de Medio Ambiente y Recursos Naturales (permit number FAUT-0243) of the Mexican Government to U.O. García-Vázquez. This work is part of the master's degree research project of U. A. García-Sotelo (Evolutionary Biology).

CO N FLI C T O F I NTE R E S T
The authors Uriel Alonso García Sotelo, Uri Omar García Vázquez, and David Nahum Espinosa Organista declare that we have no significant competing financial, professional, or personal interests that might have influenced the performance or presentation of the work described in this manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
All sequence data generated in this study will be available in the public genetic sequence database GenBank (https://www.ncbi.nlm.nih. gov/genba nk/).