Articles – PLOS Currents Tree of Life Wed, 08 Mar 2017 22:58:33 +0000 en-US hourly 1 Red Algal Phylogenomics Provides a Robust Framework for Inferring Evolution of Key Metabolic Pathways Fri, 02 Dec 2016 17:00:49 +0000 Red algae comprise an anciently diverged, species-rich phylum with morphologies that span unicells to large seaweeds. Here, leveraging a rich red algal genome and transcriptome dataset, we used 298 single-copy orthologous nuclear genes from 15 red algal species to erect a robust multi-gene phylogeny of Rhodophyta. This tree places red seaweeds (Bangiophyceae and Florideophyceae) at the base of the mesophilic red algae with the remaining non-seaweed mesophilic lineages forming a well-supported sister group. The early divergence of seaweeds contrasts with the evolution of multicellular land plants and brown algae that are nested among multiple, unicellular or filamentous sister lineages. Using this novel perspective on red algal evolution, we studied the evolution of the pathways for isoprenoid biosynthesis. This analysis revealed losses of the mevalonate pathway on at least three separate occasions in lineages that contain Cyanidioschyzon, Porphyridium, and Chondrus. Our results establish a framework for in-depth studies of the origin and evolution of genes and metabolic pathways in Rhodophyta.



Red algae (Rhodophyta) form a monophyletic lineage containing ~7,000 described species1 that exhibit a wide variety of morphological and ultra-structural forms and have complex reproductive strategies. The Cyanidiophytina (e.g., Galdieria and Cyanidioschyzon) include extremophiles that thrive in volcanic areas surrounding hot springs. In contrast, their mesophilic sisters (Rhodophytina) are globally distributed from freshwater environments to open oceans and deep oceans (>200 m) to the intertidal zone. Despite a highly reduced core gene inventory that resulted from an ancient phase of genome reduction2, red algae represent one of the few eukaryotic lineages that have evolved complex multicellularity3, typified by red seaweeds such as Porphyra and Gracilaria. Red seaweeds account for ~95% of known red algal taxa and are important sources of agricultural (e.g., nori) and industrial products (e.g., agar and carrageenan).

Studies of red algal systematics have largely relied on a handful of plastid and nuclear genes4,5,6,7,8 and focused on a broad diversity of lineages within the Florideophyceae9,10. One of the major findings of these analyses is the separation of Cyanidiophytina from the Rhodophytina4,8. Whereas Cyanidiophytina contain only two known families (Cyanidiaceae and Galdieraceae), Rhodophytina encompass six classes: Bangiophyceae, Florideophyceae, Compsopogonophyceae, Porphyridiophyceae, Rhodellophyceae, and Stylonematophyceae4. Excluding the well-supported monophyly of Bangiophyceae and Florideophyceae (hereafter, collectively referred to as red seaweeds), relationships among the remaining classes remain controversial4,5,6,7,8.

In this study, we applied phylogenomics to a rich genomic dataset to erect a robust red algal tree of life. The dataset encompassed 298 orthologous nuclear-encoded genes from all major red algal lineages. In contrast to previous phylogenies built using smaller datasets4,5,6,7,8, our results support a fundamental, ancient split between red seaweeds and non-seaweed lineages among mesophiles. We discuss the implication of this new perspective on red algal phylogeny to understanding the evolution of multicellularity in red algae, and demonstrate the utility of this phylogenetic framework to infer the evolution of the mevalonate (MVA) pathway of isoprenoid biosynthesis in Rhodophyta.


Construction of single-copy orthologous gene alignments

We created a local database that includes protein sequences (translated from EST or predicted from genome sequences) from 15 red algal taxa2,11,12,13,14,15,16 (Fig. 1A) and 3 green algae17,18,19 (Table 1, Appendix 1). This database, after removing short sequences with length <100 amino acids, was used in a self-query using BLASTp (e-value cutoff = 1e-5). The BLASTp search output was used as input for OrthoMCL20 with parameters (evalueExponentCutoff = -10, percentMatchCutoff = 40, inflation = 1.5) to construct orthologous gene families. Among these families, we searched for single-copy orthologous genes with one gene copy per species (allowing missing data in up to three red algae and in no more than one green alga). For each orthologous gene family, the corresponding sequences were retrieved and aligned using MUSCLE (version 3.8.31) under the default settings21. The alignments were then trimmed using TrimAl (version 1.4)22 in automated mode (-automated) and then ‘polished’ with T-COFFEE (version 9.03)23 to removed poorly aligned residues (conservation score ≤ 5) among the aligned blocks. A total of 298 single-gene alignments (length >150 amino acids and with ≥15 sequences) were retained for downstream analysis.

Construction of the multi-protein phylogeny

The 298 single-copy gene alignments were concatenated into a super-protein alignment. A phylogenetic tree was inferred using Phylobayes (version 3.3)24 under the CAT model25. This is a mixture model that takes into consideration site-specific evolutionary properties (such as rate and profile) within the alignment25. The CAT model generally fits data significantly better than one-matrix models such as LG and WAG. We set up two chains that ran in parallel and assessed convergence periodically using ‘bpcomp’ and ‘tracecomp’ functions. Convergence assessments were done based on sampled trees (taking one from every 10 trees) following burnin equal to 20% of the entire length of the chain. The two chains were stopped when they converged to an acceptable level that allows good qualitative measurement of the posterior consensus. According to the user instructions (, an acceptable run corresponds to a maximum discrepancy across all bipartitions (maxdiff <0.3) when monitored with the ‘bpcomp’ function, and statistical discrepancies <0.3 and effective sizes >50 for all parameters when monitored with the ‘tracecomp’ function.

Construction of coalescence model-based species trees

We built a coalescence model-based red algal phylogeny with 100 replicates following Seo’s method26. For each replicate, we randomly sampled 298 genes with replacement. For each sampled alignment, a pseudo-alignment was generated by random sampling of amino acid site from the original alignment with replacement. Only one green algal sequence (as outgroup) was retained with the priority given to Chlamydomonas reinhardtii, Chlorella variabilis, and Micromonas RCC299 in order. A ML tree was built for each pseudo-alignment using IQtree (version 0.9.6)27 under the best-fit amino acid evolutionary model selected on the fly (-m TEST). The resulting 298 ML trees, rooted with outgroup sequences, were then used for maximum pseudo-likelihood tree construction using MP-EST (version 1.4) under the default settings28. This procedure was repeated 100 times and the resulting 100 maximum pseudo-likelihood trees were summarized under majority rule using the ‘consense’ function in Phylip (

Phylogenetic analyses of mevalonate pathway genes

Galdieria sulphuraria proteins in the MVA pathway (module identifier: M00095) and the methylerythritol phosphate (MEP) pathway (module identifier: M00096) were retrieved from the KEGG database29 and used as queries against NCBI (nr) using BLASTp (e-value cutoff = 1e-5) ( The representative sequences (e.g., from Metazoa and land plants) were retrieved from Genbank. Local BLASTp searches (e-value cutoff = 1e-5) were done against our red algal database aforementioned followed by retrieval of the significant hits. Galdieria phlegrea sequences were retrieved from the previous study30. Each G. sulphuraria query, together with the homologs (from Genbank and our local database), were aligned using MUSCLE (version 3.8.31)21 under the default settings. The alignment was trimmed using trimAl (version 1.4)22 in the automated mode (-automated). ML trees were built using IQtree (version 0.9.6)27 under the best amino acid evolutionary model selected using (-m TEST) with branch support values estimated using 1,500 ultrafast bootstrap replicates (-bb 1500). The resulting trees were manually inspected. Distantly related paralogs (if any) were removed manually and the trees were rebuilt following the procedure described above.

Validation of gene losses in red algae

We searched for the G. sulphuraria MVA and MEP proteins in a red algal nucleotide database (genome and transcriptome) using tBLASTn (e-value cutoff = 1e-5). The homologous protein sequences translated from the hit nucleotide sequences were collected using an in-house script. For each query sequence, the translated proteins corresponding to the three top bit-score hits and the three top-identity (query-hit identity) hits were incorporated into the single-gene ML tree building procedure described above. Distantly related homologs were manually identified and removed. Red algal sequences that were monophyletic with G. sulphuraria were considered to be orthologs.

Results and Discussion

Red algal tree of life

We constructed single-gene alignments for a total of 298 one-to-one orthologous genes (98,494 amino acid positions in total) that are conserved in 15 red algal and 3 green algal taxa (see Methods). Analysis of the concatenated super-protein alignment under the CAT model led to a highly supported phylogenetic tree that received 1.00 posterior probability for all interior nodes (Fig. 1A). This tree confirmed the early split between Cyanidiophytina and Rhodophytina4,8 and monophyletic relationship between Bangiophyceae and Florideophyceae4,8. The relationships within Florideophyceae are consistent with previous analyses10,31 with Hildenbrandiophycidae (Hildenbrandia) in the basal position. Nemaliophycidae (Palmaria) is sister to the monophyletic group containing Corallinophycidae (Calliarthron) and Rhodymeniophycidae (Chondrus)10,31. The remaining non-seaweed mesophilic lineages formed a robust monophyletic group, with Stylonematophyceae in the basal position. Compsopogonophyceae formed a sister group to the monophyletic Porphyridiophyceae and Rhodellophyceae.

Concatenation-based analysis has previously been shown in some instances to result in inflated statistical support for incorrect topologies32 due to heterogeneity across genes and gene-specific evolution, such as gene duplication33. To minimize this problem, we used a tree summarization approach that does not rely on the concatenation of multiple single-gene alignments. This method takes a population of single-gene trees as input and estimates the species tree using a coalescence model28. This analysis led to the same tree topology (Fig. 1A) to the concatenation-based analysis with high bootstrap support for the monophyletic group comprising red seaweeds (bootstrap support = 100%) and non-seaweed mesophilic red algae (bootstrap support = 90%). The relationships among non-seaweed red algal lineages are however weakly supported (bootstrap support = 49-51%). Taken together, our phylogenomic analyses strongly support a separation between seaweeds and non-seaweed lineages at the base of mesophilic red algae (Fig. 1A).

(A) A phylogenetic tree inferred from a concatenated alignment of 298-proteins. The outgroup species are not shown. Statistical supports (separated by a back slash) for each branch are derived from the super-protein analysis (posterior probability) and from the coalescence model-based analysis (bootstrap support). (B) Schematic representation of the positions of red seaweeds and land plants (thick branches) in red algae and Viridiplantae, respectively. The phylogenies are derived from this study (panel I), Scott et al (Ref. 6, panel II) and Leliaert et al. (Ref. 35, panel III). The arrows indicate genome reduction (GR). Bangiophyceae (Bangio.), Compsopogonophyceae (Compsopogo.), Cyanidiophyceae (Cyanidio.), Florideophyceae (Florideo.), Porphyridiophyceae (Porphyridio.), Stylonematophyceae (Stylonemato.), Coleochaetophyceae (Coleochaeto.), Chlorokybophyceae (Chlorokybo.), Klebsormidiophyceae (Klebsormidio.), Mesostigmatophyceae (Mesostigmato.), Zygnematophyceae (Zygnemato.).

Fig. 1: Red algal phylogenomics


The early divergence of red seaweeds within mesophilic red algae is consistent with the antiquity of the Bangiophyceae. A putative bangiophyte (Bangiomorpha pubescens) has been found in rocks dated at ca. 1.2 billion years old34. This result suggests the existence of distinct fates for the two lineages that split from the common ancestor of mesophilic red algae. One lineage remained unicellular with the development of filaments in some species (e.g., Rhodochaete and Purpureofilum), whereas the other developed complex filamentous plant bodies leading to red seaweeds with bi- and tri-phasic life cycles. The basal position of red seaweeds among mesophiles (Fig. 1B, scenario I) contrasts with previous analyses4,5,6,8,10,34 (Fig. 1B, scenario II) and to the highly derived position of land plants in Viridiplantae35. Land plants evolved within streptophyte green algae that have simpler morphological forms (Fig. 1B, scenario III)35. Likewise, kelps (phaeophytes) are also nested among multiple unicellular stramenopile lineages36. Among mesophilic red algae, red seaweeds appear to have ‘recovered’ from the extensive genome reduction they shared with the red algal ancestor2 and that was further exacerbated in Cyanidiophytina due to their extremophilic lifestyles30. The early (>1 billion year old) emergence of the multicellular lineage is all the more remarkable when placed in context to the early evolution of red algae. Alternatively, the lack of simpler lineages (in terms of morphology) in the red seaweed clade may suggest their high rate of extinction or the existence of yet unknown species that remain to be discovered in this clade. An early emergence of a peculiar type of multicellularity (green seaweed Palmophyllophyceae) was also discovered recently in the basal position of Chlorophyta37.

Parallel losses of MVA pathway

To demonstrate the usefulness of this novel perspective on red algal phylogeny, we used the reference tree to elucidate the evolution of the isopentenyl pyrophosphate (IPP) biosynthetic pathway. IPP is the building block of isoprenoids that comprises a large diversity of lipids found in all three domains of life. In photosynthetic eukaryotes, two independent pathways exist to produce IPP, the cytosolic and peroxisome localized MVA pathway and the plastid MEP pathway38. Whereas the MEP pathway is conserved across many species, the MVA pathway has been lost in green algae (Chlorophyta)38 and in some red algal lineages such as C. merolae16 and P. purpureum12. Our analysis of red algal sequence data (see Methods) showed that the MEP pathway is present in all examined lineages. The minor gene losses that were found are most likely to be explained by missing data commonly associated with transcriptome datasets (Fig. 3, Appendix 2). In contrast, the MVA pathway is largely absent (3rd to 6th enzymes in the pathway, Fig. 2A) in most red algal lineages except the Stylonematophyceae (Rhodosorus marinus and Purpureofilum apyrenoidigerum) and G. sulphuraria. Presence of the MVA pathway in G. sulphuraria39 and Cyanidium caldarium40 is supported with genetic and biochemical evidence39,40. This result suggests that loss of MVA pathway is more widespread than previously thought. The red algal origin of the MVA genes in Stylonematophyceae is supported with phylogenetic data (see Methods). For example, in the phylogeny of HMG-CoA reductase (HMGR, Fig. 2B), R. marinus and P. apyrenoidigerum form a monophyletic group with and Galdieria species, whereas no other red algae were present in this clade. A similar pattern is found for other MVA pathway genes that were lost in most red algal species (Fig. 4, Appendix 3).


Fig. 2: MVA pathway in red algae

(A) The distribution of MVA pathway genes across red algal species. Black and open circles denote the presence and absence of the genes, respectively. For each gene, the gray boxes indicate gene presence for the corresponding classes. Arrows indicate genome reduction. Red vertical bars indicate gene losses. ACAT (acetyl-CoA acetyltransferase), HMGS (hydroxymethylglutaryl-CoA synthase), HMGR (3-hydroxy-3-methylglutaryl-CoA reductase), MVK (mevalonate kinase), PMK (phosphomevalonate kinase), MVD (mevalonate decarboxylase), IDI (isopentenyl-diphosphate delta-isomerase). (B) A ML tree of HMGR. The taxa in red color: red algae, green: Viridiplantae, orange: chromalveolates, brown: Opisthokonta.

Absence of the MVA pathway in all five sampled red seaweeds suggests it was most likely lost in their common ancestor. BLASTp searches (e-value cutoff = 10) against nucleotide databases (expressed sequence tag and transcriptome shotgun assembly) in NCBI did not return any significant hits to MVA pathway genes from Bangiophyceae and Florideophyceae. In addition, their losses in C. merolae, P. purpureum, and C. crispus that have both transcriptome and genome data available are well supported. Given the red algal phylogeny (Fig. 1A), these losses were unambiguously resulted from three parallel events (Fig. 2A). Under this scenario, the MVA pathway survived the ancient phases of genome reduction (arrows, Fig. 2A) and underwent gene loss more recently after the split of the seaweed and non-seaweed lineages. MVA pathway loss in C. merolae likely resulted from an additional phase of genome reduction specific to this lineage30 (Fig. 2A). The selective forces that led to the retention or loss of the MVA pathway across the mesophilic red algal lineages are presently unknown. Nonetheless, MVA pathway loss suggests that IPP biosynthesis is dependent on the plastid MEP pathway and requires transporters for the export of IPP from the plastid to the cytosol38. The MVA pathway was also lost in Chlorophyta (including most unicellular green algae)38 and G. sulphuraria is physiologically distinct from mesophilic species. For this reason, the discovery of possible MVA pathway-containing and -absent lineages among mesophilic red algae provides an algal model for studying the evolution of isoprenoid biosynthesis and intracellular trafficking among compartments.


Our phylogenomic analyses resulted in a well-supported red algal phylogeny that provides new insights into the evolution of red seaweeds. Our results will allow more accurate reconstruction of evolutionary events (e.g., gene family evolution2 and molecular calibration10) and provide a framework to map the distribution of red algal functions and traits. Further efforts are needed to substantiate the relationships among non-seaweed mesophilic red algae with high quality genome data from these taxa41.

Data Availability

The multi-protein alignment is available for download (ID: 20087) from TreeBASE (

Competing Interests

The authors have declared that no competing interests exist.

Corresponding Author

Huan Qiu, Department of Ecology, Evolution and Natural Resources, Rutgers University, New Brunswick, NJ 08901, USA.


Appendix 1

Table 1

Algal genome and transcriptome data used for the phylogenomic analysis

Classification Species Source Data type MMETSP ID
Seaweed Hildenbrandia rubra Ref. 2 Transcriptome
Seaweed Palmaria palmata Ref. 2 Transcriptome
Seaweed Calliarthron tuberculosum Ref. 14 Partial genome
Seaweed Chondrus crispus Ref. 13 Whole genome
Seaweed Porphyra umbilicalis Ref. 14 Transcriptome
Mesophiles Purpureofilum apyrenoidigerum Ref. 2 Transcriptome
Mesophiles Rhodochaete pulchella Ref. 2 Transcriptome
Mesophiles Rhodosorus marinus Ref. 11 Transcriptome MMETSP0315
Mesophiles Rhodella maculata Ref. 11 Transcriptome MMETSP0167
Mesophiles Compsopogon coeruleus Ref. 11 Transcriptome MMETSP0312
Mesophiles Erythrolobus australicus Ref. 11 Transcriptome MMETSP1353
Mesophiles Timspurckia oligopyrenoides Ref. 11 Transcriptome MMETSP1172
Mesophiles Porphyridium purpureum Ref. 12 Transcriptome
Extremophiles Galdieria sulphuraria Ref. 15 Transcriptome
Extremophiles Cyanidioschyzon merolae Ref. 16 Whole genome
Green algae Chlorella variabilis Ref. 17 Whole genome
Green algae Chlamydomonas reinhardtii Ref. 18 Whole genome
Green algae Micromonas pusilla Ref. 19 Whole genome

Appendix 2


Fig. 3: Distribution of the MEP pathway across red algal lineages

Black and open circles denote the presence and absence of the genes, respectively. For each gene, the gray boxes indicate the gene presence for the corresponding classes. DXS (1-deoxy-d-xylulose 5-phosphate synthase), DXR (1-deoxy-d-xylulose 5-phosphate reductoisomerase), MCT (2-C-methyl-d-erythritol 4-phosphate cytidylyltransferase), CMK (C-methyl-d-erythritol kinase), MDS (2-C-methyl-d-erythritol 2,4-cyclodiphosphate synthase), HDS (4-hydroxy-3-methylbut-2-en-1-yl diphosphate synthase), HDR (4-hydroxy-3-methylbut-2-en-1-yl diphosphate reductase), IDI (isopentenyl-diphosphate isomerase).

Appendix 3


Fig. 4: ML trees for six MVA pathway genes


]]> 0
How Really Ancient Is Paulinella Chromatophora? Tue, 15 Mar 2016 16:40:57 +0000 The ancestor of Paulinella chromatophora established a symbiotic relationship with cyanobacteria related to the Prochloroccocus/Synechococcus clade. This event has been described as a second primary endosymbiosis leading to a plastid in the making. Based on the rate of pseudogene disintegration in the endosymbiotic bacteria Buchnera aphidicola, it was suggested that the chromatophore in P. chromatophora has a minimum age of ~60 Myr. Here we revisit this estimation by using a lognormal relaxed molecular clock on the 18S rRNA of P. chromatophora. Our time estimates show that depending on the assumptions made to calibrate the molecular clock, P. chromatophora diverged from heterotrophic Paulinella spp. ~ 90 to 140 Myr ago, thus establishing a maximum date for the origin of the chromatophore.



Mitochondria and plastids evolved from free­living bacteria by symbiogenesis more than one billion years ago1. Both events boosted the evolution of eukaryotes by expanding their metabolic abilities. Primary endosymbiosis leading to organelles (i.e., mitochondria and plastids) was thought to be unique in the history of life until the recent discovery of an independent primary endosymbiosis in Paulinella chromatophora2. This thecate filose amoeba hosts in its cytoplasm photosynthetic organelles of cyanobacterial origin, called chromatophores. Phylogenetic analysis of 16S rRNA showed that the chromatophores originated from marine α­cyanobacteria from the Prochlorococcus/Synechococcus clade. In contrast, “classical” plastids of plants and algae evolved from an ancient unknown lineage of cyanobacteria3.

The genome of the cyanobacteria Synechococcus WH5701, one of the closest known free­-living relatives of the chromatophore with a sequenced genome, has 2917 protein­-coding genes4. In contrast, the chromatophores of P. chromatophora strains FK01 and M0880/a contain between 841 and 867 protein­-coding genes respectively5,6. Remaining genes in the chromatophore suggest a strong metabolic interdependence with the amoebal nucleocytoplasm. The discovery that proteins of the chromatophore photosynthetic apparatus are encoded in the host genome and imported back into the cyanobacterial­-derived compartment, reinforces the suggestion that the chromatophore is a bona fide primary organelle7.

It is likely that the origin of the chromatophore is one or two orders of magnitude more recent than the establishment of the primary plastids of plants and algae. But, how ancient is the chromatophore? Since there is no fossil record of P. chromatophora or its close relatives, it is difficult to have a precise answer to this question. However, an initial guess suggested the chromatophore has a minimum age of 60 Myr. This was based on the proposal that in Buchnera aphidicola (Aphid’s endosymbiotic bacteria) a pseudogene needs between 40-60 Myr to disintegrate completely8. Since the genome of the chromatophore has pseudogenes, the same tempo was extrapolated for the chromatophore6.

Surprisingly, the suggestion that the chromatophore in P. chromatophora has a minimum age of 60 Myr underwent a change in part of the subsequent literature indicating that P. chromatophora has ~60 Myr. Clearly, this is an assertion that requires clarification and further analysis. The recent identification by single cell genomics of non-­photosynthetic close relatives of P. chromatophora9 and novel attempts to calibrate the origin of major eukaryotic groups10,11 offer an opportunity to revisit the origin in time of this extraordinary symbiosis.


Taxon sampling. Based on previous published phylogenetic analyses we retrieved a sample of 18S rRNA sequences from rhizaria and stramenopila from SILVA database12. These included sequences from: i) the phylogenetic analysis whereby close relatives of P. cromatophora were identified9; ii) the phylogenetic analyses describing divergence times of major eukaryotic groups10,11; and iii) sequences from phylogenetic analyses of euglyphids13 and diatoms14. We did not include sequences from foraminifera due to their high rate of evolution reflected in extreme long branches. Sequences were aligned with MUSCLE15. The final alignment contained 43 sequences and 1128 aligned positions without gaps. The phylogenetic tree reconstructed from these sequences (see below) is in general terms congruent with published phylogenetic analyses.

Molecular clock analyses. BEAUTi was used to prepare xml files. Time trees were inferred by using a lognormal relaxed molecular clock as implemented in BEAST 216. To select the model of evolution we used jModelTest17. The Tamura and Nei 1993 model (TrN) plus gamma distribution and an estimated proportion of invariant sites were rated best by the bayesian information criterion (BIC). Therefore, all the analyses were conducted under this model of evolution. In addition, we used a Yule tree prior for all analyses.

To calibrate the clock, we relied upon several sources of information (Fig 2). These included organismal as well as chemical fossils. We also used previous published estimates of the time divergence between rhizaria and stramenoplia. These sources of information were organized into four different calibration schemes. All four schemes used information from non­-controversial fossil record (evidences: a, b, c and d). However, the schemes differ on the use of soft evidence from previous molecular clock analyses (evidence: e.1 and e.2); and on the use of evidence from vase­shaped microfossils (VSM) described by18. VSM were originally assigned to rhizaria. However, a posterior molecular clock analysis favored a re­interpretation of VSM as members of amoebozoa10. Thus, the four calibration schemes were: (a, b, c, d, e.1); (a, b, c, d, e.1, f); (a, b, c, d, e.2); (a, b, c, d).

Priors on the age of nodes were adjusted to reflect confidence on time divergence. For instance, we assigned a strong prior on the origin of rhizosolenid diatoms because the origin of this lineage has been dated to 91.5 ± 1.5 Myr based on chemical fossils14. The 95% probability distribution on this case is between 91.5 to 97.0 Myr. Other calibrations received less strong priors, as indicated by the larger number of years contained along the 2.5% to 9.5% quantiles of the shape column from Table 1. Priors used for the origin of pennate diatoms and for the origins of earliest diatoms were based on those used by11 and represent minimal divergence times. The prior used for the origin of Eugliphydae was based on the supposition that the group is much older than the fossil evidence as suggested by19. This prior also represents a minimal divergence time. Finally, to assign a prior to the divergence in time of rhizaria from stramenopila we first averaged BEAST time estimates published by11 and then we parametrized a normal distribution to include the variability of these estimates. We used this normal distribution as a prior. We constructed a second prior for the same divergence of rhizaria from stramenopila but now based on the inference published by10. All sequences within each prior were restricted to be monophyletic. The root of the tree was determined by making rhizaria monophyletic.

We evaluated each calibration scheme by looking at convergence and ESS > 200 in Tracer after running chains of length 1010 and sampling each 104 generations. We discarded 10% of generations as burnin. Finally, we conducted all analyses without sequences (i.e., sampling from priors) to determine the impact of priors and to test whether sequences are informative on estimated divergence dates. Trees were obtained with TreeAnnotator and visualised with FigTree. We provide BEAUTi xml files of the four calibration schemes for BEAST 2 analyses here.

Results and discussion

In Fig. 1 we show an estimate of the divergence in time of P. chromatophora by using a lognormal relaxed clock on 18S rRNA. To calibrate this tree we relied on: a) the origin of rhizosolenid diatoms, which is known with high confidence (91.5 ± 1.5 Myr)14; b) a minimal time divergence of pennate diatoms (80 Myr)11,20; c) a minimal time divergence for diatoms (133.9 Myr)11,21; d) a minimal time divergence of Euglyphidae 40 Myr19; and e.1) a time estimation of the divergence of rhizaria from stramenopila ~1232 Myr11. This gave us an estimation for the origin of P. chromatophora of 93.6 Myr (38.1 – 138.2).


Fig. 1: The origin of P. chromatophora in time according to SSU rRNA lognormal molecular clock.

Tree calibrated under the scheme (a, b, c, d, e.1). Small open red circles indicate calibrations based on fossil record (a, b, c, and d); small open blue circle indicate calibration based on previous estimation (e.1); red full circle, estimated age of P. chromatophora. Notice that the divergence in time of Plasmodiophorabrassicae (black open circle) is consistent with the proposed origin in time of forams. Species names in orange: stramenopila; species name in blue and green: rhizaria.


Fig. 2: Calibration constrains.

Evidence: used to derive priors for calibration constraints; Min: Minimal divergence time (Offset parameter in BEAUTi) in Myr; Distribution: Used to model each prior together with their respective parameter values (alpha α and beta β for gamma γ and mean and sigma σ for normal distributions); Shape: shows the median for each distribution and the values containing the 2.5% and 97.5% quantiles. Stronger priors span few years between the 2.5% and 97.5% quantiles.

However other calibration schemes are possible. For instance, we can use the report of vase­-shaped microfossils (VSM) of ~742­ to 770 Myr (calibration f) that were described tentatively as members of Euglyphida18 together with calibrations a, b, c, d and e.1, to get an estimate of 141.4 Myr (48.7 – 210.4) for the origin of P. chromatophora. Alternatively, we can follow the suggestions that VSMs belong to amoebozoa and that rhizaria diverged from stramenopila 754.1 Myr (639.8 – 903.3)10 together with calibration points a, b, c, d and e.2 to get an estimate of 55.8 Myr (25.4 – 78.6) for the origin of P. chromatophora. Finally, we can estimate the divergence of P. chromatophora based only on calibration points a, b, c and d; thus avoiding controversial VSM fossils (calibration f) and the soft information provided from previous estimates of the divergence between rhizaria and stramenopila (calibrations e.1 and e.2). By this, we get an estimate of 47.9 Myr (28.8 – 64.9).

Which of the four time estimates is closer to the true divergence time of P. chromatophora? Foraminifera offer a clue. Foraminifera belong to rhizaria and as well as diatoms have a detailed fossil record. However, their 18S sequences are too divergent to include them in our analysis. The first appearance of forams in the fossil record is dated at 545 Myr22. Recent molecular time estimates of forams suggest this group originated 770 Myr (650 – 920)23. A phylogenetic analysis based on 15 concatenated genes shows that forams are a sister group of Plasmodiophorabrassicae with the exclusion of Gromiaoviformis as an outgroup11. Therefore, the origin of the lineage leading to P. brassicae in our tree has to be at least as old as the proposed origin of forams. By looking at our time estimates, it happens that only the two proposals that are based on a time of ~1232 Myr for the origin of rhizaria (i.e., that uses calibration e.1) result in a date for the origin of P. brassicae that is congruent with the estimated origin of forams (Fig. 3).


Fig. 3: Concordance between different time divergence inferences and earliest forams in the fossil record.

Only calibration schemes assuming a time divergence of rhizaria from stramenopila of ~1232 Myr (e.1) are consistent with the fossil record and proposed time divergence of forams. The horizontal black broken line represents the date of the oldest foram fossil at 545 Myr. The blue broken line represents the estimated origin of forams 770 Myr (650 – 920).

Although the 95% confidence intervals are rater large, the estimates of 93.6 and 141.4 Myr are consistent with the original suggestion that the chromatophore in P. chromatophora has a minimum age of 60 Myr6; and with a more recent proposal of a maximum age of 200 Myr for the divergence of P. chromatophora from Euglypha rotunda24.

If our estimates are correct, the two strains of P. chromatophora diverged from each other about (45.7 – 64.7 Myr). The genomes of the chromatophores in these two strains differ in about 33 genes5. This give us a rate of ~ 0.25 to 0.36 gene inactivations per million year since the divergence of the two strains. These show that genome reduction in the chromatophore is a process that continues slowly however stately. Whether these genes are been lost by genetic drift or natural selection is still a mater of research. Further refinements on the time of origin of P. chromatophora and their chromatophores are expected as our understanding of the evolution of rhizaria and the eukaryotic tree of life improves.

Character State Reconstruction of Call Diversity in the Neoconocephalus Katydids Reveals High Levels of Convergence Fri, 11 Mar 2016 16:14:18 +0000 The katydid genus Neoconocephalus is characterized by high diversity of the acoustic communication system. Both male signals and female preferences have been thoroughly studied in the past. This study used Bayesian character state reconstruction to elucidate the evolutionary history of diverse call traits, based on an existing, well supported phylogenetic hypothesis. The most common male call pattern consisted of continuous calls comprising one fast pulse rate; this pattern is the likely ancestral state in this genus. Three lines of call divergence existed among the species of the genus. First, four species had significantly slower pulse rates. Second, five species had alternating pulse periods, resulting in a double pulse rhythm. Third, several species had discontinuous calls, when pulses were grouped into rhythmically repeated verses. Bayesian character state reconstruction revealed that the double-pulse pattern likely evolved convergently five times; the slow pulse rate also evolved four times independently. Discontinuous calls have evolved twice and occur in two clades; each of which contains reversals to the ancestral continuous calls. Pairwise phylogenetically independent contrast analyses among the three call traits found no significant correlations among the character states of the different traits, supporting the independent evolution of the three call traits.



One of the most intriguing open questions in biology is, how did the diversity of species, forms, adaptations and behaviors that we find today evolve? In particular, how do novel traits evolve, i.e. traits without obvious precursors? Acoustic communication systems with their diversity of call patterns and call recognition mechanisms (i.e., call preferences) have long served as a model for this question (reviews in 1,2)

Communication signals often account for a large part of the phenotypic diversity among closely related species3,4. The rapid diversification of communication systems5 enables them to play a major role in speciation; in turn, speciation might also be a driver for the diversification of communication (e.g. 6,7). Signal diversity may be under selection from the environment, and/or constrained by many biological factors3,8. This complexity insures that the diversification of communication systems remains a controversial and important issue.

Understanding the processes that generate and maintain phenotypic diversity requires a solid hypothesis of phylogenetic relationships and knowledge of the variation in the trait(s) under consideration (e.g. 9. Here we present the character state analysis of diverse call traits in the katydid genus Neoconocephalus (Orthoptera, Tettigoniidae), utilizing a well supported phylogenetic hypothesis10.

The number of species under consideration here is relatively small compared to similar studies (e.g. 11,12). However, a plethora of data exists regarding acoustic communication of Neconocephalus and this group has long been established as model for novel trait evolution (reviews in 13,14). In this system males produce calls by rubbing their forewings together, mainly to attract receptive females. The calls are species specific in both temporal and spectral properties. Indeed, the male calls are among the most diverse phenotypic traits in this group and are in many cases the species defining traits. Many aspects of the communication system have been studied: male call pattern (e.g. 15,16,17,18), call mechanics19,20,21,22, female preferences (e.g. 23,24,25,26), and sensory processing (e.g. 27,28,29,30,31). Additionally, a well supported phylogeny for Neoconocephalus exists, encompassing most of the Central American, Caribbean and North American species10.

Three lines of call diversity exist within this group. Two of them (call structure and pulse pattern) are discrete traits, each with two character states. The third, pulse rate, is quantitative in nature. However, based on the sensory processing of this call character, we treat it here with two character states (fast and slow, see below). Here we present Bayesian character state reconstructions for these three call traits. We then tested for potential links among the diversity of these traits.

Materials and Methods

We performed a Bayesian character state reconstruction using the Neoconocephalus total evidence tree of Snyder et al.10 that combined AFLP markers, three nuclear gene markers (ITS1, ITS2, H3) and one mitochondrial gene (CO1). The methodology underlying this tree is given in10. The total evidence tree included 17 Neoconocephalus species from North and Central America as well as the Caribbean islands. Taxon sampling included all Neoconocephalus species occurring in wide geographic ranges. Missing species encompass one endemic species from North America and species endemic to two Caribbean islands. Bucrates malivolans and Belocephalus davisi were used as out group in 10; these species were not included in the character state reconstruction of this study.

Character State Reconstruction for discrete Call Traits

For the three lines of call divergence (see below) we used 10,000 post-burnin trees from the Bayesian total evidence phylogeny to reconstruct the ancestral character states with the program Discrete, implemented in BayesTraits v1.032,33 This method reconstructs the most probable character states at each tested node, maximizing the likelihood of the character states observed in terminal taxa. Call state reconstructions were completed in batch mode using the concatenated data matrix by running MCMC, 2 million generations 100k burnin, at a 0.2 acceptance rate.

Phylogenetically Independent Contrasts for Discrete Call Traits

We used generalized estimating equations (GEE), which handle categorical data and multichotomies found in a tree34,35, to analyze discrete call traits. The compare.gee command requires gee (v.4.13-1836), nlme (v3.1-11837), and ape (v3.338), packages in R. We ran this analysis using RStudio v0.98.1091 and R v3.1.2. This command runs phylogenetic independent contrasts pair-wise between evolutionary history of the three discrete call traits and the species relationships represented in10. These results are then compared to a students t-test to look for significance; additionally, GEE calculates the phylogenetic number of degrees of freedom dfP34.

Call traits and character states of Neoconocephalus calls

The calls of all Neoconocephalus species included in this analysis have been described in the literature, as reviewed in13, with only a few details added since that time18. The published descriptions of the call traits of interest were confirmed during the taxon collection of 10 and during numerous studies of our research group. Neoconocephalus calls are produced during opening and closing of the forewings (elytra); loud sound pulses are produced during the closing, while sounds during the opening are typically much softer22 and are not required for female responses (e.g.17,18,24,25; but see 26). We therefore limit our description of the temporal pattern to the closing pulses and refer to them simply as ‘pulses.’ We describe each species’ calls for the three call traits ‘call structure’, ‘pulse pattern’, and ‘pulse rate.’ The character states for these three call traits are listed in Table 1.

Table 1: Description of Neoconocephalus calls

The character states for the three call traits of the 17 species considered. C continuous, D discontinuous (=versed) call structure; Pulse rates < 100 Hz are classified as ‘slow’, > 100Hz as fast. The call structure in N. ensiger was equivocal, indicated by a ? (see text); * the verse duration for N. ensiger as given here interprets individual pulses as chirps (see text).

Species Call Structure Verse Duration Verse Rate[Hz] Pulse Pattern Pulse Rate[Hz] Source
N. affinis C double slow 10-12 1526
N. bivocatus C double fast 155-175 1739
N. caudelianus D 300 ms 1.1 single fast 192 13
N. ensiger ? 40 ms* 10-15 single slow 10-14 13222742
N. exciliscanorus D 87 ms 3.5 single slow 83 13
N. maxillosus C double fast 265-280 1518
N. melanorhinus C single fast 140 51
N. nebrascensis D 1200-1400ms 0.5 single fast 200 172251
N. palustris C single fast 225 51
N. punctipes C single fast 245 15
N. retusiformis C single slow 49 15
N. retusus C double fast 140-170 1822
N. robustus C single fast 200 172239
N. saturatus D 50-80ms 2-4 single fast 161 15
N. spiza D 50-80ms 2-4 single fast 153 15
N. triops D 500-1000ms 1-2 double fast 180-250 152240
N. velox C single fast 195 51

The simplest and most common call structure in Neoconocephalus consists of pulses produced at a constant, fast pulse rate in the range of 150-250 Hz, produced for extended periods of time (minutes up to hours) without a rhythmic second order time structure. During such calls, the wing movement never stopped at either full closing or full opening22; accordingly opening and closing pulses were never separated by more than a few ms. We refer to this call as ‘continuous’ call structure, with ‘single pulse’ pattern and a ‘fast pulse rate’ (Fig. 1A, 3A).

Pulse Pattern: Most species of Neoconocephalus produce calls that consist of a single pulse rate, i.e. the pulses are evenly spaced (Fig. 1A). In some species, however, shorter and longer pulse periods alternate. The resulting pulse pairs or ‘double pulses’ are then regularly repeated at a ‘double pulse rate,’ which is equivalent to one half of the individual pulse rate (Fig. 1B). The wing movement during one double pulse begins with a complete opening of the wings followed by a partial closing. Wings open then fully again followed by complete closure.

The wing movement does not stop either within the pulse pair, or between pulse pairs in all Neoconocephalus species with this pattern, independently of the double pulse rate i.e. no silent intervals longer than a few ms occur in these species. Five species have double pulse pattern (Fig. 1): N. bivocatus, N. affinis, N. retusus, N. maxillosus, and N. triops15,18,39.

Pulse Rate: Calls of most Neoconocephalus species have unusually fast pulse rates for Tettigoniids (150-250/s)13. Several species have dramatically slower pulse rates at about 50-80/s (N. retusiformis, N. exciliscanorus) or 10-20/s (N. ensiger, N. affinis). The fast pulse rates above about 150/s are at or beyond the limit of the temporal resolution for the sensory system and are perceived as continuous signals without amplitude modulation17. Thus, small scale changes of pulse rates are unlikely to be important for female preferences. In contrast, the slow pulse rates below 100 pulses/s can be resolved by the sensory system17. We chose to treat this trait as discrete: classifying the calls with pulse rates ≤ 100 Hz as ‘slow pulse rate’ and all other species as ‘fast pulse rate’ (Fig. 2, Table 1). For species with double pulse pattern (see above), we counted the each individual pulse (i.e. two pulses per pulse pair) when calculating the pulse rate. For the treatment of N. ensiger, see below.

Call Structure: Several species had a distinct second order time structure in their calls: pulses were produced in distinct ‘verses’ (= chirps, echemes), which were regularly repeated after silent intervals (Fig. 3B). Durations of chirps and silent intervals ranged from 50ms/250ms to 2s/1.5s. In N. triops, the silent intervals were short at 50-100ms, compared to verse durations that varied between 500 and 1000ms among populations (Table 1). Calls with such rhythmic second order time structures were classified as ‘discontinuous’ call structures, while calls without a rhythmic verse pattern, i.e. pulses produced continuously, were classified as ‘continuous’ call structures (Fig 3C). For the treatment of N. triops and N. ensiger see below.

N. triops: The range of this species extends throughout the Neotropics into temperate North America. Throughout its range, N. triops has a discontinuous call with 500-1000 ms long verses repeated after silent intervals of 50-100ms duration15,22,40. In North America, although the summer generation has the same discontinuous call structure as all other populations, the winter generation expresses an alternative call type, which is continuous40,41. This alternative call type is the result of developmental plasticity (i.e. is environmentally induced) and the same genotypes can express both phenotypes, depending on the environmental conditions. Because the alternative phenotype is limited to extreme environmental conditions in one population, we classified the call of N. triops as ‘discontinuous.’

N. ensiger: This species produces a slow pulse rate (about 10/s) continuously. In contrast to all other Neoconocephalus species considered here, the wing movement stops completely between pulses for more than 30 ms, so that the opening and closing movement of the wings cover only about 1/2 of the wing cycle. Accordingly, silent intervals of 30-40 ms occur between pulses13. This is longer than the typical interpulse durations in this genus, but shorter than the silent gaps between verses in species with discontinuous call pattern (Fig 1C). In addition, the within-male variability of pulse rate of N. ensiger’s call was much higher (at constant temperature) than that of other Neoconocephalus species42.

The long silent intervals and the variability in pulse rate raise the question of whether the pulse rate of N. ensiger is homologous to the pulse rate of the other species. Alternatively, N. ensiger’s pulses could represent a ‘one-pulse’ verse, and the pulse rate may actually be homologous to the verse rate of discontinuous species. As this question cannot be resolved at this time, we considered the call structure of N. ensiger as ‘unknown’ (Fig. 3). The combined duration of opening and closing pulses in N. ensiger is around 30-40 ms, which is equivalent to pulse rates of 25-33 Hz. Therefore, the interpretation of this call as “slow pulse rate” is not affected by the question stated above.


Pulse Pattern: The character state analysis revealed the single pulse pattern (Fig. 1) as the likely ancestral state for this call trait (posterior probability [pp] = 0.99). The double pulse pattern appeared at five tips in the tree, supporting convergent evolution of this character. All five state changes were highly supported (≥ 0.99).

Character state reconstruction of the trait "pulse pattern" of Neoconocephalus calls. Call characters of extant species are indicated by shapes at the tips (squares = double pulse, dots = single pulse). Reconstructed character states are represented by branch colors (black = single, blue = double pulse); values at nodes represent posterior probability of character states. Inset: Oscillograms of A single pulse call, B double pulse call.

Fig. 1: Character states of ‘pulse pattern’


Pulse Rate: This call trait showed a similar pattern as the pulse pattern (Fig. 2). Fast pulse rate was revealed as the ancestral state (pp = 0.99). Slow pulse rate appears at four tips, each with a highly supported character change (pp ≥ 0.99).

Character state reconstruction of the trait "pulse rate" of Neoconocephalus calls. Call characters of extant species are indicated by shapes at the tips (squares = slow, dots = fast). Reconstructed character states are represented by branch colors (black = fast, yellow = slow); values at nodes represent posterior probability of character states. Calls with pulse rates > 140 Hz were classified as 'fast', pulse rates below 100 Hz as 'slow' (see methods). Inset: Oscillograms of calls with A fast pulse rate and B slow pulse rate. The call shown in B is from N. affinis, which has double pulse structure; the call of N. ensiger shown in Fig. 3C also has slow pulse rate.

Fig. 2: Character states of ‘pulse rate’

Character state reconstruction of the trait “pulse rate” of Neoconocephalus calls. Call characters of extant species are indicated by shapes at the tips (squares = slow, dots = fast). Reconstructed character states are represented by branch colors (black = fast, yellow = slow); values at nodes represent posterior probability of character states. Calls with pulse rates > 140 Hz were classified as ‘fast’, pulse rates below 100 Hz as ‘slow’ (see methods). Inset: Oscillograms of calls with A fast pulse rate and B slow pulse rate. The call shown in B is from N. affinis, which has double pulse structure; the call of N. ensiger shown in Fig. 3C also has slow pulse rate.

Call structure: The phylogenetic pattern of this call trait differed from that of the other two traits (Fig. 3). Discontinuous calls occurred in two distinct clades and the character state analysis indicated two separate origins with high posterior probability (pp = 0.99 and pp = 0.98). The ancestral state for Neoconocephalus was accordingly a continuous call structure (pp = 0.93). Within each clade with discontinuous calls at least one reversal to continuous calls took place. The treatment of N. triops as discontinuous caller did not affect the outcome of this character state analysis: coding this species as ‘unknown’ or continuous affected only the character state of this species and shifted the character change one node upwards (data not shown).

Character state reconstruction of the trait "call structure" of Neoconocephalus calls. Call characters of extant species are indicated by shapes at the tips (squares = discontinuous, dots = continuous). Reconstructed character states are represented by branch colors (black = continuous, red = discontinuous); values at nodes represent posterior probability of character states. The red/black dashed line indicates ambiguous character states; the red/white dashed line and ? indicate the unknown character state of the N. ensiger call (see methods). Inset: Oscillograms of A continuous call, B discontinuous call, and C N. ensiger call; note the different time scale in C

Fig. 3: Character states of ‘call structure’


Ancestral call state in Neoconocephalus: The character state analyses indicate a continuous, single pulsed call with fast pulse rate as the ancestral call pattern; support for each of these characters was >0.93. This call type occurred in 6 out of 17 extant species; for each individual call trait the ancestral character is more common than the derived character (Figs. 1-3, Table 1).

Correlation of derived traits: We conducted pairwise phylogenetically independent contrasts analyses among the three call traits. None of these comparisons produced a significant correlation between character states of the different traits (Table 2). This supports the conclusion that the three call traits evolve independently.

Table 2: Phylogenetically independent contrasts among call traits

Results of pairwise comparison using generalized estimating equations. The coefficient of regression (± s.e.m) is given. Phylogenetic number of degrees of freedom dfp = 4.68.

coefficient ± s.e.m p-value
pulse pattern vs. call pattern 0.22 ± 0.33 0.42
pulse rate vs. call pattern 0.14 ± 0.34 0.44
pulse rate vs. pulse pattern 0.3 ± 0.27 0.39


The Bayesian character state analysis revealed convergent evolution for the derived states of all three call traits. However, the patterns differed among the three call traits: slow pulse rates and double pulses occurred only at tips of the species tree while the discontinuous call structure occurred in two clades. We found no evidence for correlated evolution among the three call traits.

We identified the ancestral call state of Neoconocephalus calls as a continuous call with fast pulse rate and single pulse pattern. Ancestral call structure and pulse pattern are less complex than their derived character states. Calls are produced by rubbing the forewings together; the central pattern generator (CPG) producing the rhythmic pulse pattern was likely derived from the wing beat generator43, which typically produces a rate between 10 and 30 Hz (= wing-beat frequency during flight). The ancestral pulse rate in Neoconocephalus is derived (and highly unusual) for katydids, where pulse rates are typically below 50-60 Hz44. The closest relative of the new-world genus Neoconocephalus is the old-world genus Ruspolia; indeed, Neoconocephalus is likely a monophyletic group inside of Ruspolia13. The most prevalent call pattern in Ruspolia is similar to the ancestral Neoconocephalus call, with pulse rates around 100 Hz (i.e. slower than in Neoconocephalus, yet still fast for katydids) 44,45,46.

In four species, pulse rate reverted back towards slower values. These slow amplitude modulations (AM) are important for female preferences in N. affinis and N. ensiger26,42; N. retusiformis calls have a distinct slow AM rate, but we lack data on female preferences. Male N. exciliscanorus typically have seven opening and closing movements in each chirp22. A pulse pattern, however, is not reliably detectable in the resulting sound. Females do not require a pulse pattern for phonotaxis and respond to chirps without internal AM (KHF & JS, unpublished). The differences in motor patterns resulting in slow rates make it unlikely that similar neural or genetic changes underlie the pulse rates in these four species.

Double pulses are a common pattern in katydid calls, occurring in numerous genera across the different Tettiogniid subfamilies14,44. Double pulse pattern in Neoconocephalus is unusual in that the wing movement never stops22 and no silent gaps longer than a few ms occur. This is true even in the slow calls (pulse rate 10-12 Hz) of N. affinis. In all other double pulsing katydid species where data on wing movements are available, the wings come to a complete stop in the closed position44, resulting in longer silent gaps between double pulses. The unusual motor pattern in Neoconocephalus, and its consistency among the five species, suggests that similar neural changes underlie the independent origins of double pulses in this group. The distribution of double pulses suggests that this pre-disposition was present basally in the Neoconocephalus clade.

The discontinuous call structure represents a second order time structure modulating the pulse pattern. This most likely requires the integration of a second pattern generator into neural call production. The CPG generating the breathing rhythm is a likely candidate for modulating Orthopteran calls 47,48. Discontinuous calls occurred in two clades, most likely with independent origins. Character state analysis revealed reversals to the ancestral continuous call structure within each clade. In one clade, both species with discontinuous calls (N. saturatus, N. spiza) have very similar verse structure (verse duration 50-80ms, verse rates of 2-4 Hz). In the other clade, verse structure is much more variable (Table 1) across species. Whether similar genetic changes underlie the reorganization of the pattern generator (i.e. the inclusion of a second order time structure) is not clear at this time.

We detected no correlations among the character states of the three call traits. This suggests that the neural changes underlying these character states were independent from each other. Accordingly, there are species with a single derived call trait and others with multiple derived traits (N. affinis: double pulses & slow pulse rate; N. triops: double pulses & discontinuous calls).

Female call recognition mechanisms have been characterized by identifying which combination of call parameters were attractive for females of numerous species (review in 14). In species with ancestral call pattern, female recognition relies on the absence of detectable silent intervals: gaps longer than a few ms render the calls unattractive17. No specific AM besides the absence of gaps was required. This mechanism occurs in all species with the ancestral call (N. robustus17, N. palustris, N. punctipes; Schul, unpublished data) tested so far. Species with discontinuous calls and the ancestral fast single pulse pattern used the same recognition mechanism to detect pulse pattern, in addition to a verse pattern recognizer25,49. Surprisingly, even two species with the derived double pulse pattern (N. retusus, N. maxillosus) continue to use the ancestral mechanism to recognize the pulse pattern, and, accordingly, exhibit no preference between conspecific and ancestral (= single pulse) pattern18. Three species with double pulses had call recognition mechanisms that distinctly differed from the ancestral state and from each other17,26,40. One species with slow pulse rate (N. ensiger) had another, distinctly different pulse pattern recognizer42.

We and others1,15 found an interesting distribution of derived call traits across habitats. Up to four species occur within one habitat and can be heard at the same time; a fifth species may occur asynchronously, during a different season or time of day. Among the four synchronously signaling species is one with the ancestral call pattern, and one for each derived call state. For example, in central Missouri, N. robustus (ancestral state), N. bivocatus (double pulses), N. ensiger (slow pulse rate) and N. nebrascensis (discontinuous) signal simultaneously; N. retusus occurs and signals later in the season. In Trinidad, a similar acoustic community is composed of different species: N. punctipes (ancestral), N. triops (double pulses), N. saturatus (discontinuous), N. affinis (slow pulse rate) and N. maxillosus (different activity period, later at night). Accordingly, pairs of sibling species have similar calls, if they do not co-occur (e.g. N. palustris occurs in temperate North America, N. punctipes in Central America and Caribbean; N. melanorhinus is limited to salt marshes, where N. robustus does not occur). Groups of closely related species that occur sympatrically have diverse call types.

Replicated evolution among species communities, such as the pattern of derived call traits in Neonocephalus, is a key feature of adaptive radiations50. The pattern in Neoconocephalus suggests that the call diversity and the ability to discriminate among them limits the number of species at each locale. Despite the large quantitative differences within the derived call traits (e.g. chirp durations), two species with the same derived trait do not co-occur. The only exception that we are aware of is N. triops and N. affinis, which are both double pulsed. However their calls also differ in call structure and pulse rate (Table 1), so they do not occupy the same acoustic niche. These acoustic communities suggest that the sensory processing of the calls plays an important role in the stabilization of derived call traits.

The question of which evolutionary mechanisms drive the diversification of the communication in Neoconocephalus cannot be clearly resolved at this time. We have argued, based on physiological, ecological, and phylogenetic data, that genetic drift might play an important role in call diversification14. Data on the age of this group and the derived call traits as well as additional physiological data might inform this question in the future.

Competing Interests

The authors have declared that no competing interests exist.

]]> 0
PhyloPen: Phylogenetic Tree Browsing Using a Pen and Touch Interface Mon, 23 Nov 2015 15:30:54 +0000 Phylogenetic trees are used by researchers across multiple fields of study to display historical relationships between organisms or genes. Trees are used to examine the speciation process in evolutionary biology, to classify families of viruses in epidemiology, to demonstrate co-speciation in host and pathogen studies, and to explore genetic changes occurring during the disease process in cancer, among other applications. Due to their complexity and the amount of data they present in visual form, phylogenetic trees have generally been difficult to render for publication and challenging to directly interact with in digital form. To address these limitations, we developed PhyloPen, an experimental novel multi-touch and pen application that renders a phylogenetic tree and allows users to interactively navigate within the tree, examining nodes, branches, and auxiliary information, and annotate the tree for note-taking and collaboration. We present a discussion of the interactions implemented in PhyloPen and the results of a formative study that examines how the application was received after use by practicing biologists -- faculty members and graduate students in the discipline. These results are to be later used for a fully supported implementation of the software where the community will be welcomed to participate in its development.



Phylogenetic trees are a visual representation of the hierarchical relationships of common descent among a set of species or genes from a common ancestor. Often, these diagrams also include branch lengths to provide a time scale of the divergence of species or genes from their common ancestor. An example phylogenetic tree is shown in Figure 1.

Although there has been extensive work in the visualization and rendering of phylogenetic trees, relatively few systems have offered interaction with trees. Many existing tree rendering packages, such as FigTree, Mesquite, and others used by phylogenetics experts today, generate static renderings of tree diagrams (offering limited customization). Similarly, phylogenetic trees output from analysis programs (e.g. BEAST and Mr. Bayes) almost always have to be stylistically altered before presentation and/or publication, with branches colored and rotated, annotations added, and tip labels changed for readability. Without rendering flexibility, experts sometimes have to resort to manually adjusting the diagrams themselves using image editing software, such as Adobe Illustrator, which is a process that can be very time consuming. For these reasons, the goal of our work-in-progress application is to provide biology experts with an interactive experience that is efficient and hopefully more natural to use than existing tools.

We start by evaluating the impact and benefit of using pen and touch devices in browsing phylogenetic trees. It is for this purpose that we created an experimental novel pen and touch interface called PhyloPen to elicit formative feedback on this subject, as part of a larger NSF project called Arbor Workflows, a research effort to develop powerful, yet easy-to-use phylogenetic analysis tools for evolutionary biologists 1,2. In this paper, we discuss our preliminary work and findings on using pen and touch and provide a roadmap for our future development and user studies of the software on a full implementation scale. A unique contribution of this work is that it brings together the technologies of pen, touch, and annotation in a single system.


Fig. 1: An example phylogenetic tree of the Geospiza bird species from the Galapagos Islands rendered as a cladogram.

Related Work

Although many of the commonly used phylogenetic tree viewers used by experts today generate only static renderings, OneZoom 3 is one exception that is dynamic and interactive although it does not explore non-traditional modes of input. OneZoom uses a fractal layout to handle trees, giving the user the option of changing the tree structure between different types of fractals and allowing the user to pan and zoom with a mouse through the tree to points of interest, refocusing its rendering on that area to provide more details about a particular region as the user zooms closer. The fractal layout of OneZoom gives it the additional benefit of being able to handle very large phylogenetic trees.

As for using different types of input from the traditional mouse and keyboard, there have been a number of previous efforts for the visualization of data using pen-, touch-, and other input-based interactive user interfaces, some of which were recently outlined by Isenberg and Isenberg 4 . There have been a few touch-based interfaces involving phylogenetic trees and species data specifically for museum displays, such as Involv 5 and the DeepTree exhibit 6 , as well as for collaborative learning, such as Phylo-Genie 7 . In regards to projects more closely aligning with PhyloPen’s goals, a pen and touch interface has also been implemented that is oriented toward phylogeny experts for co-located collaborative comparison tasks of different phylogenetic trees 8 . Although the majority of the focus was on touch-oriented features, the pen was especially useful for leaving sticky note annotations and integrated, free-floating annotations across the collaborative display. We are not aware of any other such interfaces that study the use of pen for phylogenetic tree browsing.

Unlike other interfaces, instead of using only touch for phylogenetic tree browsing, PhyloPen uses pen and touch together. Prior studies have shown potential benefit by using the two together. In a study by Brandl et al. 9 , the results of experiments between touch only, pen only, and pen and touch together suggested that the combination of pen and touch is superior in terms of speed and accuracy and is generally preferred by users. Additionally, in a study by Walny et al. 10 , the division of labor between pen and touch was also determined to be beneficial with touch being more natural to users for physical metaphor gestures (e.g., resizing, translating, etc.) and pen being better for gestures without corresponding physical actions in the real world or for writing (e.g., annotation) and drawing. The study also demonstrated that explicit mode switching can be minimized by using the two inputs together, potentially reducing the time to perform actions.


Our pen and touch user interface, PhyloPen, is written in C# using the Windows Presentation Foundation (WPF), which is a vector-based toolkit for creating user interfaces in Windows applications. WPF was chosen for its extensive and readily-available pen and touch API features. A full demonstration video is available online at the following address:

Graphical User Interface (GUI)

The main components of our graphical user interface are the viewport to the left that contains the rendered phylogenetic tree and a resizable sidebar to the right, which contains various options and features related to the rendering. On the bottom, there is a status bar to specify the current context and status information regarding interactions with the tree rendering. The general layout of the user interface is illustrated in Figure 2. To operate on a specific phylogenetic tree, the tree is first loaded from an input file in standard format (e.g., Newick, phyloXML) into a NoSQL database (MongoDB). A list of the trees loaded into this database is then presented to the PhyloPen user for viewing and manipulation.


Fig. 2: The general layout of PhyloPen with the view of the phylogenetic tree to the left, an options sidebar to the right, and a status bar on the bottom.

In the PhyloPen user interface, the main interaction features can be summarized into several categories: gesture-based navigation, annotation, layout modification, ancillary data display, and focus+context coordination. First, gesture-based navigation includes features that allow the user to translate the tree, zoom in and out of the view, and zoom to a particular region of the tree. Second, annotation includes features that allow the user to associate ink or text annotations to a clade in the tree and copy it to other clades in the tree. The third category, layout modification, consists of features where the layout or structure of the tree is altered by the user. These features include rotating nodes (e.g. swapping the order of child clades from a node), expanding and collapsing particular clades in the tree, and changing the width and height of the tree diagram. Fourth, ancillary data display includes features that involve information about the clades in the phylogenetic tree. Only minimal data have been placed into the tree rendering for the sake of experimentation and evaluation, but this feature category will be expanded upon in the future. Currently, when the user hovers over a leaf node with the pen, information is displayed about the species it represents, including the name of the species and an optional image of the species, the latter of which could benefit educational applications. Finally, the last category, focus+context coordination, contains features that are meant to help the user maintain context while interacting with the phylogenetic tree rendering. This category currently only contains a World-In-Miniature (WIM) model, a miniaturized version of the overall phylogenetic tree, which has a moving rectangular outline that shows where the user’s viewport is relative to the entire phylogenetic tree diagram. Additionally, when the user touches the WIM model with a stylus or finger, the view will move its center to the specific location that is touched. This feature is illustrated in Figure 3 in the tabbed sidebar of the user interface.


Fig. 3: The World-In-Miniature (WIM) tree can be used to find the user’s current location in the tree as well as to jump to any location in the tree.


Of all the features of PhyloPen, the gesture-based ones are the most important to our work. Our user interface implements gestures for two types of input: touch and pen. The touch-based gestures are primarily navigation-type, physical metaphor gestures while the ink gestures are primarily reserved for direct tree manipulation. For physical metaphor gestures, touch input tends to be more natural for people, as suggested by Walny et al. 10 . By separating the navigation gestures from the manipulation gestures, both can be used without explicit mode switching. However, for this version, we were initially asked by the biologists to allow for the possibility of performing ink gestures with their fingers, which is also useful for devices that do not have pens. This resulted in some explicit mode switching in our application between ink gesture mode and (touch) navigation mode. Measures were taken, however, to minimize its disruptiveness to the user, such as automatically switching back to the default ink gesture mode whenever the pen touches the screen and going into navigation mode automatically when two fingers are pressed on the screen at once. To exit navigation mode and return to the default ink gesture mode using only one’s fingers, the user can double-tap the screen as an alternative to a pen touch.

Since our application allows drawing through the touch interface, we have two different explicit modes with different sets of gestures, all of which are summarized in Table 1. In navigation mode, the user is able to translate and scale (zoom in and out of) the phylogenetic tree diagram using the typical touch gestures for these functions: swipe to translate and pinch to zoom. However, in ink gesture mode, there is a more general purpose set of gestures. This set currently includes gestures for the following actions: zoom to a region of the diagram, rotate children within a clade, expand and collapse clades, view information on leaf nodes, and annotate. To zoom to a certain region of the tree diagram, the user draws a rectangle around the desired area in one stroke. The tree will then gradually translate and zoom to that region, becoming the new viewport. To rotate a clade’s child nodes (e.g., reverse the clade’s child ordering), the user draws a C-like symbol over the branches of the clade to which it applies. This feature can be especially useful as a teaching tool to demonstrate the arbitrariness of clade ordering. To expand or collapse a node, the user draws a loop around the node. To temporarily display information on a node, the user must hover over the node. To keep this information visible, the user double-taps the desired node. Also, to start annotating a node, the user simply draws an X symbol anywhere on the ink canvas.


Fig. 4: An example of the interface demonstrating the process of creating an annotation. The above annotation was just written and has yet to be assigned to a node.

Once the user signals the start of an annotation, the interface enters annotation mode, a distinct third mode. The user may write annotations using ink anywhere on the ink canvas (as shown in Figure 4) and then double-tap a node to associate with the annotation. Upon associating the ink annotation, a small color-coded annotation square appears next to the node (in the place of the node circle), which displays the annotation upon hover. Each node can have multiple annotations. Each new one is partially stacked over the older ones. An annotation can be held open, just as with node information, by double-tapping the square instead of hovering. If there are multiple annotations, a special annotation selection menu appears to give the user a choice between each color-coded box to eliminate the struggle of picking the correct entry by a small margin of tapping error. This selection menu, as well as the color-coded annotation squares, are illustrated by Figure 5. Another implemented feature, requested by our test participants, allows an annotation assigned to a specific clade to be copied to its ancestors or descendants by swiping left or right, respectively, from within the annotation square.


Fig. 5: When there are multiple annotations in the stack for a specific clade and the user wants to display one, a special menu appears to make it easier for the user to select the desired annotation.

Table 1: PhyloPen has a number of different touch- and ink-based gestures to help browse the phylogenetic tree diagram.


Besides the gestures listed in Table 1, there is one additional gesture-based feature in PhyloPen. On the user interface sidebar, there are buttons for changing the width and height of the tree by a fixed amount. However, we also include a gesture-based region for resizing the tree layout, which we refer to as the tree resize canvas. It is shown as the small rectangular ink canvas in the bottom-right corner of the GUI in Figure 2. When a line is drawn on the miniature canvas, the change in width and height of the tree depends on the direction, angle, and length of the stroke. This is illustrated in Figure 6.


Fig. 6: The tree resize ink canvas can change the width and height of the tree.

User Evaluations

We conducted a formative user evaluation of PhyloPen, involving three Ph.D. students and three faculty members from one of our biological collaborators’ institutions, to gauge the interface’s usefulness and usability. Each participant individually tested PhyloPen on a Microsoft Surface Pro 2 with pen and touch support and responded to a post-questionnaire. In the post-questionnaire, users were asked to rate each feature of the user interface on a 5-point Likert scaling according to its usefulness and usability, separately. In terms of usefulness, the ratings were labeled as hardly useful, a little useful, pretty useful, really useful, and extremely useful. The usability was expressed in similar terms (labeled as “fun”). Quantifying these ratings from 1 to 5, the average rating for each feature is shown in the bar chart in Figure 7. The users were also given a variety of other questions in the post-questionnaire to provide feedback on how to direct future development.


Fig. 7: The bar chart shows the average Likert ratings from users on features of PhyloPen in terms of usefulness and perceived usability.

For the most part, according to Figure 7, the relative usefulness of a feature, as perceived by reviewers, seems correlated with the usability of the feature. Annotation has the lowest score in terms of usefulness, and one reason may be that annotations can only be stored as strokes, not recognized text. Also, many users seemed to experience some degree of frustration with the system. One potential source of frustration for the users may be having to double-tap a tiny node circle with the pen in order to associate the annotation. Although it is possible to enlarge the circles, we feel it may prove better to implement a circular gesture to surround the node for selection. This should be less concentration-intensive for the user. There are also cases where the user may want to delete the annotation from certain nodes after passing it up or down the tree (e.g., annotating the bird clades with “can fly” and then removing it from the clades of ostriches, rheas, kiwis, etc). We are not certain how best to implement a group annotation deletion gesture. Due to these limitations, we feel the annotation system may need to be revamped to handle a more flexible approach that is more comfortable for the user.

Two other features that scored relatively low in terms of usability were node rotation and the tree resize canvas. Currently, the two children flip instantaneously after the rotate gesture is performed. However, when this flip is performed over a large area of the tree diagram, an instantaneous change results in context loss and potential confusion of the user. We hypothesize that there would be higher scores for rotating clade children if we animated the rotation process. As for the tree resize canvas, we discovered that there was some confusion between the coordinate system of the tree resize feature and the coordinate system of the main canvas. Our tree resize canvas feature uses a standard Euclidean coordinate system, but the graphical components in the main canvas are stored in a coordinate system with an inverted y-axis. Thus, when the tree resizes in height, the tree grows downward. However, we have the users drawing up in the direction opposite of growth, making the feature far less intuitive than if the tree grew in the direction the user draws.

One perplexing result was that the translate/zoom feature was rated higher in usefulness than in usability. However, we believe that this could be the result of issues we experienced with touch features using WPF where the transformations were not as smooth as would be expected or the touch input itself would occasionally stop being received by the program with seemingly no explanation.

Overall, the feedback from the biologists was positive despite a few points of frustration. The users were generally able to figure out how to use the interface after watching a short demonstration (although one individual expressed the desire for assistance with remembering the gestures, such as a short list/cheat sheet of all the ink gestures). The biologists found the touch navigation relatively natural and intuitive while other features still needed improvement, as noted earlier. Additional features were also suggested, the most common of which was to move a clade within a tree by re-attaching it to an alternative ancestor.

Moving Forward

From the start, we worked on this version of PhyloPen as an experimental build. As such, we did not necessarily intend to make it available to the general public and to support it. However, this paper serves as an invitation to the biological community to work with us, giving feedback in creating the next generation of PhyloPen, which will be made available and open source. In this way, we will not be working in isolation, separated from the community that we wish to serve, because it is the community that best knows what it wants and needs. The new version will account for the feedback given by the users already contributing to the evaluation and will be open to more changes based on what the community finds most appropriate. When ready, a new version of PhyloPen will be released on the Arbor Workflows GitHub ( For each version released, the users will be able to test the software and fill out a post-questionnaire, offering feedback on their overall experience, how to improve what is available, and what features are highest in priority that have not already been implemented. The software resulting from this development process will ultimately be given a summative evaluation with a larger group of people than the experimental version.


Our objective is to develop an interactive phylogenetic tree-browsing user interface for experts that is more dynamic and ideally more natural than commonly used browsers, ultimately allowing the user to navigate, annotate, and change the tree structure with ease. In this work-in-progress, we tested an experimental pen- and touch-based interactive browser we developed called PhyloPen. Our preliminary formative evaluations show positive results, demonstrating some usefulness of pen and touch interfaces for phylogenetic trees. However, additional work based on our formative evaluations, including the revision of existing features and the addition of new ones (e.g., cut and reattachment for clades), and a thorough summative evaluation will be needed to draw more definitive conclusions. In support of this quest, we invite the biological community to work alongside us in reaching that goal, giving their opinions and shaping what PhyloPen ultimately should be.

Competing Interests

The authors have declared that no competing interests exist.

]]> 0
EAPhy: A Flexible Tool for High-throughput Quality Filtering of Exon-alignments and Data Processing for Phylogenetic Methods Wed, 05 Aug 2015 14:45:54 +0000 Recently developed molecular methods enable geneticists to target and sequence thousands of orthologous loci and infer evolutionary relationships across the tree of life. Large numbers of genetic markers benefit species tree inference but visual inspection of alignment quality, as traditionally conducted, is challenging with thousands of loci. Furthermore, due to the impracticality of repeated visual inspection with alternative filtering criteria, the potential consequences of using datasets with different degrees of missing data remain nominally explored in most empirical phylogenomic studies. In this short communication, I describe a flexible high-throughput pipeline designed to assess alignment quality and filter exonic sequence data for subsequent inference. The stringency criteria for alignment quality and missing data can be adapted based on the expected level of sequence divergence. Each alignment is automatically evaluated based on the stringency criteria specified, significantly reducing the number of alignments that require visual inspection. By developing a rapid method for alignment filtering and quality assessment, the consistency of phylogenetic estimation based on exonic sequence alignments can be further explored across distinct inference methods, while accounting for different degrees of missing data.



High-Throughput Sequencing (HTS) has revolutionised the field of phylogenetics by enabling researchers to question the evolutionary relationships between taxa with large-scale multi-locus datasets 1,2. The development of these methods has been driven by a realisation that the inclusion of many genetic markers helps to account for stochastic coalescent histories of individual genes 3,4,5,6. Species tree inference methods use the multispecies coalescent model to estimate potential gene tree – species tree discordance and large numbers of unlinked loci represent a greater sample of the gene tree distribution underlying the true species tree 6. However, while phylogenetic estimation might improve by sequencing many loci 4,5,6,7, the requirement for high-quality sequence alignments remains unchanged and is fundamental for the correct inference of phylogenetic hypotheses. Existing alignment methods can be extrapolated for use with large-scale multi-locus datasets, but visual inspection of each alignment, the traditional approach for assessing alignment quality, is challenging with thousands of sequenced loci 8. As a consequence of the impracticality of visual inspection, the impact of missing data in large phylogenomic datasets is often nominally explored and the potential consequences of distinct alignment filtering criteria remain unknown. Nonetheless, contradicting opinions coexist 9,10,11 regarding the effect of missing data on phylogenetic inference and it is therefore advisable to quantify the sensitivity of empirical phylogenetic hypotheses to data filtering choices. Thus we need workflows that automate (as far as possible) the assessment of alignment quality and the consequences (in terms of missing data) of making different choices about filtering criteria. Ideally, such a workflow would facilitate the conversion of individual contiguous sequences (‘contigs’) into quality-filtered alignments, and help to minimise the demand for visual inspection.

The need for a high-throughput alignment filtering system emerged with the recent advance in molecular methods to target and sequence large numbers of orthologous loci. Since whole-genome sequencing is still too costly for most research labs that focus on non-model organisms, genome reduction protocols have been developed that isolate large numbers of orthologous loci across the genome of closely related and deeply divergent taxa 1,8. There are two increasingly popular genome reduction methods that specifically focus on exonic sequence regions and can generate genetic markers suitable for phylogenomic inference. Transcriptome sequencing is a cost-effective method that does not require the a-priori availability of genomic resources. RNA is extracted from the same tissue in different target species and with the expected expression of similar genes, orthologous loci are isolated and sequenced for phylogenetic comparison. An alternative method, exon-capture 12, is a target enrichment approach 13,14,15,16 that benefits from an increasing number of readily available genomic resources and enables the design of study-specific capture systems. The use of exonic sequence regions for phylogenetic inference, generated by transcriptome sequencing or exon-capture, is promising and has been successfully demonstrated at different levels of divergence across the tree of life 17,18,19,20. The tremendous increase in the scale of available exonic loci benefit inference methods, but also requires a significant investment in the development of bioinformatic resources to process such data.

Whereas several excellent bioinformatic pipelines have been constructed for processing raw sequence data and conducting sequence assembly 13,14,15,16,19, a bioinformatic scheme is needed for subsequent alignment, alignment quality assessment and alignment filtering. Most published studies still conduct visual inspection of alignment quality and account for missing data by dividing datasets into a limited number of categories manually (i.e. 14,18)or automated (i.e. 21). Recently, Misof et al.19 developed a method to assess alignment quality in an extensive study that used transcriptomes to infer the phylogeny of insects. They identified potentially erroneous alignments by calculating the BLOSUM62 distance between each amino acid sequence and the best reciprocal hit of a reference taxon. A distance calculation based on a BLOSUM matrix was warranted, due to a significant level of protein divergence between most taxa. The BLOSUM alignment score matrix values the alignment of each amino acid pair differently, representing the likelihood of amino acid substitutions, but lacks resolution when the expected level of protein divergence between two sequences is limited. However, although this has not been tested prior, it can be expected that at shallower levels of divergence subtle misalignments might actually have more significant consequences for phylogenetic estimation than when inferring relationships between distant taxa, stressing the need to identify such misalignments. When assessing alignments with limited levels of sequence divergence, the exact number of clustered amino acid changes is more likely a better indicator of alignment quality than the overall BLOSUM62 distance score. In this short communication, I describe a flexible high-throughput pipeline for quality assessment of exonic sequence alignments and subsequent filtering of missing data. The pipeline is specifically designed to be flexible and process both population and phylogenetic level data, but the method developed by Misof et al.19 will likely be more effective at deep phylogenetic scales.

EAPhy, exon alignment for phylogenetics, was developed to process exonic sequence data for phylogenetic inference, but is valuable for any type of analysis that requires high-quality filtered alignments (i.e. population genomics or molecular evolution). In this manuscript I will focus on its application for phylogenetic inference. The first objective of the pipeline is to quantify alignment quality and highlight just those loci that require visual inspection. By translating exonic nucleotide alignments into amino acid alignments, EAPhy infers the relative quality of sequences and alignments by assuming that most mutations within exons are silent. In addition, the identification of regions that harbor an excessive cluster of amino acid replacements distinct from a summary reference sequence, is used as a proxy for alignment quality. Simultaneously, insertions and deletions that result in frame shifts and the introduction of multiple stop-codons are unlikely to represent true biological events and such alignments should be addressed. The pipeline can be adapted based on the expected level of divergence between taxa by adjusting the stringency of filtering criteria. The second objective is to provide a user-friendly method to account for missing data. By enabling filtering criteria for missing data to vary, the consistency of phylogenetic estimation can be quantified across different levels of missing data. Lastly, EAPhy was designed to generate alignments of different sorts (haplotype, diplotype and SNP based), in the formats required for most commonly used inference software and facilitate the further exploration of distinct analysis methods. With the development of a high-throughput method for alignment filtering and processing, the overarching aim of this pipeline is to reduce the bioinformatic burden of data analysis involving exonic sequence alignments and ultimately promote further research into the (in)congruences between inference methods, while accounting for different degrees of missing data.

Overview of Methods

EAPhy consists of a collection of scripts that takes as input a set of unaligned sequences for an arbitrary number of species and loci. It will generate multiple sequence alignments using existing aligning software and subsequently filters these alignments for a number of user-specified criteria. The final output consists of quality filtered multiple sequence alignments, allowing different degrees of missing data as preferred, and a list of alignments that still require visual inspection. The output files are automatically exported in the input format of commonly used phylogenetic inference programs. The complete package is freely available at

EAPhy can be run on most individual computers (i.e. does not require a cluster set-up) and an individual run for a modest dataset can be completed within hours. The pipeline has been used with exon-capture datasets involving tens to hundreds of individuals and thousands of loci, and finished within six hours on a Macintosh desktop computer with a 3.1 GHz Intel i7 processor (2012) and 16 GB of RAM. It is important for the user to adopt filtering criteria suitable for the dataset (level of divergence and data quality) analyzed, but if filtering criteria have been carefully reviewed EAPhy should be able to handle larger datasets than currently tested. For EAPhy to function appropriately, I advise to run the pipeline initially with a small subset of the data and replicate the analysis with alternative filtering criteria. If filtering and flagging of alignments works well, then the analysis can be extrapolated for usage with the complete dataset. The importance of specifying appropriate filtering criteria should not be underestimated, since misspecification of filtering criteria will result in a significantly reduced dataset or alternatively a dataset that equals the input data, regardless of potential low-quality alignments.

EAPhy is not designed to identify individual sequencing errors that are often associated with HTS datasets, but will identify sequence regions with excessive non-synonymous substitutions (potential ‘low-coverage’ sequences) if these have not been filtered out beforehand and appear anomalous in the resulting alignments. Several excellent pipelines have been developed to filter raw sequence data and generate assemblies, and the starting point of this pipeline requires assembled individual contigs that start in first codon frame, for each presumed orthologous locus. A complete overview of the pipeline is outlined in Figure 1 and a general description of the most important components is provided here.

Specification of configuration script

At the onset of each EAPhy run, the system path to an align program executable and all filtering criteria for downstream analysis are specified in a single configuration script. Muscle 22 is the default aligner used by EAPhy, but should be installed by the user independently from downloading EAPhy. Alternative alignment software can be used but requires modifications of several scripts. The EAPhy pipeline is designed as a set of modules that can be executed independently or in consecutive order as a complete analysis (Fig. 1). This provides a straightforward system to reiterate specific components of the pipeline, with alternative filtering criteria for alignment quality or missing data. A complete description of all filtering parameters can be found in the manual and is part of the EAPhy package that can be downloaded from GitHub.


Fig. 1: A schematic overview of the EAPhy workflow

The user specifies filtering instructions in a single configuration script (1). EAPhy will first subset the target individuals from each locus and create new contig files (2). With an existing aligner, new alignments are created for each locus (3) and are subsequently processed and checked (4). Alignments are highlighted that do not fulfill the filtering criteria and can be visually inspected. If deemed appropriate for inclusion, manually checked alignments can be added to the filtered alignment list (5). Alternatively, EAPhy can automatically continue with the alignments that passed filtering and exclude the problematic loci. Once the complete collection of filtered loci has been identified, final alignments are generated (6). If diplotype sequence data was used and heterozygous positions coded according to IUPAC format, concatenated (7) and SNP (8) alignments can be generated if required.

Missing data – within alignments

The effect of missing data on phylogenetic inference is not well understood and contradicting opinions coexist 9,10,11. Phylogenetic estimation is likely unbiased with large numbers of loci, if there are no systematic differences in sequence length between individuals for any given locus. However, the maximum-likelihood (ML) estimation might cluster individuals by sequence length rather than sequence similarity for complete positions, when missing data are non-randomly distributed and specific individuals have systematically shorter contigs or are completely missing for specific loci. At sites with missing data, the probability of observing an ‘A’, ‘T’, ‘C’ or ‘G’, is set to 1 and ML will group taxa together for which there is more signal and less uncertainty (Stamatakis, pers. comm.). Thus the effect of missing data is not limited to small sequence datasets but should also be accounted for and characterized in large-scale datasets. With the development of EAPhy, I do not advocate to discard or include incomplete sites but rather provide the opportunity to account for missing data by generating datasets where different filtering criteria have been enforced.

Missing data within individual sequences are particularly prevalent at the beginning and end of alignments (‘jagged edges’), since individual contig sequences often differ in length. Once alignments have been constructed using an existing aligner (e.g. Muscle 22), EAPhy will first address missing data by processing alignments in accordance with stringency criteria specified in the configuration script (Fig. 2). First, potential gaps within individual sequences are removed to yield long consecutive sequences (Fig. 2.1). EAPhy converts all sequence alignments into amino acid alignments and then uses a ‘jump-sliding window’ approach to assess the presence of potential non-consecutive sequence stretches that are often prevalent at the start/end of individual sequences. A jump-sliding window approach was developed since a conventional sliding window approach would remove the complete individual sequence if the first window would contain more missing data than allowed. Each window is assessed on the presence of amino acid sequence gaps and if a window contains more gaps than allowed, the complete window is removed for that individual sequence. In-frame gaps (i.e. triplet insertions) are retained if the amount of inserted codon gaps per window does not contain more missing codons than allowed. The window then ‘jumps’ a sequence distance of half the window size plus one codon and the process reiterates. By converting a nucleotide alignment in its amino acid equivalent, EAPhy specifically takes into account the coding-codon character of exonic sequences. When nucleotide sequence data is removed by codon, the remaining sequence is still in correct frame and codon position can still be inferred for each nucleotide position.

After individual sequences have been trimmed for missing data, EAPhy then assesses missing data between individuals by evaluating the amount of missing data for each amino acid alignment column (Fig. 2.2). The algorithm used is similar to the jump-sliding window approach, but now focuses on the amount of missing data within each amino acid alignment column. The window-length of amino acid columns and the amount of missing data allowed within each column, can be specified in the configuration script. The algorithm evaluates for each amino acid column whether the amount of individuals with missing data exceeds the cut-off specified. If more than half of the columns in a given window have more missing data than allowed, the columns in the first half of the window are removed from the alignment. The window then ‘jumps’ a sequence distance of half the window size plus one codon and the process reiterates. Amino acid columns at the end of alignments are removed, if they have not been evaluated but the specified window length exceeds the number of remaining columns. When alignments have been filtered for missing data within and between individuals, EAPhy evaluates the presence of single nucleotide insertions, by assessing the frequency of sequenced individuals for each nucleotide alignment column. If the number of sequenced individuals is below a user specified cut-off, the site is assumed to be a sequencing error and removed from the alignment.


Fig. 2: An exemplary overview of the three main filtering steps conducted during alignment filtering and quality assessment

First, potential gaps within individual sequences are removed to yield long consecutive sequences (1). By converting nucleotide codons into amino-acids, the number of missing amino acids is assessed for each window using a ‘jump-sliding-window’ approach. If less amino-acids are missing for a given window than the specified ‘jump-window gap ratio’, the complete corresponding nucleotide stretch is retained (see green frames). If more amino acids are missing for a window, the complete nucleotide stretch is removed (see red frame). Secondly, the amount of missing data for each amino-acid alignment column in a given window is quantified (2). If more than half of the amino acid columns in a window miss more individuals than the specified ‘column gap ratio’, all corresponding nucleotide columns are removed (see red frame). Lastly, the quality of the resulting alignment is assessed, by comparing each individual sequence to a consensus sequence (3). Following a common sliding window approach for each individual sequence, the number of amino acids identical to the consensus is quantified. If, for each window, the number of amino acids distinct is less than the specified ‘difference ratio’, the alignment is retained (see green frames). If for any individual within an alignment, a window would fail this criterion, the alignment is flagged for visual inspection. In addition, for each sequence the number of stop-codons is quantified and if any individual sequence contains more stop-codons than a specified cut-off number (e.g. > 1), the alignment is also flagged for visual inspection.

Alignment quality

Once each alignment has been filtered for missing data, EAPhy then inspects the alignment quality by translating the nucleotide sequences and evaluating the resulting amino acid alignment (Fig. 2.3). First, if the number of stop codons for any individual sequence exceeds a user specified cut-off value (e.g. > 1), the alignment is flagged for visual inspection. Subsequently, a general consensus sequence is estimated for each alignment, and each individual sequence is compared to the consensus sequence in a ‘normal’ sliding-window approach. The window length is specified in the configuration script and each individual sequence is compared to the consensus sequence by sliding window. For each window, the number of amino acids distinct from the consensus is quantified and if greater than the proportion specified in the configuration script, the alignment is flagged for visual inspection.

Finally, phylogenetic inference is dependent on the comparison of orthologous genetic markers and comparing potential paralogous loci might yield confounded estimates of relationship. EAPhy assumes that the sequenced contigs for each locus are orthologous but has an additional option to potentially identify paralogous loci, by identifying markers with excessive levels of average individual heterozygosity. The user can inspect the distribution of average individual heterozygosity across all loci and based on this observation make an informed decision whether to exclude a certain percentage of loci with the highest level of average individual heterozygosity.

Concatenation and SNP selection

After visual inspection and filtering of flagged alignments, the collection of final high quality alignments can then be used for a variety of phylogenetic estimation methods. Gene trees can be inferred based on single alignments and a concatenated maximum likelihood tree can be estimated based on all alignments combined. Since all alignment filtering was conducted by codon, each nucleotide can still be assigned its correct codon position. PartitionFinder 23estimates the most optimal partitioning scheme across all sequence positions and appropriate substitution model for each partition. A PartitionFinder input file is automatically created with each gene and codon position of the concatenated alignment specified.

In addition to sequence-based alignments, EAPhy will also generate concatenated alignments that include polymorphic sites exclusively. SNAPP 24 is a species tree method that uses unlinked biallelic markers, instead of sequence-based alignments, and EAPhy can generate alignments with a biallelic SNP randomly sampled from each locus. It will verify whether polymorphic sites are biallelic and neglect polymorphic sites with more than two allelic states. Alternatively SNP alignments can be constructed where every single SNP is considered, regardless of allele count, or with all SNP’s across all loci concatenated. If a study is geared towards recovering population structure, such alignments can be used in analyses that model allele frequencies (e.g. 25).

Missing data – number of sequenced individuals

Sequencing success can vary among individual samples. If specific individuals are systematically underrepresented and miss data for many loci, it is possible that the phylogenetic placement of such taxa is ambiguous and the investigator would prefer to exclude these samples. Thus, the potential impact of missing individuals across loci should be accounted for. EAPhy attempts to highlight where this is likely by: a) providing alternative datasets with different numbers of missing individuals allowed and b) providing summary statistic output files quantifying the number of loci sequenced for each individual. This enables the investigator to further explore the potential effects of missing data on phylogenetic inference.

In Summary

The first objective of developing EAPhy was to provide a flexible and rigorous tool to generate reliable alignments, while minimizing the need for extensive visual inspection. Secondly, EAPhy was designed to allow filtering criteria for missing data to vary and investigate the impact of missing data on phylogenetic estimation. Lastly, EAPhy creates a large number of desired input formats for subsequent analysis, enabling the exploration of distinct inference methods. Negating the effort of manual alignment filtering and processing, EAPhy will hopefully stimulate further research into the potential consequences of applying alternative criteria for missing data and datatype, and how this might ultimately result in (in)congruent estimates of phylogenetic relationships across methods. The simultaneous development of novel molecular approaches to sequence orthologous genetic markers and bioinformatic methods to analyze such data, will ultimately provide us with the tools to generate a phylogenetic framework for all taxa across the tree of life.

Competing Interests

The author has declared that no competing interests exist

]]> 0
Visualising Geophylogenies in Web Maps Using GeoJSON Tue, 23 Jun 2015 13:00:40 +0000 This article describes a simple tool to display geophylogenies on web maps including Google Maps and OpenStreetMap. The tool reads a NEXUS format file that includes geographic information, and outputs a GeoJSON format file that can be displayed in a web map application.



The increasing number of georeferenced sequences in GenBank 1 and the growth of DNA barcoding 2 means that the raw material to create geophylogenies 3 is readily available. However, constructing visualisations of phylogenies and geography together can be tedious. Several early efforts at visualising geophylogenies focussed on using existing GIS software 4, or tools such as Google Earth 5,6,7 . While the 3D visualisations enabled by Google Earth are engaging, it’s not clear that they are easy to interpret. Another tool, GenGIS 12,13 , supports 2D visualisations where the phylogeny is drawn flat on the map, avoiding some of the problems of Google Earth visualisations. However, like Google Earth, GenGIS requires the user to download and install additional software on their computer.

By comparison, web maps such as Google Maps 15 and OpenStreetMap 16 are becoming ubiquitous and work in most modern web browsers. They support displaying user-supplied data, including geometrical information encoded in formats such as GeoJSON, making them a light weight alternative to 3D geophylogeny viewers. This paper describes a tool that makes use of the GeoJSON format and the capabilities of web maps to create quick and simple visualisations of geophylogenies.

2D layout of geophylogenies

The following discussion assumes that we have a phylogeny, and that for most (if not all) of the OTUs in that phylogeny are associated with a point locality for which we know the latitude and longitude.

In order to draw a geophylogeny on a web map we need to solve three problems. The first, relatively trivial problem is to place the the localities of the OTUs on the map (I shall refer to these as “occurrences”).

The second is to draw the phylogeny. Typically when drawing an evolutionary tree we compute x and y coordinates for a device where these coordinates have equal units and are linear in both horizontal and vertical dimensions, such as a computer screen or printer. In web maps coordinates are expressed in terms of latitude and longitude, and in the widely-used “web mercator” projection the y-axis (latitude) is non-linear. Furthermore, on a web map the user can zoom in and out, so pixel-based coordinates only make sense with respect to a particular zoom level.

The single 256 × 256 pixel tile representing the globe a zoom level 0 showing the pixel coordinates for the top left corner, the centre (corresponding to longitude 0, latitude 0), and the bottom right. Tile image map tiles by CartoDB under CC-BY 3.0 license.

Fig. 1: Web map tile


Web maps use “tiles” of a fixed size to represent the globe. Each tile is typically 256 × 256 pixels in size, and the number of tiles comprising a map is 2zoom where zoom is the zoom level. At zoom level 0 the map comprises a single tile (Fig. 1), at zoom level 1 the map comprises 4 tiles, and zoom level 2 eight tiles, and so on. To accommodate the web mercator projection, we first compute a geographic bounding box for the tree based on the bounding box that encloses the occurrences, then offset that box so that so that the tree is drawn, say below, the occurrences. We can then convert the longitude λ and latitude Ω coordinates of the bounding box to pixels x and y at zoom level 0 using the following formulae:

Note that the maximum latitude that can be displayed in the web mercator projection is 85.051129° north and south. The tree drawing is then laid out within that bounding box, with the nodes positioned in terms of pixels. Once pixel coordinates have been computed for the whole tree, they are then converted back to latitude and longitude values:

Expressing the tree in terms of latitude and longitude coordinates means that the rendering of the tree as the user zooms in and out is handled automatically by the web map application.

If we want to provide the user with a visual connection between each occurrence on the map and the location of the corresponding OTU in the phylogeny, we can draw a line connecting the two. These lines may criss-cross creating visual clutter, reducing this clutter is the third problem. To make the diagram more comprehensible, I adopt the approach used by GenGIS 12,13 to reorder the nodes in the tree to minimise the number of crossings 8. As an additional feature, if a taxon is represented by more than one occurrence, we can enclose the set of occurrences by a convex polygon to represent the range of that taxon.

Having computed a layout, we then need to render that on a web map. There are a number of different web maps available, each with their own API. Rather than tie the visualisation to a particular API, we can use a standardised output format, such as GeoJSON, to encode the layout, so that users can pick which web map they wish to use for the visualisation.


GeoJSON 17 is a format for encoding geographic data in JSON (JavaScript Object Notation). It includes various geometry types (such as Point, LineString, and Polygon), and is supported by a number of online mapping tools, including Google Maps 15 and Leaflet 18 . A GeoJSON document comprises a set of one or more features, each of which has a geometry and additional properties. Using the GeoJSON geometry types we can encode occurrences (Point), the tree (a set of LineString), and taxon distributions (Polygon) in GeoJSON, then have the entire visualisation rendered by the web application. The GeoJSON specification does not, by itself, include any information on how to display the objects encoded in a GeoJSON document (e.g., what colour to use for a line), but some informal standards have emerged, such as storing CSS styles as properties.

Input format

In order to create the visualisation we also need a way to input a phylogeny and the geographic localities. The approach taken here is to use the NEXUS format 9 , and the GEOGRAPHIC datatype introduced by the Mesquite Cartographer package 14 . While some might argue that XML represents the future of phylogenetic file formats 10 , NEXUS is easy to manually edit and hence facilitates debugging and exploring the software. Given a set of OTUs, the tool expects a NEXUS file with a TREES block describing a tree, followed by a CHARACTERS block encoding the location of each OTU. Each OTU is typically a DNA sequence. Sets of sequence may belong to the same taxon (e.g., a species or a DNA barcode BIN 2 ). Following Mesquite, this information can be stored in an ALTTAXNAMES command in a NOTES block.

Figure 2 shows a NEXUS file for the widely used Banza example 11,19

NEXUS file for Hawaiian Banza, with geographical data encoded in the CHARACTERS block.

Fig. 2: NEXUS file for Hawaiian Banza



I have implemented a NEXUS to GeoJSON converter using PHP. The code parses the NEXUS file, computes a bounding box based on the distribution of the OTUs, draws the tree, and exports the result in GeoJSON. The code is available on github Code for the examples in this article are available from A live demo can be explored at which includes examples of visualising geophylogenies using both Google Maps (Fig. 3) and Leaflet (Fig. 4).

Fig3-GoogleMaps-CC-BY-no logo

Fig. 3: Geophylogeny for South American marsupial

Geophylogeny for DNA barcodes for the marsupial Proechimys guyannensis, showing two distinct clusters that are geographically allopatric (data from BOLD, map tiles by CartoDB under CC-BY 3.0 license).

Geophylogeny for Hawaiian katydids (genus Banza) displayed using the Leaflet framework with map tiles by CartoDB under CC-BY 3.0 license.

Fig. 4: Geophylogeny for Hawaiian katydids



At present the method described here requires a middle layer (written in PHP) that resides on a web server and converts the NEXUS file to GeoJSON. An obvious extension would be to port that code to Javascript and have the entire tool function within the web-browser client.

Although lacking some of the functionality of more specialised software such as GenGIS, an advantage of a web map-based tool is that it brings phylogenies into an environment already familiar to users of biodiversity data, such as the GBIF portal. Many users will have already encountered points on maps, and layers (e.g., of environmental data, or estimated species distributions). By representing phylogeny in GeoJSON we open the way for phylogenetic information to be incorporated into these maps.

Another reason GeoJSON is attractive is that because it is a JSON document it could be stored and indexed in a document database such as CouchDB 20 , which I’ve used elsewhere for taxonomic and phylogenetic data 21 . Hence we could imagine being able to quickly build a database of geophylogenies that can be queried both taxonomically and spatially. This would be one way to tackle the challenge of Kidd’s call for a “map of life”3.

Competing interests

The authors have declared that no competing interests exist.

]]> 0
Concatenation Analyses in the Presence of Incomplete Lineage Sorting Fri, 22 May 2015 16:00:52 +0000 Incomplete lineage sorting (ILS), modelled by the multi-species coalescent, is a process that results in a gene tree being different from the species tree. Because ILS is expected to occur for at least some loci within genome-scale analyses, the evaluation of species tree estimation methods in the presence of ILS is of great interest. Performance on simulated and biological data have suggested that concatenation analyses can result in the wrong tree with high support under some conditions, and a recent theoretical result by Roch and Steel proved that concatenation using unpartitioned maximum likelihood analysis can be statistically inconsistent in the presence of ILS. In this study, we survey the major species tree estimation methods, including the newly proposed “statistical binning” methods, and discuss their theoretical properties. We also note that there are two interpretations of the term “statistical consistency”, and discuss the theoretical results proven under both interpretations.



Estimating species trees from multiple loci is commonly performed using concatenation methods, in which multiple sequence alignments from different genomic regions are concatenated into one large supermatrix, and then a tree is estimated on the supermatrix. Yet, incomplete lineage sorting (ILS) (modelled by the multi-species coalescent model) can result in different loci having topologically different phylogenies, with high probability of gene tree incongruence when the effective population size is large and the time between speciation events is small 3. Most importantly, Roch and Steel 1 recently proved that using unpartitioned maximum likelihood to estimate a species tree on a concatenated alignment from different loci can converge to a tree other than the species tree as the number of loci increases, even if the sequence length per locus is allowed to increase; in other words, unpartitioned maximum likelihood can be statistically inconsistent (and even positively misleading, a stronger statement) under the multi-species coalescent model. Furthermore, simulation studies have shown that species trees estimated using concatenation can result in incorrect trees with high support 2. Thus, both theory and empirical studies show that the use of concatenation to estimate a species tree from multiple loci can lead to incorrect phylogenetic estimations.

Because concatenation can result in incorrect trees, coalescent-based “species tree methods” have been developed that are statistically consistent under the multi-species coalescent model. However, there are two meanings for “statistical consistency under the multi-species coalescent model”, and so it is important to distinguish between them.

  • Statistical consistency under the multi-species coalescent model – weak version. The most common use of the term asserts that the tree estimated by the species tree method will converge in probability to the true species tree as the number of sites per locus and the number of loci both increase 3.
  • Statistical consistency under the multi-species coalescent model – strong version. The other use of the term asserts that the estimated tree will converge in probability to the true species tree as the number of loci increases, but limits the sequence length per locus (perhaps to a constant number of sites). The first use of the term is clearly the weaker condition, since it makes stronger assumptions.

Thus, what Roch and Steel proved 1 is that unpartitioned maximum likelihood is statistically inconsistent in both senses. However, their proof does not extend to fully partitioned maximum likelihood analyses, which allow the numeric model parameters to change between the different loci in the concatenated alignment. Indeed, it is not yet established whether unpartitioned maximum likelihood analyses are inconsistent or consistent under either of these interpretations of the meaning of statistical consistency.

On the other hand, many species tree estimation methods have been developed that are provably statistically consistent in the first (i.e., weaker) sense of the term, using a variety of algorithmic approaches. For example, some methods (e.g., *BEAST 20) co-estimate gene trees and species trees) and others (e.g., SVDquartets 15) estimate the tree using a single site from each locus within a set of unlinked loci; however, most of the commonly used methods estimate the species tree by combining estimated gene trees. These “summary methods” (such as MP-EST 12 and ASTRAL 4) are statistically consistent under the multi-species coalescent model in the first (i.e., weaker) sense of the term, and are increasingly popular due to their relative speed and ease of use.

Only a few coalescent-based species tree estimation methods have been proven to be statistically consistent in the second (i.e., stronger) sense of the term, which establishes that the species tree estimated by the method converges to the true species tree as the number of loci is allowed to increase, even when the sequence length per locus is bounded 7,16,18,19. However, it is not known whether the commonly used summary methods such as MP-EST are statistically consistent under the second sense of the term.

Simulation studies evaluating the relative performance of species tree estimation methods in comparison to concatenation analysis have had mixed results (e.g., sometimes concatenation is more accurate, sometimes coalescent-based methods are more accurate, and sometimes the differences are not statistically significant 4,5,6,8,9,10,11). In addition, coalescent-based summary methods have been shown to have reduced accuracy in the presence of substantial gene tree estimation error 4,5,6,8,10,11. Furthermore, the proofs of statistical consistency for standard coalescent-based summary methods (such as MP-EST and ASTRAL) have assumed that all the input gene trees are estimated without any error 7. Since estimated gene trees are often imperfect, the assumption of 100% accurate gene trees is unlikely to be biologically realistic 8.

Statistical Binning

In two recent papers 10,11, we presented techniques designed to improve coalescent-based species tree estimation by improving gene tree estimation when estimated gene trees have insufficient accuracy as a result of limited phylogenetic signal (e.g., because of low rates of evolution or sequence lengths that are too short). These two techniques partition the loci into sets, so that each set should contain loci that are deemed likely to have a common evolutionary tree. The first of these techniques 11 is called “unweighted statistical binning” and the second 10 is called “weighted statistical binning”; the use of weighting in statistical binning addresses a weakness we identified in unweighted statistical binning. Each type of statistical binning uses a fully partitioned maximum likelihood analysis to estimate a tree on each bin, and then these “supergene” trees are given to a summary method (such as MP-EST) to compute a species tree. Hence, statistical binning is a combination of concatenation and coalescent-based species tree estimation.

Statistical binning has the following steps:

  1. Compute gene trees with bootstrapping on each locus, using maximum likelihood.
  2. Use the bootstrap support on the edges of each gene tree to determine for each pair of genes whether they are likely to have a common gene tree topology, and build an “incompatibility graph” to represent this information (so that each node represents a gene, and two genes are connected by an edge if the topological differences between their gene trees are considered statistically significant according to the test).
  3. Partition the vertices of the graph into sets of approximately the same size, so that no two vertices in any set are adjacent; these are called the “bins”.
  4. For each bin, concatenate the alignments in the bin, and compute a fully partitioned maximum likelihood tree on the bin; these are called the “supergene trees”. If performing a weighted statistical binning, then replicate the supergene tree by the number of genes in its bin. If performing an unweighted statistical binning, then do not replicate the supergene trees.
  5. Apply the preferred summary method (e.g., ASTRAL or MP-EST) to the set of supergene trees to compute the species tree.

Note that the only difference between weighted and unweighted statistical binning is that weighted statistical binning replicates the supergene trees by the number of genes in the bin, and unweighted statistical binning does not do this. This difference is essential to the theoretical properties of the two methods. As proven in 10, pipelines based on weighted statistical binning followed by summary methods such as MP-EST or ASTRAL (that will converge to the true species tree as the number of true gene trees increases) are statistically consistent using the first definition, which allows the number of sites per locus as well as the number of loci to increase. However, unweighted statistical binning – even if followed by MP-EST or ASTRAL – is not statistically consistent under the first definition!

As shown in 10, MP-EST and ASTRAL used with weighted statistical binning were typically more accurate than when used without weighted statistical binning with respect to the estimation of species tree topologies and branch lengths, on datasets simulated under the multi-species coalescent model. However, there are some cases where using weighted statistical binning reduced accuracy, but these are limited to small numbers of taxa with very high levels of ILS.

Thus, statistical binning used with a coalescent-based summary method provides a blend of concatenation and coalescent-based methods: supergene trees are computed on concatenated alignments using fully partitioned maximum likelihood analyses, and then a coalescent-based summary method (such as MP-EST) is applied to these supergene trees to estimate a species tree. However, importantly, all maximum likelihood analyses used in statistical binning are fully partitioned, and this is critically important to the statistical properties that can be proven about these methods.

Liu et al. 17 argue against the use of statistical binning, saying “We show that approaches such as binning, designed to augment the signal in species tree analyses, can distort the distribution of gene trees and are inconsistent”. The claim that they show phylogenomic pipelines using naive binning 21 (not the same as statistical binning) are inconsistent is invalid, since they argue through the use of simulations (their own, and by reference to other simulations). However, statistical consistency or inconsistency is a property of behavior as the amount of data goes to infinity — and simulations are limited to finite data. Therefore, simulations by definition cannot prove statistical consistency or inconsistency. Hence, their paper does not establish inconsistency for naive binning, nor for any kind of binning approach.

Although their simulation study only examined naive binning 21, Liu et al. also expressed concerns about the potential for statistical binning to have a deleterious effect on species tree estimation. Therefore, their study does raise the following question: Are statistical binning pipelines statistically consistent under the multi-species coalescent model?

We can definitively answer this question with respect to the first meaning of statistical consistency: as shown in 10, phylogenomic pipelines that use weighted statistical binning followed by coalescent-based summary methods such as MP-EST will converge in probability to the true species tree as the number of loci and sites per locus both increase. We sketch the proof, to illustrate how the different aspects of the algorithm design are used to ensure statistical consistency. First, 10 showed that as the sequence length per locus increases, estimated gene trees converge to the true gene tree, and their bootstrap support values also converge to 100%; hence, the binning procedure produces bins that contain genes with the same gene tree topology with probability converging to 1. Recall that in a weighted statistical binning analysis, the loci that are placed in the same bin are then analyzed under a fully partitioned maximum likelihood analysis. Under the assumption that all the loci in the same bin have evolved down the same tree, a fully partitioned maximum likelihood analysis will converge to the tree associated to the bin as the length of the sequences for each locus increased. Then, in a weighted statistical binning analysis, this “supergene” tree would be replicated as many times as the number of loci in its bin. Finally, as the number of loci increases, the distribution of the supergene trees will converge to the true gene tree distribution defined by the species tree. Therefore, if the supergene trees are analyzed using a summary method (e.g., MP-EST or ASTRAL) that is statistically consistent given true gene trees, the pipeline described would be statistically consistent (in the first sense) under the multi-species coalescent model.

However, as shown in 10, pipelines using unweighted statistical binning (which does not replicate the supergene trees) are very different. In particular, as the number of genes increases, the estimated distribution of supergene trees will converge to the flat distribution in which every possible gene tree topology appears exactly once. Hence, no matter how these supergene trees are combined, the true species tree cannot be estimated with high probability. Therefore, pipelines that use unweighted statistical binning are not consistent even in the weak sense.

But what about the second interpretation of statistical consistency, where the number of sites per locus is fixed, but the number of loci increases? Are pipelines that use weighted statistical binning, followed by a summary method such as MP-EST or ASTRAL, statistically consistent under this second meaning (i.e., the strong sense of statistical consistency)? Does the Roch and Steel result help shed any light on this issue? Before we can answer this, we need to understand the difference between unpartitioned and fully partitioned maximum likelihood analyses.


Suppose we are given multiple sequence alignments for p different loci, and we wish to compute a maximum likelihood tree on the concatenated alignment under the Jukes-Cantor site evolution model, where a Jukes-Cantor model tree consists of a rooted binary tree T and numeric model parameters (the branch lengths of the tree). In an unpartitioned Jukes-Cantor maximum likelihood analysis we assume all the sites within the alignment evolve down a single Jukes-Cantor model tree, and we seek the tree and its numeric model parameters that is most likely to have generated the observed data.

In a fully partitioned Jukes-Cantor maximum likelihood analysis, we no longer assume that all the sites evolve down a single Jukes-Cantor model tree; instead, we assume that the different parts of the concatenated alignment each evolves down its own model tree, and the only constraint we make is that the model tree for the different parts share the same tree topology. Hence, we allow the numeric model parameters (i.e., branch lengths) to differ between the different loci. Therefore, the result of a fully partitioned Jukes-Cantor maximum likelihood analysis on a concatenated alignment from p loci is a set of p (potentially different) Jukes-Cantor model trees, each sharing the same tree topology. Equivalently, the result is a single tree topology T, but also p branch lengths for each branch in T.

We will show that fully partitioned maximum likelihood analyses and unpartitioned maximum likelihood analyses have very different theoretical properties. First, consider how maximum likelihood evaluates a sequence alignment where the sequences for the different species are all identical (i.e., an input of 100 sequences each identical to AACATAG). It is easy to see that all tree topologies are equally good under Jukes-Cantor maximum likelihood, and so cannot be distinguished under the maximum likelihood criterion. We will refer to alignments of this type as “invariant alignments” and loci that have invariant alignments as “invariant loci”. Now consider a concatenated alignment based on p loci, where all the loci but one are invariant. Under a fully partitioned Jukes-Cantor maximum likelihood analysis of such a dataset, because the p-1 invariant multiple sequence alignments fit every tree topology equally well, they do not impact the fully partitioned Jukes-Cantor maximum likelihood analysis of the concatenated alignment, and the result is identical to the Jukes-Cantor maximum likelihood analysis of the single variable multiple sequence alignment (see 10 for details) . However, as Roch and Steel showed (see below), unpartitioned maximum likelihood analyses converge (as the number of loci increases) to a maximum parsimony analysis of the single variable multiple sequence alignment. Thus, unpartitioned and fully partitioned maximum likelihood analyses behave very differently, and theoretical results (positive or negative) about one do not imply the same results for the other.

Roch and Steel’s proof that concatenation is statistically inconsistent under the multi-species model

Roch and Steel 1 prove that unpartitioned maximum likelihood is statistically inconsistent under the multi-species coalescent model, under the assumption that the gene sequence evolution model is the r-state symmetric model (i.e., models such as Jukes-Cantor, in which substitutions between every pair of distinct states are equiprobable). They establish this proof by showing that for some model species tree with very high levels of ILS and very low rates of evolution (in which all loci have the same sequence length), with high probability, nearly all loci will be invariant (i.e., their sequence alignments will have no changes on them). They then show that under the r-state symmetric model, unpartitioned maximum likelihood will converge to maximum parsimony as the number of loci increases (see Appendix). Roch and Steel then argue that under these conditions, as the number of loci increases, maximum parsimony will converge to a tree that is different from the species tree. In other words, unpartitioned maximum likelihood on the concatenated alignment will be positively misleading because maximum likelihood will be identical to maximum parsimony under some conditions (for large enough numbers of loci), and maximum parsimony will be positively misleading under the multi-species coalescent model.

Consider an input where the first p-1 loci have invariant multiple sequence alignments, and the pth locus has a variable sequence alignment. By the argument provided by Roch and Steel, as p→∞, an unpartitioned maximum likelihood analysis will converge to a maximum parsimony analysis of the concatenated alignment. Since p-1 of these loci are invariant, the unpartitioned Jukes-Cantor maximum likelihood analysis of the concatenated alignment converges to the maximum parsimony analysis of the pth locus.

For the same data, but under a fully partitioned maximum likelihood analysis, the first p-1 loci have no impact on the GTR maximum likelihood analysis, and so the result is identical to Jukes-Cantor maximum likelihood on the pth locus. Thus, for all values of p, the fully partitioned maximum likelihood analysis will always be the Jukes-Cantor maximum likelihood of the pth locus. Thus, as p→∞, the unpartitioned maximum likelihood analysis will converge to the maximum parsimony analysis of the pth locus but the fully partitioned maximum likelihood analysis will be the maximum likelihood analysis of the pth locus. Since Jukes-Cantor maximum likelihood and maximum parsimony are not identical methods (i.e, they can return different trees on the same dataset), this shows that fully partitioned and unpartitioned maximum likelihood analyses are different methods, and can return different trees.

To summarize, the theorem by Roch and Steel established that an unpartitioned maximum likelihood analysis can be inconsistent under the multi-species coalescent model, but their result does not apply to a fully partitioned maximum likelihood analysis. Most importantly, their proof uses the fact that the maximum likelihood analysis is unpartitioned to show that as the number of loci increases and nearly all the loci are invariant, the unpartitioned maximum likelihood analysis converges to maximum parsimony. This statement is explicitly not true for fully partitioned analyses, which are unaffected by invariant loci. Hence, their entire argument is restricted to unpartitioned maximum likelihood, and their proof does not apply to a fully partitioned maximum likelihood analysis.

Summary and Discussion

This article has focused on what has been established theoretically for some standard methods for estimating species trees (i.e., concatenation using maximum likelihood and summary methods such as MP-EST and ASTRAL), as well as for the newer approach of weighted statistical binning, followed by a summary method. Because the term “statistical consistency” has been used in two different ways in the literature, we have summarized what is known under each meaning: the weaker sense where both parameters (sequence length per locus and number of loci) increase, and the stronger sense where only the number of loci increases but the sequence length per locus is bounded, perhaps by a constant. We have also clarified that Roch and Steel’s theorem about concatenation using maximum likelihood being statistically inconsistent is restricted to unpartitioned maximum likelihood.

So, what do we know about statistical consistency (for either sense) of standard techniques for estimating species trees in the presence of ILS? Methods that have been proven to be statistically consistent in the first sense include the standard summary methods (e.g., MP-EST, ASTRAL, the population tree from BUCKy 22, etc.), methods that estimate trees directly from alignments, such as *BEAST and SVDquartets, and also weighted statistical binning, paired with standard summary methods. On the negative side, Roch and Steel’s theorem establishes that unpartitioned maximum likelihood can be statistically inconsistent. A separate argument shows that unweighted statistical binning is inconsistent in this first sense.

However, even for this weak sense of statistical consistency (where both the number of sites per locus and number of loci increase), the statistical consistency of many methods is still unknown. In particular, Roch and Steel’s theorem does not establish statistical inconsistency for fully partitioned maximum likelihood. Because the loci can have different tree topologies under the multi-species coalescent model, it seems likely that maximum likelihood analyses, even if fully partitioned, will be found to be inconsistent. However, the proof for an unpartitioned analysis being inconsistent provided by Roch and Steel in 1 depends on using an unpartitioned analysis, and so establishing the inconsistency (or consistency, as the case may be) of fully partitioned maximum likelihood analysis requires a different mathematical argument.

In terms of the second definition (and stronger sense) of statistical consistency, where only the number of loci increase but the number of sites per locus is bounded, very little has been established. In fact, the only established results are that unpartitioned maximum likelihood and unweighted statistical binning are both inconsistent. Finally, we do not know whether any of the standard summary methods, fully partitioned maximum likelihood analyses, or weighted statistical binning, are statistically consistent under this second definition.

In fact, to try to prove a method is statistically consistent or inconsistent under the second definition is very difficult; the only methods that have been proven to be statistically consistent under this definition (which constrains the number of sites per locus to be bounded) are explicitly designed to have this property, but all (to our knowledge) require some additional constraints (e.g., the same constant rate of evolution across all loci 19 , or a strict molecular clock 7). Proofs of inconsistency under the second definition have been established, but again only for unpartitioned maximum likelihood (proven by Roch and Steel) and unweighted statistical binning. Attempts to prove statistical consistency or inconsistency under this second definition have so far failed for any of the standard methods.

In other words, the major phylogenomic estimation methods in common use that are designed to address incomplete lineage sorting — summary methods such as MP-EST, ASTRAL, and the population tree from BUCKy, and co-estimation methods such as *BEAST — have been proven to be statistically consistent only in the first sense, where both the number of loci and number of sites increase. None of them have been proven to be consistent in the second sense, where the sequence lengths per locus are bounded. Similarly, the weighted statistical binning method is statistically consistent in this first sense, but it is not known if it is statistically consistent in the second sense.

It is also clear that Roch and Steel’s theorem (as well as the observations provided by their proof) is limited to unpartitioned maximum likelihood. Therefore, it cannot help us understand the theoretical properties of fully partitioned maximum likelihood, and hence also cannot help us understand the theoretical properties of methods that use fully partitioned maximum likelihood (such as weighted statistical binning). Therefore, any attempt to establish the statistical consistency or inconsistency of these methods will need to use other mathematical arguments.

In other words, the established theory regarding species tree estimation methods is limited. Yet, performance in practice (i.e., on simulated and biological data) suggests that many methods (including concatenation) provide good accuracy under some model conditions. Unfortunately, as discussed in 7,8,10,13, most simulation studies have explored performance only on very small datasets (e.g., with at most ten species) and under unrealistic conditions. For example, 14 examined performance on model species tree with only 5 species and very high levels of ILS, where sequence evolution was under a strict molecular clock, and used 1000 sites per locus (so that the sequence lengths per locus are too large to avoid recombination events. Instead, to understand the relative performance of coalescent-based summary methods and concatenation, more extensive analyses based on biologically realistic conditions (and hence based on short sequences, or modeling sequence evolution with recombination within loci) are needed.

Finally, despite the interest in coalescent-based species tree estimation, the current methods are still in their infancy, and new methods will need to be developed in order to obtain highly accurate species trees under realistic conditions. While the focus of this article has been on the performance of concatenation analyses, and the implications for coalescent-based summary methods that use weighted statistical binning, alternative approaches have been developed that do not require the estimation of gene trees or supergene trees (e.g., 10111516181920). Given the theoretical and empirical challenges in producing accurate gene trees, these approaches may provide the best accuracy for genome-scale phylogenomic analysis.

Statistical consistency (both senses) of standard species tree estimation techniques

Statistical consistency of some standard methods

We present the current status with respect to statistical consistency (of the first or second kind) of some standard phylogenomic estimation methods. The first column is for the first meaning of statistical consistency, which states that the species tree estimated by the method will converge to the true species tree as the number of loci and number of sites per locus both increase. The second column is for the second meaning, which states that the species tree estimated by the method will converge to the true species tree as the number of loci increases, even for bounded number of sites per locus. We also cite the paper in which the theoretical result is established.

Consistency – first kind Consistency – second kind
Unpartitioned concatenated maximum likelihood NO (1) NO (1)
Fully partitioned maximum likelihood UNKNOWN UNKNOWN
Unweighted statistical binning followed by consistent summary method (e.g., ASTRAL) NO (10) NO (10)
Weighted statistical binning followed by consistent summary method (e.g., ASTRAL) YES (10 ) UNKNOWN

Appendix 1: Quote from Roch and Steel’s Paper

An outline of the proof of the main theorem is as follows: We show that the expected proportion of sites that are constant can be made arbitrary large with low rates of evolution (the lower bounds are formalized in Claim 4) and that the empirical frequencies of site patterns is concentrated around the expected values (Claim 2). When there are a large enough number of invariable sites, it can be shown that likelihood scores and parsimony scores converge to the same answer (formalized in Claim 1). Thus trees that have better parsimony score have better likelihood under these scenarios. Therefore, it suffices to show that parsimony is not statistically consistent under arbitrary low rates of evolution.

Sebastien Roch and Mike Steel, “Likelihood-based tree reconstruction on a concatenation of sequence datasets can be statistically inconsistent”, Theoretical Population Biology 100 (2015): 56-62

Competing Interests

The authors have declared that no competing interests exist.

]]> 0
Concatenation and Species Tree Methods Exhibit Statistically Indistinguishable Accuracy under a Range of Simulated Conditions Mon, 09 Mar 2015 13:25:54 +0000


For decades, phylogeneticists have recognized the potential for discordance between a gene tree and its species tree 1,2. Maddison’s seminal insights on the topic 3 illustrate how discordance may arise from any of the following sources: gene duplication and extinction (“hidden paralogy”), lateral gene transfer, and incomplete lineage sorting (ILS). Given the computational and theoretical difficulties of adequately modeling gene descent, phylogeneticists have traditionally not attempted to reconstruct species phylogenies with models that explicitly account for sources of gene trees-species tree discord, and have instead opted to increase the total number of loci included in their concatenated data matrices under the expectation that, on average, the signal of multiple genes will recover the true species tree.

In recent years, however, molecular systematists have begun eschewing traditional supermatrix methods in favor of explicit modeling of the processes that may generate discordance. Although issues related to gene duplication and lateral gene transfer have received some attention 4,5, most efforts have concentrated on accounting for the effects of ILS by modeling the multispecies coalescent 6. The barrage of new methods now available to account for error due to ILS (Table 1) and the general transition away from supermatrix methods seem to follow the predicted emergence of a new paradigm in molecular systematics 7.

Table 1. Species tree methods accounting for incomplete lineage sorting in molecular systematics

Note that this list is not exhaustive.

Method Details
BATWING8 Uses Metropolis–Hastings algorithms to estimate population histories and mutation rates. The ancestral DNA data are inferred at each node of the genealogical tree. Different types of sequence data can be analyzed. The program estimates demographic parameters, but do not allow for the effects of either recombination or selection.
*BEAST6 Explicitly models the multispecies coalescent, coestimating multiple gene trees within a shared species tree and effective population sizes of extant and extinct populations in a Bayesian Markov chain Monte Carlo framework.
BEST9,10,11 Estimates species tree topology, divergence times, and population sizes from gene trees under a multispecies coalescent model in a Bayesian framework; unlike *BEAST, BEST estimates each gene tree individually, then uses importance sampling to infer the species tree.
BUCKy12 Uses non-parametric clustering of genes to estimate concordance factors based on proportions of genes supporting a clade; these factors provide the basis for inferring species-trees. Implemented in a Bayesian phylogenetic framework, and thus integrates over gene tree uncertainty.
GLASS13 An estimator of pairwise divergence times under the multispecies coalescent. Applies pairwise estimates by single-linkage clustering to generate species trees with branch lengths, using multiple unlinked loci and allow the use of several alleles from each population.
iGLASS14 Improves GLASS estimation of divergence times by deriving the expected waiting time until the first interspecific coalescence occurs among independent loci for a pair of taxa.
MCMCcoal15 A Bayes method for simultaneous estimation of the species divergence times and current and ancestral population sizes. Uses multiple loci, the topology of the species tree is assumed known, and a MCMC algorithm integrates over uncertain gene trees and branch lengths (or coalescence times) at each locus as well as species divergence times. It can handle any species tree and allows different numbers of sequences at different loci.
MDC16 Parsimony-based method that heuristically searches for the species tree that minimizes the number of deep coalescences implied by the gene trees, and thus does not model the coalescent process.
MP-EST17 Uses pseudo-likelihood function of the species tree to obtain maximum pseudo-likelihood estimates of species trees, with branch lengths in coalescent units. It assumes no gene flow or horizontal gene transfer, but method is robust to a small amount of the former.
NJst18 The distance between two species is defined as the average number of internodes between two species across gene trees, then the species tree is estimated by the neighbor joining tree built from the distance matrix.
SNAPP19 Uses a polynomial-time algorithm that computes the likelihood of a species tree directly from the sequences under a finite-sites model of mutation effectively integrating over all possible gene trees. The method applies to unlinked biallelic markers and it is implemented in a Markov chain Monte Carlo sampler for inferring species trees, divergence dates, and population sizes.
STAR20 Uses summary statistics of coalescence times by ordering the expected ranks of the coalescences among sequences, which is consistent with the ancestral order of populations in the species tree. It is resistant to variable substitution rates along the branches in gene trees.
STEAC20 Uses summary statistics of average coalescence times to estimate the species tree. The gene trees are generated from multilocus sequences using a consistent method for gene tree estimation (e.g. maximum likelihood) without the molecular clock assumption.
STEM21 Analytically derives the maximum likelihood species tree for a set of gene trees with branch-length information, explicitly modeling discord as a function of the coalescent process.
SVDquartets22 Infers relationships among quartets of taxa under the coalescent model using techniques from algebraic statistics; it accounts for mutational and coalescent variance. Uncertainty in the estimated relationships is quantified using the nonparametric bootstrap.
STELLS23 Uses coalescent-based maximum likelihood approach for inferring the species tree from a set of gene tree topologies (no branch lengths), which are used to estimated branch lengths in coalescent units.; the algorithm assumes free recombination between genes. It assumes the gene tree topologies are correctly inferred, therefore does not model uncertainty of gene tree estimation from sequences.

Inadequate model specification is known to introduce systematic bias into phylogenetic analyses, and thus the drive to account for potential error introduced by ILS has merit on theoretical grounds 16,24,25. However, the conditions under which traditional supermatrix approaches to molecular systematics necessarily result in less accurate topologies than methods accounting for population genetic processes remain unresolved, and recent empirical studies have reported high congruence between supermatrix and species tree methods 26,27,28.

Moreover, while theoretical studies have demonstrated the existence of species trees for which discordant gene trees are more likely than genealogies that agree with the species tree (the “anomalous gene trees” [AGTs] of Degnan and Rosenberg 24), such species histories are apparently unlikely to mislead most phylogenetic analyses. Huang and Knowles 29 have shown that, for species trees containing anomaly zones (i.e. short internal branches that result in AGTs), the probability of observing a polytomous gene tree is greater than the probability of recovering an anomalous gene tree. Thus, for species trees with short, deep branches, gene trees are more likely to be uninformative than to be outright misleading.

Given such empirical and theoretical findings, traditional supermatrix approaches may still be adequate for densely sampled data matrices and for clades evolving under non-extreme rates of change. However, simulation studies evaluating the accuracy of phylogenetic methods that explicitly treat error due to ILS do not always compare these methods to concatenation 30. Those simulation studies that do make comparisons to concatenation have not demonstrated uniformly superior accuracy of ILS-sensitive methods across all simulated conditions 31, and often use simple models of sequence evolution and/or few taxa 6,12,25, conditions that do not reflect typical empirical studies.

Furthermore, simulation studies concluding that concatenation is statistically inconsistent in the presence of anomaly zones and gene tree discordance often do not demonstrate superior accuracy by species tree methods under the same conditions 32, and are thus uninformative as to the practical merits of supermatrix and species-tree paradigms. Even when concatenation is compared to species-tree methods, it is not always specified whether apparent differences in performance are statistically significant 12; but see 33 and discussion below of other recent simulation studies that explicitly assess the accuracy of concatenation vis-à-vis species-tree methods.

Huang et al. 30 demonstrate how two sources of error – mutational variance and coalescent variance – impact the estimation of species trees using methods that account for ILS, but do not make direct comparisons to concatenation. Many empiricists still rely on concatenation for producing phylogenetic hypotheses for many taxa (e.g. ) and calibrating phylogenies with fossil information 34,35,36. Empiricists are likely to continue using concatenation for the foreseeable future, until methods that co-estimate genes trees and species trees *BEAST6 become more computationally feasible 37, and until truly synthetic and well-understood approaches to simultaneously account for the various sources of gene tree-species tree discord are developed. Thus, it is crucial that the relative accuracy of concatenation, vis-à-vis species-tree approaches, be explicitly delimited.

Herein we extend the experimental design of Huang et al. 30 to include traditional supermatrix methodologies, and explicitly assess the relative accuracies of species-tree methods and concatenation in estimating phylogeny. We find that the species-tree methods fail to outperform concatenation, and argue that, as a matter of best practices, studies involving methods sensitive to sources of gene tree-species tree discord should continue to compare results against those obtained by concatenation.


Huang et al. 30 performed simulations to quantify the effects of mutational and coalescent variance in species-tree estimation. Specifically, they compared the performance of two methods – STEM 21 and MDC 16 – in resolving the true species tree under a number of conditions, but did not present comparisons to the supermatrix approach. Our goal was to explicitly make these comparisons. The relevant steps of the simulations performed by Huang et al. 30 involved 1) generating a species tree under a uniform speciation model, 2) simulating coalescent gene trees for each species tree, 3) simulating DNA sequences under a specified model of nucleotide evolution along the branches of each gene trees, 4) estimating gene trees from the simulated DNA matrix, 5) estimating species trees from the estimated gene trees using MDC and STEM, and 6) calculating the discord between the true species tree and the estimated species tree. Our simulations followed the same trajectory, except that we do not estimate gene trees separately and instead estimate species trees directly from concatenated matrices. Rather than rerun the simulations of Huang et al. 30 in full, we obtained from the authors their calculated measures of discord between results obtained for each species tree method and the true tree under a variety of conditions. We then compared these discord measurements to values that we produced by identically parameterized simulations analyzed under concatenation.

We obtained 50 species trees of eight taxa used by Huang et al. 30, which those authors generated in the Mesquite software package 38 under a uniform speciation model. We simulated 540 coalescent gene trees for each of the 50 species tree under a neutral coalescent model, constant population size, and no migration in the program ms 39 using the script provided by Huang et al. 30. For each of the 50 species trees, the simulations produced gene trees at total tree depths of 1N and 10N generations; these values reflect a range of total lineage duration. We used only the subset of all simulated trees with a single individual per species, in accordance with the common practice for supermatrix methods. For each gene tree, sequences of 1000 base pairs were simulated using Seq-Gen 40 under an HKY model of nucleotide substitution 41, a transition-transversion rate ratio of 3.0, gamma-distributed rate heterogeneity shape parameter of 0.8, and nucleotide frequencies of A=0.3, C=0.2, T=0.3, and G=0.2 for the ancestral sequence, in accordance with the model of evolution used by Huang et al. 30. These simulations assume a single clock-like rate of sequence evolution across the tree.

Huang et al. 30 evaluate the effects of increasing the number of loci in the data matrix on the accuracy of MDC and STEM in obtaining the true tree, measuring the accuracy of trees produced with matrices of 3, 9, and 27 genes. To make our supermatrix analyses comparable to those of Huang et al. 30, the 540 genes of each of the 50 species trees were concatenated into matrices of 3, 9, and 27 genes; for each species tree at each tree depth, this procedure produced 180 matrices with 3 loci (1N/3 loci and 10N/3 loci), 60 matrices with 9 loci (1N/9 loci and 10N/9 loci), and 20 matrices with 27 loci (1N/27 loci and 10N/27 loci). These concatenated data sets were analyzed in MrBayes 3.2.2 42 under the same model of nucleotide substitution used to produce the sequences (i.e. HKY), a uniform molecular clock model, and Dirichlet distributed nucleotide frequencies. MrBayes was run for each tree estimation until the standard deviation of split frequencies fell below 0.01. Estimating species trees using the substitution model under which they were simulated ensured that model misspecification was not an additional source of error in comparing method accuracy. After discarding the first 25% of the posterior as burn-in, we sampled every 100 trees from the posterior distribution to produce a majority rule consensus tree with all compatible clades added (using the “allcompat” command in MrBayes). The use of a consensus tree for subsequent analyses of discord between methods poses a potential problem, as calculations of topological discord (described below) may be biased in trees with polytomous clades. This potentiality is not addressed by Huang et al. 30, and their measure of discord was also based on the 50% majority rule consensus with all compatible clades. However, adding all compatible clades to the majority rule consensus ensures that most trees are fully resolved, and manual inspection of a representative subset of our estimated species trees found no polytomous clades. Encouragingly, recent research indicates that the majority-rule consensus is expected to minimize the discord between the estimated and the true tree 43.

Topological discord between estimated species trees and their true species tree was measured using the Robinson-Foulds (RF) distance 44 in the R package phangorn 45. The RF statistic measures the distance between two unrooted trees as the summation of the number of bipartitions unique to the first tree and the number of bipartitions unique to the second tree. At each tree depth, mean discord was calculated by averaging RF distances for each of the 180, 60, and 20 concatenated matrices associated with each of the 50 species trees. We then took the average of the means for all 50 species trees for each matrix size, such that the mean discord for matrices consisting of 3 loci, 9 loci, and 27 loci was estimated from a total of 9,000, 3,000, and 1,000 phylogenetic analyses, respectively. We note that Huang et al. 30 ran 20 independent replicates for each species tree at each permutation of tree depth and sampling design; we did not adopt this procedure here, as a single iteration of the simulations produced sufficiently large sample sizes. ANOVA statistical tests were used to evaluate differences in mean RF distances within and across methods using the values obtained by Huang et al. 30. We did not investigate whether clade support, rather than topological disparity, differed significantly between concatenation and species tree methods, but this topic merits further research. The scripts used to perform the simulations are deposited at

We acknowledge that eight-taxon matrices are significantly smaller than the matrices of most molecular phylogenetic analyses – though still larger than those used in many simulation analyses – and also recognize that STEM and MDC are not the most sophisticated ILS-sensitive species tree methods available. Our goal, however, was to perform a general test of the importance of considering concatenation when choosing among phylogenetic estimation methods. For the sake of comparability to Huang et al. 30, who seek explicitly to make recommendations regarding method choice, we did not expand their basic approach.


Comparisons among methods of a given matrix size demonstrate that mean discord values are statistically indistinguishable, with two exceptions: 1N/3 loci and 1N/9 loci (Table 2, Fig. 1). For 1N/3 loci, MDC performed significantly worse than the other methods under consideration, while at 1N/9 loci, STEM exhibited statistically superior performance. Within methods, increasing the number of loci results in a significant increase in accuracy for both 1N and 10N simulations (Table 2, Fig. 1). Across all simulations, concatenation is as or more accurate than either STEM or MDC, except at 1N/9 loci, where STEM appears, on average, to produce slightly less discordant trees than concatenation.

Under the simulated conditions, increasing the number of loci sampled improved accuracy across all inference methods, as is expected for methods that exhibit statistical consistency. Likewise, all methods showed increased accuracy for analyses of taxa related by longer effective internal branches, i.e. branches characterized by smaller ancestral populations, less ancestral polymorphism, and/or increased generational turnover (Fig. 1). Still, the RF values were surprisingly high at 1N, reflecting the inherent difficulty in confidently resolving recent divergences.

Table 2. Results of ANOVA tests for differences in accuracy among methods (concatenation, STEM, MDC) with different numbers of loci and between methods with the same number of loci

Accuracy measured by the topological distance (Robinson-Foulds distance) between the true and the obtained tree. Significant differences denoted with asterisks. See also Figure 1.

3 loci 9 loci 27 loci Concatenation MDC STEM
1N <0.001* 0.0255* 0.6384 <0.0001* <0.0001* <0.0001*
10N 0.3601 0.8233 0.8505 <0.0002 <0.0001* <0.0011*

Fig. 1: Mean and standard deviation estimated from a total of 9,000, 3,000, and 1,000 phylogenetic analyses for the 3 loci, 9 loci, and 27 loci analyses, respectively.

The Phylogeneticist’s Toolbox

Our simulations demonstrate that concatenation can perform as well or better than methods that attempt to account for sources of error introduced by ILS (Table 2, Fig. 1). This result corroborates some recent empirical findings that concatenation and species-tree methods can produce near-identical results 27,28,33,46, and is concordant with simulation studies indicating that ILS-sensitive methods do not consistently outperform concatenation under all conditions 31,33,46,47,48,49,50. Indeed, concatenation can outperform species-tree methods when levels of ILS are low, few loci are used, or gene trees have low phylogenetic signal 33,47, conditions that likely reflect the reality of many phylogeneticists working with mitochondrial and nuclear markers generated by traditional methods of sequencing to estimate large scale phylogenetic relationships 34,35,36,51,52,53. In addition, simulations by Patel et al. 49 and analysis of mammalian empirical data 50 indicate that concatenated analyses benefit from increased power vis-à-vis “shortcut” coalescence methods (i.e. those that estimate gene trees and species trees separately) when trying to resolve deeper divergences, especially when more loci are used. In these cases, high mutational error and gene tree heterogeneity are significantly more detrimental to shortcut coalescent methods than they are to concatenation, resulting in more accurate concatenation-based trees despite the incorrect assumption that all loci share a common topology. Binning techniques that use concatenation to produce “supergene trees” as input for coalescence-based analyses may improve the accuracy of species tree estimation by alleviating the bias introduced by using numerous small loci 54.

Gene tree-species tree discordance is a well-documented phenomenon, and we do not question that it poses a significant challenge to accurate inference of phylogeny; however, the growing body of research indicating that concatenation is not uniformly worse than discordance-sensitive methods should give phylogeneticists pause. That concatenation remains reasonably accurate under a range of conditions suggests that it retains practical utility. Where results from concatenation are concordant with those obtained from discordance-sensitive methods, the phylogeneticist can be reasonably confident that the resultant phylogenetic hypothesis is robust. Where concatenation and discordance-sensitive methods disagree significantly, the disparity may indicate regions of the tree in need of focused attention and might suggest possible biological factors confounding phylogeny estimation 55, or may indicate critically violated assumptions in one or both approaches 50. In these cases, morphological, physiological or ecological evidence may provide support for one hypothesis over the other.

Though many methods have been developed to account for error due to ILS (Table 1), other biological processes (hidden paralogy and horizontal gene transfer) can confound accurate phylogenetic inference; however, empirical researchers cannot know a priori the relative contributions of these processes to the genealogical history of a lineage, or to what extent discordance among gene trees may be due to methodological or sampling error (e.g. poor taxonomic sampling; model misspecification, species misidentification, long branch misplacement, arbitrary resolution of gene trees) 50,56. While concatenation has been rightly criticized for not accounting for known sources of phylogenetic error, shortcut coalescence methods assume that all error is due to incomplete lineage sorting – that is, that gene trees are known with certainty 50. To date, no species tree method simultaneously accounts for all sources of gene discordance 57; until well-understood, synthetic approaches to co-estimating gene trees and species tree and accounting for all sources of gene tree-species tree discord are developed, the continued use of concatenation has practical merit insofar as it avoids making a priori assumptions about the extent to which various sources of uncertainty have influenced the evolutionary history of the taxa under consideration. Thus, we argue that concatenation remains a useful component of the phylogeneticist’s toolbox as a “null” methodology against which results from alternative methods of inference can be compared and interpreted.

That concatenation remains an important part of the phylogeneticist’s toolbox is borne out by recent analyses indicating that the method exhibits greater power to overcome sampling error and discrepant patterns of homoplasy 50, and that “concatalescence” methods can make explicit use of concatenation in improving the quality of input trees prior to coalescence-based estimation of the species tree 54. Moreover, estimation of time-calibrated phylogenies combining molecular and morphological data (extant and fossil) is only possible in concatenation analyses, and concatenation has been shown to be more resilient to the effects of missing data than species tree methods 58,59.

We did not test improvement in accuracy of species tree estimation by using multiple individuals of the same species. New species-tree methods that have been developed to handle phylogenomic data of hundreds of loci have also used one individual per species when testing for accuracy and performance 46, a common practice in studies aiming to infer large-scale interspecific relationships 34,36,51,52,53. Although our results and the findings of phylogenomic studies indicate concatenation exhibits reasonable accuracy under certain conditions, we recognize that phylogenetic studies using either supermatrix or species-tree methods may benefit from including multiple individuals per species, especially for recent divergences 30,49, as the additional data help account for intraspecific allelic diversity when estimating interspecific relationships.

Research is only beginning to shed light on the relative accuracies of species-tree methods and concatenation for genomic-scale data matrices with thousands of taxa. Simulations indicate that coalescent-based methods improve accuracy for genome-scale matrices 33, though adding loci for which gene tree estimates are poor may be more detrimental for ILS-sensitive methods than for concatenation 49,50. Recent empirical research suggests that rapid radiations may make some internal nodes permanently intractable, regardless of methodology or sampling (e.g. few taxa and genes 60, many taxa and few genes 61, or few taxa and many genes 27).

The phylogeneticist’s toolbox continues to grow as our understanding of the dynamics of evolving populations improves, and the last year has witnessed the development of numerous methodological tools and considerations for practicing systematists. For example, recent research has produced posterior predictive model checking for coalescent models 62 and new methods of data exploration to reveal the causes and consequences of patterns of gene tree discordance 63,64. We do not dispute that the proliferation and refinement of new methodologies is crucial to the continued success of phylogenetic research. Our purpose here is to demonstrate the importance of explicit comparison of contemporaneous and classical methodologies, and to highlight that the addition of theoretically important parameters may not always produce practical gains. We conclude that concatenation exhibits statistically comparable accuracy under a range of sampling and tree depth conditions vis-à-vis some existing species tree methods, and urge molecular phylogeneticists to thoroughly evaluate the performance of methods that model gene tree-species tree discord against concatenation.

Competing Interests

The authors have declared that no competing interests exist.

Correction: Building a Phylogenomic Pipeline for the Eukaryotic Tree of Life – Addressing Deep Phylogenies with Genome-Scale Data Fri, 06 Mar 2015 22:15:10 +0000


In table 1 an incorrect value was provided for the “Contigs” value for “Filamoeba nolandi” (Column 3, row 2). The corrected table is provided below:

Table 1: Transcriptome statistics for Amoebozoa used in case study.
Taxon ATCC Contigs SSU Bacterial Eukaryotic Unknown Genes
Filamoeba nolandi 50430 21671 17 2409 12205 7057 171
Pessonella sp. PRA29 17104 31 1783 8695 6626 199
Trichosphaerium sp. 40318 22137 15 2252 9541 10344 35
Eukaryota sp. JRG-2011 50979 14754 33 1940 7303 5478 178
Stereomyxa ramosa 50982 21354 38 2141 9854 9321 210


Grant JR, Katz LA. Building a Phylogenomic Pipeline for the Eukaryotic Tree of Life – Addressing Deep Phylogenies with Genome-Scale Data. PLOS Currents Tree of Life. 2014 Apr 2. Edition 1. doi: 10.1371/currents.tol.c24b6054aebf3602748ac042ccc8f2e9.

View Article

]]> 0
One Tree to Link Them All: A Phylogenetic Dataset for the European Tetrapoda Fri, 08 Aug 2014 10:49:46 +0000


The use of phylogenetic data into ecological analyses has grown rapidly in the last decades, giving rise to new disciplines such as community phylogenetics which incorporate information on species relatedness into the study of community structure1,2, as well as to studies of large-scale distribution of species and their phylogenetic diversity3,4. Additionally, the integration of ecological and evolutionary information holds promise to improve ecological forecasting in the current context of climate and land change and biodiversity loss5,6. Since the pioneering work of Dan Faith7, conservation biology has long recognized the importance of considering phylogenetic diversity as a relevant feature for conservation8,9. The EDGE framework is, in this regard, an important initiative that combines the evolutionary distinctiveness of species (i.e. the evolutionary contribution of a species to the tree of life) with globally endangered risk assessment to derive conservation priorities10. Recent works have also focused on how future climate and land use change could further jeopardize the tree of life in certain parts of the world11,12.

To foster the developments of these emergent fields and timely questions, detailed and broadly sampled phylogenetic hypotheses are needed to appropriately integrate evolutionary information into ecological and conservation studies. Recent phylogenomic studies have improved our understanding of the evolutionary relationships within the main Tetrapoda groups, especially at high levels such as families and orders. For instance, Roelants and colleagues13 clarified the relationships between global amphibians at the family level, while Pyron et al.14 performed a similar achievement on Squamata, sampling all families and sub-families. Concerning birds, Hackett et al.15 elucidated the inter-ordinal relationships of extant birds, and a later study16 confirmed the partly controversial results found by Hackett and colleagues. Despite these achievements, we still lack detailed species-level phylogenies for such groups. Moreover, there is a lack of phylogenies for particular regions (but see 17), as systematists mainly focus on building species-level phylogenies for entire clades. Although it is of obvious interest, research areas such as community phylogenetics and conservation planning do not specifically require complete taxonomic sampling, but rather complete spatial, or biogeographic, sampling. In order words, ecological studies that wish to integrate evolutionary data usually require a phylogenetic hypothesis for the entire species pool under study, which might be along a specific gradient18 or a continental scale assessment11,19 . For instance, incorporating phylogenetic diversity in reserve design or gap analysis only require a complete phylogenetic tree for the entire group with the region of interest (see for example 19, 85). It should however be noted that since the complete coverage only concerns Europe, estimates of phylogenetic uniqueness are therefore biased and should be accounted for in the analysis of the data (e.g. 86).

For that purpose, we here construct and provide a phylogenetic dataset for all Tetrapoda species that occur in the entire European sub-continent (including Turkey) built on relevant phylogenetic data in Genbank and consensus tree knowledge, based on a supermatrix-supertree mixed approach20 . We also check the congruence of the phylogenies obtained with previous evolutionary studies.


Squamata and Testudinae

The list of European Squamata species was extracted from Maiorano et al21. DNA sequences of 7 nuclear (BDNF, c-mos, NT3, PDC, R35, RAG-1, RAG-2) and 6 mitochondrial loci (12S, 16S, COI, cytB, ND2, ND4) were downloaded from Genbank with PHLAWD22. These regions have been shown to be useful for phylogenetic inference in previous studies of squamates according to Pyron et al14. Only 16 species of a total of 239 had no molecular data available in Genbank. In addition to Squamata species, we included 3 levels of outgroup taxa: Sphenodon punctata (closest living relative to Squamata); all 10 species of European turtles, two crocodilians (Alligator and Crocodylus) and two birds (Dromaius and Gallus); and finally two mammals (Mus and Pan). Genbank accession numbers are detailed in Table S1 (Appendix 1).

For each region, DNA sequences were aligned with MAFFT23 and checked by eye with Seaview24. Ambiguous alignment positions were trimmed with trimAl25. All the regions were concatenated in a supermatrix with FASConCAT26. The phylogenetic inference analysis was conducted with RaxML v. 7.8.127 using the GTRGAMMA model and employing the rapid hill-climbing algorithm28; we searched for 100 Maximum Likelihood trees applying a family tree constraint for squamates based on Pyron et al14. Bootstrapping was conducted with 1000 replicates to assess clade support.

The 100 ML trees were dated with penalized-likelihood as implemented in r8s29; we constrained 5 nodes based on fossil information extracted from Mulcahy et al.30: we set a minimum and a maximum age of 256 and 300 mya respectively for the most recent common ancestor (mrca) of all Reptilia31,32, a minimum and a maximum age of 239 and 250 mya respectively for the mrca of Birds and crocodilians32 , a minimum age for the mrca of Lepidosauria of 223 mya33,34, a minimum age of 111 mya for the stem branch of Amphisbaenidae34,35, and a minimum age of 93 mya for the stem branch of Alethinophidia36. The best smoothing value was determined by a cross-validation procedure, following 29 .

The data matrix and the phylogenetic tree with the highest likelihood are available in Treebase (accession number: S15708).


For Amphibians, we include here the phylogenetic tree constructed for a previous study19. The list of European Amphibian species was extracted from Maiorano et al21. We retrieved from GenBank sequences of phylogenetic informative regions that were available for at least 30% of the species: 9 mitochondrial (12S, 16S, COI, cytb, ND1, ND2, ND4, tRNA-Leu, tRNA-Val) and 2 nuclear (RAG-1, rho) regions. We found relevant molecular data for all species, but we excluded the two hybrid species Pelophylax grafi and Pelophylax hispanicus. We included Xenopeltis unicolor, Gallus gallus and Mus musculus as outgroups to root the tree. For each region, alignment was conducted with four programs (Clustal37, Kalign38, MAFFT23 , MUSCLE39). The best resulting alignment was selected based on Mumsa38, and checked visually. Ambiguous regions of each alignment were removed with trimAl25. All regions were concatenated in a supermatrix with FASconCAT. As with Squamata, we obtained 100 ML phylogenetic trees by conducting a phylogenetic inference analysis with RaxML, this time applying a family-level tree constraint based on Roelants et al13. A bootstrap analysis was conducted with 1000 replicates to assess clade support.

We dated the 100 ML trees with penalized-likelihood (r8s) using the following fossil data to constrain minimum ages for selected nodes: 155 mya for the crown-origin of salamanders40, 170 mya for Bombianura41, 250 mya for Batrachia42, 110 mya for the split of Pelobatidae and Pelodytidae families43, 145 mya for the split of Pelobatidae and Neobatrachia43, and 61 mya for the split of Plethodidae and Proteidae44. Additionally, we set a minimum and maximum age (312-330 mya) for the split between diaspid (Gallus gallus, Xenopeltis unicolor) and synapsid amniotes (Mus musculus), based on Benton and Donoghue45.

The data matrix and the phylogenetic tree with the highest likelihood are available in Treebase (accession number: S13561).


We include here 100 dated phylogenetic trees for 430 species of European breeding birds from Roquet et al20. This phylogenetic dataset was built upon sequences retrieved from GenBank for 10 mitochondrial gene regions (12S , ATP6 , ATP8 , COII , COIII , ND1 , ND3 , ND4 , ND5 , ND6) and six nuclear ones (28S , c-mos, c-myc , RAG-1 , RAG-2 , ZENK). The alignment procedure was the same as for Amphibians. We also performed 100 ML phylogenetic inference searches and standard bootstrapping (1000 replicates) with RaxML, applying a tree constraint at the ordinal level based on Hackett et al15. The 100 trees were dated with penalized likelihood (r8s) applying fossil calibrations for 14 clades (Table S2, Appendix 2). The best ML tree can be found in Treebase (study number 10770).


The phylogenetic data here included for mammals is based on the super-tree of Fritz and colleagues46; concretely, we extracted the resampled dataset of 100 fully resolved phylogenetic trees from Kuhn et al.47, where polytomies of the super tree from Fritz et al.46 were randomly resolved applying a birth-death model to simulate branch-lengths. Then, for each tree, we replaced the Carnivora clade with the update performed on a recent study48, which provides a better resolution and increases the sampling from 252 sp to all Carnivora species (286 sp). Later, we removed all non-European species. These modifications of the phylogenetic trees were done with the R package ape.

Phylogenetic inference

As stated before, for each taxon group except mammals we have conducted 100 ML inferences with RAxML. Every inference begins with a different starting tree, which is built by adding sequences one by one in random order, identifying their optimal location on the tree under the parsimony optimality criterion. Since sequences are added in random order, it is very likely that a different starting tree is generated at every search8788. RAxML searches were then performed with the method “lazy subtree rearrangement” (a variant of subtree prunning and regrafting method) under a ML framework. Like all heuristic search strategies, the Maximum Likelihood search strategy employed by RAxML is not guaranteed to find the most probable tree of the tree-space, and because of that, it is important to conduct multiple searches from different starting trees. To check if all the searches converged on trees with similar likelihoods, we performed the Shimodaira-Hasegawa test89 (SH) implemented in RAxML. In all cases, the likelihoods of the trees of a same group of taxa were not significantly different (p < 0.01). This increases our confidence on the trees found being close/similar to the most likely tree, and that the trees obtained do not result from the algorithm getting stuck in a local optima.

Supertree construction

The trees cited above were combined, after pruning the outgroups, by joining them with the R package ape; to do so we set divergence ages between these main groups based on the information retrieved in the webpage Timetree49: the divergence age between mammals and sauropsids (i.e. birds, turtles and squamates) was set to 324 mya, and the divergence age between amphibians and the rest of the groups was set to 361 mya. To build the final tetrapod tree, we randomly selected one tree from each of 100 trees available for each group. We repeated this approach 100 times to get 100 realisations of the tetrapod tree. These combinations were done randomly since the likelihoods of the trees of each group were not significantly different according to the SH test. The 100 dated trees of each group and the 100 dated supertrees for all European Tetrapoda are available from the Dryad digital repository (DOI: X).

Results and Discussion


The study of Pyron et al.14 constituted a major advance in our understanding of the phylogenetic relationships between the main lineages of Squamata. Their study had a broad taxonomic and molecular sampling: they included members of all currently recognized families and subfamilies, for which 7 nuclear and 5 mitochondrial loci were analysed. Here, we took profit of the knowledge derived from that study by incorporating a tree constraint to the family level based on their results. We also performed the analysis without the tree constraint (results not shown); the results were congruent with the first analysis, but the lack of a family-tree constraint yield low bootstrap (BS) support for the deepest nodes.

Our phylogenetic results are largely congruent with those of Pyron and colleagues14. We have similar levels of strong nodal supports except for the relationships between genera of Lacertidae; 67.8% of the nodes had a strong support (BS > 70%, Fig. 1, Appendix 3) and 13.1% of the nodes had a moderate support (BS 50-70%). In accordance with their study, we detected that some genera are not monophyletic: Ablepharus (Scincidae), Cyrtopodion (Gekkonidae), Zamenis (Colubridae). We also found strong evidence that Hierophis and Dolichophis (Colubridae) are not monophyletic genera, as D. cypriensis (which was not included in 14) is nested within Hierophis with a 100% BS.

Available dating studies on Squamata differ considerably on age estimates. For instance, a recent study30 estimated the squamate crown group to be c. 180 mya, while two other studies estimated the same group to be c. 240 mya50,51. Our estimates of divergence times are in general roughly similar to those of Kumazawa’s study50. It has been suggested30 that the use of only mitochondrial regions (which is not the case here) may bias the results towards older ages, but anyway differences in methodology and in taxon and molecular sampling make difficult to identify all the causes of those discrepancies.


The phylogenetic inference analysis for the amphibians yielded a particularly robust topology: 83.5 % of the nodes showed a strong support (BS>70%, Fig. 2, Appendix 3). Supported nodes of our ML trees were congruent with previous phylogenetic studies52,53,54,55,56,57. Concerning the divergence age estimates, we obtained younger ages for the deepest nodes compared to the work of Roelants and colleagues13, for instance, Batrachia was estimated with r8s to be c. 330 mya in that study, while we estimated it to be c. 300 mya; in contrast, we retrieved older ages for the shallowest nodes (e.g. we estimated that the divergence between Salamandra and Pleurodeles occurred 100 mya, Roelants and colleagues estimated it to have occurred c. 75 mya). These differences might be linked to the difference in molecular and taxon sampling: Roelants et al. sampled only one species per genera; and several families that were included in their work are not present in Europe and thus were not included in our supermatrix.


Supported nodes of our ML trees are congruent with previous phylogenetic studies (Anseriformes: Donne-Goussé et al. 58, Eo et al.59, Gonzalez et al.60; Galliformes: Gutierrez et al.61 , Dimcheff et al.62, Crowe et al.63, Kimball et al.64, Kriegs et al.65, Lislevand et al.66; Gruiformes: Fain et al.67; Procellariiformes: Penhallurick and Wink68; Ardeidae: Sheldon et al.69 ; Accipitridae: Lerner and Mindell70, Griffiths et al.71 ; Charadriiformes: Paton et al.72 , Thomas et al.73, Pons et al.74 , Bridge et al.75 , Paton and Baker76 , Fain and Houde77; Passeriformes: Alström et al.78, Nguembock et al.79, Treplin et al.80; Piciformes-Coraciiformes: Johansson et al.81, Benz et al.82; Strigiformes: Wink et al.83); 68.7% of the nodes had a strong BS support (BS>70%, Fig. 3, Appendix 3), and an additional 12.4% had a moderate support (BS=50-70%). Divergence age estimates were, in general, congruent with those obtained by Brown et al84.


The modification of the most recent mammals supertrees available on the literature47 with the update of Carnivora clade48 allowed to increase the phylogenetic resolution (only nine polytomies remain in the updated Carnivora clade) and to have a higher species sampling.

The importance of accounting for phylogenetic uncertainty

Phylogenetic information is sometimes incorporated in ecological analyses based on a single phylogenetic tree, assuming the tree is known without error. Any phylogenetic tree estimate will probably not be an exact representation of the true phylogeny due to possible bias or uncertainties such as molecular and taxon sampling, sequence alignment, homoplasy, or long-branch attraction9091. For all these reasons, it is important to include phylogenetic uncertainty in order to avoid overestimating our confidence in subsequent analyses (i.e. obtaining too narrow confidence intervals). This type of uncertainty can be accounted for in two ways: with a single consensus tree (in which unsupported nodes are collapsed into polytomies), or running the analyses with a range of trees and later summarising the results11, 93. The first approach (i.e. consensus tree) may not be preferable, as polytomies can influence the results of tree-based statistical analyses (e.g. see 92 for the influence of phylogenetic resolution on several community phylogenetics indices), and do not allow to fully explore the variation in ecological patterns resulting from phylogenetic uncertainty. Moreover, not all methods have been adapted to allow for polytomies, some of them require completely bifurcating trees (e.g. the EDGE index10). For all these reasons, we highly recommend to account for phylogenetic uncertainty by including a set of high-probability trees.

Data availability statement

The 100 dated supertrees for all European Tetrapoda and the 100 dated trees of each taxon group (amphibians, birds, mammals, squamates and turtles) are available from the Dryad digital repository (DOI: X).


We provide here a phylogenetic dataset constituted of 100 chronograms of European Tetrapoda species as a tool for ecological studies that aim to incorporate an evolutionary perspective, and for phylogenetic conservation assessment. This phylogenetic dataset is in general agreement with previous studies, and we expect it to be coarsely approximate with the “true” Tetrapoda evolutionary tree. Instead of providing the best ML tree for every group, we provide 100 trees (available on Dryad repository), as computing analyses with several trees allows taking in account phylogenetic uncertainty. Regarding the taxonomic sampling, the big majority of species are included (91%). On the other side, some molecular regions have low sampling, thus, this dataset will be useful until substantial amount of molecular data becomes available for a considerable number of species.

Competing Interests

The authors have declared that no competing interests exist.

]]> 1