Ancient genetic divergence in bumblebee catfish of the genus Pseudopimelodus (Pseudopimelodidae: Siluriformes) from northwestern South America

Study area and sampling

This study analyzed a total of 257 Pseudopimelodidae muscle tissues from specimens collected in seven hydrographic sub-zones of Colombia (IDEAM, 2015; Table 1, Fig. 1): (1) Magdalena River upper sector, (2) Magdalena River middle sector, (3) Magdalena River and Cauca River- lower sectors and San Jorge River, (4) Cauca River upper and middle sectors, (5) Caribbean drainage (Atrato River and Sinú River), (6) Amazon River hydrographic zone (Orteguaza River and Vaupés River), and (7) Orinoco River hydrographic zone (Negro River and Meta River).

Muscle tissues of specimens from Magdalena River, Cauca River and Atrato River, preserved in 97% ethanol, were provided by Integral S.A., through the scientific cooperation agreement CT-2013-002443. In addition, muscle tissue from Pseudopimelodus specimens from Sinú River and Orinoco-Amazon hydrographic zones, as well as members of the genera Cruciglanis, Microglanis and Rhyacoglanis, were collected, identified, and provided through the scientific agreements 00466 FUNINDES-INCODER; 0003, 0033, and 00187 FUNINDES-AUNAP; and 037-2014 FUNINDES-HUMEDALES.

Some collection sites were located in sectors with steep slopes in the Magdalena River- upper sector, upstream of the municipality of Honda, which marks the boundary between the upper and the middle Magdalena. Collection sites in the Cauca River upper and middle sectors were located upstream in the Cauca River canyon, which is the steepest margin of the Antioqueño Plateau in the northern portion of the Central Cordillera (Restrepo-Moreno et al., 2009) and marks the boundary between the middle and lower sectors of Cauca River. This landform has been considered a geographic barrier for many fish species (Dahl, 1971) and is the focus of the largest hydropower project in Colombia (Hidroituango). Additionally, sampling localities in Magdalena River and Cauca River lower sectors and the San Jorge River correspond to lowlands sites downstream in the Cauca River canyon, including flood plain “ciénagas” that are part of the Momposina depression, one of the most important wetlands in northwestern South America. Finally, collection sites in the cis-Andean hydrographic sub-zones of Orinoco (Meta River and Negro River) and Amazon (Orteguaza River and Vaupés River) are separated by the Vaupes Arch, the major drainage that divide the Llanos region of eastern Colombia (Winemiller & Willis, 2011).

DNA extraction, PCR amplification and sequencing

Genomic DNA was extracted using the PureLink^® Genomic DNA extraction kit (Invitrogen), following the manufacturer’s protocol. For phylogeographic analyses, we amplified a partial region of cox1 gene in all 257 specimens, using the universal primer cocktail VF2, FishF2, FishR2 (Ward et al., 2005), and FR1d (Ivanova et al., 2007), as previously published. Furthermore, rag2 gene was amplified from 77 samples, using the primers MHRAG2-F1 and MHRAG2-R1, as reported by Hardman & Page (2003). For both markers, polymerase chain reactions (PCR) were performed with a total volume of 30 µl, containing 2.5 µl of genomic DNA, 3 µl 10X buffer, 1.5 µl MgCl₂ (50 Mm), 0.75 µl dNTPs (10 mM), 0.75 µl of the primer cocktail (0.25 pmol/µl each), 0.15 µl of Platinum™ Taq DNA Polymerase (5 U/µl), and 20.6 µl of sterile nuclease-free water (Amresco). The thermal profile for cox1 consisted of an initial cycle at 95 °C for 3 min, 32 cycles of 94 °C for 3 min, 25 s at 60 °C, and 30 s at 72 °C, and the final extension was omitted. A similar thermal profile was used to amplify rag2, but with changes in annealing temperature (58 °C) and extension time (45 s). PCR products were confirmed by agarose gel electrophoresis and EZ vision staining and sequenced in both forward and reverse directions using an ABI sequencer 3730XL. Sequences were checked and edited using Geneious v10.0.9 (https://www.geneious.com) software, and multiple alignment to create consensus sequences was performed with the MAFTT v7.308 (Katoh & Standley, 2013). Both genes were translated into amino acids to confirm absence of stop codons or unexpected frameshift errors using the sequence translation tool Transeq, available in The European Molecular Biology Open Software Suite (EMBOSS https://www.ebi.ac.uk/Tools/st/emboss_transeq/).

Genetic diversity and phylogenetic relationships

Calculation of genetic diversity indices, such as number of haplotypes, haplotype (Hd) and nucleotide (π) diversity (Nei, 1987), and polymorphic sites, was performed using DNAsp v6 (Rozas et al., 2017). Phylogenetic analysis for each gene and concatenated dataset was conducted using Bayesian Inference, with MrBayes (MB) v3.2 (Ronquist et al., 2012), based on the best-fit model estimated in the software IQ-TREE v1.6.12 (Kalyaanamoorthy et al., 2017). Parameters included two independent Markov Chain Monte Carlo (MCMC) iterations for 20 million generations sampled every 1000 generations, discarding the first 25% sampled generations as burn-in; settings for remaining parameters were left at their default values. Convergence of the MCMC was based on the Potential Scale Reduction Factor (PSRF) approaching 1.0 and the standard deviation of split frequencies approaching 0.0. In addition to sequences obtained in this study, the phylogenetic analysis included previously published rag2 and cox1 sequences of species in all other genera of the family Pseudopimelodidae and the pimelodid species Pseudoplatystoma magdaleniatum and Pimelodus yuma as outgroups (Table 1). In addition, spatial relationships of cox1 haplotypes were explored by constructing a haplotype network based on the Median-Joining algorithm (Bandelt, Forster & Röhl, 1999), using the software PopART v1.7 (Leigh & Bryant, 2015).

Species delimitation

Mitochondrial cox1 gene sequences were used to implement three single-locus species-discovery (SLSD) methods in R v3.4.163, by using the scripts developed by Machado et al. (2018) and available at https://github.com/legalLab/publications. These methods were implemented to define a priori clusters of putative species (lineages) of Pseudopimelodus as trait values for StarBEAST2 v14 (Ogilvie, Bouckaert & Drummond, 2017) species tree analysis. The three methods implemented corresponded to (1) a divergence threshold optimizing and clustering approach (locMin; Brown et al., 2012), (2) the general mixed Yule coalescent model (GMYC; Pons et al., 2006; Fujisawa & Barraclough, 2013), and (3) a Bayesian approximation of the GMYC (bGMYC; Reid & Carstens, 2012). For all analyses, a set of ultrametric trees was generated with BEAST v1.8.4 (Drummond & Rambaut, 2007) choosing HKY+G as the best-fit substitution model, selected by IQ-TREE v1.6.12, setting a relaxed molecular clock model and a Birth-Death Process as a tree prior (Gernhard, 2008). We ran three independent MCMC chains for 20 million sampled every 18,000 generations, discarding the first 10% of trees as burn-in. Results were combined, and then we sub-sampled a total of 1,000 trees for subsequent analyses. Tracer v1.7 (Rambaut, 2014b) was used to verify convergence and ESS values (>300).

All three analyses were based on point estimates from the maximum clade credibility tree generated using TreeAnnotator v1.8.4 and calculation of confidence intervals from the posterior sample of 1,000 trees. The R packages bGMYC v1.0.2 (Reid & Carstens, 2012), splits v1.0-19 (Fujisawa & Barraclough, 2013), and ape v4.1 (Paradis, Claude & Strimmer, 2004) were used in these analyses. A conservative posterior probability of con-specificity at 0.05 were used to summarize the bGMYC posterior samples into putative species.

Species tree and divergence time estimation

We combined sequence information from both mitochondrial cox1 and nuclear rag2 markers to implement a dated species tree considering the molecular substitution model of each gene using StarBEAST2 v14 (Ogilvie, Bouckaert & Drummond, 2017) implemented in BEAST2 v2.6.2 (Bouckaert et al., 2014). For both markers we set a strict molecular clock, constant population as population size model and Yule Process for species tree prior. To estimate the time to the most recent common ancestor (TMRCA), and the corresponding credibility intervals (95% HPD) for the main lineages, two fossil calibrations points were included: one corresponds to the oldest pimelodid fossil from the Paleogene of South America (Gayet & Otero, 1999) that provide a minimum age of 55.8 million years ago (mya), for the divergence between the families Pimelodidae and Pseudopimelodidae. This single calibration age constraint was set in the clade including the two pimelodids species used as outgroups, specifying a normal distribution prior, with a standard deviation equal to 1. Furthermore, an additional fossil calibration corresponding to a cf. Cephalosilurus or Pseudopimelodus from South America middle Miocene (lognormal distribution, mean 0.1, SD 0.8, offset 11.5, range 15.9–11.5 mya; Lundberg, 1998; Lundberg et al., 2010) was incorporated on the ancestral node containing all trans- and cis-Andean Pseudopimelodus. Analysis was run for 150,000,000 generations sampling every 15,000 generations and discarding the 10% of the burn-in samples. Convergence was assessed using Tracer v1.7. The summary of all trees and the annotation of mean ages of all nodes, and the corresponding HPD ranges, was conducted with TreeAnnotator v2.6.0. Finally, the selected tree was visualized using the software FigTree v1.4.3 (Rambaut, 2014a).

We used the MODEL_SELECTION Package v1.5.2 implemented in BEAST2 v2.6.2 to test alternative species tree models for Pseudopimelodus based on the results of SLSD approach. Analyses were performed using the PathSampler icon on the Package Application Launcher, with a path sampling of 10 steps, a chain length of 20 million generations and a burn-in of 10% with no pre burn-in phase. Remaining parameters were left as default following Leaché et al. (2014). Marginal likelihood estimates (MLE) were obtained for each different model run and Bayes Factor Delimitation methods were used to rank the different species delimitation models. Bayes Factors (BF) were calculated by subtracting the MLE among these models, and then multiplying the difference by two [BF = 2 × (MLE1 − MLE0)]. A positive BF value indicates support in favor of model 1, while a negative BF value indicates support in favor of model 0 (Leaché & Bouckaert, 2018).

Sequence analysis and phylogenetic relationships

Sequence edition yielded a final alignment matrix of a 357 bp and 520 bp fragment for cox1 and rag2, respectively. No evidence of indels or stop codons was found after translation into amino acids, indicating that the amplified product corresponds to functional cox1 and rag2 sequences. For the cox1 dataset, sequence analyses revealed 33 different haplotypes out of 247 sequences analyzed (GenBank accessions MH553571–MH553603; MH800618–MH800835), defined by 72 polymorphic sites, 52 of which were parsimoniously informative. Values of haplotype and nucleotide diversities for cox1 dataset were h = 0.761 ± 0.020 and π = 0.046 ± 0.001, respectively. Additionally, sequence analysis of rag2 in the trans-Andean Pseudopimelodus revealed 26 haplotypes of 77 sequences analyzed (GenBank accessions MH595749–MH595774; MH800567–MH800617), defined by 22 polymorphic sites, 16 of which were parsimoniously informative. Values of haplotype and nucleotide diversities were Hd = 0.878 ± 0.024 and π = 0.010 ± 0.000. According to the Bayesian Information Criterion, the optimal model for cox1 and rag2 datasets were HKY + G4 and HKY respectively. Collection information, voucher numbers and GenBank accessions for all new sequences obtained in this study are provided in Tables S1 and S2.

All three SLSD methods consistently suggested five unique molecular delimitations corresponding to lineages 1–5 (Fig. 2). However, we found incongruence among the methods for the delimitation of the valid species P. mangurus (EU179816) and P. charus (EU179815). Pseudopimelodus mangurus was only correctly delimited by GMYC, while P. charus was correctly delimited by locMin and GMYC (Fig. 2). The point estimates were 13 for bGMYC, 15 for locMin, and 16 for GMYC. Confidence intervals (95%) were largest for locMin (12–16 species, mean: 14), followed by GMYC (5–19 species, mean: 16), and lowest for bGMYC (4–15 species, mean: 11, median: 12, mode: 14).

Median joining network (Fig. 3) confirmed these lineages and also showed a deeper genetic divergence between lineage 1 and the remaining Pseudopimelodus lineages, as evidenced by the largest number of mutation steps (19 point mutations). Notably, lineage 1 and lineage 5, represented by a higher number of samples, showed unique haplotypes related to a common ancestral haplotype distributed in different sampling sites (Fig. 3). Network analysis showed a lineage composed of haplotypes from the Magdalena River middle sector, Magdalena River and Cauca River lower sectors and the San Jorge River (lineage 1); haplotypes from Amazon River hydrographic zone (lineage 2) and Orinoco River hydrographic zone (lineage 3); haplotypes from Atrato River and Sinú River (lineage 4); and haplotypes from the Magdalena River upper sector and, Cauca River upper and middle sectors (lineage 5).

The coalescent-based analyses using the MODEL_SELECTION Package ranked the GMYC species delimitation model as the more likely hypothesis (MLE = −4038.037016) instead of the alternatives bGMYC (MLE = −4042.52818) or locMin (MLE = −4080.28218) models. Moreover, the five evolutionary lineages of Pseudopimelodus were also supported by the dated species tree using the concatenated rag2 and cox1 genes (Fig. 4). This analysis in conjunction to all three SLSD models and phylogenetic analysis using MrBayes (Supplementary Figs. 1–3) showed that cox1 and rag2 recovered the genus Pseudopimelodus as a monophyletic group. In addition, SLSD results and dated species tree support the idea that the genus Rhyacoglanis is a sister clade of Pseudopimelodus plus Cruciglanis.

Divergence time estimates (Fig. 4) indicated that lineage 1 diverged from all remaining members of Pseudopimelodus during Middle Miocene (16.21 mya 95% HDP 11.91–20.79). Similarly, a Middle Miocene split was detected between lineage 2 and the clade comprising lineages 3, 4 and 5 (12.37 mya HPD 8.98–16.23), while a Late Miocene divergence was detected between lineage 3 and the clade formed by lineages 4 and 5 (9.93 mya HPD 6.54–13.20), and between lineages 4 and 5 (5.28 mya HPD 1.61–8.20).