Long‐term exposure to elevated temperature leads to altered gene expression in a common bloom‐forming cyanobacterium

Cyanobacteria have a strong potential to compete well under elevated temperatures. Understanding how they acclimate and evolve under climatic stressors can help us accurately predict their response to forecasted future conditions. However, it is unclear whether increased temperature results in microevolution and/or changes in gene expression. This is the first study to investigate how long‐term exposure under increased temperature influences cyanobacterial genomes. Here, we cultivated three strains of Microcystis aeruginosa (M10, M11, and M12) under two temperature conditions, ambient (22°C) and high‐temperature (26°C) for 2 yr and subsequently sequenced the full genomes. The six genomes were then compared to a reference genome and analyzed for single‐nucleotide polymorphisms, from which the mutation rate was calculated to see if temperature influenced the prevalence of gene changes. Furthermore, we investigated how temperature impacted the gene expression of six genes involved in thermal tolerance and heat shock response. We found that M. aeruginosa exposure to high temperatures demonstrated a stronger expressional response with genes associated with heat shock and thermal tolerance due to exposure to elevated temperature. Although the functionality of many genes encoding for the carbon concentrating mechanisms, nutrient metabolism and secondary metabolites were unaffected, temperature could be a possible driver of genetic change due to enhanced mutation rates. Yet, differing patterns in M10 exposed to high temperatures suggests strain specifics components are also a factor. These patterns suggest changes in plasticity, which would allow for M. aeruginosa to respond rapidly to changes in temperature and to be resilient to environmental change.

Microcystis is one of the most common and widely studied bloom-forming cyanobacteria (Zhao et al. 2018), colonizing a wide variety of global freshwater environments (Paerl et al. 2016).However, these blooms are predicted to become more intense, widespread, and frequent, as global temperatures increase (Dick et al. 2021).With this genus most commonly responsible for toxic cyanobacterial blooms, their predicted increase could lead to enhanced water contamination and public health threats (Qu et al. 2018).Freshwater lake average temperatures are predicted to rise by as much as 4 C by 2100 (Maberly et al. 2020), so it is vital to understand how Microcystis are able to dominate in such conditions (Paerl and Huisman 2008), although empirical data on their evolutionary responses are scarce (Zhao et al. 2018).Therefore, there is a need to understand the physiological and adaptive capability of key phytoplankton taxa, such as Microcystis (Ji et al. 2017;Burson et al. 2018;Huisman et al. 2018;Drug a et al. 2019).The majority of studies completed so far used short term experiments, but this research overlooks the potential of cyanobacteria to acclimate to new environments over the long term (Ji et al. 2017;Drug a et al. 2022).Understanding the evolutionary (long-term) response of cyanobacteria to multiple environmental states, allows reliable projections to be made on how aquatic freshwater ecosystems will respond to future changes such as global warming (Thomas and Litchman 2016).
Hence, further research is needed into how increased temperatures can influence the genome and its potential for microevolutionary adaptation.
Microcystis is a cyanobacterial genus with abundant genetic diversity (Pérez-Carrascal et al. 2019).Previous studies have shown high similarities in their core genomes with a large accessory genome containing genes related to secondary metabolites and toxins (Harke et al. 2016).The Microcystis genome has strain specific genes consisting of the functional categories of transferases (5.9%), transposases (2.9%) and endonucleases (2.6%).This is combined with a large proportion of repeat sequences, restriction enzymes, and insertion sequences compared to other cyanobacteria (Frangeul et al. 2008;Meyer et al. 2017), which could play a role in the ability of the genus to respond to changing environments (Meyer et al. 2017).Increases in temperature affects Microcystis in a variety of ways including metabolic rate, chlorophyll a concentration, growth rate and toxin concentrations (Dick et al. 2021).
Recent genomic studies have shown that warming can lead to some species rapidly evolving through modulating phenotypic and genetic responses (McGaughran et al. 2021).Evolution involves changes within the genomes such as mutations or gene removal (Kutschera and Niklas 2004), including microevolution.This is a change in DNA sequences through mutation or nucleotide replacement, as a response to intra-specific variation, while macroevolution is a result of interspecific variation (Therkildsen et al. 2013).From a genetic perspective, temperature can directly increase mutation rates (through replication errors and DNA damage) and indels (insertions and deletions of one or more nucleotide) leading to increased de novo genetic diversity (Chu et al. 2018;McGaughran et al. 2021;Arromrak et al. 2022).
Increased plasticity is defined as an increase in the speed in which an organism can change its phenotype repetitively (such as upregulation/downregulation in the expression of certain genes) in response to variations in environmental conditions (Ratikarinen and Kokko 2019).These include alterations of gene expression such as the upregulation of certain heat shock, thermal tolerance, oxidative stress, and metabolism genes, which could indicate changes to plasticity (Tillich et al. 2014).Furthermore, temperature can influence the production and concentrations of secondary metabolites such as toxins (microcystins; Walls et al. 2018).Previous studies have shown a lack of consensus about whether elevated temperatures cause an increase or decline in toxin production (Bui et al. 2018;Walls et al. 2018;Wilhelm et al. 2020).Therefore, by exposing Microcystis aeruginosa (a commonly occurring, unicellular, colonyforming cyanobacterium; Yamaguchi et al. 2019) to increased temperatures, we can determine whether a strain of Microcystis evolves or undergoes changes to gene expression.
Variant calling (single-nucleotide polymorphisms [SNPs] and Indels) is a crucial technique for comparing genomes, as it can determine nucleotide differences between organisms (Olson et al. 2015).SNPs and indels can identify differences in nucleotides position relative to a reference genome (Nielsen et al. 2011).Previous research has shown that variant calling can help to determine differences between various strains of cyanobacteria (e.g., whether one strain of Synechococcus had substantially higher growth rates at elevated temperatures and enhanced light intensities, compared to a non-adapted strain due to advantageous mutations, identified with variant calling; Ungerer et al. 2018).However, the position of the SNPs and indels is vital, as those found around translation/coding and transcription start sites, can have a direct effect on gene transcription and protein translation (Neininger et al. 2019).Chu et al. (2018) determined the bias of SNPs toward coding vs. noncoding DNA in high temperatures, suggesting temperature can influence mutational heterogeneity across different genomic regions.This might be caused by the dependence of temperature on the function and target of the DNA repair mechanisms.This is the first study on the impact of long-term exposure to increased temperature on Microcystis genomes.
Understanding the drivers of genomic change in M. aeruginosa, especially with respect to temperature, has considerable socioeconomic and ecological implications.Although, previous studies have shown acclimation in phytoplankton can result from multiple drivers, and phenotypes can evolve following multiple paths (Brennan et al. 2017;Walworth et al. 2021), the information on how M. aeruginosa will respond from the genetic perspective to a warmer climate is missing.This research aimed to determine if exposure to increased temperature resulted in microevolution and changes to the genome, or in altered gene expression.First, we analyzed the expression of six genes responsible for heat shock and thermal tolerance, throughout the evolution experiment.Second, we compared the genomes of the ambient and high-temperature strains using SNP analysis to determine the differences within the 27 genes of interest (GOI) against a reference genome (NIES-843).These techniques were combined to ascertain if long-term exposure to increased temperature results in microevolution and/or variations in gene expression.

Strain collection and growth conditions of the Microcystis strains
The three M. aeruginosa strains (M10, M11, M12) used in this study were isolated from three lakes, Ciuperca (45.18 N,28.79 E),Tas ¸aul (44.35 N,28.61 E),and Gorgova (45.15 N,29.18 E), in the summer of 2018.These eutrophic freshwater lakes had an average summer water temperature of 22 C.After collection, the samples were deposited as non-axenic cultures in the Collection of Cyanobacteria and Algae at the Institute of Biological Research in Cluj-Napoca, Romania (Dragos ¸1997).The phylogenetic identities of the three strains were confirmed using the 16S rDNA gene amplified with specific primers (Frank et al. 2008).The polymerase chain reaction (PCR) fragments were sequenced by a third-party company (Macrogen Europe, Amsterdam, The Netherlands).A sequence similarity search against GenBank public database showed the affiliation of the three sequences to M. aeruginosa species.
All strains were exposed for 2 yr (from June 2018 to June 2020) in semi-batch conditions (corresponding to approx.200 generations) in ultrapure water-based BG11 medium (Allen and Stanier 1968; Supporting Information Table S1) under controlled 16 h : 8 h light : dark conditions provided by white LED lamps (25 μmol photon m À2 s À1 ).The strains were grown under two different temperature conditions, 22 C (mean lake summer temperature) and 26 C (the predicted temperature by 2100 [IPCC 2018]).To distinguish between the different M. aeruginosa strains and exposure conditions, the following labeling was used: strain_exposure temperature.Samples were grown in triplicates using 100-mL glass tubes.Each sample was bubbled daily with filtered atmospheric air (0.22 μm; Minisart, Sartorius, Göttingen, Germany) for 90 min (15 min aeration every 4 h) to prevent cells gathering.All samples were diluted every 3 weeks regardless of strain or temperature of exposure using ultrapure water-based BG11 medium (Supporting Information Table S1) to the same starting density of OD 600 = 0. 1 (approximately 1.06 Â 10 6 cells L À1 ).Previously collected data showed that the M. aeruginosa strains reach the stationary growth phase after approximately 3 weeks and, therefore using this dilution schedule, prevented nutrient depletion.All other growth conditions were kept identical for all cultures.

Growth rate measurement and thermal reaction norms
The full methodology for the growth rate measurements and thermal reaction norms (TRNs) for M11 are stated in Briddon et al. (2022).The same methodology was used for M10 and M12.The growth rate measurements were calculated based on OD, with all samples diluted to the same starting density of OD 600 = 0.1 (approximately 1.06 Â 10 6 cells L À1 ), before each growth rate and TRN experiment.During the 2-yr evolution experiment, the growth rate (μ) of both ambient and high-temperature M10, M11, and M12 were measured (during 2-week experiments) at the beginning of the exposure period (June 2018, Month 0) as well as after 8, 16, and 24 months of cultivation at 22 C and 26 C. All strains were tested in their original environment (22 C).This was used to determine if exposing M. aeruginosa to 26 C gave them a competitive advantage when grown back in the ambient temperature (22 C).Therefore, the growth rate was only measured at 22 C regardless of whether the strain was previously exposed at either temperature.The strains were not acclimated back to ambient conditions before measuring their growth rates at 22 C.This step of acclimation (of approx. 1 week) is sometimes used in such experiments; however, during a side experiment, we noticed that the growth parameters of the strains acclimated at 26 C did not reverse to the original values even after several months of being cultured back at 22 C (Supporting Information Fig. S1).Therefore, a 1-week re-acclimation step would not have resulted in significantly different growth rates.The growth rate was calculated using the equation, μ = (ln Ndln N0)/d, where N0 and Nd represent culture density at the beginning and the end of each test and d is the duration of the test in days.Number of divisions per day (Div d À1 ) was calculated as follows: Div d À1 = μ/ln2, while generation time (Gen.t)was calculated according to: Gen.t = 1/Div d À1 .One-way ANOVAs and twosample t-tests were used to determine if growth rates differed between strains, the temperature of exposure, and time.
TRNs were completed to determine the optimal and lethal temperatures for growth for the ambient and hightemperature M10, M11, and M12 strains, to see if 2 yr of exposure to increased temperature resulted in a shift in their lethal temperature.All six strains were grown at a range of temperatures from 20 C to 40 C, with intervals of 2 C. The growth rate was calculated as shown above.

Reverse transcription, quantitative PCR
The full methodology for reverse transcription, quantitative PCR (RT-qPCR) including primer design, RNA extraction, and cDNA synthesis are stated in Briddon et al. (2022).After each growth rate experiment, completed at Months 0, 8, 16 and 24, biomass was collected for the strains M10, M11, and M12 for the ambient and high-temperature conditions.All cultures (both ambient and high-temperature ones) were tested at ambient conditions (22 C), except for "Month 0_26" versions, where the cells were collected after 1 d of exposure to 26 C, to see if exposing the cultures to 26 C for just 1 d resulted in a change in gene expression.The biomass was subsequently harvested for RNA and analyzed for cDNA synthesis and RT-qPCR to determine the gene expression of six genes involved in heat shock response and thermal tolerance (cya1, pnp, pyrR, clpC1, clpC2, sigF;Tillich et al. 2014).The qPCR results were calculated as fold change gene expressions of the seven target genes (2 ÀΔΔCt ) relative to the M. aeruginosa M10, M11 and M12 rnpB (reference gene) and the initial ambient-temperature Month 0 condition.

DNA isolation
The genomes of the six M. aeruginosa were analyzed after 2 yr of exposure to different temperatures.First, the water from each of the three replicates was combined (totalling approx.240 mL), and then centrifuged to remove the excess liquid.Genomic DNA of the six Microcystis samples was extracted using the Quick-DNA™ Fecal/Soil Microbe Kits (Zymo Research) according to the manufacturer's instructions.The DNA concentration and quality were assessed using a Nan-oDrop™ 2000 Spectrophotometer.Subsequently, the genomic DNA was subjected to the whole-genome sequencing (WGS) procedure, carried out by a third-party company (Genewiz, Germany).An Illumina MiSeq 2 Â 150 bp sequencing configuration was used and with a quality guarantee of Q30 for greater than 75% of bases.

SNP analysis
The raw WGS data were cleaned with default settings in Fastp (Chen et al. 2018).The mapping of the short reads onto the M. aeruginosa NIES-843 reference genome was completed using the BWA-MEM (Li and Durbin 2010).The mapped reads were processed with SAMtools (Li et al. 2009), piping the output to GATK3 (UnifiedGenotyper; DePristo et al. 2011) for calling the SNPs and indels.To calculate the maximum likelihood (ML) phylogenetic trees, we pruned the SNP markers for a 100% genotyping rate among all the six samples.In addition, we corrected for linkage disequilibrium using the indeppairwise 100 5 0.5 parameter in PLINK 1.9 (Chang et al. 2015), which removes a pair of SNPs if the linkage disequilibrium is greater than 0.5, while considering a window of 100 SNPs, and shifting the window 5 SNPs forward before repeating the procedure.The resulting SNP markers were concatenated, converted to fasta format, and used as an input for MEGA X (Kumar et al. 2018), to calculate the best model for building a ML tree, which included the M. aeruginosa NIES-843 reference genome in the analysis.Model T92 + I was chosen, as it had the lowest Bayesian Information Criterion (BIC) score and the bootstrap consensus tree was inferred from 1000 replicates.Bedtools (Quinlan and Hall 2010) was used to intersect the SNP data set with the GFF file available for the reference genome.All the SNPs detected within the 27 GOI were used to calculate a 2 nd ML tree, using the T92 model.GC bias in the six samples was calculated using the CollectGcBiasMetrics implemented in Picard ("Picard Toolkit," 2019.Broad Institute, GitHub Repository; https:// broadinstitute.github.io/picard/;Broad Institute).The SNP density plot was created with Circos (Krzywinski et al. 2009).To determine if the pressure of selection at different temperatures influenced the frequency of gene changes, the mutation rate was calculated.First, we determined the number of SNPs for each M. aeruginosa genome (compared to the reference genome, NIES_843).Subsequently, the number of SNPs was divided by the approximate number of generations that would occur over the 2 yr of exposure.This was used to determine if exposure to two different temperature conditions resulted in conflicting mutation rates.

Genomes assembly and annotation
As the cultures were not axenic, we treated the sequence data as belonging to small metagenomes.Quality filtration of the reads, genomes assembly and binning, gene prediction and annotation was completed as described in Chiriac et al. (2022).Viral sequences were predicted with VirSorter and Vibrant, and contigs that had > 25% of their genes annotated as viral genes were eliminated from the bins (Roux et al. 2015;Kieft et al. 2020).The level of completeness and contamination were estimated with CheckM v1.0.18 (Parks et al. 2015) and the default collection of 120 bacterial marker genes set.Taxonomy was assigned to the bins using GTDB-Tk v.

Growth rate evolution (including TRNs)
The growth rate and TRN results for M11 ambient and hightemperature conditions were previously published in Briddon et al. (2022).At Month 0 (the start of the two-year evolution experiment), the growth rate of M11 was significantly higher than M10 (determined using a t-test; p < 0.001), but not M12 (Fig. 1a).However, this did not result in significantly higher growth rates for M11 over time for either condition (ambient or high temperature).There was no pattern for the strains M10 and M12 over the 2 yr of exposure, with similar growth rates in Month 24 to Month 0. At the end of the experiment (Month 24), only M11 demonstrated a significant difference between the ambient and high-temperature conditions (determined using a t-test; p < 0.001).The TRNs showed clear differences between the lethal temperatures of the ambient and hightemperature strains for M11 (Fig. 1b).The lethal temperature for the high-temperature M11 was 2 C higher compared to the ambient condition.However, for M10 and M12, the ambient and high-temperature conditions had similar growth rates at all temperatures including the lethal temperature.

qPCR gene expression
The full qPCR results for the strain M11 (including both the ambient and high-temperature conditions) are published in Briddon et al. (2022).Ambient-temperature M10 showed an increase in gene expression (upregulation) of fivefold to sixfold for the genes clpC1 and clpC2 (at Month 24 compared to Month 0).However, missing results (due to sample loss) for cya1, pnp, pyrR, and sigF means any trends in gene expression cannot be determined (Fig. 2a).Both M11 and M12 exposed to ambient temperatures showed a slight upregulation in clpC1 but either no trend or a downregulation in the remaining genes over the 2 yr (Fig. 2b,c).All high-temperature conditions (M10, M11, and M12) demonstrated a proportional increase in the gene expression of the genes clpC1, and clpC2.M11 and M12 exposed to high temperatures also showed an upregulation in pyrR, while M11_26 also had an increase in the gene expression of the sigF and pnp genes over time.

General features of the six M. aeruginosa genomes
The sequencing of the 6 genomes resulted in between 160 and 324 contigs.The genomes have an average GC content of 42.7% (AE 0.16%; Supporting Information Table S2) and a coding density of 84% (AE 0.4%).Genome completeness varied between 85.79% (M12_22) and 94.74% (M12_26) with contamination of between 0.07% (M10_26 and M12_22) and 0.51% (M11_26 and M12_26).Variant analysis results showed 88,704 SNPs in the 6 M. aeruginosa genomes compared to the references genome (NIES-843), 1117 insertions and 955 deletions (Supporting Information Table S3).The majority of SNP base changes were transitions (211,974) compared to transversions (106,980; Supporting Information Tables S4, S5).The SNPs were well spread out across the genomes with no areas of high variability (Fig. 3).The mutation rate (e.g., the number of SNPs that occurred per generation) was higher in the high-temperature condition (for M11 and M12) compared to the ambient condition, but this pattern was not observed in M10_26 (Supporting Information Table S6).

Carbon concentration mechanisms
In all six genomes, the genes of the two bicarbonate and two CO 2 uptake systems for inorganic carbon were present.These include the two sodium-dependent bicarbonate uptake systems (bicA and sbtA) and BCT1 (encoded by the genes cmpABCD), which is ATP dependent (Song and Qiu 2007).SbtB (encoded by a companion gene near sbtA) regulates the activity of sbtA (Du et al. 2014) and was present in all genomes.The gene bicA was truncated in all genomes.This is in combination with the NDH uptake systems that convert the passively diffusing CO 2 within the cell into HCO À 3 via NADPH-driven electron flow (Sandrini et al. 2015), again observed in all genomes.

Nutrient assimilation
Among the Microcystis genomes studied, the genes associated with the P and N metabolic pathway were observed in all samples (as stated in Harke and Gobler 2013).For P metabolism, the major genes include those associated with phosphorus scavenging, such as the phosphate-binding protein (pstS), phosphate import ATP-binding protein (pstB), and the phosphate transport system permease protein (pstC).Both members of the two-component regulatory system involved in the phosphate regulon gene expression were also present in all six genomes.This includes the phosphate regulon of transcriptional regulatory proteins (phoB1 and phoB2; a positive regulator for P assimilation for when phosphate is limited; Lu et al. 2019) and the membrane-associated protein (phoR) kinase that phosphorylates phoB1 in response to environmental signals.PhoB2 was absent in M10_22 and M12_22.
Analysis of genes that necessitate N metabolism (glnAB, glnBF, gltBDSX, nadABC, nrtB-D, ureA-F, urtA-E) demonstrated, that the majority of them were present in all six genomes including those involved in ammonium, nitrate/nitrite, and urea transport (Harke and Gobler 2013).An exception to this is that a majority of the nitrogenase structural genes were absent from all genomes (such as nigD, nifK, nigF; Pancrace et al. 2017).Only the nifB and nifH genes were observed in all genomes.This coincides with previous studies which highlighted that Microcystis has limited/no ability to support N 2 fixation (Lu et al. 2019).

Active secondary metabolites including toxins
Microcystis has the ability to produce a variety of secondary metabolites such as microcystins, microviridins, cyanopeptolins, aeruginosins, and cyanobactins, some of which are toxic (Filatova et al. 2021).The microcystin gene cluster encodes a single enzyme complex of 55 kb and comprises of 10 genes (mcyA-J).It consists of three peptide synthases (mcyA-C), a molecular polyketide synthase (mcyD), two hybrid enzymes comprising of peptide synthetase and polypektide synthase models (mcyF and mcyG) and enzymes involved in tailoring (mcyJ, mcyI, and mcyF) and transport (mcyH) of the toxin.Interestingly, the whole multifunctional gene complex required for microcystin biosynthesis was absent from the M. aeruginosa expoured to high temperatures yet present in those strains exposed to ambient temperatures.For M10_26 the coverage of the mcy complex was fragmented and had poor consensus.The organization of the mcy cluster was identical in the strains exposed to ambient temperatures.To confirm these results, a PCR was completed to confirm the presence of the whole mcy cluster in all genomes (Supporting Information Fig. S2).The PCR confirmed the presence of the whole mcy cluster in M10_26 and M11_26.Enzyme immunoassay (ELISA) completed demonstrated the production of microcystins in all samples with concentrations ranging from 0.153 μg L À1 (M12_26) to 2.525 μg L À1 (M11_22; Supporting Information Fig. S3), with no relationship between temperature of exposure and microcystin concentration (Supporting Information Table S6).
The cluster of genes involved in the biosynthesis of microviridins (mdnA-E; Uzunov et al. 2021) were present in all genomes, except for mdnA in M10_26, which was heavily fragmented and had poor consensus.This secondary metabolite is an important agent against the natural predators of cyanobacteria by acting as an antifeeding agent (Do Amaral et al. 2021).
The mcn gene cluster consists of seven genes (mcnA-G) with a total size of 32 kp (Chen et al. 2021).They encode cyanopeptolins, a potential neurotoxin and occur through nonribosomal peptide synthase (NRPS; Rounge et al. 2007).
McnA-E encode the NRPS modules, while mcnF has an ABC Fig. 2. The gene expression results for the six genes for M10, M11, and M12 exposed to both the ambient and high temperature conditions for Months 0, 8, 16, and 24.All strains were tested in their original environment (22 C) except for "Month 0_M10_26," "Month 0_M11_26," and "Month 0_M12_26," which refers to the expression pattern of the cultures maintained at 26 C for only 24 h.The values represent the relative abundance of the target genes based on rnpB, which was used as reference gene.These results were normalized against Month 0 ambient M10, M11, or M12 to determine if temperature has caused a fold change in gene expression (2 ÀΔΔCt ).Each test was done in three replicates and then averaged.The abbreviations on the y-axis states the strain and the condition that it was exposed to and the month of the evolution experiment the samples were analyzed (e.g., Month_8_M10_22 is the strain M10 exposed to the ambient condition (22 C) in the eighth month).The different shades of color correspond to the strength of the gene expression, blue for downregulation, white for neutral and red for upregulation during the 2-yr evolution experiment.The gray color corresponds to missing results.
transporter domain and mcnG has an unknown function (Huang and Zimba 2019).The presence of mcnBCE was unclear in M11_26 and M12_26; therefore, a PCR was carried out for the genes mcnBCE, which confirmed the presence of the three genes in all genomes (Supporting Information Fig. S4).
Aeruginosins are potent inhibitors of serine proteases with certain variants being highly toxic to invertebrates (Pearson et al. 2020).The gene cluster spans 25 kb and consists of > 8 genes (aerA-N; Chen et al. 2021).While aerK, aerD-F, and aerN were present in all genomes, aerL was absent from all.A PCR confirmed the presence of the genes aerAB (Supporting Information Fig. S4).
Cyanobactins, which are small cyclic peptides help defend against predators (Leikoski et al. 2012), and consist of the genes pirA-G.PirA, pirB, pirC, pirE1, pirE2, and pirG were present in all genomes, with only the gene pirF being absent.

Genes of interest
The 27 GOI involved a variety of different gene clusters and processes and were chosen due to their enhancement under elevated temperature (Meyer et al. 2017; Guan  -843).There are no areas with a large number of SNPs, which suggest that no regions of the genomes are highly variable (characterized by a cluster of dark blue bars).In order to calculate the SNP density, the 5.8 Mbp genomes was split in 10k bp bins and the number of SNPs in each such bin was calculated (as shown by each individual bar).The color and the length of the bars represent the number of SNPs found in each 10k bp bins, with the darker the blue and the longer a bar is, the higher the number of SNPs found.These were later plotted using the Circos software to produce the figure.The max length of the bars is 637, that is, the maximum number of SNPs per bin considered.et al. 2018;Willis and Woodhouse 2020;Dick et al. 2021).All three GOI involved in the carbon concentration mechanism (bicA, sbA, and sbtB) were present in all genomes; however, bicA was incomplete (Table 1).For P and N metabolism, the genes cphB (used for N storage) and phnD (encodes for the phosphonate-binding protein of the ABC-type phosphonate transporter found in picoplankton; Ilikchyan et al. 2009) were absent from all samples.Only nrtA (nitrate/ nitrite binding protein) had a deletion within the middle of the gene and lacked a stop codon for all strains.Most of the secondary metabolites (mcyDE, aerAB, pirA, mdnDE, and mcnBCE) were present in all genomes.Except for aerG1 and mdnA, which had fragmented coverage for the M. aeruginosa that had been exposed to high temperatures (Table 1).
Of the GOIs, only mcyD had an insertion within the gene, which was only observed in the ambient-temperature condition.The 15 bp sequence (5 0 -AGAAAAATCCTCTGT-3 0 ) is identical in all samples exposed to the ambient-temperature condition and contains direct repeats of the target site (eight bases; 5 0 -AGAAAAAT-3 0 ) and does not contain any start/stop codons.As the ELISA results (Supporting Information Fig. S3) showed the production of toxins in genomes that had been exposed to ambient temperatures, the insertion does not appear to impact functionality.No other indels were observed in the other GOIs.The start/stop codons were present for all the GOIs, which had complete sequences.The only exception to this is the mcnB and mcnC, where the stop codon (of mcnB) and the start codon (of mcnC) were absent as the two genes were merged in all genomes.Only one SNP occurred in the start/stop codon in mcnC, with the stop codon changing from TTA to TCA.As the change still resulted in a stop codon, it is unlikely to affect the functionality of the gene.

Phylogenetic tree
The phylogenetic tree using pruned SNPs (totaling 2315) from across the six M. aeruginosa genomes and one reference genome, showed clear groupings (Fig. 4a).The two M10 conditions (ambient and high temperature) were most closely associated with each other, while, the M11 and M12 strains are separated by temperature.The Microcystis reference genome (NIES-843) was most closely associated with M10_22, M10_26, M11_26, and M12_26.The phylogenetic tree created by SNPs within the GOI followed a similar pattern except the Microcystis reference genome (NIES-843), was most closely associated with M11_22 (Fig. 4b).

Discussion
Investigating how the long-term exposure to temperature of M. aeruginosa influences phenotypic plasticity and expression, can help explore how species might respond under a future warming climate and whether it results in reversible plasticity and/or microevolution (Sandoval-Castillo et al. 2020).Analysis of the qPCR gene expression of three strains of M. aeruginosa showed an upregulation in the gene clpC1 regardless of the temperature of exposure.While the M. aeruginosa exposed to high temperatures also demonstrated an upregulation clpC2, pyrR (except M10_26), and sigF and pnp (M11_26 only).Tillich et al. (2014) determined that the genes clpC1, clpC2, pyrR, sigF, and pnp gave elevated fitness advantage under increased temperatures and were paramount for the enhanced thermal tolerance in Synechocystis.However, the lack of consistency suggests that other factors apart from temperature are also influencing gene expression.Both Huang et al. (2002) and Marin et al. (2004) determined that the ClpC genes have been linked to various stress responses apart from temperature, such as high-light irradiance.Furthermore, temperature change can denature proteins and change their shape (Trotter et al. 2001).Clp (an ATPdependent Clp protease ATPase subunit) is responsible for digesting misfolded proteins.Therefore, its upregulation in all conditions could be the result of a higher proportion of misfolded proteins due to a heat-shock response, regardless of temperature of exposure (Tillich et al. 2014).Also, Eriksson and Clarke (1996) found that the level of gene expression of the ClpC proteins remained unchanged (regardless of temperature), when the gene ClpB was present in Synechococcus.Upregulation (as temperature increased) of the ClpC genes only occurred when ClpB was missing.As the gene ClpB is present in all genomes, this could explain the similar levels of gene expression of the ClpC genes for the M. aeruginosa in both conditions.This suggests that temperature can result in changes to gene expression; however, it is not necessarily a proportional response.Yet, we also found that exposure to increased temperature could have led to substantial changes in the genomes, with high-temperature M11 and M12 having a higher mutation rate compared to those samples exposed to ambient temperatures.
All genomes have a genotype I, in that they have the bicA and sbtAB genes, which is common, due to the genetic diversity among Microcystis.The presence of both transporters suggests these genomes are versatile in their growth responses to different inorganic carbon concentrations (Sandrini et al. 2014).CO 2 and inorganic C concentrations can be highly variable within lake systems.The supersaturation of lakes from high CO 2 (resulting from a variety of processes such as respiration, lake mixing and water temperature; Yang et al. 2022) can result in dissolved CO 2 concentrations reaching 10,000 ppm, yet enhanced photosynthetic activity can also deplete dissolved CO 2 , raising the pH to over 10.5 (Lazzarino et al. 2009).Furthermore, carbonate chemistry can vary more strongly within dense algal blooms (or cultures) due to distortion by biological processes.This can lead to substantial shifts in C availability (Wolf et al. 2022), however, without pH measurements it is not possible to quantify how the pH varied.Because of such fluctuations, the presence of both transporters could provide a competitive advantage depending on the CO 2 concentrations and other environmental conditions, including situations where CO 2 concentrations continue to increase (Dick et al. 2021).
The similarities between the genes used for N and P metabolism suggest that exposure in laboratory conditions to temperature has limited influence on nutrient assimilation.Constant dilutions and replenishment of nutrients could be a possible explanation for the loss of the genes cphB (used for N storage; Harke and Gobler 2015) and phnD.Ilikchyan et al. (2009) determined that the phnD gene was expressed during P limitation in the ARC-21 strain of freshwater Synechococcus sp.Furthermore, the Synechococcus sp.CC9902 lacked a phoB gene, so phnD could have been used as an alternative control (Ilikchyan et al. 2009).Therefore, high P concentrations and the presence of the phoB in all genomes, could possibly explain its disappearance.As the expression of different nutrient assimilation genes was not investigated as part of this study, it would be interesting to test whether exposure to enhanced temperatures results in upregulation/downregulation of these genes.
Interestingly, there was no relationship between microcystin concentration and temperature of exposure (Supporting Information Table S6).This is in agreement with other studies, which highlight conflicting results on whether enhanced temperatures leads to an increase or decline in the production of toxins (Bui et al. 2018;Martin et al. 2020;Dick et al. 2021).Walls et al. (2018) determined that the optimal temperature for microcystin production was 15-25 C, which declined in elevated temperatures in the toxin-producing Planktothrix agardhii.Moreover, Martin et al. (2020) determined that a temperature reduction caused a doubling of the microcystin quota over a 7-9 d period, yet a reversal led to the microcystin quota declining.This could possibly explain why the M11 and M12 that had been exposed to high temperatures, had a lower microcystin production compared to the ambient-temperature versions.However, high-temperature M10 had the highest microcystin production of all samples (2.525 μg L À1 ).Laboratory and field studies have shown that toxic strains can benefit more from warming than non-toxic ones, yet this differed between strains (Bui et al. 2018).This indicates the response to temperature is strain dependent, suggesting that the costbenefit analysis of microcystin production is not necessarily beneficial at higher temperatures.
Of all the strains exposed to high temperatures, M11 and M12 have a greater genomic similarity compared to M10 (as according to the ML phylogenetic trees; Fig. 4), suggesting strain specific factors, apart from temperature, could have caused the genetic differences.All long-term exposure conditions were identical to minimize the number of differences between the samples, therefore, temperature is most likely to be the driver of differences between genomes, however, other factors must be considered such as microbiome diversity, viral infections, and laboratory conditions.Viruses can influence the evolution of host populations such as M. aeruginosa and can drive intraspecific diversification (Dick et al. 2021).The high heterogeneity of CRISPR-Cas systems in M. aeruginosa demonstrates its evolutionary dynamics (Yang et al. 2015).Visual observation during the 2-yr exposure process suggests the presence of cyanophages and the VIRSORTER results highlighted multiple viral fragments, suggesting that viruses could have facilitated strain-specific differences.Microcystis can form attached microbiomes of heterotrophic bacteria (Dick et al. 2021).Recent studies have suggested that secondary metabolites can be shared between the Microcystis and its microbiome, resulting in only one having to encode certain genes (Pérez-Carrascal et al. 2021).Long-term association during long-term exposure process could result in the loss of redundant genes in M. aeruginosa strains (Black Queen Hypothesis; Morris et al. 2012).The two genomes with the most varied microbiomes (M10_26 and M12_22; Supporting Information Table S7), also had the most fragmented/ truncated gene coverage of secondary metabolites.Therefore, it could be hypothesized that the diversity of the microbiomes in these two samples could explain some of the differences in the secondary metabolite genes (Janssen 2019).Laboratory conditions could have also influenced the composition of secondary metabolites.The physiology for non-axenic cultures can change as a response to culturing conditions and time of cultivation (Meyer et al. 2017).However, as genomes of the three strains were not sequenced at the beginning of the exposure period, the impact of the cultivation period on genomic features is unknown.
The grouping of M11 and M12 (exposed to high temperatures) suggests temperature could be a possible driver of genetic microevolution (as shown by the phylogenetic trees; Fig. 4).Temperature has been identified as one of the dominant drivers modulating phenotypic and genetic responses (Arromrak et al. 2022).It can cause rapid evolutionary change through an increase in growth rates and elevated chance of mutations (McGaughran et al. 2021).Prokaryotes such as Microcystis have a higher rate of adaptation to warming compared to more complex organisms (Barton et al. 2020).Nevertheless, Qu et al. (2022) determined that during 2 yr of exposure to different temperature regimes, the marine diazotroph Trichodesmium displayed no capacity to adapt to warming.At the same time, Crocosphaera (also a marine N 2 -fixing cyanobacterium) underwent adaptation to enhanced temperature, through a suite of specific genetic changes, suggesting that evolutionary capacity could be considered for this organism.Furthermore, Jerney et al. (2019) found no evidence of long-term adaptation to warming in a marine dinoflagellate, which was able to adjust its growth to different temperatures and salinities, showing a broad tolerance and a high acclimation potential under a range of temperatures.These studies have shown that although longterm exposure to temperature can be a driver of genetic change, it is not universal across all microalgae, supporting our results that the selection pressure of high temperature results in alterations to gene expression rather than genetic changes in Microcystis sp.This is shown by the M11 and M12 exposed to high temperatures having the highest mutation rate per generation (of 201.6 and 201.5, respectively) compared to the ambient condition.Yet, the strain M10 (which had also been exposed to high temperatures) had a lower mutation rate (of 199), which could indicate a partial change that had not progressed as far as M11 and M12, suggesting a strain-specific effect.There was also evidence that temperature caused the upregulation of genes resulting in a fitness advantage.This is especially true for the strain M11 that was exposed to high temperatures, which had significantly higher growth rates and lethal temperature after 2 yr of heat exposure (Fig. 1A).This, combined with differences in growth rates between the strains exposed to high temperatures suggests that tolerance to h.zwarming under selection pressure differs within species.

Conclusion
In summary, this study shows that the upregulation of genes associated with thermal tolerance suggests that temperature can lead to plastic responses within M. aeruginosa strains exposed to high temperatures (compared to those exposed to ambient temperatures).The differences that we found between genomes (SNPs, higher mutation rate in the strains exposed to high temperatures) suggest intraspecies variation, but not evidence of phenotypic adaptation in reference to growth rate, due to the long-term exposure to increased temperature.Instead, it seems that changes to gene expression suggest that M. aeruginosa is able to respond plastically to long-term selection under high-temperature conditions.It would be interesting to test our observations using field experiments to determine if exposure to enhanced temperature still causes changes in gene expression and genotype composition in natural phytoplankton communities.

Fig. 1 .
Fig. 1. (A) Mean growth rate and standard deviation of Microcystis aeruginosa strains (M10, M11, and M12) measured every 8 months for the two temperature conditions used for the evolution experiment (ambient [22 C; mean lake summer temperature for the present day] and high temperature [26 C; the predicted temperature by 2100[) and (B) the thermal reaction norms with standard deviation showing the growth rate at 2 C intervals between 20 C and 40 C for both ambient and high-temperature M10, M11, and M12 after 2 yr of exposure.

Fig. 3 .
Fig. 3. Whole-genome comparisons of the six, 5.8 Mbp Microcystis aeruginosa genomes showing the location of the SNPs in comparison to the reference genome (NIES-843).There are no areas with a large number of SNPs, which suggest that no regions of the genomes are highly variable (characterized by a cluster of dark blue bars).In order to calculate the SNP density, the 5.8 Mbp genomes was split in 10k bp bins and the number of SNPs in each such bin was calculated (as shown by each individual bar).The color and the length of the bars represent the number of SNPs found in each 10k bp bins, with the darker the blue and the longer a bar is, the higher the number of SNPs found.These were later plotted using the Circos software to produce the figure.The max length of the bars is 637, that is, the maximum number of SNPs per bin considered.

Fig. 4 .
Fig. 4. ML phylogenetic trees from (A) using pruned 2315 SNPs from across the six genomes and one reference genome and (B) 1007 SNPs within the

Table 1 .
Presence/absence of genes of interest in the six genomes of Microcystis aeruginosa with, green = complete gene present, white = gene absence, and blue = fragmented coverage.
*Presence confirmed via PCR and gel electrophoresis (see Supporting Information).