MtDNA species-level phylogeny and delimitation support significantly underestimated diversity and endemism in the largest Neotropical cichlid genus (Cichlidae: Crenicichla)

Specimen sampling

This study is based on a total of 681 Crenicichla individuals (including Teleocichla) (Appendix 1, Table S1) that were collected throughout the whole distribution of the two genera in South America (Fig. 1) from 129 locations (La Plata 28, E. coast 10, W. Amazon 17, Orinoco 19, Negro 6, E. Guiana shield 24, E. Amazon 25) in 31 major river drainages (=ecoregions sensu Abell et al., 2008; La Plata 7, E. coast 6, W. Amazon 3, Orinoco 5, Negro 1, E. Guiana shield 3, E. Amazon 6).

Tissue samples for DNA analysis were taken from the right pelvic fin, stored in 95% ethanol and deposited in the tissue collection of the Department of Zoology, Faculty of Science, University of South Bohemia. Voucher specimens were photographed alive, permanently tagged and identified using a unique code (C1 to C704; Table S1), and following euthanasia (overdose of anaesthetic MS-222 Tricaine methanesulfonate and Benzocaine hydrochloride dissolved in water) preserved in 10% formalin and later transferred to 70% ethanol. Voucher specimens are deposited in the fish collection of the Department of Zoology, Faculty of Science, University of South Bohemia, and in the fish collections of Museo de La Plata (MLP), Museo Argentino de Ciencias Naturales (MACN), and Asociación Ictiológica La Plata (AI). Study of the animals was approved by the Control Commission for Ethical Treatment of Animals, University of South Bohemia, Faculty of Science, Czech Republic (17864/2005-30/300).

Specimens for the study were obtained by extensive collecting lasting several decades through our collecting and a network of dedicated cichlid hobbyists, cichlid fish importers, and colleagues who have collected and kept animals over this period. Our collecting was done under permanent collection permits to Hernán Ortega (National University of San Marcos, Lima, Peru), Javier A. Maldonado Ocampo (Pontificia Universidad Javeriana, Bogotá, Colombia), Cecilia Rodríguez-Haro (Universidad Regional Amazónica IKIAM, Tena, Ecuador), and Nadacion FUDECI (Venezuela) under bilateral agreements, and under temporary permits from the Administración de Parques Nacionales and Ministerio de Ecología y Recursos Naturales Renovables, Argentina (permits nos. NEA328 rnv3 and Res: 509/07, respectively). Specimens from Brazil and the Guianas were collected and imported by a consortium of professional aquarium fish importers coordinated by JG (Amazon Peixes Ornamentais LTDA, Manaus, Brazil; Netuno Aquarium Peixes Orn LTDA, Manaus, Brazil; Amazonstar Comercio de Peixes Ornamentais LTDA, Manaus, Brazil; FLOCAMBUS S.A.R.L, French Guiana; Oleeja Ornamental Fish Export, Suriname, South America; and Panta Rhei GmbH, Hannover, Germany) (permits nos. 0121/UVAGRO/AERO/16, 0545/UVAGRO/AERO/15, 00101/SVA-AIRJ/LPS L, 05021/SVAGRU/10, 0229 UVAGRO-ABEL 15, 0938 UVAGRO-AEROBEL 14, 0462 UVAGRO-AEROBEL 15, 0069 UVAGRO-AEROBEL 20, 011013/2016/SVAGRU, 004989/2020/VIGIAGRO, 00006331/2020-SVA-GRU, 1351 SVA-GRU 13, 2012403102A, and 201804004E). The imported specimens are kept alive and are being reared in care of JG, in a network of fish hobbyists coordinated by JG and voucher specimens are kept at the fish collection of the Department of Zoology, Faculty of Science, University of South Bohemia. Collection localities are given in Table S1 and are visualized in Fig. 1.

Species determination

We combine morphological species determination with post-hoc species delimitation using molecular mtDNA cytochrome b (cytb) and ND2 markers following the protocol of Říčan et al. (2019). Specimens were identified to species with the use of original descriptions, identification keys, and comparative material. Specimens that are new undescribed species or that could not be identified to species level are reported as “Crenicichla sp.” (possibly new/unidentified species), specimens that were identified to species only with reservation are reported as “Crenicichla cf. species” (a species that conforms to diagnosis with reservation and occurs outside its supposed distribution, possibly new) or “Crenicichla aff. species” (closely related species, revealed by DNA phylogeny, possibly new). For summary statistics, “C. aff. species” and “C. cf. species” are considered equivalent to “C. species”.

Molecular markers

For molecular phylogenetic analyses and delimitation of putative species we have used two mitochondrial (mtDNA) markers, the cytochrome b (cytb) gene and the ND2 gene following the protocol of Říčan et al. (2019). The cytb is the single most often used molecular marker in cichlid phylogenetic studies and has been one of the first molecular markers to provide well resolved phylogenies of the Neotropical cichlids that have subsequently been virtually completely confirmed by later studies including recent phylogenomic studies (e.g. Říčan et al., 2016; Ilves, Torti & López-Fernández, 2018; Burress et al., 2018a, 2018b). The ND2 gene has been found to have in Crenicichla even better and more evenly distributed support values across both shallow and deep nodes in the phylogeny than the cytb (Piálek et al., 2012).

Laboratory methods

Genomic DNA was extracted from ethanol-preserved fin tissue using the JETQUICK Tissue DNA Spin Kit (Genomed, Hannover, Germany) following standard protocol and as described in Říčan et al. (2019). The primers and reaction conditions of polymerase chain reaction (PCR) amplification are as in Říčan, Zardoya & Doadrio (2008). The products were analyzed in an ABI 3730XL automated sequencer (Applied Biosystems; Macrogen Inc., Seoul, Korea). Contiguous sequences of the gene segments were created by assembling DNA strands (forward and reverse) using GENEIOUS v. 11.0.2 (http://geneious.com, Kearse et al., 2012). Nucleotide coding sequences were also translated into protein sequences to check for possible stop codons or other ambiguities. All newly generated sequences (558) were deposited in GenBank under Accession numbers MW554711–MW584412 (Table S1). Sequences were aligned using MUSCLE v. 3.8 (Edgar, 2004), using the default settings.

Phylogenetic methods

As outgroups we have included Astronotus ocellatus and Satanoperca jurupari based on previous knowledge of cichlid relationships and studies dedicated to Crenicichla-Teleocichla and its most closely related groups (Piálek et al., 2012; Ilves, Torti & López-Fernández, 2018; Burress et al., 2018a). Maximum parsimony (MP) analysis in PAUP* 4b.10 (Swofford, 2003), maximum likelihood (ML) analyses in RAxML v8.2.4 (Stamatakis, 2014) and Bayesian inference analyses (BI) in MrBayes v3.1.2 (Huelsenbeck & Ronquist, 2001; Ronquist & Huelsenbeck, 2003) and BEAST v1.8.4 (Drummond & Rambaut, 2007) were used for phylogenetic inference following the protocol of Říčan et al. (2019). The MP phylogenetic analyses in PAUP* were run with 500 random sequence additions, 10 trees kept per addition, and a hs (heuristic) search on the saved trees to find all the shortest trees. Bootstrap analyses were done using the same approach, with five random sequence additions per one replication. Bootstrap analyses were run with 1,000 replications. The BI analyses in MrBayes (and BEAST, see below) were run with partitioning into codon positions (1st + 2nd vs. 3rd). An optimal model of evolution according to Akaike criterion was selected using MrModeltest 2.2 (Nylander, 2004) and PAUP* v. 4.0b10 (Swofford, 2003). The BI analysis using the Markov chain Monte Carlo (MCMC) simulation was run for two million generations with trees sampled and saved every 1,000 generations. Two independent analyses, each comprising two runs with eight chains, were performed to compare results of independent analyses. The analyses were run at the freely available Cipres server (https://www.phylo.org/) and in the Czech academic National Grid Infrastructure MetaCentrum (www.metacentrum.cz). The first 10% of trees from each run were discarded as burn-in. Convergence of the runs was estimated with the use of graphical visualization and diagnostics (especially the effective sample size; ESS) in Tracer v. 1.8.4 (Rambaut et al., 2018). The remaining trees were used for reconstruction of the 50% majority-rule consensus tree with posterior probability (PP) values of the branches.

Molecular clock dating analyses in BEAST

For divergence time estimation we used the Bayesian Evolutionary Analysis by Sampling Trees (BEAST) software package version v.1.8.4 (Drummond & Rambaut, 2007) with parameters as in MrBayes analyses (except for 30 million generations). We used the relaxed molecular clock model with lognormal distribution of rates and for tree prior the coalescent model with constant size. The calibration of the molecular clock was performed using secondary calibration from the study of Musilová et al. (2015) which was focused on the dating of Neotropical cichlids. Musilová et al. (2015) employed a calibration using a set of Neotropical fossil cichlid species and the study is so far the best sampled dated phylogeny of the Neotropical cichlids. For the calibration of Crenicichla we have used the basal node of the genus, estimated by Musilová et al. (2015) at 28 Ma (normal distribution; SD = 1).

The analyses were run at the freely available Cipres server (https://www.phylo.org/). Runs were checked for convergence with Tracer v.1.8.4 (Rambaut et al., 2018). Four well converged runs were combined in LogCombiner v.1.8.4 with a burnin of 10% for each of the data partition schemes. The final tree for each data partition scheme was produced from these data with TreeAnnotator v.1.8.4.

Species delimitation analyses using GMYC and PTP

We employed the General Mixed Yule Coalescent (GMYC) and Poisson tree processes (PTP) analyses for molecular species delimitation using the studied markers and following the protocol of Říčan et al. (2019). Both methods were designed for delimiting species based primarily on single molecular markers (hence where multilocus coalescent-based methods are not applicable).

The General Mixed Yule Coalescent (GMYC) model (Pons et al., 2006; Fujisawa & Barraclough, 2013) is frequently used in empirical studies (Fontaneto et al., 2007; Monaghan et al., 2009; Carstens & Dewey, 2010; Vuataz et al., 2011; Powell, 2012) and the newer Poisson tree processes (PTP) model (Zhang et al., 2013) has been shown to even outperform the GYMC method where distances between species are small. Both methods outperform OTU-picking methods (relying on simple sequence similarity thresholds) and are more robust to cases where the barcoding gap is absent (Zhang et al., 2013). Both methods model speciation (among-species branching events) via a pure birth process and within-species branching-events as neutral coalescent processes. The methods identify the transition points between inter- and intra-species branching rates by maximizing the likelihood score of the model. While the GMYC method uses time to identify branching rate transition points (hence only on a time-calibrated ultrametric tree) the PTP method directly uses the number of substitutions and does not require a time-calibrated ultrametric tree. Both methods assume that all lineages leading from the root to the transition points are different species.

PTP and GMYC analyses were run at the freely available web interface (http://species.h-its.org and https://mptp.h-its.org). Prior to analyses a haplotype dataset was created in Fabox (Villesen, 2007). The mPTP and bPTP analyses were run on the MrBayes haplotype tree and the GMYC analysis on the ultrametric haplotype BEAST tree.

Phylogenetic relationships

Prior to this study the phylogenies with the most representative sampling of Crenicichla and Teleocichla were the studies of Piálek et al. (2012) and Burress et al. (2018a). Piálek et al. (2012) included 161 specimens and 35 species in a multilocus mtDNA + nDNA phylogeny, while Burress et al. (2018a) included 57 species in a reduced genome representation nDNA ddRAD phylogeny. At the species level the relationships in Piálek et al. (2012) and Burress et al. (2018a), derived from different molecular markers, are highly congruent and are also highly congruent with our present study with a much denser sampling, the only exceptions are within the Missioneira group and parts of the Saxatilis group which lack statistical robustness of resolution in all three studies.

At the level of the species groups Teleocichla is in both Piálek et al. (2012) and Burress et al. (2018a) nested within Crenicichla as the sister group of the Wallacii, Saxatilis and Lugubris groups, in Burress et al. (2018a) with all analyses of the ddRAD data using various settings, parameters and matrix sizes giving a 100% bootstrap support for this phylogenetis position. In our mtDNA phylogeny in this study the position of Teleocichla is equivocal, either nested within Crenicichla (MP and ML analysis; Figs. S1, S2A and S2B), or as a sister group of Crenicichla (BI analysis; Fig. 2), in both cases without statistical support. The uncertainty of Teleocichla position in mitochondrial markers is likely due to saturation at this large phylogenetic distance (see Piálek et al., 2012 for details). Based on available studies (Piálek et al., 2012; Burress et al., 2018a; this study) Teleocichla is thus most probably a species group within Crenicichla, the sister group of the Wallacii, Saxatilis and Lugubris groups (including Acutirostris group) (Burress et al., 2018a).

The southern groups of Crenicichla (Missioneira, Lacustris, Scottii) form a well supported monophyletic clade in our mtDNA analyses (Figs. 2, 5, 6, S1, S2A and S2B) as well as in other studies of Crenicichla (Piálek et al., 2012; Burress et al., 2018a) with the Missioneira being the sister group of the Lacustris. The Scottii group is found strongly nested within the Lacustris group as it is in Piálek et al. (2012) and Burress et al. (2018a).

The sucessive sister groups of this southern clade are the Macrophthalma group and the Reticulata group, both relationships being well supported in our mtDNA analyses (Figs. 2, 4–6, S1, S2A, S2B) as they are in the mtDNA + nDNA multilocus study of Piálek et al. (2012), while in the nDNA ddRAD study of Burress et al. (2018a), the Reticulata group is the sister group of the southern clade followed by the Macrophthalma group (bootstrap supports for this alternative relationship have values of 60%, 79%, 86%, 88% and 99% in the various analyses of the ddRAD data in Burress et al., 2018a).

The remaining species groups (Wallacii, Saxatilis, Lugubris and Acutirostris) form the second main clade of Crenicichla in the phylogenetic BEAST analysis of our mtDNA data (Figs. 2–4), and with the exclusion of the Wallacii group also in the RAXML analysis (Figs. S1, S2A and S2B). The relationships between the Saxatilis, Lugubris and Acutirostris groups however differ between the BEAST and RAXML analyses and lack statistical support. Both Piálek et al. (2012) and Burress et al. (2018a) found Lugubris and Acutirostris as the sister group of Saxatilis followed by Wallacii, which is a different relationship from both our mtDNA phylogenies with lack of statistical support at these nodes (Figs. 2–4, S1, S2A and S2B), while the topology in both Piálek et al. (2012) and Burress et al. (2018a) is strongly supported.

The Acutirostris group (and the Scottii group; see above) is the only species group whose species composition and monophyly is questioned by molecular phylogenetic analyses. In our mtDNA analyses the Acutirostris group appears as a separate clade but without statistical support that includes C. acutirostris, C. multispinosa, C. ternetzi and also C. cincta (Figs. 2–4, S1, S2A and S2B). In the nDNA ddRAD phylogeny of Burress et al. (2018a) the group is not monophyletic but is paraphyletic to the Lugubris group (C. acutirostris at the base followed by C. multispinosa followed by C. cincta) with strong statistical support and hence probably best subsumed into the latter group, as has been proposed in Stawikowski & Werner (2004), where it is treated as a species complex within the Lugubris group.

In summary our analyses together with previous molecular studies of Crenicichla support the monophyly of the Lacustris, Lugubris, Macrophthalma, Missioneira, Reticulata, Saxatilis, and Wallacii species groups which form two main clades of Crenicichla, the Lacustris-Missioneira-Reticulata-Macrophthalma clade, and the Lugubris-Saxatilis-Wallacii clade, and we subsume the Acutirostris group into the Lugubris group and the Scottii group into the Lacustris group.

Species diversity

Our results clearly demonstrate that the species diversity in Crenicichla is significantly underestimated. While our study included 77 valid species our analyses delimited between 126 and 168 species (Figs. 3–6), at least twice as many. Extrapolation from the ratio of total valid species (105) and the number of valid species included in our study (77) (105/77) and the here delimited species (126–168) suggests the minimum number of species in Crenicichla to be 172–229 (Figs 3–6).

While the profesional ichthyological knowledge of the species diversity in the group is very limited as seen from our results, the cichlid-hobby knowledge of the Crenicichla fauna is actually much better. This is evident from the correspondence between the here delimited species and the potentially undescribed species identified in the cichlid hobby literature. The cichlid hobby literature identifies at least 85 potentially undescribed species of Crenicichla and Teleocichla (Appendix 3) and essentially all of those that were included in the present study (50) have been delimited by the molecular analyses (40–46 depending on the delimitation method; Figs. 3–6; this group of delimited taxa shown in violet). Adding the at least 85 identified potentially undescribed species to the presently valid species (105) still gives an underestimated total of 190 Crenicichla plus Teleocichla species compared to our extrapolation based on molecular species delimitation which is in the range of 172–229 species. Virtually all of these species unrecognized by the professional literature are found in the Amazon-Orinoco-Guiana core of South American fish diversity (Albert & Reis, 2011; Van der Sleen & Albert, 2018), where they in the C. lugubris and C. saxatilis groups actually form the majority of the delimited species diversity (Figs. 3–6; delimited species 31–41, 32–51 vs. included valid species 12, 19, respectively). Many of the delimited species in this category are presently subsumed under described species and especially under those with putatively large distribution areas. Our mtDNA phylogeny and species delimitations together with their various degrees of morphological distinctiveness clearly demonstate that these putatively widely distributed species are an oversimplification that does not reflect true diversity.

In the following paragraphs we provide a short review of the putative new species previously proposed in the cichlid-hobby literature and here for the first time delimited and supported in molecular analyses. We provide this review by species groups, starting with the C. lugubris and C. saxatilis groups, where this hidden diversity actually forms the majority of the species diversity.

Underestimated Neotropical fish biodiversity

The minimum percentage increase of possibly unrecognised species observed here in Crenicichla (65–121%) significantly surpasses the conclusions of Reis et al. (2016), who estimated that 34–42% of Neotropical freshwater fishes remain undescribed. The percentage increase identified here varies widely between Crenicichla species group from 0% in the Missioneira and Macrophthalma groups, through 25–40% in the Lacustris group, 50–87% in the Reticulata group, 60–80% in Teleocichla, 68–168% in the Saxatilis group, 125–200% in the Wallacii group, and 158–241% in the Lugubris group (Fig. 1). The numbers clearly indicate that the supposedly widely distributed species in the Reticulata, Wallacii, Saxatilis and especially Lugubris group have no foundation in the studied molecular markers. In agreement with Reis et al. (2016) the localization of the unrecognized diversity in Crenicichla is also in the Amazon basin, especially in its eastern portion and in the remaining two areas of the AOG core (Greater Amazonia) of South American ichthyological diversity.

The general explanations for the unrecognised diversity of Neotropical fishes stem largely from (1) lack of systematic sampling throughout the distibution of the given group and lack of sampling in areas of difficult access, (2) widespread taxa or heterogeneous taxa with insufficient or overwhelming amounts of museum material or improperly preserved material and information, and/or (3) cryptic or pseudocryptic (morphological differences apparent but overlooked) diversity in widespread species (Reis et al., 2016). Genetic data are an important instrument in uncovering cases of the latter. In Crenicichla the explanation for the underestimated diversity is the same with two additions, namely (1) limitations in available diagnostic characters due to predominant study of preserved specimens which have lost the diagnostic coloration pattern characters only present in live fishes (this is a specific cichlid problem among all other Neotropical fish groups), and (2) lack of acknowledgement in the professional ichthyological literature of the effort brought by the hobbyists, in this case the cichlid amateurs. As our study demonstrates, the majority of the genetically delineated putatively new species are not unexpected units (of these there are conservatively 25 in our study), but rather previously proposed potentially undescribed species by the cichlid hobbyists (of these there are conservatively 39 in our study).

Therefore, we stress the conclusions of Reis (2013), who demonstrated that the primary instruments for recognition, understanding, and conservation of Amazon basin fishes is increasing expertise in fish taxonomy and systematics and we additionally clearly demonstrate that DNA study and other genetic tools are powerful complementary methods for uncovering fish diversity and highlighting groups in need of taxonomic revisions.

Despite our wide ranging specimen sampling and potential for understanding of Crenicichla diversity our study has one limiting factor and that is the use of a single uniparentally inherited molecular marker (mtDNA cytb and ND2 loci). The potential limitations for species delimitation using mtDNA have been discussed in literature (e.g., Rannala & Zang, 2003; Yang & Rannala, 2010; Dupuis, Roe & Sperling, 2012; Fujita et al., 2012). The use of additional nuclear markers is generally recommended as the use of these unlinked markers has the potential to improve the accuracy of phylogenetic reconstructions and species delimitation. A study doing just that using the same specimen sampling and the nDNA reduced genome representation ddRAD method is currently in preparation and shows results that are highly congruent with the mtDNA markers employed here.

In spite of having used only mtDNA loci, the majority of molecularly delimited putative species are actually not previously unrecognized but previously synonymized or already postulated species based on coloration patterns, morphological characters and inferred endemism due to potential distribution barriers. All these additional data support and reinforce the inference based on the mitochondrial markers.