All methods were carried out in accordance with relevant guidelines and regulations, and the study was carried out in compliance with the ARRIVE guidelines. Sampling of the O. infernale individual used to generate the mitogenome presented here was accomplished with the assistance of a professional cave diver who captured the specimen using a custom-made hand net specifically designed for efficient capture and secure storage while cave diving. The sample was collected under collecting permit SGPA/DGVS/05375/19 issued by the Mexican Ministry of Environment and Natural Resources (Secretaría de Medio Ambiente y Recursos Naturales; SEMARNAT) to JA. The sampling locality is the cenote Kan-Chin (Huhí, Yucatán), located at 20°40'11"N, 89°10'6"W. The voucher specimen was euthanized with MS-222 prior to preservation in accordance with recommended guidelines for the use of fishes in research (Nickum et al. 2004), fixed in a 10% formalin solution, and subsequently transferred to 70% ethanol for long-term storage in the Colección Nacional de Peces (CNPE) of the Instituto de Biología (IB) at the Universidad Nacional Autónoma de México (UNAM), where it has been catalogued and deposited (CNPEIBUNAM 23285). A tissue sample (muscle fragment) was taken prior to specimen fixation, preserved in 95% ethanol, and eventually cryopreserved at -80 °C. High-molecular genomic DNA was extracted using the phenol-chloroform protocol (Sambrook et al. 1989). The DNA was sheared by sonication with a Bioruptor pico of Diagenode and Minichiller. Sonication was performed using six cycles of alternating 30 s ultrasonic bursts and 30 s pauses in a 4 °C water bath. For library preparation we used a DNA sample of 200 ng which was quantified using a Qubit fluorometer (Invitrogen). Library preparation was carried out using the KAPA Biosystem Hyper Kit (Kapa, Biosystem Inc., Wilmington, MA). Fragmented DNA was ligated to custom, TruSeq-style dual-indexing adapters (Glenn et al. 2016). Fragments were size selected in a ~300–500 bp range which was enriched through PCR, purified and normalized. The Illumina NextSeq v2 300 cycle kit was used for sequencing paired-end 150 nucleotide reads at the Georgia Genomics Facility, University of Georgia, Athens, USA.

The quality of the raw data was assessed with FastQC (Andrews, 2010). Good-quality sequences that did not contain ambiguous nucleotides and reads with average quality of 30Q were demultiplexed, trimmed and merged using Geneious Prime 2020.0.4 (https://www.geneious.com). Mitogenome assembly was conducted with MITObim v.1.9 (Hahn et al. 2013) using two reference mitogenomes from close relatives of O. infernale available in GenBank: Ophisternon candidum (MT436449) and Synbranchus marmoratus (AP004439). These reference mitogenomes were aligned in order to generate a consensus sequence for use during the annotation procedure. MitoFish MitoAnnotator (Iwasaki et al. 2013) and MITOS (Bernt et al. 2013) were used to identify and annotate protein-coding genes (PCGs), transfer RNAs (tRNAs), and ribosomal RNAs (rRNAs). The resultant annotated O. infernale mitochondrial genome was deposited in the GenBank database under accession number OM388306.

Nucleotide and amino acid composition, codon usage profiles of protein-coding genes (PCGs), Relative Synonymous Codon Usage (RSCU), and characterization the non-coding mtDNA control region (CR) were computed with mega X (Kumar et al. 2018). Nucleotide composition skewness was calculated with the formulas AT skew = (A − T)/(A + T) and GC skew = (G − C)/(G + C) (Perna and Kocher 1995). Prediction of tRNAs secondary structure was accomplished with tRNAScan-SE 2.0 (Chan and Lowe 2019) through the webserver http://lowelab.ucsc.edu/tRNAscan-SE/, using Infernal without HMM filter search mode and “vertebrate mitochondrial” as sequence source (Lowe and Chan 2016). Analysis and prediction of CR secondary structure in O. infernale was accomplished using the software ClustalW (Thompson et al. 2003) as implemented in mega X (Kumar et al. 2018) by comparison (via multiple sequence alignment) with reports of secondary CR structure from two other teleost fishes, namely Siniperca chuatsi (EU659698) (Zhao et al. 2006) and Cyprinion semiplotum (MN603795) (Sharma et al. 2020).

We investigated patterns of selection on PCGs on a mitogenomic scale and phylogenetic relationships among major synbranchiform lineages based on all mitogenomic comparative data for the group available on GenBank. To measure of the strength and mode of natural selection acting on PCGs, we estimated the ratio of non-synonymous (K_A) to synonymous (K_S) substitutions (K_A/K_S, also known as ω or d_N/d_S) using the HyPhy 2.5 package (Kosakovsky Pond et al. 2020) as implemented in mega X (Kumar et al. 2018) based on the newly assembled mitochondrial genome (Ophisternon infernale, OM388306) and seven additional synbranchiform mitogenomes previously available in GenBank: Ophisternon candidum (MT436449), Synbranchus marmoratus (AP004439), Monopterus albus (NC003192), Mastacembelus armatus (NC023977), Mastacembelus erythrotaenia (NC035141), Macrognathus aculeatus (KT443991), and Macrognathus pancalus (NC032080). To compare patterns of selection between synbranchiform families, we conducted two separate K_A/K_S analyses, one for synbranchids and one for mastacembelids. The taxonomic sampling for phylogenetic analyses included representatives of the synbranchiform families Synbranchidae and Mastacembelidae, as well as a representative of Indostomidae, a monogeneric family historically classified in the Gasterosteiformes on the basis of morphological evidence (Van Der Laan et al. 2014; Nelson et al. 2016) but more recently assigned to the Synbranchiformes in accordance to the results of molecular phylogenetic studies (Betancur-R. et al. 2013; Betancur-R et al. 2017). The lack of published mitochondrial genomes of fishes from the synbranchiform family Chaudhuriidae prevented us from including representatives of this taxon in our analyses. The ingroup consisted of the synbranchids Ophisternon infernale (OM388306), Ophisternon candidum (MT436449), Synbranchus marmoratus (AP004439) and Monopterus albus (NC003192), the mastacembelids Mastacembelus armatus (NC023977), Mastacembelus erythrotaenia (NC035141), Macrognathus aculeatus (KT443991) and Macrognathus pancalus (NC032080), and the indostomid Indostomus paradoxus (NC004401). The outgroup consisted of representatives of close relatives of Synbranchiformes such as the anabantiforms Channa micropeltes (NC030542) and Nandus nandus (AP006809), and the gasterosteiforms Gasterosteus aculeatus (NC041244) and Pungitius pungitius (NC011571); the last two included to test the phylogenetic position of I. paradoxus with respect to members of the Gasterosteiformes. The phylogeny was rooted at the viviparous brotula Diplacanthopoma brachysoma (AP004408). Phylogenetic relationships were inferred based on a concatenated alignment of all 13 PCGs. DNA sequence data from each PCG was independently aligned via multiple sequence alignment using the software MUSCLE (Edgar 2004) under default parameters via the “translation align” tool of the software Geneious Prime 2020.0.4 (https://www.geneious.com). The best-fit substitution model for each PCG was determined according to the corrected Akaike Information Criterion (AICc) with the software jModelTest2 (v. 2.1.10) (Darriba et al. 2012) under the following likelihood settings: number of substitution schemes = “3”; base frequencies = “+F”; rate variation = “+I and + G with nCat = 4”; base tree for likelihood calculations = “ML optimized”; and base tree search = “Best” (effectively evaluating among all 24 “classical” GTR-derived models). Individual alignments (ATP6=681 bp, ATP8=168 bp, COX1=1,539 bp, COX2=690 bp, COX3=783 bp, CYTB=1,137 bp, NAD1=975 bp, NAD2=1,053 bp, NAD3=348 bp, NAD4=1,380 bp, NAD4L=294 bp, NAD5=1,836 bp, and NAD6=525 bp) were subsequently concatenated using the software 2matrix (Salinas and Little 2014), yielding a data matrix totaling 11,409 aligned bp. Maximum Likelihood inference of phylogeny was carried out on the concatenated alignment partitioned by gene using the software RAxML-NG (v. 1.0.1) (Kozlov et al. 2019) through the CIPRES Science Gateway (Miller et al. 2010), with nodal support estimated by means of the bootstrap character resampling method (Felsenstein 1985) based on 1000 pseudoreplicates.

The complete mitochondrial genome of O. infernale presented herein (GenBank accession number OM388306) is 16,804 bp in total length (Fig. 1; Table 1), a somewhat larger size than previously published synbranchiform mitogenomes, which range from 16,493 bp (in M. erythrotaenia) (or from 16,152 bp if considering the putative synbranchiform I. paradoxus) to 16,622 bp (in M. albus). Although the mitogenome of the synbranchid S. marmoratus reported in GenBank (AP004439) is considerably shorter (15,561 bp), this significant difference in length is actually due to it missing the NAD1 gene (normally ~1,000 bp) as a result of reported technical difficulties during sequencing (Miya et al. 2003). The composition and general arrangement of mitochondrial genes in O. infernale is identical to that reported for other synbranchiforms (Li et al. 2016; Han et al. 2018; White et al. 2020) as well as for more distantly related teleosts (Miya et al. 2001, 2003; Satoh et al. 2016), and consists of a total of 37 genes divided into the following categories: 13 PCGs, 2 rRNAs, 22 tRNAs, and the non-coding control region (CR) (Fig. 1; Table 1). Twenty-eight genes (12 PCGs, 2 rRNAs, 14 tRNAs) plus CR are located on the H-strand, while the remaining nine genes (NAD6 and 8 tRNAs) are located on the L-strand (Table 1); a configuration that corresponds to those of previously reported synbrachiform mitogenomes (Li et al. 2016; Han et al. 2018; White et al. 2020). The overall base composition of the O. infernale mitogenome is T=28.7%, A=31.6%, G=13.2%, and C=26.5%, which is fairly similar to those of other synbrachiform mitogenomes (Table 2). Nucleotide composition, however, is biased toward A+T (60.4%), with O. infernale displaying the highest values of this metric among the analyzed synbranchiforms. The mitogenome of O. infernale exhibits positive AT (0.046) and negative GC (-0.277) skewness, a general pattern shared with other species of the Synbranchiformes (Table 2).

Figure 1. 

Annotated map of the mitochondrial circular genome of O. infernale. The outer ring corresponds to the H- (outermost) and L-strands, and depicts the location of PCGs (in black, except for ND6 which is encoded in the L-strand and is portrayed in red), the non-coding control region (in dark brown), tRNAs (in red), and rRNAs (in light brown). The inner ring (black sliding window) denotes GC content along the genome. Live specimen photograph taken in the Cenote Kancabchen (Homún, Yucatán), courtesy of cave diver Erick Sosa.

The 13 PCGs, altogether totaling 11,449 bp, correspond to 68.1% of the O. infernale mitogenome. These genes consist of seven regions that code for the subunits of the NADH dehydrogenase (ubiquinone) protein complex (NAD1-6, NADL4), three that code for the subunits of the enzyme cytochrome c oxidase (COX1-3), one that codes for the enzyme cytochrome b (CYTB), and two that code for the subunits 6 and 8 of the enzyme ATP synthase F_O (ATP6, ATP8). Except for COX1 and ND4L, PCGs exhibit an ATG (Met) start codon­, which is the standard in eukaryotic systems (Kozak 1983). The start codon exhibited by COX1 (GTG), however, is fairly common among vertebrates (Nwobodo et al. 2019). Conversely, an initiation-codon change from ATG (Met) to ATA (Ile) such as the one observed in ND4L is less common. Notably, of the synbranchiform mitogenomes analyzed, only that of O. infernale displays ATA as ND4L initiation codon. Most PCGs (10 out of 13) exhibit a TAA stop codon, which is a standard termination codon common in vertebrate mtDNA. However, of these 10 genes only three (NAD1, NAD4L, ATP8) display a complete codon (TAA), while the remaining seven (ATP6, COX2, COX3, NAD4, NAD5, NAD6, CYTB) contain an incomplete stop codon (either TA or T). Of the remaining three PCGs, NAD2 and COX1 have the stop codon AGA, while NAD3 has the stop codon GAC (Table 2). PCGs in the mitogenome of O. infernale exhibit levels of A+T content (60.1%) comparable to–though slightly higher than–those of other synbranchiforms, which range from 53.1% in M. armatus to 59.1% in O. candidum (Table 2). In contrast to our findings for the entire mitogenome, AT-skews in PGCs across all synbranchiform mitogenomes analyzed exhibit negative values. Conversely, and in correspondence with our whole-mitogenome results, GC-skews in PGCs also exhibit negative values and highly similar across most analyzed synbranchiforms. A total of 3816 amino acids are encoded by PCGs in the mitogenome of O. infernale, with Leu (14.7%), Ser (9.4%), Thr (7.8%), and Pro (7.7%) being the most frequent, while Met (1.1%) being the least common. RSCU values represent the ratio between the observed usage frequency of one codon in a gene sample and the expected usage frequency in the synonymous codon family, given that all codons for the particular amino acid are used equally. The synonymous codons with RSCU values > 1.0 have positive codon usage bias and are defined as abundant codons, whereas those with RSCU values < 1.0 have negative codon usage bias and are defined as less-abundant codons (Gun et al. 2018). Results from RSCU analysis of PCGs in the mitogenome of O. infernale indicate that the most frequently used codons are ACC (1.59%), AAA (1.56%), TTA, ATA, and GAA (1.49%), which code for the amino acids Thr, Lys, Leu, Met, and Glu, respectively. On the other hand, codons encoding Prol (CCG, 0.16%), Thr (ACG, 0.2%), Ala (GCG, 0.23%), Ser (TCG, 0.39%), and Leu (CTG, 0.4%; TTA, 0.48%) are the least frequent (Fig. 2; Table 3).

Figure 2. 

Results from analysis of Relative Synonymous Codon Usage (RSCU) of the mitochondrial genome of O. infernale. Codon families are plotted on the x-axis. The label for the 2, 4, or 6 codons that compose each family is shown in the boxes below the x-axis, and the colors correspond to those in the stacked columns. RSCU values are shown on the y-axis.

The mitogenome of O. infernale contains the typical 22 tRNAs usually documented for mitogenomes of other teleosts and vertebrates (Lee et al. 1995; Díaz-Jaimes et al. 2016; Satoh et al. 2016; Nwobodo et al. 2019; White et al. 2020). The genomic organization of tRNAs in O. infernale is identical to that reported for O. candidum (White et al. 2020) and other synbranchids (Li et al. 2016; Han et al. 2018). Altogether, tRNAs total 1547 bp, with individual ones ranging from 64 bp (tRNA^Ser) to 76 bp (tRNA^Thr) (Table 1). Fourteen tRNAs are encoded in the H-strand, while the remaining eight in the L-strand (Fig. 1; Table 1). Twenty-one of the 22 tRNAs fold into the canonical cloverleaf secondary structure that consists of four domains (AA stem, D arm, AC arm, and T arm) and a variable loop (Fig. 3). Notably, the tRNA^Ser (11889–11952) exhibits an unusual structure in which the D arm is missing. Although any change in tRNA secondary structure could potentially alter its amino acid recognition capability (Nwobodo et al. 2019), it has been shown that loss of the D arm does not necessarily imply reduced functionality; in fact, almost all tRNAs^Ser for AGY/N codons lack the D arm, and truncated tRNAs appear to have been compensated for by several interacting factors (Watanabe et al. 2014). Furthermore, among fishes, loss of the tRNA^Ser D arm is not unique to O. infernale, for it has been reported in several species, including chondrichthyans such as Chiloscyllium griseum (Chen et al. 2013), Triaenodon obesus (Chen et al. 2016), and Cephalloscyllium umbratile (Zhu et al. 2017), as well as teleosts such as Oreochromis andersonii and O. macrochir (Bbole et al. 2018). Although most tRNAs present the canonical 7-bp T loop, nonstandard T-loop lengths were observed in tRNA^Met (6 bp), tRNA^Phe (8 bp), and tRNA^Ser (9 bp). Other deviations from the traditional tRNA secondary structure that could affect functionality is the presence of extra loops. The tRNA^Arg (10080–10148) in the mitogenome of O. infernale exhibits a loop at the base of the AA stem, thus potentially affecting aminoacylation. The nucleotide composition in the tRNAs of the O. infernale mitogenome is T=29%, A=31.4%, G=20.2%, and C=19.3%. The genes that code for the mitochondrial 12S and 16S rRNA subunits in O. infernale are 948 bp and 1092 bp long, respectively, and are located on the H-strand separated by the tRNA^Val, just like in most teleost fishes (Lee et al. 1995; Satoh et al. 2016).

Figure 3. 

Secondary structure of the 22 tRNA genes of the mitochondrial genome of O. infernale predicted by tRNAScan-SE 2.0.

The mtDNA control region of O. infernale is 1074 bp long (15731–16804), encoded in the H-strand, and flanked by tRNA^Pro and tRNA^Phe at the 5' and 3' ends, respectively (Fig. 1; Table 1), which is consistent with our understanding of mitogenome structure and organization in fishes (Lee et al. 1995; Rasmussen and Arnason 1999; Satoh et al. 2016). Ophisternon infernale CR nucleotide composition is T=33.1%, A=36.7%, G=10.9%, and C=19.3%, with A+T content (69.8%) larger than that of the entire mitogenome but similar to that of other fishes including synbranchids (Li et al. 2016; Han et al. 2018). Like in other fishes, CR in O. infernale is divided into three domains: a central conserved domain flanked and two hypervariable domains (upstream and downstream). Three conserved sequence blocks (CSBs) were detected at the central conserved domain (CSB-F, CSB-E, CSB-D) as well as at the downstream hypervariable region (CSB1, CSB2, CSB3) (Fig. 4). Although additional CBSs have been identified for the central conserved domain (CSB-B, CSB-C) in mammals (Southern et al. 1988), the three identified herein for O. infernale are those commonly found in fishes (Broughton and Dowling 1994; Chen et al. 2012). The upstream hypervariable domain in the CR of O. infernale has a length of 256 bp and includes two copies of the motif TACAT and three copies of palindromic motif ATGTA. A change in the motif sequence (TGCAT) was observed in C. semiplotum and S. chuatsi but not in O. infernale. Compared to those from the central conserved domain, CSBs in the downstream hypervariable domain displayed larger variation across the three fish species compared. Notably, CSB2 and CSB3 were slightly more conserved than CSB1, a pattern that has been reported for other fishes (Chen et al. 2012).

Figure 4. 

Comparison (multiple sequence alignment) of the mtDNA control region of O. infernale with those of fellow teleosts Siniperca chuatsi and Cyprinion semiplotum. The alignment displays the three canonical domains distinguished by Termination Associated Sequences (TAS) of the upstream hypervariable region (in red), central conserved domain blocks (CSB-F, CSB-E, CSB-D) (in blue), and conserved sequence blocks of the downstream hypervariable region (CSB-1, CSB-2 and CSB-3) (in green).

Results from K_A/K_S analyses (Fig. 5) indicate that most mtDNA PCGs in synbranchiform fishes have evolved under strong purifying selection (K_A/K_S <<1), preventing major structural and functional protein changes. Exceptions to this general pattern were observed for COX1 and NAD6 in synbranchids (Fig. 5a) and for NAD4 and NAD6 in mastacembelids (Fig. 5b), where significant signals of positive selection were detected. Studies in different groups of animals, including cephalopods (Almeida et al. 2015), rodents (Tomasco and Lessa 2011), and humans (DeHaan et al. 2004), have linked amino acid replacements in NAD6 to adaptive selection to hypoxic conditions. Because numerous synbranchiform species are known to be fossorial and to inhabit low-oxygen waters, the observed signature of positive selection in NAD6 might be related to adaptation to decreased oxygen availability. Notably, a recent comparative mitogenomic study of the African tilapias Oreocrhomis andersonii and O. macrochir similarly uncovered a pattern of positive selection in NAD6 suggestive of adaptation in response to changing environments (Bbole et al. 2018). In contrast to the pattern observed for NAD6, selection in COX1 and NAD4 is completely conflicting between synbranchiform families. While in mastacembelids COX1–like most mitochondrial genes–has evolved under purifying selection (K_A/K_S <1), the opposite happens in synbranchids. Although speculative at this point, the fact that half of our synbranchid dataset consists of troglomorphic cave-dwelling species (O. infernale and O. candidum) (vs. none in the mastacembelid dataset) could explain the observed differences in COX1 selection patterns. Compared to surface waters, subterranean waters such as those of karst environments that harbor populations of O. infernale and O. candidum (in Mexico and Australia, respectively) contain low dissolved oxygen (Huppop 2000). Because of its role in aerobic metabolism, COX1 might therefore be a target of directional selection promoting the evolution of more metabolically efficient variants in hypogean lineages (Boggs and Gross 2021). In contrast, the observed conflicting patterns of selection in NAD4–another gene involved in cellular respiration–between mastacembelids (positive) and synbranchids (purifying), do not seem to be readily explained by ecological differences related to cave life.

Figure 5. 

Patterns of selection in mtDNA PCGs of synbranchiform fishes. Results from K_A/K_S ratio analysis on mitochondrial PCGs (x-axis) in synbranchiform fishes of the families Synbranchidae (a) and Mastacembelidae (b).

Our understanding of phylogenetic relationships in synbranchiform fishes is incipient compared to that of other teleost groups. Despite the fact that for the past two decades molecular systematics has been routinely employed to refine and update the classification of fishes and our knowledge of their evolutionary history (Betancur-R et al. 2017), a comprehensive molecular phylogeny of the Synbranchiformes has yet to be proposed. Apart from a phylogenetic study focused on Central American synbranchids (Perdices et al. 2005), no studies have investigated synbranchiform relationships using comparative DNA sequence data. Surprisingly, recent phylogenetic studies focused on higher-level relationships among major lineages of bony fishes (Betancur-R. et al. 2013; Betancur-R et al. 2017) resulted in the reassignment of armored sticklebacks (family Indostomidae, traditionally placed in the suborder Gasterosteoidei, order Scorpaeniformes) to the order Synbranchiformes. Although the classification of indostomids as gasterosteoids had been previously questioned on the basis of mitogenomic evidence (Miya et al. 2003, 2005; Kawahara et al. 2008), it was not until the phylogenetic classification of Betancur et al. (2013, 2017) that the family Indostomidae was transferred to the order Synbranchiformes. This proposal, however, was not adopted by the most authoritative contemporary standard references of fish systematics (Van Der Laan et al. 2014; Nelson et al. 2016), on the grounds of lack of morphological support and the need for further corroboration. It should be noted that previous molecular phylogenetic studies that cast doubt on the traditional placement of indostomids, whether based on “legacy” markers (Betancur-R. et al. 2013; Betancur-R et al. 2017) or complete mitochondrial genomes (Miya et al. 2003, 2005; Kawahara et al. 2008), relied on a very limited representation of synbranchiform diversity. In contrast, our phylogenetic analysis used mitogenomic data from a comparatively larger taxon sampling that included eight synbranchiform species from five genera (Ophisternon, Synbranchus, Monopterus, Mastacembelus, and Macrognathus) and two families (Synbranchidae, Mastacembelidae). Notably, our phylogenetic results (Fig. 6) corroborate the notion that indostomids are more closely related to synbranchiforms than to gasterosteoids. Nevertheless, contrary to the findings of studies that have recently challenged the traditional classification of indostomids with respect to synbranchiforms (Kawahara et al. 2008; Betancur-R. et al. 2013; Betancur-R et al. 2017), our inferred phylogenetic placement of Indostomus does not render Synbranchiformes paraphyletic. With the caveat that our sampling of synbranchiforms and closely related lineages is only partial, our results imply that indostomids are in fact the sister lineage of the order Synbranchiformes. While this phylogenetic pattern (topology) might be considered sufficient for lumping indostomids with synbranchiforms, examination of relative branch lengths (Fig. 6) suggests that Indostomus is indeed a highly divergent lineage. In order to acknowledge their genetic and morphological (Britz and Johnson 2002) distinctiveness, indostomids may in fact warrant an order of their own. Within Synbranchiformes, our results remarkably do not support the monophyly of the synbranchid genus Ophisternon, for O. infernale is resolved as more closely related to Synbranchus marmoratus than to O. candidum (Fig. 6). While at first sight this novel finding of a sister-group relationship between O. infernale and S. marmoratus is certainly unexpected, this hypothesis might not be that far-fetched from a biogeographic perspective, and when considering both the striking external morphological similarity between the two genera and the taxonomic ambiguities surrounding the classificatory history of the group (Rosen and Greenwood 1976). Synbranchus is restricted to the New World and comprises three species: S. marmoratus (Central and South America), S. madeirae (Madeira River basin, Bolivia), and S. lampreia (Pará, Brazil). Ophisternon as currently delimited exhibits an essentially Gondwanan distribution, with six valid species distributed in Middle America (O. infernale, O. aenigmaticum), Australia (O. candidum, O. gutturale), South Asia and Western Pacific (O. bengalense), and West Africa (O. afrum). Assuming that Gondwanan drift vicariance is the main process responsible for the present-day globally disjunct distribution of the genus (Rosen 1975), the split between the Mexican-endemic O. infernale and the West Australian-endemic O. candidum should be at least as old as the Middle Jurassic separation of Eastern Gondwana (Antarctica, Madagascar, India, and Australia) from Western Gondwana (South America and Africa), dated at ca 165 Ma (McLoughlin 2001). From this it follows that the split between Ophisternon and Synbranchus should be even older. Notably, the only phylogenetic study that has investigated divergence times via molecular dating in a group of synbranchiforms (Perdices et al. 2005) estimated a comparatively much younger age (< 20 Ma) for the split between Ophisternon (aenigmaticum) and Synbranchus (marmoratus). Although marine dispersal and extinction could be invoked in an attempt to reconcile biogeographic patterns with our admittedly limited knowledge of the timescale of synbranchiform diversification, the paraphyly of Ophisternon remains problematic. Our phylogenetic results coupled with the abovementioned estimates of synbranchid divergence times (Perdices et al. 2005) lead us to hypothesize that perhaps New World species of Ophisternon (O. infernale and O. aenigmaticum) are in fact more closely related to Synbranchus species than to the remaining Ophisternon species. As such, New World species of Ophisternon would have to be transferred to the genus Synbranchus. This phylogenetic scenario is also compatible with a likely very recent origin of the cave-dwelling O. infernale. Although there is virtually no information regarding the timing of origin and colonization of the fishes that inhabit the cenotes and submerged caves of the YP karstic aquifer (Arroyave et al. 2021), these aquatic habitats are supposed to be extremely young, effectively established not before 20,000 years ago, at the end of the last glacial maximum in the Northern Hemisphere, when rising sea levels eventually resulted in the flooding of karstic sinkholes and dry caves (Coke IV 2019). Such a recent origin for O. infernale is certainly much easier to explain as a result of speciation from a fellow New World lineage, such as O. aenigmaticum or S. marmoratus. Regardless of the appeal and feasibility of these hypotheses concerning the systematics of New World Ophisternon in general and the origins of O. infernale in particular, our phylogenetic findings and their interpretation need to be taken with caution because of their absolute reliance on mtDNA only. It is well known that the mitochondrial genome is effectively a single locus (Avise 2012), that individual gene and species trees are not always congruent (Maddison 1997), and that nuclear and mtDNA inheritance patterns are not always congruent either (Funk and Omland 2003). Notwithstanding these limitations, our results emphasize the pressing need for a comprehensive systematic and biogeographic study of synbranchiform fishes, ideally based on genome-wide sequence data.

Figure 6. 

Phylogenetic relationships of major synbranchiform lineages. Molecular phylogeny based on comparative mitochondrial PCGs from relevant available mitogenomes and the newly generated herein for O. infernale. Troglobitic cave-dwelling species are marked with an asterisk to distinguish them from surface-dwelling ones. Outgroup taxa not shown. Colored circles on nodes indicate degree of clade support as determined by bootstrap values.