The genomes of polyextremophilic cyanidiales contain 1% horizontally transferred genes with diverse adaptive functions

The role and extent of horizontal gene transfer (HGT) in eukaryotes are hotly disputed topics that impact our understanding of the origin of metabolic processes and the role of organelles in cellular evolution. We addressed this issue by analyzing 10 novel Cyanidiales genomes and determined that 1% of their gene inventory is HGT-derived. Numerous HGT candidates share a close phylogenetic relationship with prokaryotes that live in similar habitats as the Cyanidiales and encode functions related to polyextremophily. HGT candidates differ from native genes in GC-content, number of splice sites, and gene expression. HGT candidates are more prone to loss, which may explain the absence of a eukaryotic pan-genome. Therefore, the lack of a pan-genome and cumulative effects fail to provide substantive arguments against our hypothesis of recurring HGT followed by differential loss in eukaryotes. The maintenance of 1% HGTs, even under selection for genome reduction, underlines the importance of non-endosymbiosis related foreign gene acquisition.


Introduction
Eukaryotes transmit their nuclear and organellar genomes from one generation to the next in a vertical manner. As such, eukaryotic evolution is primarily driven by the accumulation, divergence (e.g., due to mutation, insertion, duplication), fixation, and loss of gene variants over time. In contrast, horizontal (also referred to as lateral) gene transfer (HGT) is the inter-and intraspecific transmission of genes from parents to their offspring. HGT in Bacteria (Doolittle, 1999;Ochman et al., 2000;Boucher et al., 2003) and Archaea (Nelson-Sathi et al., 2012) is widely accepted and recognized as an important driver of evolution leading to the formation of pan-genomes (Tettelin et al., 2005;Vernikos et al., 2015). A pan-genome comprises all genes shared by any defined phylogenetic clade and includes the so-called core (shared) genes associated with central metabolic processes, dispensable genes present in a subset of lineages often associated with the origin of adaptive traits, and lineage-specific genes (Vernikos et al., 2015). This phenomenon is so pervasive that it has been questioned whether prokaryotic genealogies can be reconstructed with any confidence using standard phylogenetic methods (Philippe and Douady, 2003;Doolittle and Brunet, 2016). In contrast, as eukaryote genome sequencing has advanced, an increasing body of data has pointed towards the existence of HGT in these taxa, but at much lower rates than in prokaryotes (Danchin, 2016). The frequency and impact of eukaryotic HGT outside the context of endosymbiosis and pathogenicity however, remain hotly debated topics in evolutionary biology. Opinions range from the existence of ubiquitous and regular occurrence of eukaryotic HGT (Husnik and McCutcheon, 2018) to the almost complete dismissal of any eukaryotic HGT outside the context of endosymbiosis as being Lamarckian, thus false, and resulting from analysis artefacts (Martin, 2018;Martin, 2017). HGT skeptics favor the alternative hypothesis of differential loss (DL) to explain the current data. DL imposes strict vertical inheritance (eukaryotic origin) on all genes outside the context of pathogenicity and endosymbiosis, including putative HGTs. Therefore, all extant genes have their root in LECA, the last eukaryotic common ancestor. Patchy gene distributions are the result of multiple ancient paralogs in LECA that have been lost over time in some eukaryotic lineages but retained in others. Under this view, there is no eukaryotic pan-genome, there are no cumulative effects (e.g., the evolution of eukaryotic gene structures and accrual of divergence over time), and therefore, mechanisms for the uptake and integration of foreign DNA in eukaryotes are unnecessary.
A comprehensive analysis of the frequency of eukaryotic HGT was recently done by Ku and Martin (2016). These authors reported the absence of eukaryotic HGT candidates sharing over 70% protein identity with their putative non-eukaryotic donors (for very recent HGTs, this figure could be as high as 100%). Furthermore, no continuous sequence identity distribution was detected for HGT candidates across eukaryotes and the 'the 70% rule' was proposed ('Coding sequences in eukaryotic genomes that share more than 70% amino acid sequence identity to prokaryotic homologs are most likely assembly or annotation artifacts.') (Ku and Martin, 2016). However, as noted by others (Richards and Monier, 2016;Leger, 2018), this result was obtained by categorically dismissing all eukaryotic HGT singletons located within non-eukaryotic branches as assembly/annotation artefacts, as well as those remaining that exceeded the 70% threshold. In addition, all genes that were presumed to be of organellar origin were excluded from the analysis, leaving a small dataset extracted from already under-sampled eukaryotic genomes.
Given these uncertainties, the aim of our work was to systematically analyze eukaryotic HGT using the Cyanidiales (known as Cyanidiophytina in some taxonomic schemes) as model organisms. The Cyanidiales comprise a monophyletic clade of polyextremophilic, unicellular red algae (Rhodophyta) that thrive in acidic and thermal habitats worldwide (e.g., volcanoes, geysers, acid mining sites, acid rivers, urban wastewaters, geothermal plants) (Castenholz and McDermott, 2010). With a divergence age estimated to be around 1.92-1.37 billion years , the Cyanidiales are the earliest split within Rhodophyta and define one of the oldest surviving eukaryotic lineages. They are located near the root of the supergroup Archaeplastida, whose ancestor underwent the primary plastid endosymbiosis with a cyanobacterium that established photosynthesis in eukaryotes (Reyes-Prieto et al., 2007;Price et al., 2012). In the context of HGT, the Cyanidiales became more broadly known after publication of the genome sequences of Cyanidioschyzon merolae 10D (Matsuzaki et al., 2004;Nozaki et al., 2007), Galdieria sulphuraria 074W (Schö nknecht et al., 2013), and Galdieria phlegrea DBV009 . The majority of putative HGTs in these taxa was hypothesized to have provided selective advantages during the evolution of polyextremophily, contributing to the ability of Galdieria, Cyanidioschyzon, and Cyanidium to cope with extremely low pH values, temperatures up to 56˚C, as well as high salt and toxic heavy metal ion concentrations (Castenholz and McDermott, 2010;Doemel and Brock, 1971;Reeb and Bhattacharya, 2010;Hsieh et al., 2018). In such environments, they can represent up to 90% of the total biomass, competing with specialized Bacteria and Archaea (Seckbach, 1972), although some Cyanidiales strains also occur in more temperate environments Gross et al., 2002;Ciniglia et al., 2004;Barcytė et al., 2018;Iovinella et al., 2018). The integration and maintenance of HGT-derived genes, in spite of strong selection for genome reduction in these taxa (Qiu et al., 2015) underlines the potential ecological importance of this process to niche specialization (Schö nknecht et al., 2013;Qiu et al., 2013;Raymond and Kim, 2012;Bhattacharya et al., 2013;Foflonker et al., 2018;Schö nknecht et al., 2014). For this reason, we chose the Cyanidiales as a model lineage for studying eukaryotic HGT.
It should be appreciated that the correct identification of HGTs based on sequence similarity and phylogeny is rarely trivial and unambiguous, leaving much space for interpretation and erroneous assignments. In this context, previous findings regarding HGT in Cyanidiales were based on single genome analyses and have therefore been questioned (Ku and Martin, 2016).
Many potential sources of error need to be excluded during HGT analysis, such as possible bacterial contamination in the samples, algorithmic errors during genome assembly and annotation, phylogenetic model misspecification, and unaccounted for gene paralogy (Richards and Monier, 2016). In addition, eukaryotic HGT reports based on single gene tree analysis are prone to misinterpretation and may be a product of deep branching artefacts and low genome sampling. Indeed, false claims of prokaryote-to-eukaryote HGT have been published (Boothby et al., 2015;Crisp et al., 2015) which were later corrected (Koutsovoulos et al., 2016;Salzberg, 2017).
Here, we used multi-genomic analysis with 13 Cyanidiales lineages (including 10 novel, long-read, draft genome sequences) from nine geographically isolated habitats. This approach increased phylogenetic resolution within Cyanidiales to allow more accurate assessment of HGT while avoiding many of the above-mentioned sources of error. The following questions were addressed by our research: (i) did HGT have a significant impact on Cyanidiales evolution? (ii) Are previous HGT findings in the sequenced Cyanidiales genomes an artefact of short read assemblies, limited genome databases, and uncertainties associated with single gene trees, or do they hold up with more extensive sampling? (iii) And, assuming that evidence of eukaryotic HGT is found across multiple Cyanidiales species, are cumulative effects observable, or is DL the better explanation for these results? Table 1. Summary of the 13 analyzed Cyanidiales genomes.
Genome Size (

Features of the newly sequenced cyanidiales genomes
Genome sizes of the 10 targeted Cyanidiales (Figure 1) range from 12.33  Mbp, similar to other members of this red algal lineage (Matsuzaki et al., 2004;Schö nknecht et al., 2013;Qiu et al., 2013) (Table 1). PacBio sequencing yielded 0.56 Gbp -1.42 Gbp of raw sequence reads with raw read N50 ranging from 7.9 kbp -14.4 kbp, which translated to a coverage of 28.91x -70.99x at the unitigging stage (39.46x -91.20x raw read coverage) (Appendix 1). We predicted a total of 61,869 novel protein coding sequences which, together with the protein data sets of the already published Cyanidiales species (total of 81,682 predicted protein sequences), capture 295/ 303 (97.4%) of the highly conserved eukaryotic BUSCO dataset. Each species, taken individually, scored an average of 92.7%. In spite of massive gene losses observed in the Cyanidiales (Qiu et al., 2015), these results corroborate previous observations that genome reduction has had a minor influence on the core eukaryotic gene inventory in free-living organisms (Qiu et al., 2016). Even C. merolae Soos, the species with the most limited coding capacity (4406 protein sequences), includes 89.5% of the eukaryotic BUSCO dataset. The number of contigs obtained from the Galdieria genomes ranged between 101-135. G. sulphuraria 17.91 (a strain different from the ones sequenced) was reported to have 40 chromosomes, and strains isolated from Rio Tinto (Spain), 47 or 57 chromosomes (Moreira et al., 1994). Pulsed-field gel electrophoresis indicates that G. sulphuraria 074W has approximately 42 chromosomes that are between 100 kbp and 1 Mbp in size (Weber, 2007). The genome assembly of C. merolae Soos produced 35 contigs, which approximates the 22 chromosomes (including plastid and mitochondrion) of the C. merolae 10D telomere-to-telomere assembly. Whole genome alignments indicate that a portion of the assembled contigs represent complete chromosomes.

Orthogroups and phylogeny
The 81,682 predicted protein sequences from all 13 genomes clustered into a total of 9075 orthogroups and phylogenetic trees were built for each orthogroup. The reference species tree was constructed using 2,090 OGs that contained a single-copy gene in at least 12 of the 17 taxa (Porphyra umbilicalis (Brawley et al., 2017), Porphyridium purpureum , Ostreococcus tauri RCC4221 (Blanc-Mathieu et al., 2014), and Chlamydomonas reinhardtii (Merchant et al., 2007) were added to the dataset as outgroups). As a result, the species previously named G. sulphuraria Soos and C. merolae MS1 were reannotated as G. phlegrea Soos and G. sulphuraria MS1. Given these results, we sequenced a second genome of C. merolae and a representative of the G. phlegrea lineage. The species tree reflects previous findings that suggest more biodiversity exists within the Cyanidiales ) ( Figure 2).

Analysis of HGTs
The most commonly used approach to identify HGT candidates is to determine the position of eukaryotic and non-eukaryotic sequences in a maximum likelihood tree. Using this approach, 96 OGs were identified in which Cyanidiales genes shared a monophyletic descent with prokaryotes, representing 1.06% of all OGs. A total of 641 individual Cyanidiales sequences are considered as HGT candidates ( Table 1). The amount of HGT per species varied considerably between members of the Cyanidioschyzon (33-34 HGT events, all single copy genes) and Galdieria lineages with 44-54 HGT events (52.6 HGT origins on average, 47-62 HGT gene candidates). In comparison to previous studies (Schö nknecht et al., 2013;Qiu et al., 2013), no evidence of massive gene family expansion regarding HGT genes was found because the maximum number of gene copies in HGT orthogroups was three. We note, however, that one large gene family of STAND-type ATPases that was previously reported to originate from an archaeal HGT (Schö nknecht et al., 2013) did not meet the criteria used in our restrictive Blast searches; that is the 10 À5 e-value cut-off for consideration and a minimum of three different non-eukaryotic donors. This highly diverged family requires more sophisticated comparative analyses that were not done here (Appendix 2).
Gene co-localization on raw sequence reads One major issue associated with previous HGT studies is the incorporation of contaminant DNA into the genome assembly, leading to incorrect results (Boothby et al., 2015;Crisp et al., 2015;Koutsovoulos et al., 2016;Salzberg, 2017). Here, we screened for potential bacterial contamination in our tissue samples using PCR analysis of extracted DNA with the rbcL and 18S rRNA gene markers prior to sequencing. No instances of contamination were found. Furthermore, our work relied on PacBio RSII long-read sequencing technology, whereby single reads frequently exceed 10 kbp of DNA. Given these robust data, we also tested for co-occurrence of HGT gene candidates and 'native' genes in the same read. The protein sequences of each species were queried with tblastn (10 À5 e-value, 75 bitscore) against a database consisting of the uncorrected PacBio RSII long reads. This analysis showed that 629/641 (98.12%) of the HGT candidates co-localize with native red algal genes on the same read (38,297 reads in total where co-localization of native genes and HGT candidates was observed). It should be noted that the 10 novel genomes we determined share HGT candidates with C. merolae 10D, G. sulphuraria 074W, and G. phlegrea DBV009, which were sequenced in different laboratories, at different points in time, using different technologies, and assembly pipelines. Hence, we consider it highly unlikely that these HGT candidates result from bacterial contamination. As the accuracy of long read sequencing technologies further increases, we believe this criterion for excluding bacterial contamination provides an additional piece of evidence that should be added to the guidelines for HGT discovery (Richards and Monier, 2016). . Species tree of the 13 analyzed extremophilic Cyanidiales genomes using mesophilic red (Porphyra umbilicalis, Porphyridium purpureum) and green algae (Ostreococcus tauri, Chlamydomonas reinhardtii) as outgroups. IQTREE was used to infer a single maximum-likelihood phylogeny based on orthogroups containing single-copy representative proteins from at least 12 of the 17 taxa (13 Cyanidiales + 4 Others). Each orthogroup alignment represented one partition with unlinked models of protein evolution chosen by IQTREE. Consensus tree branch support was determined by 2000 rapid bootstraps. All nodes in this tree had 100% bootstrap support, and are therefore not shown. Divergence time estimates are taken from Yang et al. (2016). Similarity is derived from the average one-way best blast hit protein identity (minimum protein identity threshold = 30%). The minimal protein identity between two strains was 65.4%, measured between g. sulphuraria SAG21.92, which represent the second most distant sampling locations (12,350 km Differences in molecular features between native and HGT-derived genes One of the main consequences of HGT is that horizontally acquired genes may have different structural characteristics when compared to native genes (cumulative effects). HGT-derived genes initially retain characteristics of the genome of the donor lineage. Consequently, the passage of time is required (and expected) to erase these differences. Therefore, we searched for differences in genomic features between HGT candidates and native Cyanidiales genes with regard to: (1) GC-content, (2) the number of spliceosomal introns and the exon/gene ratio, (3) differential transcription, (4) percent protein identity between HGT genes and their non-eukaryotic donors, and (5) cumulative effects as indicators of their non-eukaryotic origin (Danchin, 2016;Ku and Martin, 2016;Schö nknecht et al., 2013).

GC-content
All 11 Galdieria species showed significant differences (GC-content of transcripts is normally distributed, Student's t-test, two-sided, p 0.05) in percent GC-content between native sequences and HGT candidates ( Table 1). Sequences belonging to the Galdieria lineage have an exceptionally low GC-content (39%-41%) in comparison to the majority of thermophilic organisms that exhibit higher values (~55%). On average, HGT candidates in Galdieria display 1% higher GC-content in comparison to their native counterparts. No significant differences were found for C. merolae 10D and C. merolae Soos in this respect. Because native Cyanidioschyzon genes have an elevated GC-content (54%-56%), this makes it difficult to distinguish between them and HGT-derived genes (Appendix 3).

Spliceosomal introns and exon/Gene
Bacterial genes lack spliceosomal introns and therefore the spliceosomal machinery. Consequently, genes acquired through HGT are initially single-exons and may acquire introns over time due to the invasion of existing intervening sequences. We detected significant discrepancies in the ratio of single-exon to multi-exon genes between HGT candidates and native genes in the Galdieria lineage. On average, 42% of the Galdieria HGT candidates are single-exon genes, whereas only 19.2% of the native gene set are comprised of single-exons. This difference is significant (categorical data, 'native' vs 'HGT' and 'single exon' vs. 'multiple exon', Fisher's exact test, p 0.05) in all Galdieria species except G. sulphuraria SAG21.92 ( Table 1). The Cyanidioschyzon lineage contains a highly reduced spliceosomal machinery (Qiu et al., 2018), therefore only~10% of native genes are multi-exonic in C. merolae Soos and only 1/34 HGT candidates has gained an intron. C. merolae 10D has only 26 multi-exonic genes (~0.5% of all transcripts) and none of its HGT candidates has gained an intron. Enrichment testing is not possible with these small sample sizes (Appendix 4). We analyzed the number of exons that are present in multi-exonic genes and obtained similar results for the Galdieria lineage ( Table 1). All Galdieria species show significant differences regarding the exon/gene ratio between native and HGT genes (non-normal distribution regarding the number of exons per gene, Wilcoxon-Mann-Whitney-Test, 1000 bootstraps, p<=0.05). HGT candidates in Galdieria have 0.97-1.36 fewer exons per gene in comparison to their native counterparts. Because the multi-exonic HGT subset in both Cyanidioschyzon species combined includes only one multi-exonic HGT candidate, no further analysis was performed (Appendix 4).

Differential transcription
Several RNA-Seq datasets are publicly available for G. sulphuraria 074W (Rossoni, 2018) and C. merolae 10D (Rademacher et al., 2016). We aligned (Kim et al., 2015) the transcriptome reads to the respective genomes, using an identical data processing pipeline (Robinson et al., 2010) for both datasets to exclude potential algorithmic errors ( Figure 3). The average read count per gene (measured as counts per million, CPM), of native genes was 154 CPM in G. sulphuraria 074W and 196 CPM C. merolae 10D. The average read counts for HGT candidates in G. sulphuraria 074W and C. merolae 10D were 130 CPM and 184 CPM, respectively. No significant differences in RNA abundance between native genes and HGT candidates were observed for these taxa (non-normal distribution of CPM, Wilcoxon-Mann-Whitney-Test, p<0.05).
Gene function -not passage of time -explains percent protein identity (PID) between Cyanidiales HGT candidates and their non-eukaryotic donors Once acquired, any HGT-derived gene may be fixed in the genome and propagated across the lineage. The PID data can be further divided into different subsets depending on species composition of the OG. Of the total 96 OGs putatively derived from HGT events, 60 are exclusive to the Galdieria lineage (62.5%), 23 are exclusive to the Cyanidioschyzon lineage (24%), and 13 are shared by both lineages (13.5%) ( Figure 4A). Consequently, either a strong prevalence for lineage specific DL exists, or both lineages underwent individual sets of HGT events because they share their habitat with other non-eukaryotic species (which is what the HGT theory would assume). The 96 OGs in question are affected by gene loss or partial fixation. Once acquired only 8/13 of the 'Cyanidiales' (including 'Multiple HGT' and 'Uncertain') OGs and 20/60 of the Galdieria specific OGs are encoded by all species. Once acquired by the Cyanidioschyzon ancestor, the HGT candidates were retained by both C. merolae Soos and C. merolae 10D in 22/23 Cyanidioschyzon specific OGs. It is not possible to verify whether the only Cyanidioschyzon OG containing one HGT candidate is the result of gene loss, individual acquisition, or due to erroneously missing this gene model during gene prediction. The average percent PID between HGT gene candidates of the 13 OGs shared by all Cyanidiales and their non-eukaryotic donors is 41.2% (min = 24.4%; max = 65.4%) ( Figure 4B). From the HGT perspective, these OGs are derived from ancient HGT events that occurred at the root of the Cyanidiales, well before the split of the Galdieria and Cyanidioschyzon lineages. The OGs were retained over time in all Cyanidiales, although evidence of subsequent gene loss is observed. Under the DL hypothesis, this group of OGs contains genes that have been lost in all other eukaryotic lineages except the Cyanidiales. Similarly, the average PID between HGT candidates their non-eukaryotic donors in OGs Figure 3. Differential gene expression of G. sulphuraria 074W. (A) and C. merolae 10D (B), here measured as log fold change (logFC) vs transcription rate (logCPM). Differentially expressed genes are colored red (quasi-likelihood (QL) F-test, Benjamini-Hochberg, p <= 0.01). HGT candidates are shown as large circles. The blue dashes indicate the average logCPM of the dataset. Although HGT candidates are not significantly more or less expressed than native genes, they react significantly stronger to temperature changes in G. sulphuraria 074W ('more red than black dots'). This is not the case in high CO 2 treated C. merolae 10D. DOI: https://doi.org/10.7554/eLife.45017.005 exclusive to the Cyanidioschyzon lineage is 46.4% (min = 30.8%; max = 69.7%) and 45.1% (min = 27.4%; max = 69.5%) for those OGs exclusive to the Galdieria lineage. According to the HGT view, these subsets of candidates were horizontally acquired either in the Cyanidioschyzon lineage, or in the Galdieria lineage after the split between Galdieria and Cyanidioschyzon. DL would impose gene loss on all other eukaryotic lineages except Galdieria or Cyanidioschyzon. Over time, sequence similarity between the HGT candidate and the non-eukaryotic donor is expected to decrease at a rate that reflects the level of functional constraint. The average PID of 'ancient' HGT candidates shared by both lineages (before the split into Galdieria and Cyanidioschyzon approximately 800 Ma years ago ) is~5% lower than the average PID of HGT candidates exclusive to one lineage which, according to HGT would represent more recent HGT events because their acquisition occurred only after the split (thus lower divergence) ( Figure 4C). However, no significant difference between Galdieria-exclusive HGTs, Cyandioschyzon-exclusive HGTs, and HGTs shared by Only genes from the sequenced genomes were considered (13 species). A total of 60 OGs are exclusive to the Galdieria lineage (11 species), 23 OGs are exclusive to the Cyanidioschyzon lineage (two species), and 13 OGs are shared by both lineages. A total of 46/96 HGT events seem to be affected by later gene erosion/partial fixation. (B) OG-wise PID between HGT candidates vs. their potential non-eukaryotic donors. Point size represents the number of sequenced species contained in each OG. Because only two genomes of Cyanidioschyzon were sequenced, the maximum point size for this lineage is 2. The whiskers span minimum and maximum shared PID of each OG. The PID within Cyanidiales HGTs vs. PID between Cyanidiales HGTs and their potential non-eukaryotic donors is positively correlated (Kendall's tau coefficient, p=0.000747), showing evolutionary constraints that are gene function dependent, rather than time-dependent. (C) Density curve of average PID towards potential non-eukaryotic donors. The area under each curve is equal to 1. The average PID of HGT candidates found in both lineages ('ancient HGT', left dotted line) is~5% lower than the average PID of HGT candidates exclusive to Galdieria or Cyandioschyzon ('recent HGT', right dotted lines). This difference is not significant (pairwise Wilcoxon rank-sum test, Benjamini-Hochberg, p>0.05). (D) Presence/Absence pattern (green/white) of Cyanidiales species in HGT OGs. Some patterns strictly follow the branching structure of the species tree. They represent either recent HGTs that affect a monophyletic subset of the Galdieria lineage, or are the last eukaryotic remnants of an ancient gene that was eroded through differential loss. In other cases, the presence/absence pattern of Galdieria species is random and conflicts with the Galdieria lineage phylogeny. HGT would assume either multiple independent acquisitions of the same HGT candidate, or a partial fixation of the HGT candidate in the lineage, while still allowing for gene erosion. According to DL, these are the last existing paralogs of an ancient gene, whose erosion within the eukaryotic kingdom is nearly complete. DOI: https://doi.org/10.7554/eLife.45017.006 both lineages was found (non-normal distribution of percent protein identity, Shapiro-Wilk normality test, W = 0.95, p=0.002; Pairwise Wilcoxon rank-sum test, Benjamini-Hochberg, all comparisons p>0.05). Therefore, neither Cyanidioschyzon nor Galdieria specific HGTs, or HGTs shared by all Cyanidiales, are significantly more, or less, similar to their potential prokaryotic donors. We also addressed the differences in PID within the three groups. The average PID within HGT gene candidates of the 13 OGs shared by all Cyanidiales is 75.0% (min = 51.9%; max = 90.9%) ( Figure 4B). Similarly, the average PID within HGT candidates in OGs exclusive to the Cyanidioschyzon lineage is 65.1% (min = 48.9%; max = 83.8%) and 75.0% (min = 52.6%; max = 93.4%) for those OGs exclusive to the Galdieria lineage. Because we sampled only two Cyanidioschyzon species in comparison to 11 Galdieria lineages that are also much more closely related (Figure 2A), a comparison between these two groups was not done. However, a significant positive correlation (non-normal distribution of PID across all OGs, Kendall's tau coefficient, p=0.000747) exists between the PID within Cyanidiales HGTs versus PID between Cyanidiales HGTs and their non-eukaryotic donors ( Figure 4B). Hence, the more similar Cyanidiales sequences are to each other, the more similar they are to their noneukaryotic donors, showing gene function dependent evolutionary constraints.

Complex origins of HGT-impacted orthogroups
While comparing the phylogenies of HGT candidates, we also noted that not all Cyanidiales genes within one OG necessarily originate via HGT (phylogenetic trees of each HGT-OG are included in Figure 5-figure supplements 1-96). Among the 13 OGs that contain HGT candidates present in both Galdieria and Cyanidioschyzon, we found two cases ( Figure 4A, 'Multiple HGT'), OG0002305 and OG0003085, in which Galdieria and Cyanidioschyzon HGT candidates cluster in the same orthogroup. However, these have different non-eukaryotic donors and are located on distinct phylogenetic branches that do not share a monophyletic descent ( Figure 5A). This is potentially the case for OG0002483 as well, but we were uncertain due to low bootstrap values ( Figure 4A, 'Uncertain'). These OGs either represent two independent acquisitions of the same function or, according to DL, the LECA encoded three paralogs of the same gene which were propagated through evolutionary time. One of these was retained by the Galdieria lineage (and shares sequence similarity with one group of prokaryotes), the second was retained by Cyanidioschyzon (and shares sequence similarity with a different group of prokaryotes), and a third paralog was retained by all other eukaryotes. It should be noted that the 'other eukaryotes' do not always cluster in one uniformly eukaryotic clade which increases the number of required paralogs in LECA to explain the current pattern. Furthermore, some paralogs could also have already been completely eroded and do not exist in extant eukaryotes. Similarly, 6/60 Galdieria specific OGs also contain Cyanidioschyzon genes (OG0001929, OG0001938, OG0002191, OG0002574, OG0002785 and OG0003367). Here, they are nested within other eukaryote lineages and would not be derived from HGT ( Figure 5B). Also, eight of the 23 Cyanidioschyzon specific HGT OGs contain genes from Galdieria species (OG0001807, OG0001810, OG0001994, OG0002727, OG0002871, OG0003539, OG0003929 and OG0004405) which cluster within the eukaryotic branch and are not monophyletic with Cyanidioschyzon HGT candidates ( Figure 5C). According to the HGT view, this subset of candidates was horizontally acquired in either the Cyanidioschyzon lineage, or the Galdieria lineage only after the split between Galdieria and Cyanidioschyzon, possibly replacing the ancestral gene or functionally complementing a function that was lost due to genome reduction. According to DL, the LECA would have encoded two paralogs of the same gene. One was retained by all eukaryotes, red algae, and Galdieria or Cyanidioschyzon, the other exclusively by Cyanidioschyzon or Galdieria together with non-eukaryotes.
Stronger erosion of HGT genes impedes assignment to HGT or DL As already noted above, only 50/96 of the sampled HGT-impacted OGs do not appear to be affected by erosion. Dense sampling of 11 taxa within the Galdieria lineage allowed a more in-depth analysis of this issue. Here, a bimodal distribution is observed regarding the number of species per OG in the native and HGT dataset ( Figure 6C). Only 52.5% of the native gene set is present in all Galdieria strains (defined as 10 and 11 strains in order to account for potential misassemblies and missed gene models during prediction). Approximately 1/3 of the native OGs (36.1%) has been affected by gene erosion to such a degree that it is present in only one, or two Galdieria strains. In comparison, 26.7% of the candidate HGT-impacted OGs are encoded in >10 Galdieria strains, Figure 5. The analysis of OGs containing HGT candidates revealed different patterns of HGT acquisition. Some OGs contain genes that are shared by all Cyanidiales, whereas others are unique to the Galdieria or Cyanidioschyzon lineage. In some cases, HGT appears to have replaced the eukaryotic genes in one lineage, whereas the other lineage maintained the eukaryotic ortholog. Here, some examples of OG phylogenies are shown, which were simplified for ease of presentation. The first letter of the tip labels indicates the kingdom. A = Archaea (yellow), B = Bacteria (blue), E = Eukaryota (green). Branches containing Cyanidiales sequences are highlited in red. (A) Example of an ancient HGT that occurred before Galdieria and Cyanidioschyzon split into separate lineages. As such, both lineages are monophyletic (e.g., OG0001476). (B) HGT candidates are unique to the Galdieria lineage (e.g. OG0001760). (C) HGT candidates are unique to the Cyanidioschyzon lineage (e.g. OG0005738). (D) Galdieria and Cyanidioschyzon HGT candidates are derived from different HGT events and share monophyly with different non-eukaryotic organisms (e.g., OG0003085). (E) Galdieria HGT candidates cluster with non-eukaryotes, whereas the Cyanidioschyzon lineage clusters with eukaryotes (e.g., OG0001542). (F) Cyanidioschyzon HGT candidates cluster with non-eukaryotes, whereas the Galdieria lineage clusters with eukaryotes (e.g., OG0006136). DOI: https://doi.org/10.7554/eLife.45017.007 The following figure supplements are available for figure 5:     whereas 53.0% are present in less than three. The latter number might be an underestimation due to the strict threshold for HGT discovery which led to the removal of HGT candidates that were singletons. The HGT distribution is therefore skewed towards OGs containing only a few or one Galdieria species as the result of recent HGT events that occurred; for example after the split of G. sulphuraria and G. phlegrea. In spite of the strong erosion which would also lead to partial fixation of presumably recent HGT events, we analyzed whether the distribution patterns of HGT candidates across the sequenced genomes reflect the branching pattern of the species trees ( Figure 4C). This is true for all HGT candidates that are exclusive to the Cyanidioschyzon or Galdieria lineage. Either the HGT candidates were acquired after the split of the two lineages (according to HGT), or differentially lost in one of the two lineages (according to DL). In the 60 Galdieria specific OGs we found 12 OGs containing less than 10 and more than one Galdieria species ( Figure 4C). In 5/12 of the cases, the presence absence pattern reflects the species tree (OG0005087, OG0005083, GO0005479, OG0005540). Here, the potential HGT candidates are not found in any other eukaryotic species. According to HGT, they were acquired by a monophyletic sub-clade of the Galdieria lineage. According to DL, they were lost in all eukaryotes with the exception of this subset of the Galdieria lineage (e.g., OG0005280 and OG0005083 were potentially acquired or maintained exclusively by the last common ancestor of G. sulphuraria 074W, G. sulphuraria MS1, G. sulphuraria RT22, and G. sulphuraria SAG21). In the remaining OGs, the HGT gene candidate is distributed across the Galdieria lineage and conflicts with the branching pattern of the species tree. HGT would assume either multiple independent acquisitions of the same HGT candidate, or partial fixation of the HGT candidate in the lineage, while still allowing for gene erosion. According to DL, these are the last existing paralogs of an ancient gene, whose erosion within the eukaryotic kingdom is nearly complete. However, it must be considered that in some cases, DL must have occurred independently across multiple species in a brief of time after the gene was maintained for hundreds of millions of years across the lineage (e.g., OG0005224 contains G. phlegrea Soos, G. sulphuraria Azora and G. sulphuraria MS1). This implies that the gene was present in the ancestor of the Galdieria lineage and also in the last common ancestor of closely related G. sulphuraria MS1, G. sulphuraria 074W and G. sulphuraria RT22 (as well as G. sulphuraria SAG21) and the last common ancestor of closely related G. sulphuraria MtSh, G. sulphuraria Azora and G. sulphuraria YNP5587.1 (as well as G. sulphuraria 5572). A gene that was encoded and maintained since LECA, was lost independently in 6/8 species within the past few million years.

The seventy percent rule
In their analysis of eukaryotic HGT (Ku and Martin, 2016), Ku and co-authors reach the conclusion that prokaryotic homologs of genes in eukaryotic genomes that share >70% PID are not found outside individual genome assemblies (unless derived from endosymbiotic gene transfer, EGT). Hence, they are considered as assembly artifacts. We analyzed whether our dataset supports this rule, or alternatively, it is arbitrary and a byproduct of the analysis approach used, combined with low eukaryotic sampling (Richards and Monier, 2016;Leger, 2018). In addition to the 96 OGs potentially acquired through HGT, 2134 of the 9075 total OGs contained non-eukaryotic sequences, in which the Cyanidiales sequences cluster within the eukaryotic kingdom, but are similar enough to non-eukaryotic species to produce blast hits. Based on the average PID, no OG contains HGT HGT (yellow) orthogroups when compared to non-eukaryotic sequences in each OG. The red lines denote the 70% PID threshold for assembly artifacts according to 'the 70% rule'. Dots located in the top-right corner depict the 73 OGs that appear to contradict this rule, plus the 5 HGT candidates that score higher than 70%. 18/73 of those OGs are not derived from EGT or contamination within eukaryotic assemblies. (B) Density curve of average PID towards non-eukaryotic species in the same orthogroup (potential non-eukaryotic donors in case of HGT candidates). The area under each curve is equal to 1. The average PID of HGT candidates (left dotted line) is 6.1% higher than the average PID of native OGs also containing non-eukaryotic species (right dotted line candidates that share over 70% PID to their non-eukaryotic donors with OG0006191 having the highest average PID (69.68%). However, 5/96 HGT-impacted OGs contain one or more individual HGT candidates that exceed this threshold (5.2% of the HGT OGs) ( Figure 6A). These sequences are found in OG0001929 (75.56% PID, 11 Galdieria species), OG0002676 (75.76% PID, 11 Galdieria species), OG0006191 (80.00% PID, both Cyanidioschyzon species), OG0008680 (72.37% PID, 1 Galdieria species), and OG0008822 (71.17% PID, 1 Galdieria species). Moreover, we find 73 OGs with eukaryotes as sisters sharing over 70% PID to non-eukaryotic sequences (0.8% of the native OGs) ( Figure 6A). On closer inspection, the majority are derived from endosymbiotic gene transfer (EGT): 16/73 of the OGs are of proteobacterial descent and 33/73 OGs are phylogenies with gene origin in Cyanobacteria and/or Chlamydia. These annotations generally encompass mitochondrial/plastid components and reactions, as well as components of the phycobilisome, which is exclusive to Cyanobacteria, red algae, and red algal derived plastids. Of the remaining 24 OGs, 18 cannot be explained through EGT or artifacts alone unless multiple eukaryotic genomes would share the same artifact (and also assuming all gene transfers from Cyanobacteria, Chlamydia, and Proteobacteria are derived from EGT). A total of 6/24 OGs are clearly cases of contamination within the eukaryotic assemblies. Although 'the 70% rule' captures a large proportion of the dataset, increasing the sampling resolution within eukaryotes increased the number of exceptions to the rule. This number is likely to increase as more high-quality eukaryote nuclear genomes are determined. Considering the paucity of these data across the eukaryotic tree of life and the rarity of eukaryotic HGT, the systematic dismissal of eukaryotic singletons located within non-eukaryotic branches as assembly/annotation artifacts (or contamination) may come at the cost of removing true positives.

Cumulative effects
We assessed our dataset for evidence of cumulative effects within the candidate HGT-derived OGs. If cumulative effects were present, then recent HGT candidates would share higher similarity to their non-eukaryotic ancestors than genes resulting from more ancient HGT. Hence, the fewer species that are present in an OG, the higher likelihood of a recent HGT (unless the tree branching pattern contradicts this hypothesis, such as in OG 0005224, which is limited to 3 Galdieria species, but is ancient due to its presence in G. sulphuraria and G. phlegrea). In the case of DL, no cumulative effects as well as no differences between the HGT and native dataset are expected because the PID between eukaryotes and non-eukaryotes is irrelevant to this issue because all genes are native and occurred in the LECA. According to DL, the monophyletic position of Cyanidiales HGT candidates with non-eukaryotes is determined by the absence of other eukaryotic orthologs (given the limited current data) and may be the product of deep branching effects. First, we tested for general differences in PID with regard to non-eukaryotic sequences between the native and HGT datasets ( Figure 6B). Neither the PID with non-eukaryotic species in the same OG for the native dataset, nor the PID with potential non-eukaryotic donors in the same OG for the HGT dataset was normally distributed (Shapiro-Wilk normality test, p=2.2e-16/0.00765). Consequently, exploratory analysis was performed using non-parametric testing. On average, the PID with non-eukaryotic species in OGs containing HGT candidates is higher by 6.1% in comparison to OGs with eukaryotic descent. This difference is significant (Wilcoxon rank-sum test, p=0.000008).
Second, we assessed if OGs containing fewer Galdieria species would have a higher PID with their potential non-eukaryotic donors in the HGT dataset. We expected a lack of correlation with OG size in the native dataset because the presence/absence pattern of HGT candidates within the Galdieria lineage is dictated by gene erosion and thus independent of which non-eukaryotic sequences also cluster in the same phylogeny. Jonckheere's test for trends revealed a significant trend within the native subset: the fewer Galdieria species are contained in one OG, the higher the average PID with non-eukaryotic species in the same tree (Jonckheere-Terpstra, p=0.002). This was not the case in the 'HGT' subset. Here, no general trend was observed (Jonckheere-Terpstra, p=0.424).
Third, we compared the PID between HGT-impacted OGs and native OGs of the same size (OGs containing the same number of Galdieria species). This analysis revealed a significantly higher PID with non-eukaryotic sequences in favor of the HGT subset in OGs containing either one Galdieria sequence, or all 11 Galdieria sequences (Wilcoxon rank-sum test, Benjamini-Hochberg, p=2.52e-08| 3.39e-03) ( Figure 6D). Hence, the 'most recent' and 'most ancient' HGT candidates share the highest identity with their non-eukaryotic donors, which is also significantly higher when compared to native genes in OGs of the same size.

Natural habitat of extant prokaryotes with closely related orthologs
We next set out to explore the natural habitats of extant prokaryotes that harbor the closest orthologs with candidate HGTs in the Cyanidiales. To this end, we counted the frequency at which any non-eukaryotic species shared monophyly with Cyanidiales ( Table 2). A total of 568 non-eukaryotic species (19 Archaea, 549 Bacteria), from 365 different genera representing 24 divisions share monophyly with the 96 OGs containing HGT candidates. Most prominent are Proteobacteria that are sister phyla to 53/96 OGs. This group is followed by Firmicutes (28), Actinobacteria (19), Chloroflexi (12), and Bacteroidetes/Chlorobi (10). The only frequently occurring archaeal orthologs were found in Euryarchaeota (6 OGs). Interestingly, the closest orthologs often occurred in extremophilic prokaryotes that share similar (current) habitats with Cyanidiales. We hypothesize that potential non-eukaryotic HGT donors might share similar habitats because proximity is thought to favor HGT. However, we have no direct evidence of what the environment might have been at the time of HGT, or whether a third organism acted as the vector and has not been sampled in our analyses. It is worth noting that the phylogenetic data clearly demonstrate that Cyanidiales have been extremophiles for hundreds of millions of years. It is however conceivable that the HGTs may have occurred when these cells were being dispersed (they have a worldwide distribution) from one extreme site to another and would have encountered mesophilic donors at these times. Given these caveats, it is interesting to note that Sulfobacillus thermosulfidooxidans (Firmicutes), a mixotrophic, acidophilic (pH 2.0), and moderately thermophilic (45˚C) bacterium that was isolated from acid mining environments in northern Chile (where Galdieria is also present) was most prominent amongst the prokaryotic orthogroups. Sulfobacillus thermosulfidooxidans shares monophyly in 6/96 HGT-derived OGs and is followed in frequency by several species that are either thermophiles, acidophiles, or halophiles and share habitats in common with Cyanidiales ( Table 2).

Functions of horizontally acquired genes in cyanidiales
We analyzed the putative molecular functions and processes acquired through HGT. Annotations were curated using information gathered from blast, GO-terms, PFAM, KEGG, and EC.  , 2013). As such, we report a strong bias for metabolic functions among HGT candidates (Figure 7).

Metal and xenobiotic resistance/detoxification
Geothermal environments often contain high arsenic (Ar) concentrations, up to a several g/L as well as high levels of mercury (Hg), such as >200 mg/g in soils of the Norris Geyser Basin (Yellowstone National Park) and volcanic waters in southern Italy (Stauffer and Thompson, 1984;Aiuppa, 2003), both known Cyanidiales habitats (Castenholz and McDermott, 2010;Ciniglia et al., 2004;Toplin et al., 2008;Pinto, 1975). Studies with G. sulphuraria have shown an increased efficiency and speed regarding the biotransformation of HgCl 2 compared to eukaryotic algae (Kelly et al., 2007). Orthologs of OG0002305, which are present in all 13 Cyanidiales genomes, encode mercuric reductase that catalyzes the critical step in Hg 2+ detoxification, converting cytotoxic Hg 2+ into the less toxic metallic mercury, Hg 0 . Arsenate (As(V)) is imported into the cell by high-affinity P i transport systems (Meharg and MacNair, 1992;Catarecha et al., 2007), whereas aquaporins regulate arsenite (As(III)) uptake (Zhao et al., 2010). Galdieria and Cyanidioschyzon possess a eukaryotic gene-set for the chemical detoxification and extrusion of As through biotransformation and direct efflux (Schö nknecht et al., 2013). Arsenic tolerance was expanded in the Galdieria lineage through the acquisition (OG0001513) of a bacterial arsC gene, thus enabling the reduction of As(V) to As(III) using thioredoxin as the electron acceptor. It is known that As(III) can be converted into volatile dimethylarsine and trimethylarsine through a series of reactions, exported, or transported to the vacuole in conjugation with glutathione. Two separate acquisitions of a transporter annotated as ArsB are present in G. sulphuraria RT22 and G. sulphuraria 5572 (OG0006498, OG0006670), as well as a putative cytoplasmic heavy metal binding protein (OG0006191) in the Cyanidioschyzon lineage.
In the context of xenobiotic detoxification, we found an aliphatic nitrilase (OG0001760) involved in styrene degradation and three (OG0003250, OG0005087, OG0005479) Galdieria specific 4-nitrophenylphosphatases likely involved in the bioremediation of highly toxic hexachlorocyclohexane (HCH) (van Doesburg et al., 2005), or more generally other cyclohexyl compounds, such as cyclohexylamine. In this case, bioremediation can be achieved through the hydrolysis of 4-nitrophenol to 4-nitrophenyl phosphate coupled with phosophoesterase/metallophosphatase activity. The resulting cyclohexyl compounds serve as multifunctional intermediates in the biosynthesis of various heterocyclic and aromatic metabolites. A similar function in the Cyanidioschyzon lineage could be taken up by OG0006252, a cyclohexanone monooxygenase (Chen et al., 1988) oxidizing phenylacetone to benzyl acetate that can also oxidize various aromatic ketones, aliphatic ketones (e.g., dodecan-2one) and sulfides (e.g., 1-methyl-4-(methylsulfanyl)benzene). In this context, a probable multidrugresistance/quaternary ammonium compound exporter (OG0002896), which is present in all Cyanidiales, may control relevant efflux functions whereas a phosphatidylethanolamine (penicillin?) binding protein (OG0004486) could increase the stability of altered peptidoglycan cell walls. If these annotations are correct, then Galdieria is an even more promising target for industrial bioremediation applications than previously thought (Henkanatte-Gedera et al., 2017;Fukuda et al., 2018).

Cellular oxidant reduction
Increased temperature leads to a higher metabolic rate and an increase in the production of endogenous free radicals (FR), such as reactive oxygen species (ROS) and reactive nitrogen species (RNS), for example during cellular respiration (Phaniendra et al., 2015). Furthermore, heavy metals such as lead and mercury, as well as halogens (fluorine, chlorine, bromine, iodine) stimulate formation of FR (Dietz et al., 1999). FR are highly biohazard and cause damage to lipids (Ylä-Herttuala, 1999), proteins (Stadtman and Levine, 2000) and DNA (Marnett, 2000). In the case of the superoxide radical ( . O 2-), enzymes such as superoxide dismutase enhance the conversion of 2 x . O 2-, into hydrogen peroxide (H 2 O 2 ) which is in turn reduced to H 2 O through the glutathione-ascorbate cycle. Other toxic hydroperoxides (R-OOH) can be decomposed various peroxidases to H 2 O and alcohols (R-OH) at the cost of oxidizing the enzyme, which is later recycled (re-reduced) through oxidation of thioredoxin (Rouhier et al., 2008). The glutathione and thioredoxin pools and their related enzymes are thus factors contributing to a successful adaptation to geothermal environments. Here, we found a cytosolic and/or extracellular peroxiredoxin-6 (OG0005984) specific to the Cyanidioschyzon lineage and two peroxidase-related enzymes (probable alkyl hydroperoxide reductases acting on carboxymuconolactone) in the Galdieria lineage (OG0004203, OG0004392) (Chae et al., 1994). In addition, a thioredoxin oxidoreductase related to alkyl hydroperoxide reductases (OG0001486) as well as a putative glutathione-specific gamma-glutamylcyclotransferase 2 (OG0003929) are present in all Cyanidiales. The latter has been experimentally linked to the process of heavy metal detoxification in Arabidopsis thaliana (Paulose et al., 2013).

Carbon metabolism
G. sulphuraria is able to grow heterotrophically using a large variety of different carbon sources and compounds released from dying cells (Gross et al., 1995;Gross, 1998). In contrast, C. merolae is strictly photoautotrophic (De Luca et al., 1978). G. sulphuraria can be maintained on glycerol as the sole carbon source (Gross et al., 1995) making use of a family of glycerol uptake transporters likely acquired via HGT (Schö nknecht et al., 2013). We confirm the lateral acquisition of glycerol transporters in G. sulphuraria RT22 (OG0006482), G. sulphuraria Azora and G. sulphuraria SAG21 (OG0005235). The putative HGT glycerol transporters found in G. sulphuraria 074W did not meet the required threshold of two Cyanidiales sequences (from different strains) in one OG. In addition, another MIP family aquaporin, permeable to H 2 O, glycerol and other small uncharged molecules (Liu et al., 2007) is encoded by G. sulphuraria Azora (OG0007123). This could be an indication of a very diverse horizontal acquisition pattern regarding transporters. OG0003954 is the only exception to this rule, because it is present in all Galdieria lineages and is orthologous to AcpA|SatP acetate permeases involved with the uptake of acetate and succinate (Robellet et al., 2008;Sá-Pessoa et al., 2013).
The irreversible synthesis of malonyl-CoA from acetyl-CoA through acetyl-CoA carboxylase (ACCase) is the rate limiting and step in fatty acid biosynthesis. The bacterial ACCase complex consists of three separate subunits, whereas the eukaryotic ACCase is composed of a single multifunctional protein. Plants contain both ACCase isozymes. The eukaryotic enzyme is located in the cytosol and a bacterial-type enzyme consisting of four subunits is plastid localized. Three of the HGT orthogroups (OG0002051, OG0007550 and OG0007551) were annotated as bacterial biotin carboxyl carrier proteins (AbbB/BCCP), which carry biotin and carboxybiotin during the critical and highly regulated carboxylation of acetyl-CoA to form malonyl-CoA [ATP +Acetyl CoA + HCO 3-* ) ADP + Orthophosphate + Malonyl-CoA]. Whereas OG0002051 is present in all Cyanidiales and located in the cytoplasm, OG0007550 and OG0007551 are unique to C. merolae Soos and annotated as 'chloroplastic'. Prior to fatty acid (FA) beta-oxidation, FAs need to be transformed to a FA-CoA before entering cellular metabolism as an exogenous or endogenous carbon source (eicosanoid metabolism is the exception). This process is initiated by long-chain-fatty-acid-CoA ligases/acyl-CoA synthetases (ACSL) (Mashek et al., 2007) [ATP + long-chain carboxylate + CoA * ) AMP + diphosphate + Acyl-CoA]. Five general non-eukaryotic ACSL candidates were found (OG0001476, OG0002999, OG0005540, OG0008579, OG0008822). Only OG0001476 is present in all species, whereas OG0002999 is present in all Galdieria, OG0005540 in G. sulphuraria 074W and G. sulphuraria MS1, and OG0008579 and OG0008822 are unique to G. phlegrea DBV009. The GO annotation suggests moderate specificity to decanoate-CoA. However, OG0002999 also indicates involvement in the metabolism of linoleic acid, a C 18 H 32 O 2 polyunsaturated acid found in plant glycosides. ACSL enzymes share significant sequence identity but show partially overlapping substrate preferences in terms of length and saturation as well as unique transcription patterns. Furthermore, ACSL proteins play a role in channeling FA degradation to various pathways, as well as enhancing FA uptake and FA cellular retention. Although an annotation of the different ACSL to their specific functions was not possible, their involvement in the saprophytic adaptation of Cyanidioschyzon and especially Galdieria appears to be plausible.

Amino acid metabolism
Oxidation of amino acids (AA) can be used as an energy source. Once AAs are deaminated, the resulting a-ketoacids ('carbon backbone') can be used in the tricarboxylic acid cycle for energy generation, whereas the remaining NH 4 + can be used for the biosynthesis of novel AAs, nucleotides, and ammonium containing compounds, or dissipated through the urea cycle. In this context, we confirm previous observations regarding a horizontal origin of the urease accessory protein UreE (OG0003777) present in the Galdieria lineage  (the other urease genes reported in G. phlegrea DBV009 appear to be unique to this species and were thus removed from this analysis as singletons; for example ureG, OG0008984). AAs are continuously synthesized, interconverted, and degraded using a complex network of balanced enzymatic reactions (e.g., peptidases, lyases, transferases, isomerases). Plants maintain a functioning AA catabolism that is primarily used for the interconversion of metabolites because photosynthesis is the primary source of energy. The Cyanidiales, and particularly the Galdieria lineage is known for its heterotrophic lifestyle. We assigned 19/ 96 HGT-impacted OGs to this category. In this context, horizontal acquisition of protein|AA:proton symporter AA permeases (OG0001658, OG0005224, OG0005596, OG0007051) may be the first indication of adaptation to a heterotrophic lifestyle in Galdieria. Once a protein is imported, peptidases cleave single AAs by hydrolyzing the peptide bonds. Although no AA permeases were found in the Cyanidioschyzon lineage, a cytoplasmic threonine-type endopeptidase (OG0001994) and a cytosolic proline iminopeptidase involved in arginine and proline metabolism (OG0006143) were potentially acquired through HGT. At the same time, the Galdieria lineage acquired a Clp protease (OG0007596). The remaining HGT candidates are involved in various amino acid metabolic pathways. The first subset is shared by all Cyanidiales, such as a cytoplasmic imidazoleglycerol-phosphate synthase involved in the biosynthetic process of histidine (OG0002036), a phosphoribosyltransferase involved in phenylalanine/tryptophan/tyrosine biosynthesis (OG0001509) and a peptydilproline peptidyl-prolyl cis-trans isomerase acting on proline (OG0001938) (Dilworth et al., 2012). The second subset is specific to the Cyanidium lineage. It contains a glutamine/leucine/phenylalanine/valine dehydrogenase (OG0006136) (Kloosterman et al., 2006), a glutamine cyclotransferase (OG0006251) (Dahl et al., 2000), a cytidine deaminase (OG0003539) as well as an adenine deaminase (OG0005683) and a protein binding hydrolase containing a NUDIX domain (OG0005694). The third subset is specific to the Galdieria lineage and contains an ornithine deaminase, a glutaryl-CoA dehydrogenase (OG0007383) involved in the oxidation of lysine, tryptophan, and hydroxylysine (Rao et al., 2006), as well as an ornithine cyclodeaminase (OG0004258) involved in arginine and/or proline metabolism. Finally, a lysine decarboxylase (OG0007346), a bifunctional ornithine acetyltransferase/N-acetylglutamate synthase (Martin and Mulks, 1992) involved in the arginine biosynthesis (OG0008898) and an aminoacetone oxidase family FAD-binding enzyme (OG0007383), probably catalytic activity against several different L-amino acids were found as unique acquisitions in G. sulphuraria SAG21, G. phlegrea DBV009 and G. sulphuraria YNP5587.1 respectively.

One carbon metabolism and methylation
One-carbon (1C) metabolism based on folate describes a broad set of reactions involved in the activation and transfer C1 units in various processes including the synthesis of purine, thymidine, methionine, and homocysteine re-methylation. C1 units can be mobilized using tetrahydrofolate (THF) as a cofactor in enzymatic reactions, vitamin B12 (cobalamin) as a co-enzyme in methylation/rearrangement reactions and S-adenosylmethionine (SAM) (Ducker and Rabinowitz, 2017). In terms of purine biosynthesis, OG0005280 encodes an ortholog of a bacterial FAD-dependent thymidylate (dTMP) synthase converting dUMP to dTMP by oxidizing THF present in G. sulphuraria 074W, G. sulphuraria MS1, and G. sulphuraria RT22. In terms of vitamin B12 biosynthesis, an ortholog of the cobalamin biosynthesis protein CobW was found in the Cyanidioschyzon lineage (OG0002609). Much of the methionine generated through C1 metabolism is converted to SAM, the second most abundant cofactor after ATP, which is a universal donor of methyl (-CH 3 ) groups in the synthesis and modification of DNA, RNA, hormones, neurotransmitters, membrane lipids, proteins and also play a central role in epigenetics and posttranslational modifications. Within the 96 HGT-impacted dataset we found a total of 9 methyltransferases (OG0003901, OG0003905, OG0002191, OG0002431, OG0002727, OG0003907, OG0005083 and OG0005561) with diverse functions, 8 of which are SAMdependent methyltransferases. OG0002431 (Cyanidiales), OG0005561 (G. sulphuraria MS1 and G. phlegrea DBV009) and OG0005083 (G. sulphuraria SAG21) encompass rather unspecific SAMdependent methyltransferases with a broad range of possible methylation targets. OG0002727, which is exclusive to Cyanidioschyzon, and OG0002191, which is exclusive to Galdieria, both methylate rRNA. OG0002727 belongs to the Erm rRNA methyltransferase family that methylate adenine on 23S ribosomal RNA (Yu et al., 1997). Whether it confers macrolide-lincosamide-streptogramin (MLS) resistance, or shares only adenine methylating properties remains unclear. The OG0002191 is a 16S rRNA (cytidine1402-2'-O)-methyltransferase involved the modulation of translational fidelity (Kimura and Suzuki, 2010).

Osmotic resistance and salt tolerance
Cyanidiales withstand salt concentrations up to 10% NaCl (Albertano, 2000). The two main strategies to prevent the accumulation of cytotoxic salt concentrations and to withstand low water potential are the active removal of salt from the cytosol and the production of compatible solutes.
Compatible solutes are small metabolites that can accumulate to very high concentrations in the cytosol without negatively affecting vital cell functions while keeping the water potential more negative in relation to the saline environment, thereby avoiding loss of water. The G. sulphuraria lineage produces glycine/betaine as compatible solutes under salt stress in the same manner as halophilic bacteria (Imhoff and Rodriguez-Valera, 1984) through the successive methylation of glycine via sarcosine and dimethylglycine to yield betaine using S-adenosyl methionine (SAM) as a cofactor (Lu et al., 2006;Waditee et al., 2003;Nyyssola et al., 2000). This reaction is catalyzed by the enzyme sarcosine dimethylglycine methyltransferase (SDMT), which has already been characterized in Galdieria (McCoy et al., 2009). Our results corroborate the HGT origin of this gene, supporting two separate acquisitions of this function (OG0003901, OG0003905). In this context, a inositol 2dehydrogenase possibly involved in osmoprotective functions (Kingston et al., 1996) in G. phlegrea DBV009 was also found in the HGT dataset (OG0008335).

Non-Metabolic functions
Outside the context of HGT involving enzymes that perform metabolism related functions, we found 6/96 OGs that are annotated as transcription factors, ribosomal components, rRNA, or fulfilling functions not directly involved in metabolic fluxes. Specifically, two OGs associated with the bacterial 30S ribosomal subunit were found, whereas OG0002627 (Galdieria) is orthologous to the tRNA binding translation initiation factor eIF1a which binds the fMet-tRNA(fMet) start site to the ribosomal 30S subunit and defines the reading frame for mRNA translation (Simonetti et al., 2009), and OG0004339 (Galdieria) encodes the S4 structural component of the S30 subunit. Three genes functioning as regulators were found in Cyanidioschyzon, a low molecular weight phosphotyrosine protein phosphatase with an unknown regulator function (OG0002785), a SfsA nuclease (Takeda et al., 2001), similar to the sugar fermentation stimulation protein A and (OG0002871) a MRP family multidrug resistance transporter connected to parA plasmid partition protein, or generally involved in chromosome partitioning (mrp). Additionally, we found a Cyanidioschyzon-specific RuvX ortholog (OG0002578) involved in chromosomal crossovers with endonucleolytic activity (Nautiyal et al., 2016) as well as a likely Hsp20 heat shock protein ortholog (OG0004102) unique to the Galdieria lineage.

Various functions and uncertain annotations
The remaining OGs were annotated with a broad variety of functions. For example, OG0001929, OG0001810, OG0004405, and OG0001087 are possibly connected to the metabolism of cell wall precursors and components and OG0001929 (Galdieria) is an isomerizing glutamine-fructose-6-phosphate transaminase most likely involved in regulating the availability of precursors for N-and O-linked glycosylation of proteins, such as for peptidoglycan. In contrast, OG0004405 (Cyanidioschyzon) synthesizes exopolysaccharides on the plasma membrane and OG0001087 (Cyanidiales) and OG0001810 (Cyanidioschyzon) are putative undecaprenyl transferases (UPP) which function as lipid carrier for glycosyl transfer in the biosynthesis of cell wall polysaccharide components in bacteria (Apfel et al., 1999). The OGs OG0002483 and OG0001955 are involved in purine nucleobase metabolic processes, probably in cAMP biosynthesis (Galperin, 2005) and IMP biosynthesis (Schrimsher et al., 1986). A Cyanidioschyzon specific 9,15,9'-tri-cis-zeta-carotene isomerase (OG0002574) may be involved in the biosynthesis of carotene (Chen et al., 2010). Two of the 96 HGT OGs obtained the tag 'hypothetical protein' and could not be further annotated. Others had non-specific annotations, such as 'selenium binding protein' (OG0003856) or contained conflicting annotations.

Discussion
Making an argument for the importance of HGT in eukaryote (specifically, Cyanidiales) evolution, as we do here, requires that three major issues are addressed: a mechanism for foreign gene uptake and integration, the apparent absence of eukaryotic pan-genomes, and the lack of evidence for cumulative effects (Martin, 2017). The latter two arguments are dealt with below but the first concern no longer exists. For example, recent work has shown that red algae harbor naturally occurring plasmids, regions of which are integrated into the plastid DNA of a taxonomically wide array of species (Lee et al., 2016). Genetic transformation of the unicellular red alga Porphyridium purpureum has demonstrated that introduced plasmids accumulate episomally in the nucleus and are recognized and replicated by the eukaryotic DNA synthesis machinery (Li and Bock, 2018). These results suggest that a connection can be made between the observation of bacterium-derived HGTs in P. purpureum  and a putative mechanism of bacterial gene origin via longterm plasmid maintenance. Other proposed mechanisms for the uptake and integration of foreign DNA in eukaryotes are well-studied, observed in nature, and can be successfully recreated in the lab (Leger, 2018;Li and Bock, 2018).

HGT-the eukaryotic pan-genome
Eukaryotic HGT is rare and affected by gene erosion. Within the 13 analyzed genomes of the polyextremophilic Cyanidiales (Foflonker et al., 2018;Schö nknecht et al., 2014), we identified and annotated 96 OGs containing 641 single HGT candidates. Given an approximate age of 1,400 Ma years and ignoring gene erosion, on average, one HGT event occurs every 14.6 Ma years in Cyanidiales. This figure ranges from one HGT every 33.3 Ma years in Cyanidioschyzon and one HGT every 13.3 Ma in Galdieria. Still, one may ask, given that eukaryotic HGT exists, what comprises the eukaryotic pan-genome and why does it not increase in size as a function of time due to HGT accumulation? In response, it should be noted that evolution is 'blind' to the sources of genes and selection does not act upon native genes in a manner different from those derived from HGT. In our study, we report examples of genes derived from HGT that are affected by gene erosion and/or partial fixation ( Figure 4A). As such, only 8/96 of the HGT-impacted OGs (8.3%) are encoded by all 13 Cyanidiales species. Looking at the Galdieria lineage alone ( Figure 6C), 28 of the 60 lineage-specific OGs (47.5%) show clear signs of erosion (HGT orthologs are present in 10 Galdieria species), to the point where a single ortholog of an ancient HGT event may remain. When considering HGT in the Cyanidiales it is important to keep in mind the ecological boundaries of this group, the distance between habitats, the species composition of habitats, and the mobility of Cyanidiales within those borders that control HGT. Hence, we would not expect the same HGT candidates derived from the same non-eukaryotic donors to be shared between Cyanidiales and marine/freshwater red algae (unless they predate the split between Cyanidiales and other red algae), but rather between Cyanidiales and other polyextremophilic organisms. In this context, inspection of the habitats and physiology of potential HGT donors revealed that the vast majority is extremophilic and, in some cases, shares the same habitat as Cyanidiales ( Table 2). A total of 84/96 of the inherited gene functions could be connected to ecologically important traits such as heavy metal detoxification, xenobiotic detoxification, ROS scavenging, and metabolic functions related to carbon, fatty acid, and amino acid turnover. In contrast, only 6/96 OGs are related to methylation and ribosomal functions. We did not find HGTs contributing other traits such as ultrastructure, development, or behavior (Figure 7). If cultures were exposed to abiotic stress, the HGT candidates were significantly enriched within the set of differentially expressed genes (Figure 3). These results not only provide evidence of successful integration into the transcriptional circuit of the host, but also support an adaptive role of HGT as a mechanism to acquire beneficial traits. Because eukaryotic HGT is the exception rather than the rule, its number in eukaryotic genomes does not need to increase as a function of time and may have reached equilibrium in the distant past between acquisition and erosion.

HGT vs. DL
Ignoring the cumulative evidence from this and many other studies, one may still dismiss the phylogenetic inference as mere assembly artefact and overlook all the significant differences and trends between native genes and HGT candidates. This could be done by superimposing vertical inheritance (and thus eukaryotic origin) on all HGT events outside the context of pathogenicity and endosymbiosis. Under this extreme view, all extant genes would have their roots in LECA. Consequently, patchy phylogenetic distributions are the result of multiple putative ancient paralogs existing in the LECA followed by mutation, gene duplication, and gene loss. Following this line of reasoning, all HGT candidates in the Cyanidiales would be the product of DL acting on all other eukaryotic species, with the exception of the Cyanidiales, Galdieria and/or Cyanidioschyzon ( Figure 5A-C). However, we found cases where either Galdieria HGT candidates (six orthogroups), or Cyanidioschyzon HGT candidates (eight orthogroups) show non-eukaryotic origin, whereas the others cluster within the eukaryotic branch ( Figure 5E-F). In addition, we find two cases in which Galdieria and Cyanidioschyzon HGT candidates are located in different non-eukaryotic branches ( Figure 5D). DL would require LECA to have encoded three paralogs of the same gene, one of which was retained by Cyanidioschyzon, another by Galdieria, whereas the third by all other eukaryotes. The number of required paralogs in the LECA would be further increased when taking into consideration that some ancient paralogs of LECA may have been eroded in all eukaryotes and that eukaryote phylogenies are not always monophyletic which would additionally increase the number of required paralogs in the LECA in order to explain the current pattern. The strict superimposition of vertical inheritance would thus require a complex LECA, an issue known as 'the genome of Eden'.
Cumulative effects are observed when genes derived from HGT increasingly diverge as a function of time. Hence, a gradual increase in protein identity towards their non-eukaryotic donor species is expected the more recent an individual HGT event is. The absence of cumulative effects in eukaryotic HGT studies has been used as argument in favor of strict vertical inheritance followed by DL. Here, we also did not find evidence for cumulative effects in the HGT dataset. 'Recent' HGT events that are exclusive to either the Cyanidioschyzon or Galdieria lineage shared 5% higher PID with their potential non-eukaryotic donors in comparison to ancient HGT candidates that predate the split, but this difference was not significant ( Figure 4C). We also tested for cumulative effects between the number of species contained in orthogroups compared to the percent protein identity shared with potential non-eukaryotic donors under the assumption that recent HGT events would be present in fewer species in comparison to ancient HGT events that occurred at the root of Galdieria ( Figure 6D). Neither a gradual increase in protein identity for potentially recent HGT events, nor a general trend could be determined. Only orthogroups containing one Galdieria species reported a statistically significant higher protein identity to their potential non-eukaryotic donors which could be an indication of 'most recent' HGT.
Whereas the absence of cumulative effects may speak against HGT, this does not automatically argue in favor of strict vertical inheritance followed by DL. Here, the null hypothesis would be that no differences exist between HGT genes and native genes because all genes are descendants of LECA. This null hypothesis is rejected on multiple levels. At the molecular level, the HGT subset differs significantly from native genes with respect to various genomic and molecular features (e.g., GC-content, frequency of multiexonic genes, number of exons per gene, responsiveness to temperature stress) (Table 1, Figure 3). Furthermore, HGT candidates in Galdieria are significantly more similar (6.1% average PID) to their potential non-eukaryotic donors when compared to native genes and non-eukaryotic sequences in the same orthogroup ( Figure 6B). This difference cannot be explained by the absence of eukaryotic orthologs. We also find significant differences in PID with regard to non-eukaryotic sequences between HGT and native genes in orthogroups containing either one Galdieria sequence, or all eleven Galdieria sequences regarding ( Figure 6D). Hence, the 'most recent' and 'most ancient' HGT candidates share the highest resemblance to their non-eukaryotic donors, which is also significantly higher when compared to native genes in OGs of the same size. Intriguingly, a general trend towards 'cumulative effects' could be observed for native genes, highlighting the differences between these two gene sources in Cyanidiales.
Given these results and interpretations, we advocate the following view of eukaryotic HGT. Specifically, two forces may act simultaneously on HGT candidates in eukaryotes. The first is strong evolutionary pressure for adaptation of eukaryotic genetic features and compatibility with native replication and transcriptional mechanisms to ensure integration into existing metabolic circuits (e. g., codon usage, splice sites, methylation, pH differences in the cytosol). The second however is that key structural aspects of HGT-derived sequence cannot be significantly altered by the first process because they ensure function of the transferred gene (e.g., protein domain conservation, threedimensional structure, ligand interaction). Consequently, HGT candidates may suffer more markedly from gene erosion than native genes due to these countervailing forces, in spite of potentially providing beneficial adaptive traits. This view suggests that we need to think about eukaryotic HGT in fundamentally different ways than is the case for prokaryotes, necessitating a taxonomically broad genome-based approach that is slowly taking hold.
In summary, we do not discount the importance of DL in eukaryotic evolution because it can impact ca. 99% of the gene inventory in Cyanidiales. What we strongly espouse is that strict vertical inheritance in combination with DL cannot explain all the data. HGTs in Cyanidiales are significant because the 1% (values will vary across different eukaryotic lineages) helps explain the remarkable evolutionary history of these extremophiles. Lastly, we question the validity of the premise regarding the applicability of cumulative effects in the prokaryotic sense to eukaryotic HGT. The absence of cumulative effects and a eukaryotic pan-genome are neither arguments in favor of HGT, nor DL.

Materials and methods
Cyanidiales strains used for draft genomic sequencing Ten Cyanidiales strains ( Figure 1) were sequenced in 2016/2017 using the PacBio RS2 (Pacific Biosciences Inc, Menlo Park, CA) technology (Rhoads and Au, 2015) and P6-C4 chemistry (the only exception being C. merolae Soos, which was sequenced as a pilot study using P4-C2 chemistry in 2014). Seven strains, namely G. sulphuraria 5572, G. sulphuraria 002, G. sulphuraria SAG21.92, G. sulphuraria Azora, G. sulphuraria MtSh, G. sulphuraria RT22 and G. sulphuraria MS1 were sequenced at the University of Maryland Institute for Genome Sciences (Baltimore, MD). The remaining three strains, G. sulphuraria YNP5587.1, G. phlegrea Soos, and C. merolae Soos were sequenced at the Max-Planck-Institut fü r Pflanzenzü chtungsforschung (Cologne, Germany). To obtain axenic and monoclonal genetic material for sequencing, single colonies of each strain were grown at 37˚C in the dark on plates containing glucose as the sole carbon source (1% Gelrite mixed 1:1 with 2x Allen medium [Allen, 1959], 50 mM Glucose). The purity of single colonies was assessed using microscopy (Zeiss Axio Imager 2, 1000x) and molecular markers (18S, rbcL). Long-read compatible DNA was extracted using a genomic-tip 20/G column following the steps of the 'YEAST' DNA extraction protocol (QIAGEN N.V., Hilden, Germany). The size and quality of DNA were assessed via gel electrophoresis and the Nanodrop instrument (Thermo Fisher Scientific Inc, Waltham, MA).

Assembly
All genomes (excluding the already published G. sulphuraria 074W, G. phlegrea DBV009 and C. merolae 10D) were assembled using canu version 1.5 (Koren et al., 2017). The genomic sequences were polished three times using the Quiver algorithm (Chin et al., 2013). Different versions of each genome were assessed using BUSCO v.3 (Simão et al., 2015) and the best performing genome was chosen as reference for gene model prediction. Each genome was queried against the National Center for Biotechnology Information (NCBI) nr database (Geer et al., 2010) in order to detect contigs consisting exclusively of bacterial best blast hits (i.e., possible contamination). None were found.

Gene prediction
Gene and protein models for the 10 sequenced Cyanidiales were predicted using MAKER v3 beta (Cantarel et al., 2008). MAKER was trained using existing protein sequences from Cyanidioschyzon merolae 10D and Galdieria sulphuraria 074W, for which we used existing RNA-Seq (Rossoni, 2018) data with expression values > 10 FPKM (Rademacher et al., 2016) combined with protein sequences from the UniProtKB/Swiss-Prot protein database (UniProt Consortium T, 2018). Augustus (Stanke and Morgenstern, 2005), GeneMark ES (Borodovsky and Lomsadze, 2011), and EVM (Haas et al., 2008) were used for gene prediction. MAKER was run iteratively and using various options for each genome. The resulting gene models were again assessed using BUSCO v.3 (Simão et al., 2015) and PFAM 31.0 (Finn et al., 2016). The best performing set of gene models was chosen for each species.

Orthogroups and phylogenetic analysis
The 81,682 predicted protein sequences derived from the 13 genomes listed in Table 1 were clustered into orthogroups (OGs) using OrthoFinder v. 2.2 (Emms and Kelly, 2015). We queried each OG member using DIAMOND v. 0.9.22 (Buchfink et al., 2015) to an in-house database comprising NCBI RefSeq sequences with the addition of predicted algal proteomes available from the JGI Genome Portal (Nordberg et al., 2014), TBestDB (O'Brien et al., 2007), dbEST (Boguski et al., 1993), and the MMETSP (Moore Microbial Eukaryote Transcriptome Sequencing Project) (Keeling et al., 2014). The database was partitioned into four volumes: Bacteria, Metazoa, remaining taxa, and the MMETSP data. To avoid taxonomic sampling biases due to under or overabundance of particular lineages in the database, each volume was queried independently with an expect (e-value) of 1 Â 10 -5, and the top 2000 hits were saved and combined into a single list that was then sorted by descending DIAMOND bitscore. Proteins containing one or more bacterial hits (and thus possible HGT candidates) were retained for further analysis, whereas those lacking bacterial hits were removed. A taxonomically broad list of hits was selected for each query (the maximum number of genera selected for each taxonomic phylum present in the DIAMOND output was equivalent to 180 divided by the number of unique phyla), and the corresponding sequences were extracted from the database and aligned using MAFFT v7.3 (Katoh and Standley, 2013) together with queries and hits selected in the same manner for remaining proteins assigned to the same OG (duplicate hits were removed). A maximum-likelihood phylogeny was then constructed for each alignment using IQTREE v7.3 (Nguyen et al., 2015) under automated model selection, with node support calculated using 2000 ultrafast bootstraps. Single-gene trees for the referenced HGT candidates from previous research regarding G. sulphuraria 074W (Schö nknecht et al., 2013) and G. phlegrea DBV009  were constructed in the same manner, without assignment to OG. To create the algal species tree, the OG assignment was re-run with the addition of proteomes from outgroup taxa Porphyra umbilicalis (Brawley et al., 2017), Porphyridium purpureum , Ostreococcus tauri RCC4221 (Blanc-Mathieu et al., 2014), and Chlamydomonas reinhardtii (Merchant et al., 2007). Orthogroups were parsed and 2090 were selected that contained singlecopy representative proteins from at least 12/17 taxa; those taxa with multi-copy representatives were removed entirely from the OG. The proteins for each OG were extracted and aligned with MAFFT, and IQTREE was used to construct a single maximum-likelihood phylogeny via a partitioned analysis in which each OG alignment represented one partition with unlinked models of protein evolution chosen by IQTREE. Consensus tree branch support was determined by 2,000 UF bootstraps.

Detection of HGTs
All phylogenies containing bacterial sequences were inspected manually. Only trees in which there were at least two different Cyanidiales sequences and at least three different non-eukaryotic donors were retained. The singleton HGT candidates in Cyanidiales are presented in the appendix (Appendix 5) and were not analyzed further here. Phylogenies with cyanobacteria and Chlamydiae as sisters were considered as EGT and excluded from the analysis. Genes that were potentially transferred from cyanobacteria were only accepted as HGT candidates when homologs were absent in other photosynthetic eukaryotes; that is the cyanobacterium was not the closest neighbor, and when the annotation did not include a photosynthetic function, to discriminate from EGT. Furthermore, phylogenies containing inconsistencies within the distribution patterns of species, especially at the root, or UF values below 70% spanning over multiple nodes, were excluded. Each orthogroup was queried against NCBI nr to detect eukaryotic homologs not present in our databases. The conservative approach to HGT assignment used here allowed identification of robust candidates for in-depth analysis. This may however have come at the cost of underestimating HGT at the single species level. Furthermore, some of the phylogenies that were rejected because <3 non-eukaryotic donors were found may have resulted from current incomplete sampling of prokaryotes. For example, OG0001817 is present in the sister species G. sulphuraria 074W and G. sulphuraria MS1 but has a single bacterial hit (Acidobacteriaceae bacterium URHE0068, CBS domain-containing protein, GI:651323331).

Data deposit
The nuclear, plastid, and mitochondrial sequences of the 10 novel genomes, as well as gene models, ESTs, protein sequences, protein alignments, orthogroup and single gene trees, and gene annotations are available at http://porphyra.rutgers.edu. Raw PacBio RSII reads, and also the genomic, chloroplast and mitochondrial sequences, have been submitted to the NCBI and are retrievable via BioProject ID PRJNA512382.

Acknowledgements
This work was supported by the Deutsche Forschungsgemeinschaft (under Germany´s Excellence Strategy -EXC-2048/1 -Project ID: 390686111 and EXC 1028, and WE 2231/21-1 to APMW) and by the Heine Research Academy (AWR). We thank Dr. Luke Tallon and Dr. Bruno Hü ttel for the excellent technical assistance for PacBio sequencing. Dr. Marion Eisenhut for her fantastic Power Point templates. DCP and DB are grateful to the New Jersey Agricultural Experiment Station and the Rutgers University School of Environmental and Biological Sciences Genome Cooperative for supporting this research. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Data availability
The genomic, chloroplast and mitochondrial sequences of the 10 novel genomes, as well as gene models, ESTs, protein sequences, and gene annotations are available at http://porphyra.rutgers. edu. These data have also be uploaded to Dryad doi:10.5061/dryad.m06n200. Raw PacBio RSII reads, and also the genomic, chloroplast and mitochondrial sequences, have been submitted to the NCBI and are retrievable via BioProject ID PRJNA512382.
The following datasets were generated: from high evolutionary rates achieved at the subtelomeric regions. But the high evolutionary rates also made it impossible to correctly embed the aforementioned gene families in a phylogenetic tree. As such, it is not to be excluded that a similar case occurred regarding Galdieria sulphuraria's 'archaeal ATPases', although a permissive search might indicate an archaeal origin of single protein domains. Also, as only a patchy subset of the ATPases reacts to temperature fluctuations, it cannot be determined that temperature is the driving factor.
on the average value. (Right) Ranking all transcripts based upon their %GC content. Red '*' demarks HGT candidates. As the %GC content was normally distributed, students test was applied for the determination of significant differences between the native gene and the HGT candidate subset. number of exons. Red '*' demarks HGT candidates. As the number of exons was not normally distributed, transcripts were ranked by number of exons. In order to resolve the high number of tied ranks (e.g. many transcripts have two exons) a bootstrap was implied by which the rank of transcripts sharing the same number of exons was randomly assigned 1000 times. Wilcoxon-Mann-Whitney-Test applied for the determination of significant rank differences between the native gene and the HGT candidate subset. (Right) Violin plot showing the number of exons per transcript distribution across native transcripts and HGT candidates. Red '*' demarks HGT candidates. As the number of exons was not normally distributed, transcripts were ranked by number of exons. In order to resolve the high number of tied ranks (e.g. many transcripts have two exons) a bootstrap was implied by which the rank of transcripts sharing the same number of exons was randomly assigned 1000 times. Wilcoxon-Mann-Whitney-Test applied for the determination of significant rank differences between the native gene and the HGT candidate subset. (Right) Violin plot showing the number of exons per transcript distribution across native transcripts and HGT candidates.