Meta-Analysis as a Tool to Identify Candidate Genes Involved in the Fagus sylvatica L. Abiotic Stress Response

: In this study, we aimed to evaluate whether candidate genes for abiotic stresses in Fagus sylvatica L. are also candidate genes for herbaceous plants, with the purpose of better deﬁning the abiotic stress response model of F. sylvatica . Therefore, a meta-analysis was performed on published papers related to abiotic stress. Firstly, we carried out a systematic review regarding the activity of 24 candidate genes selected for F. sylvatica under abiotic stress reported in 503 articles. After choosing the inclusion criteria, 73 articles out of 503, regarding 12 candidate genes, were included in this analysis. We performed an exploratory meta-analysis based on the random-effect model and the combined effect-size approach (Cohen’s d). The results obtained through Forest and Funnel plots indicate that the candidate genes for F. sylvatica are considered to be candidate genes in other herbaceous species. These results allowed us to set up models of plants’ response to abiotic stresses implementing the stress models in forest species. The results of this study will serve to bridge knowledge gaps regarding the pathways of response to abiotic stresses in trees based on the meta-analysis. The study approach used could be extended to observe larger gene databases and different species. GTLs, GT-element-binding proteins transcription factor [113]; NIG1, basic helix-loop-helix-type [114].


Introduction
Plants, during their life cycle, are exposed to continuous chemical-physical and biological changes in their habitat or in an environment altered by humans. These changes are often unfavourable or stressful for growth and development. Therefore, plants have evolved various dynamic molecular reprogramming events to rapidly perceive changes and adapt accordingly [1][2][3]. In recent decades, significant progress has been achieved in understanding the physiological, cellular, and molecular mechanisms important for plant adaptation to certain climatic conditions [1][2][3][4][5]. Using molecular biology approaches, many abiotic stress-inducible genes have been identified, and their functions have been characterised in plants [6,7]. This has permitted finding evidence of adaptive loci at the DNA level using the candidate gene approach [8]. The development of Next-Generation Sequencing technologies has allowed the sequencing of thousands of genotypes that will give a significant advancement in the identification of genes important to adaptation [6,[8][9][10][11][12][13][14]. This large amount of available data has improved all omics disciplines with respect to the study of regulatory pathways and gene expression and their correlation to abiotic stress. Despite the variety of available data and rich resources, it is important to note that the diversity of experimental designs and types of analyses can be a limiting factor for comparing these genomic data and the development of stress models. Due to the variability present in the methods applied to induce stress in the sampled tissues and the bioinformatics tool used for data generation and analysis, this large amount of information has given rise to a new source of heterogeneity [11][12][13]. Ambiguous information from distinct resources may result in methodological artefacts derived from observational data. In addition, the heterogeneity of results present among different plant species demonstrates that efforts are necessary to transfer the assessment of gene roles from reference models to other species [14]. The challenge is to find ways to combine these resources to build reliable datasets of information that are easily manageable and testable. Systemic approaches could bring many advantages and pave the way toward comprehensive modelling. To obtain a comprehensive view of plant responses to environmental stresses, it will be important to integrate omics data with systemic biology data and to develop computational models [11,14,15]. In their review, Cramer et al. [14] explored the perspectives of systemic approaches in determining molecular responses to abiotic stresses. In the post-genomic era, comprehensive analyses using systemic or omics approaches have increased our understanding of the complex molecular regulatory networks associated with stress adaptation and tolerance. Although many of the functions of individual parts are unknown, their function can sometimes be inferred through association with other known parts, providing a better understanding of the biological system as a whole [14]. New models can be formed from the large amount of data collected and can lead to new hypotheses generated by these [15]. The most used models to describe signalling or metabolic pathways are based on theoretical models [15][16][17]. These models are based on the laws of physics and use differential or algebraic equations to represent biological processes. Therefore, these models cannot represent dynamic processes such as gene regulation. Over the years, alternative methods based on computational models have emerged, and statistical and systemic models are the dominant models in stress-response studies [18,19]. Due to the potential for hypothesis testing, prediction, and uncertainty quantification, statistical models have proven to be very useful in these studies. However, their main limitation is that they fail to incorporate the diversity present in gene expression data acquired under different conditions and experimental conditions [15]. Models based on mixed types overcome these limitations by setting fixed-and random-effects terms. These models are useful for combining omics data from highly variable studies, and meta-analysis is one of the main types of analysis that allows this combination. Meta-analysis is a statistical procedure for analysing data combined from several studies and can be an important source of concise and updated information. This type of analysis has been widely used to observe candidate genes for abiotic stress response [20][21][22][23][24]. The meta-analysis approach makes it possible to overcome the lack of characterising candidate gene datasets or differentially expressed genes, mainly caused by the heterogeneity of single studies. In their review, Tseng et al. [25] identify four principal meta-analysis approaches: combined p-value, combined effect size, combined ranks, and directly merging raw data [25]. The objective of this work is to provide further insight into the abiotic stress response and to increase the accuracy of response models. For this, we conducted a statistical analysis on some candidate gene data and related stress. We used an exploratory meta-analysis approach as a tool for the evaluation of a dataset of candidate genes regarding the abiotic stress response in Fagus sylvatica L. More specifically, we conducted a systematic quantitative review, through the method of combined effect size (ES), of a wide range of studies to increase the knowledge on the roles of candidate genes included in the initial dataset, through the acquisition of information present in the literature. The aim of this work is to improve existing models that describe the response to environmental stimuli of forest species, characterised by long life cycles. To achieve this goal, we first performed a systematic review to obtain a large database of articles concerning the study of candidate genes involved in responding to different abiotic stresses in herbaceous plants. Subsequently, we conducted an exploratory meta-analysis that provided a quantitative summary of the results.

Materials and Methods
The experimental design followed to perform the meta-analysis is shown in Figure 1. The details of each step are reported in successive paragraphs.

Candidate Genes Database
The first step was to search for and collect the scientific articles concerning the expression, sequencing, and annotation of specific genes in response to environmental stimuli in F. sylvatica, and to place them into a single dataset (Table S1). A database of 24 selected genes involved in response to abiotic stresses (drought, salt, cold and heat stress) in F. sylvatica was implemented.

Systematic Review
This step is particularly important as it allows us to decide which article can be considered as following the inclusion criteria (decision on inclusion criteria) and obtain the final papers' database (included papers database) to be used for the successive meta-analysis.
As a first step, we collected information regarding enzymes encoded by the F. sylvatica selected genes. In particular, we searched information on: (i) enzyme name, (ii) transcript name referred to Arabidopsis thaliana genome annotation TAIR, (iii) family, (iv) subfamily, (v) molecular function, (vi) biological process, (vii) cellular component, and (viii) protein class. Such information was collected through a search on major databases (KEGG, Panther, UniProt, TAIR) and collected into a table (Table S2). Successively, we conducted a systematic review regarding each gene present in the initial dataset following the guidelines proposed by the PRISMA protocol [26]. This procedure provides a standardised framework for meta-analysis and systematic reviews, thus allowing the reliability and replicability of the results obtained. The PRISMA protocol is based on identifying articles through exhaustive bibliographical research, screening of collected articles through exclusion of duplicates, title-abstract screening, evaluation of articles for eligibility, and the inclusion of studies of interest in quantitative synthesis (meta-analysis) [26,27].
We searched the available literature for each selected gene through the major databases (Scopus, Research Gate, and Spring Link) and search engines (PubMed and Google Scholar) using the following search terms: "gene ID" + "plant *" + "candidate gene *" + "abiotic stress" + "response *" + "SNPs". The terms have been used both together and in combination of two or three terms (such as 'Gene ID' + 'abiotic stress' or 'Gene ID' + 'abiotic stress' + 'response', etc.).
Based on the guidance of Nakagawa et al. [28], to overcome data independence, publication bias, and outlier issues, we limited our search to English peer-reviewed papers and excluded reviews (which could have led to duplicate information from some papers). After initial screening, we determined the inclusion criteria necessary to collect usable data: compare control to treated; quantitative PCR articles with the presence of comparable gene expression data (|log 2 fold change| ≥ 0.00012); provide the mean, a measure of dispersion (Standard Error-SE or Standard Deviation-SD), and sample size, or provide the original data to calculate mean and SD. Moreover, the minimum number of four articles available for each abiotic stress was considered necessary for meta-analysis. In the case of sub-group analysis, we needed at least four articles for each gene for each stress. Because some articles presented more than one case study, we considered those reporting results from different treatments, species, or both as independent cases to avoid cases of pseudo-replication.

Effect Size Calculation
After applying the inclusion criteria, the effect size was calculated. Effect size is a way of quantifying the size of the difference between two groups. It is particularly valuable for quantifying the effectiveness of a particular treatment relative to a control. By emphasizing the most important aspect of a treatment rather than its statistical significance, effect size promotes a more scientific approach to knowledge accumulation. Based on some other studies [29][30][31][32], we used Cohen's d as the effect size. Cohen's d describes the standardised mean difference in an effect and is useful to compare effects across studies when the dependent variables are measured differently. It is defined as the mean difference between two groups, divided by the standard deviation of both groups. Conventionally, in metaanalysis, the two groups are considered the experimental and control groups. Cohen [33] is given by Equation (1): where: The numerator is the difference between the log 2 fold change of the two groups of observations. The denominator is the pooled SD. The pooled estimate represents the mean value of the standard deviations of the treatment and the control groups. The articles collected during the systematic review were characterised by the presence of continuous outcomes. The main problem detected during the calculation of the effect size was the presence of incomplete articles of SD or other measures of dispersion. As this lack of information was present in most of the articles collected, we followed the lists multiple solutions of varying nature method reported by Weir et al. [30] to impute dispersion measures from incomplete reports. Among the listed solutions, we observed and applied the prognostic method, also described by the article by Ma et al. [29]. Therefore, the calculation of effect sizes can be divided into two steps: (i) calculation of assigned SE and the assigned SD through prognostic method, and (ii) calculation of SD pooled and then Cohen's d. To apply the prognostic method for SD imputation, we collected data from the study by Van Zhong et al. [34]. From this article, 27 SE values were selected. Through application of the prognostic method, it was possible to impute SD and calculate Cohen's d (Equation (1)).

Meta-Analysis: Random Effect Model, Forest and Funnel Plot
After conducting the systematic review of the collected papers and calculating the effect size, we performed the data analysis using the "Included papers database" (Table S4) and conducting an exploratory meta-analysis. For this analysis, we chose to follow the combined effect size approach [25], which assumes that standardised effect sizes can be combined across studies. According to Tseng et al. [25], the two most used models in this category are the fixed effects model (FEM) and the random effects model (REM). We decided to use REM given that: (i) the data analysed were continuous results and the REM is the most effective; (ii) data vary over time, the use of this model is recommended; and (iii) the REM is used when heterogeneous data characterised by a large variability of experimental conditions must be compared, without committing methodological errors. The analyses were conducted with the JASP 0.14.1 Meta-Analysis module [35]. JASP software supports a wide range of techniques commonly used for meta-analysis. These include fixed and random effects analysis, fixed and mixed effects meta-regression, Forest and Funnel plots, tests for funnel plot asymmetry, and more. The engine behind this power of analysis is the metafor R package [36]. JASP is user-friendly and freely available (https://jasp-stats.org/, last accessed on 18 December 2021) [35].

Interactive Analyses of Stress Meta-Analyses Data
To better understand the results obtained from the meta-analysis and the information present in the bibliography, we investigated the presence of possible relationships between the observed genes and the studied stresses. To do so, we illustrate simple set relationships between abiotic stresses, using the genes relevant to those stresses as commonalities. We then visualised such data through a Venn diagram, obtained using the R Venndiagram package [37]. The Venn diagram compares lists of genes from a set of experiments and identifies the genes shared between the experiments or unique to an experiment in relation to a stress.

Candidate Gene Dataset, Systematic Review, and Inclusion Criteria
From the bibliographic review carried out to collect information regarding genes involved in response to abiotic stress in Fagus sylvatica L., we obtained information regarding 24 selected genes (Table S1). Based on the results obtained, we carried out a systematic review as previously described. We collected a total of 503 articles related to the expression study of 24 selected genes in other species (Table S3). We applied the PRISMA protocol parameters for screening the collected articles, and Figure 2 shows the complete workflow of study selection and screening of eligible datasets.
Following the guidelines suggested by Nakagawa et al. [28], we limited the systematic review to English peer-reviewed papers and excluded reviews. Applying these limitations allowed us to avoid the collection of duplicate studies. During the screening phase, we had to exclude 94 articles. The remaining 409 articles were submitted to the eligibility phase of the PRISMA protocol. This involved the selection and application of the inclusion criteria, which considers both: (i) whether or not the comparison between treated and control is available, and (ii) the presence of measures of dispersion, sample size, the presence of original data, or a combination of the three to calculate the SD. We excluded 336 articles that did not meet the chosen inclusion criteria. From the remaining 73 articles, data useful for conducting the meta-analysis were derived, as described in Materials and Methods. The reduction we had to apply due to the lack of usable data to conduct the meta-analysis resulted in the exclusion of 12 selected genes represented in Table 1, which also shows the number of initial and included studies and the number of case studies analysed. For each selected gene, a database was made with the information and data collected for the: (i) selected genes, (ii) reference, (iii) species (iv), type of abiotic stress, (v) log 2 fold change of treated plant and control plant, (vi) number of treated plant and control plant, (vii) imputed SE of treated plant and control plant, (viii) imputed SD of treated plant and control plant, (ix) S pooled , and (x) effect sizes (Table S4).  [26]. Our results are shown within the brackets. Table 1. Summary of the 12 selected genes for which the number of articles included in the metaanalysis is reported. For each gene, the initial articles, included articles, and case studies are reported. Study cases represent selected cases within the same article that differ in exposure time, treatment, plant species, and type of stress.

Meta-Analysis Results
We performed the meta-analysis through the software JASP, which allowed us to obtain the cumulative effect size values for each selected gene. These values, including 95% confidence intervals, are reported in Table 2, and the graphical representation is shown through Forest and Funnel plots ( Figure S1). Through observation of the I 2 index (Table 2), we observed that for Caffeic acid o-methyltransferase, Cytosolic class I small heat shock protein, Formate dehydrogenase, Glyceraldehyde-3-phosphate dehydrogenase, Heat shock protein 70, Light harvesting complex II protein, S-adenosylmethionine synthase, and Xyloglucan endotransglusylase hydrolase protein 23 genes, the heterogeneity was high because I 2 is greater than 75% [33]. To investigate what caused the heterogeneity, we decided to carry out subgroup analysis to assess whether, by conducting a meta-analysis by single stress (drought, heat, cold, salinity), the heterogeneity decreased. As shown in Table 3, we performed subgroup analyses by subdividing the analysed studies according to the type of stress observed. We could carry out this type of analysis only for Cytosolic class I small heat shock protein, S-adenosylmethionine synthase, and Heat shock protein 70 genes because the number of articles available for each abiotic stress was sufficient (at least four articles for each gene for each stress). Table 2. Summary of results obtained through applying the random-effects model (REM). For each candidate gene observed, we report: the name of the gene, the species studied, the stress studied, the cumulative effect size values, the 95% confidence interval values, and the I 2 index. The I 2 index quantifies the degree of heterogeneity in analysis and measures the extent of true heterogeneity by dividing the difference between Q (chi-square statistic) and its degrees of freedom (k − 1) by Q, and multiplying the result by 100. After subgrouping, the heterogeneity not caused by chance of the observations for S-adenosylmethionine synthase does not seem to decrease. The heterogeneity could be derived from the variability of the experimental conditions under which the studied plant species are subjected to the stresses. However, a decrease in heterogeneity can be detected for only those studies involving the observation of cold stress response, an indicator that the experimental conditions were similar. We observed similar results for Heat shock protein 70, where only the subgroup related to heat stresses appears to decrease in heterogeneity. Subgrouping for Cytosolic class I small heat-shock protein resulted in decreased heterogeneity values for all observed stresses.

Candidate Genes
The main result of the meta-analysis is represented by a graph, called a forest plot, which depicts the extension of the effect size for each study case and the cumulative effect size obtained with the random effects model (REM). Figures 3 and 4 show the forest plots for the subgroup analysis of Cytosolic class I small heat shock protein and S-adenosylmethionine decarboxylase. Table 3. Summary of heterogeneity indexes for only those genes to which subgrouping could be applied. For each candidate gene, we report: the name of the gene, I 2 cumulative, the stress studied for subgroup, the cumulative effect size values for each subgroup, the 95% confidence interval values for each subgroup, and I 2 for each subgroup. The I 2 index represents the degrees of variability: it quantifies the degree of heterogeneity in analysis and measures the extent of true heterogeneity by dividing the difference between Q (chi-square statistic) and its degrees of freedom (k − 1) by Q and multiplying the result by 100.    Table S4.

Candidate
Observing the cumulative effect size values, shown in Table 3 and Figure 3, it can be inferred that the expression of the gene encoding for cytosolic class I small heat-shock protein is positively regulated in response to salinity, cold, and drought stress. The large heterogeneity of studies for the response to cold stress confers a lower statistical potential, observable by the model confidence interval intersecting zero (Table 3, Figure 3). Figure 4 displays the forest plots of the subgroup analysis for the gene encoding for S-adenosylmethionine synthase. Observing the cumulative effect size reported in Figure 4 and Table 3, it is possible to observe that the expression of the gene encoding for S-adenosylmethionine synthase is positively regulated during the response to cold stress and negatively regulated during the response to salinity stress. Again, the heterogeneity of the studies confers greater uncertainty. The confidence interval, in both cases, intersects 0. Thus, although a subgroup analysis can be performed, the heterogeneity of these studies still confers a degree of uncertainty. The situation is somewhat different for analyses performed in the absence of subgroup division of study cases. The variance of the studies in this type of analysis is higher than in subgroups. Despite this, it was possible to observe values of the cumulative effect size of the REM having a confidence interval that does not intersect the 0. In Figure 5, forest plots with the relative effect size of some genes analysed without the subdivision into subgroups are reported. It can be seen that the genes encoding for S-adenosyl decarboxylase and xyloglucan endotransglucosylase hydrolase protein 23 exhibit positive expression regulation during the response to abiotic stresses.

Comparison between Selected Genes and Related Stress
To assess the relationship between the genes relevant to stress response in the studied stresses, we analysed the data obtained from the meta-analysis through a Venn diagram ( Figure 6). The graph represents the genes that have been observed to be relevant to the response to one or more stresses. As can be observed, all the genes analysed are involved in the response to multiple stresses, and none of them shows a stress-specific relationship. It is interesting to note that there is the presence of only three genes (GAPDH, SAM, and LHC-II) that are relevant for the response to all the stresses observed (drought, cold, heat, and salt stress). Figure 6. Venn diagram regarding stresses and candidate genes whose expression is regulated in response to these stresses. The results were generated using Venndiagram R package [37]. The data present within and at the intersections of the sets represent the ID of the observed genes. The sets represent the type of stress studied.

Discussion
The information collected and analysed by an exploratory meta-analysis based on the random-effect model and the combined effect-size approach (Cohen's d) allowed us to identify and select 12 genes responding to abiotic stresses in F. sylvatica. As can be seen from our overall synthesis of the collected studies revealed heterogeneity between studies (Table 2, Figure S1, Figures 3 and 4), for Heat shock protein 70 ( Figure S1), cytosolic class I small heat-shock protein (Figure 3), and S-adenosylmethionine synthase (Figure 4), it was possible to conduct an analysis in subgroups, limited only to case studies concerning a specific stress (Table 3). This type of analysis, in some cases, has led to a reduction in the heterogeneity present originally (Table 3, Figure S1, Figures 3 and 4). We hypothesised that this heterogeneity might also be due to different conditions and the time of exposure to which plants are subjected to abiotic stress. Thus, we speculated that the heterogeneity depends on methodological differences between studies included in the meta-analysis. Despite this great heterogeneity, it was possible to observe the cumulative effect of the various case studies for the 12 genes analysed ( Table 2). This effect, obtained through the application of the REM, allowed us to observe the expression of genes in response to stresses.
The results obtained allowed us to set up models of plants' responses to abiotic stresses implementing the stress models in forest species reported below. As there is no specific information on a metabolic pathway in forest species, we considered herbaceous species to acquire more specific information of the selected genes to check if they are involved in abiotic stresses Using our results, we have developed a specific model for each abiotic stress (heat, cold, drought, and salt). We have followed the general scheme suggested by Harfousche et al. [3] as this is one of the most comprehensive for illustrating the main generic response pathways to abiotic stresses.

Heat Stress Model
Sudden increases in temperature trigger the stress response in many different organisms, such as bacteria, fungi, and plants. In plants, it has been observed that initial exposure from medium to high heat stress provides resistance against a subsequent lethal dose. This phenomenon is called acquired tolerance [38,39]. The heat stress response is characterised by an elevated synthesis of a specific set of heat shock proteins (Hsps) associated with the development of thermotolerance [40][41][42]. Among the genes involved in the heat stress response, we observed the gene encoding for Hsp70 to be positively regulated by heat stress (Table 2). In eukaryotes, including plants, many different members of the Hsp70 family have been observed. From the information obtained from the systematic review, we observed that members of the Hsp70 family exhibit functions attributed to heat and cold stress and seed maturation and germination ( Figure S1, Table 2). Together, these data suggest that plant Hsp70 proteins interact with diverse substrates and take part in a plethora of cellular processes in Glycine max L. and Zea mays L [43][44][45]. Additionally, concerning the heat stress response, we observed the regulation of the gene encoding for cytosolic small heat shock protein (sHsps) (Figure 3). This gene appears to be over-expressed in response to this stress. Based on the information obtained from the systematic review, we observed that sHsps are a group of proteins whose production is ubiquitous in the cell in response to heat stress. In fact, these proteins are the most dominant proteins produced in response to such stress [38]. These proteins exhibit typical chaperone-like activity; in fact, they bind to unfolding intermediates to protect them from irreversible aggregation and maintain them in a competent refolding state. Expression patterns and chaperone-like function of sHsps suggest that their production correlates with the thermotolerance status acquired by cells following mild heat stress in Nicotiana tabacum L. [46]. One of the genes that we have observed to be regulated in response to heat stress is the gene encoding for ADP-glucose pyrophosphorylase large subunit (AGP-ase) ( Figure S1). In contrast to what has been reported for the previous two genes, the expression of AGP-ase seems to be negatively regulated by heat stress conditions. Indeed, from the information obtained as a result of the systematic review, we observed that AGP-ase activity is decreased, inhibited, or both in response to heat stress ( Figure S1, Table S4) in accordance with what was reported by Shayanfar et al. [46] in Triticum aestivum L. and Mangelsen [47] in Hordeum vulgare L. Kaur et al. [48] reported that AGP-ase activity is directly correlated with starch biosynthesis in Triticum aestivum L. The decrease in AGP-ase activity during heat stress appears to be directly related to the redirection of carbon flux away from starch biosynthesis pathways. Inhibition of this gene appears to be directly related to the decrease in 3-phosphoglycerate caused by increased respiration. Another gene that we observed to be down-regulated in response to this stress is the gene encoding for Light harvesting complex II protein (LHC-II). As can be seen from Table 2, the cumulative effect size for this gene is negative (−4.79). Based on the information obtained from the systematic review, we observed how heat stress conditions could directly affect the structure of this protein (Table S3). Accordingly, Song et al. [49] and Shakeel et al. [50] have also observed a down-regulation of LHC-II in response to elevated temperatures in Populus simonii Carrière and Agave americana L. An interesting finding reported by Vayghan et al. [51] concerns the phosphorylation status of the gene in response to elevated temperatures. This epigenetic pattern is usually linked with increased gene expression. Indeed, in the study by Vayghan et al. [51], the gene encoding for LHC-II was found to be up-regulated in Lepidium sativum L. Figure 7 shows the model of response to heat stress considering the results obtained from the meta-analysis of our study.  [52]; HPT/AHP2 Histidine-containing phosphotransfer ABF [52]; CNGC6, Cyclic nucleotide gated Ca 2+ channel 6 [53][54][55][56]; MID1, Ca 2+ -permeable mechanosensitive channels mid1-complementing activity 1 [57]; phyB, Photosensory receptor B [56,58]; TMS10, TMS10L, Thermo-sensitive genic male sterile 10 [5,56,59]; CaM3, Calmodulin 3 [56,60]; CBK3, CaM-binding protein kinase 3 [56,60]; CRLK1, Calcium/Calmodulin-regulated receptor like cytoplasmic kinases 1 [5,56,61]; CRLK2, Calcium/Calmodulin-regulated receptor like cytoplasmic kinases 2 [5,61]; DREB2A, DRE-binding protein 2A [62]; HSFs, Heat stress transcription factor [62]; WRKY39, WRKY transcription factor 39 [2,62]; bZIP17, Basic doamin/leucine zipper17 Transcription factor [2,62]; bZIP28, Basic doamin/leucine zip-per28 Transcription factor [2,62]; bZIP60, Basic doamin/leucine zipper60 Transcription factor [2,62].

Cold Stress Model
Low temperature is one of the major abiotic stresses limiting the growth and development of many plant species and affecting their geographic distribution. Arabidopsis thaliana L. can increase freezing tolerance in response to low temperatures through extensive gene expression reprogramming events that result in appropriate metabolic-structural alterations [63,64]. The result of these reprogramming events is the increased levels of hundreds of metabolites with protective and signalling effects [64]. Through the analysis of the results of the systematic review, we observed the behaviour of some key genes (Cytosolic class I small heat-shock protein, Light-harvesting complex II protein, S-adenosylmethionine synthase, S-adenosyl-l-homocysteine hydrolase) for the response to these environmental conditions (Table 2, Figure 8). These include the previously mentioned genes encoding for Hsps ( Figure 3) and sHsps ( Figure S1). Hsps are commonly associated with response to high temperatures, but evidence indicates that they can also respond to low temperatures in Arabidopsis thaliana [65,66]. Renaut et al. [67] observed that two Hsp70-like proteins were up-regulated in response to cold stress in Prunus persica L. As already observed for heat stress, the role of Hsps is to prevent the aggregation of denatured proteins, facilitating their refolding under cold stress conditions [63]. Among the genes found to be regulated in response to multiple stresses, we observed the gene encoding for Light-harvesting complex II protein (LHC-II). From the observation of the results obtained from the systematic review, in fact, the expression of the LHC-II gene appears to be induced by cold stress ( Table 2). In some studies on tomatoes and tobacco, the up-regulation of this gene in response to low-temperature stress has been observed [68][69][70]. Up-regulation of LHC-II allows a reduction in ROS content, acting as an antioxidant. This effect protects the biological membrane system and protects photosystem-II from photoinhibition. After the application of the systematic review and meta-analysis, we observed cold-stress-induced regulation of two interesting genes coding for S-adenosylmethionine synthase (SAMS) and S-adenosyl-l-homocysteine hydrolase (SAHH) ( Figure S1). S-adenosyl methionine is an important protein that serves as a universal donor of methyl groups. Being the keystone of one of the most important epigenetic patterns, this protein is fundamental during the growth and development phases of plants. SAM has, in fact, a fundamental role in the methylation of DNA, RNA, and polyamine and biotin biosynthesis [71,72]. Scientific evidence has shown that SAM can interact in the regulatory pathway of environmental stimuli response through epigenetic modifications and hormonal control in plant cells, as demonstrated in Oryza sativa L. and Arabidopsis thaliana L. [73]. SAM is synthesised through the action of the enzyme S-adenosylmethionine synthetase (SAMS) that catalyses this synthesis from ATP and l-methionine. From the observation of our results, SAM seems to be up-regulated during cold stress ( Table 2, Table S4). Studies of Guo et al. [74], Mahatma et al. [75], and Heidari et al. [76] demonstrated that in plants such as Medicago sativa L. [74], Cajanus cajan L. [75], and Solanum lycopersicum L. [76], a gene encoding for SAMS was up-regulated in response to cold stress. Induction of SAMS during stress increases SAMS levels as positive effects of ethylene and polyamine biosynthesis [71]. Because it appears to correlate with hormone synthesis, induction of SAMS activity may elevate ethylene levels. Ethylene might be involved in cell wall thickening [73]. According to Poulton et al. [77], the availability of the SAM enzyme is a key prerequisite for methylation. However, these methylation reactions have by-products that can inhibit methyltransferase activity [78,79]. One of these by-products is S-adenosyl-l-homocysteine (SAH). This protein can compete with SAM for the same binding site. Therefore, during developmental stages and during the stress response, this protein must be efficiently removed. To ensure the proper methyl-transferase activity of SAM, SAH is rapidly hydrolysed by another key enzyme during the stress response, namely S-adenosyl-l-homocysteine hydrolase (SAHH), which catalyses the hydrolysis of SAH into l-homocysteine and adenosine (Ado) [79]. In fact, the results of the systematic review and meta-analysis showed that in Arabidopsis thaliana plants subjected to cold stress, there is over-regulation of SAHH ( Figure S1), in accordance with Puyaubert et al. [80]. The results obtained in our study permitted us to depict the model of response to heat stress reported in Figure 8.  [52]; HPT/AHP2 Histi-dine-containing phosphotransfer ABF [52]; CNGC6, Cyclic nucleotide gated Ca 2+ channel 6 [53][54][55][56]; MID1, Ca 2+ -permeable mechanosensitive channels mid1-complementing activity 1 [57]; COLD1, Chilling tolerance divergenece 1 [56,81]; RGA1, G-protein alpha subunit 1 [56,81]; CIPK7, CBL-interacting protein kinase 7 [56][57][58][59][60]; CRLK1, Calcium/Calmodulin-regulated receptor like cytoplasmic kinases 1 [5,56,61]; CRLK2, Calcium/Calmodulin-regulated receptor like cytoplasmic kinases 2 [5,56,61]; ICE1, Induced CBP expression 1 [82]; CBF3, C-repeat-binding factor 3 [82];DREB1A, DRE-binding protein 1A [62]; ERF, Ethylene-responsive transcription factor [82]; AP2, APETALA2/ethylene response factor [82]; RAP2.1, CBF-regulation transcription factor 2.1 [82]; RAP2.6, CBF-regulation transcription factor 2.6 [82]; STZ/ZAT10, C2H2-type zinc finger transcription factor [82]; CCA1, MYB-type transcription factor [82]; LHY, MYB-type transcription factor [82].

Salt Stress Model
The salinity of the soil, often caused by NaCl, has a negative impact on the growth and development of plants [3]. This stress is represented by ion imbalance and hyperosmotic stress, leading to disorganisation of the cell membrane, ion toxicity, and oxidative damage [3]. In addition, salt stress causes oxidative stress due to the generation of reactive oxygen species (ROS) [94]. Regarding this kind of stress, through systematic reviews, we found interesting aspects concerning the behaviour of glyceraldheyde 3-phosphate dehydrogenase (GAPDH), S-adenosylmethionine synthase (SAMS), and S-adenosyl-l-homocisteine hydrolase (SAHH) genes (Table 2, Figure 10). GAPDH is involved in glycolysis, the photosynthetic reductive pentose phosphate pathway, and in signal-transduction processes related to the perception/signalling of abiotic stresses. Therefore, to date, this enzyme is defined as a moonlighting protein [95]. It is possible that the enzyme contributes to increasing salt tolerance in plants maintaining higher recycling rates of ADP and NADP + to decrease ROS production, helping to maintain photosynthetic efficiency and plant development [96]. We observed from our results ( Table S4) that GADPH expression was positively influenced by salt stress in Thellungiella halophila, Medicago sativa, Arabidopsis thaliana, and rice [97][98][99]. Furthermore, glyceraldheyde 3-phosphate dehydrogenase participated in the polyamines (PAs)mediated salt stress responses of Arabidopsis thaliana roots [98]. Another interesting finding that we observed concerns the activity of two key genes, S-adenosylmethionine synthase (SAM) and S-adenosyl-l-homocysteine (SAHH). SAM, from our results (Table S4), appears to be up-regulated under salt stress in Lycopersicon esculentum. From the information obtained through systematic reviews (Table S3), it emerged that there is an association between lignin deposition in vascular plant tissues and up-regulation of the gene encoding for SAMS during salt stress [100]. In addition, in response to salinity stress, vascular plants activate mechanisms to increase root cell wall synthesis and modification [101] (Table S3). These mechanisms involve a strong methyltransferase activity that requires elevated SAM synthesis. As observed for cold stress, methyltransferase activity appears to be crucial to the stress response. The dynamic process involves three foundational enzymes: SAM, SAMS, and SAHH. As already observed for cold stress [80], upregulation of the gene encoding for SAHH was also observed for salt stress in Spinacia oleracea L. [102]. Figure 10 shows the model of response to salt stress implemented with the results of our study.

Conclusions
In this study, we were able to deepen the knowledge reported in some theoretical models for the response to abiotic stresses showing the usefulness of the meta-analysis approach. In fact, large-scale theoretical models, based on mathematical methods, are very useful to understand the physical and biophysical laws at the basis of stress response metabolic pathways [15,25], but cannot clearly represent the dynamic mechanisms that regulate gene expression in response to one or more stress conditions. These theoretical models present some knowledge gaps due to the lack of representation of transcription factors or, indeed, candidate genes.
The approach we used for the identification and representation of B-ketoacyl-coA synthase, Caffeic acid o-methyltransferase, Cytosolic class I small heat shock protein, Formate dehydrogenase, Glyceraldehyde-3-phosphate dehydrogenase, Heat shock protein 70, S-adenosylmethionine decarboxylase, Light harvesting complex II protein, S-adenosylmethionine synthase, Xyloglucan endotransglusylase hydrolase protein 23, Adenylate kinase, and ADP-glucose pyrophosphorylase large subunit genes (Table 1, Figures 7-10) for stress-responsive metabolic pathways, proved to be very useful and effective in collecting and analysing information to better understand their role in the response to abiotic stresses. The future direction of this type of analysis should certainly be to analyse entire databases of genes involved in abiotic stress responses in other trees to deepen the knowledge of molecular and physiological responses in relation to different environmental stimuli. Future research should lead to the development of theoretical-statistical models, where information is obtained with mathematical, physical, statistical, and biophysical models.

Data Availability Statement:
The data presented in this study are openly available in the supplementary MS Access database included with this submission: Table S3.