Divergence times of the Rhoadsia clade (Characiformes: Characidae)

The family Characidae is the most diverse group of fishes in the Neotropics with challenging systematics. The three genera Carlana , Parastremma , and Rhoadsia , formerly considered the subfamily Rhoadsiinae, are now included in the subfamily Stethaprioninae. Previous phylogenetic analyses did not include all genera of Rhoadsiinae, specifically Parastremma . Here, we estimated the phylogenetic relationships and divergence times of the genera of Rhoadsiinae (the Rhoadsia clade) relative to the most representative genera of the Characidae . We used six molecular markers from the mitochondrial and nuclear genome to estimate the phylogeny and divergence times. We confirmed the monophyly of the Rhoadsia clade. Furthermore, we estimated that the Central American genus Carlana and the western Colombian genus Parastremma diverged approximately 13 Mya (95% HPD 8.36–18.11), consistent with the early-closure estimates of the Isthmus of Panama (~15 Mya). The genus Rhoadsia , endemic to Western Ecuador and Northern Peru, was estimated to originate at around 20 Mya (95% HPD 14.35–25.43), consistent with the Andean uplift (~20 Mya).


INTRODUCTION
The Family Characidae is the most diverse family of fishes in the Neotropics (Albert, Reis, 2011;Fricke et al., 2022). Due to its great diversity, species of this large group are classified into various subfamilies. However, classifying the species into subfamilies is still challenging and constantly changing as new information becomes available. The use of molecular markers in combination with morphology has helped clarify a lot of the uncertainty within Characidae. That is the case of the former subfamily Rhoadsiinae with the three genera Rhoadsia Fowler, 1911, Parastremma Eigenmann, 1912, and Carlana Strand, 1928(Cardoso, 2003 (here also referred as the Rhoadsia clade). However, more recent total-evidence phylogenetic analysis prompted the reclassification of this group into a larger subfamily Stethaprioninae (Mirande, 2019), which is consistent with phylogenomic evidence (Betancur-R. et al., 2019).
The subfamilial recognition and membership of the Rhoadsia clade have shifted over time. For example, the Rhoadsia clade includes the genus Rhoadsia with two species recognized, R. minor Eigenmann &Henn, 1914 andR. altipinna Fowler, 1911, the genus Carlana with its only species C. eigenmanni (Meek, 1912), and the genus Parastremma with three species, P. album Dahl, 1960, P. pulchrum Dahl, 1960, and P. sadina Eigenmann, 1912. Cardoso (2003 defined the group as having a single series of teeth in the premaxilla when young and a two series when reaching adulthood (except for Carlana which does not develop an outer series). Adult specimens have two conical teeth in their outer series and five multicuspid teeth in their inner series. Mirande (2009Mirande ( , 2010 proposed the inclusion of Nematocharax in the subfamily Rhoadsiinae sensu Cardoso (2003), based on morphological phylogenetic analyses. The characters that unified the four genera were the form of teeth of inner premaxillary row with cusps aligned in straight series, without anterior concavity and five or more cusps of anterior maxillary teeth (Mirande, 2010). By contrast, a more recent phylogeny based on combined morphological and molecular data showed that these four genera are not monophyletic. That hypothesis recovered Nematocharax Weitzman, Menezes & Britski, 1986 closer to the polyphyletic genus Moenkhausia Eigenmann, 1903(Mariguela et al., 2013Mirande, 2019), a group of medium size fish (~12 cm) widely distributed in the Amazon basin and adjacent drainages, which were classified as part of a tribe Stethaprionini. On the other hand, the genera Rhoadsia and Carlana were closer to species like Pseudochalceus kyburzi Schultz, 1966and Nematobrycon palmeri Eigenmann, 1911(Mirande, 2019 mainly found in the northwestern South America; together, these were classified as members of a tribe Rhoadsiini. Other members of this Rhoadsiini tribe can also be found in the Amazon basin and Southeast Brazil. Consequently, these genera are now within the much larger subfamily Stethaprioninae (Mirande, 2019). Although Parastremma was generally assumed to be within the Rhoadsia clade with Rhoadsia and Carlana, Parastremma has not been formally included in a phylogeny until recently, where it was used as an outgroup of populations of Rhoadsia sp. along with Carlana (Cucalón et al., 2022). However, Cucalón et al. (2022) did not sample other characid taxa to test the monophyly of the Rhoadsia clade, nor did they perform a fossil-calibrated divergence time analysis to estimate when members of the group diverged from each other.
In this study, we investigated the phylogenetic relationships and divergence time of the Rhoadsia clade (Rhoadsiinae sensu Cardoso, 2003) relative to other members of the family Characidae, with the intension to provide a better understanding of the evolutionary history of the Rhoadsia clade within the Characidae.

MATERIAL AND METHODS
Data collection. Sequences for representative taxa within families of Characoidei (sans Crenuchoidea) were retrieved from GenBank using phylogenies estimated in Mirande (2019) and Terán et al. (2020) as guides to maximize phylogenetic diversity. Most extensive sampling was done within the family Characidae to achieve representative sampling of most clades within the subfamilies. In addition, clades outside of Characidae were represented with up to 6 species per family. We retrieved sequences representing three mitochondrial markers, 16S Ribosomal RNA (16S) (~600 bp), Cytochrome Oxidase I (COI) (~600 bp), and Cytochrome b (Cytb) (~600 bp), and up to three nuclear markers, Myosin Heavy Chain 6 (Myh6) (~1000 bp), Recombination Activating 1 (RAG1) (~1200 bp), and Recombination Activating 2 (RAG2) (~1200 bp), when available from GenBank. These markers were chosen since they are commonly used for phylogenetic reconstructions and were the most frequently available across taxa. Genes available from different individuals were chosen arbitrarily for each species, however whenever possible genes coming from the same individual were selected to reduce the likelihood of contamination (see Tab. S1). Mitochondrial genes were obtained from complete mitogenomes when available using a custom script ( (Palumbi, 1996), RAG1 (López et al., 2004;Li, Ortí, 2007;Oliveira et al., 2011), and RAG2 (Oliveira et al., 2011; this study) (see Tab. S2). The polymerase chain reaction (PCR) was carried out using the following volumes: For a 15 µl reaction we used 7.5 µl of GoTaq® master mix (www.promega.com), 0.3 µl of each primer at 10 µM (Tab. S2), 2 µl of DNA template, and complemented the reaction with molecular grade water. We performed a nested PCR for RAG1, and RAG2, and a regular PCR for the 16S, with the following thermocycler protocol. For the first PCR nested and the regular PCR, one cycle of denaturation for 1 min at 95°C, 30 cycles of denaturation at 95°C for 1 min, annealing (50-52°C) (Tab. S2), extension at 72°C for 2 min, and one cycle of final extension at 72°C for 10 min. The second PCR nested was similar as the first PCR except it ran for 35 cycles.

Alignment.
Each gene was aligned independently using MAFFT version 7.453 (Katoh, Standley, 2013) using the option -auto recommended when unsure which alignment strategy to use based on data size. Then, the aligned genes were concatenated and converted to NEXUS format for further analyses using the tool AMAS (Alignment Manipulation And Summary) (Borowiec, 2016).
Phylogenetic analysis. We used Maximum Likelihood (ML) to reconstruct the phylogeny of the Rhoadsia clade relative to the family Characidae. Species from the other families within Characoidei (sans Crenuchoidea) were used as outgroup taxa. The ML analysis was carried out using IQ-TREE2 v. 2.0.6 (Minh et al., 2020), using a partitioned model (Chernomor et al., 2016), with 1000 iterations of ultrafast bootstrap (Hoang et al., 2017). To determine the best partitioned substitution model for phylogenetic analysis, we implemented ModelFinder (Kalyaanamoorthy et al., 2017) to simultaneously estimate each gene's best-fit substitution model and best-fit alignment partitioning model scheme (option --merge). ModelFinder selects the best-fit model that minimizes the Bayesian Information Criterion (BIC) score. We enforced the relationship of some of the deep nodes based on phylogenomic results from Betancur-R et al. (2019) as followed: Chalceidae + Characidae, sister to Alestoidea and Erythrinoidea + Curimatoidea. The ML tree was edited for visualization using the R packages phytools v1.0-1 (Revell, 2012) and ggtree v3.2.1 (Yu et al., 2017).

Divergence time estimation.
Divergence times among clades were estimated using BEAST version 2.6.4 (Bouckaert et al., 2019). We constrained the analysis using a starting tree with estimated divergence time based on the penalized likelihood (PL) method (Cole et al., 2014) implemented in the software treePL (Smith, O'meara, 2012) following Maurin (2020). We used the rooted ML tree generated from IQ-TREE as the input tree for treePL and calibrated the nodes using 13 calibration nodes used in Kolmann et al. (2021). See Tab. S3 for detailed information about the node calibration including, taxa calibrated, mean age, minimum and maximum age, fossil species name, and tips (i.e., species in the tree) used to inform treePL which node to calibrate. Most parameters were set through the program BEAUti version 2.6.4 (Bouckaert et al., 2019) included in the BEAST2 package. The parameters were as followed: We assumed that all genes had the same molecular clock and tree topology by linking the "Clock" and "Tree" model for all partitions. The "Site model" was left unliked across partitions to use the best-fit substitution model. The substitution model was estimated using ModelFinder (Kalyaanamoorthy et al., 2017) from IQTREE2 v. 2.0.6 (Minh et al., 2020) using the option "TESTONLY" to exclude testing of the free rate model (assumed by default), since this is not currently implemented in BEAST2. The option "--merge" was also set to estimate the best substitution model scheme to optimize the number of parameters used during the analysis. The molecular clock model was set as an uncorrelated relaxed clock log-normal to allow for rate heterogeneity across branches. We used the birth-death (Gernhard, 2008) model tree prior and calibrated the nodes of the phylogeny using prior settings following the 13 calibration fossils used in Kolmann et al. (2021) (Tab. S3). For details on fossil calibration and rationale see materials and methods described in Kolmann et al. (2021). The prior distribution for each calibration was set as exponential to account for increasing uncertainty at further points in the past, except for the root node that was uniform distribution, following Kolmann et al. (2021). To fix the tree topology during the Markov Chain Monte Carlo (MCMC) chain, we modified the XML file by removing the lines that contained the operators "wideexchange," "narrow-exchange," "subtree-slide," and "Wilson-balding" following the instruction from the BEAST2 website (www.beast2.org). We ran the analysis with 100 million MCMC iterations, sampling every 5000 generations. The log file was inspected in Tracer (Rambaut et al., 2014) for convergence of the MCMC. We summarized the maximum clade credibility (MCC) tree discarding 10% burn-in in TreeAnnotator from the BEAST package version 2.6.3 (Bouckaert et al., 2019). The time calibrated (i.e., MCC tree) and ML trees were visualized in FigTree (Rambaut, 2016). We used the R package MCMCtreeR (Puttick, 2019) to visualize the posterior ages distribution of the MCC tree.

RESULTS
Sampling and phylogenetic analysis. Sequences for a total of 211 species of the Characoidei group were obtained for the phylogenetic reconstruction. The total length of the concatenated sequence alignment after trimming was 8276 bp. Accession number of the gene sequences used including the ones obtained from this study for the Rhoadsia clade can be found in Tab. S1. The best scheme substitution models selected by ModelFinder for the ML method allowing the incorporation of the free rate model were GTR+F+R5 (16S), GTR+F+I+G4 (COI), GTR+F+I+G4 (Cytb), and TIM2e+R5 (Myh6, RAG1, RAG2). For the Bayesian divergence time analysis (without allowing for free rate model), the best substitution model scheme was GTR+F+I+G4 (16S), GTR+F+I+G4 (COI), GTR+F+I+G4 (Cytb), and TIM2e+I+G4 (Myh6, RAG1, RAG2). The MCMC run reached stationarity with ESS greater than 200.
Phylogenetic relationships of the Rhoadsia clade. The ML tree placed the Rhoadsia clade (ultrafast bootstrap [BS]: 100) within the subfamily Stethaprioninae ( Fig.  1). Rhoadsia was reconstructed as monophyletic (BS: 100), and this clade was sister to a highly supported clade formed by Carlana eigenmanni and Parastremma sadina (BS: 100). The closest relative to the Rhoadsia clade was the species Pseudochalceus kyburzi (BS: 100). That clade was sister of Nematobrycon palmeri (BS: 99) and followed by Inpaichthys kerri Géry & Junk, 1977 (BS: 83). Refer to Fig. S4 for full ML tree and BS support values for all nodes.   Fig. S5 for full maximum clade credibility tree with 95% HPD for all nodes.

DISCUSSION
We herein inferred the phylogenetic relationships and divergence times of the Rhoadsia clade (Rhoadsiinae sensu Cardoso, 2003), based on six molecular markers from the mitochondrial and nuclear genome. The phylogenetic relationships of Rhoadsia and Carlana to other characids were consistent with previous reports (Mirande, 2019;Terán et al., 2020). The three genera of the Rhoadsia clade -Rhoadsia, Carlana, and Parastremma -showed a monophyletic relationship ( Fig. 1)  Relationships of the Rhoadsia clade and closest relative within the subfamily Stethaprioninae. The phylogenetic relationships of the Rhoadsia clade were consistent with the most recent phylogeny of the group based on morphological and molecular data (Mirande, 2019;Terán et al., 2020). In this study, we added Parastremma sadina to the analysis. The genus Parastremma with three valid species, P. album, P. pulchrum, and P. sadina (the latter analyzed here) has been previously designated within Rhoadsiinae sensu Cardoso (2003) based on morphology but without formal phylogenetic analysis (Cardoso, 2003;Mirande, 2010;Jimenez-Prado et al., 2015). The genus Parastremma is endemic to the Chocó-Darien ecoregion in Colombia (DoNascimiento et al., 2017). The phylogeny presented here shows Parastremma sadina sister to the monotypic Carlana from Central America, corroborating its place within the Rhoadsia clade (Fig. 1). The genus Carlana with its only species Carlana eigenmanni has been regarded previously as a junior synonym of Rhoadsia by some authors (Eigenmann, Myers, 1921;Géry, 1977) while other authors associated Carlana with the subfamily Cheirodontinae after observing that Carlana was the only member of the Rhoadsia clade keeping a single tooth series in the premaxilla in adult fish as opposed to a double series (Fink, Weitzman, 1974). However, this trait appears to be a homoplasy.
The genus Rhoadsia contains two recognized species, R. minor, and R. altipinna. This genus is endemic to drainages from the Pacific slope of the Andean mountains from northern Ecuador to northern Peru, an area known for being highly threatened (Aguirre et al., 2021). The most distinctive characteristic of Rhoadsia is a dark spot located on the side of the body below the insertion of the dorsal fin (Jimenez-Prado et al., 2015). Although their taxonomic status as two species have being questioned by some authors (Géry, 1977), recent studies based on molecular markers showed the two species appear to differ genetically and are allopatrically distributed (Aguirre et al., 2016;Cucalón et al., 2022) Aguirre et al., 2019). Nevertheless, the taxonomic distinction within the genus Rhoadsia might still require further study using a more integrative approach.
The genus Nematocharax, a Brazilian freshwater fish that was previously thought to be closely related to the Rhoadsia clade (Weitzman et al., 1986;Mirande, 2009Mirande, , 2010, has shown to be problematic in regards to its systematics due to morphological similarities with various characid genera (Weitzman et al., 1986). Weitzman et al. (1986) suggested the potential relationship of Nematocharax with the Rhoadsia clade due to the compressed, relatively deep body, long dorsal fin, and fully toothed maxilla, which are characteristics found in the Rhoadsia clade. The phylogenetic relationship was then supported in Mirande (2010), uniting Nematocharax with members of the Rhoadsia clade by the form of the teeth of the inner premaxillary tooth row with cusps aligned in 9/14 ni.bio.br | scielo.br/ni Roberto V. Cucalón and Milton Tan straight series and without anterior concavity and five or more cusps in the anterior maxillary teeth. However, the designation of Nematocharax as a member of the Rhoadsia clade was challenged after the inclusion of the nuclear and mitochondrial markers into the phylogeny, placing Nematocharax sister to species of the polyphyletic genus Moenkhausia (Mariguela et al., 2013;Mirande, 2019; this study).
Major historical processes associated with the divergence of the Rhoadsia clade. The high levels of endemism found in Western Ecuador have been attributed to species isolation due to the uplift of the Andes estimated at 20 Mya (Schaefer, 2011). In addition, the effect of the uplift of the Andes has been investigated on the diversification of birds (Benham, Witt, 2016;Hazzi et al., 2018) and plants (Luebert, Weigend, 2014). More recently, it has been associated with the high diversity of the freshwater fish family Characidae in South America (Melo et al., 2022) and South American freshwater fishes in general (Cassemiro et al., 2021). Freshwater fishes are especially prone to diversify after such geological events due to their limitation to move outside their river systems. The subfamily Stethaprioninae has being previously estimated to have diversified close to 30 Mya (Melo et al., 2022). Other species-rich families like Trichomycteridae and Loricariidae (Siluriformes) are associated with major geological events in South America including multiple marine transgressions and regressions as well as mountain formations between 55-10 Mya (Cassemiro et al., 2021). The origin of the genus Rhoadsia seems consistent with these estimations, showing a time of divergence from its closest relative (14-25 Mya) (Fig. 1), within the range of the formation of the Andes (Schaefer, 2011).
The genera of the Rhoadsia clade each inhabit adjacent, non-overlapping regions (Western Ecuador, Colombia, and Central America). The genus Carlana is found only in Central America from Panama to Nicaragua (Cardoso, 2003). The genus Parastremma inhabits Colombia's freshwater rivers on both the Caribbean slope and Pacific slope, and is the closest geographically to Rhoadsia from Western Ecuador (Eigenmann, 1912;Cardoso, 2003;DoNascimiento et al., 2017). Both Parastremma and Rhoadsia are part of the Tumbes-Chocó-Magdalena ecoregion spanning from Southeastern Panama to Northwestern Peru, but their species ranges do not overlap. Fish faunas from Western Ecuador and Western Colombia are known to differ, potentially indicating independent evolutionary histories (Eigenmann, 1921;Aguirre et al., 2017). Eigenmann (1921) suggested that Rhoadsia and Parastremma independently dispersed from the east of the Andes, although other hypotheses appear equally likely (i.e., north-south or south-north origin). It is noteworthy that species like Hoplias microlepis (Günther, 1864), which inhabit Western Ecuador like Rhoadsia, are also found in the Chagres River, Panama, while absent in Colombia, presenting a disjoint distribution (Eigenmann, 1921;Mattox et al., 2014). This may indicate that the subdivision in fish composition observed between Western Ecuador, Western Colombia, and Central America does not seem universal for all fishes of this region. However, genetic analysis is still needed to determine the level of divergence between the two disjunct populations of Hoplias Gill, 1903. The genus Parastremma and the Central American genus Carlana diverged ~13 Mya (95% HPD 8. 36-18.11). This is consistent with recent estimates of an ancient closure of the Isthmus of Panama (~15 Mya) and migration patterns of other land and freshwater organisms into Central America (Hurt et al., 2009;Bacon et al., 2015a;Thacker, 2017), rather than the traditionally accepted closure at 3.5 Mya (reviewed in Coates and Stallard, 2013). Furthermore, dispersal events prior to the full closure of the Isthmus of Panama appear to be commonly inferred across taxa (Bermingham, Martin, 1998;Thacker, 2017), and the divergence between Carlana and Parastremma seems to be more consistent with an early closure of the Isthmus as well (Coates, Stallard, 2013;Bacon et al., 2015a;Montes et al., 2015). We also observed a similar divergence time between the species Astyanax microlepis Eigenmann, 1913 found exclusively in Colombia and its sister clade that includes A. mexicanus (De Filippi, 1853), A. petenensis (Günther, 1864), A. nasutus Meek, 1907, andA. nicaraguensis Eigenmann &Ogle, 1907 from North and Central America (~20 Mya, 95% HPD 11.79-23.17) (Fig. 1). The timing of the closure of the Isthmus of Panama is still an ongoing debate (Bacon et al., 2015a,b;Hoorn, Flantua, 2015;Montes et al., 2015;O'Dea et al., 2016;Jaramillo et al., 2017;Molnar, 2017) and we cannot discard the possibility of the divergence between Carlana and Parastremma happening in South America first followed by the dispersal of an ancestral population of Carlana into Central America, followed by extinction. Given the Gondwanan origin of characiforms (Arroyave et al., 2013), it is safe to assume that the genus Carlana should have expanded up into Central America not earlier than 18 Mya at the beginning of the Miocene. It has been hypothesized that later geological events like the Pliocene high sea stand and the rise of the Central Cordillera might have contributed to the extirpation of some species in Central America, allowing others to occupy some new niches favoring allopatric speciation and the high endemism found in the lower Mesoamerican region (Smith, Bermingham, 2005).
In this study, we include a member of the genus Parastremma in a densely-sampled molecular phylogeny of the Characidae for the first time, supporting a relationship sister to Carlana and both closely related to Rhoadsia as a clade, supporting the monophyly of the Rhoadsia clade. We also estimated the divergence times of the group, with the Ecuadorian genus Rhoadsia diverging from its closest relative between 14-25 Mya, consistent with the uplift of the Andean Mountains. The Central American Carlana and Colombian Parastremma diverged between 8-18 Mya, consistent with the early closure hypothesis of the Isthmus of Panama. This study adds knowledge regarding the evolutionary history and biogeography of the Rhoadsia clade.

ACKNOWLEDGMENTS
We thank Dr. Mark Davis for granting access to the molecular laboratory and supplies for PCR work, Dr. Windsor Aguirre for providing the DNA samples used to complement the loci of members within the Rhoadsia clade, Rachel Skinner for providing a python script to retrieve multiple accessions from GenBank, and Jeffrey Haas for allowing us to use the Department of Agricultural and Biological Engineering server at UIUC, which was essential for running the analyses. We thank Dr. Marcos Mirande and an anonymous reviewer for their valuable comments and suggestions that help improve the quality of this manuscript. • Maurin KJL. An empirical guide for producing a dated phylogeny with treePL in a maximum likelihood framework. ArXiv. 2020.