Molecular evolution of globin genes in Gymnotiform electric fishes: relation to hypoxia tolerance
© The Author(s). 2017
Received: 5 July 2016
Accepted: 26 January 2017
Published: 13 February 2017
Nocturnally active gymnotiform weakly electric fish generate electric signals for communication and navigation, which can be energetically taxing. These fish mainly inhabit the Amazon basin, where some species prefer well-oxygenated waters and others live in oxygen-poor, stagnant habitats. The latter species show morphological, physiological, and behavioral adaptations for hypoxia-tolerance. However, there have been no studies of hypoxia tolerance on the molecular level. Globins are classic respiratory proteins. They function principally in oxygen-binding and -delivery in various tissues and organs. Here, we investigate the molecular evolution of alpha and beta hemoglobins, myoglobin, and neuroglobin in 12 gymnotiforms compared with other teleost fish.
The present study identified positively selected sites (PSS) on hemoglobin (Hb) and myoglobin (Mb) genes using different maximum likelihood (ML) methods; some PSS fall in structurally important protein regions. This evidence for the positive selection of globin genes suggests that the adaptive evolution of these genes has helped to enhance the capacity for oxygen storage and transport. Interestingly, a substitution of a Cys at a key site in the obligate air-breathing electric eel (Electrophorus electricus) is predicted to enhance oxygen storage of Mb and contribute to NO delivery during hypoxia. A parallel Cys substitution was also noted in an air-breathing African electric fish (Gymnarchus niloticus). Moreover, the expected pattern under normoxic conditions of high expression of myoglobin in heart and neuroglobin in the brain in two hypoxia-tolerant species suggests that the main effect of selection on these globin genes is on their sequence rather than their basal expression patterns.
Results indicate a clear signature of positive selection in the globin genes of most hypoxia-tolerant gymnotiform fishes, which are obligate or facultative air breathers. These findings highlight the critical role of globin genes in hypoxia tolerance evolution of Gymnotiform electric fishes.
Teleosts of the nocturnally-active neotropical clade Gymnotiformes produce and detect weak electric signals for the purposes of electrolocation and communication . Variation in the patterning or frequency of the electric organ discharge (EOD) plays a vital role in electrical communication during behaviors such as aggression, courtship, and mating . EODs are classified as wave- or pulse-type. Wave-type EODs are formed by regularly-emitted pulses where the pulse duration is approximately equal to the inter-pulse interval thereby approximating a sine wave. Pulse-type EODs are emitted irregularly and the EOD pulse duration is short . The order Gymnotiformes comprises six families: Hypopomidae, Rhamphichthyidae and Gymnotidae, which are all pulse-type, and Sternopygidae, Eigenmannidae and Apteronotidae, which are wave-type .
In most species EOs are composed of cells, called electrocytes, derived from muscle tissue. The EO is directly controlled by the nervous system, which commands the electrocytes to fire . Electrocytes are large cells capable of generating large ionic currents, especially sodium (Na+) currents. In order to restore the ionic gradient for Na+ following action potentials, electrocytes have large amounts of the “sodium pump” Na+/K+ ATPase; this pump uses one molecule of ATP for every three Na+ ions pumped back out of the cell. Thus, the generation of electricity requires energy. Physiological studies show that performance-related costs of EOD generation may be surprisingly high, from 10 up to 30% of routine oxygen consumption . Yet, surprisingly, the O2 consumption of gymnotiforms does not differ from that of other similarly sized teleosts  suggesting that gymnotiforms have adaptations for energy efficiency.
Gymnotiforms have their highest diversity throughout the ecologically varied lowland aquatic habitats of the vast Amazonian floodplains . Although most gymnotiforms inhabit well-oxygenated streams and rivers (including all wave-type and a few pulse-type species), several lineages (only pulse-type species) have additionally radiated within oxygen-poor, stagnant habitats. Furthermore, all species are in danger of being trapped in shrinking, hypoxic pools during the dry season. Thus, hypoxia poses additional physiological challenges, most importantly, a deficiency in oxygen to fuel oxidative metabolism. As a consequence, coping with hypoxia is a potentially daunting task for gymnotiforms.
Different gymnotiform species possess a variety of anatomical and physiological adaptations to cope with these energetic demands including some common to other fish living in hypoxic environments such as large gill surface area , the use of aquatic surface respiration (ASR, in which fish swim to the surface and take in water from the topmost heavily oxygenated layer of water over their gills), and various forms of air-breathing . In addition, some species have gymnotiform-specific means of conserving energy such as decreasing EOD amplitude (which lowers Na+ influx and, therefore, the amount of ATP needed to power Na+/K+ ATPase) during the day when they are inactive, upon encountering hypoxic conditions, or when food is scarce [10–12].
Little is known, however, about molecular adaptations for oxygen delivery or usage. Therefore, in this study, we focused on the well-known oxygen carriers, the globins. Globins are the most widespread respiratory proteins, existing in fungi, plants, and animals [13, 14]. Globins are conserved metalloproteins that typically have seven alpha helices that form a pocket containing an oxygen-binding heme . They have been investigated for over a century, and eight globin types have been identified in vertebrates including hemoglobin (Hb), myoglobin (Mb), and neuroglobin (Ngb). Hb and Mb are best known for their respiratory functions, which play critical roles in the maintenance of cellular oxygen supply in support of aerobic metabolism [16, 17].
Hb is a heterotetramer protein composed by two α- and two β-chains, which are encoded by the corresponding alpha (Hba) and beta (Hbb) globin gene family . Moreover, it is responsible for facilitating O2 from the respiratory system to the inner organs via the circulatory system . Mb is a compact and highly soluble monomer protein containing one proto-heme, which is involved in the oxygen storage and transportation within heart and skeletal muscle cells , and has a higher oxygen affinity than Hb. Ngb is a monomer and an oxygen-carrying protein essentially restricted to neurons , which plays a key role in facilitated diffusion and local storage of O2. Although its function is still uncertain, there is general agreement that Ngb is associated with mitochondria and thus oxidative metabolism, and serves a neuroprotective role during hypoxic stress . Interestingly, Hb, Mb and Ngb have also been proposed to have roles in nitric oxide (NO) metabolism and the detoxification of reactive oxygen species during hypoxia .
We gathered coding sequences of Hba, Hbb, Mb, and Ngb genes from 12 gymnotiform species to: 1) test whether these globin genes have evolved adaptively (i.e.,—show signs of positive selection) during gymnotiform origins and evolution; 2) evaluate whether gymnotiforms that inhabit hypoxic/anoxic vs. well-oxygenated water display different patterns of molecular evolution and; 3) provide a more comprehensive picture of the hypoxia tolerance in gymnotiforms.
Taxon sampling, and primary treatments of data
Characteristics of species used in this study
Positive selected genes
Genomic DNA and tissue RNA isolation and sequencing
Genomic DNA was extracted using QIAGEN DNeasy Tissue Kit (QIAGEN, Inc.). The total RNA was isolated from skeletal muscle, electric organ and brain using the QIAGEN RNeasy Mini Kit (QIAGEN, Inc.), and used as a template for single-stranded (ss) cDNA synthesis using the SuperScript III reverse transcriptase protocol (Invitrogen). Degenerate and specific primers were designed using the homologous sequences from the genome of electric fishes by using Oligo . The primers are listed in Additional file 1: Table S3. This ss cDNA (1:10 dilution) and genomic DNA were used as the template with degenerate primer pairs in the first nested PCR, and their products were templates in second nested PCR reactions with primer pairs. PCR amplifications were carried out using the following program: 95 °C 4 min, 35 cycles of 94 °C denaturation 1 min, 55–50 °C annealing 30 s, 72 °C extension 1 min, and 72 °C elongation for 15 min. PCR products were sequenced in both directions. The newly determined sequences were deposited in GenBank under the accessions KX827275-KX827283, KX833066-KX833070 and KX982253- KX982258.
Sequences were assembled and primer sequences were trimmed using MacVector . Nucleotide sequences were translated into amino acid sequences and aligned using MEGA 6.0 . The alignments were manually inspected and edited by eye (see Additional file 2: Figure S1). The best fit model of nucleotide evolution was determined by jModeltest , and Maximum Likelihood (ML) trees was reconstructed by PhyML 3.0 .
Molecular evolutionary analysis
To determine whether adaptive evolution might have occurred on the globin genes of electric fishes, we used the PAML package , which uses a maximum-likelihood (ML) approach to calculate nonsynonymous to synonymous rate ratios (ω = dN/dS). The ratios >1, =1 and <1 indicate positive selection, neutrality and negative selection, respectively. For all genes, both the species tree (Fig. 1) and putative ML tree (Additional file 2: Figure S2) were separately used as the working topology in all the analyses.
We used a site model for positive selection at individual codons in electric fish samples for each gene, i.e., M8 and M8a . Model M8a only allows codons to evolve neutrally or under purifying selection (ω <1), whereas M8 model includes a class of sites with ω > 1. Amino acids under selection for model M8 were identified using a Bayes empirical Bayes approach (BEB), and we considered as candidates sites with a posterior probability >80% . Then, we further employed a series of ML methods implemented in the Datamonkey web server (http://www.datamonkey.org), which has the advantage of improving the estimation of the dN/dS ratio by incorporating variation in the rate of synonymous substitution . The single likelihood ancestor counting (SLAC) model is a conservative test, which counts the synonymous and nonsynonymous changes at each codon position in a phylogeny. The Fixed-Effect Likelihood (FEL) calculates site-by-site dN/dS without assuming a prior distribution. The Random-Effect Likelihood (REL) assumes a prior distribution across sites. In addition, the Fast Unconstrained Bayesian AppRoximation (FUBAR) ensures robustness against model misspecification. Each module was run using the default cutoffs with p = 0.2 for SLAC and FEL, Bayes Factor = 50 for REL and posterior probability = 0.8 for FUBAR.
To test for possible heterogeneity of ω ratios along independent branches, we used the free-ratio model, which allows each branch to have a separate dN/dS. The null model is a very strict model called the one-ratio model (M0) that allows only a single ω ratio for all branches. We further executed branch-site tests to explore positive selection affected by a few sites along a specific branch . We compared modified model A, which assumes four classes of sites, especially, allowing codons under positive selection along foreground lineage with ω2 > 1, to the null hypothesis, in which fixed ω2 = 1 is allowed based on branch-site model A. For all the analyses, the nested models were compared using a likelihood ratio test (LRT) with various degrees of freedom, and all analyses were run twice to ensure convergence. Branch-site REL analysis was also performed to determine whether specific-lineage is evolving under positive selection using a web-based implementation of HyPhy package (http://www.datamonkey.org), which is based on likelihood ratio tests that identify all lineages with a proportion of sites that are evolving with dN/dS > 1, and do not require partitioning lineages into foreground and background branches .
Each gene sequence alignment was also analyzed using the program TreeSAAP , which further supports PAML and Datamonkey at the protein physicochemical level. TreeSAAP compares the magnitude of physicochemical changes of non-synonymous residues across a phylogeny and identifies specific amino acid properties that have likely been affected by positive destabilizing selection during evolution. In this study, amino acid properties were considered to be under positive-destabilizing selection if positive selection was detected in radical magnitude ranges (6–8). The number of radical changes in the amino acid properties was used as a proxy for determining the strength of positive selection. More radical changes in amino acid properties might suggest adaptive evolution. The residues that had fewer than six amino acid property changes were categorized as type I sites, whereas those that had more than six were categorized as type II sites.
To provide further insights into the underlying effects of these positively selected sites, we mapped them onto the three-dimensional (3D) structure. The crystal structures of P11748 for Hba, P11749 for Hbb, and 2NRL for Mb were taken from the Protein Data Bank (http://www.rcsb.org/pdb). Postulated functional regions or residues were searched from Uniprot (http://www.uniprot.org/uniprot). Pymol  was used to produce the images of the three-dimensional models of corresponding gene.
Tissues were removed from a single Brachyhypopomus gauderio that was previously housed under normoxic conditions in the laboratory and total RNA was extracted from brain, skeletal muscle, heart and electric organ. RNA was treated with RiboZero (Illumina, MRZH-11124) kit to remove ribosomal RNA, and cDNA libraries (200 bp, paired ends) were made. Libraries were sequenced on an Illumina HiSeq machine. We processed the raw reads with Trimmomatic v0.32  for adapter removal (IlluminaClip: TruSeq3-PE.fa:2:30:10), quality trimming (SlidingWindow:4:5, Leading:5, Trailing:5) and size filtering (MinLen:25). These are Trinity’s default settings, which are based on the work of MacManes . We performed quality control of both raw and processed reads with FastQC v0.11.3 (http://www.bioinformatics.babraham.ac.uk/) (Table S6).
We combined the processed PE reads across organs, in silico normalized, and de novo assembled them with Trinity v2.2.0 [37, 38]. Then we used BUSCO v1.1b1 , along with BLAST+ v2.2.31 , HMMER v3.1b1 , and EMBOSS v6.5.7  to measure transcriptome completeness, against the Vertebrates, Metazoans and Eukaryotes datasets. In all cases, the assembly displayed a large percentage of complete orthologs (Table S7).
We estimated levels of expression for each transcript and gene, per organ, using the Trinity-provided scripts align_and_estimate_abundance.pl and abundance_estimates_to_matrix.pl. We used both RSEM v1.2.19  and kallisto v0.42.5  quantification methods, which produced qualitatively very similar results. Reported abundances are TMM-normalized values calculated with the RSEM method. Only one Trinity gene blasted against the myoglobin gene sequence with an E-value of 0.00. This gene’s per organ abundances are the ones reported.
The abundance estimation results suggested that very few genes accounted for a large fraction of gene expression. Upon inspection of said genes, many were related with rRNA, and therefore were expected to be depleted during the library preparation process. Although our B. gauderio transcriptome assembly meets our quality standards, we don’t recommend future use of the RiboZero kit (which is designed for human, mouse and rat) when working with RNAseq from this taxon. Sequences for B. gauderio genes are available at: http://efishgenomics.integrativebiology.msu.edu/blast_search/.
We successfully amplified cDNAs for the Mb, Hba, Hbb and Ngb gene from eight Gymnotiforms species (see Additional file 1: Table S2). There are no insertion/deletion mutations or changes that result in stop codons in gene sequences, suggesting the presence of functional proteins in electric fishes. These were added to mRNA sequences previously obtained by NextGen sequencing from four other species for a total of 12 species.
We constructed phylogenetic trees using maximum likelihood (ML) performed by PhyML from nucleotide alignments of the four genes dataset. The Akaike Information Criteria (AIC) in jModeltest selected the GTR + G substitution model for all genes. Relationships of the gene trees largely reflected species relationships previously estimated with morphological data by Tagliacollo et al.  (see Additional file 2: Figure S2). For example, the phylogenetic tree placed Gymnotiformes and Siluriformes together, and they had a closer relationship with Characiformes than Cypriniformes in the Mb gene tree. However, the topology based on the genes still failed to resolve the relationships within Gymnotiformes, which is not surprising since this has been difficult to resolve even with large datasets. The bootstrap support values for these relationships were not high, which is most likely due to the short length of globin genes and small number of species that were used to reconstruct the topologies.
Adaptive molecular evolution
Evidence of positive selection identified by different methods in Hyphy
SLAC (P > 0.8)
FEL (P > 0.8)
REL (PP > 50)
Fubar (PP > 0.8)
35 (0.134), 73 (0.199), 109 (0.133)
35 (0.058), 73 (0.077), 79 (0.111), 83 (0.132), 109 (0.073)
35 (60.944), 73 (60.689), 83 (50.668), 109 (84.333)
35 (0.909), 73 (0.840), 79 (0.940), 109 (0.824)
E.electricus: p = 0.007
123 (0.197), 133 (0.154)
15 (0.093), 45 (0.146), 81 (0.127), 106 (0.164), 112 (0.162), 123 (0.072), 133 (0.067), 144 (0.100)
123 (60.003), 133 (1832.56)
15 (0.818), 81 (0.874), 111 (0.899), 123 (0.854), 133 (0.979), 144 (0.921)
B.gauderio: p = 0.003;
E.virescens: p = 0.037
27 (0.097), 98 (0.096)
27 (0.021), 92 (0.110), 98 (0.029), 108 (0.131)
27 (337.517), 98 (94.269)
27 (0.973), 44 (0.849), 92 (0.849), 98 (0.925)
E.electricus: p = 0.005
B.gauderio: p = 0.017
Selective pressure analyses of Gymnitifores by the Branch-Site Model
Positive selected sites
ω0 = 0.129, ω1 = 1, ω2 = 6.251
ω0 = 0.129, ω1 = 1, ω2 = 1
ω0 = 0.145, ω1 = 1, ω2 = 34.694
20-0.955*, 29–0.973*, 113–0.989*, 114–0.912, 132–0.981*, 133–0.888
ω0 = 0.139, ω1 = 1, ω2 = 1
ω0 = 0.137, ω1 = 1, ω2 = 15.025
9-0.987*, 18–0.964*, 98–0.940, 109–0.990**, 122–0.994**, 134–0.992**, 138–0.911
ω0 = 0.136, ω1 = 1, ω2 = 1
ω0 = 0.136, ω1 = 1, ω2 = 381.189
37-0.998**, 106–0.982*, 108–0.960*, 109–0.993**, 117–0.944
ω0 = 0.132, ω1 = 1, ω2 = 1
ω0 = 0.142, ω1 = 1, ω2 = 999
ω0 = 0.141, ω1 = 1, ω2 = 1
Last common ancestor of gymnotiform
ω0 = 0.143, ω1 = 1, ω2 = 60.065
ω0 = 0.141, ω1 = 1, ω2 = 1
ω0 = 0.063, ω1 = 1, ω2 = 189.298
ω0 = 0.063, ω1 = 1, ω2 = 1
Radical amino acid sites under positive selection detected by Datamonkey, Branch-site model and TreeSAAP Simultaneously
Radical Changes in AA Propertiesc
pK’, R a ,
N s , R a
K o , H t
N s , B r , pK’ , R a , H p , H t
F, E t
P α , P c , P
P α , P c , P
P α , N s , c , K o , pK’ , α n
P α , P c , P
P α , P c , P
N s , R a
N s , B r , c
P β , B l , P c , F , R a , P
N s , B r , c , h , p , E t
N s , B r , μ
P α , P c , P
N s , B r , c , h , p , E t
N s , B r , c
P α , N s , B r , c , h , E l , P r , p , α c , H p , E t
E t , B r , K o , h, E sm
B l , α c , H t
B l ,
N s , R a
B l , R F , P c , C a , E l , F , M v , M w , V 0 , u , R a , H t
N s , B r , pK’ , h , F , p , E t , E l , R a , H p , H t
N s , P β , B r , R F , pK’ , h , E l , F , H nc , P r , p , R a , H p , E t
Ph i , E t
P α , P c , P
pK’, R a
P α , P c , P
P β , B l , P c , F , R a , P
P α , P c , P
pK’, C, R a , P
pK’, C, R a , P
N s , B r , h , E l , F , p , R a , H p , H nc , E t
K o , E t
P α , P c , P
Structural links to protein function
Expression profile of Mb and Ngb in E. Electricus
Transcriptomes were made from some of the same tissues (brain, heart, skeletal muscle, EO) of B. gauderio (see Additional file 1: Table S6 and S7; Additional file 3). The overall pattern of Mb expression is similar to E. electricus. That is, expression is highest in the heart, also observed in brain, and low in muscle and EO. Ngb levels were too low to measure accurately.
Pervasive adaptive evolution of gymnotiform globin genes
Recent studies have shown strong evidence for positive selection in hypoxia tolerant species, like yak , hummingbirds , and cetaceans , when compared with their hypoxia intolerant relatives, as well as hypoxia-tolerant populations of humans such as Tibetans .
In this study, we surveyed all four primary oxygen-carrier globin genes in gymnotiforms for signs of positive selection and to assess whether hypoxia tolerance has influenced the evolution of these genes. Our analyses provide strong evidence that globin genes have been subjected to positive selection during gymnotiform evolution. First, neutral models of evolution were rejected for Hbb genes, and more than two ML methods identified specific codons with a high probability of being under selection for Hba, Hbb and Mb genes. Second, adaptive evolution was further supported by evidence of radical changes in positively selected amino acids in gymnotiform globins. Again, 34.6% (9/26) belong to the Type II class, suggesting robust positive selection. Finally, several of the putatively selected sites fall in, or close to, regions important for function based on structural information.
Positively selected sites are scattered pervasively along lineages of gymnotiform phylogeny (see Table 2 and Fig. 1), suggesting the contribution of the respiratory proteins for oxygen storage and transportation during adaptation to expensive oxygen consumption in gymnotiforms. Moreover, a signal of positive selection was also detected in the lineage leading to the common ancestor of gymnotiforms. This lineage represents the early evolutionary history of the gymnotiform’s evolved EOs, during which the gymnotiforms were faced with the challenges of high energetic cost for electric organ discharge generation. Although this branch was only detected in the Mb gene with three positively selected codons (Mb: N98A, G12P; Hbb: Q133C), three or more radical property changes occurred at each amino acid (see Table 4). That is to say, globin genes may have adaptively enhanced oxygen binding and transportation in accordance with the changes of high energetic cost during the early evolutionary phase of EOs in gymnotiformes.
Functional consequences of amino acid replacement
Although gymnotiform globin genes contain putatively positively selected sites, it is important to assess their functional relevance. We thus analyzed selected residues for their structural properties to predict potential functional implications. We found that all 26 radical amino acid changes residues were concordant between three methods and thus constitute robust candidates for positive selection (see Table 4). Radical substitutions of amino acids at key positions may change the properties of the molecule. For example, residue 83 in Hba is located very close to the promixal histidine, which is implicated in the iron-proximal histidine linkage . Consequently, substitution at this site seems to be essential in the maintenance of heme oxygen binding.
Furthermore, position 144 in Hbb is close to 147 His. It has been reported that His-HC3(147)β carp (Cyprinus carpio) hemoglobin plays a key role in the Root effect, which is a phenomenon associated with non-cooperative oxygen binding and decreased oxygen affinity . Hence, we suggest that positive selection acting on this site is likely to be involved in modulating hemoglobin oxygen combination and cooperation. For Mb, site 134 is mainly responsible for formation of hydrogen bonds with water in a hydrophobic environment ; therefore, amino acid changes at this site are likely to be involved in the regulation of water bonds. In spite of the evidence for selection documented here, functional studies of these candidate positively selected sites are necessary in gymnotiforms in the future.
Myoglobin and NO production
Interestingly, the Mb of Gymnarchus niloticus a member of the other independently evolved group of electric fishes, the mormyroidea, shares Cys 106 (see Fig. 5, blue) with E. electricus at a site that is conserved among other teleosts. Whereas all other (200+ species) of mormyroid fishes utilize gills and cannot breathe air under hypoxic conditions [56–58], G. niloticus is an air-breathing fish with a highly vascularized gas bladder. These fish also grow to be large (>1.5 meters) similar to E. electricus. Therefore, it appears that the shared Cys 106 is a functional convergence in two large, air-breathing electric fish.
Relationship of globin evolution to gymnotiform life histories
The two species with strongest signature of positive selection of globins—E. electricus and B. gauderio--live comfortably in hypoxic and even anoxic environments. E. electricus, which is the largest gymnotiform capable of reaching ~2 m in length, is unusual among gymnotiforms as it is an obligate air-breather. It surfaces every few minutes to gulp air, which it stores in its mouth and, if prevented from breathing air, it will drown . E. electricus is unique even among air breathing fishes as it obtains oxygen from elaborated papillae in the mouth not related to the gills; indeed, its gills are small and underdeveloped. Furthermore, its circulatory system is unlike that of other teleosts in that the oxygenated blood from the oral papillae mixes with the venous circulation before being pumped out of the heart resulting in poorly oxygenated blood. The hematocrit, Hb content, and oxygen capacity of the blood, as well as heart rate and volume of blood moved per unit time are all higher than in most teleosts and may be adaptations to overcome the poor oxygenation of mixed arterial-venous blood .
Fish of the genus Brachyhypopomus inhabit hypoxic/anoxic waters  and in the laboratory, tolerate >6 h of anoxia. They have large well-developed gills, but are also facultative air-breathers either gulping bubbles of air at the surface then descending, or “skulking” at the surface with open mouths taking in air . These two behaviors differ from ASR in that they involve taking in air, whereas ASR is merely taking in well-oxygenated water from the surface layer. Fish that hold air in their mouths must have gills that do not collapse or they would provide too little useful surface area for oxygen absorption.
Gymnotus species are also capable of tolerating >6 h of anoxia by gulping air and storing it in a highly vascularized, extended, lung-like posterior swim bladder unique among gymnotiforms . They use this swim bladder “lung” to extend their aerobic range . The swim bladder is contacted by branches of the celiac artery and hepatic portal systems as well as segmental arteries and veins that then penetrate the body musculature. Although we tried multiple primers to amplify Gymnotus omarorum Mb from muscle, heart and EO tissue, we failed to amplify this gene from RNA samples even though we obtained the Hba and Hbb genes from the same samples. Recently, Macqueen et al.  revealed that a cardiac Mb deficit has evolved repeatedly in teleosts under diverse ecological settings; in sticklebacks (Gasterosteus aculeuatus) Mb is even a pseudogene. Although we have not examined Gymnotus DNA for a possible psuedogene, we suspect that the Mb gene was lost in G. omarorum. In keeping with this, we did not note strong selection on Hbs in this genus. We suggest that their blood is adequately oxygenated in their unique vascularized swim bladder with no further selection pressure to alter Hb affinity for oxygen.
Crampton  included E. virescens among the species that are capable of inhabiting hypoxic environments. However, this species is just barely able to survive in hypoxia. They tend to avoid hypoxic water although they are capable of surviving an hour or so of anoxic conditions using ASR , but they do not gulp air. They are sensitive to hypoxia under which the amplitude of their EODs decreases rapidly . In accordance with this, we observe only minimal evidence of positive selection on Hbb in this species.
The other species in this study, which show little or only weak evidence of positive selection on globin genes, do poorly in hypoxic conditions and cannot survive anoxia . Under anoxic conditions in the laboratory they respond with ASR as long as they are able and eventually fall immobile to the bottom. These include all the wave species and some of the pulse-type fish.
Finally, the genera Brachyhypopomus and Gymnotus are speciose  with member species distributed throughout the gymnotiform’s whole range. For example, fish from these two genera are the primary species in the most Southern extent of the gymnotiform’s distribution in Uruguay [63, 64]. The evolution of different modes of air-breathing may have given them an advantage over those species unable to derive oxygen from the air.
Expression of Mb and Ngb in E. electricus and B. gauderio
As in some other teleosts, Mb is most highly expressed in heart (although this varies considerably across species ) and less so in muscle. However, as Mb is expressed in “slow” oxidative but not “fast” non-oxidative muscle fibers, and slow fibers are small and few in number in teleost muscle, its mRNA will not be abundant in a sample of epaxial muscle even if it is highly expressed in slow muscle [65, 66]. Thus, the difference in muscle Mb expression between E. electricus and B. gauderio may depend on composition of the muscle sub-types and the exact location from which samples were obtained.
Our expectation was that Mb levels would be high in the EO consistent with the ongoing EOD activity. But Mb levels are modest in the EO of both species. We did not note a difference in expression between the Sach’s organ and Hunter’s/Main EO of E. electricus. Sach’s EO constantly fires a low voltage-pulse at a low frequency (~1–10 Hz), whereas all three EOs discharge in high frequency bursts of >400 Hz when fish are capturing prey or defending themselves . B. gauderio discharges at rest at ~15 Hz and during social interactions or foraging at 30–40 Hz, with occasional high-frequency bursts (100–200 Hz) called “chirps” that are aggressive signals .
There are a few possible and not mutually exclusive explanations for the low levels of Mb in the EO. First, EO originates from fast muscle fibers, and fast muscle fibers do not express Mb; there might be a developmental constraint minimizing Mb expression in the EO. Second, some gymnotiforms possess mechanisms to decrease EOD amplitude (and therefore oxygen consumption) during anoxia . This has only been studied in a few species and it would be intriguing to compare the distribution of EO Mb expression, EOD “usage” (discharge rate), and occurrence of mechanisms for reducing EOD amplitude in the face of anoxia in the EO across gymnotiforms. Perhaps different species trade off EO Mb concentration with extent of anoxia-dependent reduction of the EOD amplitude. Third, EOs are well vascularized and it may not be necessary to have a local reserve source of O2. Fourth, it is possible that EOs switch to anaerobic metabolism during and after periods of high activity. At this point there is no evidence for or against this. It is also worth noting that Mb expression is high in heart and low in muscle in the elephant shark suggesting that this is the ancestral vertebrate pattern .
Recently, it has been recognized that Mb is expressed in the brain and kidney as well as in a few other non-muscle tissues in teleosts. Mainly, Mb is expressed in the endothelial cells of blood vessels in the brain and gills, and the epithelial cells of the tubules in the kidney . The expression of Mb in the brain in E. electricus and B. gauderio and the kidney of E. electricus is in accord with these observations.
Ngb is primarily expressed in neural tissues in mammals and teleosts. It is reported to be in the gills in zebrafish . As there are no gill transcriptomes available of E. electricus we do not know if it is expressed in the gills of this species. Ngb was expressed in E. electricus kidney but zebrafish kidney has not yet been examined. Finally, it is not expressed in zebrafish muscle nor did we observe it in E. electricus muscle. Qualitatively, the expression pattern of Mb and Ngb does not differ from non-gymnotiforms. Thus, it appears that the main effect of selection on E. electricus and B. gauderio Mb is on its sequence rather than its expression pattern in individuals housed under normoxic conditions. Further experiments would be necessary to determine whether there are additional adaptations on the regulation of globin genes under different environmental conditions.
The gymnotiforms arose in Gondwana  and diverged rapidly into lineages with differences in EO morphology and signaling. Some amino acid substitutions occurred in myoglobin genes as the EO was initially evolving, perhaps as a result of the new energy demands imposed by electrogenesis. While many gymnotiforms rely on gills for oxygen uptake, a few lineages evolved distinct modes of air-breathing. Along with this, we note variation in the strength of selection on different globin genes depending on the mode of air breathing. We anticipate that future studies will elucidate the physiology and biochemistry of these globin genes and seek other molecules under selection in the oxygen transport or metabolic pathways underlying energy usage.
The authors acknowledge Drs. Christopher Braun (Hunter College, NYC), G. Troy Smith (Indiana University), Ana Silva (Universidad de la República de Uruguay), Michael Markham (University of Oklahoma, Norman), and Hernán Lopez-Fernandez (Royal Ontario Museum) for contributing RNA sequences or tissues for PCR for this study. This work was supported by the China Scholarship Council (RT), the University of Texas (HZ), NSF grant IOS 1557857 (HZ), Nanjing Normal University (GY), NSFC grants 31325025 and 31630071 (GY), Nanjing Normal University (RT) and Priority Academic Program Development of Jiangsu Higher Education Institutions (RT).
Availability of data and materials
HZ, RT, and GY conceived and designed the study. RT and ML collected the data. YL prepared the samples and submitted the sequences. RT performed the laboratory work. RT analyzed the data. All authors contributed to writing the manuscript and approved the final version of it.
The authors declare that they have no competing interests.
Ethics approval and consent to participate
As the research did not include living or anesthetized animals, no ethical approval was necessary. All methods used to collect observational data were non-invasive, were in compliance with the regulations and guidelines of the University of Texas at Austins’ Institutional Animal Care and Use Committee.
Open AccessThis 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.
- Von der Emde G. Electroreception. In: Bullock TH, Fay CDH, Richard R, editors. Neurosciences-from molecule to behavior: a university textbook. Springer: Springer Science & Business Media; 2013. p. 409–25.Google Scholar
- Caputi AA, Carlson BA, Macadar O. Electric organs and their control. In: Electroreception. New York: Springer; 2005. p. 410–51.View ArticleGoogle Scholar
- Hopkins CD. Neuroethology of electric communication. Ann Rev Neurosci. 1988;11:497–535.View ArticlePubMedGoogle Scholar
- Tagliacollo VA, Bernt MJ, Craig JM, Oliveira C, Albert JS. Model-based total evidence phylogeny of Neotropical electric knifefishes (Teleostei, Gymnotiformes). Mol Phylogenet Evol. 2016;95:20–33.View ArticlePubMedGoogle Scholar
- Salazar VL, Krahe R, Lewis JE. The energetics of electric organ discharge generation in gymnotiform weakly electric fish. J Exp Biol. 2016;216:2459–68.View ArticleGoogle Scholar
- Julian D, Crampton WGR, Wohlgemuth SE, Albert JS. Oxygen consumption in weakly electric Neotropical fishes. Oecologia. 2003;137:502–11.View ArticlePubMedGoogle Scholar
- Crampton WGR. Gymnotiform fish: an important component of Amazonian fioodplain fish communities. J Fish Biol. 1996;48:298–301.Google Scholar
- Crampton WGR, Chapman LJ, Bell J. Interspecific variation in gill size is correlated to ambient dissolved oxygen in the Amazonian electric fish Brachyhypopomus (Gymnotiformes: Hypopomidae). Environ Biol Fish. 2008;83:223–35.View ArticleGoogle Scholar
- Val AL, Silva MNP, Almeida-Val VMF. Hypoxia adaptation in fish of the Amazon: a never-ending task. S Afr J Zool. 1998;33:107–14.View ArticleGoogle Scholar
- Crampton WGR. Effects of anoxia on the distribution, respiratory strategies and electric signal diversity of gymnotiform fishes. J Fish Biol. 1998;53:307–30.View ArticleGoogle Scholar
- Reardon EE, Parisi A, Krahe R, Chapman LJ. Energetic constraints on electric signalling in wave-type weakly electric fishes. J Exp Biol. 2011;214:4141–50.View ArticlePubMedGoogle Scholar
- Sinnett PM, Markham MR. Food deprivation reduces and leptin increases the amplitude of an active sensory and communication signal in a weakly electric fish. Horm Behav. 2015;71:31–40.View ArticlePubMedGoogle Scholar
- Burmester T, Hankeln T. Function and evolution of vertebrate globins. Acta Physiol. 2014;211:501–14.View ArticleGoogle Scholar
- Hardison RC. A brief history of hemoglobins: plant, animal, protist, and bacteria. Proc Natl Acad Sci U S A. 1996;93:5675–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Bashford D, Chothia C, Lesk AM. Determinants of a protein fold: Unique features of the globin amino acid sequences. J Mol Biol. 1987;196:199–216.View ArticlePubMedGoogle Scholar
- Rahbar S. Hemoglobin: Structure, function, evolution, and pathology. Am J Hum Genet. 1983;35:781.PubMed CentralGoogle Scholar
- Wittenberg BA, Wittenberg JB. Transport of oxygen in muscle. Annu Rev Physiol. 1989;51:857–78.View ArticlePubMedGoogle Scholar
- Dickerson RE, Geis I. Hemoglobin: structure, function, evolution and pathology. Menlo Park: Benjamin-Cummings Publishing Company; 1983.Google Scholar
- Gros G, Wittenberg BA, Jue T. Myoglobin’s old and new clothes: from molecular structure to function in living cells. J Exp Biol. 2010;213:2713–25.View ArticlePubMedPubMed CentralGoogle Scholar
- Hankeln T, Wystub S, Laufs T, Schmidt M, Gerlach F, Saaler-Reinhardt S, Reuss S, Burmester T. The cellular and subcellular localization of neuroglobin and cytoglobin clue to their function? IUBMB life. 2004;56:671–9.View ArticlePubMedGoogle Scholar
- Vinogradov SN, Moens L. Diversity of globin function: enzymatic, transport, storage, and sensing. J Biol Chem. 2008;283:8773–7.View ArticlePubMedGoogle Scholar
- Offerman JD, Rychlik W. Oligo primer analysis software. In: Krawetz SA, Womble DD (Eds.). Introduction to Bioinformatics. A Theoretical and Practical Approach. New Jersey: Humana Press Inc.; 2003. p. 345–61.Google Scholar
- Griffin AM, Griffin HG, Olson SA. MacVector: an integrated sequence analysis program for the Macintosh. Computer Analysis of Sequence Data: Part II. 1994. p. 195–201.Google Scholar
- Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: molecular evolutionary genetics analysis version 6.0. Mol Biol Evol. 2013;30(12):2725–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Posada D. jModelTest: phylogenetic model averaging. Mol Biol Evol. 2008;25:1253–6.View ArticlePubMedGoogle Scholar
- Guindon S, Dufayard JF, 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.View ArticlePubMedGoogle Scholar
- Yang ZH. PAML 4: phylogenetic analysis by maximum likelihood. Mol Biol Evol. 2007;24:1586–91.View ArticlePubMedGoogle Scholar
- Swanson WJ, Nielsen R, Yang Q. Pervasive adaptive evolution in mammalian fertilization proteins. Mol Biol Evol. 2003;20:18–20.View ArticlePubMedGoogle Scholar
- Yang Z, Wong WSW, Nielsen R. Bayes empirical bayes inference of amino acid sites under positive selection. Mol Biol Evol. 2005;22:1107–18.View ArticlePubMedGoogle Scholar
- Pond SLK, Frost SDW. Datamonkey: rapid detection of selective pressure on individual sites of codon alignments. Bioinformatics. 2005;21:2531–3.View ArticlePubMedGoogle Scholar
- Zhang J, Nielsen R, Yang Z. Evaluation of an improved branch-site likelihood method for detecting positive selection at the molecular level. Mol Biol Evol. 2005;22:2472–9.View ArticlePubMedGoogle Scholar
- 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.View ArticleGoogle Scholar
- Woolley S, Johnson J, Smith MJ, Crandall KA, McClellan DA. TreeSAAP: selection on amino acid properties using phylogenetic trees. Bioinformatics. 2003;19:671–2.View ArticlePubMedGoogle Scholar
- DeLano WL. The PyMOL molecular graphics system. San Carlos: DeLano Scientific; 2002.Google Scholar
- Bolger AM, Marc L, Bjoern U. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20.View ArticlePubMedPubMed CentralGoogle Scholar
- MacManes MD. On the optimal trimming of high-throughput mRNA sequence data. Front Genet. 2014;5:13.View ArticlePubMedPubMed CentralGoogle Scholar
- Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Trinity: reconstructing a full-length transcriptome without a genome from RNA-Seq data. Nat Biotechnol. 2011;29:644.View ArticlePubMedPubMed CentralGoogle Scholar
- Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat Protoc. 2013;8:1494–512.View ArticlePubMedGoogle Scholar
- Simão FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics. 2015;31(19):3210–2.View ArticlePubMedGoogle Scholar
- Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, Madden TL. BLAST+: architecture and applications. BMC bioinformatics. 2009;10:1.View ArticleGoogle Scholar
- Finn RD, Clements J, Eddy SR. HMMER web server: interactive sequence similarity searching. Nucleic Acids Res. 2011;39(Web Server issue):W29–37.View ArticlePubMedPubMed CentralGoogle Scholar
- Rice P, Ian L, Alan B. EMBOSS: the European molecular biology open software suite. Trends Genet. 2000;16:276–7.View ArticlePubMedGoogle Scholar
- Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC bioinformatics. 2011;12:1.Google Scholar
- Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol. 2016;34:525–7.View ArticlePubMedGoogle Scholar
- Gallant JR, Traeger LL, Volkening JD, Moffett H, Chen PH, et al. Genomic basis for the convergent evolution of electric organs. Science. 2014;344:1522–5.View ArticlePubMedGoogle Scholar
- Traeger LL, Volkening JD, Moffett H, Gallant JR, Chen PH, et al. Unique patterns of transcript and miRNA expression in the South American strong voltage electric eel (Electrophorus electricus). BMC Genomics. 2015;16:243.View ArticlePubMedPubMed CentralGoogle Scholar
- Qiu Q, Zhang G, Ma T, Qian W, Wang J, et al. The yak genome and adaptation to life at high altitude. Nat Genet. 2012;44:946–9.View ArticlePubMedGoogle Scholar
- Projecto-Garcia J, Natarajan C, Moriyama H, Weber RE, Fago A, et al. Repeated elevational transitions in hemoglobin function during the evolution of Andean hummingbirds. Proc Natl Acad Sci U S A. 2013;110:20669–74.View ArticlePubMedPubMed CentralGoogle Scholar
- Tian R, Wang Z, Niu X, Zhou K, Xu S, Yang G. Evolutionary genetics of hypoxia tolerance in cetaceans during diving. Genome Biol Evol. 2016;8:827–39.View ArticlePubMedPubMed CentralGoogle Scholar
- Yi X, Liang Y, Huerta-Sanchez E, Jin X, Cuo ZX, et al. Sequencing of 50 human exomes reveals adaptation to high altitude. Science. 2010;329:75–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Parkhurst LJ, Goss DJ, Perutz MF. Kinetic and equilibrium studies on the role of the. beta.-147 histidine in the Root effect and cooperativity in carp hemoglobin. Biochemistry. 1983;22:5401–9.View ArticleGoogle Scholar
- Friedman JM, Scott TW, Stepnoski RA, et al. The iron-proximal histidine linkage and protein control of oxygen binding in hemoglobin. A transient Raman study. J Biol Chem. 1983;258:10564–72.PubMedGoogle Scholar
- Teeter MM. Myoglobin cavities provide interior ligand pathway. Protein Sci. 2004;13:313–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Helbo S, Fago A. Allosteric modulation by S-nitrosation in the low-O2 affinity myoglobin from rainbow trout. Am J Physiol-Reg I. 2011;300:R101–8.Google Scholar
- Helbo S, Gow AJ, Jamil A, Howes BD, Smulevich G, Fago A. Oxygen-linked S-nitrosation in fish myoglobins: A cysteine-specific tertiary allosteric effect. PLoS One. 2014;9, e97012.View ArticlePubMedPubMed CentralGoogle Scholar
- Chapman LJ, Chapman CA. Hypoxia tolerance of the mormyrid Petrocephalus catostoma: implications for persistence in swamp refugia. Copeia. 1998;998:762–8.View ArticleGoogle Scholar
- Chapman LJ, Hulen KG. Implications of hypoxia for the brain size and gill morphometry of mormyrid fishes. J Zool. 2001;254:461–72.View ArticleGoogle Scholar
- Berenbrink M, Koldkjaer P, Kepp O, Cossins AR. Evolution of oxygen secretion in fishes and the emergence of a complex physiological system. Science. 2005;307:1752–7.View ArticlePubMedGoogle Scholar
- Johansen K, Lenfant C, Scmidt-Nielsen K, Petersen JA. Gas exchange and the control of breathing in the electric eel, Electrophorus electricus. Z Vergl Physiologie. 1968;61:137–63.View ArticleGoogle Scholar
- Liem KF, Eclancher B, Fink WL. Aerial Respiration in the Banded Knife Fish Gymnotus Carapo (teleostei: Gymnotoidei). Physiol Zool. 1984;57:185–95.View ArticleGoogle Scholar
- McKenzie DJ, Steffensen JF, Taylor EW, Abe AS. The contribution of air breathing to aerobic scope and exercise performance in the banded knifefish Gymnotus carapo L. J Exp Biol. 2012;215:1323–30.View ArticlePubMedGoogle Scholar
- Macqueen DJ, Johnston IA. Cardiac myoglobin deficit has evolved repeatedly in teleost fishes. Biol Lett. 2014;10:20140225.View ArticlePubMedPubMed CentralGoogle Scholar
- Richer-de-Forges MM, Crampton WGR, Albert JS. A new species of Gymnotus (Gymnotiformes, Gymnotidae) from Uruguay: description of a model species in neurophysiological research. Copeia. 2009;3:538–44.View ArticleGoogle Scholar
- Giora J, Malabarba LR. Brachyhypopomus gauderio, new species, a new example of underestimated species diversity of electric fishes in the southern South America (Gymnotiforme, Hypopomidae). Zootaxa. 2009;2093:60–8.Google Scholar
- Cossins AR, Williams DR, Foulkes NS, Berenbrink M, Kipar A. Diverse cell-specific expression of myoglobin isoforms in brain, kidney, gill and liver of the hypoxia-tolerant carp and zebrafish. J Fish Biol. 2009;212:627–38.Google Scholar
- Jaspers RT, Testerink J, Della Gaspera B, Chanoine C, Bagowski CP, van der Laarse WJ. Increased oxidative metabolism and myoglobin expression in zebrafish muscle during chronic hypoxia. Biol Open. 2014;3(8):718–27.View ArticlePubMedPubMed CentralGoogle Scholar
- Catania K. The shocking predatory strike of the electric eel. Science. 2014;346:1231–4.View ArticlePubMedGoogle Scholar
- Silva AC, Perrone R, Zubizarreta L, Batista G, Stoddard PK. Neuromodulation of the agonistic behavior in two species of weakly electric fish that display different types of aggression. J Exp Biol. 2013;216:2412–20.View ArticlePubMedGoogle Scholar
- Opazo JC, Lee AP, Hoffmann FG, Toloza-Villalobos J, Burmester T, Venkatesh B, Storz JF. Ancient duplications and expression divergence in the globin gene superfamily of vertebrates: insights from the elephant shark genome and transcriptome. Mol Biol Evol. 2015;32:1684–94.View ArticlePubMedPubMed CentralGoogle Scholar
- Fuchs C, Heib V, Kiger L, Haberkamp M, Roesner A, Schmidt M, Hamdane D, Marden MC, Hankeln T, Burmester T. Zebrafish reveals different and conserved features of vertebrate neuroglobin gene structure, expression pattern, and ligand binding. J Biol Chem. 2004;279:24116–22.View ArticlePubMedGoogle Scholar
- Lavoué S, Miya M, Arnegard ME, Sullivan JP, Hopkins CD, Nishida M. Comparable ages for the independent origins of electrogenesis in African and South American weakly electric fishes. PLoS One. 2012;7, e36287.View ArticlePubMedPubMed CentralGoogle Scholar
- Silva A, Quintna L, Galeano M, Errandonea P. Biogeography and breeding in Gymnotiformes from Uruguay. Environ Biol Fish. 2003;66:329–38.View ArticleGoogle Scholar