Towards a complete classification of the Neotropical thorny catfishes (Siluriformes: Doradidae)

Mark Henry Sabaj¹ and Mariangeles Arce H.¹

PDF: EN    XML: EN  | Supplementary:  S1  S2  S3 | Cite this article

Abstract​

---

Introduction​

---

Thorny catfishes (Siluriformes: Doradidae) form a monophyletic group of about 96 valid extant and one fossil species endemic to freshwaters of South America on both sides of the Andes Mountains. Most doradids are easily distinguished from other catfishes by having a conspicuous midlateral row of bony scutes, each one with a central, caudally directed thorn (Fig. 1). Each midlateral scute is formed by dorsal and ventral aliform expansions of a lateral-line tubule. A single enlarged pore perforates the skin in the axil of each thorn. The infranuchal scute is exceptionally composed of both an expanded tubule and an ossified ligament that runs between the nuchal region of the cranium and the rib supported by the sixth vertebra, which is the first long-formed rib. As such, the infranuchal scute represents an unambiguous synapomorphy for Doradidae (Birindelli, 2014). Another synapomorphy for doradids is the presence of Sörensen’s ligament (Fig. 2), an unossified ligament between the anterolateral rim of the Müllerian disk and an ossified tubule or scute in the tympanic region (Birindelli, 2014).

FIGURE 1 | Variation in scute morphology in cleared and stained specimens of Doradidae. A. Amblydoras nheco (ANSP 187416); B. Megalodoras uranoscopus (ANSP 178302); C. Hassar orestis (ANSP 181094); D. Leptodoras linnelli (ANSP 182791). Infranuchal scute (is), exceptionally composed of expanded lateral-line tubule and ossified ligament between nuchal region of skull and rib supported by 6th vertebra.

Adult thorny catfishes vary in standard length from about 22 mm (Physopyxis ananas Sousa & Rapp Py-Daniel, 2005) to over one meter (Oxydoras spp.). Doradids generally occupy benthic habitats in lowland lakes and rivers, although a few taxa frequent pelagic habitats, such as Nemadoras hemipeltis (Eigenmann, 1925) and Pterodoras Bleeker, 1862. Many of the smaller species are peculiar to floodplains and occupy sluggish streams and river margins during the low-water season. The larger species are restricted to the main channels of medium to large rivers. A few doradids (e.g., Rhinodoras Bleeker, 1862) often associate with large rocky rapids in rivers draining the Brazilian and Guiana shields. The propensity of thorny catfishes for large river channels and lowland floodplains, coupled with an absence from upland headwaters, makes Doradidae a prime candidate for investigating large scale shifts in Neotropical drainage patterns.

FIGURE 2 | Sörensen’s ligament (sl), unossified ligament between anterolateral rim of Müllerian disk (md) and first ossified tubule or scute (not visible) in tympanic region; stained specimen of Oxydoras sifontesi (ANSP 181069, 149.5 mm SL). gb = gas bladder, is = infranuchal scute, pcp = posterior cleithral process, pnp = posterior nuchal plate, ptsc = posttemporal-supracleithrum.

The taxonomic history of thorny catfishes includes 43 nominal genera and 146 nominal species (Fig. 3) dating back to the Linnaean (Linnaeus,1758) descriptions of Acanthodoras cataphractus and Platydoras costatus (Sabaj Pérez, 2014). Lacepède (1803) proposed the first genus (Doras), Bleeker (1858) established the family-group name Doradidae (van der Laan et al., 2014), and Higuchi et al. (2007) described the first valid subfamily, Astrodoradinae. Kner (1853, 1855) published the first detailed descriptions of doradids in his treatment of 18 species including 13 proposed as new. Although Kner recognized only one genus (Doras), his species spanned 14 of the 31 genera considered valid in the family. Eigenmann (1925) compiled a comprehensive monograph on Doradidae that is rich with figures and acute observations that continue to inform modern studies. Sabaj, Ferraris (2003) assembled an annotated checklist of doradids that clarified or highlighted a number of nomenclatural and taxonomic issues. The next fifteen years witnessed the descriptions of two new genera and 23 new species, nearly a quarter of the total species considered valid here. Although the classification of Doradidae is more or less complete to the genus level, taxonomic work remains at the species level for a number of genera, especially Acanthodoras Bleeker, 1862, Amblydoras Bleeker, 1862, Anadoras Eigenmann, 1925, Hemidoras Bleeker, 1858, Platydoras Bleeker, 1862, and Pterodoras Bleeker, 1862.

FIGURE 3 | Summary of taxonomic history of Doradidae. Each square represents a nominal valid taxon (black) or putative synonym (red) plotted against the year of its description (two nomina oblita and four replacement names not included). Continuous black line traces ratio of valid to nominal species through four time periods: burn in, naming, lapsus and modern revision. Monographs by Kner (1855) and Eigenmann (1925) mark ends of burn in and naming periods, respectively. Monograph by Sabaj, Ferraris (2003) marks beginning of modern revision.

Cladistic studies of doradids began with Higuchi (1992) who used morphology to hypothesize relationships within the family inclusive of a previously contentious member, Wertheimeria maculata Steindachner, 1877. Arce H. et al. (2013) provided robust support for alternative relationships based on phylogenetic analyses of molecular data. Birindelli (2014) assembled the most comprehensive morphological data set to date to investigate phylogenetic relationships among Doradidae and its sister family Auchenipteridae. Based on those results, Birindelli (2014) firmly diagnosed Doradidae and proposed a new subfamily, Wertheimerinae. Other recent studies have described variation in gas bladder morphology (Birindelli et al., 2009), sperm morphology (Quagio-Grassiotto et al., 2011), bioacoustics (Kaatz, Stewart, 2012; Zebedin, Ladich, 2013; Knight, Ladich, 2014), digestive tube morphology (de Melo Germano et al., 2014) musculature (Arce H., 2015) and cytogenetics (Baumgärtner et al., 2018; Takagui et al., 2019). Drawing heavily from variation in caudal-fin morphology, Birindelli, Sousa (2018) assembled a key to the 26 doradid genera inhabiting the Amazon, Orinoco and Guianas.

The primary goals of this study are to advance the classification and summarize the geographic distributions of thorny catfishes. We expanded the taxon sampling of the molecular data set analyzed by Arce H. et al. (2013) and compiled comprehensive data on the caudal skeleton for all doradid taxa. Based on our analyses of those data, we propose a revised classification of Doradidae and comment on morphological trends observed in caudal-fin evolution among doradids and other catfishes.

Material and methods

---

Molecular Data: markers and taxon sampling. Sequence data were assembled for one nuclear gene, recombination activating gene 1 (rag1), and two mitochondrial genes, cytochrome c oxidase subunit 1 (co1) and 16s ribosomal RNA (16s), from 218 specimens representing 37 outgroup taxa (44 specimens) and 100 ingroup taxa (174 specimens) (Tab. 1). The current analysis employed the same three markers used by Arce H. et al. (2013), but added 74 specimens (43 doradids and 31 outgroups) and 38 species-level taxa (14 doradids and 24 outgroups).

Outgroup taxa were selected on the basis of molecular studies (Sullivan et al., 2006; Lundberg et al., 2007; Nakatani et al., 2011; Arcila et al., 2017; Betancur-R. et al., 2017; Calegari et al., 2019) that support Diplomystidae as the sister group to Siluroidei, Cetopsidae as the sister group to all other siluroids, and Aspredinidae as sister to Auchenipteridae + Doradidae, with those three families comprising the superfamily Aspredinoidea Adams, 1854 (van der Laan, 2019:121; see also Results). The ingroup taxa represented all 31 valid genera of Doradidae, 83 of the 96 extant valid species, and 17 taxa that are undescribed species or currently unassigned to nominal ones.

TABLE 1 | List of taxa, voucher specimens and DNA sequences analyzed. *Denotes individuals sequenced in Arce H. et al. (2013). Museum codes follow Sabaj (2020). ^a Sequence data published by Sullivan et al. (2006) for voucher ANSP 180476. ^b Sequence data submitted to GenBank by Heok Hee Ng (2006) for voucher ANSP 180476 (tag 4515) from an unpublished study. ^c Sequence data published by Nakatani et al. (2011); no voucher data. ^d Genus assignment based on Calegari et al. (2019). ^e Questionably a junior synonym of Hemidoras boulengeri (Steindachner, 1915).

(click here to Table 1)

Molecular Data: DNA extraction, amplification and sequencing. Generally, tissues (e.g., fin, muscle or gill) were taken in the field and preserved in 95–100% ethanol; voucher specimens were fixed in 10% buffered formalin, then transferred to 70–75% ethanol for long-term museum storage. Ideally, the tissue sample is associated with a field tag number that is tied to the voucher specimen. Additional tissue samples were provided by generous colleagues (see Acknowledgments).

Total DNA was extracted using the Qiagen DNeasy blood and tissue kit. PCR was carried out in 20 µl reactions; primers for amplification and sequencing are listed in Arce H. et al. (2013:561, tab. 1). For co1 and 16s, the PCR reaction mixture consisted of 10 µl of Apex Taq DNA Polymerase Master Mix, 1.5 mM MgCl2 (Genesee Scientific), 0.5 µM of forward and reverse primer, 5–8 µl of distilled water and 1–4 µl of DNA template. Cycles of amplification were programmed accordingly: 95°C for 4 min (initial denaturation), 10 cycles of three steps, 50°C or 55°C for 30 sec (annealing, temperature decreased by 1°C after each cycle), 72°C for 2 min (extension) and 95°C for 1 min (denaturation); 30 cycles of three steps, 95°C for 1 min, 40°C or 44°C for 30 sec, and 72°C for 2 min; final extension step at 72°C for 10 min. Amplification of rag1 followed the protocol of Sullivan et al. (2006): 4 min at 95°C (initial denaturation), 35 cycles of three steps, 30 sec at either 50°C, 55°C or 59°C, 2 min at 72°C, and 30 sec at 95°C; final extension step for 4 min at 72°C. Amplifications were sent to Functional Biosciences, Inc. laboratories for purification and sequencing.

Molecular Data: sequence alignment and phylogenetic analyses. Sequences were edited and combined into contigs for each marker (rag1, co1, 16s) in Geneious 11.1.2 (Drummond et al., 2010). Complete gene sequences were aligned in MUSCLE 3.7 (Edgar, 2004) using default parameters. Alignments were refined manually, and sequences for the three markers were concatenated in Mesquite 3.40 (Maddison, Maddison, 2011). Translations of new sequences for co1 and rag1 were aligned in COBALT (Papadopoulos, Agarwala, 2007) to correct for frameshifts and to trim low-quality ends prior to DNA sequence alignment.

We analyzed combined nuclear and mitochondrial sequences using Maximum Parsimony (MP) and Bayesian Inference (BI), and employed the same parameters as Arce H. et al. (2013) for comparability. Analyses were performed on the combined dataset with terminals restricted to those represented by at least two loci (i.e., 218 specimens; Tab. 1). For MP analysis, the trees were generated using the “new technologies search” implemented in TNT (Goloboff et al., 2008) and performed in two steps. The first step used a combination of sectorial searches (RSS and CSS), 100 iterations of ratchet, 100 cycles of tree fusing, and 100 rounds of drift; driven was set to reach the minimum length 50 times. The second step used the trees produced in the first search to perform a traditional TBR search. Gaps were treated as missing data and all characters had equal weights. Godman-Bremer support (Goodman et al., 1982; Bremer, 1988, 1994; Grant, Kluge, 2008) was calculated for each node and plotted on the consensus tree.

For Bayesian analyses, the concatenated gene matrix was divided into eight partitions: one for 16s, one for each nucleotide position per co1 codon, one for each nucleotide position per rag1 codon, and one for the rag1 intron. Bayesian analyses were conducted in MrBayes 3.1.6 (Huelsenbeck, Ronquist, 2001; Ronquist, Huelsenbeck, 2003) using the GTR + GAMMA model. We ran three heated chains and one cold chain for 60 million generations, sampling every 10,000th generation. To ensure sampling of the posterior distribution we discarded 0.25% of the trees.

Morphological Data. Specimens examined for morphological data were designated as alc (alcohol), sk (dry skeleton) or cs (cleared and stained following the methods of Taylor, Van Dyke, 1985). Data on the caudal skeleton were taken from cleared and stained specimens and dry skeletons while immersed in 90% glycerin and 75% ethanol, respectively, and viewed under a Wild M3C stereomicroscope. Immersion facilitated the removal of residual muscle tissue and assessment of sutures. Midlateral scutes were removed from both sides to facilitate clear observations of the caudal skeleton. Observations were made on adult specimens and juveniles at stages where the caudal skeleton was already mostly ossified. In a few cases, ontogeny was used to hypothesize fusion between elements (e.g., procurrent caudal-fin rays in some astrodoradins). But for the most part, fusions between elements (e.g., hypurals, parhypural) was presumed and not directly observed via ontogenetic series (e.g., Vaz, Hilton, 2020). Museum codes follow Sabaj (2020).

For descriptions of the caudal skeleton, we employed the diural scheme which considers the last vertebra to be a compound caudal centrum formed by the fusion of the posteriormost preural centrum (PU1) plus anteriormost ural centrum (U1) (Lundberg, Baskin, 1969; Grande, Shardo, 2002; de Pinna, Ng, 2004; Bird, Mabee, 2003; Bensimon-Brito et al., 2012). In cases where a second ural centrum (U2) is visible, it is sometimes considered a fusion product of two or three originally distinct centra (Arratia, 2003; de Pinna, Ng, 2004; Bensimon-Brito et al., 2010, 2012). The compound caudal centrum (PU1+U1) supports the pleurostyle (PL), hypurals (HY) and parhypural (PH). We use the generic term pleurostyle for the elongate process that projects at an angle from the dorsal posterior corner of compound caudal centrum. Previous authors used the term uroneural (i.e., modified ural neural arch) for this process in catfishes (e.g., Lundberg, Baskin, 1969; Grande, Shardo, 2002; de Pinna, Ng, 2004); however, the homology and evolution of this process remains uncertain among ostariophysans (Cumplido et al., 2020). Hypurals are ventral bony elements separated into lower hypurals (HY1,2) and upper hypurals (HY3,4,5,6 in catfishes) by a diastema or gap for the passage of paired arterial and venous branches leading to and from the caudal fin (Desvignes et al., 2018). The parhypural represents the last haemal arch and spine, and the hypurals are considered modified haemal spines of the ural centra (Arratia, Schultze, 1992; Schultze, Arratia, 2013).

Lundberg, Baskin (1969) introduced a formula for describing various patterns of fusion and/or loss among the elements supported by the compound caudal centrum (PU1+U1). They used a plus sign (+) between adjacent elements that are presumably completely fused (e.g., PH+HY1+2), and a semicolon (;) between adjacent elements that remain separated or at least distinguishable, often by a long and continuous suture (e.g., PH; HY1; 2). Although the parhypural and ventral hypurals may appear separate and scored as such, these three elements are tightly associated or fused (continuous) proximally near their fusion to the compound caudal centrum from early developmental stages to adulthood in catfishes (Grande, Shardo, 2002; Adriaens, Vandewalle, 2003). When the sixth hypural was not distinguishable, it was presumed lost rather than fused, and thereby omitted from the formula.

For scoring individuals, we modified the formula of Lundberg, Baskin (1969) by using a hyphen (-) between elements that are only partially fused and retain features suggestive of independence such as distal or internal gaps and/or semitransparent windows of thin bone; figures in Grande, Shardo (2002) similarly employed hyphens. For scoring a taxon as a whole, a hyphen in the formula also might represent polymorphism where two elements may appear completely fused in some specimens, but separate in others. For completeness, we also included the pleurostyle (PL) and epural (EP) in the formula because those elements may fuse with each other or with the upper hypural plate in some taxa. Principal caudal-fin rays are reported as branched (Arabic numeral) or simple (lower case Roman numeral).

Character state mapping. For two characters associated with fusion patterns in the caudal skeleton, states were mapped on the Maximum Parsimony phylogeny generated in the current study for Aspredinidae, Auchenipteridae and Doradidae (i.e., Aspredinoidea). The first character was divided into two states: parahypural separate (1) or fused (2) with hypurals 1+2. The second character involved the upper hypurals (HY) and pleurostyle (PL) and exhibited three states treated as ordered: HY3+4; 5; PL (1), HY3+4+5; PL (2), and HY3+4+5-PL or HY3+4+5+PL (3). Next, each possible state was assigned to the common ancestor of the three families. Then, the number of transformations necessary to achieve the phylogenetic distribution of states in the terminal lineages was assessed by eye. The inferences from this exercise are presented in the Discussion.

Results​

---

Molecular Analyses. In our final analyses, 180 of the 218 specimens were represented by complete molecular data (all genes: rag1, co1, 16s; Tab. 1). Seven specimens were represented only by rag1 and co1 sequences, nine specimens were represented only by rag1 and 16s sequences, and 22 specimens were represented only by co1 and 16s sequences. The Maximum Parsimony (MP) analysis produced 144 most parsimonious trees of 9235 steps each. Under MP, the rag1 dataset consisted of 1861 total and 716 parsimony-informative base pairs for 196 specimens, the 16s dataset consisted of 583 total and 188 parsimony-informative base pairs for 211 specimens, and the co1 dataset consisted of 593 total and 246 parsimony-informative base pairs for 209 specimens. The combined dataset included 3037 total base pairs of which 1150 were parsimony informative for 218 terminals.

Trees produced by the Maximum Parsimony (MP) and Bayesian (BI) analyses were highly resolved and agreed on most intergeneric relationships (Figs. 4, S1, S2, S3) with a few notable exceptions. The largest disagreement between the MP and BI topologies involved the base of Doradidae. In the MP analysis, Acanthodoradinae was the first subfamily to diverge from the rest of Doradidae and Astrodoradinae was the second. BI reversed this topology with Astrodoradinae diverging first, followed by Acanthodoradinae. Relationships within Astrodoradinae also differed between the two analyses. Both identified Anadoras Eigenmann, 1925 as the first genus to diverge in Astrodoradinae. MP supported Physopyxis Cope, 1871 sister to Astrodoras + Hypodoras and Amblydoras Bleeker, 1862 sister to Scorpiodoras Eigenmann, 1925. BI placed Physopyxis sister to a clade composed of Scorpiodoras and Amblydoras (Astrodoras + Hypodoras).

FIGURE 4 | Phylogenetic relationships among all genera and subfamilies of Doradidae inferred from Maximum Parsimony analysis of rag1, 16sand co1DNA sequence data (strict consensus of 144 most parsimonious trees, each with 9235 steps).

Within the subfamily Doradinae, MP and BI differed in four major respects. In the parsimony analysis, Doraops + Pterodoras was the first group to diverge within Doradinae, followed by Oxydoras Kner, 1855. BI weakly supported (0.5 posterior probability) the reverse with Oxydoras as the first genus to split from the rest of Doradinae, followed by Doraops + Pterodoras. A second difference between MP and BI was placement of the clade Centrodoras (Lithodoras + Megalodoras). In the parsimony analysis, Centrodoras (Lithodoras + Megalodoras) was sister to the fimbriate-barbel doradids. Alternatively, BI supported a sister group relationship between Centrodoras (Lithodoras + Megalodoras) and Centrochir + Platydoras, and that clade was sister to the fimbriate-barbel doradids. Thirdly, MP supported the monophyly of Doras inclusive of Doras punctatus Kner, 1855 a species formerly assigned to Ossancora (Birindelli, Sabaj Pérez, 2011), and placed Doras sister to all other fimbriate-barbel doradids. In the BI analysis, Doras carinatus (Linnaeus, 1766; type species), D. micropoeus (Eigenmann, 1912), and D. higuchii Sabaj Pérez & Birindelli, 2008 formed a clade sister to all other fimbriate barbel taxa except D. phlyzakion Sabaj Pérez & Birindelli, 2008 and D. punctatus. Those two species, respectively, were successive sister taxa to the remaining fimbriate-barbel taxa. Finally, near the crown of the doradid tree, MP and BI disagreed on relationships within a clade composed of Hassar Eigenmann & Eigenmann, 1888, Nemadoras Eigenmann, 1925, Tennellus Birindelli, 2014 and Hemidoras + Ossancora. MP weakly supported two monophyletic clades, Nemadoras + Tennellus and Hassar (Hemidoras + Ossancora), each with a Godman-Bremer support value of 1 (Fig. S2). In the BI analysis, Nemadoras was the first genus to diverge and Tennellus + Hassar and Hemidoras + Ossancora formed reciprocally monophyletic clades (Fig. S1).

Our revised classification of Doradidae (Tab. 2; Fig. 4) is based on relationships supported by the Maximum Parsimony analysis of the DNA sequence data. The results of the Bayesian analysis are consistent with our classification except for the monophyly of Doras which is supported only by MP. Except for Astrodoradinae and Doradinae, each subfamily is diagnosed by caudal-fin or other characteristics.

TABLE 2 | Revised classification of Doradidae Bleeker, 1858. Nominal species that remain questionable as valid preceded by “?” and listed under possible senior synonym. Totals exclude species that are questionably valid, and species introduced to or questionably present in a given basin. Asterisk denotes species included in molecular phylogenetic analyses.