Specimens were collected from two sites in the Mukwadzi River, a tributary of the Manyame River, a south bank affluent of the middle Zambezi River, during surveys in 2016 and 2019. Samples were collected using a battery-powered Samus 725GN backpack electric fisher with a block net placed downstream to capture dislodged animals in the fast-flowing current. The specimens were photographed to document the live colour pattern then euthanized with clove oil. Muscle tissue from the right side of the specimens was cut out and preserved in 99% ethanol for molecular analysis. Voucher specimens for morphological studies were fixed in 10% formalin in the field and subsequently transferred to 70% ethanol for long term preservation. Additional tissue samples and voucher specimens used in the present study were obtained from the National Fish Collection at the National Research Foundation-South African Institute for Aquatic Biodiversity (NRF-SAIAB) in Makhanda (Tables 1, 2).

A total of six new COI sequences of Chiloglanis carnatus sp. nov. were generated for this study. Preparation and sequencing of genetic material was done in the Aquatic Genomics Research Platform at the NRF-SAIAB. Genomic DNA was extracted from preserved tissues using the salting-out method (Sunnucks and Hales 1996). The mitochondrial DNA cytochrome c oxidase subunit I (COI) gene was amplified by polymerase chain reaction (PCR) using the universal fish DNA barcoding primer set FishF1 and FishR1 (Ward et al. 2005). PCRs were performed with a Veriti 96 well thermal cycler (Applied Biosystems, USA) and each reaction mixture (25 µL) contained 50–100 ng) of template DNA, 6.5 µL of water, 0.5 µL of each primer (10 µM), and 12.5 µL Taq DNA polymerase 2× master mix red (Amplicon PCR enzymes and reagents, Denmark). The PCR amplification profile had an initial denaturation step of 3 min at 94 °C followed by 38 cycles of 30 sec at 94 °C, annealing at 55 °C for 30 sec, and extension at 72 °C for 50 sec, and final extension at 72 °C for 7 min. The amplicons were purified using an Exonuclease I-Shrimp Alkaline Phosphate (Exo/SAP, Thermo Fisher Scientific, USA) protocol (Werle et al. 1994), sequenced using standard fluorescent BigDye v. 3.1 (Applied Biosystems, USA) terminator chemistry in the forward direction, and analysed on a 3500 Genetic Analyser (Applied Biosystems, USA) at the NRF-SAIAB. Additional sequences were obtained from the public databases GenBank (https://www.ncbi.nlm.nih.gov/genbank/) and Barcode of Life Data Systems (BOLD) (http://www.boldsystems.org/) (Table 1).

Phylogenetic analyses included genetic sequences generated from Chiloglanis carnatus sp. nov., six of the seven nominal species from southern Africa, six candidate species of Chiloglanis identified by Chakona et al. (2018) (Chiloglanis sp. ‘roughskin’, Chiloglanis sp. ‘Zambezi’, Chiloglanis sp. ‘Nyangombe’, Chiloglanis sp. ‘Pungwe’, Chiloglanis sp. ‘Shire’, Chiloglanis sp. ‘dwarf’), and three outgroup species (Euchilichthys boulengeri Nichols & LaMonte, 1934; Euchilichthys royauxi Boulenger, 1902; Atopochilus savorgnani Sauvage, 1879) (Table 1). Genetic material for C. neumanni from its type locality and C. emarginatus could not be accessed before finalising this study. Mitochondrial DNA sequences were edited, aligned, and trimmed in MEGA-X (Kumar et al. 2016). The sequences were translated into amino acid sequences in MEGA-X to check for stop codons and gaps to ensure that they were copies of functional mitochondrial protein coding sequences. Haplotype groups were identified using DNASP 6 (Rozas et al. 2017). The most suitable model for nucleotide substitution was selected using the Akaike Information Criterion (AIC) (Akaike 1974) as implemented in the program jModelTest 0.1.1 (Darriba et al. 2012). Bayesian phylogenetic inference was performed in MrBayes 3.2.6 (Ronquist et al. 2012) using the TIM3+I+G evolutionary model identified using jModeltest. The phylogenetic tree and posterior probabilities were inferred using four Markov chain Monte Carlo (MCMC) chains which were run for 2 × 10⁶ generations with tree sampling every 1000 generations. The program Tracer 1.7 (Rambaut et al. 2018) was used to analyse the quality of the trace files generated by MrBayes and to determine the number of trees to be discarded as burn-in. The first 25% of the sampled trees for each analysis was discarded as burn-in, and the remaining trees were used to calculate a majority rule consensus tree. Maximum likelihood (ML) analysis of the same dataset was performed in RAxML v. 8 (Stamatakis 2014) through the graphical user interface raxmlGUI v. 2 (Silvestro and Michalak 2012). A total of 10 ML searches were performed in raxmlGUI and support values for the ML tree nodes were estimated by 1000 non-parametric bootstrap inferences (Felsenstein 1985). Bootstrap values equal to or higher that 70% (Hillis and Bull 1993), and posterior probability values at 0.95 or higher (Alfaro and Holder 2006), were considered strong support.

Four molecular species delimitation methods were used to delineate candidate species within the suckermouth catfishes of southern Africa using the same dataset used for the phylogenetic analysis. The first two methods, Automatic Barcode Gap Discovery (ABGD; Puillandre et al. 2012) and Assemble Species by Automatic Partitioning (ASAP; Puillandre et al. 2021) infer the barcode gap from the data to partition sequences into proposed candidate species. These methods were performed on their respective webservers (https://bioinfo.mnhn.fr/abi/public/abgd/ and https://bioinfo.mnhn.fr/abi/public/asap/asapweb.html). The intraspecific diversity priors were set at P_min = 0.001 and P_max = 0.1) for both methods. The Kimura (K80) TS/TV distance model was used and the remaining settings were left at their default parameters. The second pair of species delimitation methods included the Bayesian implementation of the Poisson Tree Processes (bPTP) (Zhang et al. 2013) and the General Mixed Yule Coalescent (GMYC) (Pons et al. 2006; Fujisawa and Barraclough 2013). Both GMYC and bPTP require a phylogenetic tree as input and from this tree they estimate rates of branching events to infer which parts of the tree are likely to follow a speciation model (interspecific variation) and which parts follow a coalescent model (intraspecific variation). The bPTP was performed on the web server (http://species.h-its.org/ptp/) using the same tree generated for phylogenetic reconstruction and a MCMC run for 1 × 10⁶ generations with 10% burn-in. For the GMYC analysis a fully resolved ultrametric tree was inferred in Bayesian evolutionary analysis by sampling trees (BEAST) 2.4.6 (Bouckaert et al. 2014) using a strict clock and Yule model and the MCMC was ran for 1 × 10⁷ generations with tree sampling every 1000 generations. The program Tracer 1.7.2 was used to analyse the quality of the trace files generated by BEAST. TreeAnnotator (Helfrich et al. 2018) was used to summarise the trees sampled by BEAST into a single maximum credibility tree with a burn-in of 25%. The species’ limits by threshold Statistics (splits) package (http://r-forge.r-project.org/projects/splits) in R 3.5.0 (R Core Team 2018) was used to identify the candidate species from the maximum credibility tree produced by TreeAnnotator. Model corrected intraspecific and interspecific genetic distances of the molecular taxonomic units identified by the species delimitation methods were calculated in PAUP* 4.0a163 (Swofford 2003).

A total of 19 specimens of Chiloglanis carnatus sp. nov. collected from the Mukwadzi River were examined in the present study. Comparative material included the lectotype of C. neumanni, holotypes of six valid species from southern Africa and five candidate species identified by Chakona et al. (2018) (Table 2). Because type material for C. fasciatus could not be accessed before finalising this study, 10 conspecific specimens collected from close to the type locality of this species and identified using the species description by Pellegrin (1936) and the key in Skelton (2001) were used as topotypes. The syntypes of C. pretoriae were severely deformed, thus only their meristic counts were included for comparison in this study, but 16 specimens collected from near the type locality of C. pretoriae and identified using the key in Skelton (2001) were used as topotypes for this species. Formulae and terminology of morphometric and meristic characters followed Schmidt et al. (2015), Friel and Vigliotta (2008), and Skelton and White (1990). A total of 49 morphometric characters were measured to the nearest 0.1 mm using digital Vernier callipers following Friel and Vigliotta (2008) (Table 3, Fig. 2A–D). External meristic counts were performed under a stereo microscope. Vertebrae counts were made from radiographs taken at the NRF-SAIAB using an Inspex 20i Digital X-ray Imaging System (Kodex Inc., New Jersey, USA). Radiographs for the lectotype and paralectotypes of C. neumanni were taken at the Royal Museum for Central Africa, Belgium (MRAC) using a VisiX-MedexLoncin (www.medex.be). A total of nine meristic characters were examined: number of mandibular teeth, pre-maxillary teeth, pectoral-fin rays, pelvic-fin rays, dorsal-fin rays, anal-fin rays, abdominal vertebrae, caudal vertebrae, and total vertebrae (Fig. 3). Following Roberts (1989), vertebrae counts excluded the Weberian structures and began from the fifth vertebra which was identified by a pair of large but slender ribs, and included the hypural complex which was counted as one vertebra. The abdominal vertebrae were defined as the vertebrae that occurred anterior to the first anal fin ray pterygiophore. Caudal vertebrae were defined as those that occurred posterior to the first anal fin ray pterygiophore and included the hypural complex which was counted as one vertebra (Roberts 1989) (Fig. 3). The genital papillae were examined to determine the sex of the specimens following Friel and Vigliotta (2008). Morphological measurements were standardised by transforming body measurements into percentages of the standard length (SL) and head measurements into percentages of the head length (HL). Principal component analyses (PCA) were performed in PAST v. 3.12 (Hammer et al. 2001) using the covariance matrix for the morphometric data in order to identify morphological characters that contributed the most to distinguishing Chiloglanis carnatus sp. nov. from the other Chiloglanis species from southern Africa.

The COI alignment of 80 sequences had 534 base pairs and 176 variable sites. A total of 47 unique haplotypes were identified. Although the Bayesian phylogenetic tree was not fully resolved, it showed genetic structuring that supported the monophyly of suckermouth catfishes from southern Africa (Fig. 4). Chiloglanis carnatus sp. nov. was recovered as an exclusive group that is genetically divergent (2.8–15.0% genetic distances) from other Chiloglanis species and lineages from southern Africa (Figs 4, 5; Table 4). With the exception of C. pretoriae and Chiloglanis sp. ‘Shire’, all recovered clades were well-supported (posterior probability > 0.95). Genetic divergences within valid and candidate species ranged from 0–1.5% and interspecific divergences ranged from 1.3–15.7% (Table 4). Chiloglanis paratus from the Phongolo River was recovered as the most basal clade that is sister species to all the southern African suckermouth catfishes included in the present study. The relationships among the remaining taxa were not well resolved as they were recovered within five polytomous clades with weak support between them (Fig. 4). The first clade contained Chiloglanis carnatus sp. nov. from the Manyame River and C. fasciatus from the Okavango River. The second clade contained C. anoterus and C. bifurcus from the Incomati River system as well as C. pretoriae from the Limpopo River system. The third clade contained Chiloglanis sp. ‘Nyangombe’ from the Nyangombe River and Chiloglanis sp. ‘dwarf’ from the Pungwe and Ruo rivers. The fourth clade contained Chiloglanis swierstrai from the Limpopo River system. The fifth clade contained the Chiloglanis sp. ‘Zambezi’, Chiloglanis sp. ‘Pungwe’, Chiloglanis sp. ‘roughskin’, and Chiloglanis sp. ‘Shire’ lineages. The Chiloglanis sp. ‘roughskin’ lineage occurs in the Buzi and Pungwe rivers, whereas Chiloglanis sp. ‘Pungwe’ is endemic to the Pungwe River. Chiloglanis sp. ‘Shire’ and Chiloglanis sp. ‘Zambezi’ lineages were found in the lower Zambezi River system with the latter lineage also occurring in the Okavango River. The phylogenetic tree inferred using the ML approach had similar topology to the Bayesian inference tree (Fig. 5).

All four molecular species delimitation methods identified Chiloglanis carnatus sp. nov., Chiloglanis sp. ‘Shire’, Chiloglanis sp. ‘Nyangombe’, C. swierstrai, C. anoterus, C. pretoriae, C. bifurcus, C. fasciatus, and C. paratus as unique molecular taxonomic units (Fig. 4). The Assemble Species by Automatic Partitioning method recovered the least number of candidate species, this method grouped Chiloglanis sp. ‘roughskin’, Chiloglanis sp. ‘Pungwe’, and Chiloglanis sp. ‘Zambezi’ into a single molecular taxonomic unit. The General Mixed Yule Coalescent method recovered the highest number of molecular taxonomic units. This method identified additional molecular taxonomic units within Chiloglanis sp. ‘roughskin’ and Chiloglanis sp. ‘dwarf’. The Automatic Barcode Gap Discovery and bPTP inferred similar molecular taxonomic units with the exception of Chiloglanis sp. ‘dwarf’ which was split into two molecular taxonomic units by the latter method.

Principal component analysis (PCA) of the morphometric characters showed that Chiloglanis carnatus sp. nov. is separated from C. swierstrai and C. anoterus along principal component 1 (PCI) (Fig. 6). This separation was associated with maxillary barbel length (Table 6). Chiloglanis carnatus sp. nov. (20.3–28.8%HL) has relatively shorter maxillary barbels compared to C. swierstrai (44.2–66.8%HL) and Chiloglanis sp. ‘Zambezi’ (31.3–37.0%HL, Table 5, Fig. 7A, B). Chiloglanis carnatus sp. nov. is separated from C. swierstrai, C. anoterus, and C. neumanni along principal component 2 (PCII) (Fig. 6). Separation along PCII is associated with the oral disc width (Table 6). Chiloglanis carnatus sp. nov. has a relatively smaller oral disc width (51.1–64.6%HL) compared to C. anoterus (69.1%HL, Table 5, Fig. 7C).

Additional scatterplots were generated to explore the characters that further distinguish the Chiloglanis carnatus sp. nov. specimens. The Chiloglanis carnatus sp. nov. specimens have a narrower mandibular tooth row width (4.6–8.1%HL) compared to C. pretoriae (16.0–25.6%HL), C. swierstrai (10.0–16.6%HL), C. neumanni (9.9–13.5%HL), C. emarginatus (9.6–13.5%HL), C. bifurcus (10.4–17.3%HL), C. anoterus (10.5%HL), Chiloglanis sp. ‘dwarf’ (13.6–25.0%HL), Chiloglanis sp. ‘Nyangombe’ (19.0–25.5%HL), Chiloglanis sp. ‘Pungwe’ (17.6–27.8%HL), Chiloglanis sp. ‘roughskin’ (11.9–22.2%HL), and Chiloglanis sp. ‘Zambezi’ (20.5–25.4%HL; Fig. 7D, E). Chiloglanis carnatus sp. nov. has an oral disc with relatively longer lower lips (18.3–26.6%HL) compared to Chiloglanis sp. ‘roughskin’ (9.6–16.8%HL, Fig. 7F). Chiloglanis carnatus sp. nov. has a relatively deeper caudal peduncle (11.3–13.2%SL) compared to C. neumanni (9.5–9.9%SL), C. paratus (9.6–9.9%SL), C. fasciatus (7.5–8.8%SL), C. swierstrai (7.2–8.7%SL), and Chiloglanis sp. ‘Zambezi’ (10.0–11.1%SL, Fig. 7G, H). A longer adipose-fin base length distinguishes Chiloglanis carnatus sp. nov. (17.0–23.3%SL) from C. bifurcus 9.2–13.6%SL) (Fig. 7I). Larger adipose fin height (4.1–6.8%SL) and shorter anal fin rays (11.7–17.9%SL) further distinguish Chiloglanis carnatus sp. nov. from C. neumanni (adipose fin height: 2.7–3.1%SL; anal fin ray length: 19.2–20.9%SL; Fig. 7J, K). A shorter distance between the anterior nares of Chiloglanis carnatus sp. nov. (9.5–15.5%HL) separates it from C. bifurcus (19.5–21.2%HL), C. emarginatus (16.5–22.4%HL), and C. swierstrai (15.7–22.4%HL, Fig. 7L). Chiloglanis carnatus sp. nov. has a relatively longer head (30.5–34.9%SL vs 24.8–28.0%SL), relatively wider body at pectoral-fin insertion (23.0–25.3%SL vs 17.6–21.5%SL), and relatively longer pre-pelvic (56.0–59.3%SL vs 49.1–54.8%SL) and pre-dorsal distances (39.9–43.7%SL vs 34.5–36.7%SL) that readily separated them from C. swierstrai (Fig. 7M–P).

Comparison of meristic characters revealed consistent differences between Chiloglanis carnatus sp. nov. specimens and the other species from southern Africa. Chiloglanis carnatus sp. nov. specimens have ten closely packed mandibular teeth that separate them from C. bifurcus, C. emarginatus, C. fasciatus, and C. neumanni that have eight mandibular teeth as well as from C. anoterus, C. pretoriae, C. paratus, and C. swierstrai that have > 10 mandibular teeth (Fig. 8A). A higher number of anal-fin rays separates Chiloglanis carnatus sp. nov. specimens (12–13) from C. paratus (9) and C. pretoriae (10) (Fig. 8B). Chiloglanis carnatus sp. nov. specimens have a higher number of total vertebrae (29–30) compared to C. neumanni (28), C. pretoriae (28), and C. paratus (27) (Fig. 8C).

The integrated approach used in this study provided genetic and morphological characters that clearly and consistently distinguish Chiloglanis carnatus sp. nov. from the known species and lineages from this region. This study has thus provided evidence that supports the description of the Chiloglanis carnatus sp. nov. as a new species.