- Research article
- Open Access
Tunicate mitogenomics and phylogenetics: peculiarities of the Herdmania momus mitochondrial genome and support for the new chordate phylogeny
- Tiratha Raj Singh†1,
- Georgia Tsagkogeorga†2,
- Frédéric Delsuc2,
- Samuel Blanquart3,
- Noa Shenkar1, 4,
- Yossi Loya1,
- Emmanuel JP Douzery2Email author and
- Dorothée Huchon1Email author
© Singh et al; licensee BioMed Central Ltd. 2009
- Received: 2 September 2009
- Accepted: 17 November 2009
- Published: 17 November 2009
Tunicates represent a key metazoan group as the sister-group of vertebrates within chordates. The six complete mitochondrial genomes available so far for tunicates have revealed distinctive features. Extensive gene rearrangements and particularly high evolutionary rates have been evidenced with regard to other chordates. This peculiar evolutionary dynamics has hampered the reconstruction of tunicate phylogenetic relationships within chordates based on mitogenomic data.
In order to further understand the atypical evolutionary dynamics of the mitochondrial genome of tunicates, we determined the complete sequence of the solitary ascidian Herdmania momus. This genome from a stolidobranch ascidian presents the typical tunicate gene content with 13 protein-coding genes, 2 rRNAs and 24 tRNAs which are all encoded on the same strand. However, it also presents a novel gene arrangement, highlighting the extreme plasticity of gene order observed in tunicate mitochondrial genomes. Probabilistic phylogenetic inferences were conducted on the concatenation of the 13 mitochondrial protein-coding genes from representatives of major metazoan phyla. We show that whereas standard homogeneous amino acid models support an artefactual sister position of tunicates relative to all other bilaterians, the CAT and CAT+BP site- and time-heterogeneous mixture models place tunicates as the sister-group of vertebrates within monophyletic chordates. Moreover, the reference phylogeny indicates that tunicate mitochondrial genomes have experienced a drastic acceleration in their evolutionary rate that equally affects protein-coding and ribosomal-RNA genes.
This is the first mitogenomic study supporting the new chordate phylogeny revealed by recent phylogenomic analyses. It illustrates the beneficial effects of an increased taxon sampling coupled with the use of more realistic amino acid substitution models for the reconstruction of animal phylogeny.
- Markov Chain Monte Carlo
- Mitochondrial Genome
- Complete Mitochondrial Genome
- Phylogenomic Study
- Halocynthia Roretzi
Mitochondrial genomes (mtDNAs) of bilaterian animals are short, circular DNA molecules of 14-16 kb in length, typically characterized by the absence of introns and the presence of only short intergenic regions, with the exception of the control region (CR), a non-coding region assumed to contain the elements for the regulation of replication and transcription of the genome [1, 2]. To date, about 1,700 complete metazoan mitochondrial genomes have been sequenced and used in comparative mitogenomics and phylogenetic studies on different taxonomic scales [3–9].
The mitochondrial gene content is highly conserved across the different bilaterian phyla, with typically 37 genes [1, 2, 10]. Among them, 13 genes encode for proteins [ATP synthase subunits 6 and 8 (atp6 and atp8), cytochrome oxidase subunits (cox1, cox2, cox3), apocytochrome b (cytb), and dehydrogenase subunits (nd1, nd2, nd3, nd4, nd5, nd6 and nd4L)]. The remaining genes encode two ribosomal subunits (srRNA and lrRNA) and usually 22 tRNAs. However, cases of duplication and loss of tRNAs have been reported within bilaterians [11, 12]. Tunicate mitochondrial genomes illustrate such exceptions and typically encode 24 tRNAs, apart from two Phallusia species which lack the tRNA-Asp , and Halocynthia roretzi which encode two tRNA-phe . The two additional tRNAs present in tunicate mitochondrial genomes when compared to vertebrates are (i) the tRNA-Gly (for AGR codons), which is necessary for the translation due to the derived tunicate mitochondrial genetic code , and (ii) the tRNA-Met (for AUA codon), whose presence might reduce the conflict between translation initiation -- which requires a tRNA-Met (for AUG) -- and translation elongation that involves AUG codons .
The mitochondrial gene order is highly conserved within Deuterostomia [10, 17], and particularly in chordate genomes. Conversely, mitochondrial gene arrangement shows an important plasticity in some animal phyla, (e.g. molluscs and nematodes [1, 2]), and especially in tunicates [2, 13, 18–22]. Tunicates, or Urochordates, are marine deuterostomes characterized by markedly diversified developmental and life history traits, and traditionally encompass three major classes: Ascidiacea (sea squirts), Thaliacea (salps) and Appendicularia (larvaceans). Ascidiacea, commonly referred to as ascidians, is the most speciose and widespread group. Several ascidian species have been identified as invasive species, such as Styela clava and Pyura praeputialis [23–26], and have a strong ecological impact on the invaded marine ecosystems. Some species are also widely used as model organisms in evo-devo studies like Ciona intestinalis and Botryllus schlosseri [27–31]. According to the traditional classification, the class Ascidiacea is subdivided into three major orders: Phlebobranchia, Aplousobranchia, and Stolidobranchia. In contrast with this taxonomic view, 18S rRNA-based phylogenies have shown that ascidians are in fact paraphyletic [32–34]. According to the 18S rRNA phylogenetic framework, Aplousobranchia, Phlebobranchia, and Thaliacea are closely related, whereas Stolidobranchia forms a distinct and monophyletic group, which might be close to Appendicularia, although the position of the latter is still debated [32, 33].
To date, complete mitochondrial genomes of tunicates are mainly available for a single representative of Thaliacea (Doliolum nationalis) and five ascidians [13, 18–22], including four phlebobranchians (Ciona intestinalis type A and B, C. savignyi, Phallusia fumigata, P. mammillata) and one stolidobranchian (Halocynthia roretzi). The available mitochondrial data suggest that several unique features characterize mitogenomic evolution in tunicates relative to other chordate phyla. Two main peculiarities can be distinguished. The first refers to the highly variable gene order observed within the group, which implies that extensive gene rearrangements have occurred even at low taxonomic levels [13, 18, 19]. However, since most available tunicate complete mtDNA sequences belong to phylogenetically-related species (except H. roretzi) according to the 18S rRNA reference [32–34], it is not possible to evaluate whether mitochondrial gene rearrangements characterize the whole order or only the Aplousobranchia + Phlebobranchia + Thaliacea clade. The second specificity is that of an accelerated evolutionary rate of tunicates, as revealed by the long branches of the group in mitogenomic topologies [35–37] and the associated composition bias . However, whether this accelerated substitution rate is restricted to protein coding genes as in snakes  or is a more general feature of the whole mtDNA of Tunicates, has yet to be investigated.
These two peculiar evolutionary features of tunicate mitochondrial genome evolution have hampered their reliable phylogenetic placement within metazoans. Analyses of mitochondrial protein-coding genes have almost always systematically placed tunicates outside Bilateria [20, 35, 36, 39]. This is in sharp contrast with recent nuclear-based phylogenomic studies that identified tunicates as the closest living relatives of vertebrates within chordates [40–43]. Only two recent mitogenomic study have found marginal support for chordate monophyly. Bourlat et al.  grouped cephalochordates with vertebrates according to the traditional Euchordata hypothesis using a concatenation of the 13 protein coding genes under a site- and time-heterogeneous mixture model in Bayesian phylogenetic reconstructions. Alternatively, Zhong et al.  recovered the new chordate phylogeny under the maximum likelihood framework when removing the fastest evolving vertebrates species and when considering only the four most conserved mitochondrial proteins.
Here, we sequenced the complete mitochondrial genome of the solitary ascidian Herdmania momus (Ascidiacea: Stolidobranchia: Pyuridae), an Indo-Pacific species that was introduced into the Mediterranean Sea via the Suez Canal . We describe the structural and compositional features of H. momus mtDNA, discuss its evolutionary dynamics with respect to other tunicate and chordate mitochondrial genomes, and provide an updated metazoan phylogeny based on probabilistic analyses of the 13 mitochondrial proteins using site- and time-heterogeneous mixture models of amino acid substitutions.
General features of H. momus mtDNA
In contrast to other deuterostome mtDNA genomes, all genes are encoded on the same DNA strand and thus mitochondrial gene transcription shows only one orientation in H. momus, as in all other tunicate genomes sequenced to date. In chordates, although the majority of genes is co-orientated on one major DNA strand (H-strand), a minimum set of tRNAs and the nd6 gene are encoded from the minor one (L-strand) [2, 17, 45]. The genome architecture of H. momus thus provides additional evidence in support of the hypothesis that gene arrangement in only one strand is a shared derived feature of tunicate mtDNAs .
With regard to genome organization, H. momus mtDNA seems substantially less compact than other available tunicate mitochondrial genomes. Adjacent genes overlap in four cases: cox1 - tRNA-Val (9 bp), nd4L - tRNA-Leu (15 bp), atp8 - tRNA-Asp (1 bp) and nd2 - tRNA-Trp (2 bp). The total gene overlap thus appears slightly lower with respect to the average ~35 bp gene overlap of the other tunicate genomes . Interestingly, there is no protein-protein gene overlap. Furthermore, H. momus shows a high percentage of non-coding (NC) sequence (ca. 10%).
Concerning base composition, the mitochondrial genome of H. momus is AT rich with a value of 63.5%. The overall base frequencies as estimated in the coding strand follow the decreasing order T, G, A and C. T is by far the most abundant nucleotide with 41.3% in the genome. Most tunicate mtDNAs present similar compositional AT-rich profiles varying from 61.1% in D. nationalis to 78.6% in C. intestinalis. The only exception concerns the mitochondrial genomes of species from the genus Phallusia, which present more balanced AT and GC compositions with 52.8% and 53.2% AT in P. mammillata and P. fumigata, respectively .
The predicted protein-coding genes in the mtDNA of H. momus present lengths that are overall similar to their orthologues in the other tunicate genomes. Nevertheless, H. momus was not predicted to carry any deletion at the C-terminal ends of nd4 and nd5 genes similar to those reported in H. roretzi .
The base composition of protein-coding genes in H. momus was measured as being A+T rich, as is the profile of its entire genome. Consequently, the codon usage follows the same pattern, with C-rich codons showing lower frequencies than those that are T- or G- rich. Extreme cases are constituted by the CGC codon for Arginine, which is never used, and by the UUU codon for Phenylalanine and UUG codon for Leucine, which are the most frequently used. The CGC codon is also never used in the two Ciona genomes, whereas it is present in other tunicate genomes [18, 22].
The preferred start and stop codons in H. momus mitochondrial genes are GTG (6 genes) and TAG (7 genes), respectively. In the other Pyuridae (H. roretzi), GTG is the preferred start codon too, but TAA is the most frequent stop codon (8 genes). Like in other tunicate genomes, incomplete T-starting stop codons are predicted at the end of the three genes coding for nd1, nd4 and nd4L. The complete TAA stop codon is probably created by polyadenylation.
Because of high sequence divergence of tunicate rRNA sequences, the boundaries of all tunicate mitochondrial rRNA sequenced so far have been inferred from the flanking genes. Likewise, the srRNA gene in H. momus is estimated at most as 673 bp long, thus being slightly shorter than its acidian orthologues which range from 687 to 738 bp, and about 30 bp longer than the srRNA of D. nationalis mtDNA. On the genome map, it is found located between the tRNA-Ala and tRNA-Glu. The lrRNA, on the other hand, appears to be at most 1,159 bp long, a length similar to that of other tunicates, and is located downstream of tRNA-Tyr and upstream of cytb (Figure 1).
Both lengths and locations of the two rRNA genes in the H. momus mitochondrial genome provide additional evidence for two unique features of mtDNA evolution in tunicates. The first is that tunicates present the shortest mitochondrial rRNA genes among chordates, and more generally among deuterostomes, with lengths varying at most from 641 to 738 bp for the small ribosomal subunit and from 1,059 to 1,279 bp for the large ribosomal subunit . The corresponding lengths for the srRNA and lrRNA genes in the other chordate lineages are estimated as higher than 844 bp and 1,367 bp, respectively. The second feature shared by all tunicates concerns the location of the two genes in the mitochondrial genome. In contrast to all other chordates, where the two rRNA coding genes are usually adjacent, srRNA and lrRNA are found rearranged and separated in all available tunicate mtDNAs . The distance between the two genes in H. momus genome is about 5 kb. Distances of the same order of magnitude (5-7 kb) were also identified in the mtDNAs of Cnidaria, some Protostomia, and some Echinodermata .
The predicted cloverleaf structures for the 24 tRNA genes of the H. momus mitochondrial genome are presented in Figure 2. As in H. roretzi, only two tRNAs are characterized by an unusual structure: tRNA-Ser (AGY) and tRNA-Asn. tRNA-Ser (AGY) appears to possess an 8 bp-long anticodon stem structure. Similar long stems (7-9 bp) have been reported in Ciona ssp., D. nationalis, and H. roretzi tRNA-Ser (AGY) [18, 20–22]. However, this structure is absent in Phallusia tRNA-Ser (AGY), which appears to have retained the classical cloverleaf conformation . Concerning the tRNA-Asn gene, two alternative cloverleaf structures have been inferred, as also previously reported for H. roretzi . The structure predicted by tRNA-scan (Figure 2) is characterized by a two-nucleotide spacer between the DHU and the anti-codon stem as well as a mismatch at the start of the T stem. The alternative conventional cloverleaf structure (Figure 2) appears less stable, in forming a mismatch at the start of the D arm. Similar two-nucleotide spacers have also been found in the tRNA-Asn structures of Phallusia ssp., H. roretzi, and D. nationalis, but not in Ciona ssp [13, 18, 20–22].
Statistics on non-coding (NC) sequences in tunicate mitochondrial genomes.
Longest NC region
nd1 - lr RNA
nd1 - tRNA Pro
cox1 - tRNA Gly (GGN)
tRNA Cys - nd4L
cox3 - cox2
nd4 - tRNA Val
nd6 - tRNA Tyr
Considering length distributions, 19 of the 28 NC regions present lengths over 20 bp, with the longest accounting for 197 bp and being located between the nd6 and tRNA-Tyr genes (Figure 1). Three additional NC regions measuring more than 100 bp have been identified in H. momus mtDNA, involving the disjunction of the following gene pairs: cytb - nd2, nd4 - tRNA-His, and tRNA-Lys - nd3. The remaining nine NC sequences are shorter intergenic spacers, mostly 3-10 bp long, distributed homogenously through the genome.
The presence of palindrome sequences was checked by sequence similarity searches and has been detected within 12 NC regions, and in particular within the three longest NC regions of the H. momus mitochondrial genome. However, no stem-loop structure similar to that of the control region of vertebrate mt-genomes has been identified in H. momus. Further studies are necessary to establish whether these regions are involved in the control of replication and transcription.
Sequence similarities searches have finally revealed the occurrence of a duplication event in the NC regions of the H. momus genome, similar to those previously determined in Ciona intestinalis and Phallusia fumigata mtDNAs . Here, a fragment of 89 bp was found duplicated in two different locations of the genome: the first, situated in the non-coding region downstream of the nd6 gene and upstream of the gene encoding for the tRNA-Tyr (NC = 197 bp); while the second involved the non-coding region downstream of the tRNA-Asn (NC = 88 bp) and 13 bp of the 5' region of the tRNA-Lys. The estimated distance between the repeated motifs is about 5 kb (Figure 1).
Thus, the mitochondrial genome organization in H. momus when compared to that of H. roretzi suggests that extensive gene rearrangements occur within Stolidobranchia, similar to those observed in Phlebobranchia . Future work should determine whether co-generic species in Stolidobranchia also present high variations in gene order, as reported for the genera Ciona and Phallusia within Phlebobranchia .
The mitochondrial genome has been repeatedly used in molecular phylogenetics of animals, principally due to several convenient features, such as its relatively small size, its cellular abundance, and its mostly uniparental mode of inheritance facilitating orthology assignment . However, the phylogenetic resolving power of mtDNAs is often reduced by pervasive lineage-specific variations of base composition and substitution rate [46, 47]. Bilaterian evolutionary relationships obtained from previous phylogenetic reconstructions based on mitochondrial genes [35, 36] are in relatively good agreement with the currently accepted view of Protostomia phylogeny . However, the situation is more complex concerning Deuterostomia, mainly because of the recurrent placement of tunicates as sister-group of either the remaining Bilateria [35, 36] or the other chordates . Such placement outside chordates is in direct contradiction to results from recent large-scale phylogenomic studies, which strongly support the clade Olfactores, i.e., the grouping of tunicates with craniates (jawed vertebrates + cyclostomes) to the exclusion of cephalochordates [40–43]. The sister position of tunicates relative to the remaining Bilateria is generally interpreted as a long-branch attraction artifact caused by the peculiarities of tunicate mitogenomic evolution in terms of both lineage-specific evolutionary rate and amino-acid composition [36, 37, 39].
The maximum likelihood tree obtained on the 13 protein-coding genes concatenation using the site- and time-homogeneous mtREV+Γ4 model of amino-acid sequence evolution places tunicates as the sister group of the remaining Bilateria, as previously reported [35, 36, 39]. However, the bootstrap support for such a position is low (BP = 61). This result is presumably the outcome of two types of systematic biases: compositional heterogeneity of the tunicate lineage compared to other chordate and animal phyla (see Figure 4) combined with accelerated evolutionary rates [37, 39]. ML phylogenetic analyses were also conducted under an empirical site-heterogeneous CAT mixture model  which has been shown to alleviate long-branch attraction artefacts in the context of animal phylogenomics . The ML tree inferred under this model shows a dramatic topological shift in placing tunicates as the sister-group of craniates within monophyletic chordates and deuterostomes, in agreement with phylogenomic studies [35, 41]. However, this phylogenetic placement is not statistically supported (BP < 50).
Bayesian posterior probabilities (PP) for alternative positions of tunicates using different models of amino-acid replacements.
Tunicata sister to the remaining Bilateria
Tunicata sister to the remaining Deuterostomia
Although removal of fast-evolving gene and taxa has been found to recover the new chordate phylogeny, such an approach is not successful to solve relationships within Tunicata . Our mitogenomic analysis provided strong evidence for intra-tunicate relationships, with the overall phylogenetic scheme being well resolved despite the high lineage-specific evolutionary rate (Figure 5). The inferred relationships for the group were notably in concordance with nuclear-based phylogenies and morphological data. More precisely, Herdmania was unambiguously retrieved in a sister-group relationship with Halocynthia (PP = 1), in agreement with the systematic classification that places these two ascidian genera within the stolidobranch family Pyuridae . Firm evidence was similarly obtained for the respective grouping of the congeneric species of Ciona and Phallusia (PP = 1), and high posterior probability is retrieved for a sister-group relationship of the two genera into Phlebobranchia (PP = 1). Finally, the thaliacean Doliolum is found to branch within ascidians with solid statistical support (PP = 0.99), favoring the hypothesis of Ascidiacea paraphyly  as also supported by 18S rRNA phylogenies [32, 34, 56].
Compared evolutionary rates of rRNA- and protein-coding genes
The complete mitochondrial genome of the ascidian H. momus shares several features with the other tunicate mtDNAs so far sequenced: (i) all genes are encoded on the same strand; (ii) two additional tRNAs, when compared with vertebrates, are present (one for Gly, one for Met); (iii) the gene order is highly rearranged, and differs from all other known tunicates, with the noticeable synteny disruption of the cox2--cytb gene block; and (iv) the evolutionary rate is elevated and involves both protein- and rRNA-coding genes. From the phylogenetic viewpoint, the analysis of the concatenated mitochondrial proteins recovers the tunicates + craniates clade within monophyletic chordates, provided that site- and time-heterogeneous mixture models of amino acid replacements are used.
Sampling and DNA extraction
The H. momus (Savigny, 1816) individual used for our study was collected in June 2006 on the Eilat-Ashkelon pipeline pier in Eilat, Red Sea 29°31' N 34°55' E at a depth of 12 m . Genomic DNA was isolated from gonads by placing ~0.5 cm3 of tissue in 1 ml of lysis buffer (100 mM Tric-HCl pH = 7.8; 100 mM EDTA; 0.5% SDS; 0.2 mg/ml proteinase K) overnight at 55°C. The digestion was followed by Phenol:Chloroform:Isoamyl alcohol extraction (25:24:1) and isopropanol precipitation .
DNA sequencing and gene annotation
Primers were designed to match conserved regions of chordate genomes in the cox1, cox2, cox3, and cytb genes (Table S1 in Additional File 1). PCR amplifications were performed using the ExTaq polymerase (TaKaRa) in 25 μL reaction mixture according to manufacturer's instructions and were conducted in two steps. A first amplification was performed with external primers followed by a re-amplification of the initial PCR product using different sets of internal primers. Based on the partial sequences obtained for the aforementioned genes, longer specific primers were subsequently designed. This allowed us to determine the cox3-cox1 region and to redesign (long) primers in the nd3 gene located in the same region.
The complete mitochondrial genome was finally amplified in two overlapping fragments using the Long and Accurate (LA) Taq polymerase (TaKaRa) in 25 μL reaction mixture according to manufacturer's instructions. The first fragment covered the cox3-cytb region (~7,700 bp) and was obtained by PCR amplification using the external primers HMO_Cox3_Long_F1 5'-ACTGTTGTTTTACTTAGTTCGGGAGTTACTGTGAC-3' and HMO_Cob_Long_R1 5'-ACATATAAGCGACCCCCAACAAAAAGAAAC-3' followed by a re-amplification with the primers HMO_Cox3_Long_F2 5'-AGTTTTGGTGGAGGGCTTATGGGATTAGTTTGGAC-3' and HMO_Cob_Long_R2 5'-GACACCAGAATAGGGCCAAAAATAAATCTTT-3'. The second fragment encompassed the cytb-nd3 region (~8,900 bp) and was similarly acquired by a first PCR amplification using the external primers HMO_Cob_Long_F1 5'-AGGGGTGCTACTTTAACGCGGTTTTACAC-3' and HMO_Cox1_Long_R2 5'-ACATTATACAACTGCCCATCTCCAATTACCTG-3' and a subsequent re-amplification using the primers HMO_Cob_Long_F2 5'-AGTGGTTTTGTCTTTGGTGCATTTGGTTTTC-3' and NADH3_R 5'-CTGAATAGAATAACCCTCGCTATCACT-3'. The two fragments overlapped over ~150 bp at 5' end of the cytb and over ~780 bp at the cox3-nd3 region. The sequencing of the long PCR products was performed by the company MACROGEN (Seoul, Korea) using shotgun sequencing (i.e., shearing of the PCR fragment, bacterial library construction, and sequencing of the clones). Contigs were assembled using both ARACHNE  and Chromaspro (Technilysium Pty. Ltd.) programs. Regions of low coverage in the assembly (i.e., fewer than three clones) were confirmed by primer-walking sequencing (~4 kb of the genome was re-sequenced) with the use of specific primers and genomic DNA as template. The complete mtDNA sequence of H. momus was deposited under the EMBL accession number FN296153.
Protein coding genes were identified using the ORF-Finder tool of NCBI  and the DOGMA server . Following Iannelli et al. , boundaries of protein coding gene were inferred in such a way that overlap between genes was minimized and similarity between tunicate proteins was maximized. ATG and non-standard initiation codons  were considered as reliable start codons.
Both the tRNAscan-SE 1.23 program and the DOGMA server were used to identify and infer the secondary structures of tRNA genes . tRNAs that were not identified by the above tools were sought by folding all putative non-coding regions using the Mfold server . All predicted cloverleaf secondary structures were manually checked and compared with other known ascidian tRNAs. The positions of the small (srRNA) and large (lrRNA) ribosomal RNA genes were identified by the DOGMA server and confirmed based on sequence similarity searches to orthologous genes in other ascidians. Finally, the boundaries of the rRNA genes were inferred from the flanking genes.
The phylogenetic dataset was built upon the taxon sampling used by Bourlat et al. . We expanded the tunicate sampling to seven by including the newly-sequenced H. momus and two recently published genomes from the genus Phallusia . The snake representative (Boa constrictor) was excluded because snake mitochondrial genomes have been shown to be the subject of adaptive evolution that can severely bias phylogenetic inference . The final dataset therefore comprises 54 taxa including 31 deuterostomes, 19 protostomes, and four non-bilaterian outgroups (cf. Additional file 1). The nucleotide sequences of all mitochondrial genes were retrieved from the Organellar Genome Retrieval (OGRe) database . For the 13 protein coding genes, sequences were translated, and aligned at the amino-acid level using MAFFT with default parameters . Ambiguously aligned sites were identified and removed from all individual genes separately, using the program Gblocks  with the following parameters: minimum number of sequences for a conserved position = 18; minimum number of sequences for a flanking position = 29; maximum number of contiguous nonconserved positions = 8; minimum length of a block = 2; allowed gap positions = with half. The concatenation of the 13 proteins yielded a phylogenetic dataset including 54 taxa and 2136 unambiguously aligned amino-acid sites, of which 1697 were variable (Additional file 2).
Phylogenetic analyses were conducted using Maximum Likelihood (ML) and Bayesian Inference (BI) reconstruction approaches on the concatenated amino acid dataset. Maximum Likelihood (ML) analyses were performed using the program PHYML 3  under the models mtREV+Γ4 and CAT+Γ4, with a 4-category Gamma (Γ4) distribution of the among-site amino acid replacement rate heterogeneity , and with the number of CAT categories set to 20 (C20) as recommended by Le et al. . The heuristic ML searches were conducted by performing Subtree Pruning and Regrafting (SPR) moves on a Neighbor-Joining (NJ) starting tree. Statistical support was estimated by Bootstrap resampling with 100 pseudo-replicates generated by the program SeqBoot of the PHYLIP package . In all replicates, ML analyses were performed using PHYML through the same heuristic search strategy. Nodal bootstrap supports (BS) were obtained from the 50% majority rule consensus of the 100 reconstructed trees using the program TREEFINDER .
Bayesian inference was conducted using the program PhyloBayes 3.1  under both homogeneous amino acid models (mtREV+Γ4 and GTR+Γ4) and variants of the site-heterogeneous CAT+Γ4 mixture model (CAT+Poisson+Γ4, CAT+mtREV+Γ4, and CAT+GTR+Γ4) . The program nhPhyloBayes  was used to performed Bayesian analysis under a site- and time-heterogeneous model which combines the break-point approach (BP) in order to model variations of amino acid replacement rates along branches and the CAT mixture model in order to account for site-wise variations of these rates. In this CAT+BP+Γ4 analysis, the number of categories of the mixture component was fixed at 60 (C60), and the biochemical profiles were those inferred by Le et al. , rather than being estimated.
In Bayesian analyses under mtREV+Γ4, GTR+Γ4, and under the different variants of CAT+Γ4 (CAT+Poisson+Γ4, CAT+mtREV+Γ4, and CAT+GTR+Γ4), three independent Markov Chain Monte Carlo (MCMC) were run in parallel, whereas five MCMC were run for analyses under CAT+BP+Γ4. Each MCMC was launched from a random initialization and for a large number of cycles with parameters and trees saved every cycle. Priors were set to values as described in Blanquart and Lartillot  concerning the CAT+BP model, and in Lartillot et al.  concerning all other applied models. Convergence of MCMC was checked by monitoring the marginal likelihood through cycles. Bayesian Posterior Probabilities (PP) were obtained from the trees sampled during the stationary phase of the different MCMC.
Branch length analysis
Two concatenated nucleotide datasets were constructed to compare evolutionary rates of protein-coding and rRNA-coding mitochondrial genes. In both cases, the different genes were individually aligned using MAFFT, and ambiguously aligned sites were removed using Gblocks before building each concatenation. This led to a 13 protein-coding gene dataset containing 8,130 nucleotide sites and a rRNA dataset including 806 sites (Additional files 3 and 4). Using the rooted topology obtained from the analysis of the amino-acid sequences under the CAT+BP+Γ4 mixture model as a reference, ML branch lengths were optimized under the GTR+Γ8 model using PAUP* 4.0b10  on both the protein-coding and rRNA-coding nucleotide concatenations. For each dataset, the root-to-tip distance was calculated for each of the 54 taxa by summing the branch lengths on the path going from the root of the tree to each of its terminals. The two sets of root-to-tip distances were then compared graphically to evaluate the lineage-specific rate correlation between the two mitochondrial partitions.
We would like to thank Tamar Feldstein and Isadora Cohen for their help at the bench, Benoît Nabholz for sharing scripts, two anonymous reviewers for their comments, and Naomi Paz for editing the text. TRS was supported by a fellowship from the Israel National Committee for Planning and Budgeting Universities (VATAT). This work has been funded by the Research Networks Program in BIOINFORMATICS of the High Council for Scientific and Technological Cooperation between France and Israel to DH, EJPD and YL, and have benefited from the ISE-M bioinformatics cluster. This publication is contribution N° 2009-109 of the Institut des Sciences de l'Evolution de Montpellier (UMR 5554 - CNRS).
- Boore JL: Animal mitochondrial genomes. Nucleic Acids Res. 1999, 27 (8): 1767-1780. 10.1093/nar/27.8.1767.PubMed CentralView ArticleGoogle Scholar
- Gissi C, Iannelli F, Pesole G: Evolution of the mitochondrial genome of Metazoa as exemplified by comparison of congeneric species. Heredity. 2008, 101 (4): 301-320. 10.1038/hdy.2008.62.View ArticleGoogle Scholar
- Boore JL: The use of genome-level characters for phylogenetic reconstruction. Trends Ecol Evol. 2006, 21 (8): 439-446. 10.1016/j.tree.2006.05.009.View ArticleGoogle Scholar
- Boore JL, Brown WM: Big trees from little genomes: mitochondrial gene order as a phylogenetic tool. Curr Opin Genet Dev. 1998, 8 (6): 668-674. 10.1016/S0959-437X(98)80035-X.View ArticleGoogle Scholar
- Curole JP, Kocher TD: Mitogenomics: digging deeper with complete mitochondrial genomes. Trends Ecol Evol. 1999, 14 (10): 394-398. 10.1016/S0169-5347(99)01660-2.View ArticleGoogle Scholar
- Telford MJ, Herniou EA, Russell RB, Littlewood DT: Changes in mitochondrial genetic codes as phylogenetic characters: two examples from the flatworms. Proc Natl Acad Sci USA. 2000, 97 (21): 11359-11364. 10.1073/pnas.97.21.11359.PubMed CentralView ArticleGoogle Scholar
- Waugh J, Huynen L, Millar C, Lambert D: DNA barcoding of animal species-response to DeSalle. Bioessays. 2008, 30 (1): 92-93. 10.1002/bies.20698.View ArticleGoogle Scholar
- Zink RM, Barrowclough GF: Mitochondrial DNA under siege in avian phylogeography. Mol Ecol. 2008, 17 (9): 2107-2121. 10.1111/j.1365-294X.2008.03737.x.View ArticleGoogle Scholar
- Saccone C, De Giorgi C, Gissi C, Pesole G, Reyes A: Evolutionary genomics in Metazoa: the mitochondrial DNA as a model system. Gene. 1999, 238 (1): 195-209. 10.1016/S0378-1119(99)00270-X.View ArticleGoogle Scholar
- Saccone C, Gissi C, Reyes A, Larizza A, Sbisa E, Pesole G: Mitochondrial DNA in metazoa: degree of freedom in a frozen event. Gene. 2002, 286 (1): 3-12. 10.1016/S0378-1119(01)00807-1.View ArticleGoogle Scholar
- Cameron SL, Dowton M, Castro LR, Ruberu K, Whiting MF, Austin AD, Diement K, Stevens J: Mitochondrial genome organization and phylogeny of two vespid wasps. Genome. 2008, 51 (10): 800-808. 10.1139/G08-066.View ArticleGoogle Scholar
- Waeschenbach A, Telford MJ, Porter JS, Littlewood DTJ: The complete mitochondrial genome of Flustrellidra hispida and the phylogenetic position of Bryozoa among the Metazoa. Mol Phylogenet Evol. 2006, 40: 195-207. 10.1016/j.ympev.2006.03.007.View ArticleGoogle Scholar
- Iannelli F, Griggio F, Pesole G, Gissi C: The mitochondrial genome of Phallusia mammillata and Phallusia fumigata (Tunicata, Ascidiacea): high genome plasticity at intra-genus level. BMC Evol Biol. 2007, 7: 155-10.1186/1471-2148-7-155.PubMed CentralView ArticleGoogle Scholar
- Gissi C, Pesole G: Transcript mapping and genome annotation of ascidian mtDNA using EST data. Genome Res. 2003, 13: 2203-2212. 10.1101/gr.1227803.PubMed CentralView ArticleGoogle Scholar
- Durrheim GA, Corfield VA, Harley EH, Ricketts MH: Nucleotide sequence of cytochrome oxidase (subunit III) from the mitochondrion of the tunicate Pyura stolonifera : evidence that AGR encodes glycine. Nucl Acids Res. 1993, 21 (15): 3587-3588. 10.1093/nar/21.15.3587.PubMed CentralView ArticleGoogle Scholar
- Xia X, Huang H, Carullo M, Betran E, Moriyama EN: Conflict between translation initiation and elongation in vertebrate mitochondrial genomes. PLoS One. 2007, 2 (2): e227-10.1371/journal.pone.0000227.PubMed CentralView ArticleGoogle Scholar
- Boore JL, Daehler LL, Brown WM: Complete sequence, gene arrangement, and genetic code of mitochondrial DNA of the cephalochordate Branchiostoma floridae (Amphioxus). Mol Biol Evol. 1999, 16 (3): 410-418.View ArticleGoogle Scholar
- Gissi C, Iannelli F, Pesole G: Complete mtDNA of Ciona intestinalis reveals extensive gene rearrangement and the presence of an atp8 and an extra trnM gene in ascidians. J Mol Evol. 2004, 58 (4): 376-389. 10.1007/s00239-003-2559-6.View ArticleGoogle Scholar
- Iannelli F, Pesole G, Sordino P, Gissi C: Mitogenomics reveals two cryptic species in Ciona intestinalis. Trends Genet. 2007, 23 (9): 419-422. 10.1016/j.tig.2007.07.001.View ArticleGoogle Scholar
- Yokobori S, Oshima T, Wada H: Complete nucleotide sequence of the mitochondrial genome of Doliolum nationalis with implications for evolution of urochordates. Mol Phylogenet Evol. 2005, 34 (2): 273-283. 10.1016/j.ympev.2004.10.002.View ArticleGoogle Scholar
- Yokobori S, Ueda T, Feldmaier-Fuchs G, Pääbo S, Ueshima R, Kondow A, Nishikawa K, Watanabe K: Complete DNA sequence of the mitochondrial genome of the ascidian Halocynthia roretzi (Chordata, Urochordata). Genetics. 1999, 153 (4): 1851-1862.PubMed CentralGoogle Scholar
- Yokobori S, Watanabe Y, Oshima T: Mitochondrial genome of Ciona savignyi (Urochordata, Ascidiacea, Enterogona): comparison of gene arrangement and tRNA genes with Halocynthia roretzi mitochondrial genome. J Mol Evol. 2003, 57 (5): 574-587. 10.1007/s00239-003-2511-9.View ArticleGoogle Scholar
- Castilla JC, Guinez R, Caro AU, Ortiz V: Invasion of a rocky intertidal shore by the tunicate Pyura praeputialis in the Bay of Antofagasta, Chile. Proc Natl Acad Sci USA. 2004, 101 (23): 8517-8524. 10.1073/pnas.0401921101.PubMed CentralView ArticleGoogle Scholar
- Locke A, Hanson JM, Ellis KM, Thompson J, Rochette R: Invasion of the southern Gulf of St. Lawrence by the clubbed tunicate (Styela clava Herdman): Potential mechanisms for invasions of Prince Edward Island estuaries. J Exp Mar Biol Ecol. 2007, 342 (1): 69-77. 10.1016/j.jembe.2006.10.016.View ArticleGoogle Scholar
- Ramsay A, Davidson J, Landry T, Arsenault G: Process of invasiveness among exotic tunicates in Prince Edward Island, Canada. Biological Invasions. 2008, 10 (8): 1311-1316. 10.1007/s10530-007-9205-y.View ArticleGoogle Scholar
- Therriault TW, Herborg LM: Predicting the potential distribution of the vase tunicate Ciona intestinalis in Canadian waters: informing a risk assessment. ICES J Mar Sci. 2008, 65 (5): 788-794. 10.1093/icesjms/fsn054.View ArticleGoogle Scholar
- Cañestro C, Bassham S, Postlethwait JH: Seeing chordate evolution through the Ciona genome sequence. Genome Biol. 2003, 4 (3): 208-10.1186/gb-2003-4-3-208.PubMed CentralView ArticleGoogle Scholar
- Cañestro C, Yokoi H, Postlethwait JH: Evolutionary developmental biology and genomics. Nat Rev Genet. 2007, 8 (12): 932-942. 10.1038/nrg2226.View ArticleGoogle Scholar
- Davidson B: Ciona intestinalis as a model for cardiac development. Semin Cell Dev Biol. 2007, 18 (1): 16-26. 10.1016/j.semcdb.2006.12.007.PubMed CentralView ArticleGoogle Scholar
- Manni L, Zaniolo G, Cima F, Burighel P, Ballarin L: Botryllus schlosseri : a model ascidian for the study of asexual reproduction. Dev Dyn. 2007, 236 (2): 335-352. 10.1002/dvdy.21037.View ArticleGoogle Scholar
- Godefroy N, Hoa C, Tsokanos F, Olofsson L, Le Goff E, Douzery EJP, Baghdiguian S, Martinand-Mari C: Identification of autophagy genes in Ciona intestinalis : A new experimental model to study autophagy mechanism. Autophagy. 2009, 5 (6): 805-815.View ArticleGoogle Scholar
- Tsagkogeorga G, Turon X, Hopcroft RR, Tilak M-K, Feldstein T, Shenkar N, Loya Y, Huchon D, Douzery EJP, Delsuc F: An updated 18S rRNA phylogeny of tunicates based on mixture and secondary structure models. BMC Evolutionary Biology. 2009, 9: 187-10.1186/1471-2148-9-187.PubMed CentralView ArticleGoogle Scholar
- Zeng L, Jacobs M, Swalla B: Coloniality and sociality has evolved once in Stolidobranch ascidians. Integr Comp Biol. 2006, 46: 255-268. 10.1093/icb/icj035.View ArticleGoogle Scholar
- Zeng L, Swalla B: Molecular phylogeny of the protochordates: chordate evolution. Can J Zool. 2005, 83: 24-33. 10.1139/z05-010.View ArticleGoogle Scholar
- Bourlat SJ, Juliusdottir T, Lowe CJ, Freeman R, Aronowicz J, Kirschner M, Lander ES, Thorndyke M, Nakano H, Kohn AB, Heyland A, Moroz LL, Copley RR, Telford MJ: Deuterostome phylogeny reveals monophyletic chordates and the new phylum Xenoturbellida. Nature. 2006, 444 (7115): 85-88. 10.1038/nature05241.View ArticleGoogle Scholar
- Bourlat SJ, Nielsen C, Economou AD, Telford MJ: Testing the new animal phylogeny: a phylum level molecular analysis of the animal kingdom. Mol Phylogenet Evol. 2008, 49 (1): 23-31. 10.1016/j.ympev.2008.07.008.View ArticleGoogle Scholar
- Bourlat SJ, Rota-Stabelli O, Lanfear R, Telford MJ: The mitochondrial genome structure of Xenoturbella bocki (phylum Xenoturbellida) is ancestral within the deuterostomes. BMC Evol Biol. 2009, 9: 107-10.1186/1471-2148-9-107.PubMed CentralView ArticleGoogle Scholar
- Castoe TA, Jiang ZJ, Gu W, Wang ZO, Pollock DD: Adaptive evolution and functional redesign of core metabolic proteins in snakes. PLoS ONE. 2008, 3 (5): e2201-10.1371/journal.pone.0002201.PubMed CentralView ArticleGoogle Scholar
- Zhong J, Zhang J, Mukwaya E, Wang Y: Revaluation of deuterostome phylogeny and evolutionary relationships among chordate subphyla using mitogenome data. J Genet Genomics. 2009, 36 (3): 151-160. 10.1016/S1673-8527(08)60102-2.View ArticleGoogle Scholar
- Delsuc F, Brinkmann H, Chourrout D, Philippe H: Tunicates and not cephalochordates are the closest living relatives of vertebrates. Nature. 2006, 439 (7079): 965-968. 10.1038/nature04336.View ArticleGoogle Scholar
- Delsuc F, Tsagkogeorga G, Lartillot N, Philippe H: Additional molecular support for the new chordate phylogeny. Genesis. 2008, 46 (11): 592-604. 10.1002/dvg.20450.View ArticleGoogle Scholar
- Putnam NH, Butts T, Ferrier DE, Furlong RF, Hellsten U, Kawashima T, Robinson-Rechavi M, Shoguchi E, Terry A, Yu JK, et al: The amphioxus genome and the evolution of the chordate karyotype. Nature. 2008, 453 (7198): 1064-1071. 10.1038/nature06967.View ArticleGoogle Scholar
- Dunn CW, Hejnol A, Matus DQ, Pang K, Browne WE, Smith SA, Seaver E, Rouse GW, Obst M, Edgecombe GD, Sørensen MV, Haddock SH, Schmidt-Rhaesa A, Okusu A, Kristensen RM, Wheeler WC, Martindale MQ, Giribet G: Broad phylogenomic sampling improves resolution of the animal tree of life. Nature. 2008, 452 (7188): 745-749. 10.1038/nature06614.View ArticleGoogle Scholar
- Shenkar N, Loya Y: The solitary ascidian Herdmania momus : native (Red Sea) versus non-indigenous (Mediterranean) populations. Biol Invasions. 2008, 10 (8): 1431-1439. 10.1007/s10530-008-9217-2.View ArticleGoogle Scholar
- Nohara M, Nishida M, Miya M, Nishikawa T: Evolution of the mitochondrial genome in Cephalochordata as inferred from complete nucleotide sequences from two Epigonichthys species. J Mol Evol. 2005, 60 (4): 526-537. 10.1007/s00239-004-0238-x.View ArticleGoogle Scholar
- Philippe H: Rodent monophyly: pitfalls of molecular phylogenies. J Mol Evol. 1997, 45: 712-715.Google Scholar
- Phillips MJ, McLenachan PA, Down C, Gibb GC, Penny D: Combined mitochondrial and nuclear DNA sequences resolve the interrelations of the major Australasian marsupial radiations. Syst Biol. 2006, 55 (1): 122-137. 10.1080/10635150500481614.View ArticleGoogle Scholar
- Telford MJ: Animal phylogeny. Curr Biol. 2006, 16 (23): R981-985. 10.1016/j.cub.2006.10.048.View ArticleGoogle Scholar
- Le SQ, Gascuel O, Lartillot N: Empirical profile mixture models for phylogenetic reconstruction. Bioinformatics. 2008, 24 (20): 2317-2323. 10.1093/bioinformatics/btn445.View ArticleGoogle Scholar
- Lartillot N, Brinkmann H, Philippe H: Suppression of long-branch attraction artefacts in the animal phylogeny using a site-heterogeneous model. BMC Evol Biol. 2007, 7 (Suppl 1): S4-10.1186/1471-2148-7-S1-S4.PubMed CentralView ArticleGoogle Scholar
- Lartillot N, Philippe H: A Bayesian mixture model for across-site heterogeneities in the amino-acid replacement process. Mol Biol Evol. 2004, 21 (6): 1095-1109. 10.1093/molbev/msh112.View ArticleGoogle Scholar
- Blanquart S, Lartillot N: A site- and time-heterogeneous model of amino acid replacement. Mol Biol Evol. 2008, 25 (5): 842-858. 10.1093/molbev/msn018.View ArticleGoogle Scholar
- Philippe H, Lartillot N, Brinkmann H: Multigene analyses of bilaterian animals corroborate the monophyly of Ecdysozoa, Lophotrochozoa, and Protostomia. Mol Biol Evol. 2005, 22 (5): 1246-1253. 10.1093/molbev/msi111.View ArticleGoogle Scholar
- Lartillot N, Philippe H: Improvement of molecular phylogenetic inference and the phylogeny of Bilateria. Philos Trans R Soc Lond B Biol Sci. 2008, 363 (1496): 1463-1472. 10.1098/rstb.2007.2236.PubMed CentralView ArticleGoogle Scholar
- Bromham LD, Degnan BM: Hemichordates and deuterostome evolution: robust molecular phylogenetic support for a hemichordate + echinoderm clade. Evol Dev. 1999, 1 (3): 166-171. 10.1046/j.1525-142x.1999.99026.x.View ArticleGoogle Scholar
- Wada H: Evolutionary history of free-swimming and sessile lifestyles in urochordates as deduced from 18S rDNA molecular phylogeny. Mol Biol Evol. 1998, 15 (9): 1189-1194.View ArticleGoogle Scholar
- Jiang ZJ, Castoe TA, Austin CC, Burbrink FT, Herron MD, McGuire JA, Parkinson CL, Pollock DD: Comparative mitochondrial genomics of snakes: extraordinary substitution rate dynamics and functionality of the duplicate control region. BMC Evol Biol. 2007, 7: 123-10.1186/1471-2148-7-123.PubMed CentralView ArticleGoogle Scholar
- Castoe TA, de Koning AP, Kim HM, Gu W, Noonan BP, Naylor G, Jiang ZJ, Parkinson CL, Pollock DD: Evidence for an ancient adaptive episode of convergent molecular evolution. Proc Natl Acad Sci USA. 2009Google Scholar
- Sambrook J, Fritsch EF, Maniatis T: Molecular cloning: a laboratory manual. 1989, New York: Cold Spring Harbor Laboratory Press, 2Google Scholar
- Batzoglou S, Jaffe DB, Stanley K, Butler J, Gnerre S, Mauceli E, Berger B, Mesirov JP, Lander ES: ARACHNE: a whole-genome shotgun assembler. Genome Res. 2002, 12 (1): 177-189. 10.1101/gr.208902.PubMed CentralView ArticleGoogle Scholar
- The NCBI ORF-Finder tool. [http://0-www.ncbi.nlm.nih.gov.brum.beds.ac.uk/projects/gorf/]
- Wyman SK, Jansen RK, Boore JL: Automatic annotation of organellar genomes with DOGMA. Bioinformatics. 2004, 20 (17): 3252-3255. 10.1093/bioinformatics/bth352.View ArticleGoogle Scholar
- Wolstenholme DR: Animal mitochondrial DNA: structure and evolution. Int Rev Cytol. 1992, 141: 173-216. 10.1016/S0074-7696(08)62066-5.View ArticleGoogle Scholar
- Lowe TM, Eddy SR: tRNAscan-SE: A program for improved detection of transfer RNA genes in genomic sequence. Nucleic Acids Res. 1997, 25 (5): 955-964. 10.1093/nar/25.5.955.PubMed CentralView ArticleGoogle Scholar
- Zuker M: Mfold web server for nucleic acid folding and hybridization prediction. Nucleic Acids Res. 2003, 31 (13): 3406-3415. 10.1093/nar/gkg595.PubMed CentralView ArticleGoogle Scholar
- Jameson D, Gibson AP, Hudelot C, Higgs PG: OGRe: a relational database for comparative analysis of mitochondrial genomes. Nucleic Acids Res. 2003, 31 (1): 202-206. 10.1093/nar/gkg077.PubMed CentralView ArticleGoogle Scholar
- Katoh K, Kuma K, Toh H, Miyata T: MAFFT version 5: improvement in accuracy of multiple sequence alignment. Nucleic Acids Res. 2005, 33 (2): 511-518. 10.1093/nar/gki198.PubMed CentralView ArticleGoogle Scholar
- Castresana J: Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis. Mol Biol Evol. 2000, 17 (4): 540-552.View ArticleGoogle Scholar
- Guindon S, Gascuel O: A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. Syst Biol. 2003, 52 (5): 696-704. 10.1080/10635150390235520.View ArticleGoogle Scholar
- Yang Z: Among-site rate variation and its impact on phylogenetic analyses. Trends Ecol Evol. 1996, 11 (9): 367-372. 10.1016/0169-5347(96)10041-0.View ArticleGoogle Scholar
- Felsenstein J: PHYLIP (PHYLogeny Inference Package). Version 3.66 edn. 2001, Seattle: Department of Genome Sciences. University of WashingtonGoogle Scholar
- Jobb G, von Haeseler A, Strimmer K: TREEFINDER: a powerful graphical analysis environment for molecular phylogenetics. BMC Evol Biol. 2004, 4 (1): 18-10.1186/1471-2148-4-18.PubMed CentralView ArticleGoogle Scholar
- Lartillot N, Lepage T, Blanquart S: PhyloBayes 3. A Bayesian software package for phylogenetic reconstruction and molecular dating. Bioinformatics. 2009, 25: 2286-2288. 10.1093/bioinformatics/btp368.View ArticleGoogle Scholar
- Blanquart S, Lartillot N: A Bayesian compound stochastic process for modeling nonstationary and nonhomogeneous sequence evolution. Mol Biol Evol. 2006, 23 (11): 2058-2071. 10.1093/molbev/msl091.View ArticleGoogle Scholar
- Swofford DL: PAUP*: Phylogenetic Analysis Using Parsimony and other methods, 4b10 edn. 2002, Sinauer, Sunderland, MAGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.