Transcriptomic responses to drought stress in the Patagonian southern beech Nothofagus alpina

Background Deciphering the genetic architecture of drought tolerance could allow the candidate genes identification responding to water stress. In the Andean Patagonian forest, the genus Nothofagus represents an ecologically relevant species to be included in different genomic studies. These studies are scarce in South American ecosystems however represent an important source of genomic data in order to interpret future climate-change environment scenarios of these emblematic forests. Here, we achieved the assemblage of the transcriptome of N. alpina while searching for key genes of activated or suppressed metabolic pathways in response to drought stress. Results De novo transcriptome assembly resulted in 104,030 transcripts. Following confirmation of drought conditions, based on reduction of leaf water potential and stomatal conductance, a differential gene expression analysis resulted in 2720 significantly expressed genes (1601 up-regulated and 1119 down-regulated). Enrichment analysis (over-representation analysis and gene set enrichment analysis) resulted in more than one hundred stress-responsive term ontologies (i.e. biological processes) and pathways. Terms such as response to abscisic acid and pathways such as plant hormone signal transduction or starch and sucrose metabolism were over-represented. Protein–protein interaction assessment resulted in networks with significantly expressed top common hub gene clusters (e.g. plant-type cell wall biogenesis among down-regulated or ABA-signalling among up-regulated). These networks evidenced important regulators at gene expression such as transcriptional factors. Conclusions Responses of N. alpina seedlings to drought stress were evidenced by the activation of several genes linked to GO biological processes and KEGG pathways, which were mainly based on over-expression of specific protein


Introduction
Climate change promotes increases in frequency and severity of extreme drought events in a large proportion of terrestrial ecosystems (Allen et al. 2010), affecting forest quality, structure, diversity and health (Senf et al. 2020;Vicente-Serrano et al. 2020;Hartmann et al. 2022).Based on large-scale transcriptomic responses, researchers have aimed to elucidate gene expression profiles in different drought stress conditions (Du et al. 2018;Chen et al. 2019;Pervaiz et al. 2021;Lin et al. 2022).Plant stress responses involve complex crosstalk between different regulatory levels (Krasensky and Jonak 2012;Woldesemayat and Ntwasa 2018), including metabolism adjustment and gene expression for physiological and morphological acclimation (Lin et al. 2022).Inducible genes protecting against environmental stress are involved in stomatal closure, osmolyte metabolism, antioxidant activity and phytohormone signalling (Yamaguchi- Shinozaki and Shinozaki 2005;Shinozaki and Yamaguchi-Shinozaki 2007;Long et al. 2019).
Drought, which impairs photosynthesis, could directly trigger mortality through carbon starvation (McDowell et al. 2008) or hydraulic failure (Blackman et al. 2019).Stomatal closure and photosynthesis reduction constitute responses to diminish water loss (Shinozaki and Yamaguchi-Shinozaki 2007) mediated mainly by a plant hormone (ABA, abscisic acid) as part of a cross-talk metabolic pathways, signalling and regulation of gene expression.This hormone plays a pivotal role in several physiological processes (Long et al. 2019).Response to stress also includes the regulation of gene expression codifying for transcription factors (TFs), protein kinases, and phosphatases (Hasegawa et al. 2000;Seki et al. 2002;Yamaguchi-Shinozaki and Shinozaki 2005;Shinozaki and Yamaguchi-Shinozaki 2007).
South American temperate forests range across a decreasing precipitation gradient and the Sub-Antarctic temperate forest occupies a wide extension from 34° 50′ S to the south of the Tierra del Fuego archipelago at 55° S (Roig and Villalba 2008).In these forests, the genus Nothofagus includes ten species dispersed along a patched distribution in mono-specific or mixed forests (Suarez 2010).Among them, N. alpina (Poepp. et Endl.)Oerst., is most extensively distributed west of the Andes, in Chile, with some eastward extension into Argentina.Because of its historical commercial importance, its populations have been persistently degraded and deforested (Donoso and Lara 1996).Additionally, future distribution models predict they will be further threatened by the consequences of climate change particularly at the eastern distribution range (Marchelli et al. 2017).Due to the high value of its wood, the species was included in domestication programs to both conserve its current natural forests and promote its cultivation.To date, genomic information on the south American Nothofagus genus is insufficient (Torales et al. 2012;El Mujtar et al. 2014;Hasbún et al. 2016;Estravis-Barcala et al. 2021), and additional knowledge will be necessary for developing management and conservation policies under climate change scenarios.Previous studies on N. pumilio (Poepp.et Endl.)Krasser from areas with contrasting precipitation regimes noticed differential adaptive capacities to drought that could be revealing local adaptation to their home environment (Soliani et al. 2021).In the case of N. alpina, a partial transcriptome previously obtained with no induction to stress (Torales et al. 2012) served as a basis for our work.Thus, the main objective of this work consisted in the first identification of drought stressrelated transcriptional profile of N. alpina to unravel potential candidate key genes and pathways associated with climate change stress-related responses.We hypothesize that Nothofagus seedlings response to soil water involves over-expression of genes related to abscisic acid and osmotic adjustment pathways, as well as repression of genes related to growth.As such, this study expands the understanding of the molecular basis of drought stress mechanisms, which have become a common concern for an increasing number of plant species, including South American forest trees, for which the genomic basis of important traits related to drought tolerance remains undeciphered.

Experimental design and eco-physiological variables
A pure stand of N. alpina at Quila Quina (Tren Tren hill; 40.17°SL, 71.44° WL, 980 m a.s.l.) with an average precipitation of 1100 mm yr −1 , located in the Lácar Lake watershed (Parque Nacional Lanín, Neuquén, Argentina) at its eastern distribution, was selected for seed collection.At this location, N. alpina occupies a separate niche than its congener N. obliqua (Mirb.)Oerst.which grows at lower altitudes.Even though natural hybridization cannot be discarded, the rates of F1 hybrids at high altitudes is below 3% and introgressive backcrosses are around 8% (El Mujtar et al. 2017).
Seeds collected from a pool of 25 trees in 2009 were used for seedling germination, resulting in plants with a representative morphology of the species.Then a greenhouse pot drought experiment was performed with two levels of water irrigation: stress and control.Both treatments tried to simulate natural conditions of drought and abundant rainfall.This experiment was carried out during two austral growing seasons: 2014-2015 (5-year-old plants) and 2015-2016 (6-year-old Varela et al. (2010).In both seasons, the seedlings used in the control treatment were the same, while those used in the drought stress treatment were replaced in season 2015-2016, to avoid any cumulative effects from the first season.
The pot volume used in the experiment was previously determined by Varela et al. (2010) by measuring mean daily transpiration rates of seedlings of N. alpina.A 3-L (20 cm height) pot size allowed us to ensure water availability above the Permanent Wilting Point (PWP, below 10% vol/vol) for at least one week without watering.This enabled us both to maintain adequate soil moisture levels in control treatments (field capacity close to 30% vol/ vol) and to generate fairly gradual soil desiccation in stress treatments (Passioura 2006;Fernández 2010).Soil volumetric water content (VWC) of the pot, plant water status through pre-dawn leaf water potential (Ψ pd ), and stomatal conductance (g s ) as a physiological proxy of CO 2 fixation were assessed every time of registration on 5 and 7-8 plants per treatment for the first and second season, respectively.VWC (% vol/vol) was registered three times a week using HydroSens CS620 with a two-rod sensor of 120 mm (Campbell Scientific, Logan, USA).Pre-dawn leaf Ψ pd (MPa) was measured in one leaf per seedling per treatment at 0, 17 and 35 days of season 1, and at 0, 31 and 42 days of season 2, using a pressure chamber (PMS 1000, PMS Instruments, Corvallis, Oregon, USA).The last day of each experimental season was considered as the maximum stress condition (Ψ pd < -1,5 MPa).Finally, g s (mmol H 2 O m −2 s −1 ) was measured using a SC-1 leaf porometer (Decagon Devices Inc., Pullman WA, USA) three times a week.The proportion of harvested leaves was low in relation to total leaf area, ensuring other physiological variables measured were not affected.In all cases, fully expanded leaves were selected from the top third of the seedlings.
To evaluate the effect of drought over time for each variable measured (i.e.VWC, Ψ pd and g s ), a two-way repeated measures ANOVA was performed.Statistics was run using the tidyverse, ggpubr and rstatix packages in R Software.Model assumptions (non-significant outliers, normality and sphericity) were checked in all cases.In concordance with a statistically significant interaction between treatment and time on VWC, Ψ pd and g s in both experimental seasons, the effect of the treatment variable was analysed at each time point.Finally, p-values were adjusted using the Holm-Sidak multiple testing correction method.Ψ pd and g s records were used for comparison with previous drought physiological studies in N. alpina (Varela et al. 2010) in order to characterize the level of drought stress observed in this study.

RNA sequencing, de novo transcriptome assembly and functional annotation
At each experimental season, fully expanded leaves were trimmed (one leaf per seedling) and compiled into sample pools (10 seedlings per treatment).Each sample pool was manually grounded with mortar and pestle under liquid nitrogen and total RNA extracted according to Estravis-Barcala et al. (2021).The integrity of the RNA was assessed in a 0.8% agarose gel, and its quantity and quality checked with NanoDrop (ThermoFisher Scientific) and BioAnalyzer 2100 Plant RNA Pico chip (Agilent) before proceeding with polyA enrichment library preparation.A stranded library was constructed with a sample pool from a control treatment for a de novo assembly, and four non-stranded libraries were prepared with one sample pool per treatment at each experimental season for gene expression analysis.Sequencing was performed on an Illumina HiSeq 2500 system at the Novogene sequencing facility in University of California, Davis, CA, USA.

Differential gene expression and functional enrichment analyses
Pre-processed reads from non-stranded libraries were mapped against the de novo transcriptome assembly for transcript per million (TPM) quantification using Salmon v.0.8.1 (Patro et al. 2017).Counts were clustered into gene-level TPMs using the R package tximport v1.0.3 and normalized using the variance stabilizing transformation (VST) method in DESeq2 v1.28 (Love et al. 2014).Consistency between biological replicates of each treatment was checked via principal component analysis using the R package plotPCA.Hierarchical clustering selecting complete linkage and Euclidean distance of differentially expressed genes (DEGs) was visualized with the R package pheatmap.Differential gene expression analysis between control and stress treatments was performed with DEseq2 considering a threshold FDR < 0.05, and a volcano plot was constructed for visualization of the resulting differential expression.The null hypothesis that the logarithmic fold change (LFC) for each gene expression equals zero between treatments was tested (e.g. the gene is not affected by drought stress).The significance of differential expression results was compared in an MA plot (i.e.M for minus in the log scale; A for average in the log scale).
Functional enrichment analysis was performed with the R package clusterProfiler (Yu et al. 2012).Arabidopsis whole genome annotation was used as the background (UniProtKB, http:// www.unipr ot.org/ downl oads, accessed 10-02-2022).Annotated transcripts with |log2FC|> 1.5 were included in an over-representation analysis (ORA) using the functions enrichGO and enrichKEGG with default arguments, and in a gene set enrichment analysis (GSEA) using the gseGO and gseKEGG functions with defaults arguments.Interactive graphs of biological processes, molecular function, cellular component, and KEGG pathways of DEGs were visualized using the R package enrichplot.
Network analysis was used to further explore the correlation between genes to identify hub genes (e.g.pathways repetitive genes) among DEGs at the system level in both groups (e.g. up and down-regulated).The Cyto-Hubba plug-in of Cytoscape and the STRING database were implemented to detect validated and predicted protein-protein interactions (PPIs, Szklarczyk et al. 2019) while detecting the presence of hub genes in different metabolism pathways within the networks (Smoot et al. 2011).The CytoHubba plug-in applies specific methods such as betweenness, bottleneck, closeness, clustering coefficient, degree, dmnc, eccentricity, epc, mcc, mnc, radiality and stress to rank nodes in networks (Chin et al. 2014).Based on the ranking by 12 different methods, common hub genes were detected.Finally, to better illustrate the interaction between hub genes and other common stress-responsive genes, a sub-network was created by filtering the neighbour nodes that connected with hub genes.

De novo transcriptome assembly, differential gene expression and enrichment analyses
A total of 34,153,157 raw reads were sequenced from the stranded library, resulting in 32,730,857 reads after trimming and eliminating low-quality reads (Q < 20).The de novo assembly resulted in 179,620 transcripts, which after elimination of redundancy totalled 104,030 transcripts.Blast assessments resulted in 57,497 annotated transcripts (55.27%).The differential gene expression analysis resulted in 2720 DEGs between treatments, being 1601 up-regulated and 1119 down-regulated.Both the volcano plot (Additional file 2: Fig. S1a) and the MAplot (Additional file 2: Fig S1b ) showed significant over-expression of genes.Furthermore, the PCA assessment of the biological replicates for each treatment (Additional file 2: Fig. S1c) resulted in 94% of the total variance, with 80% and 14% explained by axes 1 and 2, respectively.Axis 1 clearly evidenced a drought effect in terms of gene expression profile.Similarly, a hierarchical clustering of DEGs showed an expression associated to drought response represented in the heatmap (Additional file 2: Fig. S1d) as two groups of genes in two distinct clusters (control and stress).The ORA analyses of these DEGs (Fig. 2) were visualized by the representation of the top 11 enriched GO terms (Fig. 2a) as well as the top 8 KEGG pathways (Fig. 2b).Among top counts, the terms response to abscisic acid, response to water deprivation, and cell wall organization or biogenesis were observed, while for the KEGG pathways were biosynthesis of secondary metabolites, plant hormone signal transduction and starch and sucrose metabolism.
On the other hand, the GSEA on GO terms (Fig. 3a) showed six out of the top 15 significantly biological processes' enriched terms were associated with activated catabolism, while seven were linked to suppressed biogenesis/development mechanisms.Relevant biological processes' suppressed terms were photosynthesis, light reaction, fatty acid derivative metabolic process, and cell wall organization or biogenesis while others over-expressed under drought were response to abscisic acid and response to hypoxia.From the KEGG GSEA (Fig. 3b) main suppressed pathways related to carbon assimilation and photorespiration included photosynthesis, and glyoxylate and dicarboxylate metabolism.Other suppressed pathways linked mainly to cell wall and meristem biogenesis as well as to auxin-activated signalling included fatty acid elongation and ABC transporters, respectively.Conversely, a main pathway such as the plant hormone signal transduction was activated and also relevant activated metabolism linked to the aspartatefamily pathway that leads to the synthesis of amino acids were identified as cysteine and methionine metabolism and lysine metabolism.
Based on ranking results with cytoHubba, the network analysis scored 25 top common hub genes among the up-regulated genes (Fig. 4a) and 15 hub genes among the down-regulated (Fig. 4b).Establishing groups by similarity, a cluster named as ABA-signalling activated genes (e.g.ABI1, HAB2) was defined by their link to this hormone pathway.This phytohormone, that activates various other stress-associated genes such as LTI65 (Fig. 4a), that encodes a protein induced in expression in response to water deprivation such as cold, high-salt, and desiccation.Other activated cluster, the galactose metabolism (Fig. 4a) included a relevant up-regulated hub gene encoding the sucrose non-fermenting 1-related protein kinase (SnRK1) was KING1.Another cluster, osmotic adjustment/carbohydrate metabolism hub genes are the GolS1 and GolS2 (Fig. 4a) coding for enzymes (e.g.galactinol synthases) associated with osmotic adjustment to avoid oxidative damage (Shen et al. 2020).Although not a hub gene, a linked to this cluster, the RFS5 gene (Fig. 4a) involved in the synthesis of the osmolyte raffinose is linked to the cluster.On the other hand, the group of suppressed hub genes were related to the clusters plant-type secondary cell wall biogenesis, where for instance the hub gene IRX15-L (Fig. 4b) encodes a protein involved in catalysing the addition of sugars to the growing polymer in the biogenesis of the cell wall biological process.A second down-regulated cluster identified as very long chain fatty acid and related to the fatty acid elongation KEGG pathway, scored-hub are the genes of the families KCS (1 and 6) and CER (1-3) (Fig. 4b).

Discussion
The transcriptome based on RNA-seq was de novo assembled for the southern beech N. alpina.Gene enrichment was assessed for the first time in this emblematic Patagonian forest tree to infer molecular responses and unravel key genes or pathways explaining changes to drought conditions.As expected, functional enrichment analyses in seedlings allowed us to identify specific drought-linked key genes identified in biological processes as well as in molecular pathways observed previously in other trees and shrubs in response to drought stress (Taïbi et al. 2017;Moran et al. 2017;Du et al. 2018;Chen et al. 2019;Feng et al. 2021;Pervaiz et al. 2021).
To disentangle the genetic bases of the traits related to drought susceptibility, we ensured the stress condition based in two physiological responses indicators (Fig. 1).On the one hand, values below Ψ 50 (Fig. 1b), previously reported for this species, were determined (point of 50% of hydraulic conductivity loss) (Varela, et al. 2010).On the other hand, a diminished stomatal conductance under stress was also observed (Fig. 1c).These physiological responses were translated into differential patterns in gene expression between stressed and control plants, whose molecular bases allowed us to elucidate that over-expression is concentrated in a particular set of genes of specific biological processes (e.g.GOterms, Fig. 2a).Up-regulated genes involved in oxidative stress (Fig. 3a) such as the activation of metabolism of secondary metabolites, phytohormone signal transduction (e.g.abscisic acid, auxin), and catabolic processes of organic acids were identified.Conversely, down-regulated genes involved in lipid-metabolism, carbon assimilation (e.g.(GO) terms.GO1 response to alcohol; GO2 response to abscisic acid; GO3 response to acid chemical; GO4 response to water; GO5 response to water deprivation; GO6 repsonse to temperature stimulus; GO7 cell wall organization or biogenesis; GO8 secondary metabolic process; GO9 terpenoid metabolic process; GO10 disaccharide biosynthetic process; GO11 pollen sperm differentiation.b Kyoto Encyclopedia of Genes and Genomes (KEGG) terms.KEGG1 biosynthesis of secondary metabolites; KEGG2 plant hormone signal transduction; KEGG3 starch and sucrose metabolism; KEGG4 MAPK signaling pathway; KEGG5 galactose metabolism; KEGG6 fatty acid degradation; KEGG7 tyrosine metabolsim; KEGG8 tropane, piperidine and pyridine alkaloid biosythensis Fig. 3 Dot plots showing enrichment of the gene set enrichment analysis.a Gene ontology (GO) terms.GO1 isoprenoid catabolic process; GO2 terpenoid catabolic process; GO3 photosynthetic electron transport chain; GO4 photosynthesis, light reaction; GO5 photosynthesis; GO6 second-messenger-mediated signaling; GO7 fatty acid derivative metabolic process; GO8 wax biosynthesis process; GO9 plant-type secondary cell wall biogenesis; GO10 organic acid catabolic process; GO11 carboxylic acid catabolic process; GO12 monocarboxylic acid catabolic process; GO13 plant-type cell wall biogenesis; GO14 lipid biosynthetic process; GO15 cell wall biogenesis; GO16 plant-type cell wall organization or biogenesis; GO17 cell wall organization or biogenesis; GO18 small molecule catabolic process; GO19 lipid catabolic process; GO20 cell division; GO21 response to decreased oxygen levels; GO22 response to oxygen levels; GO23 response to hypoxia; GO24 cellular developmental process; GO25 cell differentiation; GO26 response to cold; GO27 response to abscicic acid; GO28 response to alcohol; GO29 response to temperature stimulus.b Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.KEGG1 tropane, piperidine and pyridine alkaloid biosynthesis; KEGG2 photosynthesis; KEGG3 phenylalanine metabolism; KEGG4 tyrosine metabolism; KEGG5 glucosinolate biosynthesis; KEGG6 fatty acid elongation; KEGG7 diterpenoid biosynthesis; KEGG8 glyoxylate and dicarboxylate metabolism; KEGG9 motor proteins; KEGG10 ABC transporters; KEGG11 galactose metabolism; KEGG12 cysteine and methionine metabolism; KEGG13 plant hormone signal transduction; KEGG14 riboflavin metabolism; KEGG15 pyrimidine metabolism; KEGG16 lysine biosynthesis; KEGG17 MAPK signaling pathway Fig. 4 Protein protein interaction networks: a Up-regulated genes.b Down-regulated genes.Genes coloured in a yellow-red rank constitute hub genes according to the gene set enrichment assessment (GSEA) analysis.Hub genes related to main metabolic enriched terms are marked with a solid line.Black stars indicate DNA transcription factors photosynthesis), energy production (e.g.glycolysis and tricarboxylic acid cycle), and cell wall development/biogenesis were observed.Unravelled key genes that were significantly enriched are reported for both, encoding regulatory proteins such as gene expression/signal transduction (Fig. 4a, b).Similarly, key genes were deciphered encoding functional proteins such as key enzymes for osmolytes (e.g.raffinose) or hormone biosynthesis that respond to and tolerate stress (Fig. 4a).Regulatory proteins were identified as several transcription factors and DNA-binding factors (Additional file 1: Table S1), e.g.ERB (ethylene responsive), bHlH (leucine zipper), HSF family (heat-shock), AFP2 (ABI5 binding protein), WRKY).ERB (ERFs family members), HSF (several heat stresses factors), and WRKY transcription factors have responded similarly in N. pumilio against heat stress (Estravis-Barcala et al. 2021).In response to drought, WRKY was already identified within a complex network composed of several TFs (Lakra et al. 2013).DNA-binding proteins were also scored in the network as hub genes (Fig. 4a) such as the homeobox-leucine zipper protein (ATHB-7), identified as a regulator in an ABA-dependent pathway, while the heat stress TF (HSF4) was linked to HS-proteins regulation (Henriksson et al. 2005;Valdés et al. 2012).On the other hand, functional proteins such as enzymes (kinases or phosphatases) or transporters were evidenced through key hub genes.For instance, the enriched gene KING1 of the osmotic adjustment/ carbohydrate metabolism (Fig. 4a) elucidates that phosphorylation of the sucrose non-fermenting 1-related protein kinase (SnRK) that functions as inactivator of the sucrose synthase during a sugar starvation event.Phosphorylation has been identified as a main process acting in response to drought as enhancers of stabilization and degradation of targeted proteins (Ajadi et al. 2020).Furthermore, many metabolism pathways such as regulation of fatty acid synthesis and assimilation of nitrogen, as well as regulation of stomatal aperture/closure, are linked to this process (Ding et al. 2015).Similarly, other kinases over-expressed in this assessment are the CIPKs (Fig. 4a), which support the ABA-signalling pathway and responds to drought (Chaves et al. 2003;Yu et al. 2016;Long et al. 2019) in concordance with abiotic stresses results (Su et al. 2020;Qiu et al. 2022).Among the enzyme phosphatase group, ABI1 (Fig. 4a) is a hub gene phosphatase protein, key in several signalling ABA-pathways such as stomatal closure, osmotic water permeability, droughtinduced resistance and rhizogenesis, response to glucose, high light stress, seed germination and inhibition of vegetative growth (Chaves et al. 2003).
Photosynthesis and cell growth are primary processes to be affected by drought (Chaves et al. 2009).One of the quickest responses is suppressing water loss through stomatal closure, evidenced by the hub gene RZPF34 (Fig. 4a) encoding E3 ubiquitin-ligase, which inhibits CO 2 assimilation and causes a cellular sugar shortage, triggering a chain response including damage proteins, lipids, and nucleic acids (Batista-Silva et al. 2019).Key functional proteins linked to photosynthesis depletion were probably evidenced with suppression of KEGG pathways like galactose metabolism, photosynthesis, pyruvate metabolism, and GO terms as photosynthesis or light reaction.For instance, pyruvate metabolism pathway included two down-regulated genes encoding cell growth proteins such as ABCB4 (encoding a ATP-binding transporter) and PID (encoding protein kinase), both part of the auxin-activated signalling pathway.Likewise, photosynthesis in optimal growing conditions, but under water stress, the alternative sources of energy are glycolysis, the tricarboxylic acid (TCA) cycle and amino acid catabolism (Avin-Wittenberg et al 2012).These metabolisms are probably a stress-response to a decreasing synthesis of proteins/regulators and the consequent decrease of the main source (sucrose) for energy fulfilment evidenced with up-regulated genes as DIN4 or F24J8.4 (Fig. 4, cluster of carbohydrates metabolism).
The metabolic pathway linked to the lipid and organic acid catabolic processes were over-expressed, while those linked to cell wall organization and biogenesis were suppressed (Fig. 3a).This suggests the repression of physiological processes such as growth.The GSEA evidenced the lipid metabolism regulation related to degradation of complex molecules (e.g.carboxylic acid catabolic processes and lipid catabolic processes).It also suggests the activation of complex molecules biosynthesis (e.g.planttype cell wall biogenesis).These responses could involve the enriched fatty acid degradation pathway in response to improve drought resistance, through the inhibition of lignin and long-chain fatty acids synthases, as observed by Gu et al. (2020).Finally, the suppressed KEGG pathway glucosinolate biosynthesis (Fig. 3a) observed was contrary to previous results that show an accumulation of these compounds in response to drought stress conditions (Arbona et al. 2013).However, this experiment provides evidence that supports the over-expression of processes related to secondary metabolism linked to drought defense responses, indicating that water deprivation may hypothetically affect plant biotic interactions.
Based upon the meta-analysis from the STRING database, the evidence depicted by PPIs showed a set of core drought-responsive genes group, encoding molecular chaperones, as heat shock (HS) proteins (e.g.J8 and HSP18.2,HSP20-like chaperones superfamily) (Fig. 4a).For instance, HSP70's protect cells against heat and oxidative stress (Usman et al. 2017).In the same network, close to these HS-protein cluster, several genes (e.g.GolS1, GolS2, SIP2, RFS5 and DIN10) (Fig. 4a) coding for enzymes (e.g.galactinol synthases) associated with osmotic adjustment to avoid oxidative damage caused by ROS via the synthesis of raffinose (Nishizawa et al. 2008;Li et al. 2020) were intimately linked.

Conclusions
In this study, new transcriptomic information is presented from a novel drought stress experiment in the seedlings of N. alpina, an ecologically important tree inhabiting Patagonian temperate forests of South America.Although Nothofagus species have evidenced gene selection signals linked to drought stress along precipitation gradients suggesting local adaptation (Soliani et al. 2021), ecological niche models revealed endangered adaptation against future climate changing scenarios (Marchelli et al. 2017).Here, candidate genes for drought stress response were identified, as well as their interactions, metabolic pathways and processes potentially involved.Main results denoted the over-expression (e.g.up-regulation) of protein kinases, phosphatases, synthases and transcription factors, suggesting signalling pathway activation through phytohormones to avoid desiccation and osmotic damage based on stomatal closure and subsequent activation of carbon metabolism.For instance, ABA-signalling pathways were activated through phosphatases, whereas photosynthesis mechanism, tricarboxylic acid cycle, and pyruvate metabolism were suppressed.This indicates a possible link between over-expressed ABA genes and the expression of metabolic-related genes such as protein phosphatase 2C 56 (ABI1).On the other hand, the suppression (i.e.downregulation) of ABC transporter B family member 4 (ABCB4), enriched as an auxin-signalling pathway and cell growth ontologies, suggests a link between hormone transport to roots and the expression of genes involved in cell growth metabolism (e.g.secondary wall biogenesis, cell wall dynamics and fatty acid biosynthesis).Simultaneously, activation of genes related to secondary metabolites such as specific amino acids could be a response to osmotic stress created by the obstruction of ion exchange through stomata.This information suggests a rapid response to minimize water loss while metabolism is activated against oxidative stress and energy depletion.Overall, our results provide new genome-wide data to understand how native non-model long-living trees of Patagonian forests would acclimate, and eventually, adapt to environmental changes.

Fig. 1
Fig. 1 Physiological variables: a Soil volumetric water content (VWC, % vol of H 2 O per vol of soil).b Pre-dawn plant water potential (Ψ pd ; MPa, solid lines) and theoretical water potential at 50% loss in hydraulic conductivity (Ψ 50 ; MPa, dashed line).c Stomatal conductance (g s , mmol H 2 O m −2 s −1 ).Green lines and dots indicate mean values (± s.d.) for plants under control treatment.Red lines and dots represent mean values (± s.d.) for plants under water stress treatment.Left side: 2014-2015 growing season.Right side: 2015-2016 growing season.Variables are plotted as a function of days since beginning of experiments.Asterisks indicate statistical differences.ns = no statistically significant differences

Fig. 2
Fig.2Dot plots showing enrichment of the over-representation analysis.X axis is the RichFactor of each term.Bottom right bar represents the enriched -log10 P value.Top right dots indicate the gene number and the colour of dots corresponds to the adjusted P-value.a Gene ontology (GO) terms.GO1 response to alcohol; GO2 response to abscisic acid; GO3 response to acid chemical; GO4 response to water; GO5 response to water deprivation; GO6 repsonse to temperature stimulus; GO7 cell wall organization or biogenesis; GO8 secondary metabolic process; GO9 terpenoid metabolic process; GO10 disaccharide biosynthetic process; GO11 pollen sperm differentiation.b Kyoto Encyclopedia of Genes and Genomes (KEGG) terms.KEGG1 biosynthesis of secondary metabolites; KEGG2 plant hormone signal transduction; KEGG3 starch and sucrose metabolism; KEGG4 MAPK signaling pathway; KEGG5 galactose metabolism; KEGG6 fatty acid degradation; KEGG7 tyrosine metabolsim; KEGG8 tropane, piperidine and pyridine alkaloid biosythensis