A starting culture of G. sulphuraria strain RT22 adapted to growth at 42°C was split into two batches, which were grown separately at 42°C (control condition) and 28°C (temperature stress) for a period spanning 8 months. Bacteria were cultured on agar plates under non-photosynthetic conditions, with glucose (50 mM) as the sole carbon source. To select for fast-growing populations, the five largest colonies of each generation were picked. The samples were propagated across generations by iteratively picking the five biggest colonies from each plate and transferring them to a new plate. The picked colonies were diluted in 1 ml Allen Medium containing 25 mmol glucose. The OD₇₅₀ of the cell suspensions was measured at each re-plating step using a spectrophotometer. Approximately 1,000 cells were streaked on new plates to start the new generation. The remaining cell material was stored at -80°C until DNA extraction. This process was reiterated whenever new colonies with a diameter of 3–5 mm became visible. During the 240 days of this experiment, a total of 181 generations of G. sulphuraria RT22 were obtained for the culture grown at 42°C, whereas 102 generations were obtained for G. sulphuraria RT22 grown at 28°C.

DNA from each sample was extracted using the Genomic-tip 20/G column (QIAGEN, Hilden, Germany), following the steps of the yeast DNA extraction protocol provided by the manufacturer. DNA size and quality were assessed via gel electrophoresis and Nanodrop spectrophotometry (Thermo Fisher Scientific, Waltham, MA, United States). TruSeq DNA PCR-Free libraries (insert size = 350 bp) were generated. The samples were quantified using the KAPA library quantification kit, quality controlled using a 2100 Bioanalyzer (Agilent, Santa Clara, CA, United States), and sequenced on an Illumina (San Diego, CA, United States) HiSeq 3000 in paired-end mode (1 × 150 bp) at the Genomics and Transcriptomics Laboratory of the Biologisch-Medizinisches Forschungszentrum in Düsseldorf, Germany. The raw sequence reads are retrievable from the NCBI’s Small Read Archive (SRA) database (Project ID: PRJNA513153).

Single nucleotide polymorphisms (SNPs) and insertions/deletions (InDels) were called separately on the dataset using the GATK software version 3.6-0-g89b7209 (McKenna et al., 2010). The analysis was performed according to GATKs best practices protocols (DePristo et al., 2011; Van der Auwera et al., 2013). The untrimmed raw DNA-Seq reads of each sample were mapped onto the genome of G. sulphuraria RT22 (NCBI, SAMN10666930) using the BWA aligner (Li and Durbin, 2009) with the –M option activated to mark shorter split hits as secondary. Duplicates were marked using Picard tools¹. A set of known variants was bootstrapped for G. sulphuraria RT22 to build the covariation model and estimate empirical base qualities (base quality score recalibration). The bootstrapping process was iterated three times until convergence was reached (no substantial changes in the effect of recalibration between iterations were observed, indicating that the produced set of known sites adequately masked the true variation in the data). Finally, the recalibration model was built upon the final samples to capture the maximum number of variable sites. Variants were called using the haplotype caller in discovery mode with -ploidy set to 1 (Galdieria is haploid) and –mbq set to 20 (minimal required Phred score) and annotated using snpEff v4.3i (Cingolani et al., 2012). The called variants were filtered separately for SNPs and InDels using the parameters recommended by GATK (SNPs: “QD < 1.0 || FS > 30.0 || MQ < 45.0 || SOR > 9 || MQRankSum < -4.0 || ReadPosRankSum < -10.0,” InDels: “QD < 1.0 || FS > 200.0 || MQ < 45.0 || MQRankSum < -6.5 || ReadPosRankSum < -10.0”).

A main goal of this analysis was implementation of a method that enabled discrimination between random variants and variants that may be connected to temperature stress (non-random variants). The following logic was implemented: All variants were transformed to binary code with regards to their haplotype toward the reference genome. When the haplotype was identical to the reference genome, “0” was assigned. Variant haplotypes were assigned “1.” Random variants were gained and lost without respect to the sampling succession along the timeline and the different temperature conditions. Consequently, a “fuzzy” pattern of, e.g., “110011| 0000,” would indicate a mutation between T₀ and T₁ in the samples taken at 28°C, which was lost in T₃ and regained after T₅. The binary sequence represents the ten samples, six “cold” and four “warm,” according to their condition (“cold | warm”) and time point of sampling (“28°C_1, 28°C _2, 28°C _3, 28°C _4, 28°C _5, 28°C _6 | 42°C_1, 42°C _3, 42°C _6, and 42°C _9”). Hence, the first six digits denote samples taken at 28°C, the latter four digits those taken at 42°C; “000000| 0101” would represent a mutation in the T₂ sample taken at 42°C that was lost in T₃ and regained in T₄, and “011010| 0101” would represent a variant that does neither with respect to the sampling succession (repeated gain and loss) nor the growth condition of the samples (mutation occurs at both temperatures). By contrast, variants that were gained and fixed in the subsequent samples of a certain growth condition were considered as “non-random variants” that may reflect significant evolutionary patterns. Thus, “111111| 0000” would indicate that a mutation between T₀ and T₁ in the samples taken at 28°C was fixed over the measured period. Similarly, “000000| 0111” would indicate a mutation between T₂ and T₃ in the samples taken at 42°C that was fixed throughout the generations. As such, it was possible to determine all possible pattern combinations for non-random evolutionary patterns. The binary sequence “111111| 1111” represented the case where all ten samples contained a different haplotype when compared to the reference genome. In this specific case, systematic discrepancies between the reference genome and the DNA-Seq reads are the cause of this pattern. Variants following the “111111| 1111” pattern were removed from the dataset.

The DNA sequencing results are described in Supplementary Table S1. The Illumina HiSeq3000 raw reads reported in this project have been submitted to the NCBI’s Sequence Read Archive (SRA) and are retrievable (FASTQ file format) via BioProject PRJNA513153 and BioSamples SAMN10697271 - SAMN10697280.

Various statistical methods were applied for the different analyses performed in this project. Culture growth was measured for at 28°C (n = 6) and at 42°C (n = 10). Both datasets failed the Shapiro–Wilk normality test (p > 0.05) and showed a visible trend over time. The difference in growth between the populations was tested using the Wilcoxon rank sum test. Further, timepoints along a timeline constitute a dependent sampling approach by which the growth performance of an earlier timepoint is likely to influence the growth performance of a later timepoint.

Trends in growth over the period of this experiment were tested for significance using Jonckheere–Terpstra’s test for trends. Enrichment of GO categories as well as k-mer enrichment was tested using Fisher’s exact test for categorical data, corrected for multiple testing according to Benjamini–Hochberg. The contingency table was set up in such way that the number of times a specific GO was affected by variants was compared with the number of times the same GO was not affected by variants. This category was compared against the “background” consisting of all other GOs affected by variants and all other unaffected GOs. The same methodology was applied for k-mer enrichment testing. Differential gene expression based on previously collected data (Rossoni et al., 2018) was calculated with EdgeR (Robinson et al., 2010) implementing the QLF-test in order to address the dispersion uncertainty for each gene (Lun et al., 2016). All samples taken at 28°C were compared against all samples taken at 42°C/46°C.

Samples grown at 42°C were re-plated 10 times during the 7 months of the experiment due to faster colony growth, whereas cultures growing at 28°C were re-plated only six times (Figure 1A). Cultures grown at 42°C achieved an average doubling time of 1.32 days, equivalent to an average growth rate of 0.81/day. Cultures grown at 28°C had an average doubling time of 2.70 days, equivalent to an average growth rate of 0.39/day. This difference in doubling time/growth rate between 28°C and 42°C was significant (non-normal distribution of growth rates, Wilcoxon rank sum test, p = 0.0002) (Figure 1B). The growth rates reported here were slightly lower than in liquid batch cultures, where growth rates of 0.9/day–1.1/day were measured for heterotrophic cultures grown at 42°C (unpublished data). The changes in growth rate over time were also compared using linear regression (Figure 1C). Although the linear regression appears to indicate increasing doubling times in samples grown at 42°C, Jonckheere’s test for trends revealed no significant trend in this dataset (Jonckheere-Terpstra, p > 0.05). By contrast, samples grown at 28°C gradually adapted to the colder environment and significantly (Jonckheere–Terpstra, p < 0.05) decreased their doubling time by ∼30% during the measured period.

A total of 470,680,304 paired-end DNA-Seq reads were generated on an Illumina HiSeq 3000 sequencer. Of these, 462,869,014 were aligned to the genome (98.30%) using BWA (Supplementary Table S1). The average concordant alignment rate was 99.71%. The average genome coverage was 444 × (min = 294× and max = 579×). At least 95.5% of the sequence was covered with a depth of >20×. GATK’s haplotype caller algorithm reported 6,360 raw SNPs and 5,600 raw InDels. The SNPs and InDels were filtered separately according to GATK’s best practice recommendations. A total of 1,864 SNPs and 2,032 InDels passed the filtering process. On average, one SNP occurs every 16,177 nt and one InDel every 44,394 nt. Overall, 66.17% of the filtered variants (2578/3896) were classified as background mutations being at variance with the genome reference (“111111| 1111”). The 1243 remaining variants (966 SNPs + 277 InDels) were sorted according to their evolutionary patterns, here called “Random,” “Hot,” and “Cold” (Figure 2); 486/1243 (36.5%) are located in the intergenic region and the other 757/1243 (63.5%) in the genic region, including 5′UTR, 3′UTR, and introns. In addition, 1202/1243 (96.7%) variants followed random gain and loss patterns that do not exhibit relevant evolutionary trajectories (Figure 2). The remaining 41 variants were gained and fixed over time, thus representing non-random, evolutionary relevant variants. Twenty-three variants were fixed at 28°C and 18 variants at 42°C. Consequently, 23 variants (1.9%) were attributed to the “Cold” pattern (11 intergenic, 12 genic) and 18 variants to the “Hot” pattern (1.4%).

The vast majority of the 757 genic variants was not fixed over time and did not follow consistent evolutionary patterns (Figure 2). However, the frequency at which genes and gene functions were affected by mutation can serve as an indicator of the physiological processes most affected by evolutionary pressure at 28°C and 42°C. Here, we analyzed the functional annotations of the 757 variants located on 429 genes (6.1% of the transcriptome) using GO-Term (GO) enrichment analysis. A total of 1602 unfiltered GOs were found within the genes affected by variants (27.3% of all GOs in G. sulphuraria RT22), of which 1116 were found at least twice in the variant dataset. Of those, 234 of the GOs were significantly enriched (categorical data, “native” vs. “HGT,” Fisher’s exact test, Benjamini–Hochberg, p ≤ 0.05).

To contextualize the function in broader categories, we manually sorted all significantly enriched GOs into the following ten categories: “Cell Cycle,” “Cytoskeleton,” “Gene Regulation,” “Membrane,” “Metabolism,” “Photosynthesis,” “Stress/Signaling,” “Transport,” “Other,” and “NA” (Figure 3A). GO terms belonging to the “NA” category were considered meaningless and excluded from the analysis [e.g., “cell part” (GO:0044464), “biological process” (GO:0008150), “binding” (GO:0005488), and “ligase activity” (GO:0016874)].

To further investigate the temperature adaptation of G. sulphuraria RT22, we selected candidate genes for closer analysis using two different approaches. First, we assumed that high mutation rates in a specific gene reflect increased selective force upon its function and regulation. To identify potential targets of temperature-dependent microevolution, we searched for “variant hotspots,” here defined as the 99th percentile of genes most affected by variants. We computed variant number-dependent Z-scores for each gene and extracted genes with a Z-Score > 2.575. This procedure led to identification of seven genes, so-called “variant hotspots,” containing at least seven independent variants per gene. Next, we extracted 41 variants that followed non-random evolutionary patterns, here defined as the gain of a variant and its fixation in the subsequent samples that was exclusive to either 28°C or 42°C (Figure 2). Eighteen variants followed an evolutionary pattern defined as “Hot” (1.36%) and 23 variants followed an evolutionary pattern defined as “Cold” (1.59%). These non-random evolutionary patterns describe the gain of a variant and its fixation over time either in the 42°C dataset (“Hot,” e.g., 000000| 0001, 000000| 0011), or in the 28°C dataset (“Cold,” e.g., 000001| 0000, 000011| 0000), respectively. The underlying assumption was that this subset represented beneficial mutations. Synonymous variants were removed from further analysis. As a result, we obtained 13 genes that followed non-random evolutionary patterns (16 non-synonymous variants). An individual functional characterization of each gene is contained in the Supplementary Material (Supplementary Listing S1A for “Variant Hotspots” and Supplementary Listing S1B for “Further Non-random Genic Variants”).

The gene function of the selected temperature-dependent gene candidates broadly replicated the results of the GO enrichment analysis. Here, we found multiple enzymes involved in cell cycle control and signaling, e.g., an oxidase of biogenic tyramine (Gsulp_RT22_67_G1995), an armadillo/beta-catenin repeat family protein (Gsulp_RT22_107_G5273), the GTPase-activating ADP-ribosylation factor ArfGAP2/3 (Gsulp_RT22_82_G3036), and a peptidylprolyl cis-trans isomerase (Gsulp_RT22_64_G1844). Other candidate genes were involved in transcription and translation, e.g., a NAB3/HDMI transcription termination factor (Gsulp_RT22_83_G3136), or in ribosomal biogenesis (Gsulp_RT22_112_G5896, 50S ribosomal subunit) and required cochaperones (Gsulp_RT22_99_G4499, Hsp40). Three candidate genes were solute transporters (Gsulp_RT22_67_G2013, Gsulp_RT22_118_G6841, Gsulp_RT22_67_G1991). Most interestingly, two genes connected to temperature stress were also affected by mutations. An error-prone iota DNA-directed DNA polymerase (01_Gsulp_RT22_79_G2795), which promotes adaptive point mutation as part of the coordinated cellular response to environmental stress, was affected at 28°C (Napolitano et al., 2000; McKenzie et al., 2001), as well as the 2-phosphoglycerate kinase, which catalyzes the first metabolic step of the compatible solute cyclic 2,3-diphosphoglycerate, which increases the optimal growth temperature of hyperthermophile methanogens (Santos and da Costa, 2002; Roberts, 2005).

Horizontal gene transfer has facilitated the niche adaptation of Galdieria and other microorganisms by providing adaptive advantages (Schonknecht et al., 2013; Schönknecht et al., 2014; Foflonker et al., 2018). Five of the total 54 HGT gene candidates in G. sulphuraria RT22 gained variants (Rossoni et al., 2019). We tested whether a more significant proportion of HGT candidates gained variants in comparison to native genes. This was not the case: HGT candidates did not significantly differ from native genes (categorical data, Fisher’s exact test, p < 0.05). Of the HGT candidates, only Gsulp_RT22_67_G2013, a bacterial/archaeal APC family amino acid permease potentially involved in the saprophytic lifestyle of G. sulphuraria, accumulated a significant number of mutations (12 variants).

We tested if the 6.1% of genes that gained variants were also differentially expressed during a temperature-sensitive RNA-Seq experiment in G. sulphuraria 074W, where gene expression was measured at 28°C and 42°C (Rossoni et al., 2018). Of the 6982 sequences encoded by G. sulphuraria RT22, 4569 were successfully matched to an ortholog in G. sulphuraria 074W (65.4%); 342 were orthologs to a variant-containing gene, representing 79.7% of all genes containing variants in G. sulphuraria RT22. The dataset is representative (Wilcoxon rank sum test, Benjamini–Hochberg, p < 0.05, no differences in the distribution of variants per gene due to the sampling size). Based on this result, 36.3% of the variant-containing genes were differentially expressed. By contrast, 40.1% of the genes unaffected by variants were differentially expressed (Figure 3B). The difference between the two subsets was not significant (categorical data, Fisher’s exact test, p < 0.05). Hence, genes affected by variance during this microevolution experiment did not react more, or less, pronouncedly to fluctuating temperature.

Mutations that affect gene expression strength and pattern are a common target of evolutionary change (Barbosa-Morais et al., 2012). Intergenic DNA encodes cis-regulatory elements, such as promoters and enhancers, that constitute the binding sites of transcription factors and, thus, affect activation and transcriptional rate of genes. Promoters are required for transcriptional initiation but their presence alone results in minimal levels of downstream sequence transcription. Enhancers, which can be located either upstream, downstream, or distant from the genes they regulate, are the main drivers of gene transcription intensity and are often thought to be the critical factors of cis-regulatory divergence (Wray, 2007). Further, epigenetic changes can lead to heritable phenotypic and physiological changes without the alteration of the DNA sequence (Dupont et al., 2009). As a consequence of its evolutionary history, the genome of G. sulphuraria is highly deprived of non-functional DNA (Qiu et al., 2015). Here, we performed variant enrichment analysis of the intergenic space based on k-mers ranging from k-mer length 1 (4 possible combinations, A| C| G| T) to k-mer length 10 (1,048,576 possible combinations) to identify intergenic sequence patterns prone to variant accumulation (Table 1). The enriched k-mers were screened and annotated against the PlantCARE (Lescot et al., 2002) database containing annotations of plant cis-acting regulatory elements. Only partial hits were found, possibly due to the large evolutionary distance between plants and red algae, more specifically the Galdieria lineage, which might explain the divergence between cis-regulatory sequences (Wittkopp and Kalay, 2012). The sequence “CG,” which is the common denominator of CpG islands (Deaton and Bird, 2011), was found enriched within the k-mer set of length 2. In addition, partial hits to the PlantCARE database with a k-mer length >5 were considered as potential hits. Using this threshold, we found three annotated binding motifs, the OCT (octamer-binding motif) (Zhao, 2013), RE1 (Repressor Element 1) (Paonessa et al., 2016), and 3-AF1 (accessory factor binding sites) (Scott et al., 1996; Rhen and Cidlowski, 2005).

In this study, we subjected two populations of G. sulphuraria RT22 to a temperature-dependent microevolution experiment for 7 months. One culture was grown at 28°C, representing cold stress, and a control culture was grown at 42°C. This experiment aimed to uncover the genetic acclimation response to persistent stress, rather than the short-term acclimation response of G. sulphuraria to cold stress (Rossoni et al., 2018). We performed genomic re-sequencing along the timeline to measure changes in the genome sequence of G. sulphuraria RT22. After 7 months, corresponding to ∼170 generations of growth at 42°C and ∼100 generations of growth at 28°C, the cold-adapted cultures decreased their doubling time by ∼30%. The control cultures maintained constant growth, although a trend to slower growth might occur (Figure 1). A similar increase in the growth rate was also observed in the photoautotrophic sister lineage of Cyanidioschyzon, where cultures of Cyanidioschyzon merolae 10D were grown at 25°C for a period of ∼100 days, albeit under photoautotrophic conditions. This study found that the cold-adapted cultures outgrew the control culture at the end of the experiment (Nikolova et al., 2017). While faster doubling times at 28°C can be attributed to gradual adaptation to the suboptimal growth temperatures, we may only speculate about the causes leading to slower growth in the control condition (42°C). Perhaps G. sulphuraria RT22, which originated from the Rio Tinto river near Berrocal (Spain), may be able to thrive at high temperatures, but not for such a prolonged period.

We identified 1243 filtered variants (966 SNPs + 277 InDels), of which 757 (63.5%) were located on the coding sequence of 429 genes and 486 (36.5%) in the intergenic region. The mutation rate was estimated to be 2.17 × 10^-6/base/generation for samples grown at 28°C and 1.10 × 10^-6/base/generation for samples grown at 42°C, which we interpret as an indication of greater evolutionary stress at 28°C. Hence, suboptimal growth temperatures constitute a significant stress condition and promote the accumulation of mutations. In comparison, mutation rates in other microevolution experiments were 1.53 × 10^-8/base/generation–6.67 × 10^-11/base/generation for the unicellular green freshwater alga Chlamydomonas reinhardtii and 5.9 × 10^-9/base/generation in the green plant Arabidopsis thaliana (Ness et al., 2012; Sung et al., 2012; Perrineau et al., 2014). The 100-fold higher evolutionary rates in comparison to C. reinhardtii might result from the selective strategy employed in this experiment (only the five biggest colonies were selected to start the next generations). Although the cold-stressed samples accumulated twice as many mutations per generation in comparison to the control condition, the number of gained variants over the same period was higher in the 42°C cultures due to faster growth rates.

The impact of temperature-driven microevolution on the cellular functions of G. sulphuraria RT22 was analyzed using GO enrichment analysis. More than 75% of the 234 significantly enriched GOs affected genes functions involved in the processes of cell division, cell structure, gene regulation, and signaling. In short, the cellular life cycle appears to be targeted by variation at any stage starting with mitosis, morphogenesis, and finishing with programmed cell death. By contrast, genes directly affecting metabolic processes were less affected by mutation and made up only 10% of the enriched GOs. These observations were also confirmed through the functional annotation of the seven genes most affected by variants (“variant hotspots”) as well as the 13 genes carrying non-synonymous variants with non-random evolutionary patterns.

Historically, intergenic DNA has frequently been considered to represent non-functional DNA. It is now generally accepted that mutations affecting intergenic space can heavily influence the expression intensity and expression patterns of genes. Variants altering the sequence of cis-regulatory elements are a common source of evolutionary change (Wittkopp and Kalay, 2012). Due to two phases of genome reduction (Qiu et al., 2015), the genome of Galdieria is highly streamlined and the intergenic space accounts for only 36% of its sequence. As a consequence, it is assumed that G. sulphuraria lost non-functional intergenic regions that are affected by high random mutation rates in other organisms. In this experiment, variants accumulated proportionally between the genic and the intergenic space, which we interpret as an indication of high relevance of the non-coding regions in Galdieria. K-mer analysis revealed significant enrichment of variants occurring in CpG islands. CpG islands heavily influence transcription on the epigenetic level through methylation of the cytosines. In mammals, up to 80% of the cytosines in CpG islands can be methylated, and heavily influence epigenetic gene expression regulation. Furthermore, they represent the most common promotor type in the human genome, affecting transcription of almost all housekeeping genes and the portions of developmental regulator genes (Jabbari and Bernardi, 2004; Saxonov et al., 2006; Zhu et al., 2008). Hence, temperature adaptation is not only modulated through accumulation of mutations in the genetic region but equally driven by the alteration of gene expression through epigenetics and mutations affecting the non-coding region.