All procedures were conducted in accordance with relevant ethical guidelines. Specimen collection complied with applicable legislation in China and followed protocols approved by the Animal Ethics Committee of Nanjing Forestry University. All experimental protocols adhered to international standards for animal welfare, including the Convention on Biological Diversity (CBD), the Nagoya Protocol, and the ARRIVE guidelines (https://arriveguidelines.org). Samples (in vitro-bred ornamental fish of Chinese origin) were collected on 1 April 2022 from an aquatic market in Qinhuai District, Nanjing City, Jiangsu Province (32.005899°N, 118.841977°E). Tail fin tissues were immediately immersed in anhydrous ethanol and stored at −20°C. Genomic DNA was extracted from the samples using the CTAB method (Porebski et al., 1997). The integrity of the extracted DNA was assessed by agarose gel electrophoresis, while its purity was determined through spectrophotometric analysis.

Genomic DNA was fragmented using a Covaris ultrasonicator (Covaris, USA) to generate randomly sheared fragments. Subsequently, PCR amplification and other standard procedures were performed for library preparation. The constructed library was diluted following preliminary quantification with a Qubit 3.0 Fluorometer (Nakayama et al., 2016), and the insert size distribution was determined using a Qsep100 Bioanalyzer (BiOptic Inc., New Taipei, Taiwan). Quantitative PCR (qPCR) was conducted to precisely measure the library’s effective concentration (with a threshold of >3 nM). Upon passing quality control, libraries were pooled at appropriate ratios based on their effective concentrations and the required sequencing depth. Following DNA library qualification, paired-end sequencing was performed on the Illumina HiSeq high-throughput sequencing platform, generating a minimum of 10 Gb of raw data.

Raw sequencing data underwent quality control using FastQC, followed by adapter trimming and quality filtering with Trimmomatic. This process retained >90% of reads, yielding trimmed sequences averaging ∼140 bp in length. The average sequencing depth exceeded 20× coverage. All library preparation and quality assessments were performed by Nanjing Qingke Biotechnology Co., Ltd. (Nanjing, China). Upon receiving the raw sequencing data from the company, we employed NOVOPlasty for sequence assembly and optimization (Dierckxsens et al., 2017). The pipeline consisted of the following steps: (1) using BLAST to identify complete mitochondrial genome sequences with high similarity to our input sequences; (2) locating and modifying the config.txt file in the NOVOPlasty directory by adjusting parameters including Type, Genome Range, and Reference sequence (default parameters were used in NOVOPlasty); (3) generating of complete mitochondrial genome sequence in FASTA format; (4) performing preliminary genome annotation through the MITOS web server (Bernt et al., 2013); (5) selecting reference sequences showing both high similarity and close phylogenetic relationships using BLAST; and (6) validating the MITOS preliminary annotations using selected reference sequences as benchmarks. Transfer RNA (tRNA) genes were identified using both MITOS2 and tRNAscan-SE v2.0. (Lowe and Eddy, 1997). Ribosomal RNA (rRNA) localization proved more challenging and was primarily inferred based on their conserved positional relationships with flanking genes.

We analyzed the base composition and codon usage frequency of the annotated mitochondrial genome sequence using MEGA software (version 11) (Kumar et al., 1994). The AT-skew and GC-skew values were calculated according to the following formulas: AT-skew = (A−T)/(A + T) and GC-skew = (G−C)/(G + C) (Perna and Kocher, 1995). Finally, we generated a circular representation of the complete mitochondrial genome using the MitoFish web server (Zhu et al., 2023).

The assembled mitochondrial genome of B. cupido, together with 42 New World cichlid mitogenomes from GenBank, was included in the phylogenetic analysis. An African cichlid mitogenome served as the outgroup. A total of 44 mitochondrial genomes were selected for phylogenetic analysis based on their mitochondrial protein-coding gene (CDS) sequences (Table 1). These 44 species mitogenomes were processed using PhyloSuite v1.2.3 (Zhang et al., 2020). Sequences were initially aligned with MAFFT (Katoh and Standley, 2013) using codon-based alignment, followed by optimization in MACSE (Ranwez et al., 2018). The concatenated dataset comprised 13 CDS sequences from all 44 species. This combined dataset was subsequently analyzed with the ModelFinder plugin for optimal data partitioning and gene-specific model selection (Kalyaanamoorthy et al., 2017). For phylogenetic reconstruction, we employed both Bayesian inference (BI) and maximum likelihood (ML) methods through the MrBayes (Ronquist et al., 2012) and IQ-TREE (Minh et al., 2020) plugins, respectively. The best-fit models selected were GTR + F + I + G4 for BI analysis and GTR + F + R5 for ML analysis.

Phylogenetic reconstruction using Bayesian inference (BI) was performed on the concatenated dataset of 13 protein-coding genes (CDS) derived from 44 selected species, utilizing the MrBayes plugin within PhyloSuite. The analysis ran two independent Markov chains, with the first 25% of generations discarded as burn-in. We performed 100 million generations in total, sampling every 100 generations. For maximum likelihood (ML) analysis, we employed the IQ-TREE plugin with 50,000 ultrafast bootstrap replicates. The resulting phylogenetic trees were subsequently visualized and annotated using the interactive Tree of Life (iTOL) online platform (Letunic and Bork, 2021).

The complete mitochondrial genome of B. cupido (GenBank Accession No. OP595705) measured 16,621 bp long, exhibiting a typical circular double-stranded structure (Figure 1; Table 2). This genome contains 37 functional genes and an 891 bp control region. Gene distribution analysis revealed that nine genes were encoded on the light (L) strand, comprising one protein-coding gene (nad6) and eight tRNA genes (trnQ, trnA, trnN, trnC, trnY, trnS2, trnE, and trnP). The remaining 28 genes resided on the heavy (H) strand, including 12 protein-coding genes, 14 tRNAs, and two rRNA genes. The genome organization showed 11 intergenic spacer regions ranging from 1 to 35 bp, with the largest spacer located between trnN and trnC. Additionally, we identified eight overlapping gene regions (1-19bp in length), with the most extensive overlap occurring between atp8 and atp6.

The base composition of the B. cupido mitochondrial genome was analyzed using MEGA11 (Table 3). The overall base frequencies are as follows: A = 27.5%, T = 26.6%, C = 30.9%, and G = 14.9%, demonstrating a significant A + T bias (54.1% combined) relative to G + C content. Skewness values further indicate compositional asymmetry: a positive AT-skew (0.017) suggests a slight preference for adenine (A) over thymine (T), while a negative GC-skew (−0.350) indicats a stronger bias toward cytosine (C) over guanine (G). For protein-coding regions (total length: 11,475 bp), the base composition was A = 25.0%, T = 28.6%, C = 31.8%, and G = 14.9%, maintaining the A + T bias (53.6%). Among codon positions, the second position exhibited the strongest A + T preference (58.3%).

The 12 CDS in B. cupido use the typical ATG as the start codon, whereas cox1 use GTG as the start codon. Termination codon usage is more variable: five CDS terminate with TAA; nad6 uses TAG, cox1 uses AGA and five genes (nad2, cox2, cox3, nad3, and nad4) end with an incomplete T, while atp6 terminates with an incomplete TA (Table 2). The incomplete termination codons are completed by the addition of 3′A to the mRNA. MEGA was used to statistically analyze the codon usage and amino acid contents of the sequences (Figures 2, 3). The most frequently noted amino acids were Leu1 (15.74%), Ala (8.62%), Thr (8.30%), Ile (7.41%), and Phe (6.55%), accounting for a total of 46.62%. Cys exhibited the lowest content (0.73%) (Figure 2). The CDS used 3,825 codons, with the most frequently used codon being CUC (224 times), with a relative synonymous codon usage (RSCU) value of 1.92, followed by CUU (192 times), with an RSCU value of 1.65. All codons encoded proteins (Figure 3).

The 22 tRNAs genes range in length from 68 to 75 bp. The longest tRNA gene is trnT (75 bp), and the shortest is trnC (68 bp). Among the 22 tRNA genes, trnQ, trnA, trnN, trnC, trnY, trnS2, trnE, and trnP are determined to be encoded on the light strand, whereas the other 14 tRNA genes are found to be encoded on the heavy strand. The total length of the tRNA genes is 1558 bp. The two rRNA genes, 12S ribosomal RNA (rrnL) is located between the trnV and trnL2 tRNAs, with a length of 1719 bp, whereas 16S ribosomal RNA (rrnS) is located between trnF and trnV, with a length of 952 bp. The A + T content of the two rRNA sequences is 52.4%, similar to the A + T content of the entire sequence but only slightly higher than the G + C content. The non-coding region (D-loop) is located between trnP and trnF, with a length of 891 bp.

For phylogenetic analysis, we selected mitochondrial genome CDS from 43 New World cichlid species representing five tribes (Astronotinae, Cichlasomatinae, Cichlinae, Geophaginae, and Retroculinae), along with one African cichlid species (subfamily Pseudocrenilabrinae) as the outgroup. Both maximum likelihood (ML) and Bayesian inference (BI) trees were constructed using the concatenated sequences of 13 protein-coding genes (Figures 4, 5). Most nodes in these phylogenetic trees showed strong statistical support. The topologies of the ML and BI trees were largely congruent, with minor discrepancies primarily attributed to variations in the phylogenetic placement of Rocio octofasciata.

As shown in Figures 4, 5, all 42 New World Cichlids (except for Heros Severus) were grouped into a monophyletic group well supported (bootstrap support = 99%; Bayesian posterior probability = 1.0). Among the five New World Cichlid tribes, Geophaginae, Astronotinae, Cichlasomatinae and Cichlinae were found to be polyphyletic; only Retroculinae was monophyletic (represented by a single species and thus not evaluated for monophyly). Both ML analysis and BI analysis supported a sister-group relationship between the tribe Cichlasomatinae (except Heros severus) and the tribe Geophaginae (bootstrap support = 58%; Bayesian posterior probability = 0.96). Additionally, both ML analysis and BI analysis supported a sister-group relationship between the tribe Retroculinae and the tribe Cichlinae (except C. regani) (bootstrap support = 82%; Bayesian posterior probability = 0.99). The phylogenetic analyses also revealed that species within the tribe Cichlini occupied a basal position in the tree, while species from the tribes Cichlasomatini (except H. severus) and Geophagini formed a clade at the apex of the evolutionary tree.

The tribe Astronotinae occupied an intermediate evolutionary position between Geophaginae and Cichlinae. Specifically, Astronotus ocellatus (Astronotinae) formed a sister group to Cichlinae, whereas Chaetobranchopsis bitaeniatus (Astronotinae) exhibited a sister relationship with Geophaginae. Similarly, with the exception of C. regani, the phylogenetic analysis indicated that the tribe Geophaginae was monophyletic. Tribe Geophaginae comprised two well-supported clades, both with maximum support values (bootstrap support = 100%; Bayesian posterior probability = 1.00). The first subclade included genera Gymnogeophagus, Geophagus, Mikrogeophagus, Dicrossus, Biotodoma and Crenicichla (from tribe Cichlinae) (bootstrap support = 96%; Bayesian posterior probability = 1.00). The second subclade consisted of genera Taeniacara and Apistogramma (bootstrap support = 100%; Bayesian posterior probability = 1.00). With the exception of H. severus, the phylogenetic analysis also showed that the tribe Cichlasomatinae was monophyletic. All six examined genera within the tribe Cichlasomatinae (Amphilophus, Andinoacara, Bujurquina, Pterophyllum, Symphysodon and Thorichthys) containing two or more species were reciprocally monophyletic. Notably, genetic similarity exceeded 99.9% between Amphilophus amarillo and A. citrinellus, Symphysodon discus and Symphysodon haraldi, and Cichla monoculus and C. ocellaris. Phylogenetic relationships of Cichlidae were reconstructed after removing Crenicichla regani and H. severus. No significant changes were observed in the position of the target species or other taxa in the phylogenetic tree compared to the original topology (see Supplementary Material, Supplementary Figure S1, S2).