Construction of differentially-expressed-gene libraries and hierarchical clustering
After screening and filtering all available abiotic-stress-related transcriptomic datasets (n = 945) from the Gene-Expression Omnibus (GEO) database, 500 individual transcriptomes (from 23 selected datasets) were analyzed for differentially expressed genes (DEGs) under 10 different stress conditions: cold, complete submergence, drought, heat, high light, osmotic, salt, partial submergence, UV (UV-B), wounding, and exogenous treatment with the ethylene precursor 1-aminocyclopropane-1-carboxylate (ACC) (Supplementary Data file 1). An overview of the complete analysis is presented in Fig. S1. Based on the kinetics of individual transcriptomic analyses, the lists of DEGs were combined in early (from 1 to 6 h) and late responses to stress conditons (from 12 to 24 h) (Supplementary Data file 2; see the “Methods” section). Prior to further analyses, a first validation step was performed. To confirm the suitability and biological relevance of each DEG library, six marker genes were selected per stress condition based on their empirically determined expression. Subsequently, the expression of these markers was assessed for each stress, taking into account temporal and spatial specificity (Supplementary Data file 3). All passed the first biological validation test with at least five out of six markers being differentially expressed, demonstrating both the suitability and the accuracy of the DEG libraries.
The number of DEGs obtained in each stress, tissue, and timepoint combination provided a first insight into the regulation of abiotic-stress responses (Fig. 1a, d). In general, the responses were balanced considering up- and down-regulation and number of DEGs between tissues and timepoints. A hierarchical clustering analysis (HCA) grouped certain stresses by DEG modules, suggesting physiological resemblances between them (e.g. salt and osmotic stress in roots and shoots; partial and complete submergence in roots; wounding and drought in shoots; and early exposure to UV and high light in shoots; Fig. 1b, c, e, and f).
a Number of DEGs detected in each of the studied conditions in shoot tissue for early and late time points. Up-regulated DEGs are colored in blue, down-regulated in red. b and c HCA of the different stress conditions by DEG modules in shoots for early b and late c time points. d Number of DEGs detected in each of the studied conditions in root tissue for early and late timepoints. Up-regulated DEGs are colored in blue, down-regulated in red. e and f HCA of the different stress conditions by DEG modules in roots for early e and late f time points. Branch length is based on cluster distance.
Analysis of stress-related transcriptional responses found clear spatial and temporal differences in gene expression. In general, root responses were more stable over time, with minimal changes in clustering between early and late responses. For example, the same stress clusters were identified during early and late stress responses in roots, suggesting consistent underlying mechanisms (Figs. 1, S2 and S3). In contrast, shoot responses exhibited more dynamic shifts between early and late stress conditions. Notably, stresses like osmotic, cold, and UV elicited stronger transcriptional responses in shoots, with a higher number of DEGs detected during the early exposure. However, late responses in shoots showed a shift, with fewer up-regulated genes and a higher presence of transcriptional repression as down-regulated DEGs, highlighting the temporal complexity of shoot stress responses (Figs. 1, S2 and S3). This supports the existence of a tissue-specific mechanism to respond to stress conditions, especially when maintained over time (for more details, see Supplementary Information file).
Support vector machine (SVM) clustering classifies genes into stress cores
To identify a set of central actors participating in stress signaling, i.e. as part of a stress gene core, we computed meta-p-values for all studied genes (see the “Methods” section). We took into account their transcriptional changes in all surveyed stress conditions in four different datasets: roots early, roots late, shoots early and shoots late, and used SVM Clustering—an unsupervised version of standard SVM. First, a frequency-based pre-classification was performed. The genes appearing as DEGs in at least five of the studied conditions were assigned to the positive class (class 1), while the rest was assigned to the negative class (class 0). Subsequently, the meta-p-value dataset containing the information for all 10 stress conditions for the complete set of around 12,000 genes, was re-classified using SVM Clustering. This analysis classifies genes depending on their distribution in the 10-dimensional space taking into account the distribution of meta-p-values (reflecting how statistically significant the expression changes are under each condition).
SVM Clustering categorized the vast majority of genes as not relevant (class 0; approximately 99%), coinciding with the frequency-based pre-classification (Fig. 2a). Around 5–30% of the genes pre-classified as relevant (32, 6, 82 and 14 genes for the early-root, late-root, early-shoot and late-shoot responses, respectively) were refuted after SVM Clustering and deemed not relevant (1 → 0). In contrast, a few genes pre-classified as not relevant(0) were included in the final SVM gene core (0 → 1) (2, 4, 7 and 0 genes from the early-root, late-root, early-shoot and late-shoot responses, respectively). The number of genes pre-classified as relevant but considered irrelevant by SVM Clustering highlights the discriminatory power of this classification approach, refining the pre-classification established on the distribution of the complete set of meta-p-values. Based on this, four sets of core genes, coined SVM gene cores, with significant transcriptional alterations in all the studied conditions, were identified: 118 genes for early responses in roots, 108 for late responses in roots, 185 for early responses in shoots and 74 for late responses in shoots (Fig. 2a).
a SVM-specific hyperparameters (cost and gamma) used for the SVM Clustering of each tissue and timepoint dataset. The number of genes forming each SVM gene core is highlighted in bold at the right. A pre-classification was performed first, whereby relevant genes appearing as a DEG in at least five of the stress conditions studied were assigned to the positive class (class 1). The results of SVM Clustering are defined as 0 → 0 (genes marked as not relevant and stated as such by the algorithm), 1 → 0 (genes marked as relevant but not stated as such by the algorithm), 0 → 1 (genes not marked as relevant but stated as such by the algorithm) and 1 → 1 (genes marked as relevant and stated as such by the algorithm). b Venn diagram representing the overlap between different SVM gene cores. c Overlap between the SVM gene cores. The value, positive (over-representation) or negative (under-representation), and the statistical significance of the overlaps are indicated. A hypergeometric test was used to detect the statistical significance of the overlaps, considering a p-value < 0.001 as significant due to the reduced number of genes in the comparisons. No adjustments were applied for statistical comparissons. d Venn diagram representing the overlap between the root and shoot gene cores.
The comparison between tissues (root versus shoot in all timepoint combinations; Fig. 2b, c), as opposed to comparisons between timepoints within a tissue (root early versus root late; shoot early versus shoot late), revealed a significant under-representation of overlapping genes. Therefore, considering gene-response composition, tissue specificity is stronger than temporal specificity, supporting the results obtained by the qualitative DEG classification (Figs. 1 and S3). For that reason, and due to the statistically non-significant differences between timepoints (Fig. 2c), the genes belonging to the SVM gene cores per tissue were combined, resulting in a final dataset of 207 genes forming the root gene core and 237 genes forming the shoot gene core after removal of duplicates (Fig. 2d). Despite the tissue specificity of the SVM gene cores, 19 genes are shared between the root and shoot cores, which encompass fundamental proteins with tissue-independent functions. These predominantly cover genes involved in cell-wall maintenance and membrane integrity (EXPANSIN A1 (EXPA1), LIPID TRANSFER PROTEIN 2 (LTP2), dehydrins, such as COLD-REGULATED 47 (COR47) and LOW TEMPERATURE-INDUCED 30 (LTI30), and BLUE COPPER BINDING PROTEIN (BCB))31,32,33 in addition to some uncharacterized genes (e.g. AT1G19380 and AT5G19875). The complete list of genes that form the different SVM gene cores as well as their overlap is in Supplementary Data file 4.
We studied the composition of the gene cores in terms of annotated biological functions (Gene Ontology (GO) enrichment analysis) and gene families. The most significant function enriched in the shoot core was the response to water (GO:0009415) and response to hypoxia appeared in roots (GO:0071456) (Figs. S2 and S4). This reflects that stressed shoots prioritize maintenance of water homeostasis, while roots mostly try to maintain normoxia. In addition, amino-acid transporters were enriched in the shoot core, while EXPANSINs, related to cell-wall remodeling, appeared on the forefront in the root core (for a complete GO and gene family analysis, see Supplementary Information file).
Protein networks related to the SVM gene cores
To further elucidate the functionality of shoot and root cores, we constructed a protein network representing both physical and functional interactions. Subsequently, a k-means clustering method was applied to obtain protein clusters based on known interactions in order to provide further evidence of their biological roles (Material and Methods).
For shoots, four clusters were identified (Fig. 3a). The blue and red clusters contained the largest number of proteins (28%). However, given its higher degree of connectivity, the blue cluster was considered to be a key cluster within the shoot core (Fig. 3a, b). To support the biological relevance of the different clusters, we performed a biological-processes GO enrichment for each cluster individually (Fig. 3c–f). As expected, the biological responses of the blue cluster largely overlapped with those of the overall shoot core (Fig. S4), with ‘response to water deprivation’ (GO:0009414) the most significant GO term (Fig. 3f). Three WRKY transcription factors appeared to act as central nodes in the interaction network, strongly interacting among themselves (WRKY33, WRKY46 and WRKY18) (Fig. 3b). In addition, MAP KINASE KINASE 9 (MKK9), a Mitogen Activated Protein Kinase (MAPK) protein, directly interacts with WRKY33, the central protein in the interactome, suggesting an important regulatory role for MKK9 as well. Furthermore, the interaction network also highlighted other known stress genes as part of the stress signaling core, including the mitochondrial ALTERNATIVE OXIDASE 1A (AOX1a), as well as 31 unannotated genes, hence uncovering their function (Supplementary Data 5).
a Protein-interaction network representing functional and physical interactions between proteins encoded by the shoot gene core. Node colors represent different clusters depending on neighbor interactions; connecting-line thickness represents the strength of the connections based on experimental data. Percentages indicate the number of proteins forming the cluster in relation to the total number of genes of the SVM gene core. b Cluster that contains the highest number of proteins of the four clusters. The name or AGI number of each gene is indicated next to the corresponding node. c–f GO terms related to biological processes obtained during the GO enrichment of each cluster. The colored square indicates which GO analysis corresponds with each cluster highlighted in (a). Statistical analysis was performed using gene set enrichment analysis (GSEA) and Benjamini–Hochberg (BH) p-value adjustment with default parameters. The parameter p.adjust represents the statistical significance of the group (p-value < 0.05 was used as cutoff). Count indicates the number of genes inside the category. GeneRatio reflects the percentage of DEGs in the complete GO category.
The red cluster contained proteins mainly related to maintenance of cell-wall integrity (GO:0042546, GO:0010411, GO:0006949) (Fig. 3c). The green cluster was marked by ‘alpha-amino acid metabolic process’ (GO:1901605) and ‘response to water deprivation’ (GO:0009414) as main GO terms, reflecting its role in both metabolism and responses to water availability (Fig. 3d). Lastly, the smallest cluster (yellow) contained proteins involved in hypoxia responses (GO:0071456), together with proteins related to other stress responses (biotic and wounding stress) (Fig. 3e). In conclusion, it is evident that shoot-stress signaling mainly mitigates alterations in water status and conserves cellular water homeostasis. In addition, MKK9 and WRKY transcription factors, specifically WRKY33, appear to be pivotal in the regulation of these responses.
A similar topology was obtained for the root-core interaction network (Fig. 4a). Of the four clusters, two contained the maximum number of proteins (representing 28% of the total number in the core) and, one of them (colored in blue), exhibited the highest number of connections within the network (Fig. 4b). As expected, the blue cluster showed the GO category that characterized the root core (GO:0071456: ‘cellular response to hypoxia’). The second-most relevant GO term (‘secondary metabolic process’, GO:0019748) covered genes involved in lignan biosynthesis, such as BCB, and phenylpropanoid biosynthesis (KISS ME DEADLY 1 and 4; KMD1/4). In addition, defense responses seemed to play an important role in this blue cluster (GO:0031347), as well as responses to external stimuli (GO:0009605) and to ethylene (GO:0009723), indicating a relevant role in environmental interactions (Fig. 4f). The MAPK protein MPK11 was found at a central position in the interaction network, possibly coordinating the activity of the remaining members of the blue cluster. Interestingly, the ethylene receptor ETHYLENE RESPONSE 2 (ETR2) and the downstream transcription factor ETHYLENE RESPONSE FACTOR 2 (ERF2), which is induced by ethylene34, were also present in this cluster, corroborating a pivotal role of this hormone in the generic stress response, at least in roots.
a Protein-interaction network representing functional and physical interactions between proteins forming the root gene core. Node colors represent different clusters depending on neighbor interactions, and connecting-line thickness represents the strength of the connections based on experimental data. Percentages indicate the number of proteins forming the cluster in relation to the total of genes of the SVM gene core. b Cluster that contains the highest number of proteins of the four clusters. The name or AGI number of each gene is indicated next to the corresponding node. c–f GO terms related to biological processes obtained during GO enrichment for each cluster. The colored square indicates which GO analysis corresponds with each cluster highlighted in (a). Statistical analysis was performed using gene set enrichment analysis (GSEA) and Benjamini-Hochberg (BH) p-value adjustment with default parameters. The parameter p.adjust represents the statistical significance of the group (p-value < 0.05 was used as cutoff). Count indicates the number of genes inside the category. GeneRatio reflects the percentage of DEGs in the complete GO category.
The green cluster (Fig. 4d) was characterized by GO terms related to oxidative-stress responses (GO:0006979), cellular metabolism of amino acids (GO:0009063) and vitamins (GO:0009110; GO:0006766), and transport of inorganic compounds (GO:0006829), revealing a potential role for the maintenance of shoot central metabolism and physiology in root responses. Finally, both red and yellow clusters showed a reduced number of proteins and interaction levels compared to the previous ones (Fig. 4c, e). The yellow cluster showed GO terms involved mainly in cell-wall homeostasis and modification (GO:0009826; GO:0016049; GO:0009828; GO:0006949; GO:0010025) while the red one included GO terms water responses (GO:0009415), and hypoxia (GO:0071456), among others.
Overall, we conclude that shoot stress responses are mostly related to the maintenance of water potential and homeostasis and, secondary, to the maintenance of normoxia levels; while in roots, the opposite trend is observed. In addition, growth regulation and metabolism as well as cell-wall homeostasis are important aspects of core stress signaling in both tissues. The complete list of the genes in the SVM gene cores classified in the four clusters (blue, yellow, red and green) is found in Supplementary Data file 5.
Role of ethylene in the SVM gene cores
Because of its key role in a multitude of stresses, and given the presence of ETR2 and ERF2 in the root core, as well as MKK9—known to play a pivotal role in the activation of MPK6 under ethylene signaling35—in the shoot core, the ethylene responsiveness of the genes within the SVM gene cores was investigated. To define a robust list of such genes, we combined the publicly available data of an ETHYLENE INSENSITIVE 3 (EIN3) ChIP-seq analysis15 with a set of DEGs under early (4 h, GSE14247) and late (24 h, GSE83573) ethylene treatment, forming an ethylene-responsiveness database (Supplementary Data file 6). More than 50% of the genes in the SVM gene core for shoots and roots were ethylene responsive (Fig. 5a), underpinning the relevance of ethylene in both gene cores. Remarkably, the number of ethylene-responsive genes increased to 77% and 62.7% in the blue cluster of shoots and roots cores, respectively, further substantiating the central role of ethylene in core stress signaling.
a Number and percentage of genes of each gene core present in the ethylene-responsiveness database. Number hits refers to the number of genes from the core present in the ethylene-responsiveness database (Supplementary Data file 6), while the percentage indicates the number of hits compared to the total number of genes present in the complete core. b Number and percentage of genes of each gene core with putative WRKY33 binding sites. Number hits refers to the number of genes from the core that present WRKY33 binding motifs, while the percentage indicates the number of hits compared to the total number of genes present in the complete core. The last column shows the overlap between the presence of WRKY33 binding sites and ethylene responsiveness. c, d STRING interaction network of the ethylene-response genes from the blue cluster of the shoot c and root d gene cores. The connection between different nodes indicates protein–protein association. Thickness of connecting lines represents the strength of the connections based on experimental data. e Transcriptional responses of key genes of the SVM gene cores under stress conditions in WT (Col-0) and the ethylene-insensitive mutant ein2-5 after 1 and 3 h of stress exposure. Fold-change expression is calculated relative to the expression levels of Col-0 control samples at both 1 and 3 h. The color scale indicates downregulation (0–1; red), no change (1; black), and upregulation (1–5; green). Source data are provided as a Source Data file. f Map of predicted WRKY33 binding-site motifs in ethylene-related genes.
The subgroup of genes in the blue clusters detected as ethylene-related genes were used to construct a protein-interaction network (hereafter defined as ethylene-related clusters). In the case of the shoots, ethylene-related proteins showed the same interconnected pattern as in the complete network (Figs. 3b and 5c). Moreover, WRKY33 still appeared as a central node in shoot stress signaling, together with MKK9. Notably, MKK9, along with MPK3 and MPK6, has been directly linked to both ethylene biosynthesis and signaling35,36. In addition, LYSINE HISTIDINE TRANSPORTER 1 (LHT1), an amino-acid importer responsible for the transport of the ethylene precursor 1-aminocyclopropane-1-carboxylate (ACC37); was also part of the ethylene-related shoot cluster, as well as the mitochondrial AOX1a, which connects the regulation of respiration to stress signaling in an ethylene-dependent manner38.
The central role of WRKY33 in the shoot gene core highlights its potential in regulating stress responses. To further explore this connection, we performed an in silico study of the presence of WRKY33 binding sites, identified by the binding motif TTGACY, which was empirically determined through ChIP-seq analysis39 (Supplementary Date file 7). Given the putative central role of WRKY33 in the gene core, it is unsurprising that 70% and 75% of genes forming the blue clusters in shoots and roots, respectively, contained WRKY33 binding motifs (Fig. 5b). Furthermore, when comparing the genes presenting this motif with the ethylene-responsive genes calculated previously, more than 54% of genes overlap between these two conditions, underscoring the close relation between ethylene signaling and the potential function of WRKY33. This connection is further supported by the analysis genes related to both ethylene biosynthesis and signaling (Fig. 5f). More than 50% of the biosynthetic genes, including SAMS4, several ACSs (ACS2, ACS5–8, and ACS11), and all ACOs (ACO1–5), as well as key genes involved in ethylene signaling, such as the receptors ETR2 and ERS2, CTR1, EIN3–EIL1, and EBF2, along with one of the primary ethylene transcription factors ERF1, are also targeted by WRKY33. This further reinforces the relationship between WRKY33 and ethylene responses.
In roots, the number of genes constituting the interaction network was reduced (Figs. 4b and 5d). Consequently, the number of interactions was also decreased. Nevertheless, the ETR2–ERF2 module again appeared at its center, coordinating other nodes of the network. In addition, the core network also contained AUXIN-REGULATED GENE INVOLVED IN ORGAN SIZE (ARGOS), which is part of a negative feedback mechanism to attenuate ethylene responses, further highlighting the importance of coordinated ethylene signaling40. Several transcription factors that have been experimentally linked to specific stresses or processes, including RAP2.6L/ERF113 (wounding), NAC047 (partial submergence), and NAC6 (leaf senescence), also appeared in this cluster, suggesting a more general function for all.
Lastly, to further demonstrate the role of ethylene in the regulation of both SVM gene cores, we studied the transcriptional responses of key genes from the stress gene cores under stress conditions in wild-type Col-0 and in the ethylene-insensitive ein2-5 mutant (Fig. 5e). We observed that most of the studied genes exhibited clear down-regulation in the ein2-5 mutant compared to Col-0 under control conditions, empirically confirming that ethylene signaling is required for the transcriptional activation of these genes, particularly at early time points (1 h). The exception was WRKY33, which remained partially unaltered, especially under cold and heat conditions, indicating functioning upstream ethylene signaling during the regulation of stress responses.
The high degree of interconnection within the blue clusters of both SVM gene cores, coupled with the presence of numerous genes related to ethylene responses, underscores the central role of ethylene in regulating stress responses in both shoots and roots. This is further supported by the transcriptional data from the ethylene-insensitive ein2-5 mutant, which shows a general down-regulation of genes in both gene cores, particularly at early time points, confirming the requirement of ethylene signaling for the activation of these stress-related genes.
A central role for EXPANSINS, AP2/ERFs, WRKYs, and MAPKs in the SVM gene cores
Certain gene families were identified as crucial players in the SVM gene cores, such as EXPs, AP2/ERFs, WRKYs and MAPKs. In addition, novel gene families were also identified by the SVM Clustering algorithm, including USPs. To provide a complete and detailed map of the function of these gene families in stress responses, we investigated the transcriptional alterations of all their family members under all conditions in our meta-analysis (Fig. 6; Supplementary Data file 8). In addition, as a second validation, experimentally validated data about specific members of the selected families corroborated our results (Supplementary Data file 9). Select families are covered in the next section; with further details are presented in the Supplementary Information file.
Shown are members of the expansin (EXP) (a), ethylene-response factors (ERF) (b), WRKYs (c), mitogen-activated protein kinases (MAPK) (d) and universal stress proteins (USP) (e) families across the conditions analyzed. Green and red represent up- or down-regulation in each condition. Bicolored cells indicate that the behavior of the gene is variable depending on early or late responses (for complete information on all members of the family, see Supplementary File 7).
EXPANSIN (EXP) superfamily
EXPs were detected as the main enriched gene family in the root core (Fig. S4). They enable cell expansion and increase cell-wall flexibility41,42. The EXP superfamily is divided into four groups: EXPA, EXPB, expansins-like A (EXPLA) and EXPLB. In our DEG database, both EXPA and EXPB subfamilies are predominantly down-regulated in several stress conditions, with tissue-specific patterns (Fig. 6a; Supplementary Data file 8). Three members of EXPA (EXPA1, EXPA8 and EXPA15) and one EXPB (EXPB3) were part of SVM stress cores. While transcriptional alterations of EXPA8, EXPA15 and EXPB3 were observed for roots under certain conditions, EXPA1 showed transcriptional alterations in all tissues and timepoints in 7 out of 10 studied stresses. These findings highlight the importance of the tissue specificity of EXPs in stress signaling and of EXPA1 as a main stress regulator. Interestingly, though most EXPs were down-regulated in response to stress, the EXPLA and EXPLB groups (with two genes present in the SVM gene cores, EXPLA1 and EXPLB1) were up-regulated under several stress conditions, suggesting a potential positive role for these subfamilies in stress responses.
Ethylene response factor (ERF) family
The ERF family of transcription factors is considered to be crucial in both growth and defense responses43. ERFs are part of the AP2/ERF superfamily, comprising 147 member divided in three sub-families: 18 AP2s, 122 ERFs and six RELATED TO ABSCISIC ACID INSENSITIVE 3/VIVIPAROUS 1 (RAVs), as well as a not-yet-classified gene (AT4G13040)44. Whereas few to no transcriptional alterations were found for the AP2 family, the RAV and ERF subfamilies revealed to be highly affected by stress, often with very distinct expression patterns (Fig. 6b; Supplementary Data file 8).
Almost all ERF subgroups displayed substantial transcriptional changes under the studied stress conditions (Supplementary Data file 8). Seven ERFs were classified as members of the SVM gene core, six in the root core (ERF2, TINY (ERF040), DREB2A (ERF045), DEWAX (ERF107), RAP2.6 (ERF108), and RAP2.6L (ERF113)), and only one (RAP2.4D (ERF058)) in the shoot core (Fig. 6b). Some of these ERFs showed tissue specificity. For instance, RAP2.6/ERF108 and RAP2.6L/ERF113 showed a similar transcriptional pattern in root tissue, yet distinct patterns in shoots. While RAP2.6/ERF108 was down-regulated under the early exposure to partial submergence, RAP2.6L/ERF113 was up-regulated under both types of submergence. In contrast, ERF2 and DREB2A/ERF045 were broadly expressed in both roots and shoots in most conditions. Since the direct relationship between the members of the ERF family and ethylene is not always clear, we investigated their responsiveness to ACC as well as to ethylene (Supplementary Data files 6 and 8). ERF2, RAP2.6/ERF108 and RAP2.6L/ERF113 were up-regulated by ACC, as well as present in the ethylene-responsiveness dataset. The other ERFs were either only up-regulated by ACC (TINY/ERF040 and DEWAX/ERF107), only ethylene responsive (DREB2A/ERF045), or not detected under ACC, nor ethylene treatment (RAP2.4D/ERF058).
Mitogen-activated protein kinase (MAPK) superfamily
MAPKs are important signaling proteins in many intracellular responses to developmental, physiological and/or environmental stimuli45. MAPK signaling cascades are typically characterized by a sequence of phosphorylation and activation events along three levels comprising members of mitogen-activated protein kinase kinase kinases (MAPKKK, MKKK or MEKK), mitogen-activated protein kinase kinases (MAPKK or MKK), and mitogen-activated protein kinases (MAPK or MPK). In A. thaliana 69 MAPKKKs, 10 MAPKKs and 20 MAPKs have been described45.
Four out of the ten members of the MAPKK group were unresponsive to any of the conditions studied (Supplementary Data file 8). However, of the stress-responsive MAPKKs, MKK9 stood out, being up-regulated by six different conditions, mainly in shoot tissue (Fig. 6d). High salinity and osmotic stress elevated MKK9 transcription in roots. Not surprisingly, MKK9 appeared as part of the shoot core (Fig. 3b), indicative for its central regulatory role in abiotic stress responses.
Out of 20 MAPKs, only six did not show a transcriptional effect under any of the stress conditions. Conversely, MPK11, MPK3, MPK5 and MPK19 stood out given transcriptional alteration under seven, six, five and four different stress conditions, respectively (Supplementary Data file 8). However, only MPK11 was retained by SVM Clustering as a stress core gene (Fig. 6d). Cold, osmotic and UV stresses up-regulated MPK11 expression in all tissues, while salt induction was root-specific. Wounding (shoots and roots), drought (roots) and heat (roots) induced MPK11 transcription predominantly at early timepoints, implying its relevance specifically during the initial stages of the stress response.
Universal stress protein (USP) superfamily
USPs are proteins involved in, as their name suggests, a broad range of metabolic processes related to stress, such as nutrient starvation, heat shock and oxidative stress46. Nevertheless, their specific roles and molecular mechanisms remain largely unknown. In A. thaliana, 41 genes encode for USP proteins, cataloged according to domain organization. From our analysis, the USP family and the single gene belonging to the double USP domain group (USPUSP) appear to be involved in all stress conditions (Supplementary Data file 8).
USP12, USP25, and USPUSP1 were detected as part of the root gene core. USP12, characterized as a gene involved in ROS modulation in anoxia conditions, is up-regulated in submergence conditions, but also in heat and osmotic conditions in both tissues (Fig. 6e; Supplementary Data file 9). USP25 and USPUSP1 appear to be tissue-specific players in the general stress response, being only expressed in roots. It is evident that the function of USPs deserves more scrutiny, given their highly specific expression patterns as well as the central role of certain family members in general stress signaling.
Experimental validation: the role of EXPAs, WRKY33, MKK9, and LHT1 in general stress responses
As part of a third validation supporting the role of members of the gene cores as central stress regulators, we first studied the transcriptional alterations of three members of the EXPA family (EXPA1, EXPA10, and EXPA14), represented in the root core, and EXPA1 as a part of both gene cores, using transgenic translational-reporter lines (pEXPA1::EXPA1–mCherry (Fig. 7a), pEXPA10::EXPA10–mCherry (Fig. 7b) and pEXPA14::EXPA14–mCherry (Fig. 7c)41. We exposed the three lines to cold, salt, and osmotic stress for 1 h. The expression patterns in control conditions corresponded with the patterns described by Samalova et al. (2023)41. Short-term exposure to stress saw changes in mCherry intensity levels, indicating a change in EXPA abundance in the studied tissue. Specifically, salt treatment dramatically increased the levels of EXPA1–mCherry and EXPA10–mCherry, while decreasing the level of EXPA14–mCherry (Fig. 7d), corroborating the expression changes after the equivalent treatment obtained in our meta-analysis (Fig. 7e). Osmotic treatment modestly increased the levels of EXPA1–mCherry and EXPA10–mCherry, without affecting the level of EXPA14–mCherry. Upon cold treatment, none of the levels of the studied EXPAs differed from those of the control samples, again confirming the expression data. In conclusion, data from the translational reporter lines matched with the transcriptional data derived from the meta-analysis, experimentally validating the robustness of our analysis, and positioning EXPAs as key regulators of multiple stress responses as evidenced by their relevance in both shoot and root gene cores.
a–c mCherry fluorescence intensity for EXPApro:EXPA–mCherry reporter lines for EXPA1 (a), EXPA10 (b), and EXPA14 (c) in the lateral root cap for EXPA1, in the root tip and a region of the cortex for EXPA10 and in a region of the cortex for EXPA14. The measured region of interest (ROI) is framed with a yellow box. Scale = 50 µm. d Fluorescence-intensity changes in relation to control conditions (n for EXPA1 = 34 (control), 17 (cold), 15 (salt), 20 (osmotic); for EXPA10 = 38 (control), 20 (cold), 20 (salt), 19 (osmotic); for EXPA14 = 28 (control), 13 (cold), 16 (salt), 12 (osmotic)). Data are presented as mean values ± SD. A p-value of p < 0.05 was assumed as statistically significant (*<0.05; **<0.01; ***<0.001). e Heatmap showing intensity fold-changes represented as relative marker-gene intensity compared to control conditions (ranging from 0 (orange) to 2 (blue)) and expression fold changes based on transcriptomic analyses (represented as the log2-fold change compared to control conditions, ranging from –1.5 (red) to 1.5 (green)) for each stress condition. Data used for constructing the heatmaps are provided in the Data Source File. Shapiro’s test and Levene’s test were used to assess normality and variance of each dataset. To compare the intensity of the different treatments with the control samples, Student’s T-test and Wilcoxon rank sum exact test were used for parametric and non-parametric testing, respectively. Complete statistical analysis is available in Supplementary Data file 15.
Secondly, we investigated the function of both WRKY33 and MKK9, putative key regulators in the shoot core (Fig. 3b). We assessed stress tolerance by comparing alterations in rosette growth between the loss-of-function mutants wrky33-2 and mkk9-1 and the wild-type Col-0 (Fig. 8a–d). To consider a representative set of stress conditions, we selected five conditions studied in our analysis, covering each of the clusters obtained by HCA: complete submergence, cold, heat, salt, and wounding (Fig. 1).
Rosette area of Col-0 and wrky33-2 (a and b), mkk9-1 (c and d) lht1-5 (e and f) plantlets. Rosette area was measured after 4 days of cold (4 °C), complete submergence, heat (42 °C), salt stress (100 mM NaCl), and wounding, and after 5 days in recovery conditions for wrky33-2 plants and 10 days for mkk9-1 and lht1-5 (n > 20 plants per sample; specific n values are specified in Supplementary Data file 15). Scale = 0.5 cm. Lighting conditions differed between the sets, and exposure of panel (e) was increased equally to 1.5 to improve the quality of the images. Source data are provided as a Source Data file. A p-value < 0.05 was assumed as statistically significant and different letters represent statistical differences between samples (specific p-values are specified in Supplementary Data file 15). Shapiro’s test and Levene’s test were used to assess normality and variance of each dataset. To compare rosette areas, a two-way ANOVA followed by post-hoc Tukey’s test were used as parametric tests, while a Kruskal–Wallis rank sum test followed by Dunn’s Multiple Comparison test were used as non-parametric tests. Complete statistical analysis is available in Supplementary Data file 15.
The wrky33-2 mutant exhibited notable phenotypic differences when exposed to various stress conditions (Fig. 8a and b). Under cold, heat, and wounding conditions, wrky33-2 was hypersensitive, evidenced by reductions in both rosette area (Fig. 8a and b) and relative biomass compared to control conditions (Fig. S5a). In contrast, under complete submergence and salt stress, wrky33-2 mutants responded differently. While Col-0 plants experienced a marked decrease in both rosette size and relative biomass, wrky33-2 mutants displayed either non-significant changes or even improvements in these parameters. On the other hand, mkk9–1 mutants behaved differently from the wild type in all tested conditions and displayed an increased resistance to the vast majority of stress conditions. mkk9-1 mutants had increased tolerance under cold, heat, and wounding stress, with rosettes statistically significantly larger than those of treated Col-0 plants under the same conditions (Fig. 8c and d). However, under salt stress, mkk9-1 exhibited hypersensitivity, showing a more pronounced reduction in rosette size compared to Col-0. Lastly, under complete submergence, mkk9-1 rosettes showed enhanced growth, being larger than both treated and untreated Col-0 plants. These observations were further supported by the relative biomass measurements, which followed the same trends as rosette size (Fig. S5b).
Lastly, we analyzed the response of the lht1-5 loss-of-function mutant to the different stress treatments (Fig. 8e and f). Both its presence as a member of the shoot core, as well as its function in amino acid and ACC transport, suggest that LHT1 could act as a vital component of general stress signaling. Indeed, plants that lack functional LHT1 are hypersensitive to cold, complete submergence, heat, and salt stress. In contrast, lht1-5 plants were less affected than Col-0 plants upon wounding. Similar to the above-mentioned results, these findings were corroborated by corresponding changes in the relative biomass compared to the control condition measure (Fig. S5c).
To gain deeper insights into the roles of WRKY33 and MKK9 in the response to the tested stress conditions and the connection with ethylene, its production was analyzed under the same treatments and compared to wild-type (Col-0) (Fig. S6). Ethylene emanation increased in response to all tested abiotic stresses except cold stress, wherein a decrease was observed compared to control conditions. In WRKY33-deficient plants, ethylene levels were statistically significantly higher under control conditions and in response to heat, wounding, and complete-submergence stress, whereas slight reductions were observed under cold and salt stress. Similarly, mkk9-1 mutants displayed hypersensitivity to heat, complete submergence, and cold stress, as reflected by altered ethylene levels, although these did not reach the higher levels produced by wrky33-2. These findings highlight the regulatory influence of WRKY33 and MKK9 on ethylene production and underscore their critical roles in stress signaling mechanisms, intertwined with the ethylene pathway.
Further investigation of the transcriptional alterations of key genes from the SVM gene cores (WRKY33, MKK9, BCB, AOX1a, AT1G55450, JAZ1, OPR3, and TCH3) and ethylene-related genes (ETR2, ERF1, and EBF2) was conducted by real-time quantitative PCR analysis under the same set of stresses (Fig. S7). In the wrky33-2 mutant, a marked down-regulation of most assayed genes was seen under stress conditions, wherein the mutant exhibited hypersensitivity, particularly cold and wounding. However, in complete submergence, whereby wrky33-2 mutants performed better in terms of rosette size and relative biomass, some genes, such as BCB, ERF1, and ETR2, were up-regulated compared to Col-0. Conversely, transcriptional changes in the mkk9-1 mutant showed the opposite trend. Under stress conditions when the mutant displayed resistance (cold, complete submergence, heat, and wounding), a substantial number of genes were up-regulated relative to Col-0. In contrast, under salt stress, wherein the mkk9-1 mutant exhibited hypersensitivity, most assayed genes were statistically significantly down-regulated. These results further support contrasting regulatory roles for WRKY33 and MKK9 in stress-specific transcriptional responses.
The results of these analyses emphasize the overarching role of ethylene in regulating plant stress responses, as well as the pivotal contributions of WRKY33 and MKK9 in modulating the transcriptional activity of the stress gene core across various conditions. These findings not only validate the biological significance of the SVM gene core but also confirm its robustness, reproducibility, and suitability as a foundational framework for understanding general stress regulatory networks in plants.
