BACKGROUND
the neotropical genus austrolebias  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  <cit> .
it includes  <dig> known species  highly variable in morphology and behavior  <cit> .
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  <cit> .
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  <cit> .
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  <cit> .
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  <cit> .
as proposed by martin and palumbi  <cit> , 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  <cit> .
in the present work we describe for the first time the complete mitochondrial genome of a. charrua , 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  <cit> .
to date,  <dig> mtdna sequences  are deposited in the mitochondrial genome database of fish , of which  <dig> 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  were analyzed in the present work.
they belonged to the same f <dig> generation obtained at the lab while the parents were collected in a temporal pond located in “la coronilla” .
aus <dig> and aus <dig> were young individuals of 3 months old, while aus <dig> 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  of freshly sacrificed animals by an over-exposure to a solution of 1‰ 2-phenoxyethanol .
the sampling protocol was approved by the cnea  from uruguay.
sequencing
sequencing libraries were generated using nextera xt kit  from 1 ng of input dna.
tagmentation of genomic dna, pcr amplification  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  using the  <dig> bioanalyzer  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 .
two sets of specific primers: acr-1f , acr-1r  and acr-2f , acr-2r  were created using primer <dig> v.  <dig> .
<dig> software.
the pcr was carried out in a 10 μl total volume of 10× buffer, mgcl <dig>  <dig>  mm, dntps  <dig>  mm,  <dig>  mm of each primer,  <dig>  u of taq dna polymerase  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  <dig> dna sequencer .
mitogenome assembly and annotation
raw reads were quality filtered and trimmed to a final length of at least 50 bp.
as indicated previously  <cit> , 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  <cit> , contigs were assembled using soapdenovo-trans  <cit>  with k-mer length of  <dig> and remaining default options.
as both nuclear and mitochondrial contigs were assembled with this strategy, a local blastn  <cit>  was ran against a database of cyprinodontiforms mitogenomes downloaded from mitofish  <cit>  , to recover those contigs of mitochondrial origin.
selected contigs were annotated with mitoannotator  <cit> , a highly accurate and automated pipeline specific for fish mitogenomes, followed by manual curation in artemis  <cit> .
transfer rna genes were annotated based on their potential cloverleaf secondary structure using mitfi  <cit> , a tool for accurately locating trna genes within mitogenomic sequences .
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  and among siblings.
this was performed with a combination of burrows-wheeler aligner   <cit>  and samtools   <cit> .
of note, both the read mapping quality  and base quality  were required to have a minimum value of  <dig> 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  <cit>  using seqinr  <cit>  and ca packages  <cit> .
besides whole mitogenome analyses, calculations were performed on multiple alignments of the  <dig> protein-coding genes present in the h strand .
in these comparative analyses, the following cyprinodontiformes were included: kryptolebias marmoratus , nothobranchius furzeri , aplocheilus panchax , cyprinodon rubrofluviatilis  and fundulus olivaceus .
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.
<cit>  topology is shown in additional file  <dig>  multiple nucleotide alignments were built guided by previous alignments of amino acid sequences using muscle  <cit>  and pal2nal  <cit> .
for each alignment, ancestral sequence reconstruction was performed with a joint-likelihood reconstruction method in the codon-state space  <cit> .
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  <cit> .
uploaded alignments were matched with the reference species tree  and f. olivaceus fixed as outgroup lineage if required.
global dn/ds values were estimated with slac  <cit> .
meme  <cit>  and bs-rel  <cit>  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 α =  <dig>  level; and iii) also showed an empirical bayes factor > <dig>  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  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 v <dig>   <cit> .
also, a “supergene” phylogeny was based on a concatenated alignment of all  <dig> protein-coding genes.
in each case, five random starting trees were used to search the tree space by nni  <cit>  and spr  <cit>  methods.
survey of branch support was conducted with a bayesian-like transformation of approximate likelihood ratio test .
nucleotide substitution model parameters were determined according to the best suited sequence substitution model selected by modelgenerator v <dig>   <cit> , following bayesian information criterion.
for most genes, the hky+i+g model  <cit>  was the chosen one, while the gtr+i+g  <cit>  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 .
all inferences were recalculated after third codon positions were removed from the alignments .
RESULTS
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  <cit> , almost the whole mitochondrial genome was reconstructed in only one contig in the samples aus <dig> and aus <dig>  while two contigs were assembled in case of aus <dig> .
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  <cit> .
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 aus <dig> and aus <dig>  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 .table  <dig> illumina sequencing and soapdenovo-trans assembly
#contigs with mitochondrial sequencesa
asimilarity was determined by blastn search against a database of mitogenomes of cyprinodontiformes; the four contigs showed above 73 % similarity and e-value close to zero
the de novo assembled mitochondrial contig of aus <dig> was manually curated and edited using this additional sequence information.
the resulting complete mitochondrial sequence of  <dig>  bp  was deposited in the genebank  and considered our reference mitogenome for a. charrua .
regarding gene content and order, it exhibits the typical pattern of vertebrate mitogenomes  <cit> .
it contains the full complement of  <dig> protein-coding genes,  <dig> trna genes, two rrna genes  and one non-coding control region.
eight trnas , trna-glu and trna-pro) and the protein-coding gene nd <dig> were encoded by the light  strand, whilst the remaining genes were encoded by the heavy  strand .
atg codons initiate  <dig> of the protein-coding genes, coi being the exception .
complete taa stop codons are present in coi, coiii, atpase <dig>  nd <dig>  nd4l, nd <dig> and cytb, while a complete tag codon is present in nd <dig>  remaining protein genes end at incomplete codons , which are converted by polyadenylation into taa after transcription and processing  <cit> .
of  <dig> transfer rna genes,  <dig> were within the size range of 66–74 bp , 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  and in the predicted secondary structure lacks the d-loop arm .
while d-armless trna-cys reports do not exist for fishes  <cit> , several independent occurrences do exist in metazoans, including mammals.
atypical trnas  would be functional  <cit> .
to note, in k. marmoratus, trna-cys is 65 bp in length and has the typical secondary structure .
further analysis of new austrolebias taxa mitogenomes could clarify whether this d-armless trna-cys is a genus or species-level feature.fig.
<dig> map of the mitochondrial genome of austrolebias charrua .
visual representation created by circos  <cit>  where the innermost circle represents g+c content per every 5 bp  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
nd1
nd2
coi
coii
coiii
nd3
nd4l
nd4
nd5
nd6
cytb
standard abbreviations are used for protein-coding genes and three letter abbreviations are given for trnas and rrnas.
intergenic nucleotide refers to the nucleotide distance between pairs of adjacent genes
fig.
<dig> putative secondary structure of the trna-cys encoded in the mtdna of  a. charrua and  k. marmoratus.
the trnas were annotated using mitfi tool run through mitoannotator  <cit>
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  and mpileup files generated as mentioned in methods.
r scripts  <cit>  and careful inspection with igv  <cit>  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 ; 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  and the features  of our data.
a total of  <dig> variants were detected in protein-coding genes of the mitogenome of the full siblings analyzed, five in the mtdna of aus <dig> and seven in aus <dig>  the diversity of changes and their putative effects on the coding sequence are shown in table  <dig>  the transversion c to a identified in the coi gene  was present in both samples, while the transition t to c detected in the cytb gene  was present in all full siblings, including the reference genome .
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  <cit> .
to note, above 28 % of the reads mapping to position  <dig>  carried the alternative nucleotide cytosine.
this heteroplasmic position corresponds to the codon  <dig> of the cytb locus, and despite of being a transition it produces a non-synonymous change from phe to leu in the translated sequence .
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 .
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 bc <dig> complex  <cit> .table  <dig> nucleotide variants in the mitogenomes of three full siblings of a. charrua
reads were mapped against the reference genome aus <dig>  the number of reads carrying the variant is shown together with total number the reads mapping to each position
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 .
kryptolebias marmoratus displays an additional control region located between trna-leu and nd <dig>  <cit>  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  <cit> .
a circos plot  <cit>  highlighting synteny and identity level between a. charrua and each of the five mitogenomes used for comparison is shown in fig.
<dig>  the most notorious feature in the plot is the lack of similarity between a. charrua and the last region of n. furzeri mitochondrial chromosome , despite homologous regions can be identified in more distant species as c. rubrofluviatilis and f. olivaceous .
considering that the mitogenome of n. furzeri was built as part of a whole genome shotgun approach resulting in a narrow  <dig>  mb draft assembly   <cit>  a putative miss-assembly should be discarded before further comparative analyses.table  <dig> nucleotide composition of all cyprinodontiformes compared in present work
austrolebias charrua
b
nothobranchius furzeri
b
kryptolebias marmoratus
aplocheilus panchax
cyprinodon rubrofluviatilis
fundulus olivaceus
a, c, g, t, gcskew and atskew refer to whole mitogenome content analysis.
gc <dig>  gc <dig>  gc <dig> and coding gc refer to the analysis of the mitochondrial coding-genes located in the h strand.
see additional file  <dig> for discriminated base composition values at first, second and third codon position of these  <dig> coding-genes located in the h strand
a
k. marmoratus cr1
bannual species
fig.
<dig> circos plot mapping the synteny and identity level of a. charrua mitogenome against each of the other five mitogenomes studied , k. marmoratus  and n. furzeri .
idiograms and red ribbons represent the similarity after pairwise blastn searches.
in a. charrua idiogram, the  <dig> coding genes are coloured in blue.
in all idiograms, control region are coloured in black.
in , n. furzeri unique features  are shown in orange.
the figure was produced using circos software  <cit> , a tool for graphical representation of genome data, through the circoletto implementation  <cit> .
ac , km , nf , ap , cr  and fo
table  <dig> 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 /) and atskew /) are negative, pointing to strand asymmetry in the nucleotide composition .
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 .
in fact, this relatively strong negative value is highly unusual in actinopterygii and even chordata  <cit> .
given the differences in length of the annotated control regions  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  <dig> protein-coding genes located in the h strand.
the g+c content is again lower in both annual killifish, being the difference  larger in the first and third codon positions .
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  <cit> .
estimation of g+c content at the internal nodes of our reference topology  confirms that the lower g+c content is a recently and independently derivative character state in both annual lineages .
correspondence analyses of nucleotide composition of the coding-genes of the h strand  and on global codon usage  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  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 , second axis correlates to atskew .
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 a <dig> and a <dig> while the second appear richer in t <dig> and t <dig>  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 .fig.
<dig> 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 , km , nf , ap , cr  and fo .
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 , a fairly common characteristic  <cit> .
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  <cit> .
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 .
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  <cit> .
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  <cit> .
while this model of mtdna replication has been challenged in recent years , asymmetric trends in base composition have been confirmed  <cit> .
other various hypotheses have tried to explain the high adenine and thymine  content of almost all mitochondrial dnas  <cit> .
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  <cit> .
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  <cit> .
as expected, codon usage pattern also translates into differences in amino acid usage .
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 .
in a landscape of widespread strong purifying selection due to functional constraints, signals of positive selection have been reported for individual codon sites  of taxa known to have high energetic demanding lifestyles such as flying in bats  <cit>  and unusual oxygen availability such as high altitude caprini and subterranean octodontoid rodents  <cit> , among other conditions  <cit> .
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  and positive selection can be inferred whenever the estimated ratio dn/ds significantly exceeds one  <cit> .
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 , greatly improving the power to detect positive selection.
more recently a mixed effects model of evolution , a highly sensitive branch-site random effects phylogenetic method, was presented  <cit> .
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  <dig> protein-coding genes are under strong purifying selection, with dn/ds values below  <dig>  .
the atpase <dig>  atpase <dig> and nd <dig> genes showed the highest dn/ds values, while components of the cytochrome oxidase complex  showed the lowest rates.
these results are in perfect agreement with a previous comparative analysis in metazoans  <cit> .
the analysis conducted with meme  <cit>  inferred four sites where some branches experienced a significant greater non-synonymous than synonymous substitution rate.
these sites are: codon  <dig> of atpase6; codon  <dig> of atpase8; codon  <dig> of cytb and codon  <dig> of nd <dig> .
molecular footprints of positive selection have been previously associated to mitochondrial nd <dig> in at least one work  <cit>  and signatures of adaptation have been repeatedly detected in cytb, atpase <dig> and atpase <dig>  <cit> .
noteworthy, in three of these sites the episodic directional selection would have happened in the ancestral branch leading to rivulids , while in the cytb site it would have occurred in the tree branch leading to n. furzeri and a. panchax .
we did not find any significant evidence of sites showing episodic selection at both annual lineages.
however, as meme authors denote  <cit> , 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  <dig> codon sites and branches under episodic and pervasive selection inferred by mixed effects model of evolution
bef bayes empirical factor, syn subs and non-syn subs number of synonymous and non-synonymous substitutions inferred for the codon site, respectively
as explained, one cannot simultaneously infer the site and the branch subject to diversifying selection  <cit> , the bs-rel model  <cit>  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  <dig> .
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  at about 1 % of the sites .
the branch of n. furzeri showed about  <dig>  % of the sites under strong selection  but in this case is not statistically significant at the chosen  <dig>  level .
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  <cit> .
however, it appears that at least n. furzeri could exhibit molecular signatures of adaptive evolution in case of a proper taxon sampling scenario .
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  <cit> .
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 atpase <dig>  atpase <dig>  nd <dig> 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  <cit> .
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 .
we did not try any mitochondrial specific dna substitution models , 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, nd <dig> and nd <dig> gene trees, while the concatenated “supergene” did not recover the reference species topology .
this is concordant with previous analyses using cytb gene to resolve molecular phylogenetic relationships among species of the genus austrolebias .
moreover, most phylogeographic studies in different species groups within austrolebias were based on this gene .
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 nd <dig> genes .
therefore, third codon positions seem to be saturated in the aforementioned coding regions at intergeneric cyprinodontiformes comparisons.fig.
<dig> maximum likelihood analysis of a concatenated alignment of all  <dig> 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 v <dig>   <cit> , 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  <cit>  and garcía et al.
<cit>  as well as in other highly related rivulidae genera  <cit> .
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.
additional files
additional file 1:
phylogenetic relationships among some cyprinodontiforms families of fishes according to [
23
].
additional file 2;
selected nucleotide substitution models after the third codon positions were removed from the codon alignments.
additional file 3:
dotplot of aus7_control region obtained with zpicture .
duplicated regions are of 180 bp  and 183 bp  with a  <dig>  % of identity.
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  <cit> .
the alignment positions are coloured by percentage of identity and the column corresponding to a. charrua cytb f  <dig> is highlighted with a red dash .
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.
additional file 6:
percentage of the g+c content considering the  <dig> 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.
additional file 7:
representation of the first two dimensions of the correspondence analysis  performed on the global codon usage of the six species of cyprinodontiforms analyzed in this study.
variables  and cases  are plotted together.
blue dots represent the species: ac , km , nf , ap , cr  and fo .
red triangles correspond to codons.
additional file 8:
correspondence analysis of per gene codon usage considering only
a. charrua
,
k. marmoratus
and
n. furzeri
.
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.
additional file 9:
representation of the first two dimensions of the correspondence analysis  performed on the amino acid frequencies of the six species of cyprinodontiforms analyzed.
variables  and cases  are plotted together.
the amino acids are represented in the usual one-letter code .
blue dots represent the species: ac , km , nf , ap , cr  and fo .
additional file 10:
non-synonymous to synonymous rate ratio of all mitochondrial protein-coding genes.
additional file 11:
bs-rel results for concatenated data set.
node  <dig>   <dig> and  <dig> as in additional file  <dig>
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) nd <dig> and n) concatenated “supergene”.
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) nd <dig> and m) nd <dig>
abbreviations
aaamino acid
atpase 6–8atp synthase subunits 6–8
aus <dig>   <dig>  7austrolebias charrua individual  <dig>   <dig> and 7
coi-iiicytochrome c oxidase subunits i-iii
crcontrol region
cytbcytochrome b
gtr+i+ggeneral time reversible model with invariable sites and gamma distribution
hky+i+ghasegawa-kishino-yano model with invariable sites and gamma distribution
mlmaximum likelihood
mtdnamitochondrial dna
nd1- <dig>  4lnadh dehydrogenase subunits 1– <dig>  4l
pcrpolymerase chain reaction
rosreactive oxygen species
rrnaribosomal rna
trnatransfer rna
tstransition
tvtransversion
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.

 the neotropical genus austrolebias 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. it includes known species highly variable in morphology and behavior. 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. 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. 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. 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. as proposed by martin and palumbi, 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. in the present work we describe for the first time the complete mitochondrial genome of a. charrua, 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. to date, mtdna sequences are deposited in the mitochondrial genome database of fish, of which 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 were analyzed in the present work. they belonged to the same f generation obtained at the lab while the parents were collected in a temporal pond located in “la coronilla”. aus and aus were young individuals of 3 months old, while aus 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 of freshly sacrificed animals by an over-exposure to a solution of 1‰ 2-phenoxyethanol. the sampling protocol was approved by the cnea from uruguay. sequencing sequencing libraries were generated using nextera xt kit from 1 ng of input dna. tagmentation of genomic dna, pcr amplification 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 using the bioanalyzer 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. two sets of specific primers: acr-1f, acr-1r and acr-2f, acr-2r were created using primer v.. software. the pcr was carried out in a 10 μl total volume of 10× buffer, mgcl  mm, dntps  mm,  mm of each primer, u of taq dna polymerase 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 dna sequencer. mitogenome assembly and annotation raw reads were quality filtered and trimmed to a final length of at least 50 bp. as indicated previously, 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, contigs were assembled using soapdenovo-trans with k-mer length of and remaining default options. as both nuclear and mitochondrial contigs were assembled with this strategy, a local blastn was ran against a database of cyprinodontiforms mitogenomes downloaded from mitofish, to recover those contigs of mitochondrial origin. selected contigs were annotated with mitoannotator, a highly accurate and automated pipeline specific for fish mitogenomes, followed by manual curation in artemis. transfer rna genes were annotated based on their potential cloverleaf secondary structure using mitfi, a tool for accurately locating trna genes within mitogenomic sequences. 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 and among siblings. this was performed with a combination of burrows-wheeler aligner and samtools. of note, both the read mapping quality and base quality were required to have a minimum value of 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 using seqinr and ca packages. besides whole mitogenome analyses, calculations were performed on multiple alignments of the protein-coding genes present in the h strand. in these comparative analyses, the following cyprinodontiformes were included: kryptolebias marmoratus, nothobranchius furzeri, aplocheilus panchax, cyprinodon rubrofluviatilis and fundulus olivaceus. 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. topology is shown in additional file multiple nucleotide alignments were built guided by previous alignments of amino acid sequences using muscle and pal2nal. for each alignment, ancestral sequence reconstruction was performed with a joint-likelihood reconstruction method in the codon-state space. 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. uploaded alignments were matched with the reference species tree and f. olivaceus fixed as outgroup lineage if required. global dn/ds values were estimated with slac. meme and bs-rel 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 α =  level; and iii) also showed an empirical bayes factor > 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 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 v. also, a “supergene” phylogeny was based on a concatenated alignment of all protein-coding genes. in each case, five random starting trees were used to search the tree space by nni and spr methods. survey of branch support was conducted with a bayesian-like transformation of approximate likelihood ratio test. nucleotide substitution model parameters were determined according to the best suited sequence substitution model selected by modelgenerator v, following bayesian information criterion. for most genes, the hky+i+g model was the chosen one, while the gtr+i+g 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. all inferences were recalculated after third codon positions were removed from the alignments.  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, almost the whole mitochondrial genome was reconstructed in only one contig in the samples aus and aus while two contigs were assembled in case of aus. 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. 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 aus and aus 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.table illumina sequencing and soapdenovo-trans assembly #contigs with mitochondrial sequencesa asimilarity was determined by blastn search against a database of mitogenomes of cyprinodontiformes; the four contigs showed above 73 % similarity and e-value close to zero the de novo assembled mitochondrial contig of aus was manually curated and edited using this additional sequence information. the resulting complete mitochondrial sequence of  bp was deposited in the genebank and considered our reference mitogenome for a. charrua. regarding gene content and order, it exhibits the typical pattern of vertebrate mitogenomes. it contains the full complement of protein-coding genes, trna genes, two rrna genes and one non-coding control region. eight trnas, trna-glu and trna-pro) and the protein-coding gene nd were encoded by the light strand, whilst the remaining genes were encoded by the heavy strand. atg codons initiate of the protein-coding genes, coi being the exception. complete taa stop codons are present in coi, coiii, atpase nd nd4l, nd and cytb, while a complete tag codon is present in nd remaining protein genes end at incomplete codons, which are converted by polyadenylation into taa after transcription and processing. of transfer rna genes, were within the size range of 66–74 bp, 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 and in the predicted secondary structure lacks the d-loop arm. while d-armless trna-cys reports do not exist for fishes, several independent occurrences do exist in metazoans, including mammals. atypical trnas would be functional. to note, in k. marmoratus, trna-cys is 65 bp in length and has the typical secondary structure. further analysis of new austrolebias taxa mitogenomes could clarify whether this d-armless trna-cys is a genus or species-level feature.fig. map of the mitochondrial genome of austrolebias charrua. visual representation created by circos where the innermost circle represents g+c content per every 5 bp 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 nd1 nd2 coi coii coiii nd3 nd4l nd4 nd5 nd6 cytb standard abbreviations are used for protein-coding genes and three letter abbreviations are given for trnas and rrnas. intergenic nucleotide refers to the nucleotide distance between pairs of adjacent genes fig. putative secondary structure of the trna-cys encoded in the mtdna of a. charrua and k. marmoratus. the trnas were annotated using mitfi tool run through mitoannotator 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 and mpileup files generated as mentioned in methods. r scripts and careful inspection with igv 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; 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 and the features of our data. a total of variants were detected in protein-coding genes of the mitogenome of the full siblings analyzed, five in the mtdna of aus and seven in aus the diversity of changes and their putative effects on the coding sequence are shown in table  the transversion c to a identified in the coi gene was present in both samples, while the transition t to c detected in the cytb gene was present in all full siblings, including the reference genome. 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. to note, above 28 % of the reads mapping to position carried the alternative nucleotide cytosine. this heteroplasmic position corresponds to the codon of the cytb locus, and despite of being a transition it produces a non-synonymous change from phe to leu in the translated sequence. 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. 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 bc complex.table nucleotide variants in the mitogenomes of three full siblings of a. charrua reads were mapped against the reference genome aus the number of reads carrying the variant is shown together with total number the reads mapping to each position 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. kryptolebias marmoratus displays an additional control region located between trna-leu and nd 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. a circos plot highlighting synteny and identity level between a. charrua and each of the five mitogenomes used for comparison is shown in fig.  the most notorious feature in the plot is the lack of similarity between a. charrua and the last region of n. furzeri mitochondrial chromosome, despite homologous regions can be identified in more distant species as c. rubrofluviatilis and f. olivaceous. considering that the mitogenome of n. furzeri was built as part of a whole genome shotgun approach resulting in a narrow  mb draft assembly a putative miss-assembly should be discarded before further comparative analyses.table nucleotide composition of all cyprinodontiformes compared in present work austrolebias charrua b nothobranchius furzeri b kryptolebias marmoratus aplocheilus panchax cyprinodon rubrofluviatilis fundulus olivaceus a, c, g, t, gcskew and atskew refer to whole mitogenome content analysis. gc gc gc and coding gc refer to the analysis of the mitochondrial coding-genes located in the h strand. see additional file for discriminated base composition values at first, second and third codon position of these coding-genes located in the h strand a k. marmoratus cr1 bannual species fig. circos plot mapping the synteny and identity level of a. charrua mitogenome against each of the other five mitogenomes studied, k. marmoratus and n. furzeri. idiograms and red ribbons represent the similarity after pairwise blastn searches. in a. charrua idiogram, the coding genes are coloured in blue. in all idiograms, control region are coloured in black. in, n. furzeri unique features are shown in orange. the figure was produced using circos software, a tool for graphical representation of genome data, through the circoletto implementation. ac, km, nf, ap, cr and fo table  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 /) and atskew /) are negative, pointing to strand asymmetry in the nucleotide composition. 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. in fact, this relatively strong negative value is highly unusual in actinopterygii and even chordata. given the differences in length of the annotated control regions 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 protein-coding genes located in the h strand. the g+c content is again lower in both annual killifish, being the difference larger in the first and third codon positions. 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. estimation of g+c content at the internal nodes of our reference topology confirms that the lower g+c content is a recently and independently derivative character state in both annual lineages. correspondence analyses of nucleotide composition of the coding-genes of the h strand and on global codon usage 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 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, second axis correlates to atskew. 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 a and a while the second appear richer in t and t 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.fig. 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, km, nf, ap, cr and fo. 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, a fairly common characteristic. 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. 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. 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. 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. while this model of mtdna replication has been challenged in recent years, asymmetric trends in base composition have been confirmed. other various hypotheses have tried to explain the high adenine and thymine content of almost all mitochondrial dnas. 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. 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. as expected, codon usage pattern also translates into differences in amino acid usage. 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. in a landscape of widespread strong purifying selection due to functional constraints, signals of positive selection have been reported for individual codon sites of taxa known to have high energetic demanding lifestyles such as flying in bats and unusual oxygen availability such as high altitude caprini and subterranean octodontoid rodents, among other conditions. 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 and positive selection can be inferred whenever the estimated ratio dn/ds significantly exceeds one. 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, greatly improving the power to detect positive selection. more recently a mixed effects model of evolution, a highly sensitive branch-site random effects phylogenetic method, was presented. 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 protein-coding genes are under strong purifying selection, with dn/ds values below. the atpase atpase and nd genes showed the highest dn/ds values, while components of the cytochrome oxidase complex showed the lowest rates. these results are in perfect agreement with a previous comparative analysis in metazoans. the analysis conducted with meme inferred four sites where some branches experienced a significant greater non-synonymous than synonymous substitution rate. these sites are: codon of atpase6; codon of atpase8; codon of cytb and codon of nd. molecular footprints of positive selection have been previously associated to mitochondrial nd in at least one work and signatures of adaptation have been repeatedly detected in cytb, atpase and atpase. noteworthy, in three of these sites the episodic directional selection would have happened in the ancestral branch leading to rivulids, while in the cytb site it would have occurred in the tree branch leading to n. furzeri and a. panchax. we did not find any significant evidence of sites showing episodic selection at both annual lineages. however, as meme authors denote, 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 codon sites and branches under episodic and pervasive selection inferred by mixed effects model of evolution bef bayes empirical factor, syn subs and non-syn subs number of synonymous and non-synonymous substitutions inferred for the codon site, respectively as explained, one cannot simultaneously infer the site and the branch subject to diversifying selection, the bs-rel model 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. 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 at about 1 % of the sites. the branch of n. furzeri showed about  % of the sites under strong selection but in this case is not statistically significant at the chosen level. 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. however, it appears that at least n. furzeri could exhibit molecular signatures of adaptive evolution in case of a proper taxon sampling scenario. 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. 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 atpase atpase nd 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. 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. we did not try any mitochondrial specific dna substitution models, 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, nd and nd gene trees, while the concatenated “supergene” did not recover the reference species topology. this is concordant with previous analyses using cytb gene to resolve molecular phylogenetic relationships among species of the genus austrolebias. moreover, most phylogeographic studies in different species groups within austrolebias were based on this gene. 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 nd genes. therefore, third codon positions seem to be saturated in the aforementioned coding regions at intergeneric cyprinodontiformes comparisons.fig. maximum likelihood analysis of a concatenated alignment of all 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 v, 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 and garcía et al. as well as in other highly related rivulidae genera.  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. additional files additional file 1: phylogenetic relationships among some cyprinodontiforms families of fishes according to [ 23 ]. additional file 2; selected nucleotide substitution models after the third codon positions were removed from the codon alignments. additional file 3: dotplot of aus7_control region obtained with zpicture. duplicated regions are of 180 bp and 183 bp with a  % of identity. 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. the alignment positions are coloured by percentage of identity and the column corresponding to a. charrua cytb f is highlighted with a red dash. 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. additional file 6: percentage of the g+c content considering the 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. additional file 7: representation of the first two dimensions of the correspondence analysis performed on the global codon usage of the six species of cyprinodontiforms analyzed in this study. variables and cases are plotted together. blue dots represent the species: ac, km, nf, ap, cr and fo. red triangles correspond to codons. additional file 8: correspondence analysis of per gene codon usage considering only a. charrua, k. marmoratus and n. furzeri. 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. additional file 9: representation of the first two dimensions of the correspondence analysis performed on the amino acid frequencies of the six species of cyprinodontiforms analyzed. variables and cases are plotted together. the amino acids are represented in the usual one-letter code. blue dots represent the species: ac, km, nf, ap, cr and fo. additional file 10: non-synonymous to synonymous rate ratio of all mitochondrial protein-coding genes. additional file 11: bs-rel results for concatenated data set. node and as in additional file 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) nd and n) concatenated “supergene”. 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) nd and m) nd abbreviations aaamino acid atpase 6–8atp synthase subunits 6–8 aus 7austrolebias charrua individual and 7 coi-iiicytochrome c oxidase subunits i-iii crcontrol region cytbcytochrome b gtr+i+ggeneral time reversible model with invariable sites and gamma distribution hky+i+ghasegawa-kishino-yano model with invariable sites and gamma distribution mlmaximum likelihood mtdnamitochondrial dna nd1- 4lnadh dehydrogenase subunits 1– 4l pcrpolymerase chain reaction rosreactive oxygen species rrnaribosomal rna trnatransfer rna tstransition tvtransversion 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. 
regarding gene content and order, the circular mitogenome of a. charrua presents the typical pattern of vertebrate mitogenomes. all protein-coding genes are under strong purifying selection and we did not find any significant evidence of nucleotide sites showing episodic selection at annual lineages. 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. among teleosts, the south american genus austrolebias includes taxa of annual fishes divided into five different species groups. given the annual life-style, their mitogenomes could be under different selective pressures. among cyprinodontiforms, the mitogenome of a. charrua exhibits the lowest g+c content 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 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. 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 nd genes. therefore, third codon positions seem to be saturated in the aforementioned coding regions at intergeneric cyprinodontiformes comparisons. notably, the trna-cys was only 57 bp in length and lacks the d-loop arm. 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. 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. it contains the full complement of proteins-coding genes, trna, rrna and one non-coding control region.
 among teleosts, the south american genus austrolebias includes 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.  regarding gene content and order, the circular mitogenome of a. charrua presents the typical pattern of vertebrate mitogenomes. it contains the full complement of proteins-coding genes, trna, 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 variable sites in seven protein-coding genes. among cyprinodontiforms, the mitogenome of a. charrua exhibits the lowest g+c content 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 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 protein-coding genes are under strong purifying selection and we did not find any significant evidence of nucleotide sites showing episodic selection 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 nd genes. therefore, third codon positions seem to be saturated in the aforementioned coding regions at intergeneric cyprinodontiformes comparisons.  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. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords high-throughput sequencingaustrolebias charruacyprinodontiformescomparative mitogenomicsmolecular markersissue-copyright-statement© the author 2015 