Skip to main content

Advertisement

First complete mitochondrial genome of the South American annual fish Austrolebias charrua (Cyprinodontiformes: Rivulidae): peculiar features among cyprinodontiforms mitogenomes

Abstract

Background

Among teleosts, the South American genus Austrolebias (Cyprinodontiformes: Rivulidae) includes 42 taxa of annual fishes divided into five different species groups. It is a monophyletic genus, but morphological and molecular data do not resolve the relationship among intrageneric clades and high rates of substitution have been previously described in some mitochondrial genes. In this work, the complete mitogenome of a species of the genus was determined for the first time. We determined its structure, gene order and evolutionary peculiar features, which will allow us to evaluate the performance of mitochondrial genes in the phylogenetic resolution at different taxonomic levels.

Results

Regarding gene content and order, the circular mitogenome of A. charrua (17,271 pb) presents the typical pattern of vertebrate mitogenomes. It contains the full complement of 13 proteins-coding genes, 22 tRNA, 2 rRNA and one non-coding control region. Notably, the tRNA-Cys was only 57 bp in length and lacks the D-loop arm. In three full sibling individuals, heteroplasmatic condition was detected due to a total of 12 variable sites in seven protein-coding genes. Among cyprinodontiforms, the mitogenome of A. charrua exhibits the lowest G+C content (37 %) and GCskew, as well as the highest strand asymmetry with a net difference of T over A at 1st and 3rd codon positions. Considering the 12 coding-genes of the H strand, correspondence analyses of nucleotide composition and codon usage show that A and T at 1st and 3rd codon positions have the highest weight in the first axis, and segregate annual species from the other cyprinodontiforms analyzed. Given the annual life-style, their mitogenomes could be under different selective pressures. All 13 protein-coding genes are under strong purifying selection and we did not find any significant evidence of nucleotide sites showing episodic selection (dN >dS) at annual lineages. When fast evolving third codon positions were removed from alignments, the “supergene” tree recovers our reference species phylogeny as well as the Cytb, ND4L and ND6 genes. Therefore, third codon positions seem to be saturated in the aforementioned coding regions at intergeneric Cyprinodontiformes comparisons.

Conclusions

The complete mitogenome obtained in present work, offers relevant data for further comparative studies on molecular phylogeny and systematics of this taxonomic controversial endemic genus of annual fishes.

Background

The Neotropical genus Austrolebias (Cyprinodontiformes: Rivulidae) is a locally endemic group of annual fishes. It is distributed from northern and northeastern Argentina, Paraguay, Uruguay and southern Brazil throughout the Paraná-Plata basin and Patos-Merín system [1]. It includes 42 known species (http://www.fishbase.org/search.php) highly variable in morphology and behavior [2, 3]. They live in temporary ponds formed during rainy seasons where each generation completes full life cycle within 1 year. The population survives dry season in the form of eggs buried in the mud. During the subsequent rainy season, the ponds refill, the eggs hatch and the larvae rapidly grow to sexual maturity and reproduce [4, 5]. According to morphological and molecular characters these taxa are divided into five species groups: A. adloffi, A. alexandri, A. bellottii, A. elongatus and A. robustus [1, 6, 7]. Molecular phylogenetic analyses based on some mitochondrial genes, supported the monophyly of the genus, but exhibited limited resolution at deeper nodes of the trees. This lack of information could be due to a saturation of nucleotide changes in the genes analyzed and/or to the mode of evolutionary diversification of the Austrolebias taxa involved [6, 8]. High levels of intrageneric sequence divergence and saturation have been detected in the cytochrome b gene, and a bias in favor of AT changes has been detected at silent sites in this gene [6]. As proposed by Martin and Palumbi [9], high rates of nucleotide substitutions observed in mitochondrial genes could be related to some physiological variables such as generation time, life span, age at first reproduction, rate of population increase, and metabolic rate. In fact, high rates of growth were found related to an increase of temperature up to 25 °C in natural and laboratory populations of A. viarius [10].

In the present work we describe for the first time the complete mitochondrial genome of A. charrua (Costa and Cheffe 2001), a member of the A. adloffi species group, highlighting its structural, compositional and evolutionary features. While a wealth of data on mtDNA sequence and gene organization of fish mitogenomes have been collected in the last years, this is the first reference mitogenome for the genus, and just the second within the Rivulidae family [11]. To date, 1750 mtDNA sequences (update June 2, 2015) are deposited in the Mitochondrial Genome Database of Fish (MitoFish: http://mitofish.aori.u-tokyo.ac.jp/), of which 15 correspond to Cyprinodontiformes. We took advantage of some of these available mitogenomes and compared them in an evolutionary framework, searching for putative molecular patterns associated to the annual life style. For the different mitochondrial coding regions, we also evaluated their performance as molecular markers for phylogenetic reconstruction at different taxonomic levels.

Methods

Sample collection and mtDNA extraction

A total of three individuals of A. charrua (AUS4, AUS5 and AUS7) were analyzed in the present work. They belonged to the same F1 generation obtained at the lab while the parents were collected in a temporal pond located in “La Coronilla” (Rocha Department, Uruguay). AUS4 and AUS5 were young individuals of 3 months old, while AUS7 was an adult of 5 months old. All tissues and voucher specimens were deposited in the Sección Genética Evolutiva, Facultad de Ciencias, Universidad de la República, Montevideo, Uruguay.

The mtDNA was isolated from muscle tissue (fixed in 95 % ethanol) of freshly sacrificed animals by an over-exposure to a solution of 1‰ 2-phenoxyethanol (Sigma). The sampling protocol was approved by the CNEA (Comisión Nacional de Experimentación Animal) from Uruguay.

Sequencing

Sequencing libraries were generated using Nextera XT kit (Illumina, USA) from 1 ng of input DNA. Tagmentation of genomic DNA, PCR amplification (12 cycles) with addition of index primers and library normalization, were performed according to manufacturer’s protocol. Quality and length of libraries were assessed with Agilent High Sensitivity DNA Kit (Agilent, USA) using the 2100 Bioanalyzer (Agilent, USA) before the final normalization step. Multiplexed libraries were sequenced on an Illumina Genome Analyzer IIx platform run for 66 cycles in single-end mode, at the Institut Pasteur of Montevideo.

Amplification and sequencing of the control region

Amplification of the control region was performed by PCR. Primers were designed at both ends of the automatically assembled CR (tRNA-Pro and tRNA-Phe). Two sets of specific primers: ACR-1F (5′- ACC TGT TCC TCT AGC ACC CA -3′), ACR-1R (5′- TGT AGG AGG CAT TTA AGG TGC A -3′) and ACR-2F (5′- TTT CTG GCC CAC AAG AGA CC -3′), ACR-2R (5′- TGC TCA TGA AAC TTT TTA GGG TTT -3′) were created using Primer3 v. 0.4.0 software. The PCR was carried out in a 10 μL total volume of 10× buffer, MgCl2 2.5 mM, dNTPs 0.2 mM, 0.4 mM of each primer, 0.5 U of Taq DNA polymerase (Invitrogen) and 10 ng of DNA, under the following conditions: one denaturation step at 94 °C for 5 min, followed by 35 cycles of 94 °C for 45 s, 62 °C for 45 s and 72 °C for 1 min, and a final elongation step at 72 °C for 7 min.

Sequencing reactions were performed on each template on the automated ABI PRISM 377 DNA Sequencer (MACROGEN Inc., Korea).

Mitogenome assembly and annotation

Raw reads were quality filtered and trimmed to a final length of at least 50 bp. As indicated previously [12], small mitochondrial genomes are present in variable copy numbers in total DNA mixed samples, therefore transcriptome-designed tools are better suited for mitogenome assembly. Following [12], contigs were assembled using SOAPdenovo-Trans [13] with k-mer length of 31 and remaining default options. As both nuclear and mitochondrial contigs were assembled with this strategy, a local blastn [14] was ran against a database of cyprinodontiforms mitogenomes downloaded from MitoFish [15] (MitoFish version 2.89, update March 2014, 1361 complete mitochondrial sequences), to recover those contigs of mitochondrial origin.

Selected contigs were annotated with MitoAnnotator [15], a highly accurate and automated pipeline specific for fish mitogenomes, followed by manual curation in Artemis [16]. Transfer RNA genes were annotated based on their potential cloverleaf secondary structure using MiTFi [17], a tool for accurately locating tRNA genes within mitogenomic sequences (MiTFi was run through MitoAnnotator [15]).

Variant analysis

Reads from each sample were mapped to the chosen reference mitogenome in order to search for the presence of nucleotide substitutions and indels within mitogenomes (i.e. to search for heteroplasmy, defined as the presence of more than one type of mitogenome within one individual) and among siblings. This was performed with a combination of Burrows-Wheeler Aligner (bwa; default options) [18] and SAMtools (samtools view with -q 30; samtools mpileup with options -B -C 0 -Q 30) [19]. Of note, both the read mapping quality (−q 30) and base quality (−Q 30) were required to have a minimum value of 30 in a Phred scoring scheme to be considered in the variant analysis.

Compositional analyses

Nucleotide composition, as well as codon and amino acid usage were calculated in the R environment [20] using seqinR [21] and ca packages [22]. Besides whole mitogenome analyses, calculations were performed on multiple alignments of the 12 protein-coding genes present in the H strand (ND6 was excluded due to its location in the L strand where there are different compositional biases). In these comparative analyses, the following Cyprinodontiformes were included: Kryptolebias marmoratus (Aplocheiloidei: Rivulidae), Nothobranchius furzeri (Aplocheiloidei: Nothobranchiidae), Aplocheilus panchax (Aplocheiloidei: Aplocheilidae), Cyprinodon rubrofluviatilis (Cyprinodontoidei: Cyprinodontidae) and Fundulus olivaceus (Cyprinodontoidei: Fundulidae). To note, while K. marmoratus is the closest relative to A. charrua, N. furzeri is the only annual killifish other than A. charrua in this species sample. Cyprinodon rubrofluviatilis and F. olivaceus were considered outgroup lineages. The phylogenetic tree for these taxa according to Pohl et al. [23] topology is shown in Additional file 1. Multiple nucleotide alignments were built guided by previous alignments of amino acid sequences using Muscle [24] and Pal2Nal [25]. For each alignment, ancestral sequence reconstruction was performed with a joint-likelihood reconstruction method in the codon-state space [26, 27].

Protein-coding evolutionary analyses

Analyses of signatures of positive selection on codon sequences were carried out using the suite of routines implemented in http://www.datamonkey.org [28]. Uploaded alignments were matched with the reference species tree (Additional file 1) and F. olivaceus fixed as outgroup lineage if required. Global dN/dS values were estimated with SLAC [27]. MEME [29] and BS-REL [30] models were used to search sites and branches under episodic and pervasive selection. For MEME, sites were posited to be under directional selection at selected branches if they: i) showed non-synonymous to sysnonymous substitution rate β+ >1; ii) were significant at α = 0.05 level; and iii) also showed an Empirical Bayes Factor >20. Both MEME and BS-REL avoid the a priori and rigid partitioning of branches into “foreground” and “background” classes. BS-REL uses a random effects likelihood framework, in which ω can take one of three values along branches (ωb- ≤ωbN ≤1 ≤ωb+) to explore every branch-site combination. Sequential likelihood ratio testing is used to identify branches with significant amounts of diversifying selection.

Phylogenetic reconstruction

Nucleotide alignments for individual protein-coding genes were used to infer maximum likelihood phylogenies with PhyML v3.1 [31]. Also, a “supergene” phylogeny was based on a concatenated alignment of all 13 protein-coding genes. In each case, five random starting trees were used to search the tree space by NNI [32] and SPR [33] methods. Survey of branch support was conducted with a Bayesian-like transformation of approximate likelihood ratio test (aBayes; [34]). Nucleotide substitution model parameters were determined according to the best suited sequence substitution model selected by ModelGenerator v0.85 [35], following Bayesian information criterion. For most genes, the HKY+I+G model [36] was the chosen one, while the GTR+I+G [37] suited best for COI, COIII and the concatenated alignment. Default values were kept for remaining PhyML options. Inferred gene trees were visualized to study their concordance to the working species tree hypothesis (Additional file 1). All inferences were recalculated after third codon positions were removed from the alignments (Additional file 2).

Results and discussion

Austrolebias charrua mitogenome

The mtDNA sequence of the annual fish A. charrua was determined by massively parallel sequencing methods. After the quality filtered and trimmed reads were assembled with SOAPdenovo-trans [13], almost the whole mitochondrial genome was reconstructed in only one contig in the samples AUS4 and AUS7, while two contigs were assembled in case of AUS5 (all contigs reached mean coverage above 30×; Table 1). The mitochondrial control region is a non-coding element vital for the initiation of both mtDNA replication and transcription. Often, it contains repetitive elements and could appear as a duplicated region as shown in the mitogenome of the closely related non-annual rivulid species K. marmoratus [11]. These features pose some difficulties during the automatic assembly of mitogenomes where duplicated elements are prone to be collapsed in the final assembly. Therefore, Sanger sequencing was used to confirm the size and sequence of the control region. Eight CR sequences obtained from the samples AUS1 and AUS7, revealed that it spans from the tRNA-Pro to the tRNA-Phe and is 1349 bp long, showing an internal duplication of 180 bp which had not been resolved during the automatic assembly (Additional file 3).

Table 1 Illumina sequencing and SOAPdenovo-trans assembly

The de novo assembled mitochondrial contig of AUS7 was manually curated and edited using this additional sequence information. The resulting complete mitochondrial sequence of 17,271 bp (the initial contig was 17,103 bp) was deposited in the GeneBank (accession number KP718940) and considered our reference mitogenome for A. charrua (Fig. 1). Regarding gene content and order, it exhibits the typical pattern of vertebrate mitogenomes [38]. It contains the full complement of 13 protein-coding genes, 22 tRNA genes, two rRNA genes (12S and 16S) and one non-coding control region. Eight tRNAs (tRNA-Gln, tRNA-Ala, tRNA-Asn, tRNA-Cys, tRNA-Tyr, tRNA-Ser(1), tRNA-Glu and tRNA-Pro) and the protein-coding gene ND6 were encoded by the light (L) strand, whilst the remaining genes were encoded by the heavy (H) strand (Table 2). ATG codons initiate 12 of the protein-coding genes, COI being the exception (GTG start codon). Complete TAA stop codons are present in COI, COIII, ATPase8, ND1, ND4L, ND5 and Cytb, while a complete TAG codon is present in ND6. Remaining protein genes end at incomplete codons (i.e. T or TA), which are converted by polyadenylation into TAA after transcription and processing [39]. Of 22 transfer RNA genes, 21 were within the size range of 66–74 bp (Table 2), and each of them folded into a typical cloverleaf secondary structure. On the other hand, tRNA-Cys was only 57 bp in length, shorter than any known annotated tRNA in Cyprinodontiformes (z-score = −5.74) and in the predicted secondary structure lacks the D-loop arm (Fig. 2a). While D-armless tRNA-Cys reports do not exist for fishes [17], several independent occurrences do exist in metazoans, including mammals. Atypical tRNAs (D or T-armless) would be functional [17, 40, 41]. To note, in K. marmoratus, tRNA-Cys is 65 bp in length and has the typical secondary structure (Fig. 2b). Further analysis of new Austrolebias taxa mitogenomes could clarify whether this D-armless tRNA-Cys is a genus or species-level feature.

Fig. 1
figure1

Map of the mitochondrial genome of Austrolebias charrua (GenBank accession number KP718940). Visual representation created by Circos [51] where the innermost circle represents G+C content per every 5 bp (darker lines are the higher G+C content) and the external one represents the two mtDNA strands. Protein-coding genes are in black; tRNAs in red and are designated by their three letter abbreviations; light brown are the rRNAs and in dark brown is the CR (D-loop)

Table 2 Organization of the mitochondrial genome of Austrolebias charrua
Fig. 2
figure2

Putative secondary structure of the tRNA-Cys encoded in the mtDNA of (a) A. charrua and (b) K. marmoratus. The tRNAs were annotated using MiTFi tool run through MitoAnnotator [15]

Mitogenome differences among full siblings

The availability of deep sequencing data allows the discovery of hidden variation and heteroplasmy. Reads were mapped against the assigned A. charrua reference genome (AUS7) and mpileup files generated as mentioned in methods. R scripts [20] and careful inspection with IGV [42] were used to determine the presence of nucleotide variants. For those mitogenome positions where the mapping coverage was at least 20×, a variant was defined if: i) there was only one alternative nucleotide; ii) the alternative base showed up in at least two reads, with at least one read mapping to each DNA strand (this requirement is aimed to discard sequencing errors linked to specific biases in DNA composition as well as spurious support due to existence of duplicated reads); and iii) alternative nucleotides were far away from the read ends and the neighborhood positions were error-free. These requirements were defined after joint consideration of reviewed references [4347] and the features (e.g. mean sequencing coverage) of our data. A total of 12 variants were detected in protein-coding genes of the mitogenome of the full siblings analyzed, five in the mtDNA of AUS4 and seven in AUS5. The diversity of changes and their putative effects on the coding sequence are shown in Table 3. The transversion C to A identified in the COI gene (position 6245) was present in both samples, while the transition T to C detected in the Cytb gene (position 15,735) was present in all full siblings, including the reference genome (Table 3). The coverage depth of this experiment would not allow the detection of de novo mutations, so the variants observed would be all putative changes already carried in the maternal mitogenomes which had changed their frequency in the different siblings as consequence of sampling effects [44, 46]. To note, above 28 % of the reads mapping to position 15,735 carried the alternative nucleotide cytosine. This heteroplasmic position corresponds to the codon 368 of the Cytb locus, and despite of being a transition it produces a non-synonymous change from Phe to Leu in the translated sequence (TTT to CTT codon change). Phe and Leu are both non-polar amino acids but the change is from an aromatic residue to a non-aromatic one. However, hydrophobic aromatic amino acids can be sometimes substituted by aliphatic residues of a similar size, as in the present case. In addition, two related findings preclude considering this variant as having a notorious phenotype in the organism. First, given a multiple alignment of metazoan Cytb sequences, the position is highly variable, being Leu the consensus amino acid (Additional file 4). Second, the involved residue is part of the H helix in the intermembrane domain, close to the C-terminus of the protein. This is one of the least conserved eight helices in the protein, lacks heme binding-sites and does not interact with other components of the cytochrome bc1 complex [48, 49].

Table 3 Nucleotide variants in the mitogenomes of three full siblings of A. charrua

Comparative mitogenomics of Cyprinodontiformes: structural and evolutionary features

The whole mitogenomes compared in present work exhibit a similar total length, but few differences have been detected regarding the CR (Table 4). Kryptolebias marmoratus displays an additional control region located between tRNA-Leu and ND1 [11] while the annual species N. furzeri shows a more complex rearrangement where a duplication of tRNA-Gln sites in place of tRNA-Glu, which in turn is located in the L strand between the control region and the tRNA-Phe. Besides, a long intergenic region is interspersed between tRNA-Thr and tRNA-Pro [50]. A Circos plot [51] highlighting synteny and identity level between A. charrua and each of the five mitogenomes used for comparison is shown in Fig. 3. The most notorious feature in the plot is the lack of similarity between A. charrua and the last region of N. furzeri mitochondrial chromosome (about 3500 pb), despite homologous regions can be identified in more distant species as C. rubrofluviatilis and F. olivaceous (Figs. 3a, c). Considering that the mitogenome of N. furzeri was built as part of a whole genome shotgun approach resulting in a narrow 5.4 Mb draft assembly (out of an estimated 1.59 – 1.92 Gb total genome size) [52] a putative miss-assembly should be discarded before further comparative analyses.

Table 4 Nucleotide composition of all Cyprinodontiformes compared in present work
Fig. 3
figure3

Circos plot mapping the synteny and identity level of A. charrua mitogenome against each of the other five mitogenomes studied (a), K. marmoratus (b) and N. furzeri (c). Idiograms and red ribbons represent the similarity after pairwise blastn searches. In A. charrua idiogram, the 13 coding genes are coloured in blue. In all idiograms, control region are coloured in black. In (c), N. furzeri unique features (tRNA-Gln, tRNA-Glu and the interpersed region between tRNA-Thr and tRNA-Pro) are shown in orange. The figure was produced using Circos software [51], a tool for graphical representation of genome data, through the Circoletto implementation [88]. AC (A. charrua), KM (K. marmoratus), NF (N. furzeri), AP (A. panchax), CR (C. rubrofluviatilis) and FO (F. olivaceus)

Table 4 shows GenBank accession numbers, total length, CR length and the compositional properties of mitogenomes from A. charrua, N. furzeri, K. marmoratus, A. panchax, C. rubrofluviatilis and F. olivaceus. Regarding the A. charrua mitogenome G+C content is 37 %, while both GCskew (GCskew: (G-C)/(G+C)) and ATskew (ATskew: (A-T)/(A+T)) are negative, pointing to strand asymmetry in the nucleotide composition (excess of cytosine and thymine in the H strand, see below). It is worth noting that annual species present the lowest G+C content and GCskew, in spite of K. marmoratus, a sister taxon of A. charrua. In the case of ATskew, values are positive and negative among considered species, having A. charrua the highest asymmetry (−0.102). In fact, this relatively strong negative value is highly unusual in Actinopterygii and even Chordata [53]. Given the differences in length of the annotated control regions (Table 4) and the inherent difficulty to sequence this region, we hypothesized that their putative incompleteness in some of the species analyzed could have biased the above results. For this reason, we repeated the compositional analyses taking into account the sequences of the 12 protein-coding genes located in the H strand. The G+C content is again lower in both annual killifish, being the difference (as expected) larger in the first and third codon positions (Table 4 and Additional file 5). In order to determine if the lower G+C content in annual killifish is an ancestral or recently derivative character, we performed a sequence reconstruction of ancestral nodes with a joint-likelihood method in the codon-state space [26, 27]. Estimation of G+C content at the internal nodes of our reference topology (Additional file 1) confirms that the lower G+C content is a recently and independently derivative character state in both annual lineages (Additional file 6).

Correspondence analyses of nucleotide composition of the coding-genes of the H strand (Fig. 4) and on global codon usage (Additional file 7) segregates annual killifishes from the remaining Cyprinodontiformes species. Codons rich in A and T at third position show the highest weight in this pattern, reflected in the first axis. The second axis segregates both aplocheiloid non-rivulids fishes (N. furzeri and A. panchax) from the remaining cyprinodontiforms and separates C from G and A from T. While organisms position in the first axis correlates with mitogenome G+C content (Spearman rank correlations were: GC: Spearman’s ρ = −1.00, p-value <0.05; GC3: Spearman’s ρ = −0.99, p-value <0.05 and GCskew: Spearman’s ρ = 0.90, p-value <0.05), second axis correlates to ATskew (ATskew: Spearman’s ρ = 0.93, p-vale <0.05). Both analyses revealed that annual species appear associated with an increase in A or T content in 1st and 3rd positions. Moreover, the second axis separates A. charrua from N. furzeri, the first species associated with A1 and A3 while the second appear richer in T1 and T3. When studying per gene codon usage considering only A. charrua, K. marmoratus and N. furzeri, there is a clear difference in per species codon usage, being both annual taxa more similar, in spite of the closer phylogenetic relationship between A. charrua and K. marmoratus (Additional file 8).

Fig. 4
figure4

Correspondence analysis of nucleotide composition of all protein coding genes located in the H strand of all Cyprinodontiformes analyzed. Blue dots represent the species: AC (A. charrua), KM (K. marmoratus), NF (N. furzeri), AP (A. panchax), CR (C. rubrofluviatilis) and FO (F. olivaceus). Red triangles indicate the nucleotides at 1st, 2nd and 3rd codon position

Furthermore, comparison of base usage among organisms shows an important difference of C over G at the three codon positions (Table 4 and Additional file 5), a fairly common characteristic [54]. However, only A. charrua displayed a net difference of T over A at 1st and 3rd positions. This feature could be explained by the fast growth and high metabolic rates present in annual fishes [10, 55]. These two processes involve multiple rounds of organelle replication per cell division with a concomitant increase of replication errors and of the concentrations of ROS in the mitochondria environment [5659]. Reactive oxygen species promote GC to AT mutations through the deamination of cytosine and the oxidative conversion of guanine to 8-oxo-guanine. Hydrolytic deamination of cytosine or adenine and oxidation of guanine are among the described causes of mutational damage [60]. Given these sources of mutations and the pattern observed in A. charrua, an increase in the rate of hydrolytic deamination of cytosine in the H strand could be the main mechanism associated to the increase in thymine. In turn, the asymmetric pattern in base composition could be explained by the long-standing strand-displacement model of mtDNA replication: in this asymmetric process a portion of the H strand would remain a period of time in single strand state, which could lead to different mutational bias in H and L strands [61]. While this model of mtDNA replication has been challenged in recent years (reviewed in [62]), asymmetric trends in base composition have been confirmed [54].

Other various hypotheses have tried to explain the high adenine and thymine (AT) content of almost all mitochondrial DNAs [63]. One is that natural selection contributed to the high AT content of mitochondrial genomes, where selection for translational efficiency and accuracy shaped the nucleotide composition of codons in organelle genes, in some cases enriching the thymine content of synonymous sites [64]. Others postulated that AT richness is an adaptation for metabolic efficiency, due the increased energetic costs of producing C vs. T and G vs. A and the varying abundance of A/T vs. G/C nucleotides during organelle DNA synthesis [59, 65, 66]. As expected, codon usage pattern also translates into differences in amino acid usage (Additional file 9). While relative differences are small, there is a higher usage of phenylalanine in both rivulids, while both annual taxa share an increment of tyrosine, lysine and methionine usage.

Selective constraints in coding regions

Variations in the subunits of the oxidative phosphorylation pathway have been linked to different life-history traits and environmental adaptations [6770]. In a landscape of widespread strong purifying selection due to functional constraints, signals of positive selection have been reported for individual codon sites (or amino acids) of taxa known to have high energetic demanding lifestyles such as flying in bats [71] and unusual oxygen availability such as high altitude Caprini and subterranean octodontoid rodents [72, 73], among other conditions [70, 74, 75].

Given the annual life-style of A. charrua and N. furzeri we hypothesized that, in comparison to other cyprinodontiforms, they have been under a different selective regime and this could have left molecular signatures at the mitochondrial genome level. In fact, as explained above, both taxa present an independent decrease in their mitogenome G+C content that could be due to their fast growth and high metabolic rates. For many years codon models have been used to estimate the ratio of non-synonymous to synonymous substitution rate (dN/dS, often denoted as ω) and positive selection can be inferred whenever the estimated ratio dN/dS significantly exceeds one [76, 77]. The early codon models have been extended to permit the intensity of selection to vary among sites within a gene, among branches within a tree or both (branch-site models), greatly improving the power to detect positive selection. More recently a mixed effects model of evolution (MEME), a highly sensitive branch-site random effects phylogenetic method, was presented [29]. As it is capable of detecting episodic adaptation, as would be the case of the two annual lineages under study, we used MEME to search for evidence of positive directional selection in our mitochondrial codon sequences dataset.

As expected for their vital role in the oxidative phosphorylation, all 13 protein-coding genes are under strong purifying selection, with dN/dS values below 0.3 (Additional file 10). The ATPase6, ATPase8 and ND6 genes showed the highest dN/dS values, while components of the cytochrome oxidase complex (COI, COII and COIII) showed the lowest rates. These results are in perfect agreement with a previous comparative analysis in metazoans [53]. The analysis conducted with MEME [29] inferred four sites where some branches experienced a significant greater non-synonymous than synonymous substitution rate. These sites are: codon 121 of ATPase6; codon 47 of ATPase8; codon 369 of Cytb and codon 220 of ND2 (positions follow A. charrua annotation; Table 5). Molecular footprints of positive selection have been previously associated to mitochondrial ND2 in at least one work [70] and signatures of adaptation have been repeatedly detected in Cytb, ATPase6 and ATPase8 [70, 71, 73, 74, 78, 79]. Noteworthy, in three of these sites the episodic directional selection would have happened in the ancestral branch leading to rivulids (ATPase6, ATPase8 and ND2), while in the Cytb site it would have occurred in the tree branch leading to N. furzeri and A. panchax (Table 5). We did not find any significant evidence of sites showing episodic selection at both annual lineages. However, as MEME authors denote [29], it is difficult to accurately identify individual positively selected branches at an individual site, thus conclusions about specific lineages under selection in particular sites must be taken with caution.

Table 5 Codon sites and branches under episodic and pervasive selection inferred by mixed effects model of evolution (MEME [29])

As explained, one cannot simultaneously infer the site and the branch subject to diversifying selection [29], the BS-REL model [30] was used in a complementary approach to analyze individual protein-coding genes, as well the concatenated “supergene”. No significant results were obtained for individual genes, but A. charrua and N. furzeri showed a small proportion of COI sites with ω values above 5 (p-values <0.1, data not shown). About the concatenated data set, the tree branch leading to the clade composed by N. furzeri and A. panchax showed evidence of strong directional selection (ω >>5; corrected p-value <<0.01) at about 1 % of the sites (out of a total of 3787 codon sites analyzed; Additional file 11). The branch of N. furzeri showed about 5.5 % of the sites under strong selection (ω >>5) but in this case is not statistically significant at the chosen 0.05 level (corrected p-value <0.1).

In summary, we were unable to identify elevated rates of adaptive evolution in annual species as expected by their different life traits and lower G+C mitochondrial genome content. For the analyses, we worked with only six orthologs sequences in each case, and it is known that such small alignments use to have low power for phylogenetic methods of adaptive sequence evolution [27]. However, it appears that at least N. furzeri could exhibit molecular signatures of adaptive evolution in case of a proper taxon sampling scenario (i.e. a broaden taxa sampling within Rivulidae, Nothobranchiidae and Aplocheilidae).

Austrolebias charrua coding genes as phylogenetic molecular markers

Mitochondrial genes are some of the most popular and widely-utilized genetic markers in phylogenetic and phylogeographic studies [80]. As all genes reside in the same non-recombinant DNA molecule, their linked nature allows the assumption that they carry the same phylogenetic signal and share it with the entire mitogenome. However, as previously described, we found evidence of positive selection in the ATPase6, ATPase8, ND2 and Cytb genes and it is know that base composition changes at a different rate between lineages and genes affect adversely the efficiency of phylogenetic inference methods [54]. In this context, we assessed the performance of single mitochondrial protein-coding genes and a concatenated “supergene” at inferring evolutionary relationships throughout Cyprinodontiformes. This allowed us: i) to identify whether each protein-coding gene possesses the same phylogenetic signal and thus it is interchangeable; and ii) to distinguish how individual genes perform compared to the concatenated “supergene”. A phylogenetic inference was conducted for each gene and “supergene” in a maximum likelihood framework and each resulting gene tree visually compared to the reference species tree (Additional file 1). We did not try any mitochondrial specific DNA substitution models (e.g. GTR3 in [54]), thus we did not assessed putative inconsistencies due to base composition features.

Regarding whole coding sequence alignments, the inferred relative positioning of the lineages was in agreement with the reference species tree only for Cytb, ND1 and ND6 gene trees, while the concatenated “supergene” did not recover the reference species topology (Additional file 12). This is concordant with previous analyses using Cytb gene to resolve molecular phylogenetic relationships among species of the genus Austrolebias [68]. Moreover, most phylogeographic studies in different species groups within Austrolebias were based on this gene [8183]. When fast evolving third codon positions were removed from alignments, the “supergene” tree recovers our reference species phylogeny (Fig. 5) as well as the Cytb, ND4L and ND6 genes (Additional file 13). Therefore, third codon positions seem to be saturated in the aforementioned coding regions at intergeneric Cyprinodontiformes comparisons.

Fig. 5
figure5

Maximum likelihood analysis of a concatenated alignment of all 13 protein-coding genes “supergene” without the 3rd codon position. Branch support was conducted with a Bayesian-like transformation of approximate likelihood ratio test. The nucleotide substitution model was TVM+G determined according to the best suited sequence substitution model selected by ModelGenerator v0.85 [35], following Bayesian information criterion

While our taxonomic sampling is reduced and far from being an exhaustive and proper sample of cyprinodontiforms lineages, evidence points out differences in the phylogenetic informativeness level among the mitochondrial protein-coding genes. Even the concatenated gene set did not appear to perform consistently either. Considering cost-effectiveness trade-offs, Cytb appears as a proper molecular marker to be chosen when studying phylogenetic relationships, at inter- and intraspecific levels in Cyprinodontiformes, as was frequently used by Murphy and Collier [84, 85] and García et al. [6, 7] as well as in other highly related Rivulidae genera [86].

Conclusions

The mitogenome of A. charrua obtained in present work, represents the first one among Neotropical annual killifishes. It contains a gene arrangement and composition similar to most vertebrate mitogenomes. Among cyprinodontiforms, the mtDNA of A. charrua codes for the shortest tRNA-Cys without the D-loop arm, exhibits the lower G+C content, and the highest strand asymmetry. Our analyses showed that annual killifishes present an independent decrease in their mitogenome G+C content pointing to a change in the regimen of evolutionary forces that interact with the mitochondrial mutational bias.

Abbreviations

AA:

Amino acid

ATPase 6–8:

ATP synthase subunits 6–8

AUS4, 5, 7:

Austrolebias charrua individual 4, 5 and 7

COI-III:

Cytochrome c oxidase subunits I-III

CR:

Control region

Cytb:

Cytochrome b

GTR+I+G:

General time reversible model with invariable sites and gamma distribution

HKY+I+G:

Hasegawa-Kishino-Yano model with invariable sites and gamma distribution

ML:

Maximum likelihood

mtDNA:

Mitochondrial DNA

ND1-6, 4L:

NADH dehydrogenase subunits 1–6, 4L

PCR:

Polymerase chain reaction

ROS:

Reactive oxygen species

rRNA:

Ribosomal RNA

tRNA:

Transfer RNA

Ts:

Transition

Tv:

Transversion

References

  1. 1.

    Costa WJ. The South American annual killifish genus Austrolebias (Teleostei: Cyprinodontiformes: Rivulidae): phylogenetic relationships, descriptive morphology and taxonomic revision. Zootaxa. 2006;1213:1–162.

  2. 2.

    Vaz-Ferreira R, Sierra B. Caracteres etológicos genéricos y específicos en los peces del género Cynolebias Steindachner, 1876. Bol Soc Zool Urug. 1973;2:22–35.

  3. 3.

    Loureiro M, De Sá RO. Osteological analysis of the killifish genus Cynolebias (Cyprinodontiformes: Rivulidae). J Morphol. 1998;238:245–62.

  4. 4.

    Wourms JP. Annual fishes. Methods Dev Biol. N Y Thomas Crowell Co Pp. 1967. p. 123–37.

  5. 5.

    Wourms JP. The developmental biology of annual fishes. III. Pre-embryonic and embryonic diapause of variable duration in the eggs of annual fishes. J Exp Zool. 1972;182:389–414.

  6. 6.

    García G, Alvarez-Valin F, Gomez N. Mitochondrial genes: signals and noise in the phylogenetic reconstruction of the annual killifish genus Cynolebias (Cyprinodontiformes, Rivulidae). Biol J Linn Soc. 2002;76:49–59.

  7. 7.

    García G, Gutiérrez V, Ríos N, Turner B, Santiñaque F, López-Carro B, et al. Burst speciation processes and genomic expansion in the neotropical annual killifish genus Austrolebias (Cyprinodontiformes, Rivulidae). Genetica. 2014;142:87–98.

  8. 8.

    García G, Wlasiuk G, Lessa EP. High levels of mitochondrial cytochrome b divergence in annual killifishes of the genus Cynolebias (Cyprinodontiformes, Rivulidae). Zool J Linn Soc. 2000;129:93–110.

  9. 9.

    Martin AP, Palumbi SR. Body size, metabolic rate, generation time, and the molecular clock. Proc Natl Acad Sci. 1993;90:4087–91.

  10. 10.

    Errea A, Danulat E. Growth of the annual fish, Cynolebias viarius (Cyprinodontiformes), in the natural habitat compared to laboratory conditions. Environ Biol Fishes. 2001;61:261–8.

  11. 11.

    Lee J-S, Miya M, Lee Y-S, Kim CG, Park E-H, Aoki Y, et al. The complete DNA sequence of the mitochondrial genome of the self-fertilizing fish Rivulus marmoratus (Cyprinodontiformes, Rivulidae) and the first description of duplication of a control region in fish. Gene. 2001;280:1–7.

  12. 12.

    Rubinstein ND, Feldstein T, Shenkar N, Botero-Castro F, Griggio F, Mastrototaro F, et al. Deep sequencing of mixed total DNA without barcodes allows efficient assembly of highly plastic ascidian mitochondrial genomes. Genome Biol Evol. 2013;5:1185–99.

  13. 13.

    Xie Y, Wu G, Tang J, Luo R, Patterson J, Liu S, et al. SOAPdenovo-Trans: De novo transcriptome assembly with short RNA-Seq reads. Bioinformatics. 2014;30(12):1660–6.

  14. 14.

    Altschul SF, Madden TL, Schäffer AA, Zhang J, Zhang Z, Miller W, et al. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 1997;25:3389–402.

  15. 15.

    Iwasaki W, Fukunaga T, Isagozawa R, Yamada K, Maeda Y, Satoh TP, et al. MitoFish and MitoAnnotator: a mitochondrial genome database of fish with an accurate and automatic annotation pipeline. Mol Biol Evol. 2013;30:2531–40.

  16. 16.

    Rutherford K, Parkhill J, Crook J, Horsnell T, Rice P, Rajandream M-A, et al. Artemis: sequence visualization and annotation. Bioinformatics. 2000;16:944–5.

  17. 17.

    Jühling F, Pütz J, Bernt M, Donath A, Middendorf M, Florentz C, et al. Improved systematic tRNA gene annotation allows new insights into the evolution of mitochondrial tRNA structures and into the mechanisms of mitochondrial genome rearrangements. Nucleic Acids Res. 2012;40:2833–45.

  18. 18.

    Li H, Durbin R. Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics. 2009;25:1754–60.

  19. 19.

    Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078–9.

  20. 20.

    Team RC: R: A language and environment for statistical computing. Vienna, Austria; 2014. URL Httpwww R-Proj Org 2015.

  21. 21.

    Charif D, Lobry JR. SeqinR 1.0-2: a contributed package to the R project for statistical computing devoted to biological sequences retrieval and analysis. In: Structural approaches to sequence evolution. Berlin: Springer; 2007. p. 207–32.

  22. 22.

    Nenadic O, Greenacre M. Correspondence analysis in R, with two-and three-dimensional graphics: the ca package. 2007.

  23. 23.

    Pohl M, Milvertz FC, Meyer A, Vences M. Multigene phylogeny of cyprinodontiform fishes suggests continental radiations and a rogue taxon position of Pantanodon. Vertebr Zool. 2015;65:37–44.

  24. 24.

    Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004;32:1792–7.

  25. 25.

    Suyama M, Torrents D, Bork P. PAL2NAL: robust conversion of protein sequence alignments into the corresponding codon alignments. Nucleic Acids Res. 2006;34 suppl 2:W609–12.

  26. 26.

    Pupko T, Pe I, Shamir R, Graur D. A fast algorithm for joint reconstruction of ancestral amino acid sequences. Mol Biol Evol. 2000;17:890–6.

  27. 27.

    Kosakovsky Pond SL, Frost SDW. Not so different after all: a comparison of methods for detecting amino acid sites under selection. Mol Biol Evol. 2005;22:1208–22.

  28. 28.

    Delport W, Poon AF, Frost SD, Pond SLK. Datamonkey 2010: a suite of phylogenetic analysis tools for evolutionary biology. Bioinformatics. 2010;26:2455–7.

  29. 29.

    Murrell B, Wertheim JO, Moola S, Weighill T, Scheffler K, Pond SLK. Detecting individual sites subject to episodic diversifying selection. PLoS Genet. 2012;8, e1002764.

  30. 30.

    Pond SLK, Murrell B, Fourment M, Frost SD, Delport W, Scheffler K. A random effects branch-site model for detecting episodic diversifying selection. Mol Biol Evol. 2011;28(11):3033–43.

  31. 31.

    Guindon S, Dufayard J-F, Lefort V, Anisimova M, Hordijk W, Gascuel O. New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Syst Biol. 2010;59:307–21.

  32. 32.

    Guindon S, Gascuel O. A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. Syst Biol. 2003;52:696–704.

  33. 33.

    Hordijk W, Gascuel O. Improving the efficiency of SPR moves in phylogenetic tree search methods based on maximum likelihood. Bioinformatics. 2005;21:4338–47.

  34. 34.

    Anisimova M, Gil M, Dufayard J-F, Dessimoz C, Gascuel O. Survey of branch support methods demonstrates accuracy, power, and robustness of fast likelihood-based approximation schemes. Syst Biol. 2011;60(5):685–99.

  35. 35.

    Keane TM, Creevey CJ, Pentony MM, Naughton TJ, Mclnerney JO. Assessment of methods for amino acid matrix selection and their use on empirical data shows that ad hoc assumptions for choice of matrix are not justified. BMC Evol Biol. 2006;6:29.

  36. 36.

    Hasegawa M, Kishino H, Yano T. Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. J Mol Evol. 1985;22:160–74.

  37. 37.

    Tavaré S. Some probabilistic and statistical problems in the analysis of DNA sequences. Lect Math Life Sci. 1986;17:57–86.

  38. 38.

    Boore JL. Animal mitochondrial genomes. Nucleic Acids Res. 1999;27:1767–80.

  39. 39.

    Ojala D, Montoya J, Attardi G. tRNA punctuation model of RNA processing in human mitochondria. 1981.

  40. 40.

    Klimov PB, OConnor BM. Improved tRNA prediction in the American house dust mite reveals widespread occurrence of extremely short minimal tRNAs in acariform mites. BMC Genomics. 2009;10:598.

  41. 41.

    Jühling F, Pütz J, Florentz C, Stadler PF. Armless mitochondrial tRNAs in enoplea (nematoda). RNA Biol. 2012;9:1161.

  42. 42.

    Robinson JT, Thorvaldsdóttir H, Winckler W, Guttman M, Lander ES, Getz G, et al. Integrative genomics viewer. Nat Biotechnol. 2011;29:24–6.

  43. 43.

    Li M, Schönberg A, Schaefer M, Schroeder R, Nasidze I, Stoneking M. Detecting heteroplasmy from high-throughput sequencing of complete human mitochondrial DNA genomes. Am J Hum Genet. 2010;87:237–49.

  44. 44.

    Goto H, Dickins B, Afgan E, Paul IM, Taylor J, Makova KD, et al. Dynamics of mitochondrial heteroplasmy in three families investigated via a repeatable re-sequencing study. Genome Biol. 2011;12:R59.

  45. 45.

    Zhidkov I, Nagar T, Mishmar D, Rubin E. MitoBamAnnotator: a web-based tool for detecting and annotating heteroplasmy in human mitochondrial DNA sequences. Mitochondrion. 2011;11:924–8.

  46. 46.

    Avital G, Buchshtav M, Zhikov I, Tuval J, Dadon S, Rubin E, et al. Mitochondrial DNA heteroplasmy in diabetes and normal adults: role of acquired and inherited mutational patterns in twins. Hum Mol Genet. 2012;21(19):4214–24.

  47. 47.

    Samuels DC, Li C, Li B, Song Z, Torstenson E, Clay HB, et al. Recurrent tissue-specific mtDNA mutations are common in humans. PLoS Genet. 2013;9, e1003929.

  48. 48.

    Degli Esposti M, De Vries S, Crimi M, Ghelli A, Patarnello T, Meyer A. Mitochondrial cytochrome b: evolution and structure of the protein. Biochim Biophys Acta BBA-Bioenerg. 1993;1143:243–71.

  49. 49.

    Xia D, Yu C-A, Kim H, Xia J-Z, Kachurin AM, Zhang L, et al. Crystal structure of the cytochrome bc1 complex from bovine heart mitochondria. Science. 1997;277:60–6.

  50. 50.

    Hartmann N, Reichwald K, Wittig I, Dröse S, Schmeisser S, Lück C, et al. Mitochondrial DNA copy number and function decrease with age in the short‐lived fish Nothobranchius furzeri. Aging Cell. 2011;10:824–31.

  51. 51.

    Krzywinski M, Schein J, Birol İ, Connors J, Gascoyne R, Horsman D, et al. Circos: an information aesthetic for comparative genomics. Genome Res. 2009;19:1639–45.

  52. 52.

    Reichwald K, Lauber C, Nanda I, Kirschner J, Hartmann N, Schories S, et al. High tandem repeat content in the genome of the short-lived annual fish Nothobranchius furzeri: a new vertebrate model for aging research. Genome Biol. 2009;10:R16.

  53. 53.

    Castellana S, Vicario S, Saccone C. Evolutionary patterns of the mitochondrial genome in Metazoa: exploring the role of mutation and selection in mitochondrial protein–coding genes. Genome Biol Evol. 2011;3:1067–79.

  54. 54.

    Gibson A, Gowri-Shankar V, Higgs PG, Rattray M. A comprehensive analysis of mammalian mitochondrial genome base composition and improved phylogenetic methods. Mol Biol Evol. 2005;22:251–64.

  55. 55.

    Cardozo V. Tasa metabólica y excreción del nitrógeno en peces anuales Cynolebias viarius (Cyprinodontiformes), MSc. Dissertation. Thesis.. Uruguay: PEDECIBA. Facultad de Ciencias. UdelaR; 1999.

  56. 56.

    Martin AP. Metabolic rate and directional nucleotide substitution in animal mitochondrial DNA. Mol Biol Evol. 1995;12:1124–31.

  57. 57.

    Murphy M. How mitochondria produce reactive oxygen species. Biochem J. 2009;417:1–13.

  58. 58.

    Shokolenko I, Venediktova N, Bochkareva A, Wilson GL, Alexeyev MF. Oxidative stress induces degradation of mitochondrial DNA. Nucleic Acids Res. 2009;37:2539–48.

  59. 59.

    Smith DR. Updating our view of organelle genome nucleotide landscape. Front Genet. 2012;3.

  60. 60.

    Tanaka M, Ozawa T. Strand asymmetry in human mitochondrial DNA mutations. Genomics. 1994;22:327–35.

  61. 61.

    Bielawski JP, Gold JR. Mutation patterns of mitochondrial H-and L-strand DNA in closely related cyprinid fishes. Genetics. 2002;161:1589–97.

  62. 62.

    McKinney EA, Oliveira MT. Replicating animal mitochondrial DNA. Genet Mol Biol. 2013;36:308–15.

  63. 63.

    Min XJ, Hickey DA. DNA barcodes provide a quick preview of mitochondrial genome composition. PLoS One. 2007;2:e325.

  64. 64.

    Morton BR. Selection on the codon bias of chloroplast and cyanelle genes in different plant and algal lineages. J Mol Evol. 1998;46:449–59.

  65. 65.

    Jukes TH, Bhushan V. Silent nucleotide substitutions and G+ C content of some mitochondrial and bacterial genes. J Mol Evol. 1986;24:39–44.

  66. 66.

    Rocha EP, Danchin A. Base composition bias might result from competition for metabolic resources. Trends Genet. 2002;18:291–4.

  67. 67.

    McNab BK. Metabolic scaling: energy constraints on carnivore diet. Nature. 2000;407:584.

  68. 68.

    Arends A, McNab BK. The comparative energetics of “caviomorph”rodents. Comp Biochem Physiol A Mol Integr Physiol. 2001;130:105–22.

  69. 69.

    Bonaccorso F, McNab B. Standard energetics of leaf-nosed bats (Hipposideridae): its relationship to intermittent-and protracted-foraging tactics in bats and birds. J Comp Physiol B. 2003;173:43–53.

  70. 70.

    Da Fonseca RR, Johnson WE, O’Brien SJ, Ramos MJ, Antunes A. The adaptive evolution of the mammalian mitochondrial genome. BMC Genomics. 2008;9:119.

  71. 71.

    Shen Y-Y, Liang L, Zhu Z-H, Zhou W-P, Irwin DM, Zhang Y-P. Adaptive evolution of energy metabolism genes and the origin of flight in bats. Proc Natl Acad Sci. 2010;107:8666–71.

  72. 72.

    Hassanin A, Ropiquet A, Couloux A, Cruaud C. Evolution of the mitochondrial genome in mammals living at high altitude: new insights from a study of the tribe Caprini (Bovidae, Antilopinae). J Mol Evol. 2009;68:293–310.

  73. 73.

    Tomasco IH, Lessa EP. Two mitochondrial genes under episodic positive selection in subterranean octodontoid rodents. Gene. 2014;534:371–8.

  74. 74.

    Foote AD, Morin PA, Durban JW, Pitman RL, Wade P, Willerslev E, et al. Positive selection on the killer whale mitogenome. Biol Lett. 2011;7:116–8.

  75. 75.

    Garvin MR, Bielawski JP, Gharrett AJ. Positive Darwinian selection in the piston that powers proton pumps in complex I of the mitochondria of Pacific salmon. PLoS One. 2011;6.

  76. 76.

    Anisimova M, Kosiol C. Investigating protein-coding sequence evolution with probabilistic codon substitution models. Mol Biol Evol. 2009;26:255–71.

  77. 77.

    Bielawski JP. Detecting the signatures of adaptive evolution in protein‐coding genes. Curr Protoc Mol Biol. 2013;Chapter 19:Unit 19.1.

  78. 78.

    Menezes AN, Viana MC, Furtado C, Schrago CG, Seuánez HN. Positive selection along the evolution of primate mitogenomes. Mitochondrion. 2013;13:846–51.

  79. 79.

    Garvin MR, Bielawski JP, Sazanov LA, Gharrett AJ. Review and meta‐analysis of natural selection in mitochondrial complex I in metazoans. J Zool Syst Evol Res. 2015;53:1–17.

  80. 80.

    Avise JC. Phylogeography: the history and formation of species. Cambridge: Harvard University Press; 2000.

  81. 81.

    García G. Multiple simultaneous speciation in killifishes of the Cynolebias adloffi species complex (Cyprinodontiformes, Rivulidae) from phylogeography and chromosome data. J Zool Syst Evol Res. 2006;44:75–87.

  82. 82.

    García G, Loureiro M, Berois N, Arezo MJ, Casanova G, Clivio G, et al. Pattern of differentiation in the annual killifish genus Austrolebias (Cyprinodontiformes: Rivulidae) from a biosphere reserve site in South America: a multidisciplinary approach. Biol J Linn Soc. 2009;98:620–35.

  83. 83.

    Garcia G, Gutierrez V, Vergara J, Calvino P, Duarte A, Loureiro M. Patterns of population differentiation in annual killifishes from the Paraná–Uruguay–La Plata Basin: the role of vicariance and dispersal. J Biogeogr. 2012;39:1707–19.

  84. 84.

    Murphy WJ, Collier GE. Phylogenetic relationships within the aplocheiloid fish genus Rivulus (Cyprinodontiformes, Rivulidae): implications for Caribbean and Central American biogeography. Mol Biol Evol. 1996;13:642–9.

  85. 85.

    Murphy WJ, Collier GE. A molecular phylogeny for aplocheiloid fishes (Atherinomorpha, Cyprinodontiformes): the role of vicariance and the origins of annualism. Mol Biol Evol. 1997;14:790–9.

  86. 86.

    De León JLP, León G, Rodríguez R, Metcalfe CJ, Hernández D, Casane D, et al. Phylogeography of Cuban Rivulus: evidence for allopatric speciation and secondary dispersal across a marine barrier. Mol Phylogenet Evol. 2014;79:404–14.

  87. 87.

    Waterhouse AM, Procter JB, Martin DM, Clamp M, Barton GJ. Jalview Version 2—a multiple sequence alignment editor and analysis workbench. Bioinformatics. 2009;25:1189–91.

  88. 88.

    Darzentas N. Circoletto: visualizing sequence similarity with Circos. Bioinforma Oxf Engl. 2010;26:2620–1.

Download references

Acknowledgements

We thank anonymous reviewers for insightful comments on a previous version of the manuscript. VG was supported by PhD Fellowship (ANII, Agencia Nacional de Investigación e Innovación). The research of GG, VG and HN was also supported by the SNI (Sistema Nacional de Investigadores, Uruguay). The work was funded by the project FCE_6784 of the ANII to GG. This work was partially funded by FOCEM (MERCOSUR Structural Convergence Fund), COF 03/11. The authors are also grateful to the Japanese government for the donation of equipment.

Author information

Correspondence to Verónica Gutiérrez.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

GG conceived the study. VG conducted the sampling and carried out laboratory work. NR and HN performed all bioinformatic analyses. VG, NR and GG contributed to the writing and improving of the manuscript. All authors read and approved the final version.

Additional files

Additional file 1:

Phylogenetic relationships among some cyprinodontiforms families of fishes according to [ 23 ]. (PDF 122 kb)

Additional file 2;

Selected nucleotide substitution models after the third codon positions were removed from the codon alignments. (PDF 7 kb)

Additional file 3:

Dotplot of AUS7_control region obtained with zpicture ( http://zpicture.dcode.org/ ). Duplicated regions are of 180 bp (from 808 to 987) and 183 bp (from 970 to 1152) with a 91.67 % of identity. (PDF 59 kb)

Additional file 4:

Jalview visualization of a multiple sequence alignment of metazoan Cytb sequences [ 87 ]. A. charrua Cytb sequence was used as query for a blastp search of metazoan homolog sequences in the NCBI’s nr database. Downloaded amino acidic sequences were added to the cyprinodontiforms sequences used in this work and the whole dataset aligned with Muscle [24]. The alignment positions are coloured by percentage of identity and the column corresponding to A. charrua Cytb F 368 is highlighted with a red dash (column 375 in the present alignment). (PDF 1903 kb)

Additional file 5:

Base composition of mitochondrial coding-genes located in the H strand of all cyprinodontiforms compared in present work. Values are discriminated according to first, second and third codon position. (PDF 16 kb)

Additional file 6:

Percentage of the G+C content considering the 12 protein-coding genes located in the H strand. In green, values for the G+C content of the sequence reconstruction of ancestral nodes with a joint-likelihood method in the codon-state space. In red, G+C content of A. charrua and N. furzeri and in black, values of the remaining cyprinodontiforms mitogenomes analyzed. (PDF 472 kb)

Additional file 7:

Representation of the first two dimensions of the correspondence analysis (COA) performed on the global codon usage of the six species of cyprinodontiforms analyzed in this study. Variables (codons) and cases (species) are plotted together. Blue dots represent the species: AC (A. charrua), KM (K. marmoratus), NF (N. furzeri), AP (A. panchax), CR (C. rubrofluviatilis) and FO (F. olivaceus). Red triangles correspond to codons. (PDF 668 kb)

Additional file 8:

Correspondence analysis of per gene codon usage considering only A. charrua (AC), K. marmoratus (KM) and N. furzeri (NF). Blue dots correspond to per gene codon usage of A. charrua, red to N. furzeri and black to per gene codon usage of K. marmoratus. (PDF 545 kb)

Additional file 9:

Representation of the first two dimensions of the correspondence analysis (COA) performed on the amino acid frequencies of the six species of cyprinodontiforms analyzed. Variables (amino acids) and cases (species) are plotted together. The amino acids are represented in the usual one-letter code (red triangles). Blue dots represent the species: AC (A. charrua), KM (K. marmoratus), NF (N. furzeri), AP (A. panchax), CR (C. rubrofluviatilis) and FO (F. olivaceus). (PDF 567 kb)

Additional file 10:

Non-synonymous to synonymous rate ratio of all mitochondrial protein-coding genes. (PDF 410 kb)

Additional file 11:

BS-REL results for concatenated data set. Node 1, 2 and 5 as in Additional file 1. (PDF 479 kb)

Additional file 12:

Maximum likelihood analysis of all protein-coding genes. a) ATPase6; b) ATPase8; c) COI; d) COII; e) COIII; f) Cytb; g) ND1; h) ND2; i) ND3; j) ND4; k) ND4L; l) ND5; m) ND6 and n) concatenated “supergene”. (PDF 298 kb)

Additional file 13:

Maximum likelihood analysis of all protein-coding genes without the 3rd codon position. a) ATPase6; b) ATPase8; c) COI; d) COII; e) COIII; f) Cytb; g) ND1; h) ND2; i) ND3; j) ND4; k) ND4L; l) ND5 and m) ND6. (PDF 619 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Keywords

  • High-throughput sequencing
  • Austrolebias charrua
  • Cyprinodontiformes
  • Comparative mitogenomics
  • Molecular markers