Single specimens of C. aeneus and C. paleatus were collected from the temple of Confucius flower and wood fish market, Nanjing city, Jiangsu province, China (32°0'27.1"N, 118°50'11.5"E) in June 2020 and identified based on their morphological characteristics, according to the latest taxonomic classification of fish (Popazoglo and Boeger 2000; Huysentruyt and Adriaens 2005a, b). Their geographic data and specific origins were unknown. All fresh tissues were immediately stored at -80 °C in 95% ethanol until DNA extraction. Total DNA was extracted from the muscle tissue using a TIANamp Marine Animals DNA Kit DP324 (Tiangen Biotech Co., Ltd., Beijing, China) according to the manufacturer’s instructions. DNA integrity and purity were evaluated by 1% agarose gel electrophoresis, and DNA purity was determined with a NanoDrop 2000 (NanoDrop Technologies, Wilmington, DE, USA). DNA concentrations were quantified using a Qubit^R 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA). To ensure the accuracy of morphological identification, COI primers were designed based on the latest DNA barcoding database (NCBI and FishBase) and were amplified, sequenced, and compared. The COI sequences are provided in the Suppl. material 1. The results of the sequence alignment verify the accuracy of the morphological identification.

Next-generation sequencing was performed to determine the complete mitogenome sequence of the two species of Corydoras. The DNA libraries were sequenced on an Illumina sequencing platform by Novogene Co., Ltd. (Beijing, China). Briefly, the total DNA genome was quantified and fragmented into 250-base pair (bp) fragments using a Covaris M220 ultrasonic crushing system (Woburn, MA, USA) followed by whole-genome shotgun sequencing. According to the manufacturer’s instructions, a library was constructed based on two indices using an Illumina TruSeq DNA PCR-Free HT kit (San Diego, CA, USA). An Illumina Novaseq 6000 platform was used for sequencing of 150 paired-end reads approximately 4 Gb in size. Clean reads were generated as previously described, and the remaining high-quality reads were assembled using SPADES V3.15.2 (Bankevich et al. 2012) (http://cab.spbu.ru/software/spades/) and SOAPDENOVO2 V2.01 (Luo et al. 2012) software. The preliminary assembly results were compared with the NT database, and looped sequences annotated as mitochondrial genomes were screened. CAP3 was used to merge the splicing results from the two software programs, and the assembly results were compared with those of related species using MUMMER v3.23 (Delcher et al. 2003). The mitogenome composition was confirmed, and a complete, high-quality map of the mitochondrial genome was obtained.

The tRNA genes were verified using tRNASCAN-SE V1.3.1 (Lowe and Eddy 1997) with default settings for the vertebrate mitochondrial genetic code. The software, which integrates multiple analysis tools, can identify 99% of the tRNA genes with a very low number of false positives and predict the secondary structure of tRNAs. Protein-coding regions were re-identified using GLIMMER V3.0 (Ingram et al. 2009), and manual comparisons were performed using the SEQMAN program of LASERGENE V7.1 (Burland 2000) (DNAStar, Inc., Madison, WI, USA) based on the PCGs of nine species of Corydoras and translated into putative proteins via GenBank. The non-coding RNAs were verified using RFAM V12.0 (Griffiths-Jones et al. 2003) and INFERNAL V1.1 (Nawrocki and Eddy 2013). The rRNA genes were assumed to extend to the boundaries of flanking genes, similar to the homologous regions of other published mitogenomes of Corydoras in GenBank. The MITOS WebServer (http://mitos2.bioinf.uni-leipzig.de/index.py) and MitoFish (Iwasaki et al. 2013) (http://mitofish.aori.u-tokyo.ac.jp/) online tools were used for the final annotation of the entire mitogenome sequence of the two species of Corydoras, and the annotated mitogenomes were compared with nine published mitogenomes of Corydoras. Base compositions, genetic distances, and relative synonymous codon usage values were determined using MEGA V7.0 (Kumar et al. 1994). A graph comparing the relative synonymous codon usage was drawn using PHYLOSUITE V1.2.2 (Zhang et al. 2020). Strand asymmetry was analyzed using the formula: AT-skew = (A – T)/(A + T). The numbers of non-synonymous (Ka) and synonymous (Ks) substitutions and the ratio of Ka/Ks and nucleotide diversity for the nine species of Corydoras were calculated using DNASP 5.1 (Librado and Rozas 2009). The MitoFish (http://mitofish.aori.u-tokyo.ac.jp/) online tool was used to generate circular mitogenome maps.

Phylogenetic trees for the eleven mitogenomes of Corydoras within the family Callichthyidae and Suborder Loricarioidei were constructed by aligning 13 PCGs and two rRNA sequences with those of 42 species of Loricarioidei, 29 species from Loricariidae, and one species from Trichomycteridae (Table 1). The mitogenomes of Pterocryptis cochinchinensis (Resende et al. 2016) and Silurus asotus (Nakatani et al. 2011) (accession no. NC_027107.1 and NC_015806.1, respectively, suborder Siluroidei) were included as outgroups to root the Loricarioidei tree. All operations were performed in PHYLOSUITE V1.2.2 (Zhang et al. 2020) software package. The nucleotide sequences of 13 PCGs from 44 mitogenomes were aligned in batches with MAFFT V7.313 (Katoh and Standley 2013) (https://mafft.cbrc.jp/alignment/server/) using the codon alignment mode. The results were optimized using MACSE V2.03 (Ranwez et al. 2018). The nucleotide sequences of two rRNAs were aligned using the online tool MAFFT with default settings. Ambiguously aligned regions were removed via GBLOCKS 0.91 b with default settings. The resulting alignments were concatenated into a single dataset with PHYLOSUITE. The best partition schemes and optimal substitution models were selected by MODELFINDER (Kalyaanamoorthy et al. 2017) with the greedy algorithm and Bayesian information criterion (Watanabe 2013). The best substitution models applied to each partition are listed in Suppl. material 1: Table S1. Phylogenetic trees were constructed using two inference methods: maximum likelihood (ML) and Bayesian inference (BI). ML analyses were performed with IQ-TREE V1.6.8 with the models selected for each partition, and 1,000 bootstrap replicates were used to estimate node reliability. Bayesian analyses were performed using MRBAYES V3.2.6 (Huelsenbeck and Ronquist 2001). One million generations of two independent runs were performed with four chains and sampling trees every 100 generations. The initial 25% of trees generated prior to reaching stable log-likelihood values were discarded as burn-in. The remaining trees were used to calculate the Bayesian posterior probabilities. The resulting phylogenetic trees and gene orders were visualized and edited using iTOL (Letunic and Bork 2016).

The complete mitogenomes of C. aeneus and C. paleatus comprising 16,604 and 16,593 bp, respectively, were submitted to GenBank (accession nos. MZ571336 and MZ571337, respectively) (Fig. 1, Table 2). The two mitogenomes were circular and contained 37 mitochondrial genes (13 PCGs, 22 tRNA genes, and two rRNA genes) and one D-loop. The position of each gene in the mitogenome was identical to that in other species of Corydoras (Saitoh et al. 2003; Moreira et al. 2016a, 2017; Liu et al. 2019a, b, c, d; Chen et al. 2020; Lv et al. 2020). One of the 13 PCGs (ND6) and eight tRNAs (tRNA-Ala, tRNA-Cys, tRNA-Glu, tRNA-Asn, tRNA-Pro, tRNA-Gln, tRNA-Ser(TGA), and tRNA-Tyr) were encoded by the light chain (-), whereas the other 28 genes, including 12 PCGs, 14 tRNAs, two rRNAs, and one D-loop, were encoded by the heavy chain (+) (Fig. 1, Table 2). The 44 mitogenomes of Siluriformes (Nakatani et al. 2011; Liu et al. 2016; Moreira et al 2016b, 2018; Resende et al. 2016; Magalhães et al. 2017; Parente et al. 2017; Parente et al. 2018; Pereira et al. 2019; Ren et al. 2019; Rocha-Reis et al. 2020; Meng et al. 2021; Zhang et al. 2021) used in this study were compared, and the gene composition and order were consistent (Suppl. material 1: Fig. S1). The nucleotide composition of the two entire mitogenomes was as follows: C. aeneus A = 5417 (32.63%), T = 4299 (25.89%), G = 2451 (14.76%), C = 4437 (26.72%) and C. paleatus A = 5380 (32.42%), T = 4282 (25.81%), G = 2481 (14.95%), C = 4450 (26.82%). The two mitogenomes (values for C. aeneus followed by values for C. paleatus) had high A+T contents of 58.52% and 58.23% (Suppl. material 1: Table S2), including 58.08% and 57.67% in PCGs, 56.97% and 57.04% in tRNA genes, 59.70% and 59.10% in 16S rRNA, 55.30% in 12S rRNA, and 67.51% and 68.21% in the D-loop, respectively, which agrees with the typical base bias of fish mitogenomes (Gadaleta et al. 1989; Manchado et al. 2007; Xu et al. 2011). The overall AT and GC skew values in the entire mitogenome of C. aeneus were 0.115 and -0.288 and in C. paleatus were 0.114 and -0.284, respectively. The GC skew value of the eleven mitogenomes of Corydoras, except for tRNA, was slightly negative (-0.014 to -0.288), showing a higher occurrence of C than of G. In contrast, AT skew value, except for the second codon position, was slightly positive (0.028 to 0.379), showing a higher content of A than of T. The K2P genetic distances of the eleven mitogenomes of Corydoras were all less than 0.12 (Suppl. material 1: Table S3). C. nattereri and C. sterbai and C. nattereri and C. trilineatus showed the largest K2P genetic distances among the eleven species of Corydoras.

Figure 1. 

Gene maps of the two newly sequenced Corydoras species.

The 13 PCGs of the two new mitogenomes and those of the previously published nine mitogenomes of Corydoras contained COI–COIII, ND1–ND6, ND4L, two ATPases, and one Cyt-b, similar to that in other Siluriformes (Nakatani et al. 2011; Liu et al. 2016; Moreira et al. 2016b; Resende et al. 2016; Magalhães et al. 2017; Parente et al. 2017; Moreira 2018; Parente et al. 2018; Pereira et al. 2019; Ren et al. 2019; Rocha-Reis et al. 2020; Meng et al. 2021; Zhang et al. 2021). The total lengths of PCGs in the eleven mitogenomes of Corydoras were 11,400–11,414 bp, accounting for 67.84–69.24% of the entire mitogenome. Similar to the mitogenomes of other species of Loricarioidei, ND5 and ATPase 8 were largest (1,827 bp) and smallest (168 bp), respectively. Most PCGs stringently start with an ATG start codon, except for all COIs, which start with GTG, C. nattereri COIII (Moreira et al. 2016a) which starts with GCA, and C. schwartzi COII (Moreira et al. 2017), which starts with CCA (Suppl. material 1: Table S4). Most PCGs are stringently terminated by the stop codon TAR (TAA/TAG) or an incomplete stop codon T, except for all COIs, which terminate with AGG and C. schwartzi ATPase 6 and C. nattereri ND3, which terminate with TA. The presence of a truncated stop codon is common among vertebrate mitochondrial genes and is thought to be introduced by posttranscriptional poly-adenylation.

Similar to most previously sequenced members of Loricarioidei, the AT-skews (0.033 to 0.052) and GC-skews (-0.268 to -0.299) of the PCGs were similar among the eleven species of Corydoras (Suppl. material 1: Table S2). Summaries of the relative synonymous codon usage and the number of amino acids in the annotated PCGs are presented in Suppl. material 1: Figs S2, S3. The PCGs of the eleven mitogenomes of Corydoras (Saitoh et al. 2003; Moreira et al. 2016a, 2017; Liu et al. 2019a, b, c, d; Chen et al. 2020; Lv et al. 2020) translate into 3,798–3,802 codons and showed very similar codon usage, excluding the stop codons (26–28 bp). Ile (310.82 ± 2.69 codons), Thr (312.64 ± 2.27 codons), Ala (312.73 ± 3.08 codons), and Leu1 (CUN) (475.45 ± 12.89 codons) were the four most predominant codon families and may be associated with the coding function of the chondriosome. In contrast, Cys (24.91 ± 0.79 codons) and Ser1 (AGN) (52.18 ± 0.83 codons) had the smallest number of codons. A/T rather than G/C bias was observed in the third position, as almost all frequently used codons ended with A/T. The synonymous codon preferences for the eleven species of Corydoras were conserved, possibly because of the close relationships among members of this genus.

To reveal the evolutionary pattern of the PCGs, the Ka/Ks, nucleotide diversity, and K2P genetic distance across all mitogenomes of Corydoras were calculated for each aligned PCG. The K2P genetic distances of 13 PCGs were all less than 0.12 (Fig. 2a). Among the PCGs detected, ND4 and ATPase 8 showed the largest K2P genetic distance among the eleven species of Corydoras, followed by ND2 and ND3. The nucleotide diversity of the 13 PCGs was less than 0.11 (Fig. 2b). ND4 showed the highest nucleotide diversity, whereas COII showed the lowest diversity. To investigate the selective pressure across species of Corydoras, the Ka/Ks ratios of the PCGs of each mitogenome were estimated (Fig. 2c). The Ka/Ks value was highest for ND6, followed by ND2; the lowest values were observed for COI, COIII, ND1, and ND4L. All 13 PCGs showed Ka/Ks << 1, suggesting that all PCGs of Corydoras evolved under purifying selection.

Figure 2. 

K2P genetic distance a nucleotide diversity b Ka/Ks ratio c analyses of protein-coding genes among the eleven Corydoras mitogenomes.

The total lengths of the 22 tRNA genes ranged from 1,438 (C. schwartzi) to 1,561 bp (C. arcuatus and C. panda), whereas individual tRNA genes typically ranged from 58 to 75 bp. All tRNA genes displayed the expected cloverleaf secondary structures with normal base pairing, except for tRNA-Ser(GCT), which lacked the DHU stem (Suppl. material 1: Fig. S4), forming a loop commonly found in other vertebrates (Ojala et al. 1981; Gadaleta et al. 1989; Wolstenholme 1992). The A+T contents of these tRNAs were 56.55–57.58%. All AT-skew and GC-skew values were slightly positive, indicating a slight bias toward the use of A and G in the tRNAs (Suppl. material 1: Table S2). These rRNA genes are between tRNA-Phe and tRNA-Leu(TAA) and are separated by tRNA-Val. The average total size of the two rRNAs was 2,614 bp, and the average A+T content was 57.89%. Like the tRNAs, all AT-skew values were positive, whereas all GC-skew values were negative, indicating that rRNAs favor C compared to tRNAs in Corydoras.

The control region (D-loop), also known as the A+T rich region that contains hypervariable non-coding sequences and regulates the replication and transcription of mitochondrial DNA, is the largest non-coding region and is located between tRNA-Pro and tRNA-Phe in these mitogenomes. Compared with PCGs, the D-loop displayed a higher mutation rate and the highest variation throughout the mitogenome; thus, this region is dominant and can be used to evaluate intraspecies variations. The D-loops in the eleven species of Corydoras were 718‒1,218 bp. Compared with the other four regions (entire genome, PCGs, tRNAs, and rRNAs), the control region showed the highest A+T content, ranging from 66.77% to 71.87%. Like the rRNAs, all AT-skew values were positive, and all GC-skew values were negative.

To determine the phylogenetic relationships within the suborder Loricarioidei and family Callichthyidae, we obtained the concatenated nucleotide sequences of 13 PCGs and two rRNAs from 42 species of Loricarioidei. Phylogenetic analyses based on both ML and BI methods revealed same topologies, which also generally agreed with those presented in previous studies (Alexandrou et al. 2011; Lujan et al. 2015; Moreira et al. 2017; Roxo et al. 2019) (Figs 3, 4). These analyses confirmed that the genus Corydoras was part of the monophyletic family Callichthyidae.

Figure 3. 

Phylogenetic trees of 44 Siluriformes species using concatenated nucleotide sequences of 13 protein-coding genes and two rRNAs using the maximum likelihood method. Numbers in the ML tree represent SH-aLRT support/ultrafast bootstrap support values.

Both Callichthyidae and Loricariidae were recovered as monophyletic with very high support values (BI posterior probabilities, PP = 1; ML bootstrap, BS = 100). The 44 species of Siluriformes were divided into four major clades corresponding to the families Siluridae Callichthyidae, Trichomycteridae, and Loricariidae. The target species C. aeneus and C. paleatus were clustered into two clades (C. aeneus + C. rabauti) and (C. paleatus + C. nattereri) with a high nodal support value (PP = 1; BS = 100). The eleven species of the genus Corydoras clustered together quite well [((C. aeneus + C. rabauti) + (C. schwartzi + C. agassizii)) + (C. arcuatus + (C. panda + (C. duplicareus + (C. sterbai + C. trilineatus))))] + [(C. paleatus + C. nattereri)]. Corydoras trilineatus and C. sterbai have short, almost non-existent branch lengths; thus, they are likely the same species. The K2P genetic distances of these two species are 0.000 (Suppl. material 1: Table S3), which verifies that they are the same species. This may be caused by incorrect identification, taxonomic problems (these two species are, in fact, synonymous), and/or introgressive hybridization. Moreover, in the family Loricariidae, the genera Ancistrus and Loricariichthys were clustered into monophyletic clades [(A. cryptophthalmus + A. multispinis) + A. temminckii] and ( L. castaneus + L. platymetopon) with a high nodal support value (PP = 1; BS = 100). There was a paraphyletic relationship between the genera Hypostomus and Pterygoplichthys, [H. francisci + (H. ancistroides + H. affinis), P. pardalis + (H. plecostomus + (P. anisitsi + P. disjunctivus))]. Our results demonstrate that the concatenated nucleotide sequences of the 13 PCGs and two rRNAs were useful for determining the phylogenetic relationships of the order Siluriformes. These results can be used to improve classification of the families Callichthyidae and Loricariidae.

Figure 4. 

Phylogenetic tree of 44 Siluriformes species using concatenated nucleotide sequences of 13 protein-coding genes and two rRNAs via the Bayesian interference method. Applicable posterior probability values are shown.