Discovery and evolution of novel hemerythrin genes in annelid worms
© The Author(s). 2017
Received: 21 November 2016
Accepted: 10 March 2017
Published: 23 March 2017
Despite extensive study on hemoglobins and hemocyanins, little is known about hemerythrin (Hr) evolutionary history. Four subgroups of Hrs have been documented, including: circulating Hr (cHr), myohemerythrin (myoHr), ovohemerythrin (ovoHr), and neurohemerythrin (nHr). Annelids have the greatest diversity of oxygen carrying proteins among animals and are the only phylum in which all Hr subgroups have been documented. To examine Hr diversity in annelids and to further understand evolution of Hrs, we employed approaches to survey annelid transcriptomes in silico.
Sequences of 214 putative Hr genes were identified from 44 annelid species in 40 different families and Bayesian inference revealed two major clades with strong statistical support. Notably, the topology of the Hr gene tree did not mirror the phylogeny of Annelida as presently understood, and we found evidence of extensive Hr gene duplication and loss in annelids. Gene tree topology supported monophyly of cHrs and a myoHr clade that included nHrs sequences, indicating these designations are functional rather than evolutionary.
The presence of several cHrs in early branching taxa suggests that a variety of Hrs were present in the common ancestor of extant annelids. Although our analysis was limited to expressed-coding regions, our findings demonstrate a greater diversity of Hrs among annelids than previously reported.
KeywordsBlood pigments Respiratory proteins Transcriptome Annelida
Metabolism in metazoans requires oxidation of organic molecules. Thus natural selection has presumably favored proteins that can reversibly bind and transport oxygen . Such oxygen-binding proteins likely originated from enzymes whose primary function was to protect the body from oxygen toxicity, and, secondarily, these enzymes acquired the ability to carry oxygen molecules . Several different classes of oxygen-carrying proteins, or respiratory pigments, are found across animal life. Although these molecules can reversibly bind oxygen, their binding affinities and evolutionary origins differ. In animals, oxygen-binding proteins are usually divided into two main groups: proteins that use iron to bind oxygen, including hemoglobins and hemerythrins, and hemocyanins that use copper . Although hemoglobins and hemocyanins have been extensively investigated [4–8], knowledge on the evolutionary history of hemerythrins is limited . Interestingly, medical sciences have increasingly been taking advantage of oxygen-binding proteins as blood substitutes [10, 11] or as carrier proteins for synthetic vaccines, (e.g., cancer vaccines; [12, 13]) making further study of oxygen binding protein diversity and evolution appealing.
Hemerythrins (Hrs) are a non-heme oligomeric protein family within the ‘four-helical up-and-down bundle’ fold and ‘all alpha proteins’ class according to the Structural Classification of Proteins database (SCOP) . Oxygen-binding Hr proteins contain approximately 120 amino acid residues in a single domain and transport oxygen with the aid of two iron ions that bind directly to the polypeptide chain. Residues involved in iron binding include histidines (His) in positions 26 56, 75, 79, and 108, glutamic acid residue (Glu) in position 60 and aspartic acid residue (Asp) in position 113 (position numbers from Themiste zostericola in ). Presence of these signature residues indicates putative respiratory function for Hrs . Functional Hr subunits usually form a homooctamer, although dimeric, trimeric, or tetrameric Hrs have been observed in some sipunculid species, including Phascolosoma arcuatum, P. agasizii, and Siphonosoma funafuti [17–19]. The crystal structure of Hrs consists of a bundle of four antiparallel α-helices (A, B, C, and D) formed by polypeptides: an A α-helix formed by 19 amino acid residues from position 19 to 38, B α-helix with 23 amino acids residues from position 43 to 65, C α-helix formed by 16 amino acids residues from position 72 to 88, and D α-helix formed by 20 amino acids residues from position 98 to 118, using T. zostericola as the reference sequence . The core of active sites contains two iron atoms bridged by two carboxylate groups from aspartate and glutamate residues and an oxygen-containing ligand [20, 21]. Binding of oxygen apparently requires other currently unknown cellular factors since purified Hr, by itself, usually does not bind oxygen [22, 23]. Observed oxygen binding capacity is about 25% greater in Hrs than heme-based proteins, including hemoglobins .
Although Hr-like proteins have also been reported in prokaryotes [24–26] oxygen binding Hrs have only been reported from marine invertebrates belonging to Annelida (which include sipunculids; ) Brachiopoda, Priapulida, Bryozoa, and a single species of both Cnidaria (Nematostella vectensis) and Arthropoda (Calanus finmarchicus) [9, 18, 28, 29]. Given this phylogenetic breadth of animals, whether all metazoan Hrs share a common origin is debated [28, 29]. Overall, Hr proteins exhibit variation in their quaternary structure, and four groups have been reported based mainly on their primary structure and location within animal bodies . Specifically, hemerythrins found in vascular tissue, referred to here as circulating hemerythrins (cHrs), and muscle-specific myohemerythrins (myoHr), have been better characterized compared to the other two hemerythrin groups, ovohemerythrins (ovoHr) and neurohemerythrins (nHr) [28, 30]. cHrs are polymeric intracellular proteins that occur inside nucleated cells, hemerythrocytes or pink blood cells located in coelomic fluid or vascular systems of Hr-bearing organisms . In contrast, myoHrs are monomeric cytoplasmic proteins present in muscle cells of annelids . The main difference between these groups is the presence of a five-codon insertion found in myoHr immediately before the D α-helix. Expression of myoHrs seems not to be restricted to cHrs-bearing organisms, considering that some annelid species possess both myoHrs and hemoglobins . The other two groups of Hrs, ovohemerythrin (ovoHr) and neurohemerythrin (nHr), are also intracellular and non-circulating. ovoHr was identified in oocytes of the leech Theromyzon tessulatum and its presence during oogenesis possibly suggests a complex function in iron storage and detoxification [32, 33]. On the other hand, nHr was recently discovered in neural and non-neuronal tissues from the body wall of the leech Hirudo medicinalis, and it exhibits upregulation in response to septic injury . Nevertheless, Vanin et al. suggested that nHr of leech may in fact be a myoHr . Such diversification in Hr function may have involved gene duplications resulting in new proteins via neo- or subfunctionalization . Moreover, ovoHrs and nHrs have only been reported in the literature a few times, and more studies are required to understand their function and evolution.
Annelids have the greatest diversity of oxygen-binding proteins among metazoans  and it is the only phylum from which all subtypes of Hr proteins have been documented . While Hrs of annelids have been studied since the middle of the 20th century, until the 1990s, Hrs were recorded only from sipunculids and from a single polychaete family, Magelonidae [18, 36]. Later, Vanin et al.  found Hrs in a nereid and a leech and Bailly et al.  discovered Hrs genes in seven annelid species, suggesting Hrs are broadly distributed in annelids. Given the diversity of lifestyles among annelids known to have Hrs [29, 37], and the lack of information about Hrs in general , the occurrence and diversity of these molecules may be higher than currently recognized. Thus, to examine a wide diversity of annelid taxa for Hrs and to further understand how different forms of Hrs are evolutionarily related to each other, we employed approaches to survey Hrs from a diverse array of annelid transcriptomes in silico. We identified Hrs in 44 taxa and further describe the molecular diversity and evolution of Hrs in the light of annelid phylogeny [27, 38, 39]. Along with this, we assess whether described Hr subtypes consist of evolutionary lineages or result of independent adaptations to different organismal tissues.
List of taxa, including collection site, total number of base pairs sequenced, total number of contigs after assembly, number of putative Hr genes, and GenBank accession numbers. GenBank accession numbers for each Hr copy is indicated in Additional file 2
Total contigs number
Hr genes number
N 33° 07.’ W 076° 06.4’
Paralvinella palmiformis Desbruyères & Laubier, 1986
N 47° 56.9’ W 129° 05.9’
Paramphinome jeffreysii (McIntosh, 1868)
N 63° 30.8’ E 10° 25.0’
KY007424 to KY007431
Aphrodita japonica Marenzeller, 1879
N 48° 28.6’ W 122° 58.7’
KY007279 to KY007287
Aspidosiphon laevis de Quatrefages, 1865
N 09° 22.6’ W 82° 18.1’
KY007297 to KY007300
Heteromastus filiformis (Claparede, 1864)
N 41° 41.5’ W 070° 37.6’
KY007364 to KY007371
Chaetopterus variopedatus (Renier, 1804)
N 41° 41.5’ W 070° 43.5’
KY007307 and KY007309
Mesochaetopterus taylori Potts, 1914
N 48° 29.0’ W 123° 04.3’
KY007383 to KY007386
Arichlidon gathofi Watson Russell, 2000
N 33° 59.6’ W 76° 42.1’
KY007288 to KY007296
N 66° 33.2’ E 33° 06.7’
KY007310 to KY007313
Cossura longocirrata Webster & Benedict, 1887
N 33° 29.4’ W 074° 48.0’
Enchytraeus albidus Henle, 1837
N 59o 5.2' E 16o 03.3'
Eunice pennata (Müller, 1776)
N 39° 47.2’ W 70° 46.3’
KY007332 to KY007336
Euphrosine capensis Kinberg, 1857
S 34° 09.9’ E 18° 26.0’
KY007337 to KY007339
Poeobius meseres Heath, 1930
N 36° 41.2’ W 122° 02.0’
KY007441 and KY007442
Glycera dibranchiata Ehlers, 1868
N 41° 54.1’ W 070° 00.4’
KY007350 to KY007352
Thysanocardia nigra (Ikeda, 1904)
N 48° 28.6’ W 122° 58.7’
KY007480 to KY007485
Glycinde armigera Moore, 1911
N 36° 23.0’ W 121° 57.9’
KY007353 to KY007356
Delaya leruthi (Hrabe, 1958)
N 43° 0.8’ E 01° 2.5’
KY007320 to KY007326
Oxydromus pugettensis (Johnson, 1901)
N 48° 34.3’ W 123° 10.1’
KY007421 and KY007422
N 35° 29.4’ W 074° 48.0’
KY007409 to KY007413
Magelona berkeleyi Jones, 1971
N 36° 22.8’ W 121° 58.1’
KY007372 to KY007382
Clymenella torquata (Leidy, 1855)
N 41° 42.7’ W 070° 19.7’
KY007314 to KY007314
Nicomache venticola Blake & Hilbig, 1990
N 47° 57.0’ W 129° 05.9’
KY007405 to KY007408
Nephtys incisa Malmgren, 1865
N 40° 53.0’ W 070° 25.0’
KY007396 to KY007401
Alitta succinea (Leuckart, 1847)
N 41° 54.0’ W 070° 00.3’
KY007402 to KY007404
N 39° 54.1’ W 070° 35.1’
KY007327 to KY007330
Oenone fulgida (Savigny in Lamarck, 1818)
S 34° 37.0’ E 19° 21.4’
KY007414 to KY007416
Naineris laevigata (Grube, 1855)
S 34° 35.0’ E 19° 20.9’
KY007392 to KY007395
Galathowenia oculata (Zachs, 1923)
N 66° 33.2’ E 33° 6.7’
KY007340 to KY007349
Pectinaria gouldii (Verrill, 1874)
N 41° 37.9’ W 070° 53.3’
KY007432 to KY007440
Phascolosoma agassizii Keferstein, 1866
N 48° 31.2’ W 123° 01.0’
KY007443 to KY007453
Ancistrosyllis groenlandica McIntosh, 1879
N 40° 27.3’ W 070° 47.6’
KY007276 to KY007278
Halosydna brevisetosa Kinberg, 1856
N 48° 29.5’ W 123° 01.1’
KY007357 to KY007363
S 14o 40.0' E 145o 27.1'
KY007454 to KY007459
Myxicola infundibulum (Montagu, 1808)
N 48° 28.6’ W 122° 58.7’
KY007387 to KY007391
Schizobranchia insignis Bush, 1905
N 48° 33.3’ W 122° 56.5’
N 40° 50.3’ W 92° 5.3’
KY007461 to KY007468
Boccardia proboscidea Hartman, 1940
N 48° 29.2’ W 123° 04.1’
KY007301 to KY007307
Sternaspis scutata Ranzani, 1817
N 48° 29.1’ W 123° 04.3’
KY007469 to KY007471
Syllis cf. hyalina Grube, 1863
S 34° 37.0’ E 19° 21.4’
Themiste pyroides (Chamberlin, 1919)
N 48° 21.4’ W 123° 43.4’
KY007473 to KY007479
N 36° 41.2’ W 122° 02.0’
KY007486 to KY007488
Oligochaeta gen. sp. (unidentified Crassiclitellata - Place Kabary 2)
S 12° 59.0’ E 49° 17.4’
KY007417 to KY007420
Following trimming alignment of translated transcripts possessed 132 residue positions, with nearly all sequences, the exception being Aphrodita japonica, starting with a methionine residue. We decided to keep the apparently incomplete sequence from A. japonica due to its high similarity with the remaining sequences. All sequences in the alignment contained signature residues involved in iron binding, indicating putative respiratory function for these putative Hrs . For the 214 sequences, 100 were unique and 114 identical for at least two species at the amino acid level.
Sequences were assembled into a final dataset containing 225 sequences being 214 new, two Hr sequences from Lingula reevii, a brachiopod, and nine annelid sequences previous used as “queries”, with 396 aligned nucleotides positions. Of these, 209 sequences contained a five-codon insertion between the C and D α-helices described for myoHr, but not cHrs (Fig. 1) . Datasets supporting conclusions of this article are available in the Figshare repository, under DOIs: https://doi.org/10.6084/m9.figshare.3505883.v1 and https://doi.org/10.6084/m9.figshare.3505886.v1.
Our analysis found 13 putative cHr sequences lacking the characteristic five-codon insertion before the D α-helix that define myoHrs, distributed across nine different families. Besides well-known records of cHrs in sipunculids, such as Themiste pyroides and Phascolosoma agassizii, we discovered cHrs in the sipunculid Thysanocardia nigra as well as six annelid families; Amphinomidae Aproditidae, Capitellidae, Oweniidae, Sabellidae, and Spionidae (Fig. 2; blue clade).
The topology of the Hr gene tree did not mirror recent phylogenies of Annelida based on phylogenomic datasets [27, 38, 39]. For example, we found 10 Hr sequences identical at the nucleotide level (Fig. 2, purple clade) belonging to distant annelid families indicating a strong conservation among those orthologs. Several of these sequences were prepared and sequenced at different times, making cross contamination unlikely. Those 10 identical sequences differed 28.54% (nucleotide level) from the consensus of all others myoHr sequences and the majority of nonsynonymous substitutions are concentrated in A and B α-helices.
Regarding paralogs multiple copies of Hr genes were found for several species, including two paralogs from both Chaetopterus variopedatus and Mesochaetopterus taylori, with these paralogs forming a monophyletic clade (p.p. = 0.85; Fig. 2; green clade). Both species were from Chaetopteridae suggesting a recent paralogous duplication . Given the presence of multiple Hrs, apparently early annelids already contained several copies of Hr genes, with some paralogs arising later (as in C. variopedatus and M. taylori).
As demonstrated here expression of Hr genes among annelids is much more common than previously reported, revealing unrecognized diversity of these genes in this phylum. All 44 of the examined species possessed actively transcribed myoHrs, while cHrs were less frequently recovered from these transcriptomes. Although this diversity and wide distribution of Hrs in annelids could, in part, be explained by the need to carry oxygen (Hr have approximately 25% greater oxygen affinity than hemoglobins; ) secondary functional specializations could also be important for driving diversification. For example, Hrs participate in iron storage, metal detoxification, and immunity in some annelids (e.g., Theromyzon tessulatum, Hirudo medicinalis and Neanthes diversicolor) [33, 34, 43]. Our findings build on recent publications demonstrating that sequence diversity among Hr-bearing species was larger than traditionally suspected [9, 28, 29]. However, those studies were based predominately on genomic data. In this case, use of transcriptomes shows that Hrs genes were present and expressed.
Bayesian phylogenetic reconstruction recovered monophyletic clades for annelid cHrs sequences as well as myoHrs sequences with strong support corroborating Vanin et al.’s  previous findings. The presence of multiple copies of myoHr genes across the annelid phylogeny implies these proteins have undergone several instances of gene duplication during their evolution, as previously reported for other bacterial, archeal, and eukaryotic taxa . Moreover, the unexpected diversity of myoHrs could be associated with functional diversification of this gene, as observed for myoglobins  and also for Hrs involved in heavy metal detoxification and aspects of innate immunity .
Classification of specific Hr subtypes [28, 32–34] was not validated by the gene genealogy. Although our analyses used whole organisms (including reproductive and nerve tissues), our results failed to recover Hr proteins that corresponded to ovoHrs or nHrs. These categories, however, were described based on limited differences in amino acid sequence and do not reflect distinct monophyletic subgroups within the Hr gene family. Given this, recognizing only two primary types of Hrs, circulating Hrs (cHr) and non-circulating Hrs (ncHrs), is perhaps more appropriate. Although tentative, this reinterpretation deserves further consideration.
Our findings demonstrate that sequence diversity among Hr-bearing annelid species is much greater than traditionally suspected and that many of these Hrs are actively expressed. There are two primary types of Hrs circulating Hrs (cHr) and non-circulating Hrs (ncHrs), instead of the four subtypes reported in the literature. Incongruence between our gene genealogy relative to current knowledge of annelid evolutionary history indicates that Hrs have a complex history. Our findings indicate that both cHrs and myoHrs were likely present in the ancestor of annelids, as both subtypes occur in all lineages near the base of the annelid tree.
Information on species employed herein is provided in Table 1. Transcriptomes of these species were collected as part of the WormNet II project to resolve annelid phylogeny and were collected with a variety of techniques, including intertidal sampling, dredge and box cores. Upon collection, all samples were either preserved in RNALater or frozen at −80 °C.
Data collection & sequence assembly
RNA extraction, cDNA preparation and high-throughput sequencing generally followed Kocot et al.  and Whelan et al. . Briefly, total RNA was extracted from either whole animals (for small specimens) or the body wall and coelomic region (for larger specimens). After extraction, RNAs were purified using TRIzol (Invitrogen) or the RNeasy kit (Qiagen) with on-column DNase digestion, respectively. The SMART cDNA Library Construction Kit (Clonetech) was utilized to reverse transcribe single stranded RNA template. Double stranded cDNA synthesis was completed with The Advantage 2 PCR system (Clontech). Libraries were barcoded and sequenced with Illumina technology by The Genomic Services Lab at the Hudson Alpha Institute (Huntsville, Alabama, USA). Because transcriptomic sequencing was conducted from 2012–2015, Paired End (PE) runs were of 100 bp or 125 bp lengths, utilizing either v3 or v4 chemistry on Illumina HiSeq 2000 or 2500 platforms (San Diego, California). To facilitate sequence assembly, paired-end transcriptome data were digitally normalized to an average k-mer coverage of 30 using normalize-by-median.py  and assembled using Trinity r2013-02-25 with default settings .
Data mining and gene identification
Two approaches were utilized in searching transcriptomic data for putative Hr genes in silico. The first approach employed BLASTX  at an e-value cutoff of 10−6 to compare each assembled transcriptome contig (‘queries’) to a protein database composed of 21 protein sequence from the National Center for Biotechnology (NCBI) database (Additional file 3: Table S2) of at least 110 amino acid residues and previously identified as annelid Hrs (n = 11), myoHrs (n = 9), or nHr (n = 1). Sequences of ovoHrs were not included since there were only two relatively short (18 amino acid residues each) sequences available from NCBI. The BLASTX approach, rather than a tBLASTn, assured that any transcriptome contig with a significant ‘hit’ to an Hr would be further evaluated in the pipeline. Contigs recovered from initial BLAST searches were then utilized in BLASTX searches against the NCBI protein database (minimum e-value of 10−10) and only top hits longer than 300 nucleotides retained. These were considered putative cHr, myoHr or nHr genes, as appropriate.
The second approach processed the transcriptomic data through the Trinotate annotation pipeline (http://trinotate.github.io/) , which utilizes a BLAST-based approach to provide, among others, GO (The Gene Ontology Consortium) annotation . Recovered sequences were verified by the functional annotation they received. Transcripts annotated as hemerythrins, using the 10−6 e-value cutoff obtained by using BLASTX, were also considered putative hemerythrin-like gene orthologs.
Contigs putatively identified as Hr genes by the two approaches were subsequently translated into amino acids using TransDecoder with default settings . Since TransDecoder can produce multiple open reading frames (ORFs), all translations were additionally subject to a Pfam Domain evaluation using the EMBL-EBI database with an e-value cutoff of 10−5. Translations returning an Hr-like Pfam domain and that were longer than 100 amino acids residues were retained for subsequent analyses. Moreover, we manually evaluated the presence of residues involved in iron binding, which are: histidine residues (His) in positions 26, 56, 75, 79, and 108; glutamic acid residue (Glu) in position 60; and aspartic acid residue (Asp) in position 113, numbered by reference sequence T. zostericola. Presence of these signature residues indicates putative respiratory function for Hrs. Transcripts passing the criteria described above were considered Hr genes (Table 1). ovoHr has been just reported once for a single species and there is no complete sequence available at GenBank, so we are not able to investigate if ovoHr is present in our dataset.
The protein dataset consisted of 225 sequences, including two Hr sequences from Lingula reevii, a brachiopod, nine annelid sequences previous used as “queries” (Additional file 3: Table S2), and a remaining 214 sequences from translated transcripts (Additional file 4). All sequences were initially aligned with MAFFT using the “accurate E-INS-i” algorithm , followed by visual inspection and manual curation in order to remove spuriously aligned sequences based on similarity to the protein alignment as a whole. Subsequently, ends of aligned sequences were manually trimmed in Geneious 8.1.6  to exclude 5’residues leading to the putative start codon and 3’ residues following the first seven amino acids subsequent to the end of the D α-helix.
In order to employ nucleotide sequences in the phylogenetic analysis, an alignment from corresponding aligned protein sequences (Additional file 4) was performed using PAL2NAL . The resulting nucleotide alignment was used for all subsequent analyses (Additional file 5).
JModelTest2 was applied to carry out statistical selection of best-fit models of nucleotide substitution for the dataset using the Akaike and Bayesian Information Criteria (AIC and BIC, respectively) methods . Bayesian phylogenetic inference was performed with MrBayes 3.2.1  using the GTR + G substitution model. Two independent runs with four Metropolis-coupled chains were run for 107 generations, sampling the posterior distribution every 500 generations. In order to confirm if chains achieved stationary and determine an appropriate burn-in, we evaluated trace plots of all MrBayes parameter output in Tracer v1.6 . The first 25% of samples were discarded as burn-in and a majority rule consensus tree generated using MrBayes. Bayesian posterior probabilities were used for assessing statistical support of each bipartition.
Two alternative approaches were used to root the Hr gene genealogy. Firstly, the tree inferred in MrBayes was rooted using two Hr sequences from Lingula reevii (AAB32422.1 and AAB32423.1) as outgroup. We also inferred the root of the tree using using BEAST 1.8.3  to infer a rooted tree of Hrs under the strict molecular clock. This was done because we were unable to decisively rule out the possibility that Hr sequences from Lingula reevii were closely evolutionarily related to one of the annelid Hr lineages. The strategy using BEAST is similar to midpoint rooting, although the Bayesian implementation in BEAST allows for a more flexible treatment of the evolutionary rate via a normal prior with mean and standard deviation equal to 1. Moreover, tree topology was jointly estimated. In BEAST, we adopted the same substitution model settings used in MrBayes and the MCMC algorithm was run for 50,000,000 generations and sampled every 1,000th generation, with 50% of the run discarded as burn-in. Trees was summarized in TreeAnnotator 1.8.3  and Markov chain stationarity was assessed in Tracer by ESS values > 1,000. Since the inferred root node separated Lingula reevii sequences from annelid sequences with 0.86 posterior probability, the results are reported using the gene genealogy rooted with outgroup.
Evolutionary rate analyses
The protein alignment (Additional file 4) was used in DIVERGE  to examine site-specific shifted evolutionary rates and assesses whether there has been a significant change in evolution rate after duplication or speciation events by calculating the coefficient of divergence (θD) and determining if the null hypothesis of no functional divergence between myoHrs and cHrs could be statistically rejected. We employed a cutoff of 0.7 for detection of site-specific shifted evolutionary rates.
Additionally, we used RELAX  to detect changes in selection intensity between cHRs and myoHrs based on aligned nucleotide data (Additional file 5). Rather than just looking in dN/dS ratios averaged across branches using a ‘branch-site’ model approach, RELAX examines dN/dS ratios along each branch at a given site by drawing values from a discrete distribution of dN/dS ratios independent of other branches. In addition to assessing values in a branch independent manner, RELAX results are analyzed within a phylogenetic framework.
We thank Fernando Avila Queiroz for valuable help with the figures and Antonio Sole-Cava, Joana Zanol, and Francisco Prosdocimi for critically reading this manuscript. We are grateful to Christer Erséus and Sam James for oligochaete samples used in transcriptomes. A Fellowship for the first author was provided by CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior, Brazil). The findings and conclusions in this article are those of the authors and do not necessarily represent the views of the U.S. Fish and Wildlife Service. Use of SkyNet computational resources at Auburn University is acknowledged. This paper is part of the D. Sc. Requirements of Elisa Maria Costa-Paiva at Biodiversity and Evolutionary Biology Graduate Program of the Federal University of Rio de Janeiro. This is Molette Biology Laboratory contribution 61 and Auburn University Marine Biology Program contribution 153.
This work was funded by National Science Foundation grant DEB-1036537 to K. M. Halanych and Scott R. Santos and OCE-1155188 to K. M. Halanych.
Availability of data and material
The alignment datasets supporting the conclusions of this article are available in the Figshare repository, under DOI: https://doi.org/10.6084/m9.figshare.3505883.v1 and https://doi.org/10.6084/m9.figshare.3505886.v1. Individual sequences have been deposited to GenBank under accessions KY007275-KY007488 and are listed in Additional file 2: Table S1.
EMC and KMH conceived this study. EMC performed the in silico analysis helped by DSW. EMC, NVW and CGS participate in phylogenetic analysis. EMC, CGS and KMH interpreted the data. EMC and KMH drafted the manuscript and all others revised critically. All the authors read and approved the version to be published.
The authors declare that they have no competing interests.
Consent for publication
Ethics approval and consent to participate
This work used non vertebrate animals only and as such there were no ethic approval requirements.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
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.
- Schmidt-Rhaesa A. The Evolution of Organs Systems. New York: Oxford University Press; 2007.
- Terwilliger NB. Functional adaptations of oxygen-transport proteins. J Exp Biol. 1998;201:1085–98.PubMedGoogle Scholar
- Terwilliger RC, Terwilliger NB, Schabtach E. Comparison of chlorocruorin and annelid hemoglobin quaternary structures. Comp Biochem Phys A. 1976;55(1):51–55.
- Burmester T. Origin and evolution of arthropod hemocyanins and related proteins. J Comp Physiol B. 2002;172(2):95–107.View ArticlePubMedGoogle Scholar
- Burmester T. Evolution of respiratory proteins across the Pancrustacea. Integr Comp Biol. 2015;55(5):765–70.View ArticleGoogle Scholar
- Lecomte JT, Vuletich DA, Lesk AM. Structural divergence and distant relationships in proteins: evolution of the globins. Curr Opin Struc Biol. 2005;15(3):290–301.View ArticleGoogle Scholar
- Decker H, Hellmann N, Jaenicke E, Lieb B, Meissner U, Markl J. Minireview: Recent progress in hemocyanin research. Integr Comp Biol. 2007;47(4):631–44.View ArticlePubMedGoogle Scholar
- Vinogradov SN, Hoogewijs D, Bailly X, Arredondo-Peter R, Gough J, Dewilde S, Moens L, Vanfleteren JR. A phylogenomic profile of globins. BMC Evol Biol. 2006;6(1):31.View ArticlePubMedPubMed CentralGoogle Scholar
- Vanin S, Negrisolo E, Bailly X, Bubacco L, Beltramini M, Salvato B. Molecular evolution and phylogeny of sipunculan hemerythrins. J Mol Evol. 2006;62:32–41.View ArticlePubMedGoogle Scholar
- Rousselot M, Delpy E, Rochelle CD, Lagente V, Pirow R, Rees JF, Hagege A, Guen D, Hourdez S, Zal F. Arenicola marina extracellular hemoglobin: A new promising blood substitute. Biotechnol J. 2006;1(3):333–45.View ArticlePubMedGoogle Scholar
- Zal F, Rousselot M. Use of a haemoglobin for the preparation of dressings and resulting dressings. U.S. Patent 9220929B2, December 28, 2015.
- Jones L. Recent advances in the molecular design of synthetic vaccines. Nat Chem. 2015;7(12):952–60.View ArticlePubMedGoogle Scholar
- Zal, F. Use of haemoglobin of annelids for treating cancer. U.S. Patent 0374796A1, December 31, 2015.
- Andreeva A, Howorth D, Brenner SE, Hubbard TJ, Chothia C, Murzin AG. SCOP database in 2004: refinements integrate structure and sequence family data. Nucleic Acids Res. 2004;32:226–9.View ArticleGoogle Scholar
- Mangum CP. Physiological function of the hemerythrins. In: Mangum CP, editor. Advances in Comparative & Environmental Physiology Vol 13 – Blood and Tissue Oxygen Carriers. Berlin: Springer; 1992. p. 173–92.View ArticleGoogle Scholar
- Thompson JW, Salahudeen AA, Chollangi S, Ruiz JC, Brautigam CA, Makris TM, Lipscomb JD, Tomchick DR, Bruick RK. Structural and molecular characterization of iron-sensing hemerythrin-like domain within F-box and leucine-rich repeat protein 5 (FBXL5). J Biol Chem. 2012;287(10):7357–65.View ArticlePubMedPubMed CentralGoogle Scholar
- Addison AW, Bruce RE. Chemistry of Phascolosoma lurco hemerythrin. Arch Biochem Biophys. 1977;183:328–32.View ArticlePubMedGoogle Scholar
- Klippenstein GL. Structural aspects of hemerythrin and myohemerythrin. Am Zool. 1980;20:39–51.View ArticleGoogle Scholar
- Klotz IM, Kurtz Jr DM. Binuclear oxygen carriers: hemerythrin. Accounts Chem Res. 1984;17(9):16–22.View ArticleGoogle Scholar
- Stenkamp RE. Dioxygen and hemerythrin. Chem Rev. 1994;94:715–26.View ArticleGoogle Scholar
- Wirstam M, Lippard SJ, Friesner R. Reversible dioxygen binding to hemerythrin. J Am Chem Soc. 2003;125(13):3980–7.View ArticlePubMedGoogle Scholar
- Wells RMG. Respiratory characteristics of the blood pigments of three worms from an intertidal mudflat. New Zeal J Zool. 1982;9:243–8.View ArticleGoogle Scholar
- Kurtz Jr DM. Molecular structure/function of the hemerythrins. In: Mangum CP, editor. Advances in Comparative & Environmental Physiology Vol 13 – Blood and Tissue Oxygen Carriers. Berlin: Springer; 1992. p. 151–71.View ArticleGoogle Scholar
- French CE, Bell JML, Ward FB. Diversity and distribution of hemerythrin-like proteins in prokaryotes. FEMS Microbiol Lett. 2008;279(2):131–45.View ArticlePubMedGoogle Scholar
- Li X, Tao J, Hu X, Chan J, Xiao J, Mi K. A bacterial hemerythrin-like protein MsmHr inhibits the SigF-dependent hydrogen peroxide response in mycobacteria. Front Microbiol. 2015;5(Jan):1–11.Google Scholar
- Alvarez-Carreño C, Becerra A, Lazcano A. Molecular evolution of the oxygen-binding hemerythrin domain. PLoS ONE. 2016;11(6):e0157904.View ArticlePubMedPubMed CentralGoogle Scholar
- Weigert A, Helm C, Meyer M, Nickel B, Arendt D, Hausdorf B, Santos SR, Halanych KM, Purschke G, Bleidorn C, Struck TH. Illuminating the base of the annelid tree using transcriptomics. Mol Biol Evol. 2014;31(6):1391–401.View ArticlePubMedGoogle Scholar
- Bailly X, Vanin S, Chabasse C, Mizuguchi K, Vinogradov SN. A phylogenomic profile of hemerythrins, the nonheme diiron binding respiratory proteins. BMC Evol Biol. 2008;8:244.View ArticlePubMedPubMed CentralGoogle Scholar
- Martín-Durán JM, De Mendoza A, Sebé-Pedrós A, Ruiz-Trillo I, Hejnol A. A broad genomic survey reveals multiple origins and frequent losses in the evolution of respiratory hemerythrins and hemocyanins. Genome Biol Evol. 2013;5:1435–42.View ArticlePubMedPubMed CentralGoogle Scholar
- Klippenstein GL, Cote JL, Ludlam SE. The primary structure of myohemerythrin. Biochemistry. 1966;15(5):1128–36.View ArticleGoogle Scholar
- Ward KB, Hendrickson WA, Klippenstein GL. Quaternary and tertiary structure of hemerythrin. Nature. 1975;257:818–21.View ArticlePubMedGoogle Scholar
- Coutte L, Slomianny MC, Malecha J, Baert JL. Cloning and expression analysis of a cDNA that encodes a leech hemerythrin. Biochim Biophys Acta. 2001;1518(3):282–6.View ArticlePubMedGoogle Scholar
- Baert JL, Britel M, Sautière P, Malécha J. Ovohemerythrin, a major 14-kDa yolk protein distinct from vitellogenin in leech. Eur J Biochem. 1992;209:563–9.View ArticlePubMedGoogle Scholar
- Vergote D, Sautière PE, Vandenbulcke F, Vieau D, Mitta G, Macagno ER, Salzet M. Up-regulation of neurohemerythrin expression in the central nervous system of the medicinal leech, Hirudo medicinalis, following septic injury. J Biol Chem. 2004;279(42):43828–37.View ArticlePubMedGoogle Scholar
- Mangum CP. Major events in the evolution of the oxygen carriers. Amer Zool. 1998;38(1):1–13.View ArticleGoogle Scholar
- Manwell C, Baker CMA. Magelona haemerythrin: tissue specificity, molecular weights and oxygen equilibria. Comp Biochem Phys B. 1988;89(3):453–63.Google Scholar
- Rouse G, Pleijel F. Polychaetes. New York: Oxford University Press; 2001.
- Struck TH, Golombek A, Weigert A, Franke FA, Westheide W, Purschke G, Bleidorn C, Halanych KM. The evolution of annelids reveals two adaptive routes to the interstitial realm. Curr Biol. 2015;25(15):1993–9.View ArticlePubMedGoogle Scholar
- Weigert A, Bleidorn C. Current status of annelid phylogeny. Org Divers Evol. 2016;16:1–18.View ArticleGoogle Scholar
- DeSalle R. Can single protein and protein family phylogenies be resolved better? J Phylogenetics Evol Biol. 2015;3:e116.View ArticleGoogle Scholar
- Gu X, Zou Y, Su Z, Huang W, Zhou Z, Arendsee A, Zeng Y. An update of DIVERGE software for functional divergence analysis of protein family. Mol Biol Evol. 2013;30:1713–9.View ArticlePubMedGoogle Scholar
- Wertheim JO, Murrell R, Smith MD, Kosakovsky Pond SL, Scheffler K. RELAX: Detecting relaxed selection in a phylogenetic framework. Mol Biol Evol. 2014;32:820–32.View ArticlePubMedPubMed CentralGoogle Scholar
- Demuynck S, Li KW, Van der Schors R, Dhainaut-Courtois N. Amino acid sequence of the small cadmium-binding protein (MP-II) from Nereis diversicolor (Annelida, Polychaeta) – evidence for a myohemerythrin structure. Eur J Biochem. 1993;217:151–6.View ArticlePubMedGoogle Scholar
- Koch J, Lüdemann J, Spies R, Last M, Amemiya CT, Burmester T. Unusual diversity of myoglobin genes in the lungfish. Mol Biol Evol. 2016;33(12):3033–41.View ArticlePubMedGoogle Scholar
- Coates CJ, Decker H. Immunological properties of oxygen-transport proteins: hemoglobin, hemocyanin and hemerythrin. Cell Mol Life Sci. 2017;74:293.View ArticlePubMedGoogle Scholar
- Liu J, Ou Q, Han J, Li J, Wu Y, Jiao G, He T. Lower Cambrian polychaete from China sheds light on early annelid evolution. Naturwissenschaften. 2015;102(5):1–7.Google Scholar
- Kocot KM, Cannon JT, Todt C, Citarella MR, Kohn AB, Meyer A, Santos SR, Schander C, Moroz LL, Lieb B, Halanych KM. Phylogenomics reveals deep molluscan relationships. Nature. 2011;477:452–6.View ArticlePubMedPubMed CentralGoogle Scholar
- Whelan NV, Kocot KM, Moroz LL, Halanych KM. Error, signal, and the placement of Ctenophora sister to all other animals. Proc Natl Acad Sci U S A. 2015;112:5773–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Brown CT, Howe A, Zhang Q, Pyrkosz AB, Brom TH. A reference-free algorithm for computational normalization of shotgun sequencing data. arXiv:1203.4802 [q-bio.GN]. 2012.
- Grabherr MG, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29:644–52.View ArticlePubMedPubMed CentralGoogle Scholar
- Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403–10.View ArticlePubMedGoogle Scholar
- Gene Ontology Consortium. The Gene Ontology (GO) database and informatics resource. Nucleic Acids Res. 2004;32(1):D258–61.View ArticleGoogle Scholar
- Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, Couger MB, Eccles D, Li B, Lieber M, MacManes MD, Ott M, Orvis J, Pochet N, Strozzi F, Weeks N, Westerman R, William T, Dewey CN, Henschel R, LeDuc RD, Friedman N, Regev A. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat Protoc. 2013;8(8):1494–512.View ArticlePubMedGoogle Scholar
- Katoh K, Standley DM. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol. 2013;30(4):772–80.View ArticlePubMedPubMed CentralGoogle Scholar
- Kearse M, Moir R, Wilson A, Stones-Havas S, Cheung M, Sturrock S, Buxton S, Cooper A, Markowitz S, Duran C, Thiever T, Ashton B, Meintjes P, Drummond A. Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics. 2012;28(12):1647–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Suyama M, Torrents D, Bork P. PAL2NAL: robust conversion of protein sequence alignments into the corresponding codon alignments. Nucleic Acids Res. 2006;34 suppl 2:609–12.View ArticleGoogle Scholar
- Darriba D, Taboada GL, Doallo R, Posada D. jModelTest 2: more models, new heuristics and parallel computing. Nat Methods. 2012;9(8):772.View ArticlePubMedPubMed CentralGoogle Scholar
- Ronquist F, Huelsenbeck JP. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics. 2003;19(12):1572–4.View ArticlePubMedGoogle Scholar
- Rambaut A, Suchard MA, Xie D, Drummond AJ. Tracer v1.6. Available from http://beast.bio.ed.ac.uk/Tracer. 2014. Accessed 20 July 2016.
- Drummond AJ, Suchard MA, Xie D, Rambaut A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol. 2012;29:1969–73.View ArticlePubMedPubMed CentralGoogle Scholar
- Yang J, Yan R, Roy A, Xu D, Poisson J, Zhang Y. The I-TASSER Suite: Protein structure and function prediction. Nat Methods. 2015;12:7–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Kocot KM, Struck TH, Merkel J, Waits DS, Todt C, Brannock PM, Weese DA, Cannon JT, Moroz LL, Lieb B, Halanych KM. Phylogenomics of Lophotrochozoa with consideration of systematic error. Syst Biol. 2016. doi:10.1093/sysbio/syw079.