Monophyly of gregarines and eugregarines, and polyphyly of apicomplexan parasites

Maximum likelihood analysis with the LG+C60+F+G4+PMSF model in IQ-TREE (non-parametric and UFBoot2 supports) and PhyloBayes analysis with the CAT+GTR model (ten independent runs) produced congruent topologies that were fully resolved at most internal branches (Figure 1A). The protococcidian Eleutheroschizon duboscqi was a sister of eucoccidians, and rhytidocystids branched as basal coccidiomorphs, the group containing most medically important apicomplexans. The blastogregarine Siedleckia was strongly related to the archigregarine Selenidium pygospionis. Gregarine apicomplexans (sensu lato including blastogregarines but excluding Digyalum; see below) and eugregarines were both monophyletic, in contrast to most interpretations based on ribosomal RNA genes (rDNA) (Cavalier-Smith, 2014; Leander, 2008; Rueckert et al., 2011). To test relationships between major apicomplexan lineages, we analyzed datasets in which the two longest branches in the tree, Cephaloidophora and Gregarina, were excluded either individually or both together (LG+C60+F+G4+PMSF model with UFBoot2 supports). The resulting trees had congruent topologies with all internal branches fully supported except for two eugregarine subclades and the position of cryptosporidians (Figure 1—figure supplement 1B). Similarly, seven statistical tests on 105 tree topologies representing all possible relationships between coccidiomorphs, cryptosporidians, eugregarines, archigregarines, and blastogregarines rejected all alternative topologies at p=0.01 except those differing in the placement of cryptosporidians (Supplementary file 3). The relationship of cryptosporidians therefore requires additional support, although their sister position to gregarines, which was unambiguously recovered in all trees including the ten PhyloBayes runs, is a preferred hypothesis. Two apicomplexan-like parasites branched outside the main apicomplexan clade (Figure 1A,B). Digyalum oweni, a formally described archigregarine was fully resolved as a sister lineage to all apicomplexans and chrompodellids. A previously undescribed parasitic symbiont of the annelid Scoloplos with apicomplexan-like traits, named ‘Symbiont X’, was specifically related to Chromera velia. Thus, apicomplexan parasites in the traditional sense are polyphyletic (see Discussion). We keep using the name ‘Apicomplexa’ for the clade of parasites sensu stricto hereinafter (Figure 1A).

Plastids in Digyalum and deep-branching apicomplexans, and their multiple losses

We next explored the existence of plastids in gregarines, Digyalum, and other parasites, where none have been known. Searching their sequence data with sequences of known plastid-localized proteins (Janouškovec et al., 2015; Ralph et al., 2004; Seeber and Soldati-Favre, 2010) revealed broad presence of plastidial pathways (Materials and methods; Supplementary file 4). Digyalum, Selenidium, Siedleckia, rhytidocystids, and Eleutheroschizon all contain near-complete plastidial biosynthesis of isoprenoid precursors, heme, and fatty acids, ferredoxin redox system, iron-sulfur cluster synthesis, and plastid genomes (Figure 2A). The eugregarine Lankesteria unusually contains fatty acid biosynthesis as the only plastidial pathway, whereas Symbiont X contains only the isoprenoid pathway, similar to piroplasms (Lizundia et al., 2009). Both Lankesteria and Symbiont X appear to lack plastid genomes (Figure 2A). The distribution of control, signature plastid genes involved in polypeptide import, folding, and DNA replication in the plastid (cpn60, sDer-1, PREX), matches the presence of plastid metabolism and genomes (Figure 2A). Maximum likelihood phylogenies of all individual proteins allowed us to readily distinguish the apicomplexan sequences from bacteria and other contaminants in the datasets (Materials and methods). In most phylogenies, the apicomplexan sequences cluster with algal plastid forms, confirming that they came from the plastid endosymbiont rather than the eukaryotic host. The phylogeny is different in several genes that are either derived by horizontal gene transfer from bacteria or in fact localize outside of the plastid in Plasmodium (in heme biosynthesis; see below). N-terminal regions of plastid sequences from our new transcriptomes often carry signal peptides typical for targeting to the plastid (the other have incomplete N-termini or lack targeting signatures by default, such as most triose phosphate translocators). The signal peptides are often followed by transit peptide-like regions, although these are more difficult to predict computationally (Supplementary file 5). In Digyalum, transit peptides have a net positive charge similar to other plastid leaders (they are low in acidic and high in basic residues; Figure 2—figure supplement 1A). Digyalum transit peptides are compositionally similar to transit peptides in Plasmodium, and likewise lack the phenylalanine motif at the first position after the signal peptide cleavage site (Figure 2—figure supplement 1B) (Patron and Waller, 2007). Predicted localizations for plastidial proteins and mitochondrial ALAS correspond closely to experimental evidence in Plasmodium and Toxoplasma (Figure 2A). The only exceptions to this pattern are the last three enzymes in heme biosynthesis, which are predicted to be plastidial in some parasites (see Discussion). The reconstructed plastid pathways also reflect known dependencies between their modules (Figure 2B). Iron-sulphur cluster assembly and ferredoxin system are widely required as co-factors for isoprenoid and fatty acid synthesis, whereas heme biosynthesis can be lost independently of other modules – likely the case in Lankesteria and Symbiont X. The synthesis of 3-phosphoglycerate (GAPDH-II and PGK-II) is present selectively. SufA was not identified in the Chromera genome. Pyruvate dehydrogenase and fatty and lipoic acid synthesis protein sequences in Digyalum are notably divergent, including three that are more closely related to bacterial than plastid sequences (Figure 2A). The cytosolic mevalonate pathway for isoprenoid precursor synthesis is absent in all species, but the mitochondrial cysteine desulphurase (IscS) and cytosolic fatty acid synthase (FASI) and elongase (ELO) pathways are present, as expected (Materials and methods) (Dellibovi-Ragheb et al., 2013; Ramakrishnan et al., 2012; Zhu et al., 2004). No plastid genes were identified in five eugregarine lineages other than Lankesteria, including in the draft genome of Gregarina (Figure 2A; see Discussion about the PREX fragment in Ascogregarina). This result is suggestive of at least two losses of plastids in eugregarines, which were independent of the one in Cryptosporidium (Discussion).

Open in a separate window

Figure 2.

Core plastid metabolism in apicomplexans and their relatives.

(A) Presence of genes and pathway modules (top; abbreviations in Supplementary file 4) in representative genomes (G) and transcriptomes (left. Each gene (box) is color-coded as to its evolutionary origin, as determined by a maximum likelihood phylogeny (plastid-encoded rps, rpl, rpo, trn genes were not analyzed). Empty boxes indicate gene absence in completed genomes and blank spaces indicate absence in transcriptomes. Intracellular localization of corresponding proteins is shown by a circle inside the box and summarizes known experimental data (Supplementary file 4) or de novo prediction in silico by SignalP v4.1 (Supplementary file 5); it is missing in proteins with incomplete N-termini. Note that only some enzymes of the heme pathway (HEM) are localized in the plastid (*1) and that signal peptides in FAS:ER and ELO were not predicted (*2). (B) Dependence network of plastid protein modules for the biosynthesis of key metabolites – isoprenoid precursors IPP and DMAP, fatty acids and heme – which underlie dependency on the plastid organelle in Apicomplexa. Colored regions contain modules specific to one pathway: fatty acid (pink), isoprenoid precursor (pale green) and heme biosynthesis (yellow). Interactions are reconstructed from the literature and substrates are shown near arrows (PYR = pyruvate; AcCoA = acetyl coenzyme A, Lip:E2 = lipoylation on PDHE2; C8:ACP = octanoyl:acyl carrier protein; [FE-S]=iron sulphur cluster; PEP = phosphoenolpyruvate; e-=electron reductive power; GA3p=glyceraldehyde-3-phosphate; 3PGA = 3 phosphoglycerate; ALA = δ-aminolevulinic acid; ?=uncertainty).

Figure 2—figure supplement 1.

Open in a separate window

Transit peptides in Digyalum plastid proteins have positive charge but lack a conserved phenylalanine in the first position.

(A) Amino acid composition of plastid transit peptides in Digyalum in comparison with downstream sequences of mature plastid proteins and transit peptides in other complex plastids as reported by Patron and Waller (2007) (only the first 14 residues were analyzed to enable consistent comparison) (Patron and Waller, 2007). Average proportions of individual amino acids and their classes based on charge (colored bars) are shown: acidic (D + E) in red, basic (H + K + R) in blue, polar uncharged (C + N + Q + S + T + W + Y) in green and hydrophobic (A + F + G + I + L + M + P + V) in black. Number of proteins analyzed are indicated in square brackets behind species/group names. (B) WebLogo3 diagram of the relative composition of the first 20 amino acids in 41 target peptides of Digyalum oweni. Letter size indicates the probability of an amino acid occurring at a specific position, weighted to account for the relative abundance of each amino acids natural occurrence and unrepresentative results arising from a small sample size.

Figure 2—figure supplement 2.

Open in a separate window

Maximum likelihood phylogeny of the triose phosphate translocator, TPT (LG+F+R7 model in IQ-TREE with 10000 UFBoot2 replicates;≥80 are shown; black dots indicate full support).

Apicomplexans, chrompodellids and Digyalum are shown on orange background and other eukaryotes on gray background. Selected characterized proteins are highlighted in bold.

Apicomplexan plastids are widespread, and their genomes are highly divergent

The discovery of plastids with genomes in deep-branching apicomplexans raises questions about whether they correspond to the undescribed apicomplexan-related lineages (ARLs) as defined by plastidial 16S rDNA (Janouškovec et al., 2012). Here, we discovered ten novel ARLs among environmental 16S rDNAs in GenBank and clustered 16S rDNAs obtained from the VAMPS database (Huse et al., 2014) by a phylogenetic sorting approach similar to the one used previously (Janouškovec et al., 2012) (Materials and methods). We also discontinue the use of two ARL-X and ARL-XI described recently (Mathur et al., 2018), which we find to be members of the Vitrella clade (ARL-I; Figure 3A) and bacterial contaminants, respectively (Materials and methods). This brings the total number of unidentified ARLs to 15, in addition to three ARLs represented by Chromera, Vitrella, and corallicolid clades (see Supplementary file 6 for reference sequences for all ARLs). A global, maximum likelihood phylogeny of 16S rDNAs of bacteria, plastids and ARLs (including representative sequences of known ARLs and all GenBank and VAMPS centroid sequences of novel ARLs) readily illustrates that many plastids in apicomplexans and their relatives are yet to be discovered (Figure 3A). It also shows that some new ARLs (ARL-XV, ARL-XVI, ARL-XVII) are comparatively diversified and abundant.

Open in a separate window

Figure 3.

Plastid genomes in apicomplexans and their relatives are widespread and highly divergent.

(A) Maximum likelihood phylogeny of plastid-encoded 16S ribosomal RNA genes (rDNA) reveals 10 novel apicomplexan-related lineages (ARL-IX and ARL-XII to ARL-XX; Supplementary file 6). Tree was computed with the best-fit TVMe+R5 model in IQ-TREE with UFBoot2 supports at branches (10000 replicates;≥80 are shown; black dots indicate 100 support). Environmental sequences (dark blue) are derived from GenBank or VAMPS (97% identity cluster centroids; numbers of reads are shown in square brackets where > 1). Plastid 16S rDNA transcripts of newly sequenced species are shown in red. Note that sequences in the tree vary greatly in their AT content and substitution rates, which can induce a misleading topology - deep relationships in the tree should therefore be interpreted with caution. The fast-evolving sequences of peridinin dinoflagellates were not included. (B) Extremely high AT content in rhytidocystid plastid genomes (arrows). AT content of representative species from part A and Balanophora laxiflora and B. reflexa parasitic plants (Su et al., 2019), all abbreviated to first three letters, is shown for 16S rDNA and plastid genomes. Plastid genomes in the newly sequenced species are only partially reconstructed from transcripts (red color; see Supplementary file 7). Altered genetic codes are indicated. (C) Euler diagram of plastid genome contents in apicomplexans, dinoflagellates (ancestral gene sets for each), Digyalum, Chromera, and Vitrella. Genes on the green background are associated solely with photosynthesis. Small RNA genes are not shown.

Despite that ARLs have now more than doubled in number, none of them corresponds to the plastid 16S rDNAs of Digyalum, Eleutheroschizon, Siedleckia, Selenidium or Rhytidocystis. Instead, the five genera have some of the fastest-evolving and most AT-rich 16S rDNAs of all apicomplexans, and they cluster with compositionally similar sequences of the more distantly related hematozoans (compare Figure 1A and Figure 3A). Such artificial grouping of highly divergent sequences is a well-known phylogenetic artifact, and deep-level relationships in the tree thus ought to be interpreted with caution. The 16S rDNAs of Digyalum, Eleutheroschizon, Siedleckia, Selenidium or Rhytidocystis were recovered among a set of AT-rich transcriptomic contigs, which encode other genes typical of apicomplexan plastid genomes (Supplementary file 7; Materials and methods). In Rhytidocystis species 1 and 2, the AT content of 16S rDNAs (84% and 86%, respectively) and plastidial transcripts (88% and 91%, respectively) is among the highest of all plastid genomes described to date (Figure 3B) (Su et al., 2019). Digyalum, Siedleckia, Selenidium, and Rhytidocystis sp. one use a non-canonical genetic code in which UGA encodes for tryptophan; UGA is absent in the fragmentary plastid DNA of Rhytidocystis sp. two and encodes for STOP codons in Eleutheroschizon (Figure 3B). The Digyalum plastid encodes six genes never identified in apicomplexan plastids and lacks nine genes they do contain; only one gene (rps18) was relocated to the nucleus in parallel in both lineages (Figure 3C).

Relationships and morphological transitions in early apicomplexan evolution

Key findings of multiprotein phylogenies are that eugregarines and gregarines are unequivocally monophyletic (Figure 1 and Figure 1—figure supplement 1). This broadly supports traditional morphological classifications (Grassé, 1953) in contrast to many recent proposals based on rDNA phylogenies, which regard both groups as polyphyletic (Cavalier-Smith, 2013). Indeed, protein sequences allow for building larger phylogenetic matrices and have more even substitution rates than rDNAs, some of which are notoriously divergent and phylogenetically unstable. Although the sampling of eugregarine diversity is incomplete, our phylogeny contains six of their seven main lineages at the superfamily level (Simdyanov et al., 2017) – one of them being Polyrhabdina, which is apparently not a lecudinoid (Figure 1A and unpublished data). The eugregarine monophyly provides the first unambiguous phylogenetic support for their two candidate synapomorphies: the ultrastructure of the epimerite, and the ultrastructure of epicytic crests (Simdyanov et al., 2017). It also partially resolves the little understood relationships between eugregarine superfamilies (Cavalier-Smith, 2014; Simdyanov et al., 2017; Simdyanov et al., 2015). The blastogregarine Siedleckia groups strongly with the archigregarine Selenidium. The sampling of both groups is incomplete, especially in the archigregarines, but they do share a combination of characteristics that are rare or absent in other apicomplexans, namely active bending and twisting movement in trophozoites, pellicular folds running along the body length in many species, and one or more layers of longitudinal microtubules underlying the pellicle (Schrével et al., 2016; Simdyanov et al., 2018). Both lineages also feed by myzocytosis, which is known in some relatives of apicomplexans (Foissner and Foissner, 1984; Mylnikov and Mylnikova, 2008) but not in eugregarines, cryptosporidians or coccidiomorphs. Because the blastogregarine ultrastructure and life cycle share additional similarities with coccidians (Simdyanov et al., 2018), confirming their phylogenetic position is important for reconstructing character evolution in Apicomplexa as a whole. The monophyly of gregarines sensu lato in our trees (Figure 1 and Figure 1—figure supplement 1) provides phylogenetic support for a tentative synapomorphy of this group, the forming of a gametocyst during the life cycle (Simdyanov et al., 2017), although this characteristic is absent in blastogregarines. The position of cryptosporidians is still not fully resolved, but the group is consistently recovered next to gregarines (Figure 1 and Figure 1—figure supplement 1). Eleutheroschizon is clearly important to understanding the origin of eucoccidians (Figure 1A) (Valigurová et al., 2015), although its classification as a protococcidian is apparently incorrect (unpublished data). Two Rhytidocystis species are positioned as a basal group in the class Coccidiomorpha (Figure 1A), which supports traditional views on their classification closer to coccidians (Levine, 1979; Porchet Hennere, 1972) rather than to gregarines (Cavalier-Smith, 2014). The bigger group that includes rhytidocystids and related parasites of annelids and molluscs (Kristmundsson et al., 2011) therefore provides a stepping stone for understanding the evolution of coccidiomorphs, including the majority of medically important apicomplexans.

Widespread maintenance but multiple losses of plastids in apicomplexa

Evidence for endosymbiotic genes with plastid targeting signals and plastid genomes shows that plastids are common in deep-branching apicomplexans, despite not having been previously recognized (Figure 2, Figure 3 and Supplementary file 5). Their plastid metabolic networks are similar to those in Toxoplasma and Plasmodium and generally serve to produce three essential metabolites: isoprenoid precursors for the synthesis of carotenoids, ubiquinone and prenyls; heme as a protein co-factor; and fatty acids as constituents of cellular lipids (Imlay and Odom, 2014; Ramakrishnan et al., 2012; van Dooren et al., 2012). Heme biosynthesis is partially localized in the cytosol and mitochondria, and fatty acids can be extended by non-plastidial elongases (Figure 2A) (Ramakrishnan et al., 2012; Zhu et al., 2004), but de novo biosynthesis of all three metabolites outside of the plastid appears to be absent. This creates dependence on the plastid for the production of isoprenoid precursors, heme, and fatty acids, which can be overcome only by metabolite uptake from the environment. Salvage of host isoprenoids, heme, and fatty acids is common in parasites, but obtaining them in sufficient amounts across the life cycle is apparently challenging. Most parasites retain all three plastidial pathways, and only two lineages have lost plastid organelles altogether (Figure 5): cryptosporidians (Abrahamsen et al., 2004), and the dinoflagellates Hematodinium and Amoebophrya (Gornik et al., 2015; John et al., 2019). In free-living organisms, plastids are always retained, perhaps because replacing their metabolism by salvage is even more difficult (Janouškovec et al., 2017). Most consistently required is the biosynthesis of isoprenoid precursors (Figure 5), which is the only indispensable plastidial pathway in piroplasms, the Plasmodium blood stage, and possibly in Symbiont X (Figure 2A) (Lizundia et al., 2009; Yeh and DeRisi, 2011). A fully unprecedented situation exists in the eugregarine Lankesteria, where fatty acid biosynthesis appears to be the only retained plastidial anabolic pathway (the transcriptome is fragmentary but enzymes for the biosynthesis of isoprenoid precursors and heme are absent). Interestingly, the existence of plastids in Selenidium and Rhytidocystis is also consistent with ultrastructural reports of multimembrane organelles in some of their species, which correspond to non-photosynthetic plastids in other apicomplexans by size, appearance, and their position near the nucleus (Figure 5) (Leander and Ramey, 2006; Porchet Hennere, 1972; Schrével, 1971; Schrével et al., 2016; Wakeman et al., 2014). Five eugregarines lack any evidence for plastid presence; these include Gregarina where a draft genome is available (Figure 2A). Since Lankesteria is phylogenetically nested among them, eugregarines must have lost plastids at least twice (Figure 5). Ascogregarina contains a fragment of PREX (Figure 2A and Figure 4A), but because the species lacks evidence for other endosymbiont-derived genes we conservatively consider it encodes for an ex-plastidial protein or is a contaminant. Our results indicate that the two plastid losses in eugregarines are independent of the one in cryptosporidians (Figure 1A), unlike previous assumptions (Cavalier-Smith, 2014; Toso and Omoto, 2007). This increases total known losses of plastid organelles in apicomplexans to three and in all eukaryotes to four (Figure 5).

Open in a separate window

Figure 5.

Plastid evolution in apicomplexans and their relatives.

Plastid-related characteristics (A–D), origin of parasitism, and predicted losses of photosynthesis, plastid genomes, and plastid organelles were mapped on the updated phylogeny. The phylogeny is fully resolved except one branch where more support is needed (*1). Predicted core plastid anabolic capabilities, presence of plastid genomes, and transmission electron microscopy evidence (TEM) for multimembrane organelles corresponding to plastids by their size, appearance, and position within cells are shown on the right.

Plastids in apicomplexans and their relatives were inherited vertically from a common ancestor

Plastids in Chromera and Vitrella share several unique characteristics with the apicomplexan plastid (Janouskovec et al., 2010) but their common origin was doubted in the past (Bodył et al., 2009) and other early endosymbiotic events could underlie plastid evolution in the broader group including dinoflagellates (Waller and Kořený, 2017). Testing whether plastids in apicomplexans, chrompodellids and Digyalum derive from the same endosymbiont is therefore relevant to understanding plastid evolution and frequency of plastid losses in general. We find three lines of evidence that support a common origin. Firstly, phylogenies of individual plastid genes repeatedly group Digyalum with apicomplexans and chrompodellids, and their topologies either reiterate nuclear phylogenies directly or are unresolved without consistently supporting alternatives (Figure 4A,B, Figure 4—figure supplement 2 and Figure 4—figure supplement 3). Secondly, unidentified lineages that carry plastid-encoded 16S rDNAs related to apicomplexans and chrompodellids (ARLs) are growing in number and diversity (Figure 3A). Some are apicomplexans (Janouškovec et al., 2013; Kwong et al., 2019) and others may be free-living, but they altogether support the idea that plastids in the lineage are widespread (Janouškovec et al., 2012). Finally, unique evolutionary innovations in the plastidial PREX, ispH, and hemH link together the plastids in apicomplexans, chrompodellids, and Digyalum (bacterial fabG and split rpoC2 further link plastids in some of them). It is difficult to see how endosymbiosis would move plastids in or out of the lineage without changing the distribution of these genes. Altogether, the evidence suggests that plastids in apicomplexans, chrompodellids, and Digyalum were inherited vertically from a common ancestor, are widely distributed in the group, and are most likely retained by default, particularly in free-living representatives (Janouškovec et al., 2015; Janouskovec et al., 2010).

Metabolic evolution in the plastid

Core plastid metabolism in Digyalum, Symbiont X, and apicomplexans has been remarkably conserved across long time scales and reveals only a few outstanding variations. In Plasmodium falciparum, triose phosphate sugars are imported by two triose phosphate translocators (TPTs) residing in the outermost (PfoTPT) and innermost (PfiTPT) plastid membranes, respectively (Mullin et al., 2006). Toxoplasma contains only one TPT phylogenetically corresponding to PfoTPT. Orthologs of the PfoTPT lack N-terminal targeting signatures and are ubiquitous in apicomplexans, chrompodellids, and Digyalum (Figure 2—figure supplement 2). PfiTPT orthologs are less common (Gile and Slamovits, 2014), but they possess N-terminal signal peptides in Chromera, Vitrella, and Rhytidocystis sp. two compatible with their targeting to the inner plastid membrane (Figure 2A and Figure 2—figure supplement 2). The split between the two forms therefore likely predated Apicomplexa and the loss of the PfiTPT form in the lineage leading to Toxoplasma is apparently a derived evolutionary state. One enzyme that processes the TPT substrate dihydroxyacetone phosphate (DHAP) is triose phosphate isomerase (TPI-II), and the failure to identify TPI-II in piroplasms once led to the proposition that their plastids import glyceraldehyde-3-phosphate, the TPI-II product (Fleige et al., 2010). We find that TPI-II, although highly divergent in piroplasms, is present in all apicomplexans with plastids, and it frequently possesses N-terminal signal peptides for plastid targeting (Figure 2A). This suggests that the import and conversion of DHAP is conserved across apicomplexan plastids.

Analysis of heme biosynthesis enzymes suggests that the pathway consistently starts in the mitochondrion, as in heterotrophic eukaryotes – the algal C5-pathway must have been lost prior to the Digyalum-Apicomplexa divergence. Delta-aminolevulinic acid is likely imported to the plastid and processed to coproporphyrinogen III (Ralph et al., 2004; Sato et al., 2004). The last three enzymes in some apicomplexans sequenced by us are unexpectedly predicted to be also plastidial, similarly to Chromera and Vitrella (Koreny et al., 2011). In Plasmodium, and perhaps also in Toxoplasma and Digyalum, the last three heme biosynthesis enzymes localize in the cytosol and mitochondria, where heme is most needed (Varadharajan et al., 2004). It would be interesting to explore these contrasting localization predictions experimentally, including the possibility that some enzymes are dually targeted, or their isoforms (where present; Figure 2A) are differentially targeted to different cellular compartments.

Highly divergent plastid genomes and their multiple losses in chrompodellids

Unlike metabolism in the plastid, plastid genome structure shows unexpected variations across Apicomplexa. Plastid genomes in some of the newly sequenced taxa, as partially reconstructed from transcripts, have unusually AT-rich and fast-evolving sequences (Figure 3A). Completing the fragmentary plastid genome of Rhytidocystis species 2 (Supplementary file 7) is of particular interest because it has the most AT-rich 16S rDNA, and potentially the most AT-rich genome, ever recorded among plastids (Figure 3B) (Su et al., 2019). Plastid genome reduction in Digyalum and Apicomplexa from their common ancestor is primarily explained by the loss of photosynthesis (Figure 3C). This was accompanied by relatively modest transfer of genes to the nucleus, which involved remarkably different gene sets in the two lineages (only rps18 was relocated in both in parallel).

Underlying the very existence of plastid DNA in Plasmodium and Toxoplasma is sufB, which encodes one of only two broadly conserved, plastid-encoded proteins with function other than transcribing and translating the genome itself. Of the newly sequenced species with plastids, only Symbiont X lacks any evidence for the sufB gene or plastid genome (Figure 3A). It would be expected that Symbiont X sufB is encoded in the nucleus and was relocated there in its common ancestor with Chromera, Voromonas, Colpodella, and Alphamonas (Janouškovec et al., 2015), thus allowing heterotrophs in this lineage to lose plastid genomes at least three times independently (Figure 5).

Parasite sampling, transcriptome sequencing and assembly

Parasite cells (1 to approximately 70 individuals) were isolated from marine annelid, mollusc and barnacle hosts collected at the White Sea Biological Station of Moscow State University during August of 2016 (Supplementary file 1). Cells were hand-picked by using a glass micropipette (30% ethanol was used to detach the cells of Digyalum and Eleutheroschizon), washed 1 to 4 times in clean seawater and transferred into a clean tube in the final volume of 2–5 ul of seawater (excess seawater removed if necessary after a brief spin). Care was taken to avoid contamination by animal host cells and gut contents when isolating parasite cells but we were not able to fully prevent it (host contaminants were observed in the sequence data but they were clearly distinguishable from the parasites; see below). A 100 ul of Lysis buffer (RNAqueous-Micro Total RNA Isolation Kit; Ambion/ThermoFisher, cat no. AM1931) was added into the tube with parasite cells and the samples were stored at −80C for several weeks. Total RNA was extracted from samples by RNAqueous-Micro Total RNA Isolation Kit according to manufacturer instructions but without the DNase I digest step. In Digyalum and Eleutheroschizon, RNA from two independent cell isolations was combined before further processing (Supplementary file 1). Cells of Digyalum oweni WS3-2017 was stored in RNALater and used for reverse-transcription without extracting RNA. RNA was reverse-transcribed by using SMART-Seq v4 Ultra Low Input RNA Kit (Takara; cat no. 634888), however20-30 amplification cycles; optimization was done as in SMARTer Pico PCR cDNA Synthesis Kit). Indexed TruSeq libraries were built at Edinburgh Genomics and sequenced as paired-end 150 bp reads in two multiplexed lanes on the Illumina HiSeq 4000 machine. For the Digyalum oweni WS3-2017 sample, paired-end 150 bp HiSeq 4000 reads were produced independently of other samples. Demultiplexed reads were processed by cutadapt v1.8.3 to remove low quality nucleotides (-q 20 setting), polyA tails, SMARTer adapters and leftover Illumina barcodes (on 3' ends of shorter sequence fragments). Reads of minimum length of 100 nucleotides were assembled in Trinity v2.4.0 by using the default settings. Predicted proteomes were generated as six-frame translations of transcriptomic contigs. Analysis of assembled reads revealed minor cross-contamination between samples due to adapter swapping, but contaminants were well distinguishable by read coverage. Raw sequence reads and our assemblies were deposited in NCBI Bioprojects PRJNA557242 and PRJNA556465. Photographs of studied parasites (Figure 1B) were made using Leica DM 2500 microscopes equipped with DIC optics and a digital camera (Leica, Germany) at the White Sea Biological Station of Moscow State University and the Marine Biological Station of St. Petersburg State University.

Single protein sequence alignments for nuclear phylogeny

We built our dataset from 339 protein sequence alignments previously used in stramenopile phylogenies and tested for orthology (Derelle et al., 2016). Firstly, slow-evolving stramenopile sequence query for each alignment (mostly the Phytophthora sequence) was used to retrieve five best BLASTP hits (e-value cutoff of 1e-5) and select the closest ortholog in translated alveolate, stramenopile, and rhizarian genomes (Supplementary file 2). Secondly, a new alveolate protein sequence query (primarily from Vitrella) was used in the same way to expand the dataset with our newly generated and existing transcriptomes. Closest hits from animal, fungal, microsporidian and bodonid genomes were also included to distinguish sequences derived from animal host and other contaminants in the samples. Orthologous sequences were identified by multiple rounds of alignment clean-up as guided by Maximum Likelihood phylogenies. Initial rounds of phylogenies were focused on removing out-paralogs and contaminants by using default alignment in MAFFT v7.402 (Katoh and Standley, 2013), alignment trimming in BMGE v1.12 (Criscuolo and Gribaldo, 2010) with the -b 3 g 0.4 setting and phylogeny in Fasttree v2.1.10 (Price et al., 2010) with the -lg -gamma setting. Later rounds were focused on resolving difficult cases of paralogy, horizontal gene transfer, and the orthology of fast-evolving sequences: localpair, --linsi, alignment in MAFFT; -b 4 g 0.4 settings in BMGE; and IQ-TREE v1.6.5 (Nguyen et al., 2015) phylogeny with the LG+I+G4+F model and 1000 UFBoot2 supports. Of the original 339 single protein sequence alignments (see Supplementary file 1 in Derelle et al., 2016), 43 were excluded in the process: 10 were absent in apicomplexans and dinoflagellates (alignments #72, #290, #291, #292, #294, #295, #300, #301, #304 and #315) and 33 were sparsely sampled or strongly incongruent with known organismal relationships (alignments #29, #57, #70, #73, #82, #149, #151, #155, #157, #159, #163, #164, #167, #171, #183, #184, #191, #193, #194, #198, #241, #246, #260, #265, #270, #283, #285, #289, #302, #303, #305, #330, #331). In the remaining 296 alignments, paralogs were reduced to the most slowly evolving sequence with non-conflicting phylogenetic placement (else both were removed) and multiple gene isoforms were reduced to the most complete sequence. Cross-contaminant sequences were identified based on read coverage and excluded. Divergent regions of sequences resulting from frame shift errors or imperfect gene models were removed from alignments manually. Adjacent protein sequence fragments that were apparently derived from the same gene fragmented by wrong genome annotation were merged.

Multiprotein phylogeny and topology tests

The final 296 protein sequence datasets were realigned by the localpair algorithm in MAFFT, trimmed by using the -b 4 g 0.4 settings in BMGE and concatenated in Scafos v1.25 (Roure et al., 2007). During the latter step, sequences derived from different strains of the same species (in Colpodella angusta, Sarcocystis neurona, Voromonas pontica, and unidentified actinocephaloid parasites of Helicoverpa armigera and Helicoverpa assulta) were merged into single operational taxonomic units (Supplementary file 2). The initial phylogenetic matrix contained 54 species, 99948 sites and 14.4% missing data (Figure 1—figure supplement 1A). To further reduce missing data we merged sequences of two Oxyrrhis marina strains and of two distinct variants found in Siedleckia nematoides transcriptomes; the latter are possibly derived from different cryptic species (variant one was preferentially identified and/or found to be more complete in the transcriptome of the WS1 strain; variant two in the transcriptome of the WS2 strain). We also merged sequences of three representatives of the superfamily Actinocephaloidea with low sequence presence: Ascogregarina and two unidentified parasites of the insects Teleopsis and Helicoverpa, which contaminate the host transcriptomes (Borner and Burmester, 2017). The 18S rDNA of the Teleopsis parasite branches among other actinocephaloids (data not shown) and both unidentified group with Ascogregarina in the multiprotein phylogeny (Figure 1—figure supplement 1A). Merging all these taxa produced the main phylogenetic matrix with 50 species, 99908 positions and 10.6% missing data (Figure 1—source data 1). Three additional matrices were created by excluding all sequences of Gregarina, Cephaloidophora, or both species (Figure 1—figure supplement 1B). Final matrices were first analyzed in IQ-TREE by using the LG+I+G4+F model and the best tree was used as a guide tree in a more thorough analysis with the LG+G4+F+C60+PMSF model: 1000 UFBoot2 replicates were computed for all datasets (Figure 1A and Figure 1—figure supplement 1) and 100 non-parametric bootstrap were computed for the main matrix (Figure 1A). Seven statistical tests of 105 tree topologies corresponding to all possible relationships between coccidiomorphs, cryptosporidians, eugregarines, archigregarines, and blastogregarines in Figure 1A were calculated in IQ-TREE v1.6.5 with the LG+I+G4+F model and 10000 replicates (Supplementary file 3). The main phylogenetic matrix (Figure 1A) was also analyzed by 10 independent PhyloBayes runs (either standard version 4.1 c or MPI version 1.7b) (Lartillot et al., 2009) with the GTR+CAT model and constant sites removed (-dc setting). Each chain was run for 1000 cycles, of which the initial 250 were discarded and the remaining (10 × 750) were combined to compute a consensus tree (maxdiff = 0.23, meandiff = 0.00057).

Analysis of plastidial metabolism and horizontally acquired genes

We next searched the orthologs of apicomplexan plastid proteins with focus on core pathways (Figure 2A) in transcriptomes and genomes of apicomplexans and their relatives (including Perkinsus and dinoflagellates) by using similar approach as for the nuclear genes above. Five best BLASTP hits at the e-value threshold of 1e-5 were retrieved for each species by using a comparatively slow evolving sequence query (typically Vitrella or Chromera). These hits were included in plastid protein alignments created previously (Janouškovec et al., 2017; Janouškovec et al., 2015), or in new alignments created here by retrieving representative outgroup sequences from GenBank and the local database by BLASTP searches with the same query sequences (Supplementary file 4). This process included non-plastidial genes of the heme biosynthesis pathway, and three genes representing controls for non-plastidial pathways: iscS (mitochondrial iron-sulfur biosynthesis),enoyl reductase domain (FASI, cytosolic fatty acid synthesis), and ELO (endoplasmic reticulum-localized fatty acid elongation). Alignments were reduced by an iteration process analogous to that used for nuclear genes. Final phylogenetic matrices were prepared by localpair alignment in MAFFT and -b 4 g 0.4 trimming in BMGE, and analyzed in IQ-TREE by using the built-in ModelFinder to select the best model with the preset LG substitution matrix (for example, LG+F+R7 model for TPT in Figure 2—figure supplement 2). This approach allowed to distinguish true apicomplexan sequences from contaminants in our transcriptomes and confirmed their origins to be in the plastid endosymbiont (grouping with other plastidial sequences) or eukaryotic host (grouping with eukaryotic cytosolic or mitochondrial sequences; Figure 2A). Protein phylogenies of six horizontally acquired genes (Figure 4A,B, Figure 4—figure supplement 2 and Figure 4—figure supplement 3) were built by expanding previous datasets (Janouškovec et al., 2015) and computed by using the same method as plastid phylogenies above.

Analysis of leader sequences in plastid proteins

Proteins of core plastidial pathways (Figure 2A) derived from our transcriptomes and those of Plasmodium falciparum, Babesia microti and Toxoplasma gondii were scanned for the presence of identifiable N-terminal signal peptides in SignalP v4.1 (Petersen et al., 2011). All methionines downstream of the predicted protein start were also tested for targeting signals. Positive sequences were checked for the presence of N-terminal extensions – truncated proteins were filtered out as false positives. The remaining signal peptide-positive proteins were recorded in Figure 2A and further checked for identifiable transit peptides in ChloroP 1.1 (Emanuelsson et al., 1999). The SignalP and ChloroP statistics were listed in Supplementary file 5. Apicomplexan proteins that have known experimental localization (as described in primary literature and Apiloc3) or that were classified as high-confidence plastid-localized proteins in Plasmodium falciparum by BioID (Boucher et al., 2018), were recorded in Figure 2A and listed in Supplementary file 4. To characterize transit peptides in Digyalum oweni, we expanded its identified plastid protein set. Additional plastidial proteins in the WS1+two isolates were searched by BLASTP with apicomplexan queries, comprising known Toxoplasma and Plasmodium plastid proteins primarily involved in transcript and protein processing (e.g., stromal processing peptidase, clpP chaperone, histone-like protein HU, tRNA-Met ligase, etc.). Matched sequences were verified by reverse BLASTP searches on KEGG and NCBI BLAST websites (aided by distance tree phylogenies at the latter site). Pooling positive hits with those identified previously (Figure 2A) and filtering for sequences with N-terminal extension carrying a signal peptide (as identified by SignalP v4.1) produced 41 plastid proteins. The first 14 transit peptide residues downstream of the signal peptide cleavage site (an approach to allow comparison with results in Patron and Waller, 2007 were merged and analyzed altogether for their amino acid composition by ‘Protein Stats’ script at the Sequence Manipulation Suite website. Composition of the 41 mature proteins (without transit peptides) was also analyzed – because transit peptide cleavage sites are difficult to predict in silico the initial 50 residues were removed arbitrarily (Figure 2—figure supplement 1A). The amino acid frequency at the first 20 positions across the 41 Digyalum transit peptides was analyzed at the WebLogo3 website (Figure 2—figure supplement 1B).

We thank Ross F Waller and reviewers and editors in Elife for providing helpful comments and suggestions about the draft manuscript. We thank the staff of the White Sea Biological Station of Lomonosov Moscow State University and the Marine Biological Station of St. Petersburg State University for assistance during field sampling. The authors acknowledge the use of the UCL Myriad High Performance Computing Facility (Myriad@UCL), and associated support services, in the completion of this work. Digyalum oweni WS3 transcriptome sequencing was funded by the Russian Science Foundation (grant no. 18-14-00123).