Biodiversity is declining due to human activity throughout the globe, resulting in the degradation of ecosystem services and the loss of important components of our planet's natural heritage (Wilson, 1988; Chapin et al., 2000; Chivian and Bernstein, 2008; Newbold et al., 2015). However, it is difficult to address the biodiversity crisis when the diversity, geographic boundaries, and ecological functions of species remain unknown. Biodiversity and threats to biodiversity are heterogeneously distributed, leading some to argue that conservation efforts should focus strategically on regions where results per unit effort can be maximized. Areas that are highly threatened and have high levels of endemic species are known as biodiversity hotspots (Myers et al., 2000). One such biodiversity hotspot is the North Andean Pacific Slopes–Rio Atrato ecoregion of South America (FEOW, 2021), which includes the Chocó-Darien (Gómez N. et al., 2013; Rincon-Sandoval et al., 2019). Within this ecoregion, western Ecuador is particularly diverse and imperiled. Ecuador is a small country that harbors tremendous biodiversity, much of which is endemic. Rates of endemism are particularly high in western Ecuador because of its geographical restriction between the Andes Mountains and Pacific Ocean, disjunct drainage basins, heterogeneous topography, and climatic variability (Albert and Reis, 2011; Barriga, 2012; Jiménez-Prado et al., 2015). Studies of the evolutionary relationships and population genetic structure of species in this region are scarce despite the widespread environmental degradation (Damanik-Ambarita et al., 2016; Ribeiro et al., 2017).

The Rainbow Characin genus Rhoadsia (Fig. 1) is a charismatic freshwater fish that is endemic to the Chocó, with a range restricted to western Ecuador and northwestern Peru. Rhoadsia has long been grouped with only three other genera, Carlana, Nematocharax, and Parastremma, as part of the small subfamily Rhoadsiinae, which is characterized by having a single tooth series in the premaxilla when young and a double series when adults, with two conical teeth comprising the outer series and five multicuspid teeth comprising the inner series (Cardoso, 2003). However, current integrative DNA + morphology-based classification schemes place these genera within the much larger subfamily Stethaprioninae (Mirande, 2019).

Fig. 1

Live photos of Rhoadsia spp. (A) Male Rhoadsia minor, ROM 93672, from the Silanchi River, Esmeraldas River drainage; (B) female Rhoadsia minor, ROM 93696, from the Leilia River, Esmeraldas River drainage; and (C) male Rhoadsia altipinna, ROM 93791, from the Santa Rosa River. All scale bars = 1 cm.

Rhoadsia currently includes two species: Rhoadsia altipinna (Fowler, 1911) and Rhoadsia minor (Eigenmann and Henn, 1914; Fig. 1). Rhoadsia minor was described from specimens collected at high elevations in the Esmeraldas River drainage near Mindo, Esmeraldas Province, northwestern Ecuador (Barriga, 2012; Jiménez-Prado et al., 2015), while R. altipinna was originally described from specimens collected at lower elevations in the Guayas River basin of southwestern Ecuador, and occurs south to the coastal drainages of northwestern Peru (Barriga, 2012; Ortega et al., 2012; Jiménez-Prado et al., 2015). Populations of Rhoadsia are very similar morphologically, having a distinctive black spot on their flanks (Fig. 1), relatively deep bodies, and exhibiting strong sexual dimorphism. Males in both species are much larger, have bright red, yellow, and blue iridescent colors, and have larger dorsal, anal, and pectoral fins than females (Jiménez-Prado et al., 2015; Fig. 1).

Some authors have questioned the validity of Rhoadsia minor (Böhlke, 1958; Géry, 1977). For example, Géry (1977) suggested that R. altipinna and R. minor may constitute a single species separated into morphologically variable populations adapting to local conditions at different elevations. The main distinction between the type specimens of both nominal species is body depth, with types of R. minor being more streamlined and types of R. altipinna having a deeper body. However, body depth and vertebral number are now known to vary with elevation in river basins throughout the range of the genus (Aguirre et al., 2016, 2019; Malato et al., 2017), such that low-elevation populations from the northern Esmeraldas drainage inhabited by populations presumed to be R. minor are morphologically similar to southern, low-elevation populations presumed to be R. altipinna (Böhlke, 1958). Previous genetic studies based on limited sampling using the relatively conservative mitochondrial cytochrome oxidase I (COI) gene showed little genetic divergence among populations of Rhoadsia and no reciprocal monophyly of haplotypes (Aguirre et al., 2016; Malato et al., 2017).

The main goal of this study is to conduct a phylogeographic analysis to clarify the taxonomic status of the currently recognized species of Rhoadsia and understand how populations of these species are evolutionarily related in western Ecuador. We also assess whether patterns of genetic divergence in Rhoadsia are concordant with previous hypotheses of the aquatic zoogeographic regions in western Ecuador (Barriga, 2012). We build on previous efforts by expanding the geographic sampling for this genus, sequencing the COI gene for specimens from these new collections, and sequencing the more rapidly evolving mitochondrial cytochrome b (Cyt-b) gene and 12 microsatellites isolated by Loh et al. (2014).

MATERIALS AND METHODS

Sample collection.—A combination of new and existing samples of Rhoadsia from eight drainages in western Ecuador collected at elevations between 20 and 1260 m were included in this study (Table 1, Fig. 2). Tissue samples collected in 2008 and 2014 by Malato et al. (2017) were included, as were samples collected in 2012 and available at the Royal Ontario Museum (ROM; Supplemental Table 1; see Data Accessibility), and samples of up to 24 individuals per population collected in the summer of 2017 (Table 1). Sites sampled in 2017 were shallow (∼500 mm depth) running streams with bottoms composed mainly of cobble, boulder, gravel, and sand. For these samples, a Smith Root model LR24 electrofishing backpack was used for capturing the specimens along with sienes and dip nets. The general outputs set in the equipment were 400 to 850 volts, a frequency of 25–30 Hz, and a standard pulse type. The setting options of the equipment were adjusted depending on the environmental conditions of the sampling area. A cast net was used for deeper water sites or where electrofishing could not be used due to high conductivity. Specimens were euthanized in MS-222 or ice water. Fin clips were stored in 95% ethanol, while the rest of each specimen was fixed in 10% formaldehyde for ≥ 48 hours, washed, and transferred to 70% ethanol for long term storage.

Fig. 2

Map of Ecuador showing the sites and drainages in western Ecuador sampled for this study. The drainages are colored in gray and labeled in bold. The sites are represented with different symbols referring to their respective drainage. Streams related to the sites sampled are highlighted. Meaning of site codes can be seen in Table 1. Scale of map, 1:390000.

Table 1

Sampling sites and sample sizes for each molecular marker. N (total number of individuals), Cyt-b (cytochrome b), COI (cytochrome oxidase I), Mstl (microsatellites), ROM (Royal Ontario Museum), STRI (Smithsonian Tropical Research Institute Panama), CR (Costa Rica), Nic (Nicaragua), Pan (Panama). * Sequences retrieved from GenBank. Further details of samples from museums provided in Supplemental Table 1 (see Data Accessibility). ** Parastremma sadina (Par) and Carlana eigenmanni (Car), species used as an outgroup for the phylogeographic analysis.

Samples of two closely related genera in the subfamily Rhoadsiinae were included in the analyses as outgroups. Tissue samples of Parastremma sadina from Colombia were obtained from the ROM and of Carlana eigenmanni from the Smithsonian Tropical Research Institute in Panama (Supplemental Table 1; see Data Accessibility).

DNA extraction and mitochondrial gene amplification (Cyt-b and COI).—DNA extraction was performed using a phenol: chloroform:isoamyl alcohol method following Aguirre et al. (2016). Up to 24 specimens per site were sampled for Cyt-b. COI sequences from Malato et al. (2017) were retrieved from GenBank and included in this study (accession KY440344–KY440356). Up to eight individuals per sample from the new sites collected in this study were sequenced for the COI gene (Table 1). Smaller sample sizes were used for COI because previous work indicated that this gene varies relatively little among populations of Rhoadsia (Malato et al., 2017).

The polymerase chain reaction (PCR) was used to amplify the Cyt-b and COI genes. Two sets of primers were used to amplify Cyt-b. The primers CB1-5 (forward), 5′–CCATCCAACATCTCAGCATGATGAAA–3′ and CB3-3 (reverse) 5′–GGCAAATAGGAAATATCATTC–3′ (Palumbi, 1996) were used to amplify some samples of Rhoadsia spp. and all samples of Carlana eigenmanni and Parastremma sadina. A more specific set of primers for Rhoadsia spp. was created using invariant sites of the DNA sequences obtained with the first set of primers: Rh_Cytb_F2 (forward), 5′–GACATTTCTTTAGCTTTTTCC–3′ and Rh_Cytb_R1 (reverse), 5′–GGTTTAATATGAGGTGGTGT–3′. PCR reactions were run in a total volume of 30 µl, which included the following reagents: 1X PCR buffer, 0.2 mM MgCl₂, 12 U BSA, 0.25 µM of each dNTP, 0.35 µM of each primer, and 0.75 U Taq Polymerase. The thermocycler protocol implemented for the amplification was: one cycle of 95°C for 1 min 45 s, 52°C for 45 s, and 72°C for 45 s, 3 cycles of 94°C for 45 s, 52°C for 45 s, and 72°C for 45 s, 29 cycles of 92°C for 30 s, 52°C for 45 s, and 72°C for 45 s, and 1 cycle of 72°C for 7 min for the final extension. PCR and thermocycler conditions for COI amplification were similar to Cyt-b except for the annealing temperature, which was 50°C. The primers used for COI were: FISH-BCL (forward), 5′–TCAACYAATCAYAAAGATATYGGCAC–3′ and FISH-BCH (reverse), 5′–TAAACTTCAGGGTGACCAAAAAATCA–3′ (Baldwin et al., 2009). Ten microliters of the final product were run on a 1.5% agarose gel to confirm PCR amplification, and the remaining 20 µl were sent to the University of Arizona Genetics Core (DNA sequencing) for purification and Sanger sequencing.

Microsatellite amplification.—Up to 24 specimens per site were used for the analysis of the microsatellites. Samples consisting of fewer than ten individuals were not included because of the potential inaccuracies in estimates of allele frequencies. Microsatellite data were not collected for Carlana or Parastremma either given that the great genetic distance between them and Rhoadsia combined with the high rates of microsatellite evolution would likely result in substantial homoplasy. We used Schuelke's (2000) method for adding universal fluorescent tags (FAM, HEX, or NED) to the PCR products. Twelve microsatellite primer sets from Loh et al. (2014) were used based on their reported heterozygosity and ease of amplification: RA-06, RA-08, RA-09, RA-10, RA-13, RA-16, RA-20, RA-21, RA-22, RA-26, RA-27, and RA-30. PCR reactions were run in 10 µl volumes with conditions following Loh et al. (2014) except that the annealing temperatures for RA-06 and RA-08 were 54°C and 52°C, respectively, to improve PCR amplification. PCR products were sent to the University of Arizona Genetics Core (fragment analysis) for genotyping.

Sequence editing and genotyping.—Editing of the sequences of Cyt-b and COI was done using Geneious v11.0.5 ( https://www.geneious.com). The final alignment of the sequences was carried out using the MUSCLE algorithm. The alignment was inspected manually to correct possible errors made by the program. A total of 357 sequences for Cyt-b and 217 for COI including the outgroups were trimmed from the edges to obtain the same length for all sequences, 555 nucleotides for Cyt-b and 588 for COI. The sequences generated were uploaded to GenBank (MZ818343–MZ818699 and MZ820694–MZ820769).

Chromatograms for the microsatellites were scored in Geneious v11.0.5 using the plug-in for microsatellites and then manually edited to reduce errors. A total of 12 loci from 343 specimens corresponding to 16 populations of Rhoadsia were analyzed. Twenty-four samples were randomly selected and analyzed twice to confirm the accuracy of the scoring method. The percentage of missing data attributable to either amplification failure during PCR or uncertainty during the scoring process varied among loci. Loci RA-20 and RA-08 were the loci with the most missing data: 21% and 19%, respectively. Other loci had 10% or less missing data, with an average of 4% across all loci. Sites J1, J2, and J3 presented the greatest amount of missing data (12–14%). The percentage of missing data was under 8% for the other sites.

There was no evidence of linkage disequilibrium between any pair of microsatellite loci after applying Bonferroni correction for multiple test (Chi-squared, P > 0.01, α = 0.00075). The exact test for Hardy Weinberg Equilibrium (HWE) indicated that At and J2 significantly deviated from HWE for the locus RA-09 (P < 0.0001, α = 0.0003). Excluding these from the overall analyses did not fundamentally alter the results, so these were left in the analyses with microsatellites. The microsatellite data are available at:  https://github.com/aguirrelab/Rhoadsia_microsatellite_data_2021.

Linkage disequilibrium, Hardy Weinberg equilibrium, and neutrality tests.—Genepop on the Web (Raymond and Rousset, 1995; Rousset, 2008) was used to test the probability of linkage disequilibrium between each pair of loci along with the Hardy Weinberg Equilibrium exact test. The tests were corrected for multiple comparisons using the Bonferroni correction, α/n, where n is the total number of comparisons. Tajima's D (Tajima, 1989) and Fu's Fs (Fu, 1997) generated in Arlequin v3.5.2.2 (Excoffier et al., 2005) were used to test the assumption that the mutations are generated randomly in the population. Deviations could result from demographic events such as recent population expansions or selective sweeps (negative values), and recent population contractions or balancing selection (positive values). Zero is interpreted as neutral evolution.

Genetic differentiation.—Genetic divergence was calculated based on polymorphic sites among clades using DnaSP v6. Arlequin v3.5.2.2 (Excoffier et al., 2005) was used to conduct an Analysis of Molecular Variance (AMOVA). Three different AMOVAs were performed (Supplemental Table 2; see Data Accessibility). The first AMOVA evaluated the divergence among drainages, divergence among populations within drainages, and divergence within populations. The second AMOVA used the results from the phylogenetic analysis carried out with the mtDNA markers to test the divergence between groups of the main clades (northern sites vs. southern sites), among populations within these clades, and within populations. Finally, the hierarchical levels of the third AMOVA were based on the number of clustered populations estimated in the STRUCTURE analysis (i.e., K = 10). Pairwise F_ST was calculated in Arlequin to examine the level of genetic divergence between each pair of populations. F_ST values were interpreted as: little genetic differentiation ≤ 0.05, moderate genetic differentiation = 0.05–0.15, great genetic differentiation = 0.16–0.24, and very great genetic differentiation ≥ 0.25 (Hartl and Clark, 1997).

Mantel test and genetic diversity.—A Mantel test was used to test whether there was a significant correlation between genetic and geographic distances. Pairwise F_ST values generated in Arlequin were used for the genetic distance matrix. The pairwise geographic distance matrix was based on straight line distances between sites and was created using the program Geographic Distance Matrix Generator v1.2.3 (Ersts, 2021). The Mantel test was performed in Arlequin v3.5.2.2 (Excoffier et al., 2005).

Common genetic diversity indices (Table 2) based on microsatellites were calculated using GenAlEx v6.5 (Peakall and Smouse, 2006, 2012), while for the mitochondrial genes, the effective number of alleles (N_e) was calculated following Frankham et al. (2002) and haplotype diversity (HD) was calculated following Nei (1987).

Table 2

Genetic diversity indices of populations of Rhoadsia in western Ecuador based on mitochondrial genes cytochrome b (Cyt-b), cytochrome oxidase I (COI), and nuclear microsatellites (Mstl). N (number of individuals); mN (mean number of individuals analyzed); H (number of haplotypes); Na (mean number of alleles); Ne (effective number of alleles); mNe (mean of effective number of alleles across loci); HD (haplotype diversity); Ho (observed heterozygosity); He (expected heterozygosity); %PA (percentage of private alleles).

Phylogenetic analysis and haplotype networks.—Phylogenetic relationships of populations of Rhoadsia were tested using maximum likelihood (ML) and Bayesian analyses of the mitochondrial genes. A total of 61 and 64 sequences were used for Cyt-b and COI, respectively, after removal of duplicates, and were analyzed separately. Phylogenetic analysis of the concatenated mitochondrial data gave a similar result to the analyses conducted on the two genes separately but is not included because the reduced sampling of the COI gene resulted in loss of some of the Cyt-b data. The nucleotide substitution model selected for both the ML and Bayesian methods and both genes was the HKY + G model (Hasegawa et al., 1985) as suggested by the program jModelTest (Posada, 2008) using the Bayesian Information Criterion (BIC; Bhat and Kumar, 2010). PhyML (Guindon and Gascuel, 2003) implemented in SeaView v4.5.4 (Gouy et al., 2010) was used for the estimation of the ML tree using the approximate likelihood ratio test, Shimodaira-Hasegawa-like (aLRT SH-like; Guindon et al., 2010) for branch support. Default values were used for the rest of the parameters. Bayesian analysis parameters were set using the program Beauti v1.8.4 included in the package of BEAST v1.8.4 (Drummond et al., 2012). The parameters used were: a strict molecular clock which assumes uniform rates across branches, random starting tree and a Markov Chain Monte Carlo (MCMC) length of 50 million iterations with resampling every 5000 samples, discarding 5 million iterations as burn-in. The coalescent model used was Bayesian Skyride (Minin et al., 2008). The effective sample size (ESS), which measures the independence of sampled trees, was inspected using the program Tracer v1.6 (Rambaut et al., 2014), then summarized in TreeAnnotator v1.8.4 (included in Beast Package). The maximum clade credibility (MCC) tree produced was then visualized in FigTree v1.4.2 (Rambaut, 2016). Haplotype networks were inferred using the median-joining method for the Cyt-b and COI genes in the program PopART (Leigh and Bryant, 2015).

Population structure analysis.—STRUCTURE software v2.3.4 (Pritchard et al., 2000) was used to investigate the population structure based on the microsatellite markers. This method estimates the most likely number of genetically distinct clusters—known as K—where each individual is assigned to a distinctive cluster based on a Bayesian algorithm. The STRUCTURE analysis provides valuable insight into the number of genetically distinct units of Rhoadsia that should be recognized in western Ecuador. Twenty replicates of each possible number of K (1 to 16, the latter being the number of sites sampled) were run with 1,000,000 MCMC iterations after 250,000 of burn-in. An admixture model was used, which assumes that the populations might have multiple origins and have some genetic admixture. The allele frequencies correlated model was selected. The most likely value of K was estimated according to the mean Log-normal probability of K (LnPK) provided directly by STRUCTURE (Pritchard et al., 2000) and also by calculating ΔK (Evanno et al., 2005) using the program Structure Harvester (Earl and VonHoldt, 2012). Cluster Markov Packager Across K (CLUMPAK) was used to analyze the 20 independent runs selected after calculating the “best” value of K. It generates a consensus solution of the clustering using a Markov Clustering algorithm which relies on the similarity of the matrix between replicate runs (Kopelman et al., 2015).

RESULTS

Genetic divergence between Rhoadsia and the outgroups.—Rhoadsia differed substantially from the outgroups Carlana eigenmanni (Central America) and Parastremma sadina (Colombia), forming a clade relative to them for both mitochondrial genes. Sequences of C. eigenmanni formed two well-supported clades, one represented by samples from Costa Rica and Nicaragua and the other by samples from Panama. The four sequences of P. sadina were identical to each other. For the Cyt-b gene, both outgroups combined showed a divergence of 28% (153 mutations) when compared to populations of Rhoadsia. Samples from C. eigenmanni from Central America alone diverged about 23% (128 mutations) from Rhoadsia, while P. sadina from Colombia diverged about 20% (115 mutations). For the COI gene, the total number of mutations for the outgroups was 134 (23% divergence). Carlana eigenmanni diverged about 17% (101 mutations), while P. sadina diverged about 16% (96 mutations) from Rhoadsia.

Genetic variation in Rhoadsia.—There was substantial genetic variation documented in Rhoadsia for both the mitochondrial and nuclear markers, especially in the larger drainages (Table 2, Figs. 3–6). There were 34 polymorphic sites (6.13%) with 36 mutations total in the Cyt-b sequence alignment for specimens of Rhoadsia without the outgroups. Thirty-one of these mutations were synonymous and five were non-synonymous. The non-synonymous mutations were distributed in one individual from G1 (1 mutation), one individual from E2B (1), three individuals from SP (1), and in all individuals sequenced from Su and At (2 mutations). The COI sequence alignment showed 20 polymorphic sites (3.4%) and 21 mutations total. All mutations were synonymous. The neutrality test showed no statistically significant deviation from zero for Tajima's D and Fu's Fs in both mitochondrial genes, suggesting that these genes are evolving neutrally in Rhoadsia.

Fig. 3

Frequency of cytochrome b haplotypes (Cyt-b, left) and cytochrome oxidase I (COI, right) per site. Haplotypes colored in white correspond to private haplotypes (PH; i.e., haplotypes that are unique to a particular site). The bottom left corner graphs (a and b) present the haplotype frequency showing the number of individuals sampled in each drainage. Su-At (Sua and Atacames combined), Esm (Esmeraldas), Gua-Ta (Guayas and Taura combined), SP (San Pablo), Jub (Jubones), and SR (Santa Rosa).

Fig. 4

Haplotype network for the mitochondrial genes cytochrome b (Cyt-b, top) and cytochrome oxidase I (COI, bottom) color coded by site. The size of the circles is proportional to the haplotype frequency. The number of mutations between the haplotypes are represented by hatch marks. The populations within drainages are represented by different color shades (see legend). Haplotypes fall into two groups: the northern (N, enclosed by blue dashed line) and the southern group (S, enclosed by red dashed line).

Fig. 5

Bayesian phylogenetic trees based on the mitochondrial genes cytochrome b (Cyt-b, left) and cytochrome oxidase I (COI, right). Nodes supported by posterior probabilities greater than 0.7 are displayed with black dots unless specified (arrow). Clade I (blue) comprised populations from the northern region, while Clade II (red) is represented by populations from the southern region. The format for the label of the ingroup is: Site code_Drainage code_Haplotype number_Year_(number of individuals with identical haplotype). For the outgroup, it is: Species code_Country code_Sample ID_(Number of individuals with identical haplotype).

Fig. 6

Structure analysis of populations of Rhoadsia in Ecuador based on 12 microsatellites. Neighbor-joining of clusters (K = 10) inferred in STRUCTURE (top left). Proportion of membership of each cluster per site (top right). Proportion of membership of each cluster assigned to each individual (bottom). Clusters fell into three main groups (I, II, and III) based on the distance matrix generated in STRUCTURE. We tentatively recommend recognition of the ten clusters from the STRUCTURE analysis as evolutionarily significant units (ESUs).

Populations from Esmeraldas, Guayas, San Pablo, and Santa Rosa tended to have greater genetic diversity than populations from the Sua, Atacames, and Jubones drainages for both Cyt-b and microsatellite data (Table 2). For COI, populations from the Jubones drainage had the greatest genetic diversity. All the genetic markers examined concur that the small drainages in northwestern Ecuador, the Atacames and Sua, which are geographically close to the Esmeraldas drainage, had the lowest genetic diversity. There was a tendency for a decline in genetic diversity with elevation within drainages but it was not universal. Private haplotypes for the mitochondrial genes were relatively common but were more frequent for Cyt-b than for COI and their frequency varied among drainages (Table 2, Fig. 3).

Phylogenetic analysis of Rhoadsia mitochondrial haplotypes.—Haplotype networks for Cyt-b and COI were similar in general topology and clustered haplotypes into northern and southern groups roughly consistent with the distributions of R. minor and R. altipinna (Fig. 4). Haplotypes from the northern region were separated from the southern haplotypes by three mutations for Cyt-b and one mutation for COI, with an average sequence divergence of 6% and 3.4% of polymorphic sites for Cyt-b and COI, respectively. For Cyt-b, H9, one of the most common haplotypes (freq. 0.29), was the only haplotype shared between sites from the northern and southern regions (Supplemental Table 3; see Data Accessibility). Atacames and Sua, small isolated drainages in northwestern Ecuador, had private haplotypes (H12 and H22). Haplotypes from the northern Guayas (i.e., G1, G2, G6, and G7), such as H15, H18, and H20, had more similarities with haplotype H9 (most common haplotype in the northern group) than with haplotypes from the southern sites. The southern populations had 18 haplotypes, most of them connected to H5. Haplotypes from the southern Guayas (i.e., G3, G4, and G5) were mostly private to this drainage and were present at lower frequencies. For the COI gene, H3 was the most common haplotype present in the northern region (Supplemental Table 4; see Data Accessibility). Haplotype H1 was also shared by one site in the Esmeraldas (E8) and the northernmost sites in the Guayas (i.e., G1, G2, G6, and G7). H4, which was present almost exclusively in southern populations (and one individual from E6), was the closest southern haplotype to H3, separated by one mutation. H6 was another southern haplotype that was present in one individual from Esmeraldas (E8) along with J1, J4, and SR.

Table 3

AMOVA of populations of Rhoadsia in western Ecuador grouped by drainage, clades inferred from mtDNA genes cytochrome b (Cyt-b) and cytochrome oxidase I (COI), and clusters estimated in STRUCTURE analysis (K = 10) based on the nuclear microsatellites (Mstl). See Supplemental Table 2 (see Data Accessibility) for the specific assignment of the sites into the groups of each version of AMOVA per molecular marker. df (degrees of freedom). P-value was very significant (<0.0001) for all analyses unless otherwise specified. * P < 0.01.

Bayesian phylogenetic analysis of the mitochondrial genes identified two main clades within Rhoadsia, clustering sites from northern and southern Ecuador (Fig. 5). The maximum likelihood analysis also showed a well-defined northern clade; however, southern populations did not form a discrete clade (not shown). Topologies inferred using Cyt-b and COI data yielded similar results and the northern–southern segregations of haplotypes were generally consistent with the occurrence of two closely related species of Rhoadsia, R. minor in northwestern Ecuador and R. altipinna in southwestern Ecuador. However, the separation of haplotypes did not correspond with the currently recognized geographic ranges of R. minor and R. altipinna. Specimens from the northern Guayas sites G1, G2, G6, and G7, the same drainage from which R. altipinna was described, mostly consisted of haplotypes from the northern Clade I, which presumably is R. minor. This was true for both mitochondrial genes. There was also phylogenetic structure within the major clades that corresponded to patterns of geographic isolation, like the grouping of specimens from the small At and Su drainages as a subclade within clade I. Some patterns of phylogenetic structure were surprising, like the separation of haplotypes from the Jubones drainage, in which haplotypes from the high elevation sites J4 and J5 grouped with haplotypes from the neighboring San Pablo drainage, while haplotypes from the low elevation sites J1–J3 grouped with haplotypes from the SR and Guayas sites.

Population structure analysis with the microsatellites.—Population structure analysis of the 12 microsatellites suggested the presence of ten genetic clusters of Rhoadsia in western Ecuador (Fig. 6). Subdivision of the microsatellite dataset into ten clusters or populations (i.e., K = 10) was supported by both methods used. The Evanno method showed the greatest ΔK (ΔK = 66.27) at K 10 (Supplemental Figure 1; see Data Accessibility), and the LnPK value calculated using the Pritchard method formed a plateau at K 10 (Supplemental Figure 2; see Data Accessibility). The Structure analysis suggested some degree of population admixture in some sites. The most conspicuous example was observed between the northernmost site in the Guayas drainage, G1, and sites from the neighboring Esmeraldas drainage and southern Guayas sites. Although clustering generally correlated with the distance among sites, Cluster 8 (SP, J4, and J5) linked populations from different drainages with no modern fluvial connection. In contrast, site At (Cluster 1) was relatively highly divergent given its close geographic proximity to the neighboring sites E2B and E4.

A neighbor-joining tree based on the ten clusters estimated by STRUCTURE was created to group them based on similarity (Fig. 6). The ten clusters formed three main groups. The first group (I—the northern group) was represented at all sites from the Esmeraldas and Atacames drainages (clusters 1–4) and encompassed the type locality of R. minor. The second group (II—central-southern group) contained all sites from the Guayas drainage plus SP and the two highest sites in the Jubones drainage (J4 and J5, clusters 5–8). The geographic distribution of this group encompasses the type locality of R. altipinna. Finally, the third group (III—the southernmost group) includes SR (cluster 9) and J1, J2, and J3 from low elevation sites in the Jubones drainage (cluster 10), and also corresponds geographically to populations traditionally recognized as R. altipinna (Jiménez-Prado et al., 2015).

Partitioning of genetic variation in Rhoadsia.—Populations of Rhoadsia were grouped in three different ways for the AMOVA: by drainage, according to the results of the phylogenetic analysis of the mtDNA genes, and by clusters estimated in the STRUCTURE analysis of microsatellite data (Supplemental Table 2; see Data Accessibility). Among-groups comparisons were an important source of variation in all analyses of the mitochondrial genes, but not of the microsatellites. The AMOVA with the highest percentage of variation attributed to the among-groups component was from the clusters estimated in the STRUCTURE analysis. For the Cyt-b and COI mitochondrial genes, the among-groups component was 75% and 53%, respectively, while the within-group component corresponded to only 29.8–37.4% for the COI data and 14.4–19.3% for the Cyt-b data. In contrast, the among-groups component based on the microsatellite data was relatively low in all AMOVAs (5.0–16.5%), and the within-population component showed the highest percentage of variation, ≈80% in all three AMOVAs (Table 3).

Pairwise F_ST values and Mantel tests indicated that genetic divergence was related—in most cases—to geographical distances between and within drainages (Supplemental Figure 3; see Data Accessibility). F_ST values for COI were mostly low between sites within drainages (Supplemental Table 5; see Data Accessibility), whereas cytochrome b F_ST values showed considerably higher genetic divergence among populations between and within drainages (Supplemental Table 6; see Data Accessibility). Pairwise F_ST values for microsatellites were less correlated with geographic distances (Supplemental Figure 3, Supplemental Table 7; see Data Accessibility). The F_ST data supported several of the patterns of divergence documented previously in other analyses. The Atacames site was quite divergent from all other sites based on the Cyt-b and microsatellite data. Both the mitochondrial and nuclear markers indicated an apparent north–south divergence within the Guayas drainage such that two groups formed: a northern group consisting of G1, G2, G6, and G7 (which were similar to the geographically adjacent E8), and a southern group consisting of G5, G3, and G4. According the Cyt-b and microsatellite data, G1 from the northern Guayas drainage exhibited genetic similarities to other drainages including E8 and E2B from the Esmeraldas drainage. Also, G5 and G3 from the southern Guayas group were closely related to individuals from the neighboring site SP and high elevation Jubones sites J4 and J5. The Santa Rosa, the southernmost drainage sampled in this study, was genetically highly differentiated from the rest of the sites except its closest neighbors J1 and J2. Finally, the low elevation Jubones sites J1, J2, and J3 were genetically divergent from the high elevation sites J4 and J5, forming two genetically distinct groups segregated by elevation within the same drainage.

DISCUSSION

We conducted a phylogeographic analysis of Rhoadsia in western Ecuador and addressed two primary issues. First, we examined the validity of the two currently recognized species in this genus. The phylogeographic analyses based on both mtDNA genes (Cyt-b and COI) and the STRUCTURE analysis with microsatellites confirm the existence of two evolutionary lineages respectively identifiable as R. minor and R. altipinna, with introgression of alleles from R. minor into R. altipinna in the northern Guayas drainage, and significant geographic population structuring in R. altipinna (Figs. 5, 6). The microsatellite data gave a finer structuring pattern to populations of Rhoadsia with the identification of ten genetic clusters. The second issue was how well patterns of genetic divergence in populations of Rhoadsia correspond with previous hypotheses of aquatic zoogeographic regions (Barriga, 2012). The microsatellite data showed that the ten genetic clusters estimated for Rhoadsia are divided into three main groups with two major gaps separating clusters of populations, as well as a few minor gaps. One major gap was located between the Esmeraldas and the Guayas drainages, while the second is in the southern region between the high- and low-elevation Jubones sites. The first gap (Esmeraldas–Guayas) seems to be permeable given the significant introgression of R. minor from the Esmeraldas drainage into the northern Guayas drainage based on the discordance between the mitochondrial and microsatellite data. The second separates clusters of populations at a point within the same drainage (Jubones), a region that is geographically close to the gap that Barriga (2012) hypothesized for the boundary between the Guayas and Catamayo ichthyohydrographic zones.

Species delimitation and evolutionary relationships of populations of Rhoadsia in Ecuador.—Delimiting species is a difficult task considering the large number of species concepts that exist, and their limitations (Aldhebiani, 2018). This study employs a genetic approach in which different types of markers (i.e., mtDNA and microsatellites) were used to identify evidence of independently evolving metapopulations (i.e., multiple populations of a species that coalesce into a common ancestor; De Queiroz, 2007; Reeves and Richards, 2011). This approach has been quite useful when dealing with genetically divergent populations and cryptic species (Unmack et al., 2012; Arteaga et al., 2016; Bagley et al., 2016; Fennessy et al., 2016; Schield et al., 2018).

Populations of Rhoadsia in Ecuador were split into two main lineages allopatrically distributed in northern and southern regions of western Ecuador based on the analysis of mtDNA (Fig. 5). The two phylogroups encompassed the respective type localities of Rhoadsia minor (north) and R. altipinna (south), with some admixture near the limits of the species ranges. The distinctive clade formation was more apparent with Cyt-b than with COI, the latter presenting considerably lower node support particularly in the southern clade, which could be attributed to its slower rate of evolution (Kress et al., 2015). The STRUCTURE analysis of the microsatellite data provided a finer subdivision of Rhoadsia into ten population clusters organized into three broad groups allopatrically distributed from north to south.

Given the pattern of genetic divergence, we suggest that R. minor and R. altipinna should be retained as valid, allopatrically distributed sister species, although they are very closely related genetically and it is not clear how they should be diagnosed morphologically. Body shape, which has historically been used to distinguish between them, is now known to change with elevation in parallel in both species, such that low elevation R. minor in the Esmeraldas drainage overlap in body shape with high elevation R. altipinna from the Guayas and other southern drainages (Jiménez-Prado et al., 2015; Malato et al., 2017). Rhoadsia minor is still distinctively more streamlined at high elevations sites near its type locality in the Esmeraldas drainage (Jiménez-Prado et al., 2015; Malato et al., 2017); however, these populations are genetically closely related to the low elevation populations in the Esmeraldas and do not form a distinct clade. The two recognized species appear to show evidence of introgressive hybridization, according to the mitochondrial genes, with haplotypes common in R. minor also present in northern populations of R. altipinna (Fig. 5). The microsatellite data split populations in a matter that was more consistent with the traditionally recognized species boundaries, although there was also evidence of admixture in the northern Guayas populations, suggesting historical or ongoing gene flow between the species (Fig. 6).

We acknowledge that single genes like Cyt-b and COI may not precisely represent the species tree of Rhoadsia and should not be used as sole evidence of species delimitation (White et al., 2014). The genetic divergence between species of Rhoadsia is also relatively small. However, Rhoadsia altipinna and R. minor were described in the early twentieth century and have been recognized as distinct species since then. Rhoadsia minor also exhibits relatively strong divergence in body shape at high elevations of the Esmeraldas basin that has not been observed in R. altipinna, suggesting the potential for some adaptively important differences between species (Malato et al., 2017). The combination of the mtDNA and microsatellite data showing genetic divergence that roughly coincides with the geographic boundaries for the two species suggests that R. altipinna and R. minor should be retained as two closely related but distinct species for now.

Identifying evolutionarily significant units (ESU) of Rhoadsia in western Ecuador.—The term ESU is applied to populations that are genetically distinctive enough to be considered separate populations with distinct evolutionary histories (Ryder, 1986; Moritz, 1994; Palsbøll et al., 2007). When the isolated population seems rare and/or endangered, the evolutionarily significant unit could also be treated as a management unit in which conservation efforts can be enforced (Palsbøll et al., 2007). Rhoadsia appears to encompass multiple genetically isolated populations that may be considered as ESUs within each named species given the interesting patterns of morphological divergence that have been documented among populations (Malato et al., 2017; Aguirre et al., 2019) and the threat that freshwater fishes in western Ecuador increasingly face (Dodson and Gentry, 1991; Jiménez-Prado et al., 2015; Aguirre et al., 2021). Below we make recommendations for recognizing potential ESUs of Rhoadsia, pending confirmation through future studies of ecological, morphological, and adaptive divergence.

Samples of Rhoadsia minor in northern western Ecuador fell into group I (Fig. 6). Sua (Su) and Atacames (At) populations, which inhabit small geographically adjacent drainages, formed one of the most distinctive genetic clusters within Group I. Their genetic similarity may be due to the existence of a passage that allows for gene flow between the populations (perhaps during sporadic flooding), recent separation of the drainages, or recent headwater capture between the drainages. These populations also have two unique, non-synonymous mutations in Cyt-b that were not found in any other populations of Rhoadsia. The Su-At subclade also had the lowest levels of genetic diversity, perhaps due to the small size of the drainage basins they inhabit or historical bottlenecks (Table 2, Fig. 3). Given the relatively high level of genetic differentiation and low genetic diversity of these populations, the small size of the drainage basins that they inhabit, and the degradation that these basins are suffering (Jiménez Prado et al., 2020; Aguirre et al., 2021; Jiménez-Prado and Aguirre, 2021), we suggest designating this genetic cluster as an ESU that should be given conservation priority.

The remaining samples from Group I were from the Esmeraldas drainage basin and fell into three groups based on the microsatellite data, although these are similar enough to suggest sustained gene flow or recent divergence among populations within the drainage. The greatest divergence occurred among populations located at high elevations, with the two highest sites, Mindo (E6; Cluster 3) and Transito (E8; Cluster 4), located at 1260 m above sea level (m.a.s.l.) and 1093 m.a.s.l., respectively, being significantly divergent according to the Structure analysis conducted on the microsatellites (Fig. 6). In contrast, Teaone (E2B) and Hermoso (E4), located at lower elevations (64–282 m.a.s.l.), were similar and fell within Cluster 2. Given the identification of three distinct groups within the Esmeraldas basin and extrapolating to neighboring sites, we suggest designating three distinct ESUs, one for the low elevation populations in the Esmeraldas basin(< 300 m), one for the high elevation populations in the Toachi River (> 600 m.a.s.l.; i.e., sites E7 and E8), and one for the high elevation populations in the Blanco River (> 600 m.a.s.l.; i.e., sites E5 and E6). The ESU in the northern Blanco River sub-basin deserves particular attention because it includes the locality of Mindo, which harbors morphologically distinctive highly streamlined populations and is the type locality of R. minor (Malato et al., 2017).

Individuals of Rhoadsia altipinna from the northernmost part of Guayas, Chiguilpe (G6), Otongo (G7), Palenque (G1), and Jauneche (G2), displayed haplotypes characteristic of both R. minor and R. altipinna (Fig. 5). It is possible that both species diverged in isolation forming northern and southern species and there has been introgressive hybridization in the northern reaches of the Guayas drainage. The Esmeraldas and Guayas River basins come in very close geographical proximity to one another around the city of Santo Domingo, and a similar pattern of intermediate phenotypes between otherwise allopatrically distributed northern and southern sister species has been observed in this part of the Guayas in the suckermouth armored catfish genus Transancistrus (Lujan et al., 2015). Sites from the northern Guayas seem to genetically differ from sites from the southern Guayas (i.e., Caluma [G5], Chague Grande [G3], Chimbo [G4], and Taura [Ta]) according to all genetic markers, especially COI. Although the factors causing such genetic breaks are unclear, we suggest dividing the Guayas basin into three distinct ESUs, one corresponding to the northernmost Guayas populations (G1, G6, G7), one for G2 which fell in its own cluster in the STRUCTURE analysis, and one for the southern Guayas populations (G3, G4, G5) associated with the Babahoyo River sub-basin.

South of the Taura drainage, there are various small basins that drain east to west into the Pacific Ocean. Despite their current isolation, historical connections may have allowed for migration between these neighboring basins given their similar species communities (Barriga, 2012). A population from the San Pablo (SP) basin separated by about four basins from the Guayas, showed significant genetic divergence from populations in the Guayas basin. The SP population was closely related to sites J4 and J5 in the upper Jubones drainage basin according to the Cyt-b and microsatellites (Supplemental Tables 6, 7; see Data Accessibility). Interestingly, J4 and J5 are located about six basins south at 909 m and 1095 m of elevation, respectively, in the Jubones drainage, while the SP population was sampled at 50 m of elevation. We suggest that this cluster formed by J4, J5, and SP be designated an ESU.

The Jubones drainage, located in southern Ecuador, is segregated into two distinct genetic clusters based on elevation. As indicated above, the two highest sites (J4 and J5) formed a cluster with the site SP (Cluster 8). The three lowest sites sampled (J1, J2, and J3) formed their own cluster (Cluster 9) according to the STRUCTURE analysis of the microsatellite data and grouped most closely with the Santa Rosa River sample. This pattern was also observed with the Cyt-b data but was less obvious with the COI data, for which some haplotypes seem to be shared between sites J3, J4, and J5 (Figs. 3, 5). The level of differentiation between the low and high elevation samples from the Jubones River is surprising (Supplemental Tables 5–7; see Data Accessibility) given that they inhabit the same relatively small river basin. It is worth noting that there was a relatively large gap in elevation between the sites J3 (251 m.a.s.l.) and J4 (909 m.a.s.l.). The sampling area between these sites was difficult to access and when it was possible, there were no Rhoadsia present. It is possible that there are natural barriers like the presence of waterfalls or steep slopes with rapid water that are isolating fish populations (Silva et al., 2016). It is also worth noting that there is evidence that the lower Jubones River historically changed course. Maps of the Jubones River from the nineteenth century had it draining south of the city of Machala into the Estero de Jambelí, while today it drains well north of Machala (Wolf, 1892). If the Jubones River indeed changed course, the previous drainage route would have put low elevation populations of the Jubones River closer geographically to those of the Santa Rosa River, with which they share some genetic similarity. Given their genetic distinctiveness, we suggest recognizing the low elevation populations of the Jubones River (J1, J2, and J3) as an ESU and designating the SR population as another distinct ESU given that it forms its own group in the STRUCTURE analysis.

Aquatic zoogeographic zones in western Ecuador.—Fluctuations in topography and environmental conditions have contributed to the isolation and diversification of many species in western Ecuador (Albert and Reis, 2011; Rodríguez Tribaldos et al., 2017). Unfortunately, there are few population genetic studies on aquatic species in western Ecuador (Aguirre et al., 2013; Cucalón and Bajaña, 2015).

Our phylogeographic analysis of Rhoadsia provides an opportunity to examine how well previous hypothesis for aquatic biogeographic zones in western Ecuador hold up. Barriga (2012) divided western Ecuador into five ichthyohydrographic zones. The first is an intertidal zone representing coastal freshwater streams subject to tidal influence throughout western Ecuador with a large proportion of gobies, eleotrids, and other estuarine tolerant species. The other four zones from north to south are the Santiago-Cayapas, Esmeraldas, Guayas, and Catamayo zones, with the first three roughly corresponding to their respective drainage basins and the Catamayo zone starting at the Jubones River and extending south. Our samples were collected in the latter three zones. Although there seem to be historical records for Rhoadsia in the Santiago-Cayapas basin (Barriga, 2012), it has not been collected in recent years (Pedro Jiménez-Prado, pers. comm.) and its status there is uncertain. Nonetheless, the turnover of freshwater fish assemblages between the Santiago-Cayapas and Esmeraldas drainages indicate that the Santiago-Cayapas is a distinct biogeographic zone for freshwater fishes in western Ecuador with greater similarity to drainages of the Colombian Chocó (Jiménez-Prado et al., 2015; Navarrete et al., 2021).

The distribution of genetic variation in Rhoadsia follows Barriga's (2012) ichthyohydrographic zones relatively well, with some caveats. We indeed find an important break between the Esmeraldas and Guayas drainage basins although introgression of R. minor haplotypes into the northern Guayas basin indicates opportunities for gene flow between basins, which may result in a more complex scenario for the break that varies among species. It is worth noting that the Guayas and the Esmeraldas drainages almost meet near the city of Santo Domingo. The Toachi River in this locality has been reported to be the place where both Esmeraldas and Guayas drainages start northward and southward, respectively (Ministerio de Ambiente del Ecuador, 2012). Moreover, local news agencies in Ecuador have reported flooding by the Toachi River in the area that may periodically connect the drainages (Guerrero, 2016). Opportunities for movement between drainages is consistent with the relatively large number of species shared between these basins and with phenotypic patterns observed in a recent taxonomic revision of the armored catfish genus Transancistrus (Jiménez-Prado et al., 2015; Lujan et al., 2015). Our data on Rhoadsia also indicates that there is likely significant population genetic structuring within the Guayas River basin, possibly providing the opportunity for local adaptation in different parts of the basin. Given its size, habitat heterogeneity, and the distances between some tributaries, we suspect that genetic divergence among populations within the Guayas drainage may not be unusual, especially in less migratory species. However, a similar phylogeographic study on Hoplias microlepis, a larger, low-elevation predatory fish, showed that most populations sampled within the Guayas drainage shared very similar genetic makeup with no major divergence pattern (Aguirre et al., 2013; Cucalón and Bajaña, 2015). Thus, this likely varies among fishes with different attributes.

South of the Guayas, we found a third genetic cluster based on the microsatellite data, indicating that populations in the Jubones and the Santa Rosa Rivers differ significantly from those in the Guayas River. This coincides with Barriga's (2012) Catamayo zone, which splits from the Guayas zone at the Jubones River. We actually detected significant divergence from the Guayas samples in the SP sample, which is significantly north of the Jubones, indicating that the break between the Guayas and the Catamayo zone may occur farther north than the Jubones River in some species. Regardless of the precise location of the gap, our data for Rhoadsia are consistent with the existence of a distinct aquatic biogeographic zone in southwestern Ecuador close to Peru, with the potential for harboring locally adapted populations or possibly even endemic species. Although our study did not include samples from Peru, we speculate that the northwestern Peruvian rivers Chira and Tumbes are also part of this aquatic biogeographic zone.

Conclusion.—This study contributes to knowledge of the evolutionary relationships of populations of freshwater fishes in western Ecuador, a region known for having many endemic species that are facing grave environmental degradation caused by human activity (Aguirre et al., 2021), and an overall scarcity of phylogeographic studies of aquatic species (Damanik-Ambarita et al., 2016; Ribeiro et al., 2017). Our results provide a baseline for future research on the biogeography, evolution, and conservation of the aquatic biodiversity of western Ecuador, as well as a phylogenetic framework to examine the genetic basis of adaptation of morphologically divergent populations of Rhoadsia inhabiting Andean mountain streams in the region.

DATA ACCESSIBILITY

Supplemental material is available at  https://www.ichthyologyandherpetology.org/i2020092. The mtDNA sequences generated were uploaded to GenBank (MZ818343–MZ818699 and MZ820694–MZ820769). The microsatellite data are available at:  https://github.com/aguirrelab/Rhoadsia_microsatellite_data_2021. Unless an alternative copyright or statement noting that a figure is reprinted from a previous source is noted in a figure caption, the published images and illustrations in this article are licensed by the American Society of Ichthyologists and Herpetologists for use if the use includes a citation to the original source (American Society of Ichthyologists and Herpetologists, the DOI of the Ichthyology & Herpetology article, and any individual image credits listed in the figure caption) in accordance with the Creative Commons Attribution CC BY License.