Systems Biology of Cold Adaptation in the Polyextremophilic Red Alga Galdieria sulphuraria

Abstract

Rapid fluctuation of environmental conditions can impose severe stress upon living organisms. Surviving such episodes of stress requires a rapid acclimation response, e.g., by transcriptional and post-transcriptional mechanisms. Persistent change of the environmental context, however, requires longer-term adaptation at the genetic level. Fast-growing unicellular aquatic eukaryotes enable analysis of adaptive responses at the genetic level in a laboratory setting. In this study, we applied continuous cold stress (28°C) to the thermoacidophile red alga G. sulphuraria, which is 14°C below its optimal growth temperature of 42°C. Cold stress was applied for more than 100 generations to identify components that are critical for conferring thermal adaptation. After cold exposure for more than 100 generations, the cold-adapted samples grew ∼30% faster than the starting population. Whole-genome sequencing revealed 757 variants located on 429 genes (6.1% of the transcriptome) encoding molecular functions involved in cell cycle regulation, gene regulation, signaling, morphogenesis, microtubule nucleation, and transmembrane transport. CpG islands located in the intergenic region accumulated a significant number of variants, which is likely a sign of epigenetic remodeling. We present 20 candidate genes and three putative cis-regulatory elements with various functions most affected by temperature. Our work shows that natural selection toward temperature tolerance is a complex systems biology problem that involves gradual reprogramming of an intricate gene network and deeply nested regulators.

Introduction

Small changes in average global temperature significantly affect the species composition of ecosystems. Indeed, 252 Ma years ago up to ∼95% of marine species and ∼70% of terrestrial vertebrates ceased to exist (Benton, 2008; Sahney and Benton, 2008). This event, known as the Permian–Triassic extinction, was triggered by a sharp increase in worldwide temperature (+8°C) and CO₂ concentrations (+2000 ppm) during a period spanning 48,000–60,000 years (McElwain and Punyasena, 2007; Shen et al., 2011; Burgess et al., 2014). In comparison, atmospheric CO₂ has increased by ∼100 ppm and the global mean surface temperature by ∼1°C since the sinking of the Titanic in 1912, a little more than 100 years ago. Anthropogenic climate change and its consequences have become a major evolutionary selective force (Palumbi, 2001). Higher temperatures and CO₂ concentrations result in increased seawater acidity, increased UV radiation, and changes in oceanwide water circulation and upwelling patterns. These rapid changes represent dramatically accelerating shifts in the demography and number of species, leading to loss of habitats and biodiversity (Hendry and Kinnison, 1999; Stockwell et al., 2003). A global wave of mass extinction appears inevitable (Kolbert, 2014). In this context, it is relevant to assess the effects of temperature change on genome evolution. Aquatic unicellular eukaryotes are particularly well-suited to addressing this question due to their short generation time and straightforward temperature control of their growth environment.

Microorganisms rapidly acclimate and subsequently adapt to environmental change (López-Rodas et al., 2009; Huertas et al., 2011; Romero-Lopez et al., 2012; Osundeko et al., 2014; Foflonker et al., 2018). These adaptations are driven by natural selection and involve quantitative changes in allele frequencies and phenotype within a short period of time, a phenomenon known as microevolution. The Galdieria lineage 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, and geothermal plants) where they represent up to 90% of the total biomass, competing with specialized Bacteria and Archaea (Seckbach, 1972; Castenholz and McDermott, 2010). Accordingly, members of the Galdieria lineage can cope with extremely low pH values, temperatures above 50°C, and high salt and toxic heavy metal ion concentrations (Doemel and Brock, 1971; Castenholz and McDermott, 2010; Reeb and Bhattacharya, 2010; Hsieh et al., 2018). Some members of this lineage also occur in more temperate environments (Gross et al., 2002; Ciniglia et al., 2004; Qiu et al., 2013; Barcytė et al., 2018; Iovinella et al., 2018).

Our work systematically analyzed the impact of prolonged exposure to suboptimal (28°C) and optimal (42°C) growth temperatures on the systems biology of Galdieria sulphuraria for a period spanning more than 100 generations. We chose G. sulphuraria as the model organism for this experiment due to its highly streamlined haploid genome (14 Mb, 6800 genes) that evolved out of two phases of strong selection for genome miniaturization (Qiu et al., 2015). In this genomic context, we expected maximal physiological effects of novel mutations, thus possibly reducing the fraction of random neutral mutations. Furthermore, we expected a smaller degree of phenotypical plasticity and hence a more rapid manifestation of adaptation at the genetic level.

Materials and Methods

Experimental Design and Sampling

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 Extraction and Sequencing

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).

Read Mapping and Variant Calling

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”).

Evolutionary Pattern Analysis

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.

Data Accession

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.

Statistical Analysis

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.

Results

Culture Growth

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.

FIGURE 1

Variant Calling

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%).

FIGURE 2

GO Enrichment-Based Overview of Cellular Functions Most Affected by Variants

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)].

FIGURE 3

Cell Cycle

Functions related to the cell cycle accounted for 79/234 (33.8%) of the enriched GOs. Mitosis was affected at every stage: initiation (e.g., “positive regulation of cell proliferation,” GO:0008284, p = 0.0024; “re-entry into mitotic cell cycle,” GO:0000320, p = 0.0405), DNA replication (e.g., “DNA replication, removal of RNA primer,” GO:0043137, p = 0.0001; “ATP-dependent 5′-3′ DNA helicase activity,” GO:0043141, p = 0.014), prophase (“preprophase band,” GO:0009574, p = 0.0270), metaphase (e.g., “attachment of mitotic spindle microtubules to kinetochore,” GO:0051315, p = 0.0142), anaphase (e.g., “mitotic chromosome movement toward spindle pole,” GO:0007079, p = 0.0134), and telophase (“midbody,” GO:0030496, p = 0.0404). Mutations also accumulated in genes controlling cell cycle checkpoints of mitosis (e.g., “positive regulation of mitotic metaphase/anaphase transition, GO:0045842, p = 0.0270; “mitotic spindle assembly checkpoint,” GO:0007094, p = 0.0441).

Genes with functions involved in cell differentiation and maturation of Galdieria were also affected significantly by microevolution during organellogenesis (e.g., “regulation of auxin mediated signaling pathway,” GO:0010928, p = 0.0012; “phragmoplast,” GO:0009524, p = 0.0405; “xylem and phloem pattern formation,” GO:0010051, 0.0012), cell polarity (e.g., “establishment or maintenance of epithelial cell apical/basal polarity,” GO:0045197, p = 0.0012; “growth cone,” GO:0030426, p = 0.0096), and subcellular compartmentalization and localization (“Golgi ribbon formation,” GO:0090161, p = 0.0404; ”establishment of protein localization,” GO:0045184, p = 0.0124). Interestingly, some transcriptional regulators of cell growth seem to be conserved across the eukaryotic kingdom. GOs such as “branching involved in open tracheal system development” (GO:0060446, p = 0.0012) and “eye photoreceptor cell development” (GO:0042462, p = 0.0093) were also found, indicating high amino acid sequence similarity within this category. Further, temperature stress altered genes with functions involved in cell death (e.g., “cell fate determination,” GO:0001709, p = 0.0012; “Wnt signalosome,” GO:1990909, p = 0.0405).

Gene Regulation

Maintenance of steady and balanced reaction rates across cellular systems is essential for cell survival and poses a major challenge when an organism is confronted with changes in temperature (D’Amico et al., 2002). In this context, the second largest category within the enriched GO terms (49/234, 20.9%) was related to gene regulation. Besides cell cycle control, thermal adaptation and evolution was orchestrated predominantly through mutations in genes involved in controlling the expression profiles of other genes (“gene expression,” GO:0010467, p = 0.0118). Also, we found a significant proportion of mutations affecting genes linked to the epigenetic control of gene expression, which can occur through methylation of DNA (“hypermethylation of CpG island,” GO:0044027, p = 0.0086), as well as modulation of chromatin density and histone interactions that change the accessibility of whole genomic regions to transcription (“H4 histone acetyltransferase activity,” GO:0010485, p = 0.0040) (Jenuwein and Allis, 2001; Bird, 2002). Further, variants may have altered RNA polymerase efficiency (e.g., “RNA polymerase II transcription factor binding,” GO:0001085, p = 0.0020), mRNA processing (e.g., “regulation of RNA splicing,” GO:0043484, p = 0.0025), post-transcriptional silencing (e.g., “RNA interference,” GO:0016246, p = 0.0093) as well as alteration of ribosome structure components (e.g., “structural constituent of ribosome,” GO:0003735, p = 0.0336) and rRNA methylation components (e.g., “rRNA methylation,” GO:0031167, p = 0.0036). In this regard, GO terms linked to posttranslational protein modification were also enriched (“positive regulation of peptidyl-threonine phosphorylation,” GO:0010800, 0.0086; “N-terminal peptidyl-methionine acetylation,” GO:0017196, p = 0.0007).

Cytoskeleton

Microtubule are long polymers of tubulin that are constituents of the cytoskeleton of every eukaryote. They play a central role in intracellular organization, stability, transport, organelle trafficking, and cell division (Brouhard and Rice, 2018). Because they associate spontaneously, microtubular assembly (e.g., “microtubule nucleation,” GO:0007020, p = 0.0039) and disassembly are mostly driven by tubulin concentrations at the beginning and the end of microtubules once a critical microtubule size is reached (Voter and Erickson, 1984). However, the first steps of microtubule assembly are kinetically unfavorable. Cells solve this issue by using γ-tubulin ring complex as a template (e.g., “tubulin complex,” GO:0045298, p = 0.0031). The reaction equilibrium between tubulin polymerization and monomerization is temperature-dependent and requires accurate regulation (e.g., “tau-protein kinase activity,” GO:0050321, 0.0086). Shifting temperatures from 37°C to 25°C leads to massive microtubular dissociation in homoeothermic species (Himes and Detrich, 1989). Additionally, tubulin adaptations toward lower temperatures have been observed at the level of DNA sequence as well as at the epigenetic level in psychrophilic organisms (Detrich et al., 2000). Microtubule metabolism and its role in cellular physiology accounted for 13/234 (5.6%) of the enriched GOs.

Membranes and Transport

Another major component that is also influenced by temperature is cell integrity with regards to membrane fluidity (5/234 enriched GOs, 2.1%) and transport (25/234 enriched GOs, 10.7%). Cell membranes are selectively permeable and vital for compartmentation and electric potential maintenance. In this context, Galdieria is able to maintain near-neutral cytosolic pH against a 10⁶-fold H⁺ gradient across its plasma membrane (Gross, 2000). Membranes maintain a critical range of viscosity to be able to incorporate molecules and transport substrates and nutrients. The fluidity of a membrane is mainly determined by its fatty acid composition. Changes in temperature lead to changes in fatty acid composition, which in turn affect hydrophobic interactions as well as the stability and functionality of membrane proteins and proteins anchored to membranes. Here, we measured a significant enrichment in genes with functions connected to membrane lipid bilayers (e.g., “membrane,” GO:0016020, p = 0.0002; “mitochondrial inner membrane,” GO:0005743, p = 0.0023) as well as membrane-associated proteins (e.g., “integral component of membrane,” GO:0016021, p < 0.0001), transporters (e.g., “amino acid transmembrane transporter activity,” GO:0015171, p = 5.79568E-06), and transport functions (e.g., “transmembrane transport,” GO:0055085, p = 0.0001; “cation transport,” GO:0006812, p = 0.0028). Furthermore, temperature imposes significant restrictions to vesicles, which play a central role in molecule trafficking between organelles and in endocytosis. Vesicle formation in particular appears to be affected by temperature (e.g., “vesicle organization,” GO:0016050, p = 0.0012; “clathrin-coated endocytic vesicle membrane,” GO:0030669, p = 0.0025).

Stress and Signaling

Cell signaling comprises the transformation of information, such as environmental stress, to chemical signals that are propagated and amplified through the system where they contribute to the regulation of various processes (e.g., “response to stress,” GO:0006950, p = 0.0051; “hyperosmotic response,” GO:0006972, p = 0.0039; and “ER overload response,” GO:0006983, p = 0. 0040). Here, we found a total of 18/234 GOs (7.7%) derived from genes involved in cell signaling upon which temperature changes appeared to exhibit significant evolutionary pressure driving the accumulation of variants. A broad array of receptors (G-protein coupled, tyrosine kinases, and guanylate cyclases) performs signal transduction through phosphorylation of other proteins and molecules. The signal acceptors, in turn, influence second messengers and further signaling components that affect gene regulation and protein interactions. GO annotations indicate involvement of temperature in genes coding for receptors (e.g., “activation of protein kinase activity,” GO:0032147, p = 4.95227E-05; “protein serine/threonine/tyrosine kinase activity,” GO:0004712, p = 0.00045547; and “protein autophosphorylation,” GO:0046777, p = 1.53236E-06) as well as in genes coding for the signal acceptors (“stress-activated protein kinase signaling cascade,” GO:0031098, p = 6.45014E-06; “cellular response to interleukin-3,” GO:0036016, p = 5.05371E-06; and “regulation of abscisic acid-activated signaling pathway,” GO:0009787, p = 0.006712687).

Metabolism

Maintaining metabolic homeostasis is paramount for organism survival. The efficiency, speed, and equilibrium of metabolic pathways are modulated by enzymes and the specific kinetics of each reaction. Whereas microorganisms are not capable of controlling the amount of free enthalpy in their systems (chemical equilibriums are temperature-dependent, ΔG = -RT lnk), they are able to actively adjust their metabolic rates by regulating the amount of available enzyme (“Gene Expression”). Passively, mutations can alter enzyme structure, thereby adjusting the affinity of enzymes toward ligands. Variants affecting the genetic code of genes attributed to this category influence a broad variety of metabolic pathways (e.g., “cellular aromatic compound metabolic process,” GO:0006725, p = 0.0030 and “amine metabolic process,” GO:0009308, p < 0.0001) in both anabolism (e.g., “peptidoglycan biosynthetic process,” GO:0009252, p = 0.0015 and “glycerol biosynthetic process,” GO:0006114, p = 0.0086), and catabolism (e.g., “glycosaminoglycan catabolic process,” GO:0006027, p = 0.0011). In spite of pronounced changes in gene expression of metabolic enzymes during short-term cold stress in G. sulphuraria 074W (Rossoni et al., 2018) and Cyanidioschyzon merolae 10D (Nikolova et al., 2017), microevolution of genes directly involved in metabolic steps appeared to play a minor role in long-term temperature adaptation (34/234 GOs, 14.5%).

Photosynthesis

The majority of photosynthetic light reactions are catalyzed by enzymes located in the photosynthetic thylakoid membranes. Hence, photosynthesis is based upon temperature-dependent proteins located in temperature-dependent membranes (Yamori et al., 2014). Abnormal temperatures affect the electron transport chain between the various components of the photosynthetic process (Hew et al., 1969). If the electron transport chain between photosystem I (PSI) and photosystem II (PSII) is uncoupled, electrons are transferred from PSI to oxygen instead of PSII. This process is also known as PSII excitation pressure and leads to a boost of reactive oxygen species. Long-term microevolution did not appear to significantly affect the photosynthetic apparatus of G. sulphuraria RT22 (3/234, 1.3%), likely because the experiment was performed under heterotrophic conditions in continuous darkness.

Variant Hotspots and Non-random Genic Variants

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).

HGT Candidates Are Not Significantly Involved in Temperature Microevolution

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).

Genes Involved in Differential Expression Were Not Targeted by Mutation

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.

Intergenic Variant Hotspots

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).

Discussion

Growth Rates Adapt to Temperature

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 Cyanidioschyzonmerolae 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.

Cultures Grown at 28°C Accumulate Twice the Number of Mutations as Compared to Controls

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.

Cell Cycle and Transcription Factors Are the Main Drivers of Temperature Adaptation

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.

The Intergenic Space in Galdieria Is Equally Affected as Coding Regions

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.

Conclusion

We show here that the significant growth enhancement of samples grown at 28°C over more than 100 generations was driven mainly by mutations in genes involved in the cell cycle, gene regulation, and signal transfer, as well as mutations that occurred in the intergenic regions, possibly changing the epigenetic methylation pattern and altering the binding specificity to cis-regulatory elements. Our data indicate the absence of a few specific “key” temperature switches. Rather, it appears that the evolution of temperature tolerance is underpinned by a systems response which requires the gradual adaptation of an intricate gene expression network and deeply nested regulators (transcription factors, signaling cascades, and cis-regulatory elements). Our results also emphasize the difference between short-term acclimation and long-term adaptation with regard to temperature stress, highlighting the multiple facets of adaptation that can be measured using different technologies. The short-term stress response of G. sulphuraria and the long-term stress response in C. merolae were quantified using transcriptomic and proteomic approaches, respectively (Nikolova et al., 2017; Rossoni et al., 2018). At the transcriptional and translational levels, both organisms reacted toward maintaining energetic and metabolic homeostasis by increased protein concentrations, adjusting the protein folding machinery, changing degradation pathways, regulating compatible solutes, remodeling of the photosynthetic machinery, and tuning the photosynthetic capacity. SNP and InDel calling revealed underlying regulators mostly affected by variation which are potential drivers of altered transcript and protein concentrations and ultimately determine physiology and phenotype. Some issues, however, remained unresolved. Is the observed growth phenotype permanent, or is it mostly derived from epigenetic modification which could be quickly reversed? We also did not investigate the temperature-dependent differential splicing (Bhattacharya et al., 2018; Qiu et al., 2018) apparatus in Galdieria, or the impact of non-coding RNA elements, both of which may provide additional layers for adaptive evolution (van Bakel et al., 2010).

References

• 1

  Barbosa-MoraisN. L.IrimiaM.PanQ.XiongH. Y.GueroussovS.LeeL. J.et al (2012). The evolutionary landscape of alternative splicing in vertebrate species.Science3381587–1593. 10.1126/science.1230612

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 2

  BarcytėD.ElsterJ.NedbalováL. (2018). Plastid-encoded rbcL phylogeny suggests widespread distribution of Galdieria phlegrea (Cyanidiophyceae, Rhodophyta).Nord. J. Bot.36e01794. 10.1111/njb.01794

  • CrossRef
  • Google Scholar

• 3

  BentonM. J. (2008). When Life Nearly Died: The Greatest Mass Extinction of all Time.London: Thames & Hudson.

  • Google Scholar

• 4

  BhattacharyaD.QiuH.LeeJ.YoonH. S.WeberA. P. M.PriceD. C. (2018). When less is more: red algae as models for studying gene loss and genome evolution in eukaryotes.Crit. Rev. Plant Sci.3781–99. 10.1080/07352689.2018.1482364

  • CrossRef
  • Google Scholar

• 5

  BirdA. (2002). DNA methylation patterns and epigenetic memory.Genes Dev.166–21. 10.1101/gad.947102

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 6

  BrouhardG. J.RiceL. M. (2018). Microtubule dynamics: an interplay of biochemistry and mechanics.Nat. Rev. Mol. Cell. Biol.19451–463. 10.1038/s41580-018-0009-y

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 7

  BurgessS. D.BowringS.ShenS.-Z. (2014). High-precision timeline for Earth’s most severe extinction.Proc. Natl. Acad. Sci. U.S.A.1113316–3321. 10.1073/pnas.1317692111

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 8

  CastenholzR. W.McDermottT. R. (2010). “The Cyanidiales: ecology, biodiversity, and biogeography,” inRed Algae in the Genomic Age, ed.ChapmanD. J. (Berlin: Springer), 357–371. 10.1007/978-90-481-3795-4_19

  • CrossRef
  • Google Scholar

• 9

  CingolaniP.PlattsA.WangL. L.CoonM.NguyenT.WangL.et al (2012). A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3.Fly680–92. 10.4161/fly.19695

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 10

  CinigliaC.YoonH. S.PollioA.PintoG.BhattacharyaD. (2004). Hidden biodiversity of the extremophilic Cyanidiales red algae.Mol. Ecol.131827–1838. 10.1111/j.1365-294X.2004.02180.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 11

  D’AmicoS.ClaverieP.CollinsT.GeorletteD.GratiaE.HoyouxA.et al (2002). Molecular basis of cold adaptation.Philos. Trans. R. Soc. Lond. B Biol. Sci.357917–925. 10.1098/rstb.2002.1105

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 12

  DeatonA. M.BirdA. (2011). CpG islands and the regulation of transcription.Genes Dev.251010–1022. 10.1101/gad.2037511

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 13

  DePristoM. A.BanksE.PoplinR.GarimellaK. V.MaguireJ. R.HartlC.et al (2011). A framework for variation discovery and genotyping using next-generation DNA sequencing data.Nat. Genet.43491–498. 10.1038/ng.806

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 14

  DetrichH. W.IIIParkerS. K.WilliamsR. C.Jr.NogalesE.DowningK. H. (2000). Cold adaptation of microtubule assembly and dynamics. Structural interpretation of primary sequence changes present in the alpha- and beta-tubulins of Antarctic fishes.J. Biol. Chem.27537038–37047. 10.1074/jbc.M005699200

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 15

  DoemelW. N.BrockT. (1971). The physiological ecology of Cyanidium caldarium.Microbiology6717–32. 10.1099/00221287-67-1-17

  • CrossRef
  • Google Scholar

• 16

  DupontC.ArmantD. R.BrennerC. A. (2009). Epigenetics: definition, mechanisms and clinical perspective.Semin. Reprod. Med.27351–357. 10.1055/s-0029-1237423

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 17

  FoflonkerF.MollegardD.OngM.YoonH. S.BhattacharyaD. (2018). Genomic analysis of picochlorum species reveals how microalgae may adapt to variable environments.Mol. Biol. Evol.352702–2711. 10.1093/molbev/msy167

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 18

  GrossW. (2000). Ecophysiology of algae living in highly acidic environments.Hydrobiologia43331–37. 10.1023/A:1004054317446

  • CrossRef
  • Google Scholar

• 19

  GrossW.OesterheltC.TischendorfG.LedererF. (2002). Characterization of a non-thermophilic strain of the red algal genus Galdieria isolated from Soos (Czech Republic).Eur. J. Phycol.37477–483.

  • Google Scholar

• 20

  HendryA. P.KinnisonM. T. (1999). Perspective: the pace of modern life: measuring rates of contemporary microevolution.Evolution531637–1653. 10.1111/j.1558-5646.1999.tb04550.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 21

  HewC. S.KrotkovG.CanvinD. T. (1969). Effects of temperature on photosynthesis and Co2 evolution in light and darkness by green leaves.Plant Physiol.44671–677. 10.1104/pp.44.5.671

  • CrossRef
  • Google Scholar

• 22

  HimesR. H.DetrichH. W. (1989). Dynamics of Antarctic fish microtubules at low-temperatures.Biochemistry285089–5095. 10.1021/bi00438a028

  • CrossRef
  • Google Scholar

• 23

  HsiehC. J.ZhanS. H.LiaoC. P.TangS. L.WangL. C.WatanabeT.et al (2018). The effects of contemporary selection and dispersal limitation on the community assembly of acidophilic microalgae.J. Phycol.54720–733. 10.1111/jpy.12771

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 24

  HuertasI. E.RoucoM.Lopez-RodasV.CostasE. (2011). Warming will affect phytoplankton differently: evidence through a mechanistic approach.Proc. R. Soc. B Biol. Sci.2783534–3543. 10.1098/rspb.2011.0160

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 25

  IovinellaM.ErenA.PintoG.PollioA.DavisS. J.CennamoP.et al (2018). Cryptic dispersal of Cyanidiophytina (Rhodophyta) in non-acidic environments from Turkey.Extremophiles22713–723. 10.1007/s00792-018-1031-x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 26

  JabbariK.BernardiG. (2004). Cytosine methylation and CpG, TpG (CpA) and TpA frequencies.Gene333143–149. 10.1016/j.gene.2004.02.043

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 27

  JenuweinT.AllisC. D. (2001). Translating the histone code.Science2931074–1080. 10.1126/science.1063127

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 28

  KolbertE. (2014). The Sixth Extinction: An Unnatural History, 1st Edn. New York, NY: Henry Holt and Company.

  • Google Scholar

• 29

  LescotM.DehaisP.ThijsG.MarchalK.MoreauY.Van de PeerY.et al (2002). PlantCARE, a database of plant cis-acting regulatory elements and a portal to tools for in silico analysis of promoter sequences.Nucleic Acids Res.30325–327. 10.1093/nar/30.1.325

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 30

  LiH.DurbinR. (2009). Fast and accurate short read alignment with burrows-wheeler transform.Bioinformatics251754–1760. 10.1093/bioinformatics/btp324

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 31

  López-RodasV.CostasE.ManeiroE.MarváF.RoucoM.DelgadoA.et al (2009). Living in Vulcan’s forge: algal adaptation to stressful geothermal ponds on vulcano island (southern Italy) as a result of pre-selective mutations.Phycol. Res.57111–117. 10.1111/j.1440-1835.2009.00527.x

  • CrossRef
  • Google Scholar

• 32

  LunA. T. L.ChenY.SmythG. K. (2016). “It’s DE-licious: a recipe for differential expression analyses of rna-seq experiments using quasi-likelihood methods in edgeR,” inStatistical Genomics: Methods and Protocols, edsMatheE.DavisS. (New York, NY: Springer), 391–416. 10.1007/978-1-4939-3578-9_19

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 33

  McElwainJ. C.PunyasenaS. W. (2007). Mass extinction events and the plant fossil record.Trends Ecol. Evol.22548–557. 10.1016/j.tree.2007.09.003

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 34

  McKennaA.HannaM.BanksE.SivachenkoA.CibulskisK.KernytskyA.et al (2010). The genome analysis toolkit: a map reduce framework for analyzing next-generation DNA sequencing data.Genome Res.201297–1303. 10.1101/gr.107524.110

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 35

  McKenzieG. J.LeeP. L.LombardoM.-J.HastingsP.RosenbergS. M. (2001). SOS mutator DNA polymerase IV functions in adaptive mutation and not adaptive amplification.Mol. Cell7571–579. 10.1016/s1097-2765(01)00204-0

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 36

  NapolitanoR.Janel-BintzR.WagnerJ.FuchsR. (2000). All three SOS-inducible DNA polymerases (Pol II, Pol IV and Pol V) are involved in induced mutagenesis.EMBO J.196259–6265. 10.1093/emboj/19.22.6259

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 37

  NessR. W.MorganA. D.ColegraveN.KeightleyP. D. (2012). Estimate of the spontaneous mutation rate in Chlamydomonas reinhardtii.Genetics1921447–1454. 10.1534/genetics.112.145078

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 38

  NikolovaD.WeberD.ScholzM.BaldT.ScharsackJ. P.HipplerM. (2017). Temperature-induced remodeling of the photosynthetic machinery tunes photosynthesis in the Thermophilic Algae Cyanidioschyzon merolae.J. Plant Physiol.17435–46. 10.1104/pp.17.00110

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 39

  OsundekoO.DeanA. P.DaviesH.PittmanJ. K. (2014). Acclimation of microalgae to wastewater environments involves increased oxidative stress tolerance activity.Plant Cell Physiol.551848–1857. 10.1093/pcp/pcu113

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 40

  PalumbiS. R. (2001). Evolution - Humans as the world’s greatest evolutionary force.Science2931786–1790. 10.1126/science.293.5536.1786

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 41

  PaonessaF.CriscuoloS.SacchettiS.AmorosoD.ScarongellaH.Pecoraro BisogniF.et al (2016). Regulation of neural gene transcription by optogenetic inhibition of the RE1-silencing transcription factor.Proc. Natl. Acad. Sci. U.S.A.113E91–E100. 10.1073/pnas.1507355112

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 42

  PerrineauM. M.GrossJ.ZelzionE.PriceD. C.LevitanO.BoydJ.et al (2014). Using natural selection to explore the adaptive potential of Chlamydomonas reinhardtii.PLoS One9:e92533. 10.1371/journal.pone.0092533

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 43

  QiuH.PriceD. C.WeberA. P.ReebV.YangE. C.LeeJ. M.et al (2013). Adaptation through horizontal gene transfer in the cryptoendolithic red alga Galdieria phlegrea.Curr. Biol.23R865–R866. 10.1016/j.cub.2013.08.046

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 44

  QiuH.PriceD. C.YangE. C.YoonH. S.BhattacharyaD. (2015). Evidence of ancient genome reduction in red algae (Rhodophyta).J. Phycol.51624–636. 10.1111/jpy.12294

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 45

  QiuH.RossoniA. W.WeberA. P. M.YoonH. S.BhattacharyaD. (2018). Unexpected conservation of the RNA splicing apparatus in the highly streamlined genome of Galdieria sulphuraria.BMC Evol. Biol.18:41. 10.1186/s12862-018-1161-x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 46

  ReebV.BhattacharyaD. (2010). “The thermo-acidophilic cyanidiophyceae (Cyanidiales),” inRed Algae in the Genomic Age, ed.ChapmanD. J. (Berlin: Springer), 409–426. 10.1007/978-90-481-3795-4_22

  • CrossRef
  • Google Scholar

• 47

  RhenT.CidlowskiJ. A. (2005). Antiinflammatory action of glucocorticoids - New mechanisms for old drugs.N. Engl. J. Med.3531711–1723. 10.1056/NEJMra050541

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 48

  RobertsM. F. (2005). Organic compatible solutes of halotolerant and halophilic microorganisms.Saline Systems1:5.

  • Google Scholar

• 49

  RobinsonM. D.McCarthyD. J.SmythG. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data.Bioinformatics26139–140. 10.1093/bioinformatics/btp616

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 50

  Romero-LopezJ.Lopez-RodasV.CostasE. (2012). Estimating the capability of microalgae to physiological acclimatization and genetic adaptation to petroleum and diesel oil contamination.Aquat. Toxicol.12227–237. 10.1016/j.aquatox.2012.08.001

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 51

  RossoniA. W.PriceD. C.SegerM.LyskaD.LammersP.BhattacharyaD.et al (2019). The genomes of polyextremophilic Cyanidiales contain 1% horizontally transferred genes with diverse adaptive functions.bioRxiv [Preprint]. 10.1101/526111

  • CrossRef
  • Google Scholar

• 52

  RossoniA. W.SchönknechtG.LeeH. J.RuppR. L.FlachbartS.Mettler-AltmannT.et al (2018). Cold Acclimation of the Thermoacidophilic Red Alga Galdieria sulphuraria - Changes in gene expression and involvement of horizontally acquired genes.Plant Cell Physiol.60702–712. 10.1093/pcp/pcy240

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 53

  SahneyS.BentonM. J. (2008). Recovery from the most profound mass extinction of all time.Proc. R. Soc. B Biol. Sci.275759–765. 10.1098/rspb.2007.1370

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 54

  SantosH.da CostaM. S. (2002). Compatible solutes of organisms that live in hot saline environments.Environ. Microbiol.4501–509. 10.1046/j.1462-2920.2002.00335.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 55

  SaxonovS.BergP.BrutlagD. L. (2006). A genome-wide analysis of CpG dinucleotides in the human genome distinguishes two distinct classes of promoters.Proc. Natl. Acad. Sci. U.S.A.1031412–1417. 10.1073/pnas.0510310103

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 56

  SchonknechtG.ChenW. H.TernesC. M.BarbierG. G.ShresthaR. P.StankeM.et al (2013). Gene transfer from bacteria and archaea facilitated evolution of an extremophilic eukaryote.Science3391207–1210. 10.1126/science.1231707

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 57

  SchönknechtG.WeberA. P.LercherM. J. (2014). Horizontal gene acquisitions by eukaryotes as drivers of adaptive evolution.Bioessays369–20. 10.1002/bies.201300095

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 58

  ScottD. K.MitchellJ. A.GrannerD. K. (1996). The orphan receptor COUP-TF binds to a third glucocorticoid accessory factor element within the phosphoenolpyruvate carboxykinase gene promoter.J. Biol. Chem.27131909–31914. 10.1074/jbc.271.50.31909

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 59

  SeckbachJ. (1972). On the fine structure of the acidophilic hot-spring alga Cyanidium caldarium: a taxonomic approach.Microbios5133–142.

  • Pubmed Abstract
  • Google Scholar

• 60

  ShenS. Z.CrowleyJ. L.WangY.BowringS. A.ErwinD. H.SadlerP. M.et al (2011). Calibrating the end-permian mass extinction.Science3341367–1372. 10.1126/science.1213454

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 61

  StockwellC. A.HendryA. P.KinnisonM. T. (2003). Contemporary evolution meets conservation biology.Trends Ecol. Evol.1894–101. 10.1016/s0169-5347(02)00044-7

  • CrossRef
  • Google Scholar

• 62

  SungW.AckermanM. S.MillerS. F.DoakT. G.LynchM. (2012). Drift-barrier hypothesis and mutation-rate evolution.Proc. Natl. Acad. Sci. U.S.A.10918488–18492. 10.1073/pnas.1216223109

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 63

  van BakelH.NislowC.BlencoweB. J.HughesT. R. (2010). Most “dark matter” transcripts are associated with known genes.PLoS Biol.8:e1000371. 10.1371/journal.pbio.1000371

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 64

  Van der AuweraG. A.CarneiroM. O.HartlC.PoplinR.Del AngelG.Levy-MoonshineA.et al (2013). From FastQ data to high confidence variant calls: the genome analysis toolkit best practices pipeline.Curr. Protoc. Bioinformatics4311.10.1–11.10.33. 10.1002/0471250953.bi1110s43

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 65

  VoterW. A.EricksonH. P. (1984). The kinetics of microtubule assembly. Evidence for a two-stage nucleation mechanism.J. Biol. Chem.25910430–10438.

  • Pubmed Abstract
  • Google Scholar

• 66

  WittkoppP. J.KalayG. (2012). Cis-regulatory elements: molecular mechanisms and evolutionary processes underlying divergence.Nat. Rev. Genet.1359–69. 10.1038/nrg3095

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 67

  WrayG. A. (2007). The evolutionary significance of cis-regulatory mutations.Nat. Rev. Genet.8206–216. 10.1038/nrg2063

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 68

  YamoriW.HikosakaK.WayD. A. (2014). Temperature response of photosynthesis in C-3, C-4, and CAM plants: temperature acclimation and temperature adaptation.Photosynth. Res.119101–117. 10.1007/s11120-013-9874-6

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 69

  ZhaoF. Q. (2013). Octamer-binding transcription factors: genomics and functions.Front. Biosci.18:1051–1071.

  • Google Scholar

• 70

  ZhuJ.HeF.HuS.YuJ. (2008). On the nature of human housekeeping genes.Trends Genet.24481–484. 10.1016/j.tig.2008.08.004

  • Pubmed Abstract
  • CrossRef
  • Google Scholar