Gene Expression Profiles Suggest a Better Cold Acclimation of Polyploids in the Alpine Species Ranunculus kuepferi (Ranunculaceae)

Alpine habitats are shaped by harsh abiotic conditions and cold climates. Temperature stress can affect phenotypic plasticity, reproduction, and epigenetic profiles, which may affect acclimation and adaptation. Distribution patterns suggest that polyploidy seems to be advantageous under cold conditions. Nevertheless, whether temperature stress can induce gene expression changes in different cytotypes, and how the response is initialized through gene set pathways and epigenetic control remain vague for non-model plants. The perennial alpine plant Ranunculus kuepferi was used to investigate the effect of cold stress on gene expression profiles. Diploid and autotetraploid individuals were exposed to cold and warm conditions in climate growth chambers and analyzed via transcriptome sequencing and qRT-PCR. Overall, cold stress changed gene expression profiles of both cytotypes and induced cold acclimation. Diploids changed more gene set pathways than tetraploids, and suppressed pathways involved in ion/cation homeostasis. Tetraploids mostly activated gene set pathways related to cell wall and plasma membrane. An epigenetic background for gene regulation in response to temperature conditions is indicated. Results suggest that perennial alpine plants can respond to temperature extremes via altered gene expression. Tetraploids are better acclimated to cold conditions, enabling them to colonize colder climatic areas in the Alps.


Introduction
Temperature stress is affecting several developmental processes in the life cycle of flowering plants [1] and is considered a key constraint to the geographical distribution of species. Aside from distributional ranges, changes in average temperature can affect the phenology, defense capacity, growth, and development of plants e.g., [2][3][4][5][6].
Flowering plants evolved various adaptation strategies to survive and reproduce under adverse temperature conditions, such as plastic responses, which are observed to be positively correlated with colonization of novel habitats e.g., [20,21] and subsequent selection of fitting phenotypes over many generations [22][23][24][25][26]. In recent decades, a growing number of studies have been focused on phenotypic plasticity, one component of which is thought to be changes in gene expression patterns, and its evolutionary aspects e.g., [27][28][29][30][31][32][33]. A reliable subset of them focuses on the triggering role of environmental autotetraploid plants are facultative apomicts (aposporous), with varying proportions of sexual and asexual seeds, and colonize previously glaciated areas, i.e., the northern, central and eastern Alps as well as the northern Apennines and Corsica [84,85,88,91]. Tetraploid populations occur at higher elevations in the European Alps than diploids and exhibit a pronounced niche shift towards colder temperatures [87,88]. High alpine habitats are characterized by short growth periods and cold spells, eventually with nocturnal frost during flowering time (spring or summer). It was suggested that the niche differentiation between the cytotypes is associated with a combination of climatic conditions and reproduction mode, with the asexual taxa having a distributional advantage towards cooler conditions [92].
Previous studies revealed that the tetraploid cytotype originated 10-80 kyears ago [92], presumably by multiple and recurrent autopolyploidization events [86,93]. The genetic differentiation between cytotypes, regardless of reproduction mode, is very low and on a similar level within cytotypes (F st values~0.3 for both cytotypes) [89]. The epigenetic background of the species proposes a differentiation of the cytotypes, but also a correlation of abiotic environmental conditions with the epigenetic variation in natural populations and in experimental conditions [94,95]. Thus, a putative epigenetic background of the niche shift of the tetraploids in the Alps is indicated, with epigenetic variation being associated with elevation in natural populations and higher persistence under cold treatments [94,95]. Syngelaki et al. [96] highlighted the potential of phenotypic plasticity, with changes of growth parameters linked to DNA methylation patterns, for acclimation to environmental conditions. These experiments confirmed the different niche preferences of cytotypes in natural populations, as diploids grow better under warm conditions, while tetraploids perform better in cold treatments.
Herein, we employed diploid and tetraploid plants already used by Syngelaki et al. [95,96] and exposed them to different controlled temperature treatments, to assess the gene expression profiles of individuals. We aim to investigate whether cold temperature stress influences gene expression and try to gain further insights into the expression dynamics. A temperature-sensitivity of gene expression profiles is speculated, which could act as a rapid response towards stressful environments. We also investigated the differentiation of the gene expression profiles according to the ploidy level of the individuals. We hypothesize that the observed niche shift of the tetraploid cytotype has the physiological background of a better cold acclimation. Finally, we associated the gene expression results with DNA methylation, on a transcriptomic level, as a correlation between them would explain the potential of tetraploid R. kuepferi to adapt to cold conditions at higher altitudes during the postglacial recolonization of the European Alps.

Plant Material and Experimental Design
Diploid and tetraploid plants of R. kuepferi used in the present study were part of a long-term temperature stress experiment, as described in Klatt et al. [97] and Syngelaki et al. [95,96]. Plants were collected throughout the distribution range of the species in the European Alps [87] during the growing seasons of 2013 & 2014. Subsequently, the plants were re-potted in garden soil at the Old Botanical Garden of Göttingen University, where their ploidy level was defined via Flow Cytometry of silica gel dried leaf material collected in the field [88]. For the scope of the current experimental design, which was conducted from 2014 onwards, the plants were exposed to different temperature conditions during the sprouting and flowering period, while the rest of the parameters were kept equal. The settings for the transcriptome study were as in Syngelaki et al. [95,96] (Table 1). The temperature conditions were simulating the natural environment of the species in the Alps, including freezing during some nights. Altogether 262 individuals were categorized into four groups corresponding to their ploidy level and treatment: cold diploids (CD, 63 plants), warm diploids (WD, 79 plants), cold tetraploids (CT, 71 plants), and warm tetraploids (WT, 49 plants). At the peak of the flowering season in beginning of summer 2019, leaf tissue was collected simultaneously from three individuals per group, each originating from different natural populations, immediately frozen in liquid nitrogen and stored in −80 • C. Sample collection localities can be found in Supplementary Table S1.

RNA Extraction and Sequencing
Frozen leaf tissue was pulverized in liquid nitrogen and a maximum of 100 mg of powder was used for RNA isolation with the RNAeasy Plant ® Mini Kit (QIAGEN, Hilden, Germany) following the provided protocol. RNA quantity and quality were determined with Nanodrop, a Qubit TM with the RNA HS Assay Kit (ThermoFisher Scientific, Waltham, MA, USA) and an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Library preparation and sequencing with HiSeq 4000 (Illumina, San Diego, CA, USA) was conducted at the Transcriptome and Genome Analysis Laboratory of the Microarray & Deep-Sequencing Core Facility (UMG, Georg-August-Universität, Göttingen, Germany) producing 50 bp single end reads.

Bioinformatics
The quality of reads was assessed using FastQC v.0.11.4 [98]. Raw reads were trimmed with CutAdapt v.2.3 [99], removing adapters and bases with a phred score below 30 and removing reads shorter than 30bp after trimming. As there is no available genomic reference for R. kuepferi, the transcriptomes from all diploid individuals were pooled for a de novo assembly with Trinity v.2.8.6 and default parameters, except for max. memory was set to 50Gb [100,101] to produce a pseudoreference. The quality of the resulting assembly was verified with BUSCO v.3.0.2 [102] (Supplementary Figure S1). TransDecoder v.5.5.0 [103] was used to identify the longest Open Reading Frame per assembled contig. Coding sequences were annotated using the blastp algorithm under NCBI-BLAST v.2.10.0 [104] and the December 2020 release of Uniprot as reference. Annotation reports were produced using Trinotate v.3.2.1 [105]. Trimmed reads of each sample were mapped against the annotated pseudoreference individually using Bowtie2 v.2.3.5.1 with default parameters [106].
ClusterProfiler v.4.0.2 [113] was used for Gene Set Enrichment Analysis. This package currently only accepts single organism references via AnnotationDbi v.1.55.1 [114]. Of the currently accepted reference genomes, A. thaliana is most closely related to R. kuepferi [115]. Consequently, the pseudoreference was annotated again as described above using the updated TAIR 10 release [The Arabidopsis Information Resource (TAIR), http://www. arabidopsis.org; accessed on 11 March 2021] and the resulting annotations were employed in ClusterProfiler for a separate run for each cytotype, checking all subontologies and with several 1 Mio. permutations. For these analyses, the warm treatment was set as the control condition. Dot plots were generated with enrichplot v.1.13.1 [116].

qRT-PCR
Quantitative real-time RT-PCR was conducted to validate the differential gene expression revealed by bioinformatic analyses. The annotation reports of the pseudoreference were screened for possible genes of interest (GOIs), which are stated to be related, directly or indirectly, to DNA methylation and gene regulation [12,60,61,82,[117][118][119][120][121][122][123], as well as housekeeping genes [124,125]. Primers were designed for two methyltransferases (CAMT3, PMT2), two demethylases (JM706, JMJ25), the AGO4B gene, which is participating in the RNA-directed DNA methylation (RdDM) pathway in Arabidopsis and rice [121] and the housekeeping gene Actin (Supplementary Table S2). Primer specificity was tested with a touch-down PCR; products were sequenced and compared to the respective gene in the pseudoreference.
Complementary DNA synthesis and qRT-PCR was performed using QuantiTect Reverse Transcription Kit (QIAGEN, Hilden, Germany) and the Rotor-Gene SYBR Green PCR kit (QIAGEN, Hilden, Germany) in QIAGEN Rotor-Gene Q cycler, equipped with Q-Rex Software and following the instructions of the manufacturer for two cycling steps and a 45 cycles PCR program for three technical replicates per sample. To evaluate the differential gene expression, the amplification performance as the log of fold change was calculated with the ∆∆Ct method [126] in Excel 2016, using the housekeeping gene Actin as endogenous control for normalization and warm diploid individuals were considered the reference.

Pseudoreference and Mapping
In the current study, gene expression profiles of six diploid and six tetraploid individuals of R. kuepferi, under cold (stress) and warm (control) temperature treatments were explored. Sequencing of the samples resulted in a mean of 31,918,319 raw reads per sample with a mean of 31,707,039 reads retained after trimming. The assembly of the pseudoreference resulted in 71,444 transcripts, with 15,224 of them functionally annotated. Through Bowtie2 mapping, we obtained an average mapping rate of 93.08% per sample ( Table 2).

Differential Gene Expression
A total of 2617 significantly differentially expressed genes were identified between the four predefined groups. Among all groups, more genes were found to be down-regulated compared to up-regulated (Table 3). Gene expression was correlated with temperature, while the ploidy level of the plants under the same environmental conditions did not affect the gene expression strongly ( Figure 1). One WD individual (WD2443) seems to present gene expression patterns different from all the other samples. The same individual seems to be an outlier also in the qRT-PCR analysis.  Table  S1).

Gene Set Enrichment Analyses
Gene set enrichment was successfully assigned and statistically significant for 59 pathways in the diploids and 20 pathways in the tetraploids. Enriched gene sets with higher GeneRatios, which in ClusterProfiler are defined as 'count/setSize,' where 'count' is the number of genes that belong to a given gene set, while 'setSize' is the total number of genes in the gene set, are presented in Figures 2 and 3. Regarding Figure 2, an extract of all resulted enriched pathways is shown due to graphical purposes and a dotplot of all pathways is provided in Supplementary Data (Supplementary Figure S2). Overall, 25 and 13 pathways have been activated in diploid and tetraploid individuals, respectively. Pathways which are linked to the plasma membrane e.g., 'anchored component of plasma membrane', '(cation) transmembrane transporter activity' and the cell wall e.g., 'cell wall organization or biogenesis', 'plant-type cell wall', as well as the 'cold acclimation' and 'hydrolase activity' pathways, are activated in both cytotypes (Figures 2 and 3). However, pathways related to ion/cation homeostasis and enzymic activity, such as '(cellular) metal ion homeostasis' and 'protein serine/threonine kinase activity' are enriched only in diploids, while pathways such 'histone/chromatin modification' and 'lipid transport' are enriched only in tetraploid individuals.  Table S1).

Gene Set Enrichment Analyses
Gene set enrichment was successfully assigned and statistically significant for 59 pathways in the diploids and 20 pathways in the tetraploids. Enriched gene sets with higher GeneRatios, which in ClusterProfiler are defined as 'count/setSize,' where 'count' is the number of genes that belong to a given gene set, while 'setSize' is the total number of genes in the gene set, are presented in Figures 2 and 3. Regarding Figure 2, an extract of all resulted enriched pathways is shown due to graphical purposes and a dotplot of all pathways is provided in Supplementary Data (Supplementary Figure S2). Overall, 25 and 13 pathways have been activated in diploid and tetraploid individuals, respectively. Pathways which are linked to the plasma membrane e.g., 'anchored component of plasma membrane', '(cation) transmembrane transporter activity' and the cell wall e.g., 'cell wall organization or biogenesis', 'plant-type cell wall', as well as the 'cold acclimation' and 'hydrolase activity' pathways, are activated in both cytotypes (Figures 2 and 3). However, pathways related to ion/cation homeostasis and enzymic activity, such as '(cellular) metal ion homeostasis' and 'protein serine/threonine kinase activity' are enriched only in diploids, while pathways such 'histone/chromatin modification' and 'lipid transport' are enriched only in tetraploid individuals.

Genes of Interest and qRT-PCR
A total of 38 genes of interest, identified by their involvement in DNA methylation, were significantly differentially expressed among the four predefined groups (Table 4). Table 4. List of differentially expressed genes related to DNA methylation in R. kuepferi leaves. Annotation is according to UniProt. The genes which are highlighted are employed as genes of interest (GOIs) for the real-time qRT-PCR.

Genes of Interest and qRT-PCR
A total of 38 genes of interest, identified by their involvement in DNA methylation, were significantly differentially expressed among the four predefined groups (Table 4).  The correlation of differential gene expression analysis with epigenetics was validated via qRT-PCR. For almost all individuals the expression of all five of the selected genes was down-regulated (Supplementary Table S3), corresponding to the differential expression analysis results of DESeq2. Only for one individual, an up-regulation of all the genes was detected. This individual is the same that separates from the rest of the samples under cold treatment (Figure 1).

Discussion
Plant acclimation to cold stress induces various cellular processes through a cascade of change in gene expression and protein synthesis e.g., [127,128]. It is estimated that between 4% and 12% of the transcriptome of A. thaliana changes after a respective period of several hours, days or weeks of exposure to chilling temperatures [129,130]. This differentiation in gene expression combined with the observation of different gene clusters being up-regulated during different times of the stress exposure indicates a hierarchy in the functional response, with signaling of harmful conditions or increasing freezing tolerance comes first [129,131,132], while circadian clocks being hypothesized to play an important role in general regulation [117,133]. Polyploidy may affect cold acclimation as polyploids are thought to perform better under cold conditions e.g., [54]. We analyzed here for the first time gene expression of a perennial alpine plant under different temperature conditions, and evaluated effects of different ploidy levels on the response to cold stress.

Ploidy Effects on Gene Expression and the Distribution Pattern
The ploidy level of the individuals per se does not appear to be a significant contributor to the observed differential gene expression (Figure 1). This differs from previous studies on DNA methylation patterns, mode of reproduction, and morphological traits of R. kuepferi which revealed significant ploidy effects [95][96][97]. The congeneric species of the Ranunculus auricomus complex also showed strong ploidy effects in gene expression profiles of ovules [134]. In our study, the main effects in gene expression changes are due to treatments, not to ploidy. In autopolyploids, transcription profiles can be influenced by a multitude of factors, which are caused by genome duplication e.g., dosage effects due to the presence of additional copies of genes [78]. However, autotetraploid rice does seemingly not show a genome-wide dosage effect on gene expression; likely because subfunctionalization maintains the functional pleonasm of duplicated genes and avoids any energy waste [135][136][137][138]. Regarding the WD individual, which is an outlier for both the gene expression profiles and qRT-PCR analysis, there was no indication of lower RNA quality during wet and dry lab manipulation of the specimen. Probably, the conditional plasticity and the gene-environment interaction [39,139] of the individual, as well as the micro site of origin of the natural population, which belongs to the sympatric zone of the two cytotypes [88], could play the major role for its observed phenotype.
To investigate further how the two cytotypes cope with stress conditions, gene set enrichment for each ploidy level was assessed. Several key regulatory pathways and their interactions have been documented previously e.g., [37,117,[140][141][142]. In R. kuepferi, diploids respond more intensively to temperature treatments compared to the tetraploids, as a higher number of gene sets is significantly differentially expressed (Figures 2 and  3). Hence, diploids appear to be more stressed by cold conditions than tetraploids. The cold stress treatments of our experiments are quite similar to the natural habitat conditions of the tetraploid R. kuepferi plants in the Alps [95]. As the cold adapted genotypes are thought to have a distinct advantage over non-adapted ones in frost-prone environments, such as high mountain areas [9], we hypothesize that the postglacial colonization of the Alps by autotetraploid populations of R. kuepferi [88] did rely to a large extent on the phenotypic variation towards the climatic conditions, pronounced also as differential gene expression [143]. Thus, the geographical parthenogenesis scenario, which interrelates the cytotypes with different ecological backgrounds [88,92], is further supported. Results support the simulation study of postglacial recolonization in the Alps, which identified the acclimation/adaptation of tetraploids to a colder climatic niche in higher and more northern parts of the Alps as one decisive character of the geographical parthenogenesis pattern [92]. Our results here suggest that this niche shift has a direct physiological background of cold tolerance and is less likely due other characteristics of alpine habitats (like a lower pathogen pressure, among others).

Functional Aspects of Gene Expression Related to Cold Acclimation
Although some of the differentially expressed gene sets overlap between ploidy levels, there are some characteristics for either group. In diploids, most pathways are linked to ion/cation homeostasis and activity and are suppressed. Additional pathways are connected to the plasma membrane, cell wall, and hydrolase activity. Similar pathways are present in tetraploid individuals, with a greater focus on membranes and cell wall. In contrast to diploids, the tetraploids suppressed only six pathways, mostly related to chromatin and histone modification. In both cytotypes the cold acclimation pathway is activated.
As access point of the cell, membranes are injured by adverse environmental conditions, yet their stability contributes to cell behavior and activity maintenance [16,117,144]. Pathways related to cell lipid composition, such as the ones activated in the tetraploids, play an important role in the maintenance of plasma membrane functionality e.g., [145,146] and are affecting the downstream expression of genes linked to resilience to lower temperature [147]. Maintenance of plasma membrane functionality is especially important for freezing tolerance in alpine plants, as a fluid membrane allows transfer of water from the protoplast into the intercellular space, where extracellular ice nucleation takes place, leaving the protoplast unfrozen [13]. We suppose that our short freezing treatments (−1 • C during three nights per week) induced these expression changes. However, we observed no apparent damage of leaves in cold treatments [96], and hence both cytotypes are tolerant to short-term freezing and thawing during the day.
Decreasing membrane fluidity, coupled with its interaction with the cell wall, is considered to be one of the first cold sensors [148][149][150]. The connection of membrane rigidification to cytoskeletal rearrangements, calcium influxes and kinases, acts as trigger for the subsequent low-temperature response e.g., [151,152]. Changes in the composition of the cell wall can strongly affect plant stress resistance [153] as stress can up-regulate the expression of e.g., expansins and xyloglycan-modifying enzymes, which are acting on cell wall remodeling [154]. The hydrolase activity pathways activated in both ploidy levels in R. kuepferi under cold stress may indicate their unique roles in cell wall modification [155,156]. The cytoskeleton is also affected by cold stress [150,157] and its interactions with membranes and the cell wall play a distinctive role in cold stress tolerance [18,47,158]. The 'cell wall organization or biogenesis' pathway, which is activated in both diploid and tetraploid R. kuepferi individuals (Figures 2 and 3) is relevant to the cytoskeleton and its modifications under cold stress.
The pathways of ion/cation homeostasis and activity, as well as the serine/threonine protein kinase activity are thought to be temporal and spatial events downstream of Ca 2+ signaling e.g., [159,160]. A change in intracellular calcium ion concentration is one of the earliest signaling events triggering the response of plants to cold stress [130,161,162] with Ca 2+ dynamics being detected within 40s after cold stress exposure [163]. Presumably Ca 2+ channels are primary sensors for cooling rate and Ca 2+ efflux transporters are absolute temperature sensors [164,165], while oscillations of Ca 2+ are linked to stomatal closure in Arabidopsis thaliana, as a response to cold stress [166]. Furthermore, the serine/threonine protein kinases are plasma membrane receptor-like kinases (RLKs), several of which are calcium-moderated and promote the expression of cold-responsive genes (COR) through the activation of a mitogen-activated protein kinase (MAPK) cascade [153,[167][168][169].
The intriguing suppression of the latter pathways, especially in diploid R. kuepferi, corroborates the hypotheses of tangled information encoded through Ca 2+ kinetics, which are constantly changing as a complex mechanism of stress response and are also influenced by 'cold memory', i.e., former exposure to cold stress conditions [130,170]. We hypothesize that suppressed pathways of ion/cation homeostasis in diploids are probably linked to stomatal closure and consequently reduction of CO 2 uptake/carbon gain [171]. This would explain the lower growth performance of diploids under cold conditions, as observed in the experiments of Syngelaki et al. [96].

Gene Expression Related to the Epigenetic Mechanism of DNA Methylation
Deciphering the epigenetic background of plants, which are exposed to abiotic stress, is a constantly developing field e.g., [43]. DNA methylation is correlated with histone proteins and their post-transcriptional modifications, as the conversion of these modifications to DNA methylation profiles is often thought to be more stable e.g., [75,172]. These interactions are associated with gene expression profiles and gene transcription in general in response to cold stress [173,174] and could encompass changes in chromatin structure and accessibility [175,176]. Chromatin remodeling has a putative function as a plant thermometer, representing a relatively direct connection between cold and gene expression [177,178]. In the present study, the pathways of histone and chromatin modifications are suppressed in the tetraploid individuals ( Figure 3). Additionally, several genes, directly or indirectly related to DNA methylation, are significantly differentially expressed in both ploidy levels (Table 4). This result corroborates previous results of methylation-sensitive AFLP screenings that the cytotypes exhibit different methylation profiles [94,95].
The differential expression of genes correlated with DNA methylation is further validated by the qRT-PCR results, where the two methyltransferases (CAMT3, PMT2), two demethylases (JM706, JMJ25) and the argonaute protein AGO4B were all found to be down-regulated, as expected from the bioinformatical analyses of the transcriptomes. This fits the overall pattern of loss of methylated MS-AFLP fragments after dramatic temperature changes [95]. Methylation patterns may conserve the transgenerational epigenetic memory of response to cold treatments and hence differential acclimation and adaptation of cytotypes [94,95].
To summarize, the present study demonstrates the responses of diploid and tetraploid R. kuepferi plants towards cold stress in their gene expression patterns. Although both ploidy levels activate genes related to cold acclimation, the gene set pathways differ between cytotypes, suggesting a better cold acclimation of tetraploids than diploids. Consequently, our results strongly support the hypothesis of a physiological background of the observed ecological and geographical differentiation patterns between cytotypes. Altogether, cold stress induces differentially expressed gene profiles and several gene set pathways are involved in the response, either being activated or suppressed. Seemingly, these parallel mechanisms invoke energy conservation and help individuals to survive in novel and/or extreme environments. Lastly, DNA methylation is indicated to contribute to the regulation of gene expression and may preserve a different epigenetic memory for the two cytotypes.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/genes12111818/s1, Table S1. List of individuals used for the laboratory work. Listed are the sample and group ID, ploidy level, treatment, as well as country, province and altitude of origin [88,97]; Table S2. List of primers used for qRT-PCR validation; Table S3. Results of qRT-PCR and Delta-delta analyses; Figure S1. BUSCO (Benchmarking Universal Single-Copy Orthologs) plot of the R. kuepferi transcriptome pseudoreference; Figure S2. All gene set pathways significantly enriched from the analysis of differentially expressed genes in diploid individuals of R. kuepferi.