The role of Central American barriers in shaping the evolutionary history of the northernmost glassfrog, Hyalinobatrachium fleischmanni (Anura: Centrolenidae)

Tissue sampling

Genetic material was obtained across the entire distribution area of the species, from both museum collections and fieldwork (Fig. 1). Fieldwork was performed during the rainy season, in which at least three individuals were collected at each locality. Specimen collection permits were provided by the Ministerio de Medio Ambiente, Colombia (Resolution 120 of 24 August 2015) and the Secretaría del Medio Ambiente y Recursos Naturales, Mexico (office number 00947/16). Captured specimens were euthanized with a 20% lidocaine hydrochloride (Xylocaine) injection, and all efforts were made to minimize suffering. Liver or muscle tissues were collected in the field and were stored in an RNAlater solution until their use in the laboratory. Specimens were fixed with 10% formalin, stored in 70% ethanol and deposited in biodiversity collections at public research institutions in each country.

Molecular techniques

DNA was obtained from muscle and liver tissues following the phenol-chloroform extraction protocol (Sambrook & Russell, 2006). The quantity and quality of the DNA were verified in 1% agarose gels and by measuring absorbance using a NanoDrop spectrometer (Thermo Fisher Scientific Inc., Wilmington, DE, USA). Amplification of the mitochondrial COI (658 bp), 16S (895 bp), and ND1 (961 bp) genes was performed following the protocols described by Guayasamin et al. (2008). PCR products were visualized with agarose gels and purified according to the EXO-SAP protocol (GE Healthcare, Chicago, IL, USA). DNA sequences were obtained with the BigDye Terminator Cycle Sequencing kit (Applied Biosystems). Sequences from H. tatayoi from Venezuela (Guayasamin et al., 2008; Castroviejo-Fisher, Ayarzagüena & Vila, 2007) were also included. The sequences were assembled, manually edited and aligned with Geneious 9.1.2., and the three genes were concatenated by using the cbind function in the ape package in R (Paradis, Claude & Strimmer, 2004).

We are aware that analyses based on mitochondrial DNA (mtDNA) provide a limited view of species evolution (i.e., matrilineal inheritance), but the available genetic material from nuclear Central America is very scarce since there has been limited or no sampling in the region (due to a lack of funding as well as logistical and even bureaucratic difficulties) compared to that in northern or southern countries. Therefore, we chose to compile a database including gene sequences available from GenBank and BOLD databases with our own sequences to perform our analyses. In this way, our study encompasses the entire geographic distribution of H. fleischmanni, allowing us to evaluate the events shaping the evolutionary history of the species.

Additionally, mtDNA is a robust indicator of population structure and historical biogeographic barriers and can be used to detect cryptic diversity (Gehara et al., 2014; Dudoit et al., 2018; Lait & Hebert, 2018). Mitochondrial DNA markers have been used in a variety of Central American studies to decipher genetic structuring (for a review, see Bagley & Johnson, 2014), allowing us to make direct comparisons with previous studies from this region.

Data analysis

Phylogenetic analyses

Since phylogeographic breaks can be detected in the form of phylogenetic splits between mostly distinct geographical lineages (Bagley & Johnson, 2014), the sequences of all genes were concatenated, and a phylogenetic tree was estimated at the intraspecific level by implementing both likelihood analysis in RAxML (Stamatakis, 2006) and a Bayesian inference approach in MrBayes (Ronquist & Huelsenbeck, 2003). We rooted our phylogeny using the species Hyalinobatrachium carlesvilai, H. mondolfii, H. chirripoi and H. colymbiphyllum as outgroups. A list of the specimens and GenBank accession numbers included in this study is presented in Table S1. The best evolutionary model for each noncoding region (16S) and for the coding genes (COI and ND1) was evaluated using PartitionFinder 2 software (Lanfear et al., 2016). Maximum likelihood analysis was conducted using 10,000 rapid bootstrap analyses, the GTR+Γ evolution model and summarized support for the best tree. For Bayesian inference, we ran two independent analyses for 12 million generations, sampling trees and parameter values every 1,000 generations. Burn-in was set to 25%, and the first 3 million generations were therefore discarded.

Spatial clustering analysis

The Nei’s pairwise F_ST results showed significant differences between adjacent regions for all barriers except for the HE and the Eastern Cordillera (Table 1). The Geneland map of population membership for the concatenated genes revealed nine clusters that partially coincided with our hypothesized barriers, depicting the IT, the MPJ fault system, the Talamanca Range and the WPI as barriers. Additionally, three non-considered barriers were indicated between samples from west of the IT and on the Chortis block, between samples from both sides of the Darien and between the Ecuadorian and Colombian samples (Fig. 2A).

sPCA performed on individual genotypes revealed a significant, geographically correlated genetic structure for all three genes (nper = 999, P = 0.001). Eigenvalues indicated a higher spatial genetic structure for the main axis, related to the global structure. The positive axis of the first sPCA (regional scale) exhibited the greatest variation in genetic distance in relation to the distance network (Fig. 2B). The samples separated by the WPI, the Talamanca Range and the MPJ fault system did not overlap based on the two sPCA first components, while the samples separating the IT and HE exhibited partial overlap (Fig. 2C).

Phylogenetic patterns, times of divergence, and demography

The PartitionFinder output indicated that HKY+I+G, GTR+G, and GTR +G+I were the best models for 16S, COI, and ND1 respectively. The phylogenetic relationships based on the Bayesian and maximum likelihood approaches for the concatenated genes indicated three main, well-supported clades, although their relative positions were not fully resolved (Fig. 3). The first clade, designated the “Northern clade”, was divided into two lineages: a large lineage containing all samples from the SMB and Maya regions (pp = 1) and a smaller one from the western Chortis region (pp = 1). The second clade, the “Central clade”, was comprised of samples from the eastern Chortis and Chorotega regions (pp = 1). The third clade, the “Southern clade”, consisted of samples from the Chocó and SAP regions, including the species H. tatayoi from Venezuela (pp = 1). The Southern clade did not show any structure, displaying a polytomic topology.

The divergence time estimation results showed a pattern of divergence among the three main clades occurring during the Pliocene (∼3.40 Ma, HPD = 2.25–4.56 Ma; Fig. 3). With respect to the Northern clade, the split between the lineage from West Chortis and the remaining samples also occurred during the Pleistocene in the Gelasian age (∼2.19 Ma, HPD = 1.38–3.53), while separation between samples from the Pacific and Gulf+Maya regions occurred at the beginning of the Calabrian age (∼1.76 Ma, HPD = 0.86–2.23). Finally, the split between lineages from the Gulf and Maya regions occurred near the end of the Calabrian age (∼1.51 Ma, HPD = 0.62–1.73). The divergence between the Central and Southern clades occurred at the end of the Pliocene (∼2.64 Ma, HPD = 1.50–3.68), while splits within each clade began at the end of the Calabrian age for the Central clade (∼0.81 Ma, HPD = 0.12–0.71) and during the Gelasian age (∼1.64 Ma, HPD = 0.63–1.89) for the Southern clade (Fig. 3).

The 95% CIs of the effective population size (BSP results) overlapped along the entire time period in the Northern clade (Fig. 4A). However, the ancient and recent effective population sizes of the Southern clade did not overlap (by 95% CI of the Bayesian posterior probability), providing significant support for a change in population size and suggesting that the clade exhibited a constant population size and posterior expansion at approximately 0.1–0.3 Ma (Fig. 4B).

The results regarding Bayesian spatiotemporal diffusion (Fig. 5) highlight the Chorotega and West Chortis region as the most likely ancestral geographic area for H. fleischmanni, suggesting two subsequent dispersal events, in which the MRCA of the Northern clade was distributed in the environs of the Chortis and Maya regions, whereas that of the Southern clade (stem node) was distributed around the Chorotega and Chocó regions. Our results also reflect independent dispersal for samples west of the MPJ fault system and later divergence of the three remaining clusters in the Northern clade around the IT. For the Central clade, an initial range in the Chorotega region south of the Talamanca range was observed, with posterior dispersal towards eastern Chortis. Furthermore, an ancestral range was detected in the Chocó block within the Southern clade with subsequent dispersal towards the south and east, reaching the lowlands east of the Andes range to the south, while two lineages dispersed northwards independently, reaching the southern limit of the Chocó region.

Genetic diversity and structure

Our results provide evidence of four of the hypothesized barriers (WPI, Talamanca Range, MPJ fault system and IT) and identified two additional barriers: one east of the IT and one within the Chortis block. Hence, we identified seven genetically homogeneous regions: North American Pacific, Gulf of Mexico, Maya block, western Chortis, eastern Chortis-South Talamanca Range, North Talamanca Range and Chocó-SAP. A minimum genetic distance (Table 2) was exhibited between regions for North Talamanca Range and Chocó-SAP (K2P = 1.6%), while the maximum was observed between eastern Chortis and North Talamanca Range (K2P = 6.0%). Significant Nei’s F_ST indices between regions were obtained for almost all combinations (F_ST >0.4), except between North Talamanca Range and Chocó-SAP (Table 2).

The Southern clade

The Southern clade encompasses samples from the Chocó and SAP regions, including the species H. tatayoi from Venezuela. Samples ranging from Panama to Venezuela and Ecuador grouped together, with no clear phylogenetic separation among them. However, haplotype networks and spatial clustering analyses for coding genes allowed us to identify a partition on both sides of the Northern Andes (Supplemental Information 1 and Fig. S3). The lack of significant structure for the Southern clade is remarkable, considering that the geographic distances between populations reach 1,600 km. In addition, the CPI and the Northern Andes, which are widely recognized as speciation drivers for both highland and lowland species (Bagley & Johnson, 2014; Mendoza et al., 2015), do not seem to have exerted any effect on the genetic structure of this clade. Notably, the genetic distance observed on both sides of the Northern Andes (COI K2P = 1.9%; Table S1.1) is lower than the distances reported for lowland species with a higher dispersal capacity from similar localities, such as the hummingbird Amazilia amabilis (COI K2P = 2.06%, Mendoza et al., 2016). These results contrast with previous knowledge of the ecology of glassfrog species, which have been found to be characterized by site fidelity (Valencia-Aguilar, Castro-Herrera & Ramírez-Pinilla, 2012), low mobility and restricted gene flow even at local scales (Delia, Bravo-Valencia & McDiarmid, 2017; Robertson, Lips & Heist, 2008). However, most of the previous research on this frog has focused on calling males and reproductive behavior (Delia, Bravo-Valencia & McDiarmid, 2017), while the dispersal capability of females and tadpoles, which can have a significant impact on mtDNA genetic structure, is still unknown. Thus, it is possible that the Chocó and SAP regions present adequate conditions for tadpole dispersal, allowing range expansion. However, this hypothesis needs to be evaluated based on additional demographic studies and a greater sample size per site.

Our results do not support our initial hypothesis that the species was originally from South America and then dispersed through the Isthmus of Panama. Indeed, the Southern clade is rather young, and the various populations it encompasses differentiated during the last million years (middle Pleistocene). Our results show that this clade has experienced a recent population expansion during the last 100,000–300,000 years, reaching a relatively final stable population size, exhibiting a dispersal route from Central Panama to South America (Fig. 5). Based on the genetic and phylogenetic results for this clade, we suggest that its dispersal towards South America and on both sides of the Northern Andes occurred very recently, likely as a consequence of climatic oscillations during glacial periods (Smith, Amei & Klicka, 2012). Under this scenario, there has been insufficient time for effective genetic differentiation to occur, and the phylogenetic reconstruction therefore failed to resolve the divergence detected in the spatial analysis.

In particular, the lack of differentiation between sequences from east of Northern Andes (in the Magdalena Valley) and H. tatayoi in Serranía del Perijá is remarkable: sequences were identical with extremely low divergence (F_ST = 0.0543, Table 1), supporting a recent dispersal through the Eastern Cordillera. In contrast, we detected a slight structure for the Ecuadorian record (Geneland; Fig. 2B) based on a single specimen (QCAZ-22303). We acknowledge that a more detailed genetic structure analysis of the Colombian and Ecuadorian populations is needed, including samples from the southern departments in Colombia, and the use of additional genetic markers (e.g., microsatellites or SNPs) is suggested.

The Northern clade

Unlike the Southern clade, the Northern clade shows significant genetic structure in four different lineages, ranging from western Chortis (Central America) to lowland forests in Veracruz and Guerrero (Mexico). One remarkable finding was the split between samples from either side of the MPJ fault system, where individuals separated by only 60 km exhibit great genetic distances (K2P = 4.5% for all three genes, Table 1), even reaching the limit of the 6% barcode gap for Neotropical amphibians in COI (5.2% for COI, Table S1.2; Lyra, Haddad & Azeredo-Espin, 2017). The calibration results showed that samples from these localities have been isolated since the Gelasian age (early Pleistocene; ∼2.19 Ma, HPD = 1.38–3.53). The MPJ fault system has been recognized as the main barrier to dispersal for multiple species ranging from the Maya to the Chortis blocks (Ornelas, Ruiz-Sánchez & Sosa, 2010; Barrera-Guzmán et al., 2012; Rovito et al., 2015). Our results indeed support the hypothesis that MPJ has effectively acted as a barrier for H. fleischmanni dispersal.

On the other hand, the Geneland and phylogenetic results showed a more complex scenario, in agreement with the presence of three well-defined lineages within the Northern clade, one of which is located in the north of the SMB, another southward the SMB, and the last in the Maya block (Fig. 2B). This genetic structure is very similar to that observed for the brush-finch Arremon brunneinucha, distributed in humid montane forests (Navarro-Sigüenza et al., 2008). Samples from the Gulf of Mexico and Pacific are clearly isolated by the highlands of the Sierra Madre Oriental and the Sierra Madre Occidental (Fig. 1), indicating a divergence pattern that is frequently detected for lowland species (Mulcahy, Morrill & Mendelson, 2006; Rivera-Ortíz et al., 2016) and species groups (Streicher et al., 2014; Palacios et al., 2016).

Within this clade, we found that samples from SMB do not group as a single lineage but instead display a paraphyletic position in relation to the lineage from the Maya region. These two regions present the lowest K2P among all comparisons performed in this analysis, except for samples separated by the WPI. Importantly, unlike the known impact of the WPI and CPI on lowland species, the effect of the IT has mostly been defined in relation to montane species (Bryson, García-Vázquez & Riddle, 2011; Jiménez & Ornelas, 2016), for which it represents a dispersal-limiting barrier. Currently, the IT is occupied by dry, scrubby coastal plains that are very different from the moist areas on either side (Rodríguez-Gómez, Gutiérrez-Rodríguez & Ornelas, 2013), so one plausible explanation for this finding is that the dry forests along the IT did not always act as a barrier to H. fleischmanni. Instead, considering that those clades were likely isolated between 1.76 to 2.19 Ma (HPD = 0.86–3.53), it is likely that the successive glacial cycles during the Pleistocene could have enabled the montane forests to reach lower elevations through the impact of climate change (Barber & Klicka, 2010), thus creating a temporal corridor and allowing the species to disperse through the region. Hence, for species like H. fleischmanni, the IT probably acted as both a corridor and a barrier, which is supported by our Bayesian spatiotemporal diffusion results (Fig. 5).

The Central clade

Our results revealed a Central clade without any deep geographic structure expanding across the HE and separated from the Southern clade around the end of the Pliocene (∼2.64 Ma). We did not find any evidence suggesting differentiation between samples from the Chortis and Chorotega regions, leading us to reject the hypothesis that the HE acted as a barrier. However, we must consider the small sample sizes from Honduras and Nicaragua (n = 10 samples), which likely limits detailed structural evaluation for this region. Most phylogeographic studies performed in the region known as nuclear Central America face similar problems, with limited or null sampling from northern Honduras (Castoe et al., 2003; Crawford & Smith, 2005; Mulcahy, Morrill & Mendelson, 2006) or sampling that is biased towards the dry forests of the Pacific coast (Parkinson, Zamudio & Greene, 2000; Hasbún et al., 2005; Vázquez-Miranda, Navarro-Sigüenza & Omland, 2009; Poelchau & Hamrick, 2011), where H. fleischmanni has not been recorded. Our Bayesian spatiotemporal diffusion results showed rapid dispersal from the Chorotega to Chortis blocks, with no apparent impact on the genetic structure of these populations (Fig. 5). Nevertheless, additional work is needed to confirm whether the main geographic features present in this region have driven the dispersal of low-mobility species such as H. fleischmanni in humid forests.

Phylogeographic patterns

The three main clades that we identified for H. fleischmanni show deep intraspecific divergence, with genetic distance values (K2P) greater than 3% (Table 2). Indeed, the landscape analysis, Bayesian spatiotemporal diffusion analyses, and estimated divergence times revealed interesting patterns that allowed us to reconstruct the historical biogeography of these frogs and identify the impact of different geographic barriers on the genetic structure and phylogeographic patterns of H. fleischmanni. Although the main phylogenetic topology and the three major clades were well supported in Bayesian analyses, the maximum likelihood phylogenetic reconstruction did not resolve these relationships (see Supplemental Information 1).

The Bayesian spatiotemporal diffusion analysis did not support the hypothesis proposed by Castroviejo-Fisher et al. (2014) of a South American origin and subsequent dispersal to Central America. In contrast, H. fleischmanni appears to have originated in the region encompassing the Chorotega and eastern Chortis elements. Interestingly, we found that H. fleischmanni has undergone two dispersal events: one southwards towards the Chocó region and one northwards, reaching the Maya region, followed by vicariance events driven by the effect of the Chortis highlands and the Talamanca Range. Considering that all species that are closely related to H. fleischmanni are endemic to South America and that the divergence times among the three clades are similar (i.e., the isolation of the Northern clade occurred ca. 3.40 Ma (HPD = 2.25–4.56 Ma), while that between the two other clades occurred 2.64 Ma (HPD = 1.50–3.68)), it is likely that the MRCA ancestor of H. fleischmanni and H. carlesvilai was located in South America around 7.65 Ma (CI = 5.93–9.63). Later, the ancestor of H. fleischmanni arrived in Central America during the Pliocene, soon after the closure of the Isthmus of Panama (Montes et al., 2015), whereas the current South American populations are descendants of a second migration from Central to South America. Accordingly, the dispersal-vicariance events among the main clades potentially occurred simultaneously or over a very short time, which might explain why the position of the clades and their internal structure were not consistent between the Bayesian and maximum likelihood approaches.

Regarding the vicariance events for the Central and Southern clades, multiple elements need to be revised. The samples from each cluster that are geographically closest are located on opposite sides of the Talamanca Cordillera in the Chorotega block. The time of divergence of the MRCA for these clades (3.28 Ma, HPD = 1.59–3.86) coincides with the estimated age of the intervening mountains (1–2 Ma; Denyer, Alvarado & Aguilar, 2000; Marshall et al., 2003, which are recognized as a main driver of speciation (Savage, 2002). The time of divergence also coincides with the rise of the sea level during the mid-late Pliocene (∼3.5–3 Ma), which generated a seaway, likely reinforcing the WPI break and therefore acting as a barrier across the Pacific region (Cronin & Dowsett, 1996; Bagley, Hickerson & Johnson, 2018). Hence, the central mountain ranges on Costa Rica and Panama and the eustatic sea levels around the WPI might have increased divergence, as documented for multiple spatial divergence patterns of amphibian species (Crawford, Bermingham & Polanía-S, 2007; Wang, Crawford & Bermingham, 2008; Bagley, Hickerson & Johnson, 2018).

The isolation of the Northern clade from the other two does not entirely correspond to our hypothesis of geographical barriers. Populations from both the Northern and Central clades are distributed throughout the Chortis block, indicating that the MPJ fault system was not the main driver of divergence between clades. On the other hand, our structure (Geneland) results suggest a frontier at the center of the Chortis block, near northeastern Honduras (Fig. 2B). Similar divergence patterns have been observed between two water-dependent subspecies of Caiman crocodilus (Venegas-Anaya et al., 2008), in agreement with the eastern limit of the Chortis highlands (Morrone, 2014; Townsend, 2014). Here, the complex topography resulting from multiple volcanic activities along the Chortis highlands during the last 2 Ma and the presence of dry habitats in the Pacific region (Savage, 2002) may have isolated the H. fleischmanni populations during the late Pliocene. This hypothesis is in agreement with the high species endemism recognized for the region (Anderson et al., 2010; Townsend et al., 2012), in which intensive study is required to evaluate the underestimated regional taxonomic diversity (Townsend & Wilson, 2016).

Taxonomic implications

Previous studies have suggested that H. fleischmanni is a paraphyletic species in relation to its sister species H. tatayoi (Castroviejo-Fisher et al., 2011; Delia, Bravo-Valencia & Warkentin, 2017). Here, we confirm the paraphyly of the species, as the H. tatayoi samples are grouped within the Southern clade, lacking significant genetic differences from the western Andes samples. Furthermore, our overall Bayesian topology agrees with the results obtained by Delia, Bravo-Valencia & Warkentin (2017) for 12S sequences. We identified three main isolated lineages with large genetic distances that likely include cryptic diversity within Central America. To confirm the existence of cryptic—and potential candidate—species, different lines of evidence must be obtained; in fact, we are conducting a follow-up study of the integrative taxonomy of this species complex that includes genetic, morphologic and acoustic data (Mendoza et al., 2018, unpublished data) to describe lineage divergence and determine the identities and geographic distributions of all valid species (Padial et al., 2010).