Revisiting the phylogeny of phylum Ctenophora: a molecular perspective

Introduction

Several phylogenetic hypotheses of the phylum Ctenophora based on morphological data¹, ribosomal markers^2,3, protein-coding markers⁴, have been proposed, all of them through different approaches.

Due to the poor fossil record of the group, the morphological data of fossil taxa have been not enough to help to resolve this question, because it’s impossible to determine which characteristics arose first between the ctenophores, for several reasons such as the poor conservation state of the available fossils^5,6 or the lack of shared characteristics between the extant ctenophores and the extinct ctenophores⁷. Also some morphological characteristic have been demonstrated to be homoplastic^1,3. These situations slow down the process of reconstructing the phylogeny of this phylum.

There are also difficulties establishing an appropriate outgroup, due to the unknown position of this phylum inside Metazoa, many hypotheses have been suggested^8–11. This uncertainty could bias the phylogenetic analyses if a distant outgroup is chosen (eg. highly saturated sequences, not homologous sequences available, etc), affecting directly the support values and the topology of the reconstructed tree¹².

In this study we reconfirm the paraphyly of order Cyddipida similarly to all previous studies, also we confirm the order Lobata is paraphyletic, the reasons are exposed in the results section. Nevertheless, our data don’t support the paraphyly of Beroidae. Due to the fairly wide taxonomic sampling of this study resulting from the fusion of the protein-coding and ribosome data, some interesting relationships are suggested, such as the placement of Cestida and Thalassocalycida orders inside Lobata.

Methods

The ribosomal sequences were obtained from public data available on GenBank and automatically downloaded, then they were classified using python scripts. The NFP sequences was also obtained from GenBank and for certain taxa supplied by Steve Haddock via e-mail, from the Supplementary data of the reporting study of the marker⁴, we only included sequences from Ctenophora. The Tyrosine aminotransferase sequences and the HLH domain-containing protein were also obtained from public data of GenBank. The accession numbers of the sequences used during this study are presented in Table 1.

Previous of the concatenated final analysis, we tested several markers such as:

We execute a single locus analysis for all these markers. The ribosomal sequences were aligned by MAFFT v.7.7¹³ with the option –auto. The proteing coding sequences (NFP, Tyrosine aminotransferase and HLH domain-containing protein) were aligned using RevTrans2 (http://www.cbs.dtu.dk/services/RevTrans-2.0/web/)¹⁴.

Models for single locus analyses were selected with two programs: jModelTest 2.1.10¹⁵ for nucleotide datasets (ribosomal data), and ProtTest v 3.4.2 for protein markers¹⁶.

Single locus analyses were performed by partitions obtained with Gblocks 0.91b (http://molevol.cmima.csic.es/castresana/Gblocks_server.html)^17,18.

For phylogenetic reconstruction using the concatenated ribosomal dataset the pipeline PhyPipe was used¹⁹ (available at: https://gitlab.com/cibiop/phypipe/).

The Bayesian inference analyses (BI) were all performed by MrBayes 3.2.6²⁰ and the Maximum Likelihood (ML) analyses were performed by GARLI 2.01²¹ and RAxML 8.0.0²².

NFP was previously introduced as an ortholog by 4, which function is still unknown. So we took this marker along the ribosomal data as the backbone of the alignment matrix (and the study).

Concatenation of the sequences for the study (NFP,HLH,Tyr, 18S, 5.8S, ITS1, ITS2) was performed by 2matrix 1.0 (https://github.com/nrsalinas/2matrix)²³. This script allows the automatic concatenation of a heterogeneous matrix, and also convert the concatenated matrix to the files input of the Maximum Likelihood and Bayesian Inference programs.

Partition and selection of the models for the concatenated matrix were performed by PartitionFinder2²⁴, separately from Ribosomal data and protein-coding markers. The best scheme files are available inside the Supplementary data.

The models used for protein-coding markers in all analyses are: For NFP, JJT+I+G²⁵; for Tyrosine amino-transferase, LG+G²⁶, and for HLH, VT+I²⁷.

The set of parameters for Bayesian Inference analysis are reported in the Supplementary data inside the NEXUS file, this analysis were performed in CIPRES²⁸, we used 8 MCMC with 10’000.000 generations by duplicate, this allows an optimal performance of the analysis. For 18S and 5.8S the analysis was performed with HKY+I+G²⁹, for ITS1 and ITS2 with SYM+G³⁰.

RAxML analysis was performed in CIPRES²⁸ with 20 independent maximum likelihood analyses and 10.000 bootstrap iterations (pseudoreplicates) for nucleotide partitions, and the model used in this analysis was GTR+G+I³¹. The importance of invariant proportion executing an analysis with RAxML in this specific dataset is explained in Discussion section.

GARLI analysis was performed with 10 independent maximum likelihood analyses and 1004 bootstrap iterations (pseudoreplicates) for nucleotide partitions we used more specific models. For 18S we used TrN+I+G³²; for 5.8S, TIMef+I+G; for ITS1, TIMef+G, and for ITS2, TVMef+G. These models were indicated by jModeltest 2.1.10¹⁵.

The models were all selected by BIC criterion. In this study we did not used an outgroup.

In order to obtain a more complete matrix, we fused sequences from few species of the same genus into a single record, for example, we fused Hormiphora plumosa and Hormiphora californiensis into a single Hormiphora sp.; Bolinopsis sp. and Bolinopsis infundibulum into a single Bolinopsis sp. and Lampea lactea and Lampea pancerina into a single Lampea sp. These fused species allow an improvement of the alignment matrix and the phylogenetic reconstruction. Also few species were duplicated, such as Beroe forskalii, Bathyctena chuni, Hormiphora sp., because the HLH marker presented variation amongst the species and was not possible to obtain a consensus. We confirm the monophyly of this variation through a single locus analysis, as mentioned before.

Results

Combined dataset

The tree reconstructed from the combined dataset (protein + ribosomal DNA) is presented on Figure 1. The results from both Maximum Likelihood analyses (RAxML and GARLI) for the combined dataset (protein + ribosomal DNA) are similar, except in the specific relationships between Eurhamphaeidae + Cestida + Leucotheidae + Bolinopsidae. RAxML results matches with Mr-Bayes results, but three nodes of the analysis have low posterior probability (BI) or low bootstrap values (ML). RAxML analysis shows a clade composed by Eurhamphaeidae and Leucotheidae, and other clade composed by Cestida and Bolinopsidae. Whereas GARLI analysis shows Eurhamphaeidae as sister taxa of Leucotheidae, Cestida and Bolinopsidae. RAxML results are similar to 3. All of our analyses show Cestida within Lobata with high bootstrap values and posterior probability, defining Lobata as a clade composed by Leucotheidae, Eurhampaheidae, Bolinopsidae and Ocyropsidae. Bathocyroidae and Lampoctenidae families have an uncertain position between Lobata and the clade composed by Beroe sp. and Haeckelidae.

Figure 1. Unrooted tree reconstructed with RAxML for the combined dataset (NF, HLH, Tyr, 18S, 5.8S, ITS1, ITS2).

Support values of nodes (bootstrap values) are shown on tree branches.

Our analysis support a clade including Thalassocalycida and Lampoctenidae but the position of this clade remains controversial due to the lack of high bootstrap values. The family Bathocyroidae forms a clade with Dryodoridae family, but this clade has a low posterior probability and low bootstrap values, so for now it is not accurate to set hypothesis around this result, this results is similar to obtained by 4. Also the position of Dryodora glandiformis is still undetermined, in ML analysis this family could group with even Pleurobrachidae, due to the low bootstrap values of nodes between Lobata and Pleurobrachidae and with BI analysis with all lobates. Undescribed species T forms a good supported clade with Bathocyroidae. Further studies may focus in describing this taxon for morphological purposes. We executed a rogue taxa analysis through RogueNaRok³³, we found that Beroe ovata, Beroe cucumis, Beroe gracilis, Lampocteis cruentiventer, Dryodora glandiformis, UCS4, Llyria B, and Lyrocteis sp. were rogue taxa during this analysis, the low bootstrap values could be related to this.

The relationship of Bathyctenidae family (Represented by Bathyctena chuni) with the Mertensiidae family and the Platyctenida order remains unclear, this family shows affinity to this clade, also in several times forms a clade with two undescribed Mertensiids (A9 and undescribed sp3), this two taxa are excluded of the family Mertensiidae (Represented by Mertensia ovum, Charistephane fugiens).

The identity of spB remains unclear. species spC forms a good supported clade with Lampeidae by both methods.

Discussion

Next steps for the resolving of the phylogeny of this group is to determine who is the most basal branch inside the Ctenophora, making possible and reconstruction of ancestral characters. Inside the upper clade formed by Pleurobrachiidae, Haeckelidae, Beroida, Thalassocalycida, Dryodoridae(with low support values), Lobata and Cestida, could be very crucial the reconstruction of ancestral characters for the understanding of the plasticity of the characteristic inside this group. This could only achieved by setting a good outgroup for this group.

We strongly recommend for further studies, the identification and posterior description for undescribed species. Also an enrichment of the aligment matrix produced by this study, through sequencing crucial markers such as 18s, ITS1 and NFP, which played an important role in the reconstruction of the phylogeny presented in this study.

Also, we recommend for further studies an extensive sampling of groups like Pleurobrachidae, in an attempt to collapse the long branches presented in previous studies³ and the present study, genus as Tinerfe, which present morphological similarities with families as Haeckeliidae¹. More sampling outside groups such as Lobata would allow an improvement for further studies.

The proportion of invariant sites, plays an important role in the analysis of Ribosomal data. During the analysis for this paper, we noticed that the absence of this feature in the analysis in RAxML, forms a clade composed by Beroe and Pleurobrachidae as sister taxa of Lobata, Cestida and Thalassocalycida, this clade was of course with an extremely low bootstrap value; the presence of this feature presents Pleurobrachidae as the sister taxa of Beroe, Lobata, Cestida and Thalassocalycida (Presented in suplementary data). So the absence or presence of this feature during the analysis should be relevant.

Data availability

The raw data used for this project are available in Zenodo, DOI 10.5281/zenodo.838689 (Arteaga-Figueroa et al., 2016).