New Phylogenomic Analysis of the Enigmatic Phylum Telonemia Further Resolves the Eukaryote Tree of Life

Abstract

Introduction

The eukaryote tree of life is one of the oldest working hypotheses in biology. Over the past 15 years, this tree has been extensively reshaped based on phylogenomic studies of an ever-increasing diversity of protists (i.e., microbial eukaryotes) (Burki et al. 2007; Rodriguez-Ezpeleta, Brinkmann, Burger, et al. 2007; Brown et al. 2013; Yabuki et al. 2014; Katz and Grant 2015; Burki et al. 2016; Brown et al. 2018). The current tree places most lineages into a few redefined “supergroups”—a taxonomic rank-less concept that was devised as inclusive clades reasonably supported by molecular and/or morphological evidence (Simpson and Roger 2002,, 2004; Keeling et al. 2005). As few as eight supergroups are recognized, sometimes arranged in even higher-order divisions assuming the root of the tree falls outside of these mega-clades: Diaphoretickes, including Sar, Haptista, Cryptista, Archaeplastida; Amorphea, including Amoebozoa and Obazoa; Excavata, including Discoba and Metamonada (Burki 2014; Simpson et al. 2017; Adl et al. 2019). Importantly, this view of the tree remains a working hypothesis and whilst most groups have been robustly and consistently supported, others have not and still await confirmation (this is mainly the case for Archaeplastida and Excavata) (Burki et al. 2016; Heiss et al. 2018).

The backbone of the tree also contains a series of long-known and newly discovered enigmatic protist taxa that do not fit within any established supergroup based on ultrastructure, ribosomal DNA phylogenies, or even protein-coding genes (Brugerolle et al. 2002; Walker et al. 2006; Heiss et al. 2011; Yabuki et al. 2012; see Adl et al. 2018 for a more exhaustive list of incertae sedis taxa within eukaryotes). Interestingly, phylogenomics suggested that many of these “orphan” taxa are themselves supergroup-level lineages, likely representing early diverging evolutionary intermediates between the established supergroups (Zhao et al. 2012; Brown et al. 2013; Cavalier-Smith et al. 2015; Janouškovec et al. 2017; Heiss et al. 2018; Lax et al. 2018). Individual orphan lineages can also turn out to be evolutionarily related, forming new larger deeply branching eukaryotic clade as in the case of the recently proposed “CRuMs” supergroup (Brown et al. 2018). Taken together, orphans can profoundly impact our understanding of ancient eukaryote evolution by defining new major branches in the tree and help establishing the branching pattern among all main eukaryotic groups.

The genus Telonema is one such orphan taxon that has proven difficult to place in eukaryote-wide phylogenies. Known for over a century, this genus contains only two named species, Telonema subtile (syn. Telonema subtilis) and Telonema antarcticum (Griessmann 1913; Klaveness et al. 2005). However, environmental sequencing of the 18S rDNA has shown that the genetic diversity of Telonema is relatively high in marine and freshwater environments and in fact likely represent several genera (Shalchian-Tabrizi et al. 2007; Bråte et al. 2010; Adl et al. 2018). Telonema and related sequences are frequently detected in these aquatic environments, and although they are usually not very abundant, they are thought to play significant ecological roles as biflagellated heterotrophic predators of bacteria and small-sized eukaryotes (Klaveness et al. 2005). The ultrastructure of T. subtile and T. antarcticum revealed unique morphological features, such as a multilayered subcortical lamina and complex cytoskeleton, that distinguish Telonema from all other eukaryotes (Klaveness et al. 2005; Shalchian-Tabrizi et al. 2006; Yabuki et al. 2013). Yet, Telonema also possesses structures that are characteristics of other groups: for example, the peripheral vacuoles that are similar to the cortical alveoli in Alveolata, tripartite tubular hairs on the long flagellum that are similar to the mastigonemes in Stramenopila, or wheel-like structures resembling extrusomes in some Rhizaria (Klaveness et al. 2005; Shalchian-Tabrizi et al. 2006; Yabuki et al. 2013). This combination of unique and shared ultrastructural features makes Telonema highly distinctive from other eukaryotic groups and has supported the notion that it forms a phylum-level lineage—Telonemia—without close relatives (Shalchian-Tabrizi et al. 2006). Accordingly, the first molecular phylogenies based on up to four genes provided only few clues as to the position of the group, other than it was not closely related to any known groups (Shalchian-Tabrizi et al. 2006). Later, more substantial analyses based on larger data sets containing 127–351 genes also failed to pinpoint with statistical significance the origin of Telonemia, but this was likely due at least in part to the small number of genes sampled from only one species (T. subtile) (Burki et al. 2009, 2012,, 2016; Lax et al. 2018).

In this study, we tested whether better gene coverage and taxon sampling for Telonemia using a phylogenomic approach can help to infer with confidence the phylogenetic position of this important group of eukaryotes. To this effect, we established cultures for three telonemid species, one of which corresponded to the described T. subtile, and generated high-quality transcriptomes for all. The availability of different data sets from closely related species allowed us to robustly identify the correct copies from contaminant sequences, a task that is otherwise not easy with single data sets for orphan lineages. Our analyses unambiguously place Telonemia in a sister position to the supergroup Sar, forming a new mega-assemblage that further resolves the tree of eukaryotes.

Results

Improved Data Set and Taxon Selection

To place Telonemia in the tree of eukaryotes, we assembled a 248-genes data set with a broad taxon sampling. All genes were carefully inspected by reconstructing gene phylogenies to evaluate orthology and detect contamination (see Materials and Methods). Given the taxonomic distribution of available genomic data sets and the tentative placement of Telonemia in previous phylogenomic works (Burki et al. 2009,, 2012,, 2016), special attention was made to include representatives of Diaphoretickes lineages with sufficient data. For the telonemids, at least two sequences from different data sets were required to group together in gene trees in order to be retained. This requirement became possible by the addition of the new telonemid data sets and greatly lowered the risk that unidentified contaminants were wrongly assigned to Telonemia. Overall, the sequence coverage for the three telonemid operational taxonomic units (OTUs) ranged between 58% and 99% (supplementary table 1, Supplementary Material online), which constitutes a large improvement over previous studies (e.g., Burki et al. 2016). Sequences corresponding to the bodonid prey Procryptobia sorokini, easily identifiable in gene trees from the presence of other related taxa (e.g., Neobodo designis, Bodo saltans), were reassigned accordingly. In total, 293 taxa were initially analyzed (supplementary fig. S1, Supplementary Material online), but this number was reduced to 109 OTUs to allow computational-heavy analyses. This final taxon sampling was selected based on the taxonomic affiliation and data coverage of the OTUs, excluding faster-evolving lineages when the selection of shorter branches was possible (fig. 1, Supplementary fig. S1, Supplementary Material online, and supplementary table 1, Supplementary Material online).

Fig. 1.

Open in new tabDownload slide

ML tree showing the phylogenetic position of Telonemia within eukaryotes. The tree corresponds to the best ML tree inferred from a concatenated alignment of 248 protein-coding genes using the LG+C60+G+F-PMSF model. Support at nodes is derived from 100 BP replicates and Bayesian PPs. Black circles represent maximum support in both analyses. A number appended to a black circle indicates full support in only one of the two analyses (i.e., BP/PP). Nodes with dashes (–) indicate an alternative branching in the Phylobayes tree (see supplementary figs. S2 and S3, Supplementary Material online). The ambiguous position of Haptista with this data set is highlighted with an exclamation mark (see text for more information).

Tree Inference and Testing

The 248 genes/109 OTUs supermatrix was analyzed by maximum likelihood (ML) and Bayesian tree reconstructions using site-heterogeneous mixture models of evolution, specifically LG+C60+F+Γ4 with posterior mean site frequency (PMSF) profiles and CAT+GTR+Γ4, respectively (see Materials and Methods). These models have been shown to be more robust against heterogeneity in the amino acid substitution process across sites, which are important to take into consideration especially along the most ancient branches of the eukaryotic tree (Lartillot and Philippe 2008; Quang et al. 2008; Wang et al. 2018). As shown in figure 1, both models recovered a similar overall topology, with most currently accepted supergroups receiving maximal support (e.g., Sar, Cryptista, Haptista, Discoba, Obazoa, and Amoebozoa). Most interestingly, in both analyses the position of Telonemia was identical and unequivocal, placed within Diaphoretickes as sister to the Sar group (bootstrap [BP] = 100%, posterior probability [PP] = 1.0). Cryptista was robustly inferred in a clade including the members of Archaeplastida, but the branching pattern in this clade was unclear: whilst the LG+C60+F+Γ4+PMSF tree remained unresolved (BP = 56% for the grouping of Cryptista with Rhodophyta), the CAT+GTR+Γ4 model recovered a monophyletic Archaeplastida with 1.0 PP (fig. 1). The position of Haptista (and its possible sister taxon, Ancoracysta twista; see Janouškovec et al. 2017) also remained ambiguous showing a sister relationship to the Archaeplastida + Cryptista clade with the LG+C60+F+Γ4+PMSF model (BP = 100%) but corresponding to an unconverged node in the CAT+GTR+Γ4 analysis (two chains supporting the named relationship, and two chains supporting a sister relationship to Telonemia + Sar; for details, see supplementary figs. S2 and S3, Supplementary Material online).

The relationships within Diaphoretickes show two strongly supported blocks, Telonemia + Sar and Cryptista + Archaeplastida (whose monophyly remains ambiguous), as well as a third group, Haptista, whose position is unclear (fig. 1). In order to test for the consistency of these results, we performed two alterations of the main 248 genes/109 OTUs supermatrix. First, we evaluated the persistence of the phylogenetic signal for the Telonemia + Sar and Cryptista + Archaeplastida groupings when progressively decreasing the number of genes used to infer the trees (fig. 2A). Random subsamples of the total 248 genes were generated without replacement in 20% decrements (80% to 20%), with a replication scheme that ensured a 95% probability of sampling every gene in each subset (Brown et al. 2018). Each replicate was then subjected to ML tree reconstruction under the LG+C60+F+Γ4 model and the group support monitored with ultrafast bootstrap (UFBoot) approximation (Hoang et al. 2018). Both groups (i.e., Telonemia + Sar and Cryptista + Archaeplastida) showed remarkable stability until as few as 60% of the genes were subsampled, corresponding to an alignment length of 149 genes (in average 34,588 AA positions), but started to crumble with fewer genes (fig. 2A). None of the other tested groups were as robust against the decreasing number of genes, with the exception of Sar (used as control) that remained maximally supported throughout the experiment (fig. 2A).

Fig. 2.

Open in new tabDownload slide

Gene subsampling and fast-evolving sites removal analyses. (A) UFBoot approximation for selected nodes of interest using the full supermatrix (100%) and subsets of randomly sampled genes (20–80%; see Materials and Methods). The variability of support values is shown by Box-and-Whisker plots. (B) UFBoot approximation UFBoot for the same nodes as in (A) using the full supermatrix (0 sites removed) and subsets from which fast-evolving sites were removed in 5,000 increments (see Materials and Methods). The gray box indicates the reduced data set which recovered a different position for Haptista. Arch, Archaeplastida; Crypt, Cryptista; Glau, Glaucophyta; Hapt, Haptista; Telo, Telonemia; Viri, Viridiplantae.

The second alteration was the removal of the fastest-evolving sites from the full supermatrix, which are more likely to be homoplasic and thus leading to systematic errors due to model misspecification (Rodriguez-Ezpeleta, Brinkmann, Roure, et al. 2007). Sites were sorted from fastest to slowest and progressively removed in 5,000 sites increments until only 3,469 sites were left. Each sub-data set was again used in ML tree reconstruction (LG+C60+F+Γ4 model) with UFBoot approximation (Hoang et al. 2018). The Telonemia + Sar relationship remained strongly supported until 25,000 sites were removed, after which the support was progressively lost but not replaced by an alternative supported relationship (fig. 2B). Similarly, the Cryptista + Archaeplastida group retained strong support until the removal of 35,000 sites with no other supported grouping emerging when more sites were removed. In contrast, the position of Haptista was markedly altered by the removal of fast-evolving sites. Starting from a maximally supported sister relationship to the Cryptista + Archaeplastida group (figs. 1 and 2B), Haptista shifted to a sister position to the Telonemia + Sar group after the removal of 25,000 sites (82% UFBoot), before reverting again to its original position with further site removal albeit with no real statistical significance (fig. 2B). In order to further assess the strength of the Telonemia + Sar + Haptista relationship, we ran a Bayesian analysis with the CAT+F81 model (Lartillot et al. 2013) using the 25,000 sites-removed data set. Interestingly, the two converged chains were fully consistent with the ML tree and the position of Haptista as sister to Telonemia + Sar received maximal PP (fig. 3 and supplementary fig. S4, Supplementary Material online).

Fig. 3.

Open in new tabDownload slide

Simplified Phylobayes consensus tree inferred from a reduced concatenated alignment (25,000 fast-evolving sites removed) of the 248 protein-coding genes under the CAT+F81+Γ4 model (for a full version showing all taxa, see supplementary fig. S4, Supplementary Material online). Node support is given by Bayesian PPs and UFBoot approximation inferred from a ML analysis (LG+C60+G+F model) using the same subalignment (i.e., 25,000 sites removed; see text and fig. 2). Nodes with black circles are fully supported in both analyses. A number appended to a black circle indicates full support in only one of the two analyses (i.e., PP/UFBoot). Note the complete congruence between the Bayesian and ML analyses (all nodes are recovered by both methods), including the position of Haptista.

Discussion

The position of Telonemia in the tree of eukaryotes has been debated ever since the group was first described (Vørs 1992; Klaveness et al. 2005; Shalchian-Tabrizi et al. 2006). The difficulty in placing Telonemia is due to a combination of several factors, including the ancient origin of the group, the sparse sampling of molecular data for group members, and the lack of closely related outgroup lineages. Earlier phylogenomic studies often placed Telonemia in the vicinity of Sar or close to cryptist or haptist lineages but never convincingly (Burki et al. 2009, 2012, 2016; Katz and Grant 2015; Lax et al. 2018). Here, we tested whether improving the gene coverage in our phylogenomic supermatrix and adding more telonemid representatives help to infer the origin of this group with more confidence. The results of the phylogenetic analyses show that this is the case: for the first time, the placement of Telonemia in the global phylogeny of eukaryotes receives strong and consistent support across multiple analyses (ML and Bayesian analyses using site-heterogeneous models, fig. 1). Telonemia was inferred as the most closely related lineage to the Sar supergroup. This position is robust against variations in the gene sets although it requires a substantial amount of data to be recovered (typically above 150 genes), and is stable when the most saturated sites are removed. It is thus unlikely to be due to convergence from undetected multiple substitutions (fig. 2). To acknowledge the strength of this relationship and facilitate further discussions, we introduce the constructed name “TSAR” to refer to the Telonemia + Sar grouping. We propose that TSAR forms a provisional supergroup of eukaryotes.

The hypothesis of TSAR has multiple implications. Most importantly, it confirms that Telonemia is a bona fide eukaryotic phylum, at the same phylogenetic level as other major phyla such as the Sar members (Stramenopila, Alveolata, and Rhizaria) (Burki et al. 2007; Hackett et al. 2007; Rodriguez-Ezpeleta, Brinkmann, Burger, et al. 2007). A unique combination of ultrastructural cellular characteristics suggested that Telonemia represents a deeply branching lineage without clear evolutionary affinities (Shalchian-Tabrizi et al. 2006; Yabuki et al. 2013). However, several morphological features, for example the peripheral vacuoles, extrusomes, or tubular tripartite hairs, were also suggested to resemble cell structures found in specific lineages such as Stramenopila, Alveolata, or Rhizaria (Shalchian-Tabrizi et al. 2006; Yabuki et al. 2013). Here, we show that Telonemia is indeed closely related to the Sar group, indicating that these morphological similarities might be homologous. If true, this could explain the origin of important synapomorphies in members of Sar from a common ancestor of the group, namely the alveoli in Alveolata from peripheral vacuoles and the mastigonemes in Stramenopila from similar flagellar hairs. Telonemia thus appears crucial for understanding the early evolution of Sar, a globally distributed and highly diverse assemblage of eukaryotes likely containing the majority of protist lineages (del Campo et al. 2014; de Vargas et al. 2015; Grattepanche et al. 2018). More generally, our study identifies Telonemia as crucial to infer the evolution of other important characteristics of eukaryotes, for example their cytoskeleton that underlies cell morphology and motility (Dacks et al. 2016). With its unique, highly complex multilayered cytoskeleton (Yabuki et al. 2013), it is possible that Telonemia has retained characteristics of an ancestral complex organization before the cytoskeleton was differentially reduced in other eukaryotic groups. A better understanding of the cytoskeleton evolution and other cell systems will emerge as more data are gathered from diverse eukaryotes, especially from deep-branching lineages like Telonemia.

In addition to TSAR, our analyses strongly recovered another major clade in Diaphoretickes, including Archaeplastida members (not consistently monophyletic) and Cryptista (AC clade; fig. 1). This association has been recovered before in taxon-rich phylogenomics, and its implication for eukaryote evolution discussed (Burki et al. 2016; Cenci et al. 2018). However, the AC clade still awaits confirmation because other recent studies displayed alternative relationships, albeit with no conclusive support and poorer taxon sampling for the relevant groups (Yabuki et al. 2014; Cavalier-Smith et al. 2015, 2018; Katz and Grant 2015). Our analyses also failed to consistently place Haptista in the tree. Based on the full data set, ML and some unconverged Bayesian chains favored a close relationship to the AC clade (fig. 1). But after the removal of the 25,000 fastest-evolving sites, both ML and converged Bayesian analyses supported the placement of Haptista and TSAR together (fig. 3). The position of Haptista has varied dramatically with conflicting confidence in support (e.g., Yabuki et al. 2014; Katz and Grant 2015; Burki et al. 2016; Cavalier-Smith et al. 2018; Cenci et al. 2018). Hence, moving forward to improve the deep structure of Diaphoretickes, it will be crucial to specifically address the causes for the lack of consensus for the position of Cryptista, Haptista, and the monophyly of Archaeplastida.

Despite these lasting uncertainties surrounding the relationships within Diaphoretickes, we present here consistent results that bring us one step closer to a fully resolved eukaryotic tree of life. Telonemia was once a notorious rogue lineage, but using new transcriptome data for this important but neglected phylum, we show that it now occupies a robust position in the tree as part of an expanded supergroup (TSAR). We thus advise that Telonemia should be routinely included in future works addressing the broad-scale evolution of eukaryotes, in the same way as other major lineages are included. More generally, this study follows a series of recent phylogenomic investigations that successfully placed in the tree several orphan taxa (Brown et al. 2013, 2018; Burki et al. 2016; Janouškovec et al. 2017; Lax et al. 2018). This demonstrates that obtaining adequate genomic data sets (including transcriptomics) for orphans can go a long way towards filling important evolutionary gaps. With the maturation of culture-free genomic solutions, we expect that inaccuracy due to poor taxon sampling will continue to reduce in phylogenomic investigations. Thus, even uncultured long-standing mystery taxa will soon fill the tree of life and contribute to a better understanding of ancient eukaryote evolution.

Materials and Methods

RNA Extraction, Sequencing, and Assembling

Clonal cultures of three telonemid species were isolated from Arctic surface waters of the Kara Sea, N 75.888, E 89.508, total depth: 52 m, water temperature: 2.3°C. Isolates were propagated on the bodonid Procryptobia sorokini Frolov, Karpov, and Mylnikov, 2001 strain B-69 grown in Schmalz-Pratt’s medium by using the bacterium Pseudomonas fluorescens as food (Tikhonenkov et al. 2014).

For cDNA preparation, cells were harvested following peak abundance after eating most of the prey. Cells were collected by centrifugation (2,000 × g, room temperature) on the 0.8 µm membrane of Vivaclear Mini columns (Sartorium Stedim Biotech Gmng, Germany, Cat. No. VK01P042). Total RNA was extracted using the RNAqueous-Micro Kit (Invitrogen, Cat. No. AM1931) and was converted into cDNA prior to sequencing using the SMART-Seq v4 Ultra Low Input RNA Kit (Takara Bio USA, Inc., Cat. No. 634888). Sequencing was performed on the Illumina MiSeq platform with read lengths of 300 bp using the NexteraXT protocol (Illumina, Inc., Cat. No. FC-131-1024) to construct paired-end libraries.

Reads were trimmed with Trimmomatic v. 0.36 (Bolger et al. 2014) using the following parameters: LEADING: 5, TRAILING: 5, SLIDINGWINDOW: 5: 16, MINLEN: 80. Potential vector sequences were detected (and removed if applicable) by BLAST searches against the UniVec database (e-value: 1e-4). Cleaned reads were assembled and translated into amino acids with Trinity v. 2.4.0 and TransDecoder, respectively (default settings Haas et al. 2013). In addition to the reads obtained in our study, publicly available reads of other taxa were trimmed, assembled, and translated as described above but with slightly different settings for the trimming procedure (LEADING: 3, TRAILING: 3, SLIDINGWINDOW: 4: 15, MINLEN: 36).

Phylogenomic Data Set Construction

A previously published data set comprising 263 genes and 234 taxa was used as starting point for this study (Burki et al. 2016). A broad range of new taxa were added by the following successive steps (sources of data: ensemblgenomes.org, imicrobe.us/#/projects/104, ncbi.nlm.nih.gov, onekp.com): 1) Protein sequences for each taxon were clustered with CD-HIT (Fu et al. 2012) using an identity threshold of 85%; 2) Homologous copies were retrieved by BLASTP searches using the 263 genes as queries (e-value: 1e-20; coverage cutoff: 0.5); 3) Phylogenetic trees were constructed and carefully inspected in order to detect and remove contaminants and select orthologous copies. Sequences were removed if their positions in gene trees conflicted with an initial species tree based on an intermediate concatenated alignment (i.e., after first round of cleaning, see below); this initial species tree globally corresponded to a consensus phylogeny of eukaryotes (e.g., Burki 2014). An automated coloring scheme was developed for easier visual inspection of the gene trees in FigTree v. 1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/; Last accessed February 6, 2019). Three rounds of inspection were performed, first to detect deep paralogs and then progressively refined to remove contamination and select orthologs. Sequences were aligned with MAFFT v. 7.310 (Katoh and Standley 2013) using the -auto option (first round) and with MAFFT L-INS-i using the default settings (second and third round), and filtered with trimAL v. 1.4 (Capella-Gutiérrez et al. 2009) using a gap threshold of 0.8 (all three rounds). FastTree v. 2.1.10 (Price et al. 2010) employing -lg -gamma plus options for more accurate performances (first round) and RAxML v. 8.2.10 (Stamatakis 2014) employing PROTGAMMALGF and 100 rapid BP searches (second and third round) were used to infer single gene ML trees. Three genes (FTSJ1—RNA methyltransferase 1, tubb—beta-tubulin, and eif5A—eukaryotic translation initiation factor 5A) were excluded from the base data set due to suspicious clustering of major groups that appeared after the inclusion of new taxa: in FTSJ1, a duplication of the entire Sar group was detected; in tubb, ambiguous relationships were observed due to duplications within Stramenopila and Opisthokonta; in eif5A, an unexpected clustering of Cryptista and Rhizaria was recovered. Twelve other genes were removed because they lacked telonemid sequences.

The initial data set comprised 248 genes and 293 taxa. An unaligned version of these genes was subjected to PREQUAL v. 1.01 (Whelan et al. 2018), a new software to remove or mask stretches of nonhomologous characters in individual sequences. A PP threshold of 0.95 (-filterthresh) was used. MAFFT G-INS-i with the VSM option (–unalignlevel 0.6) was then used for alignment, and ambiguously aligned sites were filtered out with BMGE v. 1.12 (Criscuolo and Gribaldo 2010) with the parameters -g 0.2, -b 5, -m, and BLOSUM75. Partial sequences belonging to the same taxon that did not show evidence for paralogy or contamination on the gene trees were merged. Of the initial 293 taxa, a reduced set of 148 taxa was selected to allow computationally intensive analyses based on 1) knowledge of sequence divergence (slower-evolving species were retained when possible), 2) to represent all major eukaryotic groups with sufficient data (amino acid positions in supermatrix >55%); of all major eukaryotic lineages with genomic or transcriptomic data, only Picozoa was deemed of insufficient coverage, and 3) sequence completeness among the genes. Furthermore, in order to reduce missing data, several strains and species complexes were combined when monophyly inferred from a preliminary concatenated alignment was unambiguous (supplementary fig. S1, Supplementary Material online; see also supplementary table 1, Supplementary Material online for a list of 109 taxa and chimera that were used as OTUs. Finally, all 248 genes were concatenated into a supermatrix (58,469 amino acid positions) with SCaFos v. 1.25 (Roure et al. 2007).

Phylogenomic Analyses

ML trees were inferred from the supermatrix using IQ-TREE v. 1.6.3 (Nguyen et al. 2015). Following the Akaike information criterion and the Bayesian information criterion, the best fitting model was the LG+C60 mixture model (Quang et al. 2008; Minh et al. 2013) with stationary amino acid frequencies optimized from the data set (F) and four discrete gamma (Γ) categories for rate heterogeneity (LG+C60+F+Γ4). The tested models were site-homogeneous LG +/− F, I, Γ4; site-heterogeneous LG +/− F, Γ4 and C10–C60. The best ML tree for the full data set was used as guide tree under the same model to estimate the PMSFs (Wang et al. 2018) model, which was used for further tree reconstructions (LG+C60+F + Γ4+PMSF). Support values for the best LG+C60+F + Γ4+PMSF tree were obtained using nonparametric BP with 100 replicates.

Bayesian analyses were inferred using PHYLOBAYES-MPI v. 1.8 (Lartillot et al. 2013) with the CAT+GTR+Γ4 model (Lartillot and Philippe 2004). Constant sites were removed to decrease running time (-dc). Four independent Markov Chain Monte Carlo (MCMC) chains were run for >6,300 generations (all sampled). For each chain, the burnin period was chosen after monitoring the evolution of the log-likelihood (Lnl) at each sampled point, removing the generations before stabilization (generally 1,250) prior to computing a 50% majority-rule consensus tree of all chains using bpcomp. Global convergence between chains was assessed by the maxdiff statistics measuring the discrepancy in PPs. Unfortunately, global convergence was never achieved in any combination of two chains; this is a known issue of PHYLOBAYES, which is especially acute with taxon-rich data sets and has been reported on many occasions (e.g., Katz and Grant 2015; Burki et al. 2016; Kang et al. 2017). The discrepancies between chains concerning nodes under scrutiny in this study were labeled as unresolved.

Gene Subsampling and Fast-Evolving Sites Removal

From the 248 gene data set, random subsampling without replacement of 80% to 20% (in 20% increments) of genes were performed with replications (3, 5, 6, and 14 replicates in each subset, respectively); this allowed a >95% probability that every gene is represented at least once (following an approach described in Brown et al. 2018). The genes in each subset were then concatenated with SCaFoS. Fast-evolving sites were estimated with IQ-TREE v. 1.6.5 (option -wsr) giving the PMSF tree as fixed tree topology and using the LG+C60+F + Γ4 model. Sites were then removed in 5,000 increments using a Perl script. For all the subsamples with reduced number of genes or sites, phylogenetic trees with UFBoot (1,000 replicates) were inferred using IQ-TREE (LG+C60+F+Γ4) with the -bb, -wbtl, and -bnni flags employed. One sub-alignment (25,000 sites removed) was used in a Bayesian tree inference using the site-heterogeneous CAT+F81+Γ4 model as implemented in PHYLOBAYES-MPI v. 1.8 (Lartillot et al. 2013). We did not use here the CAT+GTR + Γ4 model for computational tractability. Two independent MCMC chains were run for >6,300 cycles. Before computing the majority-rule consensus tree from the two chains combined, the initial 2,500 trees in each MCMC run were discarded as burnin after checking for convergence in likelihood and in clade PPs (maxdiff <0.3).

Supplementary Material

Supplementary data are available at Molecular Biology and Evolution online.

Acknowledgments

This work was supported by a grant from Science for Life Laboratory available to F.B., which covered salaries of J.F.H.S. and M.J., and experimental expenses. Microscopical identification, cell isolation and culturing, and generation of material for sequencing were supported by the Russian Science Foundation (grant no. 18-14-00239) available to D.V.T. D.V.T. thanks Dr Natalia Kosolapova (IBIW RAS) for help with sample collection. Computational resources were provided in part by the Swedish National Infrastructure for Computing through UPPMAX under project SNIC 2018/8-192.

Data deposition: All new transcriptomic data have been submitted to the SRA database under accession number SRP150904 (https://www.ncbi.nlm.nih.gov/sra/SRP150904; Last accessed February 6, 2019). All single gene alignments, gene trees, and the supermatrix are available online in supplementary file TSAR_Sequence_data_and_trees.zip.

References