Skip to main content

Advertisement

We’d like to understand how you use our websites in order to improve them. Register your interest.
Download PDF

A phylogenomic framework and timescale for comparative studies of tunicates

BMC Biology volume 16, Article number: 39 (2018) Cite this article

Abstract

Background

Tunicates are the closest relatives of vertebrates and are widely used as models to study the evolutionary developmental biology of chordates. Their phylogeny, however, remains poorly understood, and to date, only the 18S rRNA nuclear gene and mitogenomes have been used to delineate the major groups of tunicates. To resolve their evolutionary relationships and provide a first estimate of their divergence times, we used a transcriptomic approach to build a phylogenomic dataset including all major tunicate lineages, consisting of 258 evolutionarily conserved orthologous genes from representative species.

Results

Phylogenetic analyses using site-heterogeneous CAT mixture models of amino acid sequence evolution resulted in a strongly supported tree topology resolving the relationships among four major tunicate clades: (1) Appendicularia, (2) Thaliacea + Phlebobranchia + Aplousobranchia, (3) Molgulidae, and (4) Styelidae + Pyuridae. Notably, the morphologically derived Thaliacea are confirmed as the sister group of the clade uniting Phlebobranchia + Aplousobranchia within which the precise position of the model ascidian genus Ciona remains uncertain. Relaxed molecular clock analyses accommodating the accelerated evolutionary rate of tunicates reveal ancient diversification (~ 450–350 million years ago) among the major groups and allow one to compare their evolutionary age with respect to the major vertebrate model lineages.

Conclusions

Our study represents the most comprehensive phylogenomic dataset for the main tunicate lineages. It offers a reference phylogenetic framework and first tentative timescale for tunicates, allowing a direct comparison with vertebrate model species in comparative genomics and evolutionary developmental biology studies.

Background

Large-scale phylogenetic analyses of tunicate genomic data from a handful of model species have identified this marine chordate group as the closest relative of vertebrates [1,2,3,4,5]. This discovery has had profound implications for comparative genomics and evolutionary developmental biology (evo-devo) studies aimed at understanding the origins of chordates and vertebrates [6,7,8]. Indeed, the new chordate phylogeny implies that the tunicate body plan is evolutionarily derived and has become secondarily simplified from that of more complex chordate ancestors [2, 3].

The key phylogenetic position of tunicates within chordates has prompted the selection of model species such as Ciona robusta (formerly Ciona intestinalis type A [9]), for which a full genome has been sequenced early in the history of comparative genomics to provide insight into vertebrate-specific whole genome duplications [10]. Since then, genome sequences have been assembled for additional species that are widely used as models in comparative genomics and evo-devo [11] including Ciona savignyi [12], Oikopleura dioica [5], Botryllus schlosseri [13], Molgula occidentalis, M. occulta, and M. occulata [14], Phallusia mammillata [15], and Halocynthia roretzi [15]. The available genomic data have notably revealed a stunning contrast in the evolutionary rate of nuclear protein-coding genes between tunicates and vertebrates [3, 16]. This accelerated evolution of tunicate genes is also coupled with extensive structural rearrangements observed in their genomes [5, 17, 18]. This contrast is even more pronounced for mitochondrial genomes, which are particularly fast evolving and highly rearranged in tunicates with respect to other deuterostomes, in which they are widely conserved [5, 19, 20]. The reasons behind the rapid rate of genomic evolution in tunicates remain unclear [16, 21, 22] and contrast with the unusual conservation level of embryonic morphologies between all ascidian species studied so far [7].

Despite renewed interest in tunicate evolution, phylogenetic relationships among the major tunicate lineages remain uncertain. Previous molecular phylogenetic studies relying on 18S rRNA [23,24,25,26] and mitogenomes [20, 27, 28] have proposed first delineations of major tunicate clades, revoking the traditional nineteenth century classification into the three classes Appendicularia (larvaceans), Thaliacea (salps, doliolids, and pyrosomes), and Ascidiacea (phlebobranchs, aplousobranchs, and stolidobranchs). Indeed, these studies found unanimous support for the paraphyly of Ascidiacea (ascidians) owing to the inclusion of thaliaceans in a clade also containing two main ascidian lineages (phlebobranchs and aplousobranchs) to the exclusion of stolidobranch ascidians (molgulids, pyurids, and styelids). Nevertheless, the resolving power of these standard markers — nuclear ribosomal RNA and mitochondrial protein-coding genes — appeared to be limited regarding the relationships among the three newly proposed main clades: (1) Appendicularia, (2) Stolidobranchia, and (3) Phlebobranchia + Thaliacea + Aplousobranchia. Notably, the relationships within the latter group were left unresolved, with the position of thaliaceans relative to phlebobranchs and aplousobranchs still being debated [25, 27].

The phylogenetic position of thaliaceans is key for understanding the evolution of developmental modes within tunicates [29]. Compared to their closest relatives, which are mostly solitary and sessile, the three groups of thaliaceans (salps, doliolids, and pyrosomes) are pelagic with complex life cycles including solitary and colonial phases. Their unique lifestyle also seems to be associated with spectacular differences in their embryology, such as the loss of a well-developed notochord in the larva of most thaliaceans, with the exception of only a few doliolid species [29]. Based on our current understanding of tunicate evolution, thaliaceans may have evolved from a sessile ascidian-like ancestor and therefore can serve as a model to understand how the transition from a benthic to a pelagic lifestyle has led to drastic modifications in the morphology, embryology, and life cycle of these tunicates [29]. Coloniality is another remarkable feature of the thaliaceans, which shows some similarities with the coloniality in ascidians, even though this trait probably evolved independently in the two groups [29]. It is noteworthy that doliolids have polymorphic colonies [30], a trait that is absent in colonial ascidians. A reliable phylogeny positioning thaliaceans with regard to colonial ascidians is thus necessary to understand the evolution of these unique features.

Outstanding questions in chordate evolution include the identification of the determinants of the rapid rate of genome evolution in tunicates and the emergence of vertebrates [11, 31]. A prerequisite to addressing these issues is to reconstruct a reliable phylogenetic framework and timescale to guide future comparative evolutionary genomic and evolutionary studies of chordate development. Moreover, given that the fossil record of tunicates is deceptively scarce and controversial [32,33,34], a molecular timescale for chordates would allow one to compare tunicate evolution to that of the well-calibrated vertebrates [35] for the first time. A phylogenetic and timing framework is notably critical for the identification and interpretation of both conserved and divergent developmental features of tunicates compared to model vertebrate species in the context of their fast rate of genomic evolution [11].

Here, we use new transcriptomic data obtained through high-throughput sequencing technologies (Roche 454 and Illumina HiSeq) to build the first tunicate phylogenomic dataset including all major tunicate groups. This dataset consists of 258 orthologous nuclear genes for 63 taxa including representative deuterostome species and all major chordate lineages. Using phylogenetic analyses based on the best-fitting site-heterogeneous CAT mixture model of amino acid sequence evolution, we inferred well-resolved phylogenetic relationships for the major clades of tunicates. Our molecular dating analyses based on models of clock relaxation accounting for variation in lineage-specific evolutionary rates provide a first tentative timescale for the emergence of the main tunicate clades, allowing a direct comparison with vertebrate model systems.

Methods

Transcriptome data collection

Live tunicate specimens were ordered from Gulf Specimen Marine Laboratories, Inc. (Panacea, FL, USA) and the Roscoff Biological Station (Roscoff, France) services and collected in Villefranche-sur-Mer (France) and Blanes (Spain). One single run of Roche 454 GS-FLX Titanium was conducted at GATC Biotech (Konstanz, Germany) on multiplexed total RNA libraries that were constructed for Clavelina lepadiformis, Cystodytes dellechiajei, Bostrichobranchus pilularis, Molgula manhattensis, Molgula occidentalis, Phallusia mammillata, Dendrodoa grossularia, Polyandrocarpa anguinea, and Styela plicata. Complementary RNA-seq data were acquired with paired-end 100-nt Illumina reads at Beijing Genome Institute (Shenzhen, China) for the thaliaceans Salpa fusiformis (mix of two blastozooids) and Doliolum nationalis (mix of 15 phorozooids), and with single-end 100-nt Illumina reads at GATC Biotech (Konstanz, Germany) for Clavelina lepadiformis and Cystodytes dellechiajei (mix of several individuals) [36]. Previously obtained 454 transcriptomic data for Microcosmus squamiger [16] were also considered. De novo assemblies were conducted with Trinity [37] for 454 reads and ABySS [38] for Illumina reads using the programs’ default parameters. For both kinds of libraries, we confirmed the sample taxonomic identifications by assembling the mitochondrial cytochrome c oxidase subunit 1 (CO1) and nuclear 18S rRNA barcoding genes and reconstructing maximum likelihood trees with available comparative data. Additional tunicate sequences were collected in public databases from various sequencing projects: Botryllus schlosseri, Halocynthia roretzi, and Diplosoma listerianum (expressed sequence tags (ESTs)), Molgula tectiformis (complementary DNAs), and Ciona robusta, Ciona savignyi, and Oikopleura dioica (genomes). Detailed information on biological specimens, basic statistics, and accession numbers of newly sequenced transcriptomes can be found in Additional file 1: Table S1.

Phylogenomic dataset assembly

We built upon a previous phylogenomic dataset [39] to select a curated set of 258 orthologous markers for deuterostomes. Alignments were complemented with sequences from the National Center for Biotechnology Information (NCBI) databases using a multiple best reciprocal hit approach implemented in the newly designed Forty-Two software [40]. Because 454 DNA sequence reads are characterized by sequencing errors typically disrupting the reading frame when translated into amino acids, alignments were verified by eye using the program ED from the MUST package [41]. Ambiguously aligned regions were excluded for each individual protein using Gblocks with medium default parameters [42] with a few subsequent manual refinements using NET from the MUST package to relax the fact that this automated approach is sometimes too conservative. This manual refinement step restored only 418 amino acid sites (i.e. 0.6% of the total alignment length). Potential environmental contaminations and cross-contaminations between our samples were also dealt with at the alignment stage by performing Basic Local Alignment Search Tool (BLAST) searches of each sequence against a taxon-rich reference database maintained for each curated gene alignment and were further sought by a visual examination of each individual gene phylogeny.

The concatenation of the resulting 258 amino acid alignments was constructed with SCaFoS [43] by defining 63 deuterostomian operational taxonomic units (OTUs) representing all major lineages. The taxon sampling included 18 tunicates, 34 vertebrates, and one cephalochordate, with seven echinoderms, two hemichordates, and one xenoturbellid as more distant outgroups. When several sequences were available for a given OTU, the slowest evolving one was selected by SCaFoS, according to maximum likelihood distances computed by TREE-PUZZLE [44] under a WAG+F model. The percentage of missing data per taxon was reduced by creating some chimerical sequences from closely related species (i.e. Eptatretus burgeri/ Myxine glutinosa, Petromyzon marinus/Lethenteron japonicum, Callorhinchus milii/ C. callorynchus, Latimeria menadoensis/L. chalumnae, Rana chensinensis/ R. catesbeiana, Alligator sinensis/ A. mississippiensis, Chrysemys picta/ Emys orbicularis/ Trachemys scripta, Patiria miniata/ P. pectinifera/ Solaster stimpsonii, Apostichopus japonicus/ Parastichopus parvimensis, Ophionotus victoriae/ Amphiura filiformis) and by retaining only proteins with at most 15 missing OTUs. The tunicate Microcosmus squamiger was excluded at this stage due to a high percentage of missing data resulting from the low number of contigs obtained in the assembly. The final alignment comprised 258 proteins and 63 taxa for 66,593 unambiguously aligned amino acid sites with 20% missing amino acid data.

Phylogenetic analyses

Bayesian cross-validation [45] implemented in PhyloBayes 3.3f [46] was used to compare the fit of site-homogeneous (LG and GTR) and site-heterogeneous (CAT-F81 and CAT-GTR) models coupled with a gamma distribution (Γ4) of site-rate heterogeneity. Ten replicates were considered, each one consisting of a random subsample of 10,000 sites for training the model and 2000 sites for computing the cross-validation likelihood score. Under site-homogeneous LG + Γ4 and GTR + Γ4 models, 1100 sampling cycles were run and a burn-in of 100 samples was used, and under site-heterogeneous models CAT-F81 + Γ4 and CAT-GTR + Γ4, 3100 sampling cycles were run and the first 2100 samples were discarded as burn-in.

Bayesian phylogenetic reconstruction under the best-fitting CAT-GTR + Γ4 mixture model [47] was conducted using PhyloBayes_MPI 1.5a [48]. Two independent Markov chain Monte Carlo (MCMC) simulations starting from a randomly generated tree were run for 6000 cycles with trees and associated model parameters being sampled every cycle. The initial 1000 trees sampled in each MCMC run were discarded as burn-in after checking for convergence in both likelihood and model parameters, as well as in clade posterior probabilities using bpcomp (max_diff < 0.3). The 50% majority-rule Bayesian consensus tree and the associated posterior probabilities (PPs) were then computed from the remaining combined 10,000 (2 × 5000) trees using bpcomp.

We further assessed the robustness of our phylogenomic inference by applying a gene jackknife resampling procedure [3]. A hundred jackknife replicates constituting 130 alignments drawn randomly out of the total 258 protein alignments were generated. The 100 resulting jackknife supermatrices were then analysed using PhyloBayes_MPI under the second best-fitting CAT-F81 + Γ4 model instead of the best-fitting CAT-GTR + Γ4 and for 2000 sampling cycles in order to reduce the computational burden. After removing the first 200 sampled trees of each chain as the burn-in, a majority-rule consensus tree was obtained for each replicate using the 1800 trees sampled from the posterior distribution. A consensus tree was then obtained from the 100 jackknife-resampled consensus trees. The support values displayed by this Bayesian consensus tree are thus gene jackknife support (JS) percentages. High values indicate nodes that have high posterior probability support in most jackknife replicates and are thus robust to gene sampling. We verified convergence of MCMCs in each jackknife replicate by checking that varying the burn-in value did not affect the JS percentages obtained in the final consensus.

Molecular dating

Molecular dating analyses were performed in a Bayesian relaxed molecular clock framework using PhyloBayes 3.3f [46]. In all dating calculations, the tree topology was fixed to the majority-rule consensus tree inferred in previous Bayesian analyses (Fig. 1). Dating analyses were conducted using the best-fitting site-heterogeneous CAT-GTR + Γ4 mixture model and a relaxed clock model with a birth-death prior on divergence times combined with soft fossil calibrations following Lartillot et al. [46]. Given the lack of trustable fossils within tunicates, we used 12 calibration intervals defined within vertebrates [49, 50] and one within echinoderms [51]: (1) Chordata (Max. Age 581 Mya, Min. Age 519 Mya), (2) Olfactores (Max 581 Mya, Min 519), (3) Vertebrata (Max 581 Mya, Min 461), (4) Gnathostomata (Max 463 Mya, Min 422), (5) Osteichthyes (Max 422 Mya, Min 416), (6) Tetrapoda (Max 350 Mya, Min 330), (7) Amniota (Max 330 Mya, Min 312), (8) Diapsida (Max 300 Mya, Min 256), (9) Batrachia (Max 299 Mya, Min 200), (10) Clupeocephala (Max 165 Mya, Min 150), (11) Mammalia (Max 191 Mya, Min 163), (12) Theria (Max 171 Mya, Min 124), and (13) Echinoidea (Min 255 Mya). The prior on the root of the tree (Deuterostomia) was set to an exponential distribution of mean 540 Mya.

Fig. 1
figure1

Phylogenetic relationships of 63 chordates highlighting the major tunicate groups inferred from 66,593 amino acid sites of 258 proteins. The Bayesian consensus phylogram has been inferred by PhyloBayes_MPI under the CAT-GTR + Γ4 mixture model. Values at nodes indicate Bayesian posterior probabilities (PPCAT-GTR) obtained under CAT-GTR + Γ4, and jackknife support (JS) percentages, respectively. Circles at nodes pinpoint branches with maximal support from both methods. Species with newly obtained data are indicated in bold. The branch leading to the fast-evolving Oikopleura dioica has been halved for graphical purposes

Full size image

In order to select the best-fitting clock model, we compared the autocorrelated log-normal (LN) relaxed clock model [52] with the uncorrelated gamma (UGAM) model [53] and a strict molecular clock (CL) model. These three clock models were compared against each other using the same prior settings (see above) in a cross-validation procedure as implemented in PhyloBayes following Lepage et al. [54]. However, to reduce the computational burden, the CAT-F81 + Γ4 mixture model was used instead of CAT-GTR + Γ4. The cross-validation tests were performed by dividing the original alignment in two subsets with 90% of sites for the learning set (59,934 sites) and 10% of sites for the test set (6659 sites). The overall procedure was repeated over 10 random splits for which an MCMC simulation was run on the learning set for a total 4000 cycles sampling posterior rates and dates every cycle. The first 3000 samples of each MCMC were excluded as the burn-in for calculating the cross-validation scores averaged across the 10 replicates.

The final dating calculations were conducted under both LN and UGAM relaxed clock models and the CAT-GTR + Γ4 mixture model of sequence evolution by running MCMCs for a total 25,000 cycles sampling posterior rates and dates every 10 cycles. The first 500 samples of each MCMC were excluded as the burn-in after checking for convergence in both likelihood and model parameters using readdiv. Posterior estimates of divergence dates and associated 95% credibility intervals were then computed from the remaining 2000 samples of each MCMC using readdiv. Additional dating calculations using the same sampling scheme were also conducted under the LN relaxed clock model but using the less computationally intensive CAT-F81 + Γ4 mixture model.

Results and discussion

A reference phylogenetic framework for model tunicates

The evolutionary relationships of tunicates have long been a matter of debate, mainly because tunicates are characterized by an overall accelerated rate of evolution in their nuclear and mitochondrial genomes compared to other deuterostome species. Moreover, the large lineage-specific variation in evolutionary rates among tunicates [16] could result in long-branch attraction (LBA) artefacts, which hamper the reliable reconstruction of their phylogenetic relationships [55,56,57]. Another contributing factor to our limited understanding of tunicate evolution is the uneven availability of genome data across different tunicate lineages. To address these limitations, we used: (1) a wider taxon sampling encompassing all major tunicate lineages including two divergent thaliaceans, (2) numerous nuclear genes to reduce stochastic error, and (3) powerful site-heterogeneous models that generally offer the best fit to phylogenomic data and have the advantage of being least sensitive to LBA and other potential phylogenetic artefacts [39, 58, 59]. Accordingly, the results of our Bayesian cross-validation tests showed that the CAT-GTR + Γ4 mixture model offered the best statistical fit to the data (ΔlnL = 1506 ± 98 compared to LG + Γ4), followed by the CAT-F81 + Γ4 mixture model (ΔlnL = 817 ± 112 compared to LG + Γ4) and the GTR + Γ4 model (ΔlnL = 266 ± 41 compared to LG + Γ4).

The majority-rule consensus tree obtained using Bayesian phylogenetic reconstruction under the best-fitting CAT-GTR + Γ4 site-heterogeneous mixture model is thus presented in Fig. 1. This well-supported phylogenetic tree has been rooted between Xenambulacraria (Xenoturbellida + Ambulacraria) and Chordata following the results of Philippe et al. [39] showing that Xenacoelomorpha (acoelomorphs + xenoturbellids) were related to Ambulacraria (hemichordates + echinoderms) within Deuterostomia. These results have been recently challenged by two studies claiming support for a more external position of Xenacoelomorpha as a sister group to Nephrozoa (Protostomia + Deuterostomia) [60, 61]. However, this newly proposed position is still debated, as it might be the result of an LBA artefact caused by the very long branches of acoelomorphs in phylogenomic trees [62, 63]. Hence, we have chosen to root our trees according to Philippe et al. [39], which in any case does not affect the phylogenetic relationships of chordates.

The inferred topology unambiguously recovered the monophyly of chordates (PP = 1.0; JS = 100) and grouped the reciprocally monophyletic tunicates and vertebrates into Olfactores to the exclusion of cephalochordates (PP = 1.0; JS = 100) in accordance with the newly established chordate phylogeny [1, 3, 4]. Within tunicates, the appendicularian Oikopleura dioica was the sister group of all other included taxa (PP = 1.0; JS = 100). Within the latter, there was a well-supported split (PP = 1.0; JS = 100) between Stolidobranchia on one side, and Phlebobranchia, Aplousobranchia, and Thaliacea on the other side. The monophyletic Stolidobranchia included two main clades, the first corresponding to the family Molgulidae (PP = 1.0; JS = 100), and the second grouping the families Pyuridae and Styelidae (PP = 1.0; JS = 100). Within molgulids, Bostrichobranchus pilularis was the sister group of the three species within the genus Molgula (PP = 1.0; JS = 100), while M. occidentalis was the sister group of M. manhattensis + M. tectiformis (PP = 1.0; JS = 100). Lastly, the four styelids Styela plicata, Botryllus schlosseri, Polyandrocarpa anguinea, and Dendrodoa grossularia constituted a monophyletic group (PP = 1.0; JS = 100) with respect to the single species here representing pyurids (Halocynthia roretzi). Within styelids, S. plicata diverged first (PP = 1.0; JS = 98) followed by B. schlosseri as the sister group of P. anguinea + D. grossularia (PP = 1.0; JS = 100). On the other side of the tree, Thaliacea branched with maximum statistical support (PP = 1.0; JS = 100) as the sister group of the clade Phlebobranchia + Aplousobranchia. The traditional class-level taxon Ascidiacea — currently considered to embrace the orders Aplousobranchia, Phlebobranchia and Stolidobranchia [64] — therefore refers to a paraphyletic assemblage. An alternative classification scheme based on gonad position (not commonly used nowadays) recognized two orders within ascidians: Enterogona (corresponding to Phlebobranchia + Aplousobranchia) and Pleurogona (= Stolidobranchia) [30, 65]. These alternative order-level taxa are recovered as monophyletic in our analyses. The three aplousobranchs analysed here unambiguously formed a monophyletic clade (PP = 1.0; JS = 100) with Clavelina lepadiformis being the sister group of Diplosoma listerianum and Cystodytes dellechiajei (PP = 1.0; JS = 100). The phlebobranchs appeared as a paraphyletic group with the two Ciona species branching closer to the aplousobranchs than to the other phlebobranch species (Phallusia mammillata), although with no statistical support from the gene jackknife resampling analysis (PP = 100; JS = 42).

The results from this first phylogenomic study including all tunicate lineages were in line with recent studies [20, 25,26,27,28] demonstrating that ascidians (class Ascidiacea) form a paraphyletic group. Our results showed that phlebobranchs and aplousobranchs are undoubtedly closer to thaliaceans than to stolidobranchs (Fig. 1), and that a thorough taxonomic revision of the tunicate classes is necessary. It seems clear that the use of the Ascidiacea class should be abandoned in favour of more meaningful classification schemes. Even though the position of Thaliacea was not always statistically supported, it consistently appeared as the sister group of phlebobranchs + aplousobranchs in previous studies [20, 24,25,26, 28], except for a recent genome-scale study in which the positioning of Salpa thompsoni most likely suffered artefactual LBA attraction towards the fast-evolving appendicularians [66]. The robust phylogenetic position of thaliaceans found here indicates that they likely evolved from a sessile ancestor, and their study can provide valuable information on the morphological transformations associated with the transition to the pelagic lifestyle [29].

The monophyly of the clade uniting phlebobranchs and aplousobranchs has never been challenged, and thus we suggest to re-use the term Enterogona to define this group, as originally proposed by Perrier [65] and subsequently redefined by Garstang [67]. The close relationship between thaliaceans and enterogones has also been supported by all previous molecular studies, as well as by morphological observations. The gonad position and the shared paired ontogenetic rudiment of the atrial cavity and opening might constitute two of their anatomical synapomorphies [68]. Lastly, we also confirmed the previously reported monophyly of stolidobranchs (= Pleurogona), with molgulids being the sister group to styelids + pyurids.

Finally, our phylogenomic study casts new light on two recurring issues in tunicate phylogenetics. First, phlebobranchs have been repeatedly found to be paraphyletic, albeit usually with no statistical support [25,26,27,28, 69], and the phylogenetic affinities among its members remain unclear. Notably, the traditional position of Ciona as a phlebobranch ascidian was challenged by Kott [70], who placed the genus within aplousobranchs on the basis of morphological characters. More recently, Turon and López-Legentil [69] and Shenkar et al. [28] found that Ciona was closer to aplousobranchs than to other phlebobranchs using mitochondrial DNA. These results are in agreement with the tree topology obtained in the present study, although it was not statistically supported. The positioning of the model Ciona genus and the phylogenetic relationships of phlebobranchs need to be the focus of additional phylogenomic studies including a denser taxon sampling. Second, although the position of appendicularians as sister clade to all other tunicates was well supported here and in all previous tunicate phylogenomic studies [2, 3], the extremely long branch of Oikopleura dioica coupled with our current inability to completely alleviate a potential LBA artefact — even with complex site-heterogeneous mixture models (see [59]) — prevent us from considering this species phylogenetic position as conclusive. The long appendicularian branch should be subdivided with the inclusion of additional divergent species in future phylogenomic analyses to definitively settle this point.

Evolutionary rate variations and molecular clock models

As observed in previous phylogenomic studies of chordates [2,3,4], the Bayesian phylogram estimated under the best-fitting CAT-GTR + Γ4 mixture model revealed marked branch length heterogeneity (Fig. 1). The tunicate branch lengths not only were much longer than those of all the other deuterostome clades, but they also displayed strong variations within tunicates. From the ancestral node of Olfactores, the tunicate median branch length was of 1.53 amino acid substitutions per site compared to the vertebrate median branch length, which was 0.65. From the ancestral vertebrate node, the average of branch lengths is 0.35 ± 0.05 amino acid replacements per site. In contrast, from the ancestral node of tunicates — excluding the super fast-evolving Oikopleura dioica — the average of branch lengths was 0.69 ± 0.19. For the proteins combined here for phylogenomic purposes, tunicates (with the exception of O. dioica) displayed on average a twice-higher number of amino acid substitutions than vertebrates.

Such substitution rate variation among lineages — within tunicates, and between tunicates and other deuterostomes — needs to be accounted for in molecular dating analyses by using models of clock relaxation [52]. The selection of the clock model is often arbitrary and appears mostly dependent on the software choice, with an overwhelming majority of studies relying on the BEAST software [71] using an uncorrelated gamma (UGAM, also known as UCLN) model of clock relaxation. However, it has been shown that autocorrelated rate models, such as the autocorrelated LN model, often provide a better fit with phylogenomic data [54, 72, 73]. Consequently, we compared the fit of both the UGAM and LN models to the fit of a strict CL model for our dataset using cross-validation tests under the CAT-GTR + Γ4 model. As expected given the large lineage-specific rate variation, both relaxed clock models largely outperformed the strict clock model (UGAM vs. CL: ΔlnL = 4068 ± 125; LN vs. CL: ΔlnL = 4057 ± 118). Among relaxed clock models, UGAM and LN were statistically equivalent in offering a very similar fit to our data (UGAM vs. LN: ΔlnL = 11 ± 38).

The use of a relaxed clock model allowed us to perform evolutionary rate comparisons in terms of number of substitutions per site per million years for the 63 terminal taxa considered (Fig. 2). The box plots clearly showed that tunicates evolved faster than other groups, especially compared to vertebrates, which were the slowest evolving. On average, tunicates evolved 6.25 times faster than vertebrates (two-tailed t test; t = 4.542, p < 0.001), 2.08 times faster than cephalochordates (two-tailed t test not applicable with only one cephalochordate), and 2.45 times faster than the outgroups (two-tailed t test; t = 1.711, p = 0.099) included here. The evolutionary rate variation was also much more pronounced within tunicates than within other groups, even when the very fast evolver Oikopleura dioica was excluded. For instance, the colonial species Diplosoma listerianum and Salpa fusiformis evolved considerably faster than the solitary species Ciona spp. and Styela plicata. This confirmed earlier observations based on a reduced number of taxa and substitution rate estimations on 35 housekeeping genes [16], once again underlining the peculiar genomic evolution of tunicates that might find its root in elevated mutation rates and pervasive molecular adaptation [21, 22].

Fig. 2
figure2

Evolutionary rate variation across sampled species. The bar plots represent average rate estimates (in number of substitutions per site per million years) obtained for the 63 terminal taxa regrouped by taxonomy. The rates were calculated using a rate-autocorrelated log-normal (LN) relaxed molecular clock model under the CAT-GTR + Γ4 mixture model with a birth-death prior on the diversification process and 13 soft calibration constraints. Data points are plotted as open circles with n = 10, 1, 18, 34 sample points in each taxonomic categories. Centre lines show the medians, crosses represent sample means, and box limits indicate the 25th and 75th percentiles with whiskers extending 1.5 times the interquartile range from the 25th and 75th percentiles. The width of the boxes is proportional to the square root of the sample size. This figure was made with BoxPlotR [81]

Full size image

Even though the difference in fit between the two relaxed clock models was not significant for our dataset, in general LN provided more consistent dating estimates than UGAM with respect to the mean divergence dates of numerous vertebrate groups reported in the latest phylogenomic study of jawed vertebrates [35]. Notably, as observed in a previous phylogenomic study of tetrapods [74], the application of the UGAM relaxed clock model provided unrealistically recent estimates with respect to the maximum node age for the origin of turtles (LN mean age + SD: 180 ± 19 Mya [95% credibility interval 220–146]; UGAM: 59 ± 41 Mya [173–16]) (Table 1, Fig. 3, and Additional file 2: Figure S1). The UGAM model also tended to systematically provide much wider 95% credibility intervals than LN, with several of them actually spanning hundreds of millions of years (Table 1, Fig. 1, and Additional file 2: Figure S1). Given the uncertainty associated with the dating results obtained using the UGAM model of clock relaxation, we focused our discussion below on results obtained with the more robust autocorrelated LN model, which we considered as our currently most reliable dating estimates.

Table 1 Molecular estimates of divergence dates (in Mya)
Full size table
Fig. 3
figure3

A molecular timescale for tunicates within chordates. The Bayesian chronogram has been obtained using a rate-autocorrelated log-normal (LN) relaxed molecular clock model using PhyloBayes under the CAT-GTR + Γ4 mixture model, with a birth-death prior on the diversification process and 13 soft calibration constraints. Node bars indicate the uncertainty around mean age estimates based on 95% credibility intervals. Plain node bars indicate nodes used as a priori calibration constraints. Numbers at nodes refer to Table 1

Full size image

A tentative timescale for tunicate evolution within chordates

The Bayesian chronogram obtained using the LN relaxed molecular clock model and the site-heterogeneous CAT-GTR + Γ4 mixture model of amino acid sequence evolution is presented in Fig. 3. This phylogenomic timescale showed that major tunicate clades appeared early in chordate evolutionary history. The earliest split between appendicularians and all other tunicates was dated back to ca. 450 Mya (mean age + SD: 447 ± 20 Mya [95% credibility interval 484–411]), followed by the divergence between stolidobranchs and the clade grouping thaliaceans + phlebobranchs + aplousobranchs ca. 390 Mya (389 ± 32 Mya [449–333]) and the separation of stolidobranchs into Molgulidae and Styelidae + Pyuridae ca. 350 Mya (350 ± 36 Mya [416–292]) (Table 1 and Fig. 3). Even more recent divergences such as the ones between congeneric species within Ciona and Molgula occurred more than 100 Mya.

Given the relative uncertainty on the phylogenetic position of Xenoturbella, complementary LN relaxed molecular clock analyses were also conducted using Xenoturbella as an outgroup. As the dating results previously obtained with the CAT-F81 + Γ4 and CAT-GTR + Γ4 models with the original rooting were extremely similar (linear regression on mean dates R2 = 0.99), we performed these additional analyses under the less computationally intensive CAT-F81 + Γ4 model. With the new rooting configuration, the inferred mean divergence dates between the two alternative rooting schemes were globally highly correlated within chordates (linear regression R2 = 0.89). An almost exact correspondence was found for vertebrates that contain most of the calibration points (linear regression R2 = 1.00). For tunicates, within which there is unfortunately no available calibration constraint, the correlation remained very strong (linear regression R2 = 0.95). The divergence dates within tunicates were on average older with the Xenoturbella rooting, while they remained in their vast majority within the original 95% credibility intervals (Additional file 3: Figure S2). An alternative rooting by Xenoturbella thus does not affect our main conclusions that divergence dates among the major tunicate lineages are ancient.

Our estimated divergence dates in tunicates were nevertheless associated with fairly large 95% credibility intervals, probably because of the lack of internal fossil calibrations within tunicates, in contrast to the well-calibrated vertebrates. It has recently been pointed out that, given the uncertainty associated with molecular dating estimates, building evolutionary narratives would be premature for early animal evolution [75]. In our case, we argue that in the absence of a trustable tunicate fossil record [33], our tentative molecular timescale constitutes the first and only currently available approach to provide a much-needed relative comparison of divergence times between the major lineages of tunicates and vertebrates. Such a comparison is subject to considerable uncertainty, but it has nevertheless revealed several deep divergences occurring at comparable geological times between the two groups (Fig. 3 and Table 2). For instance, these occur between tunicates (Ciona/ Oikopleura) and gnathostomes (Homo/ Callorhinchus) around 450 Mya; between thaliaceans (Salpa/ Doliolum) and lepidosaurs (Sphenodon/ Anolis) around 240 Mya; and between stolidobranchs (Molgula/ Botryllus) and tetrapods (Homo/ Xenopus) around 350 Mya.

Table 2 Parallel divergences between model tunicates and vertebrates
Full size table

The relatively ancient origins of the different tunicate lineages revealed by our molecular dating estimates have two broader implications. First, there seems to be a larger gap than previously thought between tunicate and vertebrate taxonomic ranks, which exacerbates the inadequacy of their direct comparison. For example, when a vertebrate genus usually spans less than 40 million years [76], a tunicate genus (e.g. Molgula) can span up to two hundred million years (Fig. 3 and Table 1). The meaning of Linnean categorical ranks and their temporal inconsistencies among clades have been largely discussed [76], as recently illustrated by the debate around the taxonomic status of the main chordate lineages [77,78,79]. The parallel we draw here between tunicates and vertebrates should nevertheless help tunicate developmental biologists to interpret their results in light of the large divergences that might exist between tunicate model species despite their classification in the same genus. Second, the ancient age of their major divergence events can heavily complicate orthology assessment among tunicates, as well as between tunicates and vertebrates, thus reducing the quality of genome annotations. Indeed, the fast-paced molecular evolution of tunicates prevents the identification of some genes by simple similarity methods (e.g. BLAST), even when orthologs do exist in databases. For instance, in terms of evolutionary depth, a comparative study of the genus Molgula is roughly equivalent to a comparative study among turtles representing about 180 million years of evolution. In terms of amino acid sequence divergence, the differences are much more pronounced between Molgula occidentalis / M. tectiformis (88.1% similarity) than between Phrynops/ Chrysemys (98.0% similarity; see Fig. 3 and Table 1).

From an evo-devo perspective, the phylogenetic framework and tentative timescale presented here lead to an apparent paradox. Like most nematodes [80], the embryos of each ascidian species develop in a stereotyped manner, based on the use of invariant cell lineages [7]. Unlike nematodes however, ascidian stereotyped cell lineages are shared between evolutionarily distant species such as Ciona robusta (Enterogona) and Halocynthia roretzi (Pleurogona) [11]. The extreme morphological conservation of ascidian embryogenesis therefore contrasts with the high rates of protein divergence observed in their genomes. This paradox raises questions about the underlying mechanisms involved in developmental regulation of these animals with highly dynamic genomes. In this context, our reference phylogenetic tree and divergence date estimates among tunicate lineages could be used as an evolutionary framework to select model species sufficiently close to one another (i.e. retaining sufficient phylogenetic information) for future comparative genomic analyses assessing orthology by gene tree reconciliation and estimating evolutionary rate variations among genes belonging to different ontology categories.

Conclusions

This study represents the first large-scale phylogenomic analysis including all major tunicate lineages based on transcriptomic data. The resulting phylogenetic framework and tentative timescale constitute a necessary first step towards a better understanding of tunicate systematics, genomics, and development, and in a broader context, of chordate evolution and developmental biology.

References

  1. 1.

    Bourlat SJ, Juliusdottir T, Lowe CJ, Freeman R, Aronowicz J, Kirschner M, et al. Deuterostome phylogeny reveals monophyletic chordates and the new phylum Xenoturbellida. Nature. 2006;444:85–8.

    CAS  Article  PubMed  Google Scholar 

  2. 2.

    Delsuc F, Brinkmann H, Chourrout D, Philippe H. Tunicates and not cephalochordates are the closest living relatives of vertebrates. Nature. 2006;439:965–8.

    CAS  Article  PubMed  Google Scholar 

  3. 3.

    Delsuc F, Tsagkogeorga G, Lartillot N, Philippe H. Additional molecular support for the new chordate phylogeny. Genesis. 2008;46:592–604.

    Article  PubMed  Google Scholar 

  4. 4.

    Putnam NH, Butts T, Ferrier DEK, Furlong RF, Hellsten U, Kawashima T, et al. The amphioxus genome and the evolution of the chordate karyotype. Nature. 2008;453:1064–71.

    CAS  Article  PubMed  Google Scholar 

  5. 5.

    Denoeud F, Henriet S, Mungpakdee S, Aury J-M, Da Silva C, Brinkmann H, et al. Plasticity of animal genome architecture unmasked by rapid evolution of a pelagic tunicate. Science. 2010;330:1381–5.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  6. 6.

    Schubert M, Escriva H, Xavier-Neto J, Laudet V. Amphioxus and tunicates as evolutionary model systems. Trends Ecol Evol. 2006;21:269–77.

    Article  PubMed  Google Scholar 

  7. 7.

    Lemaire P, Smith WC, Nishida H. Ascidians and the plasticity of the chordate developmental program. Curr Biol. 2008;18:R620–31.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  8. 8.

    Holland LZ. Genomics, evolution and development of amphioxus and tunicates: the Goldilocks principle. J Exp Zoolog B Mol Dev Evol. 2015;324:342–52.

    CAS  Article  Google Scholar 

  9. 9.

    Brunetti R, Gissi C, Pennati R, Caicci F, Gasparini F, Manni L. Morphological evidence that the molecularly determined Ciona intestinalis type A and type B are different species: Ciona robusta and Ciona intestinalis. J Zool Syst Evol Res. 2015;53:186–93.

    Article  Google Scholar 

  10. 10.

    Dehal P, Satou Y, Campbell RK, Chapman J, Degnan B, De Tomaso A, et al. The draft genome of Ciona intestinalis: insights into chordate and vertebrate origins. Science. 2002;298:2157–67.

    CAS  Article  PubMed  Google Scholar 

  11. 11.

    Lemaire P. Evolutionary crossroads in developmental biology: the tunicates. Development. 2011;138:2143–52.

    CAS  Article  PubMed  Google Scholar 

  12. 12.

    Small KS, Brudno M, Hill MM, Sidow A. Extreme genomic variation in a natural population. Proc Natl Acad Sci. 2007;104:5698–703.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  13. 13.

    Voskoboynik A, Neff NF, Sahoo D, Newman AM, Pushkarev D, Koh W, et al. The genome sequence of the colonial chordate, Botryllus schlosseri. elife. 2013;2:e00569.

    Article  PubMed  PubMed Central  Google Scholar 

  14. 14.

    Stolfi A, Lowe EK, Racioppi C, Ristoratore F, Brown CT, Swalla BJ, et al. Divergent mechanisms regulate conserved cardiopharyngeal development and gene expression in distantly related ascidians. elife. 2014;3:e03728.

    Article  PubMed  PubMed Central  Google Scholar 

  15. 15.

    Brozovic M, Martin C, Dantec C, Dauga D, Mendez M, Simion P, et al. ANISEED 2015: a digital framework for the comparative developmental biology of ascidians. Nucleic Acids Res. 2016;44:D808–18.

    CAS  Article  PubMed  Google Scholar 

  16. 16.

    Tsagkogeorga G, Turon X, Galtier N, Douzery EJP, Delsuc F. Accelerated evolutionary rate of housekeeping genes in tunicates. J Mol Evol. 2010;71:153–67.

    CAS  Article  PubMed  Google Scholar 

  17. 17.

    Seo H-C, Edvardsen RB, Maeland AD, Bjordal M, Jensen MF, Hansen A, et al. Hox cluster disintegration with persistent anteroposterior order of expression in Oikopleura dioica. Nature. 2004;431:67–71.

    CAS  Article  PubMed  Google Scholar 

  18. 18.

    Holland LZ, Gibson-Brown JJ. The Ciona intestinalis genome: when the constraints are off. BioEssays. 2003;25:529–32.

    CAS  Article  PubMed  Google Scholar 

  19. 19.

    Gissi C, Pesole G, Mastrototaro F, Iannelli F, Guida V, Griggio F. Hypervariability of ascidian mitochondrial gene order: exposing the myth of deuterostome organelle genome stability. Mol Biol Evol. 2010;27:211–5.

    CAS  Article  PubMed  Google Scholar 

  20. 20.

    Singh T, Tsagkogeorga G, Delsuc F, Blanquart S, Shenkar N, Loya Y, et al. Tunicate mitogenomics and phylogenetics: peculiarities of the Herdmania momus mitochondrial genome and support for the new chordate phylogeny. BMC Genomics. 2009;10:534.

    Article  PubMed  PubMed Central  Google Scholar 

  21. 21.

    Tsagkogeorga G, Cahais V, Galtier N. The population genomics of a fast evolver: high levels of diversity, functional constraint, and molecular adaptation in the tunicate Ciona intestinalis. Genome Biol Evol. 2012;4:852–61.

    Article  PubMed Central  Google Scholar 

  22. 22.

    Berna L, Alvarez-Valin F. Evolutionary genomics of fast evolving tunicates. Genome Biol Evol. 2014;6:1724–38.

    Article  PubMed  PubMed Central  Google Scholar 

  23. 23.

    Swalla BJ, Cameron CB, Corley LS, Garey JR. Urochordates are monophyletic within the deuterostomes. Syst Biol. 2000;49:52–64.

    CAS  Article  PubMed  Google Scholar 

  24. 24.

    Zeng L, Swalla BJ. Molecular phylogeny of the protochordates: chordate evolution. Can J Zool. 2005;83:24–33.

    CAS  Article  Google Scholar 

  25. 25.

    Tsagkogeorga G, Turon X, Hopcroft RR, Tilak M-K, Feldstein T, Shenkar N, et al. An updated 18S rRNA phylogeny of tunicates based on mixture and secondary structure models. BMC Evol Biol. 2009;9:187.

    Article  PubMed  PubMed Central  Google Scholar 

  26. 26.

    Govindarajan AF, Bucklin A, Madin LP. A molecular phylogeny of the Thaliacea. J Plankton Res. 2011;33:843–53.

    CAS  Article  Google Scholar 

  27. 27.

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

    Article  PubMed  PubMed Central  Google Scholar 

  28. 28.

    Shenkar N, Koplovitz G, Dray L, Gissi C, Huchon D. Back to solitude: solving the phylogenetic position of the Diazonidae using molecular and developmental characters. Mol Phylogenet Evol. 2016;100:51–6.

    Article  PubMed  Google Scholar 

  29. 29.

    Piette J, Lemaire P. Thaliaceans, the neglected pelagic relatives of Ascidians: a developmental and evolutionary enigma. Q Rev Biol. 2015;90:117–45.

    Article  PubMed  Google Scholar 

  30. 30.

    Berrill NJ. The Tunicata with an account of the British species [Internet]. London: Quaritch; 1950. https://catalog.hathitrust.org/Record/001500718

    Google Scholar 

  31. 31.

    Holland LZ. Tunicates. Curr Biol CB. 2016;26:R146–52.

    CAS  Article  PubMed  Google Scholar 

  32. 32.

    Shu D-G, Chen L, Han J, Zhang X-L. An Early Cambrian tunicate from China. Nature. 2001;411:472–3.

    CAS  Article  PubMed  Google Scholar 

  33. 33.

    Chen J-Y, Huang D-Y, Peng Q-Q, Chi H-M, Wang X-Q, Feng M. The first tunicate from the Early Cambrian of South China. Proc Natl Acad Sci U S A. 2003;100:8314–8.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  34. 34.

    Fedonkin MA, Vickers-Rich P, Swalla BJ, Trusler P, Hall M. A new metazoan from the Vendian of the White Sea, Russia, with possible affinities to the ascidians. Paleontol J. 2012;46:1–11.

    Article  Google Scholar 

  35. 35.

    Irisarri I, Baurain D, Brinkmann H, Delsuc F, Sire J-Y, Kupfer A, et al. Phylotranscriptomic consolidation of the jawed vertebrate timetree. Nat Ecol Evol. 2017;1:1370–8.

    Article  PubMed  PubMed Central  Google Scholar 

  36. 36.

    Romiguier J, Gayral P, Ballenghien M, Bernard A, Cahais V, Chenuil A, et al. Comparative population genomics in animals uncovers the determinants of genetic diversity. Nature. 2014;515:261–3.

    CAS  Article  PubMed  Google Scholar 

  37. 37.

    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.

    CAS  Article  PubMed  Google Scholar 

  38. 38.

    Simpson JT, Wong K, Jackman SD, Schein JE, Jones SJM, Birol I. ABySS: a parallel assembler for short read sequence data. Genome Res. 2009;19:1117–23.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  39. 39.

    Philippe H, Brinkmann H, Copley RR, Moroz LL, Nakano H, Poustka AJ, et al. Acoelomorph flatworms are deuterostomes related to Xenoturbella. Nature. 2011;470:255–8.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  40. 40.

    Baurain D. Forty-Two software [Internet]. https://bitbucket.org/dbaurain/42/downloads. Accessed 30 July 2017.

  41. 41.

    Philippe H. MUST, a computer package of Management Utilities for Sequences and Trees. Nucleic Acids Res. 1993;21:5264–72.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  42. 42.

    Castresana J. Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis. Mol Biol Evol. 2000;17:540–52.

    CAS  Article  PubMed  Google Scholar 

  43. 43.

    Roure B, Rodriguez-Ezpeleta N, Philippe H. SCaFoS: a tool for Selection, Concatenation and Fusion of Sequences for phylogenomics. BMC Evol Biol. 2007;7:S2.

    Article  PubMed  PubMed Central  Google Scholar 

  44. 44.

    Schmidt HA, Strimmer K, Vingron M, von Haeseler A. TREE-PUZZLE: maximum likelihood phylogenetic analysis using quartets and parallel computing. Bioinformatics. 2002;18:502–4.

    CAS  Article  PubMed  Google Scholar 

  45. 45.

    Stone M. Cross-validatory choice and assessment of statistical predictions. J R Stat Soc Ser B Methodol. 1974;36:111–47.

    Google Scholar 

  46. 46.

    Lartillot N, Lepage T, Blanquart S. PhyloBayes 3: a Bayesian software package for phylogenetic reconstruction and molecular dating. Bioinforma Oxf Engl. 2009;25:2286–8.

    CAS  Article  Google Scholar 

  47. 47.

    Lartillot N, Philippe H. A Bayesian mixture model for across-site heterogeneities in the amino-acid replacement process. Mol Biol Evol. 2004;21:1095–109.

    CAS  Article  PubMed  Google Scholar 

  48. 48.

    Lartillot N, Rodrigue N, Stubbs D, Richer J. PhyloBayes MPI: phylogenetic reconstruction with infinite mixtures of profiles in a parallel environment. Syst Biol. 2013;62:611–5.

    CAS  Article  PubMed  Google Scholar 

  49. 49.

    Benton MJ, Donoghue PCJ, Asher RJ, Friedman M, Near TJ, Vinther J. Constraints on the timescale of animal evolutionary history. Palaeontol Electron. 2015;18:1–106.

    Google Scholar 

  50. 50.

    Pyron RA. Divergence time estimation using fossils as terminal taxa and the origins of Lissamphibia. Syst Biol. 2011;60:466–81.

    Article  PubMed  Google Scholar 

  51. 51.

    Smith AB, Pisani D, Mackenzie-Dodds JA, Stockley B, Webster BL, Littlewood DTJ. Testing the molecular clock: molecular and paleontological estimates of divergence times in the Echinoidea (Echinodermata). Mol Biol Evol. 2006;23:1832–51.

    CAS  Article  PubMed  Google Scholar 

  52. 52.

    Thorne JL, Kishino H, Painter IS. Estimating the rate of evolution of the rate of molecular evolution. Mol Biol Evol. 1998;15:1647–57.

    CAS  Article  PubMed  Google Scholar 

  53. 53.

    Drummond AJ, Ho SYW, Phillips MJ, Rambaut A. Relaxed phylogenetics and dating with confidence. PLoS Biol. 2006;4:e88.

    Article  PubMed  PubMed Central  Google Scholar 

  54. 54.

    Lepage T, Bryant D, Philippe H, Lartillot N. A general comparison of relaxed molecular clock models. Mol Biol Evol. 2007;24:2669–80.

    CAS  Article  PubMed  Google Scholar 

  55. 55.

    Delsuc F, Brinkmann H, Philippe H. Phylogenomics and the reconstruction of the tree of life. Nat Rev Genet. 2005;6:361–75.

    CAS  Article  PubMed  Google Scholar 

  56. 56.

    Philippe H, Brinkmann H, Lavrov DV, Littlewood DTJ, Manuel M, Wörheide G, et al. Resolving difficult phylogenetic questions: why more sequences are not enough. PLoS Biol. 2011;9:e1000602.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  57. 57.

    Philippe H, Roure B. Difficult phylogenetic questions: more data, maybe; better methods, certainly. BMC Biol. 2011;9:91.

    Article  PubMed  PubMed Central  Google Scholar 

  58. 58.

    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.

    Article  PubMed  PubMed Central  Google Scholar 

  59. 59.

    Simion P, Philippe H, Baurain D, Jager M, Richter DJ, Di Franco A, et al. A large and consistent phylogenomic dataset supports sponges as the sister group to all other animals. Curr Biol CB. 2017;27:958–67.

    CAS  Article  PubMed  Google Scholar 

  60. 60.

    Cannon JT, Vellutini BC, Smith J, Ronquist F, Jondelius U, Hejnol A. Xenacoelomorpha is the sister group to Nephrozoa. Nature. 2016;530:89.

    CAS  Article  PubMed  Google Scholar 

  61. 61.

    Rouse GW, Wilson NG, Carvajal JI, Vrijenhoek RC. New deep-sea species of Xenoturbella and the position of Xenacoelomorpha. Nature. 2016;530:94.

    CAS  Article  PubMed  Google Scholar 

  62. 62.

    Ruiz-Trillo I, Paps J. Acoelomorpha: earliest branching bilaterians or deuterostomes? Org Divers Evol. 2016;16:391–9.

    Article  Google Scholar 

  63. 63.

    Telford MJ, Copley RR. Zoology: war of the worms. Curr Biol. 2016;26:R335–7.

    CAS  Article  PubMed  Google Scholar 

  64. 64.

    Shenkar N, Gittenberger A, Lambert G, Rius M, Moreira da Rocha R, Swalla BJ, et al. Ascidiacea World Database. [Internet]. World Register of Marine Species. 2017. http://www.marinespecies.org/aphia.php?p=taxdetails&id=1839. Accessed 30 July 2017.

  65. 65.

    Perrier E. Note sur la classification des Tuniciers. CR Acad Sci Paris. 1898;124:1758–62.

    Google Scholar 

  66. 66.

    Jue NK, Batta-Lona PG, Trusiak S, Obergfell C, Bucklin A, O’Neill MJ, et al. Rapid evolutionary rates and unique genomic signatures discovered in the first reference genome for the Southern Ocean salp, Salpa thompsoni (Urochordata, Thaliacea). Genome Biol Evol. 2016;8:3171–86.

    Article  PubMed  PubMed Central  Google Scholar 

  67. 67.

    Garstang W. Memoirs: The morphology of the Tunicata, and its bearings on the phylogeny of the Chordata. J Cell Sci. 1928;2:51–187.

    Google Scholar 

  68. 68.

    Stach T, Turbeville JM. Phylogeny of Tunicata inferred from molecular and morphological characters. Mol Phylogenet Evol. 2002;25:408–28.

    CAS  Article  PubMed  Google Scholar 

  69. 69.

    Turon X, López-Legentil S. Ascidian molecular phylogeny inferred from mtDNA data with emphasis on the Aplousobranchiata. Mol Phylogenet Evol. 2004;33:309–20.

    CAS  Article  PubMed  Google Scholar 

  70. 70.

    Kott P. The Australian Ascidiacea part 2, Aplousobranchia (1). Mem Qld Mus. 1990;29:1–266.

    Google Scholar 

  71. 71.

    Drummond AJ, Suchard MA, Xie D, Rambaut A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol. 2012;29:1969–73.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  72. 72.

    Rehm P, Borner J, Meusemann K, von Reumont BM, Simon S, Hadrys H, et al. Dating the arthropod tree based on large-scale transcriptome data. Mol Phylogenet Evol. 2011;61:880–7.

    Article  PubMed  Google Scholar 

  73. 73.

    dos Reis M, Gunnell GF, Barba-Montoya J, Wilkins A, Yang Z, Yoder AD. Using phylogenomic data to explore the effects of relaxed clocks and calibration strategies on divergence time estimation: primates as a test case. Syst Biol [Internet]. 2018. https://academic.oup.com/sysbio/advance-article/doi/10.1093/sysbio/syy001/4802240. Accessed 5 Feb 2018.

  74. 74.

    Chiari Y, Cahais V, Galtier N, Delsuc F. Phylogenomic analyses support the position of turtles as the sister group of birds and crocodiles (Archosauria). BMC Biol. 2012;10:65.

    Article  PubMed  PubMed Central  Google Scholar 

  75. 75.

    dos Reis M, Thawornwattana Y, Angelis K, Telford MJ, Donoghue PCJ, Yang Z. Uncertainty in the timing of origin of animals and the limits of precision in molecular timescales. Curr Biol. 2015;25:2939–50.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  76. 76.

    Avise JC, Liu J-X. On the temporal inconsistencies of Linnean taxonomic ranks. Biol J Linn Soc. 2011;102:707–14.

    Article  Google Scholar 

  77. 77.

    Satoh N, Rokhsar D, Nishikawa T. Chordate evolution and the three-phylum system. Proc R Soc B. 2014;281:20141729.

    Article  PubMed  PubMed Central  Google Scholar 

  78. 78.

    Lambertz M, Perry SF. Chordate phylogeny and the meaning of categorial ranks in modern evolutionary biology. Proc R Soc B. 2015;282:20142327.

    Article  PubMed  PubMed Central  Google Scholar 

  79. 79.

    Giribet G, Hormiga G, Edgecombe GD. The meaning of categorical ranks in evolutionary biology. Org Divers Evol. 2016;16:427–30.

    Article  Google Scholar 

  80. 80.

    Schulze J, Schierenberg E. Evolution of embryonic development in nematodes. EvoDevo. 2011;2:18.

    Article  PubMed  PubMed Central  Google Scholar 

  81. 81.

    Spitzer M, Wildenhain J, Rappsilber J, Tyers M. BoxPlotR: a web tool for generation of box plots. Nat Methods. 2014;11:121–2.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We thank Nicolas Galtier for giving early access to the Clavelina and Cystodytes Illumina data. We are thankful to Stefano Tiozzo, Hector Escriba, Sébastien Darras, Frédérique Viard, and the diving staffs of the marine stations of Villefranche-sur-Mer, Banyuls, and Roscoff for their help in sample collection. We also thank two anonymous referees for their comments. This is contribution ISEM 2018-024 of the Institut des Sciences de l’Evolution de Montpellier.

Funding

This work was supported by the Centre National de la Recherche Scientifique and the Agence Nationale de la Recherche (Contract ANR-13-BSV2–0011-01) to P.L and E.J.P.D. and by the Labex TULIP (ANR-10-LABX-41) to H.P. Computations were made on the Montpellier Biodiversity Bioinformatics (MBB) platform of the Labex CeMEB, and on the Mp2 and Ms2 supercomputers from the Université de Sherbrooke, managed by Calcul Québec and Compute Canada. The operation of this supercomputer is funded by the Canada Foundation for Innovation (CFI), the ministère de l’Économie, de la science et de l’innovation du Québec (MESI), and the Fonds de recherche du Québec - Nature et technologies (FRQ-NT).

Availability of data and materials

All data generated or analysed during this study are included in this published article and its additional files. Raw sequencing reads have been deposited under NCBI Bioproject PRJNA414754. Transcriptome assemblies, alignments, and trees are included in Additional file 4 and are also available from a GitHub repository (https://github.com/psimion/SuppData_Delsuc_BMCBiol_2018_Dating_Tunicata).

Author information

Author notes

    Affiliations

    1. ISEM, Université de Montpellier, CNRS, IRD, EPHE, Montpellier, France
      • Frédéric Delsuc
      • , Georgia Tsagkogeorga
      • , Paul Simion
      • , Marie-Ka Tilak
      •  & Emmanuel J. P. Douzery
    2. Centre for Biodiversity Theory and Modelling, UMR CNRS 5321, Station d’Ecologie Théorique et Expérimentale, Moulis, France
      • Hervé Philippe
    3. Département de Biochimie, Centre Robert-Cedergren, Université de Montréal, Montréal, Canada
      • Hervé Philippe
    4. School of Biological and Chemical Sciences, Queen Mary University of London, London, UK
      • Georgia Tsagkogeorga
    5. Center for Advanced Studies of Blanes (CEAB, CSIC), Girona, Spain
      • Xavier Turon
    6. Department of Biology and Marine Biology, Center for Marine Science, University of North Carolina Wilmington, Wilmington, NC, USA
      • Susanna López-Legentil
    7. Centre de Recherche en Biologie cellulaire de Montpellier, UMR 5237, CNRS, Université de Montpellier, Montpellier, France
      • Jacques Piette
      •  & Patrick Lemaire
    Authors
    1. Frédéric Delsuc
      View author publications
      You can also search for this author in
    2. Hervé Philippe
      View author publications
      You can also search for this author in
    3. Georgia Tsagkogeorga
      View author publications
      You can also search for this author in
    4. Paul Simion
      View author publications
      You can also search for this author in
    5. Marie-Ka Tilak
      View author publications
      You can also search for this author in
    6. Xavier Turon
      View author publications
      You can also search for this author in
    7. Susanna López-Legentil
      View author publications
      You can also search for this author in
    8. Jacques Piette
      View author publications
      You can also search for this author in
    9. Patrick Lemaire
      View author publications
      You can also search for this author in
    10. Emmanuel J. P. Douzery
      View author publications
      You can also search for this author in

    Contributions

    FD, HP, and EJPD designed the study. FD, GT, XT, SL-L, and JP collected and prepared biological material. FD, GT, M-KT, JP, PL, and EJPD organized sequence data collection. HP and GT constructed the supermatrix. FD, HP, PS, and EJPD analysed data. FD, HP, PS, and EJPD drafted the manuscript. All authors contributed to the final version of the manuscript and gave final approval for publication.

    Corresponding author

    Correspondence to Frédéric Delsuc.

    Ethics declarations

    Ethics approval and consent to participate

    Not applicable.

    Competing interests

    The authors declare that they have no competing interests.

    Publisher’s Note

    Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

    Additional files

    Additional file 1:

    Table S1. Origin of biological samples, sequencing, assembly statistics, and accession numbers of tunicate transcriptomes. (DOCX 105 kb)

    Additional file 2:

    Figure S1. Bayesian chronogram obtained using an uncorrelated gamma (UGAM) relaxed molecular clock model using PhyloBayes under the CAT-GTR + Γ4 mixture model, with a birth-death prior on the diversification process and 13 soft calibration constraints. Node bars indicate the uncertainty around mean age estimates based on 95% credibility intervals. Plain black node bars indicate nodes used as a priori calibration constraints. Numbers at nodes refer to Table 1. (PPTX 97 kb)

    Additional file 3:

    Figure S2. Bayesian chronogram obtained using an autocorrelated log-normal (LN) relaxed molecular clock model using PhyloBayes under the CAT-F81 + Γ4 mixture model, with a birth-death prior on the diversification process, 13 soft calibration constraints, and an alternative rooting by Xenoturbella. Node bars indicate the uncertainty around mean age estimates based on 95% credibility intervals. Plain black node bars indicate nodes used as a priori calibration constraints. Numbers at nodes refer to Table 1. (PPTX 111 kb)

    Additional file 4:

    Supporting data. Transcriptome assemblies, alignments, and trees. (ZIP 38502 kb)

    Rights and permissions

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

    Reprints and Permissions

    About this article

    Verify currency and authenticity via CrossMark

    Cite this article

    Delsuc, F., Philippe, H., Tsagkogeorga, G. et al. A phylogenomic framework and timescale for comparative studies of tunicates. BMC Biol 16, 39 (2018). https://doi.org/10.1186/s12915-018-0499-2

    Download citation

    Keywords

    • Tunicata
    • Thaliacea
    • Molecular dating
    • Transcriptomes
    • Phylogenomics
    • Evo-devo

    BMC Biology

    ISSN: 1741-7007

    Contact us