Karyotype Evolution of Talking Thorny Catfishes Anadoras (Doradidae, Astrodoradinae): A Process Mediated by Structural Rearrangements and Intense Reorganization of Repetitive DNAs

Neotropical freshwater fishes represent one of the most diversified animal groups on earth due to their remarkable morphological and chromosomal diversity, and they are considered an excellent model for evolutionary studies [Machado et al., 2018; de Moraes et al., 2021; Favarato et al., 2021; Sassi et al., 2021; Takagui et al., 2021]. Over the past years, comparative cytogenetic studies have generated relevant information that contributes to the understanding of phylogenetic relationships, taxonomy, and diversification of some Brazilian freshwater fish lineages [Ferreira et al., 2017, 2020; Prizon et al., 2017; Usso et al., 2018]. However, a major obstacle in using this type of approach is the lack of cytogenetic data for some groups, or when they are available, they are either fragmented or restricted to certain clades. Thus they are poorly representative in relation to the total diversity of a given group, as seen for instance within Doradidae thorny catfishes [Baumgärtner et al., 2018; Takagui et al., 2019]. This family has upwards of 96 valid species [Fricke and Eschmeyer, 2021] widely distributed along South America, where they are popularly known as armaus, abotoados, rique-rique, serrudo, among others. These species are easily recognized by a row of retracted spines formed from the hypertrophy of the lateral line tubules and a cephalic shield and well-developed spines in the dorsal and pectoral fins [Ferraris, 2007].

Currently, the Doradidae family is considered a monophyletic group by both morphological and molecular approaches, with phylogenetic relationships relatively well-resolved. Six subfamilies are recognized: Wertheimerinae, Astrodoradinae, Doradinae, Achantodoradinae, Agamyxinae, and Rhinodoradinae [Sabaj and Arce, 2021]. The subfamily Astrodoradinae was first proposed by Higuchi et al. [2007] and actually comprises 15 species [Fricke and Eschmeyer, 2021] allocated in 6 genera: Anadoras, Amblydoras, Scorpiodoras, Physopyxis, Hypodoras, and Astrodoras [Sabaj and Arce, 2021]. The speckled Anadoras thorny catfishes have a wide geographic distribution, however only 2 species are formally recognized: Anadoras grypus (plus Anadoras insculpitus), which occurs in the Upper and Middle Amazon, from the large Andean tributaries located in Ecuadorian and Peruvian territory to Central Amazonia in Brazil, and Anadoras weddellii (plus Anadoras regani), described from Mamoré, Guaporé, Madeira, Jari, Oiapoque Rivers, Lower Amazon to the Marajo island as well as to the Pantanal basin [Birindelli and Sousa, 2017]. A third candidate species that was temporarily named as Anadoras sp. “araguaia” occurs exclusively in the Upper Araguaia River basin, one of the most important Brazilian hydrographic systems that houses a unique fish fauna [Claro-García and Shibatta, 2013; Jardulli et al., 2014].

The present study consists of a comparative cytogenetic analysis (e.g., diploid number, heterochromatin pattern, and physical mapping of rDNAs and telomeric probes) and includes all Anadoras species from the Amazon basin and Brazilian Pantanal. Here, we add new pieces to the complex evolutionary puzzle of the thorny catfishes by analyzing the karyotypes of the oldest divergent genera within the Astrodoradinae subfamily [Sousa, 2010; Arce et al., 2013]. Additionally, a species delimitation analysis based on a 600-bp fragment of the cytochrome oxidase C subunit 1 (CO1) gene was performed for a better molecular diagnosis of the 3 Anadoras species, currently recognized only by morphological taxonomy.

The following species were analyzed: A. grypus (5 females and 2 males) from Catalão Lake, Manaus City, Amazonas State, Brazil (3°09′49.8′′S, 59°54′47.5′′W); 2 populations of A. weddellii from temporary lagoons of Miranda River (4 females and 5 males), Corumbá City, Mato Grosso do Sul State, Brazil (19°31′24.96′′S, 57°02′25.51′′W) and Cuiabá River, Cuiabá, Mato Grosso State, Brazil (16°12′46.12′′S, 55°59′09.05′′W); and Anadoras sp. “araguaia” (5 females and 3 males) from Córrego do Medo, a small tributary of Araguaia River basin, in São Miguel do Araguaia City, Goiás State, Brazil (3°07′06.6′′S, 50°26′55.9′′W) (Fig. 1). The collection of the specimens was performed under the authorization of ICMBio (Instituto Chico Mendes de Conservação da Biodiversidade). The specimens were deposited in the Museum of Zoology of the Londrina State University (MZUEL) under the voucher numbers MZUEL 17811 (A. grypus) and MZUEL 11085 (A. weddellii). Specimens of Anadoras sp. “araguaia” were deposited in the Museum of Zoology of the Universidade de São Paulo (MZUSP) under the voucher number MZUSP 110809.

Chromosome preparations were obtained from kidney cells according to the air-drying technique [Bertollo et al., 1978]. Previously, a mitotic stimulation with 300 µL of Broncho-vaxom (bacterial lysate 0.2 g/mL of water) was injected to trigger an inflammatory process and increase the number of renal cells in mitotic division [Molina et al., 2010]. To determine the diploid number, we analyzed at least 25 metaphase plates per individual in all Anadoras species using a Leica DM 2000 microscope equipped with a digital camera Moticam Pro 282B. The best metaphases were photographed using the software Motic Images Advanced 3.2. Thereafter, we measured each chromosome using the Easyideo software, and the karyotype formulas were determined according to the arm ratio parameters initially proposed by Levan et al. [1964]. The heterochromatin pattern was detected in accordance with the protocol described by Sumner [1972] with modification in the staining [Lui et al., 2012].

We investigated the distribution of 3 repetitive DNA sequences: 18S rDNA, 5S rDNA, and telomere motifs (TTAGGG)_n. The 18S rDNA probe was obtained by Mini-Prep (i.e., extraction of plasmid DNA) from Prochilodus argenteus Spix and Agassiz, 1829 [Hatanaka and Galetti, 2004] and labeled by nick translation (Roche®) with biotin-16-dUTP according to the manufacturer’s instructions. The 5S rDNA probe was obtained by Mini-Prep from Megaleporinus elongatus Valenciennes, 1850 [Martins and Galetti, 1999] and labeled by nick translation (Roche) using digoxigenin-11-dUTP according to the manufacturer’s instructions. Oligonucleotide probe containing telomere sequences (TTAGGG)_n was amplified without DNA template as described by Ijdo et al. [1991] and labeled by PCR using biotin-16-dUTP (Roche).

FISH was performed under high stringency conditions on metaphase chromosome spreads of all Anadoras species, as detailed in Pinkel et al. [1986]. Briefly, 20 µL of the hybridization mixture (2.5 ng/μL probes, 2 μg/μL salmon sperm DNA, 50% deionized formamide, 10% dextran sulfate) was dropped on the slides, and the hybridization was performed for 24 h at 37°C in a moist chamber containing distilled water. The hybridization signals of 5S rDNA probes were detected using anti-digoxigenin-rhodamine (Roche), whilst the 18S rDNA and telomere probe signals were detected using avidin fluorescein isothiocyanate (Sigma). The fluorochrome 4′,6-diamidino-2-phenyl-indol (DAPI) was used as chromosome counterstain in a concentration of 1.2 μg/mL in an antifading solution (Vector). High-resolution FISH on stretched DNA (fiber-FISH) was performed according to the protocol described by Barros et al. [2011].

Genomic DNA was extracted from muscle, using the phenol-chloroform method described by Almeida et al. [2001]. A 600-bp fragment of the mitochondrial CO1 gene was amplified by PCR with the primers described by Ward et al. [2005]. Each PCR had a final volume of 10 µL containing 1 μL of DNA template (5 ng/mL), 0.2 μL of each primer (20 μM), 5 μL of master mix, and 3.6 μL of deionized water. The amplification cycles comprised 35 cycles of 1 min at 94°C, 30 s at 54°C, and 1.5 min at 72°C; a final extension of 5 min at 72°C and a cooling period at 4°C. PCR products were observed in a 1% agarose gel.

The PCR products were checked by electrophoresis in 1% agarose gel after staining with bromophenol blue and GelRed® (Biotium, USA) at a ratio of 3:1. Subsequently, the successfully amplified products were purified in 20% polyethylene glycol and washed in 80% ethanol to perform the bidirectional labeling reaction using the BigDye Terminator v3.1 Cycle Sequencing Ready Reaction kit (Applied Biosystems/Life Technologies) according to the manufacturer’s instructions. The PCR conditions consisted of 25 cycles at 96°C for 30 s, 15 s at 50°C, and 4 min at 60°C. After the reaction, the products were precipitated and sequenced in an ABI PRISM 3500 XL Genetic Analyzer (Applied Biosystems).

In this study, 11 sequences were generated and deposited on the Barcode of Life Data Systems (BOLD) platform in the project “Species delimitation of thorny catfishes (Doradidae-Astrodoradinae): an integrated analysis based on Chromosomal and DNA barcoding - DORBC.” It was automatically assigned to a barcode index number, BIN (group of sequences corresponding to a single taxon), following the analytical procedures of Ratnasingham and Hebert [2013]. In addition, we also included in the analysis sequences of Anadoras and other Astrodoradinae species deposited on BOLD, as well as Wertheimeria maculata (Wertheimerinae subfamily) and Platydoras armatulus (Doradinae subfamily) as outgroups (online suppl. Table S1; for all online suppl. material, see www.karger.com/doi/10.1159/000523747). The dataset with 604 bp was aligned on BioEdit Sequence Alignment Editor 7.1.9 [Hall, 1999] with Clustal W tool [Thompson et al., 1994]. In MEGA v.6 program [Tamura et al., 2011], the sequences were translated and intra-/interspecific genetic distance matrices were generated.

The neighbor-joining (NJ) tree was constructed based on the nucleotide substitution model Kimura-2-Parameters (K2P) [Kimura, 1980] and 1,000 bootstrap replicates. To verify the groupings between sequences and their reliability degree, phylogenetic reconstructions were carried out with Bayesian inference (BI) and maximum likelihood (ML) in the software Mr. Bayes 3.2.6 [Ronquist and Huelsenbeck, 2003] and RAxML-HPC BlackBox 8.2.10, respectively. For BI reconstruction, 2 independent runs with 4 Markov chains were utilized and 10,000,000 generations, sampling occurring every 1,000 generations, with 10% of burn-in; all other parameters were left at default. Nucleotide evolution model HKI + I was obtained in the jModelTest v2.1.6 [Darriba et al., 2012]. The analyses were implemented in the CIPRES Science Gateway 3.3 (http://www.phylo.org/index.php/portal/). The convergence of this analysis was verified in the Tracer v1.6 [Rambaut et al., 2019], with ESS values for each parameter exceeding 200. The posterior distribution trees were used to construct a maximum clade credibility tree in TreeAnnotator v1.8.4. The final tree was edited in FigTree v.1.4.2 [Rambaut and Drummond, 2016].

The barcodes obtained in the present study were submitted to the analysis of 4 algorithms based on a single locus commonly used for the identification of potential species: (1) Automatic Barcode Gap Discover – ABGD [Puillandre et al., 2012]; (2) Barcode Index Numbers – BIN [Ratnasingham and Hebert, 2013]; (3) Single General Mixed Yule-coalescent analyses – sGMYC [Pons et al., 2006]; and (4) Bayesian Poisson Tree Process – bPTP [Zhang et al., 2013]. The ABGD analyses were carried out with a matrix of distance pairwise, based on K2P, obtained in the MEGA v.6 [Tamura et al., 2011], on a free online platform (http://wwwabi.snv.jussieu.fr/public/abgd/). The BIN analysis was performed by comparing the taxa present in the input files with all others that share identical BINs, including those uploaded by different users based on the refined single linkage (RESL) algorithm. This system divides the CO1 sequences uploaded in the BOLD systems into molecular operational taxonomic units (MOTUs) independent of their predefined taxonomic classification. Therefore, the BINs represent a useful method to evaluate whether the DNA barcodes and species designation are in agreement or not [Ratnasingham and Hebert, 2013]. The sGMYC analysis was performed for mtDNA lineage delimitation [Fujisawa and Barraclough, 2013], avaliable on https://species.h-its.org/gmyc/, and using an ultrameric tree as input, obtained in BEAST 1.8.2 [Drummond et al., 2012]. Finally, as an input file for bPTP [Zhang et al., 2013], a maximum likelihood phylogenetic tree was constructed in the RAxML-HPC BlackBox 8.2.10, implemented in CIPRES. The bPTP adds to the branches of the input tree Bayesian inference support values, making the result even more reliable, being available on the free platform https://species.h-its.org/.

The diploid number found for all specimens was 56 chromosomes consisting of 24 metacentrics, 10 submetacentrics, and 22 subtelo-/acrocentrics (24m + 10sm + 22st/a) for both sexes (Fig. 2a). The C-banding technique revealed heterochromatic blocks in metacentric pair 2 (terminal position of the q and p arms), in submetacentric pairs 15 (p arm) and 17 (terminal q arm), and acrocentric pair 28 (terminal on the q arm) (Fig. 2b). FISH revealed 2 sites of 18S rDNA in subterminal position on the q arm of the acrocentric pair 28 and 2 sites of 5S rDNA on the p arm of the submetacentric pair 15 (Fig. 2c). The sample from Cuiabá River also showed 2n = 56, the same karyotype formula, and 5S rDNA sites on pair 15. The NOR sites also occurred on the long arm of pair 28, however in a terminal position (Fig. 2c). Regarding the C-banding pattern, the individuals collected in Cuiabá River had heterochromatin in pairs 15, 17, and 28, however subtle divergences were detected: (1) chromosome pair 2 did not exhibit heterochromatin blocks and (2) pair 28 showed an interstitial block, adjacent to the NOR sites (online suppl. Fig. S1).

This “putative species” presented 2n = 56 chromosomes and possesses the same karyotype formula as described for A. weddellii populations (Fig. 2d). C-banding revealed few heterochromatic blocks only in submetacentric pair 15 (in the p arm), subtelocentric pair 20 (in the p arm), and acrocentric pair 28 (terminal in the q arm) (Fig. 2e). The 18S rDNA, like in A. weddellii, was found to occupy the terminal position in the q arm of pair 28, and the 5S rDNA in the p arm of pair 15 (Fig. 2f).

All individuals had 2n = 56 chromosomes and a karyotype formula of 20m + 14sm + 22st/a for both sexes (Fig. 2g). Heterochromatic regions were identified in the terminal position of both p and q arms in the metacentric pairs 1 and 6 and submetacentric pair 13, terminally on the q arm of the metacentric pairs 2, 3, 5, and 9 and subtelocentric pairs 19 and 20, terminally on the p arm of meta- and subtelocentric pairs 10 and 22, as well as pericentromeric blocks in submetacentric pairs 11 and 15 (Fig. 2h). FISH with rDNA probes revealed the coexistence of 2 arrangements named pattern A and B. Pattern A was observed in 5 individuals which showed three 18S rDNA sites – 2 in terminal position of the p arms of pair 13, syntenic with interstitial 5S rDNA sites, and the other in terminal position on the q arm of only one chromosome 28. Pattern B was described for 2 individuals that had only syntenic rDNA sites on the p arm of pair 13 (Fig. 2i). Among the individuals with pattern A, in the majority of the metaphases and nuclei, it was possible to detect an association between the chromosomes bearing the 18S rDNA sites, and often, one of the chromosomes of pair 13 appeared associated with the chromosome pair 28 (online suppl. Fig. S2). Fiber-FISH with rDNA probes confirmed the adjacent organization of 18S and 5S rDNA sites on the p arm of pair 13 and revealed a short interstitial region where both rDNA multigenic families are interspaced (online suppl. Fig. S2).

The physical mapping by FISH using telomeric probes, evidenced in all Anadoras species here analyzed only terminal sites in both arms of all chromosomes, i.e., no ectopic telomeric sequences were observed (online suppl. Fig. S3).

The sequencing of the CO1 gene from 11 Anadoras specimens generated a dataset composed of 604 bp where no deletions, insertions or stop codons were found, which means that this region corresponds to a functional portion, useful for DNA barcoding. The NJ tree as well as the phylogenetic reconstructions based on ML and BI methods produced a similar topology evidencing 6 major clades with high support values (above 90% bootstrap and 0.9 posterior probability). The major clades are: Wertheimeria and Platydoras (outgroups), Scorpiodoras, Hypodoras + Astrodoras, Physopyxis, and Anadoras. The last one is subdivided into 3 well-supported clades named Anadoras sp. “araguaia”, A. weddellii, and A. grypus (Fig. 3). Among the Anadoras species, the interspecific distance values ranged from 5.5% (Anadoras sp. “araguaia” and A. weddellii), 6.3% (A. weddellii and A. grypus) to 6.9% (Anadoras sp. “araguaia” and A. grypus). The intraspecific genetic distance values were 0.36% among populations of A. weddellii from Miranda, Cuiabá rivers, and Lower Amazon basin, while for A. grypus from Catalão Lake (single NORs + multiple NORs), Upper Amazon, and Negro River the intraspecific genetic distance was 0.56% (online suppl. Table S2).

The species delimitation algorithms show different results regarding the number of MOTUs detected. While the ABGD and BIN identified 20 MOTUs, bPTP detected 21 MOTUs, and GYMC recognized 22 MOTUs. Despite this numeric difference, it is important to highlight that all algorithms recognized in Anadoras 3 MOTUs and confirm that Anadoras sp. “araguaia” is a valid species (Fig. 3).

The earliest chromosome studies performed in species of Doradidae showed little karyotype variation: most species have 2n = 58 chromosomes, few heterochromatin blocks, and simple terminal NOR sites. Based on this scenario of karyotype homogeneity, Eler et al. [2007] and Milhomem et al. [2008] proposed a model of chromosome diversification mediated by pericentric inversion events, due to the constancy of 2n = 58 and great variation in karyotype formulas. Recently, this scenario changed drastically, revealing numeric and structural variations in several thorny catfishes, besides supernumerary chromosomes and a unique ZW sexual system [Takagui et al., 2017, 2021; Baumgärtner et al., 2018]. According to Takagui et al. [2021], the ancestral diploid number for the Doradidae family still remains as an open question, but several lines of evidence support the 2n = 58 chromosomes as a plesiomorphic feature in the family and as the likely ancestral karyotype configuration. Furthermore, this diploid number also occurs in most species of Auchenipteridae [Lui et al., 2009, 2010, 2013a, 2013b, 2015, 2021; Felicetti et al., 2021; Machado et al., 2021; Santos et al., 2021], the sister group to Doradidae family [Arce et al., 2013; Birindelli, 2014; Sabaj and Arce, 2021].

Assuming that the 2n = 58 chromosomes indeed represent the ancestral diploid number for thorny catfishes, our results show that the 2n = 56 described in all Anadoras species would be a derived condition, likely attributed to fusion events, as already reported in Trachydoras paraguayensis [Baumgärtner et al., 2016, 2018]. We believe that this chromosome reduction does not represent the ancestral condition for these 2 taxa, and it is likely to have arisen from independent events. Such a hypothesis can be supported by at least 3 lines of evidence: (1) Anadoras spp. and T. paraguayensis are species that are located in different subfamilies in Doradidae, more related to species with 2n = 58 than to each other; (2) the karyotypes are distinct; for instance, in Anadoras spp. several subtelo-/acrocentric pairs are reported, which do not occur in T. paraguayensis, whose karyotype is composed only of bi-armed chromosomes; (3) the probable fusion event that led to a reduction in the chromosome number in Anadoras and Trachydoras spp. did not involve the same chromosome pairs, likely to have arisen from different mechanisms. Interestingly, in T. paraguayensis, interstitial telomeric repeats were detected in the largest submetacentric chromosome pair, suggesting that this chromosome pair was involved in centric fusion events. On the other hand, in Anadoras spp. no species exhibited traces of interstitial telomeric sequences, supporting that the fusion event either eliminated the telomeres completely or changed/reduced the number of these repeats to a point that could not be detected by FISH.

The physical mapping of the rDNA probes revealed simple systems for both 18S rDNA and 5S rDNA sites, located in distinct chromosomes in A. weddellii and Anadoras sp. “araguaia,” which reinforces their phylogenetic proximity at a cytogenetic point of view. It is noteworthy that, even though the 2 populations of A. weddellii present the NOR sites in the same chromosome pair (pair 28), these sites are arranged in different positions. This subtle divergence probably occurred due to a paracentric inversion that changed the location of the NORs (subterminal/terminal) and part of the adjacent heterochromatic block (online suppl. Fig. S1B). On the other hand, both 5S and 18S rDNA are organized in a quite specific pattern in A. grypus, in which the majority of specimens from the Catalão Lake had three 18S rDNA sites, 2 of them in synteny with the 5S rDNA on the p arm of pair 13, and the other one on the q arm of only one acrocentric chromosome 28. However, 2 specimens had only the syntenic pair 13 carrying these sequences (online suppl. Fig. S2A). Fiber-FISH revealed that the rDNA sites are adjacent to each other, with the 18S rDNA sites in terminal position and the 5S rDNA sites near the centromeres, and both are interspersed in a small intermediate region (online suppl. Fig. S2B).

The existence of individuals of the same species harboring different NOR patterns in sympatry is usually related to different chromosomal rearrangements, mainly pericentric inversions and translocations [Mariotto et al., 2006; Porto et al., 2011; Konerat et al., 2015] and even with activity of moving transposable elements [Porto et al., 2014a, b]. We believe that the origin of numerical polymorphism in the distribution of the rDNA sites in A. grypus occurs as a remnant of nonreciprocal translocations between terminal segments of the pairs 13 and 28. During the interphase, the chromosomes are less condensed and occupy specific nuclear domains due to the existent interactions between telomeres and the nuclear envelope. This conformation, named Rabl’s model, guarantees a closer arrangement among nonhomologous chromosomes and may promote translocations of terminal segments [Cremer et al., 1982; Schweizer and Loidl, 1987; Cremer and Cremer, 2010]. A strong evidence supporting this hypothesis is the recalcitrant association of one of the chromosomes of pair 28 to the pair bearing the syntenic sites of rDNA; such association was also confirmed by the analysis of interphase nuclei submitted to double FISH (online suppl. Fig. S2A).

Our cytogenetic approach revealed great similarity between the karyotypes of A. weddellii and Anadoras sp. “araguaia”; both have 2n = 56 chromosomes and the same karyotype formula, as well as the rDNA sites on the same chromosome pairs. The sole difference between these 2 species was detected in the amount of heterochromatin, with A. weddellii presenting additional heterochromatic blocks to those detected in Anadoras sp. “araguaia.” A. grypus has a great level of chromosomal difference when compared to its congeners: different karyotype formula, distribution pattern of the rDNA sites, and mainly in the amount of heterochromatin. Thus, we believe that the karyotype diversification of Anadoras occurred through the following processes: (1) initially from a centric fusion, which reduced the diploid number of 58 chromosomes to 2n = 56; (2) pericentric inversions and nonreciprocal translocations, which promoted shifts in karyotype formulas and the rise of the multiple NOR system in A. grypus; (3) intense reorganization of heterochromatic segments, perhaps as a result of the jumping activity of certain transposable elements; and finally (4) paracentric inversions, which are rearrangements hardly detected by conventional cytogenetic analysis, but they modified the position of the NOR sites in the long arm of pair 28 of A. weddellii populations.

The genomic diversity that results from the interaction among certain repetitive sequences and chromosomal rearrangements may become a driver of diversification in several organisms [Rieseberg, 2001; Raskina et al., 2008; Faria and Navarro, 2010; Fernandes et al., 2017]. Their accumulation in specific genomic regions, often referred to as “fragile sites,” may induce chromosome breakages, deletions, inversions, and amplifications [Eichler and Sankoff, 2003; Pevzner and Tesler, 2003; Barros et al., 2017]. Several studies demonstrate that pericentric inversions can promote genetic incompatibilities among species, an important mechanism of postzygotic isolation [Noor et al., 2001; Kirkpatrick and Barton, 2006]. Pericentric inversions have been frequently detected in many groups of fish, especially among the Siluriformes where they are considered the major rearrangements in the karyotype diversification in Loricariichthys [Takagui et al., 2014], Rineloricaria [Venturelli et al., 2021], Trachelyopterus [Santos et al., 2021], Tatia [Lui et al., 2013a], Pimelodus [Girardi et al., 2018], and Trichomycterus [Oliveira et al., 2016]. Among Doradidae, multiple pericentric inversions events associated with heterochromatin movements occurred independently in several clades and played a key role in the diversification of Wertheimerinae [Eler et al., 2007; Takagui et al., 2019] and Doradinae subfamilies, especially in the large species allocated in the nonfimbriated barbels clade [Milhomem et al., 2008; Baumgärtner et al., 2018; Takagui et al., 2021] and in the widespread species T. paraguayensis [Baumgärtner et al., 2016].

In the light of traditional taxonomy, the 2 Anadoras species formally described (A. grypus and A. weddellii) can be easily distinguished by external and internal morphological features, however, Anadoras sp. “araguaia” has fewer diagnostic characters identified for its validation. Moreover, this putative new species is morphologically distinct from A. grypus, but shares many traits with A. weddellii, which differs only with respect to the number of gill traces and in few osteological characters [Sousa, 2010]. As previously mentioned, cytogenetic data also point to a greater similarity between Anadoras sp. “araguaia” and A. weddellii, which together with little morphological differentiation and lack of data on the CO1 gene, led us to question whether the species described for the Araguaia River is really a valid taxonomic unit or a population of A. weddellii with morphological and chromosomal variations. To elucidate this question, we performed species delimitation analysis based on sequencing of fragments (600 bp) of the CO1 gene. This methodology is based on the fact that, in order to discriminate apparent similar species, intraspecific genetic variation must be significantly lower than interspecific variation [Hebert et al., 2004]. This lack of overlap between intra- and interspecific variation creates the barcode gap, leading to the distinction of the taxa. In fact, intra- and interspecific average distances within Anadoras were 0.45% and 5.9% and demonstrate a clear gap to an interspecific differentiation 13 times greater than the intraspecific variation. Therefore, DNA barcode analysis confirms indeed that Anadoras sp. “araguaia” represents a genetically distinct species, as corroborated by the phylogenetic reconstructions (NJ, BI, ML) and delimitation algorithms (ABGD, BIN, bPTP, and GYMC).

Considering the karyotype data, sequences [Arce et al., 2013], and morphology [Sousa, 2010], we believe that the morphological similarity between A. weddellii and Anadoras sp. “araguaia” is a classic example of adaptive convergence, as already reported in Hemiodontichthys [Carvalho et al., 2018], Rineloricaria [Costa-Silva et al., 2015], and Pimelodella species [Souza-Shibatta et al., 2013]. According to Sousa [2010], Anadoras species rarely occur in the main channels of large rivers and prefer to inhabit marginal lakes associated withfloating or riparian vegetation. Thus, we believe that A. weddellii and Anadoras sp. “araguaia” still evolved convergent morphologies due to the selective pressure that the lacustrine environment exerts mainly on the body shape, color pattern, swim bladder, and on some fin structures. In fact, the DNA barcoding system is based on the premise that the distinction of different species is possible because interspecific genetic variation is higher than intraspecific variation [Hebert et al., 2003], like already established in several animal groups, including fish [Barrett and Hebert, 2005; Ward et al., 2005]. However, this short region of the genome evolves fast enough to provide resolution at species level and assist their delineation. In caveats, it is noteworthy that the use of only one kind of molecular marker or gene may not reflect the real evolutionary history of the organisms [Edwards and Beerli, 2000]. Thus, the use of different markers can result in different genealogies and allows the reduction of the effects of stochasticity [Brito and Edwards, 2009]. Therefore, combining mitochondrial and nuclear markers is recommended in order to perform a more accurate confirmation and provide an independent estimate of the evolutionary history of Anadoras species.

In the present study, we recovered several cytogenetic synapomorphies for Anadoras. These data are relevant for better understanding the evolutionary history of the thorny catfishes, especially in Astrodoradinae, since Anadoras is the clade with the oldest divergence and sister group of the other genera included in this subfamily. The karyo-evolution process of Anadoras seems to be related to structural rearrangements and an intense dispersion of repetitive sequences, which explain the origin of the high distinctive karyotype reported by the Amazonian species A. grypus. The other species, A. weddellii and Anadoras sp. “araguaia” share many karyotypic features and show subtle differences only in the C-banding pattern. Our species delimitation analysis based on the 600-bp fragment of the CO1 gene confirmed previous morphological studies and identified in Anadoras 3 MOTUs, and recognized Anadoras sp. “araguaia” as a valid species with a high genetic divergence value in comparison to the others. Thus, our data encourage the description of this new species and reaffirm the importance of cytogenetic and DNA barcoding for biodiversity assessment, mainly in endemic regions threatened by anthropogenic activities, such as the Upper Araguaia River.

The authors are grateful to Jansen Zuanon (INPA) for the assistance in collecting some of the studied fishes. We thank the Instituto Nacional de Pesquisas da Amazônia (INPA) and Universidade Estadual de Londrina (UEL), Centro de Ciências Biológicas (CCB), Departamento de Biologia Geral, for providing the laboratory structure; Coordenação de Aperfeiçoamento de Pessoal de Nível Superior CAPES (Finance code 001) and Fundação Araucária for their financial support through a doctoral grant to Fabio Hiroshi Takagui; Conselho Nacional de Desenvolvimento Científico e Tecnológico CNPq for their support through productivity grant to Lucia Giuliano-Caetano (process 302872/2018-3); and Instituto Chico Mendes de Conservação da Biodiversidade (ICMBio) for permitting the collection of biological material.

The experiments followed ethical conducts according to the Ethics Committee for Animal Use of the Universidade Estadual de Londrina, under the protocol number 60/2017.

The authors have no conflicts of interest to declare.

Coordenação de Aperfeiçoamento de Pessoal de Nível Superior CAPES (Finance code 001) and Fundação Araucária provided financial support through a doctoral grant to Fabio Hiroshi Takagui. Conselho Nacional de Desenvolvimento Científico e Tecnológico CNPq supported through a productivity grant to Lucia Giuliano-Caetano (process 302872/2018-3). The expeditions for fish collection were supported by resources provided by the Postgraduate Program in Genetics and Molecular Biology of Londrina State University and Instituto Nacional de Pesquisas da Amazônia (INPA).

Fabio Hiroshi Takagui conducted all cytogenetic and molecular experiments, collected the species, analyzed the data, conceived and designed the study, and wrote the manuscript. Paulo Cesar Venere, Lucia Giuliano-Caetano, Roberto Laridondo Lui, Patrik Viana, and Orlando Moreira-Filho collected the specimens. José Luis Olivan Birindelli helped to identify the specimens. Lucas Baumgärtner and Patrik Viana assisted in obtaining chromosomal preparations, in fluorescence in situ hybridizations, and wrote the manuscript. Jamille de Araujo Bitencourt and Moema Cristina Costa Lima assisted in the DNA barcoding experiments and analysis. Roberto Laridondo Lui, Fernanda Almeida Simões, and Lucia Giuliano-Caetano provided laboratorial structure for some cytogenetic analyses, helped in designing the study, and wrote the manuscript. All authors read and approved the final version.

All data generated or analyzed during this study are included in this article and its online supplementary material. Further inquiries can be directed to the corresponding author.