1 INTRODUCTION

Cavefish are unique hypogean species whose evolutionary histories are associated with caves and other subterranean waters (Ma et al., 2019). There are two kinds of cavefish: stygobionts (also called typical or troglomorphic cavefishes) and stygophiles (also called atypical or non-troglomorphic cavefishes) (Zhao et al., 2011). Stygobionts are obligate cavefish whose adaptive traits, such as the degeneration or loss of eyes and reduced pigmentation, both enable and restrict them to complete their entire life cycles within cave environments. Stygophiles are not limited to hypogean aquatic environments but rely heavily on cave habitats during their life cycles (Ma et al., 2019).

Globally, more than 350 species of cavefish live in caves and other subterranean habitats, including at least 199 species of stygobionts (Niemiller et al., 2019). Most cavefish species are found in the subtropics and tropics, and 75% of known populations are located in Southeast Asia or Central and South America within karst landscapes (Ma et al., 2019).

There are 148 described species of cavefish in China (accounting for 39% of such species worldwide), most of which reside within the 620 000 km² of karst landscape in the southwest (Huang et al., 2008). Karst landscapes feature complex surfaces and underground systems and provide diverse habitat to aquatic and terrestrial species (Deharveng et al., 2009; Luo et al., 2016). However, most species are limited to one independent cave or subterranean stream (Shu et al., 2013). For example, Heminoemacheilus hyalinus is only distributed in a single cave in Du'an in the Guangxi Zhuang Autonomous Region. Karst habitats comprise variable and unique geological structures that can harbor communities of highly endemic and endangered species (Luo et al., 2016); these habitats are extremely fragile and sensitive to environmental change (Zhao et al., 2011; Palandačić et al., 2012). Caves are natural laboratories for ecological and evolutionary research, and the ability of cavefish to survive in extreme environments makes them excellent models for adaptation studies (Riddle et al., 2018). However, difficulties in accessing and sampling remote caves have limited study of these ecosystems (Shu et al., 2013). As such, many cave species remain poorly studied, particularly with regard to their genome-level characteristics.

Cave environments are characterized by permanent darkness. The lack of autotrophic primary producers in dark environments can lead to restrictions in food supply; hence, cavefish have likely developed strategies to enhance their energy utilization to survive. Mitochondria play vital roles in energy metabolism. Cavefish mitochondrial DNA (mtDNA) may have undergone positive directional selection in response to the extreme environmental pressures in cave habitats (Boore, 1999; Xu et al., 2007; Tomasco and Lessa, 2011). Generally, mtDNA is characterized by its simple structure, small size, maternal inheritance, and rapid evolution (Avise et al., 1984) and has been extensively studied in relation to taxonomy, genetics, phylogenetics, and evolution (Satoh et al., 2016; Liu et al., 2017). However, few fully sequenced cavefish mitochondrial genomes (mitogenomes) are available with which to study adaptive evolution to subterranean life.

Nemacheilidae is the most diverse cavefish family (53 species) (Niemiller et al., 2019), which may have arisen due to the strong selective pressures found in isolated cave habitats (Zhao and Zhang, 2009). In southwest China, nemacheilid cavefish exhibit rich species diversity and some species, such as Paranemachilus pingguoensis, have relatively large population sizes (Lan et al., 2013). Nemacheilid fish can be valuable for systematic study of cavefish evolution and adaptation and can serve as references for the study of karst development in the Yunnan-Guizhou area. However, despite their great scientific value, Nemacheilidae cavefish have been much less studied than their surface-dwelling counterparts, with very few studies to date (Lan et al., 2013).

In the current research, the complete mitogenomes of three nemacheilid cavefish species, i.e., Paranemachilus pingguoensis, Oreonectes polystigmus, and Heminoemacheilus longibarbatus, were sequenced. In morphology, O. polystigmus can be distinguished from P. pingguoensis and H. longibarbatus by caudal fin rounded and tip of anterior nostril elongated and barbel-like. H. longibarbatus and P. pingguoensis share closely set nostrils but can be separated from each other by cheek scales present in Paranemachilus vs. absent in Heminoemacheilus. The genomic characteristics of these cavefish were compared with those of closely related species within Nemacheilidae to investigate genetic adaptations to cave habitats. The evolutionary history of the nemacheilid cavefish was explored and found to be linked to geological events that altered the karst environments.

2 MATERIAL AND METHOD 2.1 Sample collection

Paranemachilus pingguoensis (tag number: PG-20190701; GPS: 23.33°N, 107.32°E), Oreonectes polystigmus (tag number: DB-20190701; GPS: 25.03°N, 110.31°E), and Heminoemacheilus longibarbatus (tag number: CX-20190701; GPS: 24.11°N, 108.10°E) specimens used in this study were collected in Guangxi Zhuang Autonomous Region, China in July 2019 using a cage net (mesh size: 0.70 cm, length: 20 m, width: 45 cm, height: 33 cm). After collection, specimens were euthanized in tricaine mesylate (MS-222), and then dissected. Muscle tissues were preserved in 95% ethanol in the field and transferred to a -80-℃ freezer for long-term storage in a laboratory at the College of Marine Sciences, South China Agricultural University (SCAU), China. The experiments were approved by the Animal Ethics Committee of SCAU (SYXK-2019-0136).

2.2 DNA extraction and mitogenome sequencing

Genomic DNA was extracted from the muscle tissue samples using an EZNA^® Tissue DNA kit (Omega Bio-tek, USA) following the manufacturer's instructions. The DNA was evaluated and quantified using 1.0% gel electrophoresis followed by ultraviolet spectroscopy (NanoDrop 2000, Thermo Fisher Scientific, USA). The resulting high-quality genomic DNA was isolated and used for library construction (TruseqTM RNA Sample Prep Kit, Illumina, USA), with a 350-bp insert used for paired-end (150 bp) sequencing on the HiSeq4000 platform (Illumina). The mitogenomes then underwent next-generation sequencing following the manufacturer's instructions.

2.3 Analysis of sequence data

Raw sequencing data were retrieved and trimmed using Cutadapt v1.16 software with the parameters -q 20 -m 20 (http://cutadapt.readthedocs.io/) to obtain clean data. The clean paired-end reads were then used for de novo assembly using NOVOPlasty v3.8 software (Dierckxsens et al., 2017).

The mitogenomes were annotated using Mitofish v3.51 (Iwasaki et al., 2013) and ORF Finder (https://www.ncbi.nlm.nih.gov/orfinder/). The preliminary results were compared with the coding proteins and rRNA of mitogenomes from related species, using Blastp and Blastn (https://blast.ncbi.nlm.nih.gov/Blast.cgi), respectively, to confirm and correct the accuracy of the results. The complete mitochondrial DNA sequences were deposited in GenBank under accession numbers MW043720 (P. pingguoensis), MW043721 (H. longibarbatus), and MW043722 (O. polystigmus).

Mitochondrial codon usage biases can affect gene expression and reflect the evolutionary relationships among species (Wei et al., 2014). Therefore, the relative synonymous codon usage (RSCU) for each genome was analysed (Sharp and Li, 1987). Repeat sequences were isolated using Tandem Repeats Finder v4.04. The skew associated with nucleotide bases was calculated using the following formulae: adenine (A) and thymine (T) skew=(A T)/(A+T) and guanine (G) and cytosine (C) skew=(G C)/(G+C). Non-synonymous/synonymous mutation (Ka/Ks) ratios were calculated and applied to investigate the influence of selective pressures on protein-coding genes (PCGs) in cavefishes and related species; this was performed using the seqinr package in R software (Charif and Lobry, 2007).

2.4 Phylogenetic analysis

The mitogenome sequences of 61 nemacheilid species were subjected to phylogenetic analyses, with Leptobotia pellegrini (NC_031602) rooted as the outgroup. We aligned the individual genes for all species using MUSCLE v.3.8.31 (Edgar, 2004) according to the selected coding sequence. Following alignment, the genes of each species were then connected into a series that was integrated into a corresponding PCG sequence set for the following analyses.

Phylogenetic analyses were conducted using two approaches: maximum likelihood (ML) and Bayesian inference (BI). The DNA models were tested using jModelTest (Posada, 2008) based on selected sequences; the best model was selected based on the minimum Akaike information criterion. For the ML phylogenetic analysis, RAxML v8.1.5 software (Stamatakis, 2014) was applied using 1 000 bootstrap replicates. The number of replicates was determined automatically using the protein substitution model with the highest score. For the BI phylogenetic analysis, MrBayes v3.2.6 (Ronquist et al., 2012) was used, and results were considered reliable if the split frequency standard deviations were < 0.01.

A molecular clock was constructed using MEGA X v10.1.7 software (Kumar et al., 2018) following phylogenetic and corresponding sequence analyses using ML. The differentiation times among species were queried using TimeTree (Kumar et al., 2017) with the following reference differentiation times: 11.70million years ago (Ma) 45.42 Ma for Triplophysa stenura and Schistura fasciolata, 1.50 Ma 11.69 Ma for Lefua nikkonishe and Lefua costata and 48 Ma 103 Ma for Oreonectes platycephalus and Triplophysa stenura.

3 RESULT 3.1 Mitogenome organization and composition

The complete mitogenomes of Oreonectes polystigmus, Paranemachilus pingguoensis, and H. longibarbatus were 16 568, 16 562, and 16 569 bp in length, respectively. They each contained the canonical 13 PCGs, two rRNA genes (12S and 16S), 22 tRNA genes, and a D-loop control region. The majority of the genes are encoded by the heavy strand (H-strand), except for nicotinamide adenine dinucleotide hydride (NADH) dehydrogenase subunit 6 (ND6) and eight tRNA genes for glutamine, alanine, asparagine, cysteine, tyrosine, serine, glutamate, and proline, which were all located on the light strand (L-strand) (Fig. 1 and Table 1). The two rRNAs, 16S and 12S, were located between tRNA-Leu (TAA) and tRNA-Val (TAC) and between tRNA-Val (TAC) and tRNA-Phe (CAA), respectively (Fig. 1).

3.2 Base composition and base bias

The nucleotide composition (H-strand) of the three cavefish in descending order was 30.13% A, 27.69% C, 25.75% T, and 16.43% G, with an AT content of 55.88% (54.91%–56.44%).

Most PCGs had a negative GC skew, except for ND6, which exhibited a positive GC skew (Fig. 2). Furthermore, the AT skews of most PCGs were around zero, except for ND6, which had a negative AT skew (-0.48, -0.53, and -0.50 for Oreonectes polystigmus, Paranemachiluspingguoensis, and Heminoemacheilus longibarbatus, respectively).

3.3 PCGs and codon usage

The coding regions of the three species were 11 427 bp (3 800 amino acids) in length, accounting for approximately 69% of the total mitogenome. The initiation codon for most PCGs was the canonical methionine sequence (ATG), with the exception of the start codon for COXI, which was GTG. We also identified T, TA, TAA, and TAG termination codons in the cavefish mitogenomes, with TAA and T being the most prevalent. The TAG termination codon was only found in the ND1, ND5, and ND6 genes of Heminoemacheilus longibarbatus (Table 1).

For the three examined species, all amino acids were encoded by two or four codon variations, except for leucine and serine, which were encoded by six (Fig. 3). RSCU analysis revealed that the following codon sequences had higher frequencies: stop (TAA), arginine (CGA), and leucine (CTA). According to RSCU values, codons with A or C in the third position occurred more frequently than those with T or G (Fig. 3).

The PCG Ka/Ks ratios were all < 1 (Fig. 4). The Ka/Ks ratios for ATP8 and ND2 were the highest (Fig. 4).

3.4 Genes encoding tRNAs and rRNAs

In total, 22 tRNA genes were identified in the three cavefishes, ranging from 66 bp (tRNA-Cys) to 75 bp (tRNA-Leu and tRNA-Lys) in length. In general, codons and anticodons exhibited a one-to-one relationship. However, two anticodons were present for both serine (TGA and GCT) and leucine (TAG and TAA). The dihydrouridine (DHU) arm of tRNA-Ser (GCT) was absent, whereas the secondary structures of the other 21 tRNAs exhibited typical clover structures.

3.5 Compositional biases between cavefishes and other nemacheilid species

The genomic sequences of other cavefishes (stygobionts and stygophiles) and other nemacheilid species were compared to investigate differences in nucleotide composition. Absolute value of GC and AT skews were higher in the cavefishes than in other nemacheilid species (Fig. 5); this was the case for the entire mitogenome, as well as for individual elements such as PCGs, rRNAs, and tRNAs. Principal component analyses (PCA) suggested that the overall codon and amino acid usage patterns of cavefishes were notably different from those of other fish (Fig. 5e & f).

3.6 Phylogenetic and divergence time analyses of cavefishes

A phylogenetic tree was constructed based on the mitogenomic PCGs of 61 nemacheilid species, including three species of stygobionts and five species of stygophiles, with Leptobotia pellegrini as the outgroup (Fig. 6). Both ML and BI approaches produced trees with similar topologies and branch lengths, indicating a well-resolved and highly supported phylogeny of Nemacheilidae. The results show that the 61 nemacheilid species could be divided into four branches or clades. Clade Ⅰ included all species of Triplophysa and Barbatula. Clade Ⅱ was sister to Clade Ⅰ and comprised Homatula, Schistura, and Tuberoschistura. Clade Ⅲ included the genera Heminoeacheilus, Micronemacheilus, Oreonectes, Troglonectes, Paranemachilus, Lefua, and Traccatichthys. The other eight genera, Petruichthys, Schistura, Nemaceilus, Oxynoemacheilus, Nemachilichthys, Hedinichthys, Acanthocobitis, and Aborichthys, were clustered into Clade Ⅳ, which was positioned at the base of the Clade Ⅲ + (Clade Ⅱ + Clade Ⅰ) clade.

The cavefishes were distributed across five genera (Triplophysa, Heminoemacheilus, Oreonectes, Troglonectes, and Paranemachilus) in Clade Ⅰ and Clade Ⅲ. In Clade Ⅰ, the genus Triplophysa appeared to diverge from the other species around 24.83 Ma (Fig. 6). In Triplohysa, the divergence of two stygobionts (Triplophysa xiangxiensis and Triplophysa rosa) from other non-cave species was estimated to have occurred around 20.93 Ma during the Miocene epoch. In Clade Ⅲ, the cavefishes (Paranemachilus jinxiensis, Paranemachilus pingguoensis, Troglonectes daqikongensis, and Troglonectes furcocaudalis) diverged from their sister group Lefua around 24.65 Ma during the Chattian age of the late Oligocene epoch. The stygophile Oreonectes polystigmus diverged from its surface counterpart (Oreonectes platycephalus) circa 6.40 Ma during the Messinian age of the late Miocene epoch, and the stygophile Heminoemacheilus longibarbatus diverged from its surface counterpart (Y. cruciatus) circa 18.91 Ma during the Burdigalian age of the Miocene epoch. In terms of divergence times among cavefishes, the occurrence range was 7.71 Ma (Triplophysa xiangxiensis and Triplophysa rosa) to 13.76 Ma (Troglonectes daqikongenis and Troglonectes furcocaudalis) (Serravallian to Messinian Age, Miocene Epoch).

4 DISCUSSION

Previous research has shown that certain teleostean fish possess only one or similar arrangements of mitochondrial genes (Miya et al., 2003). We found similar gene arrangements in the mitogenomes of the three cavefish species, which are also teleostean fish. Gene arrangement and distribution without subsequent gene rearrangement is a common feature in vertebrates (Broughton et al., 2001). The genomes of all examined nemacheilid species exhibited higher AT content relative to GC, consistent with the high AT and low GC contents displayed in most vertebrates (Cheng et al., 2012; Luo et al., 2017; Kundu et al., 2018; Huang et al., 2019). However, negative AT skews were found in ND6 in the three cavefish, which are uncommon in mitochondrial genes, but have been detected in other vertebrates such as Marmota himalayana (Chao et al., 2014), Cetonurus globiceps (Shi et al., 2016), and Ophichthus brevicaudatus (Lü et al., 2019). These animals inhabit extreme environments, a feature shared by nemacheilid cavefish. The nucleotide skews may have been the result of mutational and selective pressures during DNA replication and transcription (Perna and Kocher, 1995; Wei et al., 2010a, b). Thus, the unusual composition of ND6 in the cavefish may be related to selective pressure from their extreme habitat. Nevertheless, further investigation is required to understand how nucleotide skews influence replication direction or gene orientation, and what selective pressures cause skews to arise.

The PCG Ka/Ks ratios were all < 1 (Fig. 4), suggesting that the PCGs evolved under purifying selection, during which harmful mutations with negative selection coefficients were naturally selected against (Yang and Bielawski, 2000). Under purifying selection, it is difficult for disadvantageous mutations to be established, and the long-term stability of biological structures is maintained. Thus, low Ka/Ks ratios for COXI, COXIII, and CYTB indicate strong purifying selection and evolutionary conservation in these genes. Differences in Ka/Ks ratios can indicate differences in selection pressures and evolutionary mechanisms among genes (Lynch et al., 2006). The Ka/Ks ratios for ATP8 and ND2 were highest (Fig. 4), suggesting that ATP8 and ND2 evolved at a faster rate than the other mitochondrial PCGs in the nemacheilid cavefish. More functional genes are present in the H-strand of the mitogenome; for example, ATP8 and ND2 both reside on the H-strand of the mitogenome, implying that selective pressures may have exerted stronger influence on the heavy chain than on the light chain in nemacheilid cavefish mtDNA.

Based on principal component analysis (PCA), the overall codon and amino acid usage patterns of the cavefish differed significantly from those of other fish (Fig. 5e & f), further supporting the nucleotide composition and skew results. Mitochondria are the centres of energy metabolism in eukaryotic organisms and produce more than 95% of all cellular energy via oxidative phosphorylation (Xu et al., 2007). In cave habitats, the lack of autotrophic primary producers due to constant darkness can lead to limitations in food supply (Zhao and Zhang, 2009). The differences in nucleotide compositional biases between cavefish and other nemacheilid species suggest that mtDNA may have been subjected to evolutionary selection to improve energy use efficiency in response to food limitations in cave environments; however, further investigations are required for this to be confirmed.

To explore the evolutionary history between cavefish and their surface-dwelling relatives, we constructed a phylogenetic tree based on the mitogenome PCGs of 61 nemacheilid species, including three species of stygobionts and five species of stygophiles, with Leptobotia pellegrini as the outgroup (Fig. 6). Both maximum-likelihood (ML) and Bayesian inference (BI) approaches produced trees with similar topologies and branch lengths, indicating a well-resolved and highly supported phylogeny of Nemacheilidae. The data revealed that Heminoemacheilus longibarbatus and Micronemacheilus cruciatus clustered into one separate group, while Paranemachilus jinxiensis grouped with Paranemachilus pingguoensis. This corroborates the findings of Du et al. (2021) that Paranemachilus jinxiensis (formerly Yunnanilus jinxiensis) belongs to the genus Paranemachilus. Our results also support Troglonectes furcocaudalis and Troglonectes daqikongensis (formerly Oreonectes furcocaudalis and Oreonectes daqikongensis, respectively) as species of Troglonectes (Huang et al., 2020). Convergent evolution is a common phenomenon in extreme cave environments, and the present study challenges the results of previous morphological phylogenetic analyses (Wilkens and Strecker, 2003; Protas et al., 2006; Trontelj, 2019). Therefore, comprehensive taxonomic revisions of these species, using both morphological and genomic data, are needed.

The separation of cave and non-cave species (24.83 Ma–24.65 Ma–20.93 Ma–6.40 Ma) occurred concurrently with the early stages of the Miocene uplift of the southern Qinghai-Tibet Plateau (QTP), which was characterized by intensive intraplate tectonism, magmatism, and metallogeny (Li, 2010). Differentiation among cavefish occurred during the late stage of the Miocene uplift (13.76 Ma 9.08 Ma– 7.71 Ma), by which time the tectonism and magmatism that had defined the earlier stages were less prominent (Li, 2010). These splitting events coincided with the intraplate orogeny of the southern QTP, which mainly occurred 25 Ma 7 Ma (Li, 2010). Throughout the Miocene uplift of the southern QTP, large-scale karst development occurred in the limestone regions (Shi et al., 1999). Uplift of the QTP significantly affected both the climate and habitat of the area (An et al., 2001; Xiao et al., 2005; Zhao and Zhang, 2009). Thus, the isolation of cave habitats, along with the unique environmental effects of the QTP uplift process, are considered key reasons for cavefish speciation (Zhao and Zhang, 2009; Yang et al., 2016). Our results provide further evidence that the speciation and diversification of nemacheilid cavefish are correlated with uplift of the QTP.

5 CONCLUSION

The complete mitogenomes of three cavefish (Order Cypriniformes, Family Nemacheilidae) were sequenced, annotated, and characterized. The mitogenome contained 37 genes (13 PCGs, 22 transfer RNAs (tRNAs), and two ribosomal RNAs (rRNAs)) and a single control region. The Ka/Ks ratios indicate that the PCGs of the three cavefish evolved under purifying selection. The low Ka/Ks ratios for COXI, COXIII, and CYTB indicate strong purifying selection and evolutionary conservation in these genes, whereas the high Ka/Ks ratios for ATP8 and ND2 indicate that these genes evolved at a faster rate than the other mitochondrial PCGs in nemacheilid cavefish. PCA revealed that the mitogenomes of the three nemacheilids exhibited nucleotide biases within the PCGs, tRNAs, rRNAs, and full genomes. These results imply that the mtDNA underwent evolutionary selection to improve energy use efficiency in response to food limitations in cave environments. Divergence time and evolutionary history analyses indicated that speciation and diversification of the cavefish coincided with the Miocene uplift of the southern QTP, when the cave habitats were undergoing significant change. This study has helped to elucidate the mitogenomes and phylogeny of fish in Nemacheilidae and provided insights into the evolutionary history of cavefish.

6 DATA AVAILABILITY STATEMENT

The complete mitochondrial DNA sequences were deposited in GenBank under accession numbers MW043720 (P. pingguoensis), MW043721 (H. longibarbatus), and MW043722 (O. polystigmus). The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

7 ACKNOWLEDGMENT

We are very grateful to Jiahu LAN from Du'an Fishery Technique Popularization Station for specimen collection and photographing on the three species.