# Landscape Genetics and Species Delimitation in the Andean Palm Rocket Frog (Aromobatidae,Rheobates)

Authors: Gabrielle Genty; Carlos E. Guarnizo; Juan P. Ramírez; Lucas S. Barrientos; Andrew J. Crawford
Journal: Journal of Zoological Systematics and Evolutionary Research
 (2022)

---

## Abstract

The complex topography of the species-rich northern Andes creates heterogeneous environmental landscapes that are hypothesized to have promoted population fragmentation and diversification by processes such as vicariance or local adaptation. Previous phylogenetic work on the palm rocket frog (Anura: Aromobatidae:Rheobates spp.), endemic to midelevation forests of Colombia, suggested that valleys were important in promoting divergence between lineages. In this study, we first evaluated previous hypotheses of species-level diversity, then fitted an isolation-with-migration (IM) historical demographic model, and tested two landscape genetic models to explain genetic divergence within Rheobates: isolation by distance and isolation by environment. The data consisted of two mitochondrial and four nuclear genes from 24 samples covering most of the geographic range of the genus. Species delimitation by Bayesian Phylogenetics and Phylogeography recovered five highly divergent genetic lineages within Rheobates, among which few to no migrants are exchanged according to IM. We found that isolation by environment provided the only variable significantly correlated with genetic distances for both mitochondrial and nuclear genes, suggesting that local adaptation may have a role in driving the genetic divergence within this frog genus. Thus, genetic divergence in Rheobates may be driven more by variation among the local environments where these frogs live rather than by geographic distance.


---

## Body

## 1. Introduction

The South American landscape and environments have experienced enormous historical changes, such as major geological upheavals and environmental fluctuations related to glacial cycles [1], which have created a complex backdrop for species diversification. The tropical Andes is a global biodiversity hotspot, home to the highest density of species per unit area in the world [2, 3]. The Andean Cordillera has a dynamic history, progressively uplifting from south to north during the Late Cretaceous [4–6]. This orogenic uplift and its concomitant environmental heterogeneity have promoted the diversification of Andean species [7, 8].One way in which the orogeny of the Andes promotes species isolation and diversification is by creating a complex topography that may fragment lowland populations or limit dispersal across high elevation regions separated by valleys [9]. The Andean uplift created valleys that may have promoted the diversification of a variety of organisms, for example, the wax palm genus Ceroxylon [10], the Ithomiini tribe of butterflies [11, 12], the montane forest subspecies of the three-striped warbler bird, suggested to have occured by allopatric divergence [13], and the Adelomyia hummingbirds, whose divergence was coupled to Andean orogeny [14]. The rise of Andean peaks promoted lowland speciation in some groups, such as the Dendrocincla woodcreepers, through vicariance during the uplift [15], yet many lowland tetrapods appear to have been largely unaffected by Andean vicariance [16].Another way that Andean orogeny can promote diversification is by providing novel environments at different elevations, which in turn could promote local adaptation and ecological divergence [17]. Previous studies on birds support the hypothesis that environmental gradients along elevation belts have promoted trait evolution [18]. Hence, elevation can promote differentiation among populations but whether this leads to speciation is controversial. While the process of speciation along an environmental gradient is well documented in some systems [19], in the Neotropics, sister species of vertebrates rarely occur as elevational replacements of each other, as predicted by the gradient hypothesis of montane speciation [20]. A recent study of Andean plants also found little support for the gradient hypothesis [21].While environmental heterogeneity may be a key factor in explaining primary divergence, when invoking the environment as a barrier, there are two factors to consider: one is the geographic distance together with the variation in environments in intervening areas that isolates pairs of populations, while the second possibility is the distinctiveness of the two environments that may characterize the respective populations. The latter scenario is referred to as isolation by environment and suggests that populations inhabiting contrasting environments have limited dispersal due to either selection against immigrants or due to individuals preferring to remain in a particular environment [22], as has been found to occur in the European common frog, Rana temporaria, in Scotland [23]. In other words, the isolation by environment model is blind to the environmental landscape conditions found in between a pair of populations and therefore contrasts sharply with the more traditional isolation by distance or isolation by resistance models [24].Tropical amphibians are ideal organisms to evaluate how the geography and environment are associated with diversification in montane regions due to their restricted dispersal, strong site fidelity, and spatially isolated breeding habitat [25]. The palm rocket frog (Anura: Aromobatidae: Rheobates) is a useful system to investigate how Andean topography has promoted biological diversification. Rheobates is endemic to Colombia and is restricted to small streams and ponds at midelevation sites and is commonly found between 1000 m and 2400 m ([26, 27]; minimum and maximum elevations are 250 m and 2520 m, respectively). Rheobates is currently composed of two named species, R. palmatus [28], found in the Eastern and Central Cordilleras of the Andes, and R. pseudopalmatus [29], found only in the northern part of the Central Cordillera. However, the taxonomic status of these species merits further scrutiny [30] due to the apparent absence of diagnostic morphological characters. Furthermore, recent molecular phylogenetic studies have recovered R. palmatus from the southern part of the Central Cordillera (Departamentos Caldas and Tolima) as being more closely related to specimens from areas close to the type locality of R. pseudopalmatus than R. palmatus from the Eastern Cordillera [31, 32].In a previous phylogenetic study, using two mitochondrial genes and one nuclear gene, Muñoz-Ortiz et al. [32] found that Rheobates was composed of three highly divergent clades. The lineage located in the town of Santa María, Boyacá, on the eastern side of the Eastern Cordillera, was isolated from the rest of the genus, with the remaining populations forming two reciprocally monophyletic groups separated by the Magdalena River Valley Figure 1. The authors, however, did not employ statistically robust, multilocus species delimitation methods to evaluate the hypothesis of Rheobates being composed of more than one lineage. Moreover, analyses which utilize genealogically independent loci are needed to confirm phylogeographic histories [33]. In the present study, we combine three new nuclear markers with the previously published mitochondrial and nuclear DNA data and add two new samples from Santa María (Table S1). We use a multilocus, coalescent-based species delimitation approach to test the previous hypotheses that Rheobates consists of two species (as per current taxonomy), three species [32], or five species (potentially supported by phylogeographic structuring, see below). Because species delimitation methods assume no migration among clades, we used a model-based historical demographic analysis to evaluate migration rates within the genus Rheobates.Figure 1
Map of the Eastern and Central Andean cordilleras of Colombia indicating the sampling localities. The yellow line between the populations of Eastern North and Eastern Central cordilleras represents the Chicamocha Canyon. The star shows the location of Bogotá, the capital of Colombia. The colors of the sampling localities match the haplotypes in Figure2 and geographic clades in Figure 3.To test the hypotheses about the role of physical geography and environmental variation as potential drivers of genetic differentiation among populations ofRheobates, we implemented a landscape genetics approach. Specifically, we tested whether the spatial pattern of divergence in Rheobates fits a model of isolation by distance [34] or isolation by environment [22, 35]. The support for a pattern of isolation by distance would indicate that the geographic distance among populations explains genetic divergence within Rheobates [34]. The support for isolation by environment, however, may indicate local adaptation that limits gene flow because allochthonous alleles are less fit [35].


## 2. Methods

### 2.1. Sampling

Genomic DNA was extracted from 24 tissue samples ofRheobates obtained through museum donations from the Museo de Historia Natural C.J. Marinkelle at the Universidad de los Andes, Bogotá (see Table S1 in Supporting Information), using a DNeasy Blood & Tissue Kit (QIAGEN, Valencia, CA, USA) following the manufacturer’s protocol. The geographic sampling covered most of the known distribution of the genus (Figure 1). While Rheobates can be locally abundant, they can be very rare in some regions, such as the eastern slope of the Eastern Cordillera of the Andes. This geographically patchy distribution, coupled with a history of social conflict in Colombia, has impeded more dense sampling.For historical analyses, we included as outgroups two aromobatid frogs,Allobates femoralis and A. cf. juanii [31]. We used the previously published sequences [32] of two mitochondrial genes, cytochrome oxidase I (COI) [36] and 16S ribosomal RNA (16S). We also included sequences from the same specimens used in Muñoz-Ortiz et al. [32] and two new samples from the Santa María population. These new sequences were successfully amplified by polymerase chain reaction (PCR), making in total 22 samples for COI and 23 samples for 16S (see Table S2 and the references therein). Additionally, we included the previously published sequences [32] of the nuclear gene proopiomelanocortin (POMC) and sequenced new samples, making in total 22 samples for POMC. We also sequenced three new, rapidly evolving nuclear loci SF232, SF328, and SF412 [37], for a total of 22, 23, and 24 additional DNA sequence fragments, respectively (the differences in the numbers of genes or individuals sequenced across loci were due to difficulties with the PCR amplification). GenBank accession numbers, field and museum voucher codes, and GPS coordinates are provided in Supporting Table S1. See Supporting Table S2 for PCR primers and molecular protocols. The purified fragments were Sanger sequenced in both directions and assembled using Geneious version 6.0 [38]. Sequences were aligned assuming the Q-INS-i strategy for the 16S loci and the default strategy for the remaining loci, as implemented in the web server (http://mafft.cbrc.jp/alignment/server/) in MAFFT version 7 [39]. For protein-coding sequences, potential misalignments were checked against amino acid translations in Mesquite version 3.04 [40]. The final aligned matrix includes a total of 2800 characters (nucleotide sites) and 26 terminals.To visually inspect gene trees for signs of incongruence among loci in terms of genetic distances among frogs within and between sampling locations, we constructed a median-joining haplotype network [41] for each mitochondrial and nuclear gene independently with the software PopART [42]. Before we ran the haplotype network on the nuclear genes, we resolved the haplotypes from our diploid frogs using PHASE [43] and SeqPHASE [44], available at https://eeg-ebe.github.io/SeqPHASE/index.html. Mitochondrial DNA sequences showed no heterozygous sites, as expected for haploid genomes, so no phasing was necessary. Networks of genes concatenated by genome are shown in Figure 2.Figure 2
Haplotype networks ofRheobates spp. representing the (a) concatenated mitochondrial sequences (COI and 16S) and (b) concatenated nuclear sequences (POMC, SF232, SF328, and SF412). Circle area is proportional to the number of individuals in possession of the particular haplotype sequence. Each red number on a branch represents the mutational steps between haplotypes. Branches without a number have one mutational step between the haplotypes.


### 2.2. Phylogenetics and Species Delimitation

PartitionFinder 2 [45] was used to select the best-fitting partition scheme for each alignment according to the corrected Akaike information criterion (AICc). We defined partitions a priori by gene identity and by codon position for protein-coding genes (all genes except the ribosomal gene, 16S, and the intron marker, SF328). The best partition scheme for each subset is reported in Supporting Table S3.To characterize the evolutionary history of populations and potential species withinRheobates, we performed a maximum likelihood (ML) phylogenetic inference as implemented in RAxML version 8.0.19 [46]. We sought to test whether the inclusion of multiple new nuclear gene sequences and a species-tree approach would corroborate the phylogenetic results obtained previously (e.g., if previous results were biased towards the mitochondrial DNA (mtDNA) tree [32]. For this analysis, we used the standard GTRGAMMA substitution model for each of the partitions recovered by PartitionFinder 2 on the data matrix of all genes concatenated. We performed a RAxML with the -f a option, using 1000 replicate searches to find the optimal ML tree and autoMRE option, which automatically determines the sufficient number of bootstrap replicates to perform. The resulting tree was rooted using the two outgroups from the genus, Allobates (see above). RAxML analyses were also repeated separately for the concatenated mitochondrial genes and concatenated nuclear genes.To corroborate the results of the concatenated ML analysis with a multispecies coalescent approach, we also estimated a species tree using the software StarBEAST2 [47]. This software simultaneously estimates the gene trees for each locus, together with a species tree by fitting a multispecies coalescent process. The input file for StarBEAST2 was created with BEAUti 2 [48]. As an initial condition, we assigned our samples to the five major clades recovered in our concatenated ML tree (see below). The bModelTest package of StarBEAST2 was applied to explore the substitution model space according to the partitioning scheme described above. bModelTest allows the simultaneous estimation of both the substitution model parameters and the phylogeny [49]. Details about the selection of priors can be found in the Supporting Information. The sampling used two independent MCMC chains of 50 million generations each, saving one sampled tree per 1000 generations and with the first 10% of trees discarded as burn-in.The posterior sample of parameters from StarBEAST2 was checked for convergence and effective sample size (ESS>200) of the estimated parameters using Tracer version 1.7 [50]. The resulting species tree inference was visualized with DensiTree version 2.2.6 [51] and summarized with TreeAnnotator version 2.6.0 [48] as a majority 50% consensus rule tree.For the species delimitation analysis, we used Bayesian Phylogenetics and Phylogeography (BPP) version 3.2 [52]. BPP takes into consideration ancestral polymorphism and incomplete lineage sorting while implementing a multispecies coalescent model. As one assumption of BPP is the absence of gene flow between populations [52], we used the MCMC coalescent simulator, IMa2 [53], to estimate migration rates among population pairs under an isolation-with-migration model [54] to evaluate symmetric migration rates between sister clades within Rheobates (see Results). Runs were conducted with 100000 MCMC generations and a thinning of 100 with a burn-in of 10%. Following the recommendations in Yang [52], we assumed a population mutation rate θ=4Nμ=10 and symmetrical migration rates (m0⟶1=m1⟶0) between pairs of locations. Given the absence of gene flow, we proceeded with our species delimitation analysis using BPP. We based our guide tree on the topology that we previously estimated with RAxML, which clustered individuals into five reciprocally monophyletic groups, not counting the outgroup, three of which showed relatively low genetic divergence (see Results). All runs consisted of 100000 generations, sampled every tenth generation, with a burn-in period of 8000 steps. Further details of the BPP analysis can be found in Supporting Information.


### 2.3. Landscape Genetics

We evaluated the potential abiotic drivers of genetic divergence by estimating the relationship between geographic and environmental distances against genetic distances among all pairs of sampled individuals within the genusRheobates, i.e., regardless of whether samples were currently assigned to R. palmatus or R. pseudopalmatus (Table S1). We estimated genetic distances between individuals with one concatenated dataset for mitochondrial genes (since mtDNA genes behave as a single locus) and each nuclear gene independently. The model selection and genetic distances for each dataset were estimated in MEGA [55], in which HKY+Γ+I was the best model for the concatenated mitochondrial genes, JC for SF412, Kimura 2 parameter for SF328 and POMC, and Tamura 3 parameter for SF232.To spatially visualize the pairwise genetic distances withinRheobates, we used the software MAPI version 1.0.1 [56]. MAPI uses a spatial network in which samples are linked by ellipses and grids of hexagonal cells encompassing the study area. This method averages genetic distances into ellipses. The study area is then divided into hexagonal cells, and the ellipses that intersect the cell follow the assumption where the bigger the ellipse, the smaller its contribution to the cells [56]. Thus, more dissimilar hexagonal cells indicate that geographically closer individuals are more genetically different. To test whether observed spatial genetic discontinuities were greater than expected by chance, we generated null expectations through a randomization procedure.. We used a beta of 0.25, 10000 permutations, and an alpha set to 0.05. Geographic distances (to test for isolation by distance) were estimated between all pairwise localities using the software Geographical Distance Matrix Generator, version 1.2.3 [57], which takes into account the Earth’s curvature.Pairwise environmental distances (to test for isolation by environment) were estimated by extracting from each sampling site the value of each of the 19 bioclimatic variables available in the WorldClim database (http://www.bioclim.org) at 30 arcseconds resolution, using the function dist in the software R version 3.6.0 [58]. This function estimates Euclidean pairwise distances in the multidimensional space between rows on a multivariate matrix, which indicates how different the 19 bioclimatic variables are between sampling sites, where smaller distances indicate more similar abiotic environments. We did not test the effect of resistance distances, a common metric estimated in landscape genetics studies, because we do not know if we are testing a single species or multiple species within the genus Rheobates and the resistance matrices that needed to calculate these distances are species dependent.We tested which of the previous spatial variables (geographic and environmental) had the strongest association with genetic distances by using a multiple regression approach. Because they were on different scales, we standardized the estimated distances by subtracting the mean and dividing by the standard deviation. We then related the possible associations of the genetic distance with the four standardized matrices of abiotic variables using a multiple matrix regression with randomization analysis (MMRR) in R [59]. Given the nonindependence of the spatial variables, the MMRR analysis uses random sampling without replacement to generate null distributions. We conducted the MMRR analysis with 10000 permutations, using the function provided by Wang [22].


## 2.1. Sampling

Genomic DNA was extracted from 24 tissue samples ofRheobates obtained through museum donations from the Museo de Historia Natural C.J. Marinkelle at the Universidad de los Andes, Bogotá (see Table S1 in Supporting Information), using a DNeasy Blood & Tissue Kit (QIAGEN, Valencia, CA, USA) following the manufacturer’s protocol. The geographic sampling covered most of the known distribution of the genus (Figure 1). While Rheobates can be locally abundant, they can be very rare in some regions, such as the eastern slope of the Eastern Cordillera of the Andes. This geographically patchy distribution, coupled with a history of social conflict in Colombia, has impeded more dense sampling.For historical analyses, we included as outgroups two aromobatid frogs,Allobates femoralis and A. cf. juanii [31]. We used the previously published sequences [32] of two mitochondrial genes, cytochrome oxidase I (COI) [36] and 16S ribosomal RNA (16S). We also included sequences from the same specimens used in Muñoz-Ortiz et al. [32] and two new samples from the Santa María population. These new sequences were successfully amplified by polymerase chain reaction (PCR), making in total 22 samples for COI and 23 samples for 16S (see Table S2 and the references therein). Additionally, we included the previously published sequences [32] of the nuclear gene proopiomelanocortin (POMC) and sequenced new samples, making in total 22 samples for POMC. We also sequenced three new, rapidly evolving nuclear loci SF232, SF328, and SF412 [37], for a total of 22, 23, and 24 additional DNA sequence fragments, respectively (the differences in the numbers of genes or individuals sequenced across loci were due to difficulties with the PCR amplification). GenBank accession numbers, field and museum voucher codes, and GPS coordinates are provided in Supporting Table S1. See Supporting Table S2 for PCR primers and molecular protocols. The purified fragments were Sanger sequenced in both directions and assembled using Geneious version 6.0 [38]. Sequences were aligned assuming the Q-INS-i strategy for the 16S loci and the default strategy for the remaining loci, as implemented in the web server (http://mafft.cbrc.jp/alignment/server/) in MAFFT version 7 [39]. For protein-coding sequences, potential misalignments were checked against amino acid translations in Mesquite version 3.04 [40]. The final aligned matrix includes a total of 2800 characters (nucleotide sites) and 26 terminals.To visually inspect gene trees for signs of incongruence among loci in terms of genetic distances among frogs within and between sampling locations, we constructed a median-joining haplotype network [41] for each mitochondrial and nuclear gene independently with the software PopART [42]. Before we ran the haplotype network on the nuclear genes, we resolved the haplotypes from our diploid frogs using PHASE [43] and SeqPHASE [44], available at https://eeg-ebe.github.io/SeqPHASE/index.html. Mitochondrial DNA sequences showed no heterozygous sites, as expected for haploid genomes, so no phasing was necessary. Networks of genes concatenated by genome are shown in Figure 2.Figure 2
Haplotype networks ofRheobates spp. representing the (a) concatenated mitochondrial sequences (COI and 16S) and (b) concatenated nuclear sequences (POMC, SF232, SF328, and SF412). Circle area is proportional to the number of individuals in possession of the particular haplotype sequence. Each red number on a branch represents the mutational steps between haplotypes. Branches without a number have one mutational step between the haplotypes.


## 2.2. Phylogenetics and Species Delimitation

PartitionFinder 2 [45] was used to select the best-fitting partition scheme for each alignment according to the corrected Akaike information criterion (AICc). We defined partitions a priori by gene identity and by codon position for protein-coding genes (all genes except the ribosomal gene, 16S, and the intron marker, SF328). The best partition scheme for each subset is reported in Supporting Table S3.To characterize the evolutionary history of populations and potential species withinRheobates, we performed a maximum likelihood (ML) phylogenetic inference as implemented in RAxML version 8.0.19 [46]. We sought to test whether the inclusion of multiple new nuclear gene sequences and a species-tree approach would corroborate the phylogenetic results obtained previously (e.g., if previous results were biased towards the mitochondrial DNA (mtDNA) tree [32]. For this analysis, we used the standard GTRGAMMA substitution model for each of the partitions recovered by PartitionFinder 2 on the data matrix of all genes concatenated. We performed a RAxML with the -f a option, using 1000 replicate searches to find the optimal ML tree and autoMRE option, which automatically determines the sufficient number of bootstrap replicates to perform. The resulting tree was rooted using the two outgroups from the genus, Allobates (see above). RAxML analyses were also repeated separately for the concatenated mitochondrial genes and concatenated nuclear genes.To corroborate the results of the concatenated ML analysis with a multispecies coalescent approach, we also estimated a species tree using the software StarBEAST2 [47]. This software simultaneously estimates the gene trees for each locus, together with a species tree by fitting a multispecies coalescent process. The input file for StarBEAST2 was created with BEAUti 2 [48]. As an initial condition, we assigned our samples to the five major clades recovered in our concatenated ML tree (see below). The bModelTest package of StarBEAST2 was applied to explore the substitution model space according to the partitioning scheme described above. bModelTest allows the simultaneous estimation of both the substitution model parameters and the phylogeny [49]. Details about the selection of priors can be found in the Supporting Information. The sampling used two independent MCMC chains of 50 million generations each, saving one sampled tree per 1000 generations and with the first 10% of trees discarded as burn-in.The posterior sample of parameters from StarBEAST2 was checked for convergence and effective sample size (ESS>200) of the estimated parameters using Tracer version 1.7 [50]. The resulting species tree inference was visualized with DensiTree version 2.2.6 [51] and summarized with TreeAnnotator version 2.6.0 [48] as a majority 50% consensus rule tree.For the species delimitation analysis, we used Bayesian Phylogenetics and Phylogeography (BPP) version 3.2 [52]. BPP takes into consideration ancestral polymorphism and incomplete lineage sorting while implementing a multispecies coalescent model. As one assumption of BPP is the absence of gene flow between populations [52], we used the MCMC coalescent simulator, IMa2 [53], to estimate migration rates among population pairs under an isolation-with-migration model [54] to evaluate symmetric migration rates between sister clades within Rheobates (see Results). Runs were conducted with 100000 MCMC generations and a thinning of 100 with a burn-in of 10%. Following the recommendations in Yang [52], we assumed a population mutation rate θ=4Nμ=10 and symmetrical migration rates (m0⟶1=m1⟶0) between pairs of locations. Given the absence of gene flow, we proceeded with our species delimitation analysis using BPP. We based our guide tree on the topology that we previously estimated with RAxML, which clustered individuals into five reciprocally monophyletic groups, not counting the outgroup, three of which showed relatively low genetic divergence (see Results). All runs consisted of 100000 generations, sampled every tenth generation, with a burn-in period of 8000 steps. Further details of the BPP analysis can be found in Supporting Information.


## 2.3. Landscape Genetics

We evaluated the potential abiotic drivers of genetic divergence by estimating the relationship between geographic and environmental distances against genetic distances among all pairs of sampled individuals within the genusRheobates, i.e., regardless of whether samples were currently assigned to R. palmatus or R. pseudopalmatus (Table S1). We estimated genetic distances between individuals with one concatenated dataset for mitochondrial genes (since mtDNA genes behave as a single locus) and each nuclear gene independently. The model selection and genetic distances for each dataset were estimated in MEGA [55], in which HKY+Γ+I was the best model for the concatenated mitochondrial genes, JC for SF412, Kimura 2 parameter for SF328 and POMC, and Tamura 3 parameter for SF232.To spatially visualize the pairwise genetic distances withinRheobates, we used the software MAPI version 1.0.1 [56]. MAPI uses a spatial network in which samples are linked by ellipses and grids of hexagonal cells encompassing the study area. This method averages genetic distances into ellipses. The study area is then divided into hexagonal cells, and the ellipses that intersect the cell follow the assumption where the bigger the ellipse, the smaller its contribution to the cells [56]. Thus, more dissimilar hexagonal cells indicate that geographically closer individuals are more genetically different. To test whether observed spatial genetic discontinuities were greater than expected by chance, we generated null expectations through a randomization procedure.. We used a beta of 0.25, 10000 permutations, and an alpha set to 0.05. Geographic distances (to test for isolation by distance) were estimated between all pairwise localities using the software Geographical Distance Matrix Generator, version 1.2.3 [57], which takes into account the Earth’s curvature.Pairwise environmental distances (to test for isolation by environment) were estimated by extracting from each sampling site the value of each of the 19 bioclimatic variables available in the WorldClim database (http://www.bioclim.org) at 30 arcseconds resolution, using the function dist in the software R version 3.6.0 [58]. This function estimates Euclidean pairwise distances in the multidimensional space between rows on a multivariate matrix, which indicates how different the 19 bioclimatic variables are between sampling sites, where smaller distances indicate more similar abiotic environments. We did not test the effect of resistance distances, a common metric estimated in landscape genetics studies, because we do not know if we are testing a single species or multiple species within the genus Rheobates and the resistance matrices that needed to calculate these distances are species dependent.We tested which of the previous spatial variables (geographic and environmental) had the strongest association with genetic distances by using a multiple regression approach. Because they were on different scales, we standardized the estimated distances by subtracting the mean and dividing by the standard deviation. We then related the possible associations of the genetic distance with the four standardized matrices of abiotic variables using a multiple matrix regression with randomization analysis (MMRR) in R [59]. Given the nonindependence of the spatial variables, the MMRR analysis uses random sampling without replacement to generate null distributions. We conducted the MMRR analysis with 10000 permutations, using the function provided by Wang [22].


## 3. Results

### 3.1. Phylogenetics and Species Delimitation

The haplotype networks inferred for each individual gene showed no major conflicts and agreed that a locality from the eastern side of the Eastern Cordillera (Santa María) was the most divergent (Figure2). The concatenated ML tree of the combined data (Figure 3), as well as with mitochondrial and nuclear genes concatenated separately (Supporting Figures S1 and S2), identified five reciprocally monophyletic lineages within the genus Rheobates and suggested that the lineage from Santa María was sister to the rest of Rheobates, although with no meaningful bootstrap support. The inferred sister clade to Santa María was composed of four well-separated and largely allopatric lineages identifiable by their geographic distribution. The Central Cordillera clade was separated from those of the Eastern Cordillera, and the later was divided into three clades. Two of these three clades along the western flank of the Eastern Cordillera are found north versus south of the Chicamocha Canyon, and the third is located further south (Figure 3). Statistical support for relationships among the two lineages of the western side of the Eastern Cordillera received low bootstrap support of just 54% (Figure 3). These phylogenetic results correspond well with those of Muñoz-Ortiz et al. [32], although nodal support tended to be higher in the previous study, even though it was based on a reduced number of genes. Results from the Bayesian multispecies coalescent approach show even weaker support for relationships among lineages within Rheobates (displayed graphically in the DensiTree of Figure 3 and the species tree in Supporting Figure S3, which provides statistical support values). Historical demographic analyses of coalescent models of isolation-with-migration using IMa2 revealed for most pairwise comparisons among the five phylogenetic lineages, a migration rate (Nem) indistinguishable from zero. The highest of these low-migration rate estimates was between Eastern South and Eastern Central Cordillera with Nem=0.178 (95% confidence interval of 0 to 0.716) migrants per generation (Table 1). Low migration rates (<0.2 migrants per generation) likely do not impede species delimitation [60].Figure 3
DensiTree based on posterior samples of species trees showing the five major clades named according to geographic range and displayed underneath the maximum likelihood tree inferred for the concatenated mitochondrial and nuclear genes using the software RAxML with the bootstrap support values for each node. The corresponding species tree is provided in Figure S3.Table 1
The effective number of migrants (Nem) between pairs of population within Rheobates. 95 % credible intervals (CI) were calculated from posterior distribution of parameter estimates assuming an isolation-with-migration model of historical demography in the Bayesian MCMC software, IMa2. Asymmetric migration rate estimates are shown here for migration from population 1 to 2 and from 2 to 1.
Population pairsMigration (95% CI)1–>22–>1Santa María vs. Central Cordillera0.078 (0.001–0.278)0.000 (0–0.131)Central Cordillera vs. Eastern Central0.024 (0–0.500)0.000 (0–0.292)Eastern Central vs. Eastern South0.000 (0–0.336)0.178 (0–0.716)Eastern Central vs. Eastern North0.087 (0–0.647)0.000 (0–0.774)Eastern North vs. Eastern South0.000 (0–0.578)0.000 (0–0.601)The BPP method of species delimitation showed that the five groups identified here have a probability of 99.98% of representing distinct, unconfirmed candidate species, i.e., Santa María, the Central Cordillera, plus the three clades distributed along the western flank of the Eastern Cordillera that formed potentially allopatric replacements of each other along a north-south transect (Figure3).


### 3.2. Landscape Genetics

The visualization made with MAPI (Figure4) indicated that Rheobates is characterized by strong spatial variation in the pairwise genetic distances among individuals, displaying several genetic discontinuities not expected by chance, and a gradual pure isolation by distance model was rejected by the mtDNA and by each nuclear fragment, with the mtDNA data displaying the strongest genetic differences. In all cases, the south-east localities (surrounding the southern part of the Eastern Cordillera) displayed higher interindividual dissimilarity values compared with the south-west, north-east, and north-west regions (Figure 4).Figure 4
MAPI spatial structure outputs. Each plot corresponds to the minimum convex polygon that delineatesRheobates sampled localities depicted in the (g) map. Polygons with darker contours correspond to areas with significantly higher interindividual dissimilarity than expected by chance. Colors towards blue indicate higher genetic dissimilarity and towards red lower genetic dissimilarity. Genetic distances were estimated with the (a) concatenated mitochondrial dataset and (b–e) each nuclear gene. (f) Current annual temperature (bio01) as a reference point.The MMRR full model revealed that the geographic variables included in this study explained 63% of the genetic variation of the concatenated mitochondrial genes inRheobates (r2=0.6334), 68% of the variation of POMC (r2=0.6825), 40% of SF232 (r2=0.3984), 20% of SF328 (r2=0.2027), and 48% of SF412 (r2=0.4805). Regarding the effect of individual variables, we found that for every gene fragment, the beta coefficients were at least twice as high for environmental distances than for geographic distances, with the former ranging between 0.37 and 0.7 (Table 2). While geographic distances also showed a significant association with genetic distances in all gene fragments except for SF412, again, the effect size was half or much less (lower betas) (Table 2). Wang et al. [22] suggest that models with greater regression coefficients are explaining more of the variation in genetic distances, and thus, we give more importance to isolation by environment, even though both environment and isolation by distance are statistically significant (Supporting Figure S4).Table 2
Regression coefficient (β), t statistic, and significance (p value) of multiple regression with randomization (MMRR) analysis on the association between genetic distances and geographic and environmental distances, for each DNA fragment in Rheobates.
Gene fragmentGeographic variableBeta coefficientt statisticp valueMitochondrial (concatenatedCOI+16S)Geographic distance0.32008.57230.00009Environmental distances0.628416.83340.00009POMC (nuclear)Geographic distance0.150133.918670.00019Environmental distances0.7963920.03000.00009SF232 (nuclear)Geographic distance0.224924.32780.00009Environmental distances0.496639.30580.00009SF328 (nuclear)Geographic distance0.1384192.470140.01379Environmental distances0.3743316.730000.00009SF412 (nuclear)Geographic distance−0.02130−0.479280.63713Environmental distances0.7001315.75330.00009


## 3.1. Phylogenetics and Species Delimitation

The haplotype networks inferred for each individual gene showed no major conflicts and agreed that a locality from the eastern side of the Eastern Cordillera (Santa María) was the most divergent (Figure2). The concatenated ML tree of the combined data (Figure 3), as well as with mitochondrial and nuclear genes concatenated separately (Supporting Figures S1 and S2), identified five reciprocally monophyletic lineages within the genus Rheobates and suggested that the lineage from Santa María was sister to the rest of Rheobates, although with no meaningful bootstrap support. The inferred sister clade to Santa María was composed of four well-separated and largely allopatric lineages identifiable by their geographic distribution. The Central Cordillera clade was separated from those of the Eastern Cordillera, and the later was divided into three clades. Two of these three clades along the western flank of the Eastern Cordillera are found north versus south of the Chicamocha Canyon, and the third is located further south (Figure 3). Statistical support for relationships among the two lineages of the western side of the Eastern Cordillera received low bootstrap support of just 54% (Figure 3). These phylogenetic results correspond well with those of Muñoz-Ortiz et al. [32], although nodal support tended to be higher in the previous study, even though it was based on a reduced number of genes. Results from the Bayesian multispecies coalescent approach show even weaker support for relationships among lineages within Rheobates (displayed graphically in the DensiTree of Figure 3 and the species tree in Supporting Figure S3, which provides statistical support values). Historical demographic analyses of coalescent models of isolation-with-migration using IMa2 revealed for most pairwise comparisons among the five phylogenetic lineages, a migration rate (Nem) indistinguishable from zero. The highest of these low-migration rate estimates was between Eastern South and Eastern Central Cordillera with Nem=0.178 (95% confidence interval of 0 to 0.716) migrants per generation (Table 1). Low migration rates (<0.2 migrants per generation) likely do not impede species delimitation [60].Figure 3
DensiTree based on posterior samples of species trees showing the five major clades named according to geographic range and displayed underneath the maximum likelihood tree inferred for the concatenated mitochondrial and nuclear genes using the software RAxML with the bootstrap support values for each node. The corresponding species tree is provided in Figure S3.Table 1
The effective number of migrants (Nem) between pairs of population within Rheobates. 95 % credible intervals (CI) were calculated from posterior distribution of parameter estimates assuming an isolation-with-migration model of historical demography in the Bayesian MCMC software, IMa2. Asymmetric migration rate estimates are shown here for migration from population 1 to 2 and from 2 to 1.
Population pairsMigration (95% CI)1–>22–>1Santa María vs. Central Cordillera0.078 (0.001–0.278)0.000 (0–0.131)Central Cordillera vs. Eastern Central0.024 (0–0.500)0.000 (0–0.292)Eastern Central vs. Eastern South0.000 (0–0.336)0.178 (0–0.716)Eastern Central vs. Eastern North0.087 (0–0.647)0.000 (0–0.774)Eastern North vs. Eastern South0.000 (0–0.578)0.000 (0–0.601)The BPP method of species delimitation showed that the five groups identified here have a probability of 99.98% of representing distinct, unconfirmed candidate species, i.e., Santa María, the Central Cordillera, plus the three clades distributed along the western flank of the Eastern Cordillera that formed potentially allopatric replacements of each other along a north-south transect (Figure3).


## 3.2. Landscape Genetics

The visualization made with MAPI (Figure4) indicated that Rheobates is characterized by strong spatial variation in the pairwise genetic distances among individuals, displaying several genetic discontinuities not expected by chance, and a gradual pure isolation by distance model was rejected by the mtDNA and by each nuclear fragment, with the mtDNA data displaying the strongest genetic differences. In all cases, the south-east localities (surrounding the southern part of the Eastern Cordillera) displayed higher interindividual dissimilarity values compared with the south-west, north-east, and north-west regions (Figure 4).Figure 4
MAPI spatial structure outputs. Each plot corresponds to the minimum convex polygon that delineatesRheobates sampled localities depicted in the (g) map. Polygons with darker contours correspond to areas with significantly higher interindividual dissimilarity than expected by chance. Colors towards blue indicate higher genetic dissimilarity and towards red lower genetic dissimilarity. Genetic distances were estimated with the (a) concatenated mitochondrial dataset and (b–e) each nuclear gene. (f) Current annual temperature (bio01) as a reference point.The MMRR full model revealed that the geographic variables included in this study explained 63% of the genetic variation of the concatenated mitochondrial genes inRheobates (r2=0.6334), 68% of the variation of POMC (r2=0.6825), 40% of SF232 (r2=0.3984), 20% of SF328 (r2=0.2027), and 48% of SF412 (r2=0.4805). Regarding the effect of individual variables, we found that for every gene fragment, the beta coefficients were at least twice as high for environmental distances than for geographic distances, with the former ranging between 0.37 and 0.7 (Table 2). While geographic distances also showed a significant association with genetic distances in all gene fragments except for SF412, again, the effect size was half or much less (lower betas) (Table 2). Wang et al. [22] suggest that models with greater regression coefficients are explaining more of the variation in genetic distances, and thus, we give more importance to isolation by environment, even though both environment and isolation by distance are statistically significant (Supporting Figure S4).Table 2
Regression coefficient (β), t statistic, and significance (p value) of multiple regression with randomization (MMRR) analysis on the association between genetic distances and geographic and environmental distances, for each DNA fragment in Rheobates.
Gene fragmentGeographic variableBeta coefficientt statisticp valueMitochondrial (concatenatedCOI+16S)Geographic distance0.32008.57230.00009Environmental distances0.628416.83340.00009POMC (nuclear)Geographic distance0.150133.918670.00019Environmental distances0.7963920.03000.00009SF232 (nuclear)Geographic distance0.224924.32780.00009Environmental distances0.496639.30580.00009SF328 (nuclear)Geographic distance0.1384192.470140.01379Environmental distances0.3743316.730000.00009SF412 (nuclear)Geographic distance−0.02130−0.479280.63713Environmental distances0.7001315.75330.00009


## 4. Discussion

Rheobates is a montane genus of frogs distributed across two of the three main mountain ranges of the Colombian Andes. In this study, we explored how these complex Neotropical landscapes may drive genetic divergence within the genus. We find evidence that supports more lineages (i.e., undescribed candidate species) than what were previously expected. Moreover, our strongest genetic signal suggests that individuals from sites with more similar environments are more genetically similar, which supports an isolation by environment scenario in Rheobates.In terms of the lineages within the genus, species delimitation using BPP supports with high confidence thatRheobates is composed of five monophyletic lineages that may correspond to distinct species, contrasting with the three species suggested by Muñoz et al. [32]. Supporting the high genetic divergence found in the ML analysis (Figure 3), these five lineages displayed migration rate point estimates of Nem<0.006 migrants per generations that were statistically indistinguishable from zero, further supporting their status as unnamed species. Preliminary analyses performed by one of us (JPR) suggest that the external morphology of adults cannot differentiate specimens from the five major clades, despite the strong genetic divergences uncovered here (unpublished data). Sukumaran and Knowles [61] argue that BPP is an inadequate tool to delimit species, as it tends to overestimate species diversity. Although BPP has been shown previously to encounter problems with delimiting recently diverged species, especially those affected by geographic distances and large population sizes [62], BPP is still one of the best tools for identifying lineages of interest. In the case of Rheobates, we have a recently diverged species affected by geographic distances, so further studies based on additional sources of evidence (e.g., more genetic loci, larval traits, or male advertisement calls) are needed to test the hypothesis that these five unconfirmed candidate species might be confirmed and perhaps described [63]. Until such time, we agree that BPP has not necessarily revealed species but distinct genetic lineages that certainly warrant further study.At first glance, the geographic position of each sampling locality seems to be related to its degree of divergence, suggesting that the geographic distance is the main mechanism that promotes divergence in the genus. For example, the most southeastern (Santa María) and most northwestern (Central Cordillera) localities are the most divergent ones. If, however, the geographic distance is the main factor promoting genetic divergence inRheobates, we would have found strong support for isolation by distance. Instead, the strongest support was for isolation by environment. Thus, rather than (or in addition to) the intervening areas between populations (geographic distance), the local environment seems to be promoting divergence independently of geographic distance between populations.Palm rocket frogs show relatively high tolerance of considerable anthropic alterations of their environment [64, 65]. We are not aware of any studies looking explicitly at geographic variation in the autecology of Rheobates, and we are thus unsure of exactly which abiotic differences among regions might be driving a pattern of isolation by environment. Based on a study of a few localities in the Eastern Cordillera, Rheobates apparently reproduces all year long [66], making potential differences in phenology between sites difficult to observe. Our estimates of isolation by environment are derived from macroclimatic variables obtained from WorldClim. Finding such a strong association between macroclimate and genetic distances is somewhat unexpected, since these frogs are tiny and we might expect specific microclimates to correlate better with genetic distances. Perhaps, macroclimate variables average out fine-scale or short-term variability in microclimatic characteristics, so they are matching the genetic distances that are built up over the long term. An interesting follow-up study would be to look at how microclimates might predict genetic variation at the genomic level in Rheobates, given the previous reports of how microclimate variation predicts genetic variation in organisms such as epiphytic bryophytes [67], forest birds [68], and frogs of the genus Pristimantis [69].We cannot reject outright the role of geographic barriers in promoting genetic differentiation among populations, especially considering that topographic barriers and environmental barriers could be correlated. For example, the eastern slope of the Eastern Cordillera is more humid than the western slope [70]. This is mostly due to the rain shadow effect caused by the trade winds that come from the east [71]. In this case, the geographic variation in environment, i.e., higher humidity along the eastern slope relative to the western slope, is due to a potential geographic barrier, the Eastern Cordillera (but see [16]). Assuming uniformitarianism that the same processes have been acting on Rheobates across its 35 million year history (the estimated stem age reported by Muñoz-Ortiz et al. [32]), a potential scenario would be that initial divergences are driven by local adaptation to environmental variation while obvious landscape features may arise secondarily, such as through the formation of valleys or mountain peaks. Such a hypothesis would fit with what has been observed in Neotropical birds, where the organism-environment interaction plays a more important role than landscape features [72].Our multiple regression analysis revealed that isolation by distance was supported by the mitochondrial loci but not by the nuclear loci, while isolation by environment was supported by the nuclear loci but not the mitochondrial loci. Evolutionary processes (e.g., natural selection and genetic drift) can act on mitochondrial and nuclear genomes differently. Discordance between mitochondrial and nuclear gene regions has been frequently encountered among animals including amphibians [73, 74]. The discordances among loci in phylogeographic structure may result from incomplete lineage sorting [73], different effective population sizes [75], natural selection [76], life-history traits such as sex-biased dispersal [77], coalescent variance [78], selective sweep on a specific area of the DNA [79], or historic isolation and secondary contact [80]. Future genomic level analyses and the use of microclimatic variables may help resolve the effect that each of these geographic features is having on the diversification of this genus of frogs.Historically, much attention has been focused on the effects of vicariance and geographic distance as main factors promoting diversification in Neotropical mountains and the case of the palm rocket frog would seem to be a prime example [32]. In the present study, however, we demonstrate that, even in Rheobates, environmental factors and ecological adaptation, as revealed by the isolation by environment model, could be more important than landscape features and vicariance [72] in promoting diversification. Our results thus lend support to the growing idea in phylogeographic studies that the geographic component of genetic diversity must be driven by the interaction between an organism and its environment, an approach referred to as trait-based phylogeography [81]. The challenge going forward will be to associate particular traits with reduced migration rates across a heterogeneous environment, such as through landscape genomics [82]. More data on the habitats, life history, and genomes of Neotropical frogs will revolutionize our understanding of the origins of diversity.



---

*Source: 6774225-2022-02-09.xml*