Molecular Plant Responses to Combined Abiotic Stresses Put a Spotlight on Unknown and Abundant Genes

Environmental stresses such as drought, heat, and salinity limit plant development and agricultural productivity. While individual stresses have been studied extensively, much less is known about the molecular interaction of responses to multiple stresses. To address this problem, we investigated molecular responses of Arabidopsis to single, double, and triple combinations of salt, osmotic, and heat stresses. A metabolite profiling analysis indicated the production of specific compatible solutes depending on the nature of the stress applied. We found that in combination with other stresses, heat has a dominant effect on global gene expression and metabolite level patterns. Treatments that include heat stress lead to strongly reduced transcription of genes coding for abundant photosynthetic proteins and proteins regulating the cell life cycle, while genes involved in protein degradation are up-regulated. Under combined stress conditions, the plants shifted their metabolism to a survival state characterized by low productivity. Our work provides molecular evidence for the dangers for plant productivity and future world food security posed by heat waves resulting from global warming. We highlight candidate genes, many of which are functionally uncharacterized, for engineering plant abiotic stress tolerance.


Introduction
Scientific and intergovernmental reports warn of dangerous consequences of global warming for various aspects of human life, with considerable threats to plant productivity (Lobell and Field, 2007;Battisti and Naylor, 2009;Lobell et al., 2011;Deryng et al., 2014;Lesk et al., 2016;IPCC, 2018). Current temperatures are approximately 1 °C above pre-industrial Downloaded from https://academic.oup.com/jxb/article/71/16/5098/5842162 by Universitaetsund Landesbibliothek Duesseldorf user on 20 November 2020 levels; a further increase of global temperatures by 0.5 °C would increase the associated risks. An increase to 2 °C above pre-industrial temperatures might lead to the complete loss of many ecosystems (IPCC, 2018). The frequency of more extreme climatic events, such as droughts and heatwaves, will increase with additional elevations of global temperatures (IPCC, 2014). Drought and high temperatures significantly reduced cereal production by about 10% between 1964 and 2007 (Lesk et al., 2016). The ever-increasing threat of extreme weather that is concomitant to global warming-in particular the increased occurrence of very high temperatures-might lead to catastrophic decreases in crop productivity and result in worldwide famine (Bita and Gerats, 2013).
One-third of the earth's landmass is dominated by arid and semiarid weather conditions (Huang et al. 2010). Considering the foreseen climate changes and global warming effects, these arid areas are expected to expand massively (Park et al., 2018); consequently, many crop plants are expected to suffer from a wide range of new stress combinations. Arid climate is characterized by high temperatures. These not only lead to high evaporation rates and hence cause drought, but also favor the accumulation of salts in the surface soil layers close to the root zones of crop plants, causing salinity stress. Thus, plants growing in arid areas are frequently subjected to a combination of heat, drought, and salt stresses; based on the climate and nature of the soil in a specific geographical region, these stresses may occur individually or in combination. In order to develop successful cultivation strategies that secure the world's food demands under the foreseen global warming, it is crucial to understand plant responses to complex environmental challenges (Battisti and Naylor, 2009;Lobell et al., 2011;Deryng et al., 2014).
A plethora of studies have investigated plant responses to individual stress conditions at the physiological and molecular levels (Schenk et al., 2000;Kilian et al., 2007;Matsui et al., 2008;Shaik and Ramakrishna, 2013). This type of study has significantly improved our understanding of the mechanisms of plant responses to environmental stresses. To apply this knowledge to conditions outside of the laboratory, we urgently need to elucidate the interactions that delineate plant responses to natural, multi-stress environments (Mittler, 2006;Atkinson and Urwin, 2012;Prasch and Sonnewald, 2015;Pandey et al., 2015;Zhang and Sonnewald, 2017).
Few studies exist on plant responses to combined stresses. Fortunately, several groups have started to overcome the technical issues that are encountered in designing the experimental set-up of combined stress factors under controlled conditions in the laboratory. Thus, in addition to two pioneering studies by Rizhsky et al. (2002Rizhsky et al. ( , 2004 on the effects of combined heat and drought stresses on gene expression profiles in tobacco and Arabidopsis, several recent genome-wide studies have made a shift to investigate plant responses to combined stresses (Rasmussen et al., 2013;Prasch and Sonnewald, 2013;Sewelam et al., 2014;Barah et al., 2016;Georgii et al., 2017;Shaar-Moshe et al., 2017;Osthoff et al., 2019). These studies suggested that plant molecular response to different individual stresses could not be used to predict the plant responses to combined stress treatments.
Here, we investigated molecular responses of the model plant Arabidopsis to single, double, and triple combinations of salt, osmotic, and heat stresses. Our results reveal a dominant effect of heat stress over salinity and osmotic stresses with regard to the global changes in gene expression and relative metabolite levels. These findings turn more attention towards the profound effects of heatwaves that are unavoidable in the current global warming era. We found that plants experiencing combined abiotic stresses reprogram their transcriptional machinery to down-regulate the expression from the most abundant genes, most probably as a strategy to save resources for defense and survival. Relative metabolite levels determined during the single stresses and their combinations highlighted the coordinated modulation of specific amino acid levels, especially in treatments including heat stress. Furthermore, we found that different compounds are used as compatible solutes depending on the nature of the abiotic stress.

Experimental design and stress treatment scheme
Seeds of wild-type Arabidopsis ecotype Columbia-0 (Col-0; WT) were sown on soil. After stratification at 4 °C for 2 d, the seeds were transferred to a long-day (16 h light-8 h dark) growth chamber using a mix of Spectralux®Plus NL 36W/840 (Radium) and Fluora L 36W/77 (Osram) bulbs with a light intensity of approximately 100 µE m −2 s −1 and a temperature regime of 22 °C day-18 °C night temperatures. Twoweek-old seedlings were transplanted to new soil. After another 14 d of growth, plants were subjected to the stress treatments (Fig. 1A).
To simulate the effects of complex field environments on Arabidopsis plants, we used well-established analog artificial stress conditions. These conditions lead to low variability and allowed us to tightly control stress level and onset and to grow plants in limited space (Verslues et al., 2006;Lawlor, 2013). As a proxy for drought stress, we applied mannitol, an osmoticum that lowers the water potential of the soil, making it difficult for the plants to extract water (Verslues et al., 2006). To simulate salt stress, we watered the plants with NaCl (Munns and Tester, 2008). For heat stress we exposed the plants to a defined elevated temperature that does not cause heat shock in Arabidopsis. For the salt stress treatment, plants were watered with 150 mM NaCl solution 16 h before sampling. The osmotic stress treatment was applied by watering plants with 200 mM mannitol 16 h before sampling. For the heat treatment, plant trays were transferred to a growth cabinet pre-adjusted at 35 °C for 4 h before sampling. We selected a sampling time of 16 h after stress initiation for salt and mannitol treatments because the plant response to these treatments peaks between 12 and 24 h. In contrast, 4 h was chosen for heat, as plant response to heat is reported to peak from 2 to 4 h (Kilian et al. 2007, Sewelam et al., 2014. To avoid the circadian clock effects that could interfere with the interactions between our treatments, we started the heat treatment 3 h after light onset, and we collected all samples at midday when the plant's metabolism is at a steady state. Double treatment of salt and mannitol stress was applied by watering plants with a mixed solution containing final concentrations of 150 mM NaCl and 200 mM mannitol. For double treatments of heat with NaCl or mannitol, a fraction of the plants that were watered with NaCl or mannitol were moved after 12 h to a growth cabinet pre-adjusted at 35 °C for an extra 4 h, resulting in a total stress time of 16 h before sampling (Fig. 1B).
For multiple stress treatments, a patch of the plants treated with the double treatment of NaCl 150 mM and 200 mM mannitol were transferred after 12 h to a growth cabinet pre-adjusted at 35 °C for an extra 4 h, resulting in a total stress time of 16 h before sampling (Fig. 1A). To apply a multiple stress treatment similar to field conditions as much as possible, Arabidopsis plants were treated with salt and mannitol 1 h before the light was turned off, and were kept for 12 h before starting the heat treatment. After 3 h of light, plant trays were transferred to a 35 °C growth cabinet for an extra 4 h, resulting in a total treatment time of 16 h (Fig. 1A). Rosettes were harvested and immediately frozen in liquid nitrogen. Tissue was ground into fine powder in the presence of liquid nitrogen, and then kept at −80 °C until use. Frozen samples were used for RNA extraction and metabolite analysis. For malondialdehyde (MDA) quantification, samples were collected 4 d after the application of the treatments. All experiments were performed in biological triplicates.

MDA measurements
Lipid peroxidation was measured by determining the amount of MDA. MDA concentration was estimated by the method of Heath and Packer (1968). A sample of 0.5 g fresh leaves was extracted in 10 ml 5% (w/v) trichloroacetic acid and the homogenate was clarified by centrifugation at 1500 g for 10 min. The supernatant (2 ml) was mixed with 2 ml of 0.67% (w/v) thiobarbaturic acid, incubated at 95 °C in a water bath for 20 min, and cooled immediately. Absorbance was read at 532 and 600 nm. MDA concentration (µM g −1 fresh weight) was calculated using the extinction coefficient 155 mM −1 cm −1 . The calibration curve was constructed in the concentration range of 0.2-2 µM.
The results were analysed using one-way analysis of variance (ANOVA). Significant differences among the means were identified using the Tukey HSD test.

RNA extraction and sequencing
For RNA extraction, rosette leaves from wild-type plants were harvested at the sampling times indicated in Fig. 1A and frozen in liquid nitrogen. Total RNA was extracted with the SV Total RNA Isolation System (Promega, Madison, WI, USA). The integrity of RNA was analysed by gel electrophoresis, and its concentration was measured using a NanoDrop (ND-1000) spectrophotometer. Libraries were prepared using the TruSeq RNA Sample Prep Kit v2 (Illumina, San Diego, CA, USA) and quantified with a Qubit 2.0 (Thermo Fisher Scientific, Waltham, MA, USA). Samples were multiplexed with 12 libraries per lane and sequenced in paired-end mode (Rapid Run, 100 bp read length) on an Illumina HiSeq 2500 platform.

Read mapping and mRNA-Seq data analysis
After successful quality control with the Fast QC software (v0.11.5), Illumina reads were quantified by mapping against the Arabidopsis reference transcriptome (Araport 11 representative CDS, Araport11_ genes.201606.cds.fasta.gz, retrieved from www.araport.org) using Kallisto v 0.45.1 (Bray et al., 2016) in default mode with 50 times bootstrapping (option, b 50). Kallisto provides normalized transcript abundance as transcripts per million (TPM). Differential transcript abundance comparing each of the treatments against the control were determined using the likelihood ratio test implemented in sleuth v0.30.0 (Pimentel et al., 2017). A significance threshold α of 0.05 was chosen after adjusting P-values (thereafter termed 'q-value') for multiple hypothesis testing via the Benjamini-Hochberg correction (Benjamini and Hochberg, 1995) as implemented in sleuth. The 'b value' supplied by sleuth was used as adapted log2-fold change between the respective treatment and control. k-means clustering was done with function k-means in R choosing eight clusters (k=8), as the clustering analyses that included more than eight clusters showed no new expression patterns. The inputs for k-means were average TPM of each treatment (n=3) after filtering for significantly differentially expressed transcripts (q≤0.05) and using a minimum foldchange threshold of 1. All statistical analyses were performed in the R environment for statistical computing (v3.6.1 provided by the CRAN project; R Core Team, 2018). Gene annotation and categorization was retrieved from publicly accessible databases, including Araport (www. araport.org/data/araport11), gene ontology (GO) annotation (www. arabidopsis.org), Mapman categorization (August 2012, www.plabipd. de/portal/web/guest/mapman), transcription factors (www.plntfdb.bio. uni-potsdam.de/v3.0), and genes involved in lipid metabolism (www. aalip.plantbiology.msu.edu/downloads).

Gas chromatography-mass spectrometry analysis
Extraction and analysis by gas chromatography-mass spectrometry (GC-MS) was performed using the same equipmental set-up and exact same protocol as described in Lisec et al. (2006). Briefly, frozen ground material was homogenized in 300 μl of methanol at 70 °C for 15 min and 200 μl of chloroform followed by 300 μl of water was added. The polar fraction was dried under vacuum, and the residue was derivatized for 120 min at 37 °C (in 40 μl of 20 mg ml −1 methoxyamine hydrochloride (Sigma-Aldrich, cat. no. 593-56-6) in pyridine followed by a 30 min treatment at 37 °C with 70 μl of N-methyl-N-(trimethylsilyl)trifluoroacetamide (MSTFA reagent;Macherey-Nagel, cat. no. 24589-78-4). An autosampler Gerstel Multi-Purpose system (Gerstel GmbH & Co.KG, Mülheim an der Ruhr, Germany) was used to inject 1µl of the samples in splitless mode to a chromatograph coupled to a time-of-flight mass spectrometer system (Leco Pegasus HT TOF-MS; Leco Corp., St Joseph, MI, USA). Helium was used as carrier gas at a constant flow rate of 2 ml s −1 and GC was performed on a 30 m DB-35 column (capillary column, 30 m length, 0.32 mm inner diameter, 0.25 μm film thickness, PN: G42, Agilent). The injection temperature was 230 °C and the transfer line and ion source were set to 250 °C. The initial temperature of the oven (85 °C) increased at a rate of 15 °C min −1 up to a final temperature of 360 °C. After a solvent delay of 180 s, mass spectra were recorded at 20 scans s −1 with m/z 70-600 scanning range. Chromatograms and mass spectra were evaluated using Chroma TOF 4.5 (Leco) and TagFinder 4.2 software. Metabolites were annotated based on a retention index calculation with deviation <5% and compared with the reference data of the Golm Metabolome Database, http://gmd.mpimp-golm.mpg.de (Luedemann et al., 2008). The results were statistically analysed using Student's t-test, using a threshold of P<0.05 between the samples of the control and treated plants.

Heat map and principal component analysis
Heat maps were created by MeV software (http://mev.tm4.org), using the log2 of the metabolite fold changes with respect to the control. Each cell in the heatmap represents the log2 average fold change value. Hierarchical clustering was applied to rows and columns using Pearson's correlation coefficient (>0.80). Principal component analysis was performed with the function prcomp implemented in R.

Plant phenotypes and stress level under individual and combined stress treatments
To study the interacting effects of different combinations of environmental stresses, 4-week-old Arabidopsis plants were subjected to single, double, and triple combinations of salt, osmotic, and heat stresses as shown in Fig. 1A. The applied treatments were selected and designed to simulate the environmental conditions in arid areas, where high temperatures cause high evaporation rates, leading to water shortage and accumulation of salts. We applied mild stress doses to impose physiological changes close to those occurring in natural plant habitats. Five days after the stress treatments, plants morphology (Fig. 1B) and the level of lipid peroxidation (Fig. 1C) were assessed. Among the single treatments, mannitol caused the highest damaging effect on plant growth and the highest lipid peroxidation level. The double treatments led to damaging effects stronger than any of the individual stresses, in which all plants showed reduced growth. After the single heat treatment, the plants recovered and grew normally, but the heat dose intensified the negative effects of salinity and osmotic stress on plant growth. Out of the double treatments the combination of salt and osmotic stresses presented the most pronounced damaging effects. The triple combined treatment caused the most sever effects compared with all other stress treatments ( Fig. 1B, C).

The genome-wide transcriptomic data indicate a dominant effect of heat stress
We investigated genome-wide changes in the Arabidopsis transcriptome as a result of the interacting effects of abiotic stresses by RNA-seq. A principal component analysis (PCA) showed that the first principal component (PC1) separated the expression patterns of the treatments that include heat stress from the other treatments, explaining 46% of the gene expression variation ( Fig. 2A). The second principal component (PC2), which explained an additional 13.4% of the variation, separated the patterns after treatments that include mannitol stress from all the others ( Fig. 2A). Out of 27 655 Arabidopsis genes represented in our study, 19 244 genes were differentially expressed (up-or down-regulated; P<0.05) by at least one of the seven treatments compared with the control (see Supplementary Dataset S1 at JXB online). Comparisons of qRT-PCR data with results from genome-wide studies indicated that relative expression changes of <2-fold change (FC) can be unreliable (Lee et al., 2005), and thus we set a cut-off value of log 2 of 1, which corresponds to 2-FC in expression with regards to the control.
At a cut-off value of 2-FC we found that 3529 genes were differentially expressed by at least one stress treatment compared with control (see Supplementary Dataset S1). Stress treatments including heat significantly affected the expression of the largest numbers of genes (Fig. 2B). Heat treatment significantly induced the expression of 1107 genes and reduced the expression of 697 genes. Mannitol treatment led to the induction of expression of 229 genes, and repressed the expression of 91 genes. Salt stress induced the expression of 63 genes, and repressed the expression of 22 genes. The effects of the single stresses could be viewed at first as having synergistic effects, in which the numbers of genes induced or repressed by the combined treatments were considerably higher than those of the individual treatments. The triple combination of stresses affected the expression level of the largest number of genes, with 1222 up-and 1071 down-regulated genes (Fig. 2B).

Different stress combinations cause significant reprograming of gene expression profiles: k-means clustering
To find specific expression patterns that may untangle the interactions between different stress treatments, we clustered the data by k-means clustering over the 3529 genes that were differentially expressed ≥2-FC ( Fig. 3; Supplementary  Table S1). We obtained eight clusters that pinpoint candidates of the abiotic stress response with potential for engineering stress-tolerant plants.

Cluster 1: osmotic stress and heat have antagonistic effects on gene expression
Cluster 1 included 281 genes, of which overall expression was not changed by salt treatment, while osmotic treatment caused induction of the majority of these genes, and heat repressed the expression of a large number of them (Fig. 3, Supplementary Table S1). This pattern suggests an antagonistic interaction of the osmotic and heat treatments. The double stress combinations show that salt treatment does not change the expression patterns induced by the individual osmotic or heat treatments. This suggests a lack of interaction between the effects of salt with either mannitol or heat on the expression profiles of the genes in this cluster. It is noteworthy that the combined osmotic and heat treatments largely reprogramed the gene expression patterns observed by the individual treatments, confirming the antagonistic interaction of these stresses on a large set of genes (120 genes, Supplementary Table S1). The triple stress treatment further validated this conclusion (Supplementary Table S1). These differences in gene expression may explain the opposing physiological effects between drought and heat (Rizhsky et al., 2002). As an example, under drought stress plants tend to close their stomata to avoid   (Rizhsky et al., 2002(Rizhsky et al., , 2004. The top 50 genes antagonistically affected by osmotic and heat treatments included many genes coding for late embryogenesis abundant (LEA) proteins, receptor-like kinases (RLK), glutathione S-transferases and a number of WRKY transcription factors (Fig. 4A). Many of these proteins are involved in responses to drought (Goyal et al., 2005;Tunnacliffe and Wise 2007;Tolleter et al. 2010). The decline of expression of these protective genes by heat may be, at least in part, responsible for the enhanced damage observed in plants treated with the combined osmotic and heat stresses with respect to the individual conditions (Fig. 1B, C).
Clusters 2, 3 and 5: a large number of genes are specifically repressed by heat Genes included in clusters 2 (912 genes), 3 (240 genes), and 5 (407 genes) showed repressed expression under all treatments that included heat stress ( Fig. 3; Supplementary Table S1). Cluster 3 showed the most pronounced repression. Single salt and osmotic treatments and their combination had minor effects on the expression levels of genes in cluster 2, while a slight repression was observed in clusters 3 and 5.
Cluster 2 represented about 26% of the total considered genes in our analysis. The most represented gene categories included those involved in DNA synthesis, secondary metabolism, biotic stress responses, and signaling events mediated by receptor-like kinases, calcium and G-proteins (see Supplementary Table S1). A large number of genes in this cluster (258 genes) were annotated as having unknown functions, making them candidates for further molecular analyses. Clusters 3 and 5 showed similar enrichment patterns to cluster 2. Especially cluster 5 contained many genes involved in cell wall metabolism (Supplementary Table S1), indicating that growth is negatively affected by heat.

Cluster 4: a small number of genes are induced by all treatments, and salt and osmotic stresses show synergistic interactions
Cluster 4 included a relatively small number of genes (55 genes) that were induced by all of the applied treatments (see Supplementary Table S1). Compared with the individual salt and osmotic treatments, the combined treatment of both showed a synergistic effect that changed the expression levels of most of the genes in this cluster (Fig. 3, cluster 4). Also, compared with the single and double stress treatments, the triple stress treatment caused a pronounced increase in the levels of a large number of genes. Salinity has two effects, early osmotic and relatively late ionic phases (Munns and Tester, 2008). Most probably, the genes of cluster 4 are involved in the osmotic responses of salinity as these genes are common between salt and osmotic stresses.
Clusters 6, 7, and 8: almost half of the differentially expressed genes, which include many oxidative phosphorylation components, are specifically induced by treatments including heat Genes in clusters 6, 7 and 8 were specifically up-regulated by treatments with a heat stress component (Fig. 3, see Supplementary Table S1). The only difference between these three clusters was the degree of the effect. Cumulatively, these three clusters included 1634 genes that represented more than 46% of the total considered genes in our clustering analysis. A large number of these genes (around 25%) are of unknown function; others are mostly involved in protein degradation or encode HSPs (Supplementary Table S1).
We found that 71 genes out of the 122 genes of the Arabidopsis mitochondrial genome were specifically induced by heat treatment and its combinations (see Supplementary  Table S1). Most of these genes encode for subunits of complex I, cytochrome oxidase, and the ATPase,

Heat shock transcription factors and heat shock proteins are mainly induced by heat and combinations of stresses including heat
Heat shock transcription factors (HSFs) and HSPs represent master regulatory and protective elements involved in the response to heat stress. Many genes coding for these proteins have already been incorporated in molecular breeding programs for sustainable agriculture (Ritossa, 1962;Swindell et al., 2007;Scharf et al., 2012;Niu and Xiang, 2018). As our transcriptome analysis showed that heat has a dominant effect on global gene expression, we specifically analysed the expression patterns of all annotated Arabidopsis HSFs (23 genes) and HSPs (156 genes) included in our study (see Supplementary  Table S2) under all stress treatments.
We found that 11 HSFs were significantly induced by treatments including heat stress (Fig. 4B). In addition, HSFA7A was significantly repressed by salt and osmotic stresses. Notably, HSFA6A was strongly induced by all treatments that included salinity and osmotic stress. HSFA6B showed specific induction by the triple combinations of stresses. Three HSF genes, HSFA1E, HSF3, and HSFA5, were repressed by treatments including a heat stress component. The remaining four HSF genes, HSFB3, HSFA1A, HSFA9, and AT4G18870, showed no responsiveness to any of the applied stress treatments (Fig. 4B). The HSF genes that were repressed by heat and/or other abiotic stress represent interesting targets for functional characterization not only to provide new insights into the plant response to abiotic stresses, but also for their promising incorporation in molecular breeding programs to enhance plant fitness.
In our analysis we also found that most members of all HSP families-in particular sHSPs and HSP70-were specifically induced by treatments that included heat stress (Fig. 4C), indicating that high temperature is the most severe environmental stress causing serious damage to cellular proteins. Notably, BIP3, an endoplasmic reticulum HSP70, was the only HSP gene repressed by heat. Instead, we found that BIP3 was induced by salt stress (Fig. 4C); BIP3's involvement in the unfolding protein response caused by salinity was recently discussed (Henriquez-Valencia et al., 2015;Guan et al., 2018).

Abiotic stress combinations induce intensive reprograming of the expression profiles of various transcription factor families
In all living organisms, transcription factors (TFs) represent the master control components that dynamically alter the cellular transcriptome, which in turn adapts metabolism and phenotype to the surrounding environment (Riaño-Pachón et al., 2007;Mitsuda and Ohme-Takagi, 2009). Here, we investigated the interacting effects of various stresses on the expression profiles of the 1187 TF genes annotated in the Arabidopsis genome (see Supplementary Table S3). As genes encoding TFs frequently do not show high expression changes compared with non-regulatory genes (Cheng et al., 2007), we monitored the changes in the expression profile of the TFs at a cut-off value of 1.5-FC. With this threshold we found that 701 TF genes were significantly differentially expressed by at least one stress treatment compared with the control (Supplementary Table S3).
The AP2-EREBP, WRKY, bHLH, and MYP/MYP related TF families showed the highest responsiveness to the applied stress combinations. We found that the members of a given TF family did not exhibit similar expression patterns under a given stress condition (see Supplementary Table S3). As a case study, we analysed the changes in the expression patterns of the gene members of the WRKY TF family. These TFs are assembled into three main groups (I, II, and III) based on the number of WRKY domains and the features of their zinc-finger-like motif (Eulgem et al., 2000;Rushton et al., 2010). In our study, 29 out of 72 genes annotated as WRKY TFs were differential expressed by at least one of the seven stress treatments compared with control (Fig. 5A). After categorizing the genes into induced or repressed sets, we found that there is no specific expression trend among the members of the same groups or even subgroups of this TF family under the applied stress treatments (Fig. 5A). This analysis indicates that the categorization of TF genes based on their binding motifs does not correlate with their expected biological functions inferred from their responsiveness to different environmental conditions. This raises the historic question related to genome-wide studies, of whether proteins that are identified according to their sequence similarities with known proteins are actually correctly categorized (Riechmann, 2002).

Genes involved in cell cycle are down-regulated especially under heat stress
Various events of the cell cycle in eukaryotes are controlled by cyclin-dependent kinases and their associated activators, cyclins, and cell division cycle proteins (Komaki and Sugimoto, 2012). Our data showed that the expression of many genes coding for these proteins was reduced mainly by the stress treatments that included heat, and particularly under the triple stress treatment (see Supplementary Table S4). Furthermore, we found that many genes involved in DNA synthesis and repair were also down-regulated mainly by heat and combinations of stresses including heat, and many genes involved in cell division were specifically repressed by heat (Supplementary Table S4).

Almost 40% of the differentially expressed genes under stress treatments are of unknown function
In the Arabidopsis genome, 38% (10 415 genes) of the total number of genes are of unknown function. Here, we found that out of these, 1312 genes are significantly (P<0.05) differentially expressed ≥2-FC by at least one of the seven stress treatments compared with control (see Supplementary Table S5). Interestingly, this number represents 37% of the total number (3529) of the differentially expressed genes.
Ten genes attracted our attention as they show specific expression patterns under various stress treatments (Fig. 5B). Three genes, AT4G19430, AT4G23493, and AT2G38465, were specifically highly induced by heat, suggesting that they may play an important role in thermotolerance. Two genes, AT5G53710 and AT1G62515, showed high induction by all individual stress treatments. Two genes, AT4G18253 and AT1G65500, were specifically induced by osmotic stress. Finally, AT5G01015, AT2G27402, and AT3G22235, showed significant repression by all treatments, with very strong repression under the combined triple stress treatment.

The transcription of most abundant genes is downregulated under stress
Abundantly expressed genes are usually ignored or underestimated in analyses that only use FC to evaluate changes in gene expression. The analysis of genes with transcript abundances higher than 100 TPM, regardless their FC expression levels (see Supplementary Table S6 ), indicated that around 70% of these genes were repressed by the stress treatments. This observation most probably reflects a strategy of the stressed plants to save resources for other urgent tasks required for defense and survival. Out of the considered 1397 abundantly expressed genes, 336 genes were differentially expressed by salinity (out of them 289 genes were repressed), 296 genes were differentially expressed by osmotic stress (out of them 226 were repressed), 1173 genes were differentially expressed by treatments that included heat (out of them 731 were repressed), and 1332 were differentially expressed by the triple stress treatment (out of them 962 were repressed) (Supplementary Table S6).
When we only considered the relative FC in expression levels, we inferred that the salt treatment has a negligible effect on expression patterns of Arabidopsis genome (see Supplementary  Table S6). Nevertheless, by considering the TPM we found a hidden effect of salinity on the expression patterns of abundantly expressed genes. Among the salinity-responsive genes, photosynthetic genes were the most significantly repressed (Supplementary Table S6; Fig. 5C).
The expression of the most abundant photosynthetic genes is highly down-regulated under stress Out of the considered abundant genes, around 100 genes are involved in photosynthesis (see Supplementary Table S6). If only the relative FC expression levels are considered, 90 of those photosynthetic genes would be disregarded (Supplementary  Table S6). The changes in expression levels of these genes are actually significant, with a general trend of down-regulation. Analysis of the top 10 most abundant photosynthetic genes (Fig. 5C) showed that most of these genes are repressed by heat and an intensification of this response is observed when other stresses simultaneously act with heat.
The most abundant photosynthetic genes encode the subunits of ribulose bisphosphate carboxylase/oxygenase (Rubisco) and the light harvesting complexes (LHC) (see Supplementary Table S6). Indeed, on Earth, Rubisco is the most abundant soluble protein (Ellis, 1979;Raven, 2013;Bar-On and Milo, 2019), and light harvesting a/b binding proteins are suggested to be the most abundant membrane-bound proteins (Xu et al., 2012). The expression level of the Rubisco gene RBCS-1A was repressed by the seven stress treatments, with the highest repression occurring under the combined triple treatment (Fig. 5C). In accordance with the reduction in the transcriptional expression of Rubisco, its activity was also found to be inhibited under drought and heat stresses in different plant species (Fahad et al., 2017;Haworth et al., 2018). We also found that LHCII genes are also repressed by stress treatments, particularly by heat and its combinations (Fig. 5C). In addition, our results showed that two LCHII genes, OHP1 and OHP2, are abundant and repressed by the stress combinations (Supplementary Table S6). Arabidopsis mutant lines of both genes, ohp1 and ohp2, showed severe reduction in growth, decrease in pigmentation and abnormal thylakoid architecture (Hey and Grimm, 2018;Li et al., 2019). Probably, the downregulation of these important genes explains, at least partially, the plant growth reduction observed under abiotic stresses (Fig. 1B).

Treatments that include heat stress cause downregulation of most ribosomal genes
Among the abundantly expressed genes are those involved in the synthesis of ribosomal proteins. We found that 234 abundant genes coding for ribosomal proteins were significantly repressed by treatments involving heat, especially by the triple stress treatment (see Supplementary Table S6). Most of these genes are predicted to localize to the cytosol (Supplementary  Table S6), where the bulk of the protein synthesis occurs. It was reported that stress does not dramatically modify the basal expression patterns of most genes encoding cytoplasmic ribosomal protein (Sormani et al., 2011). However, considering the high expression levels of the ribosomal genes, small changes in the relative expression levels represent big changes in protein amounts, which most probably is transduced in significant functional effects. This observed down-regulation of expression of the abundant ribosomal proteins parallels the downregulation of expression of the most abundant proteins in plant cells, Rubisco and LHC (Supplementary Table S6; Fig. 5C).

The profile of specific metabolites is affected by the abiotic stress treatments
We further investigated the effects of the different stress treatments on the profiles of several metabolites in leaves of Arabidopsis using GC-MS (see Supplementary Table S7). A PCA revealed that the PC1 separated the samples of the treatments that include heat stress from all other treatments, explaining 42.8% of the total variance within the dataset (Fig. 6A). Other components of the analysis did not provide any further clear difference between the treatments.
Our results revealed that specific amino acids were coordinately affected by treatments that included heat stress. Glycine, alanine, valine, isoleucine, lysine, tyrosine, and phenylalanine clustered together in a Pearson correlation heatmap and show increased levels under treatments with a heat stress component (Figs 6B, 7). We also observed that the levels of the osmoprotective amino acid proline were specifically higher following salt treatment during single or combined application with other stresses (Fig. 7). Treatment with mannitol also enhanced the levels of proline but with less intensity (Fig. 7). Accumulation of proline was also recorded in different plant species under salinity and drought stresses (Szabados and Savouré, 2010).
Soluble sugars are highly sensitive to environmental stresses, as they represent the supply of carbohydrates from source to sink organs and are also involved in osmoprotection. Our results showed that raffinose, galactinol, and trehalose, which all have osmoprotective properties, were highly increased by all treatments that included heat stress, while other treatments had a positive but much lower influence (Figs 6B, 7). Our transcriptomic data indicate that the expression of three genes encoding enzymes involved in the biosynthesis of galactinol and raffinose (Go1S1, -2, and -3) are enhanced under the imposed stress treatments (see Supplementary Table S1, cluster 6). In line with these results, transgenic Arabidopsis overexpressing GolS2 have enhanced levels of galactinol and raffinose and are more tolerant to drought and salinity stresses (Taji et al., 2002). Galactinol and raffinose also protect plants from oxidative damage (Nishizawa et al., 2008). We also found significant increment in the levels of fructose, a free sugar with osmoprotective function, specifically under treatments that include osmotic stress (Fig. 7).

Heat has a dominant effect on global gene expression
This study focused on molecular plant responses to combined global warming-associated abiotic stresses. The data presented here provide evidence that high temperature strongly affects global gene expression and metabolite profiles and governs the plant response even in combination with salinity and osmotic stresses.
The observed dominant effects of heat on the expression of a large number of genes (Figs 2, 3) can be mainly ascribed to the hypersensitivity of plant membranes to high temperature, which rapidly affects the structures and functions of biological membranes (Niu and Xiang, 2018). As most of the biological functions of the plants require molecular interactions through cellular and organellar membranes (Suzuki et al., 2012;Pospíšil and Prasad, 2014), the plant seems to quickly reorganize its transcriptome to rapidly and efficiently cope with the deleterious effects of heat stress. In addition to this, high temperature is mostly responsible for protein misfolding and enhanced protein degradation (McLoughlin et al., 2019). In line with this, we found that most members of all HSP families are specifically induced by treatments that include heat stress, indicating that high temperature is the most severe environmental stress causing serious damage to cellular proteins. Longer heat waves accompanying the global climate changes are expected to dramatically exacerbate these effects (IPCC, 2014;Haworth et al., 2018).
We also found that the mitochondrial activity is highly responsive to temperature stress (Supplementary Table S1). The increased abundance of the mitochondrial transcripts observed could represent an up-regulation of gene expression to replace damaged and degraded proteins. In accordance with this, a proteome analysis of cauliflower mitochondria indicated that the abundance of some oxidative phosphorylation components specifically increases in response to high temperature (Rurek et al., 2018). Changes in complex I composition most probably influence the redox homeostasis of mitochondria, further impacting the expression of antioxidant defense components and stress-responsive genes in a retrograde fashion. Furthermore, we found that genes involved in controlling cell cycle, cell division, DNA synthesis, and DNA repair are downregulated especially under heat stress (Supplementary Table  S4). These observations suggest that under high temperatures, plants tend to arrest their cell cycles. While this response can lead to reduced yield in crops, it has an important adaptive significance as a survival strategy (Kitsios and Doonan 2011).

Abundant and unknown genes in the spotlight
In our study, we draw attention to two categories of genes that are usually disregarded in genome-wide studies: those annotated as of unknown function and those having high transcript levels under normal growth conditions.

As a tradeoff strategy, plants down-regulate the transcription of most abundant genes under abiotic stress
Our analysis indicated that abundantly expressed genes, such as those encoding photosynthetic and ribosomal proteins, are down-regulated specifically in response to heat stress (Supplementary Table S6; Fig. 5C); an intensification of this response was observed when other stresses acted simultaneously with heat. These responses most probably reflect a strategy of the plants to save resources and energy that can instead be invested into defense mechanisms. Maintaining an efficient use of limited resources, such as nitrogen, through the reduction of expression of photosynthetic proteins represents an acclimation strategy to changing environmental factors (Pao et al., 2019). The reduction in photosynthesis and growth rates of plants observed in many studies (Fahad et al., 2017;Haworth et al., 2018;Huang et al., 2019) may be, at least partly, due to this plant strategy, and may not only be a direct consequence of the damaging effects of the applied stresses. However, it cannot be ruled out that the decreased expression of abundant transcripts for photosynthetic proteins could result from feedback control of these genes, as carbon fixation is inhibited by closed stomata and a reduced demand for photosynthates due to reduced growth. In addition to the down-regulation of genes coding for photosynthetic proteins, our results indicate that the expression of a number of genes involved in the degradation of chloroplastic proteins, such as FtsH6, FtsH8, and FtsH12, was up-regulated by heat and its stress combinations (see Supplementary Table S8). FtsH6 was reported to be involved in degradation of the PSII light harvesting proteins Lhcb1 and Lhcb3 under senescence and high-light conditions (Zelisko et al., 2005;Sedaghatmehr et al., 2016). Indeed, highly regulated degradation of the LHCII represents a crucial action to avoid photochemical damage to chloroplasts and other cellular constituents under harsh conditions (Zelisko et al., 2005). This adaptive down-regulation of the photosynthetic machinery, especially under prolonged heat waves, will inevitably inhibit the overall plant performance, diminishing plant growth.
Some functionally uncharacterized genes are candidates for engineering plant abiotic stress tolerance Interestingly, 38% of the Arabidopsis genome is of unknown function, and 37% of the genes that were significantly affected by stress treatments were also of unknown function (Supplementary Table 5). The expression of a group of these genes of unknown functions was highly modified under specific treatments applied (Fig. 5B). This highlights them as candidate genes for functional characterization and for breeding programs seeking to increase stress tolerance in crops.

Different compounds act as compatible solutes depending on the abiotic stress applied
Our metabolite analysis indicated the accumulation of specific compatible solutes depending on the nature of the stress applied (Figs 6B, 7). This most probably reflects the differential activation of specific metabolic pathways and the necessity to avoid cellular toxic effects of given osmoprotectants.
These results therefore highlight the importance of analysing stress-associated metabolites under specific stresses in programs aimed at engineering stress tolerance.
We also found a coordinated increase in the levels of some free amino acids under heat stress (Fig. 6B), which most probably reflects the catabolism of proteins to ensure the survival of the whole organism at the expense of other cellular processes. In line with this scenario, catabolic pathways for several amino acids were found to be induced under abiotic stress conditions, suggesting their use as a respiratory substitute to compensate the down-regulation of photosynthesis (Araújo et al., 2011;Batista-Silva et al., 2019). It is also possible that the enhanced levels of the branched-chain and aromatic amino acids reflect an increased production of secondary metabolites during the stress responses (Vetter, 2000;Dixon, 2001). Nevertheless, high temperatures can also negatively impact secondary metabolite production. This is illustrated by the case of fruits of Vitis vinifera and Ribes nigrum, in which reduced downstream polyphenol metabolism at high temperatures results in increased levels of phenylalanine (Mori et al., 2007;Allwood et al, 2019).

Concluding remarks
The general down-regulation of highly transcribed genes and cell cycle genes, combined with increased protein degradation, suggests that under unfavorable environmental conditions, Arabidopsis enters an emergency state with arrested growth and enhanced activity of molecular mechanisms for survival, causing a major reduction of plant development and productivity (Fig. 8). It can be argued that these responses are particularly important for a wild plant in a natural ecosystem, where the effects of environmental stresses are likely to be more dramatic than for a crop plant. Nevertheless, it is well documented that environmental stresses also significantly reduce crop productivity (Bita and Gerats, 2013;Lesk et al., 2016). Our findings provide molecular evidence for the serious effects of global warming on world food security, and the necessity of developing strategies to sustain future crop productivity.

Supplementary data
Supplementary data are available at JXB online.
Dataset S1. Quantitative information and annotation for all reads mapped onto the reference genome of Arabidopsis. Table S1. List of all differentially expressed genes (3529 genes) at 2-FC analysed by k-means clustering. Table S2. Expression of Arabidopsis HSP genes (156 genes) under seven stress treatments compared with control. Table S3. Expression of all Arabidopsis TF genes (1187 genes) under seven stress treatments compared with control. Table S4. Expression of Arabidopsis genes involved in cell division and DNA repair/synthesis at 2-FC under seven stress treatments compared with control. Table S5. Expression of 1312 genes annotated as having unknown function, differentially expressed by at least one treatment. Table S6. Expression of 1397 genes with transcript abundances higher than 100 transcripts per million. Table S7. Profiles of several metabolites in leaves of Arabidopsis using GC-MS. Table S8. Expression of genes involved in protein degradation with >2-FC.