Salt Stress Tolerance in Casuarina glauca: Insights from the Branchlets Transcriptome

Climate change and the accelerated rate of population growth are imposing a progressive degradation of natural ecosystems worldwide. In this context, the use of pioneer trees represents a powerful approach to reverse the situation. Among others, N2-fixing actinorhizal trees constitute important elements of plant communities and have been successfully used in land reclamation at a global scale. In this study, we have analyzed the transcriptome of the photosynthetic organs of Casuarina glauca (branchlets) to unravel the molecular mechanisms underlying salt stress tolerance. For that, C. glauca plants supplied either with chemical nitrogen (KNO3+) or nodulated by Frankia (NOD+) were exposed to a gradient of salt concentrations (200, 400, and 600 mM NaCl) and RNA-Seq was performed. An average of ca. 25 million clean reads was obtained for each group of plants, corresponding to 86,202 unigenes. The patterns of differentially expressed genes (DEGs) clearly separate two groups: (i) control- and 200 mM NaCl-treated plants, and (ii) 400 and 600 mM NaCl-treated plants. Additionally, although the number of total transcripts was relatively high in both plant groups, the percentage of significant DEGs was very low, ranging from 6 (200 mM NaCl/NOD+) to 314 (600 mM NaCl/KNO3+), mostly involving down-regulation. The vast majority of up-regulated genes was related to regulatory processes, reinforcing the hypothesis that some ecotypes of C. glauca have a strong stress-responsive system with an extensive set of constitutive defense mechanisms, complemented by a tight mechanism of transcriptional and post-transcriptional regulation. The results suggest that the robustness of the stress response system in C. glauca is regulated by a limited number of genes that tightly regulate detoxification and protein/enzyme stability, highlighting the complexity of the molecular interactions leading to salinity tolerance in this species.


Introduction
Climate change is unequivocally driving major environmental struggles associated with greenhouse gas (GHG) emissions and their impact on global warming [1,2] and sealevel rise [3]. Together with the accelerated rate of population growth, these changes are imposing a substantial loss of biodiversity and arable land [4,5]. Salinization is one of the major global challenges for agricultural systems. Besides the rising seawater levels, prolonged drought events cause secondary soil salinization as a result of an imbalance between water input (irrigation or rainfall) and use (transpiration) [6]. In this context, the
NOD + plants expressed a higher number of genes compared to KNO 3 + plants, with a total of 20,278 genes at 0 mM NaCl and an average of 19,780 genes in samples exposed to salinity (Table 1). From those, an average of 16,023 genes was common to control and salinity-stressed plants. Similar to KNO 3 + plants, a decreasing number of common genes and an increasing number of DEGs were observed with increasing salinity. The percentage of significant DEGs at all salinity conditions was slightly lower than in KNO 3 + plants: 0.03% (6), 0.5% (104), and 1.3% (254), respectively in 200, 400, and 600 mMNOD + . At 200 mM NaCl, the number of up-and down-regulated DEGs was similar. At 400 and 600 mM NaCl, the percentage of downregulated vs. upregulated DEGs was 60% (62) vs. 40% (42), and 85% (217) vs. 15% (37), respectively.
Heatmaps were used to analyze the regulation pattern of the top 10 salt-responsive DEGs from 600 mM NaCl relative to the control. For each treatment, the normalized read counts of these genes were also plotted to allow the comparison of gene regulation over the salinity gradient ( Figures S4 and S5). Two main clusters were formed, one with samples from 0 mM and 200 mM NaCl, and another with samples from 400 mM NaCl and 600 mM NaCl. The top downregulated DEGs from KNO 3 + plants exposed to 600 mM NaCl included an association with response to hypoxia ( (Table 3; Figure S4). In NOD + plants, 43% (136 out of 313) DEGs were uniquely mapped to known proteins, four of which uncharacterized ( Figure S3B). Considering each salt condition, four (67%), 49 (47%) and 111 (44%) of mapped DEGs were identified in 200-, 400-, and 600-NOD + , respectively (Table 1). Among the annotated DEGs from 200-NOD + , the most upregulated gene was PEPR1 (log2FC: 3.59) which is involved in PAMP-triggered immunity (PTI) signaling, followed by GDH2 (log2FC: 3.32) that catalyzes the production of GABA, and CYP78A5 (log2FC: 2.89) required for the promotion of leaf and floral growth and the prolongation of the plastochron. The only downregulated gene was GLPR2.8 (log2FC: −3.25), involved in light signal transduction and calcium homeostasis via the regulation of calcium influx into cells. In 400 NOD + , the most responsive DEGs were GDH2 (log2FC: 5.09) encoding a glutamate dehydrogenase, and RL5 (log2FC: −7.79) related to DNA binding. In 600 NOD + , the most up-or downregulated DEGs were PEPR1 (log2FC: 4.59) encoding a Leucine-rich repeat receptor-like kinase involved in PTI signaling, and RL5 (log2FC: −7.79), respectively. The top 10 most responsive DEGs in the NOD + group are summarized in Tables 4 and 5 and the full list of DEGs can be found in Table S8. Regarding the top 10 up-or downregulated DEGs, two main clusters were also formed in NOD + , one including the 0 mM and the 200 mM NaCl groups, and the other including the 400 mM and 600 mM NaCl groups. However, unlike the situation in KNO 3 + plants, a clear separation of control and 200-NOD + was observed in the first cluster. Similar to KNO 3 + plants, the top downregulated DEGs in NOD + plants exposed to 600 mM NaCl were associated with response to hypoxia, DNA-and mannose binding, cell death and cell wall biogenesis and organization (Table 3 vs. Table 5; Figure S4 vs. Figure S5). Additionally, these DEGs were also related to defense responses (GO:0006952), leaf morphogenesis (GO:0009965), response to oxidative stress (GO:0006979), cellular response to amino acid stimulus (GO:0071230), phosphorylation (GO:0016310), defense responses to bacteria (GO:0042742), signal transduction (GO:0007165) and metabolic processes (GO:0008152) ( Table 5; Figure S5). Also comparable to KNO 3 + plants, in NOD + plants the top upregulated DEGs at 600 mM NaCl were related to abscisic acid-activated signaling (GO:0009738), alkaloid metabolism (GO:0009820) and regulation of stomatal opening (GO:1902456) ( Table 3 vs.  Table 5; Figure S4 vs. Figure S5). Additionally, they were also associated with response to sulphur starvation (GO:0010438), jasmonic acid (GO:0009753) and wounding (GO:0009611), protein transport (GO:0015031), chaperone cofactor-dependent protein refolding, transmembrane transporter activity (GO:0022857), carbohydrate metabolism (GO:0005975), cell growth (GO:0016049), and lumen of the endoplasmic reticulum (GO:0005788) ( Table 5; Figure S5). Table 4. Top 10 most down and upregulated differentially expressed genes (DEGs) in Casuarina glauca nodulated by nitrogen-fixing Frankia casuarinae Thr (NOD + ), grown at 400 mM NaCl relative to the control (0 mM NaCl), ordered by log2 fold-change (log2FC). NBCI gene symbol, protein name, log2FC and respective Gene Ontology (GO) molecular function(s) and biological process(es) retrieved from UniprotKB database. Protein names coloured in blue are also among the top 10 most down or upregulated DEGs, respectively, at 400 mM NaCl in non-nodulated plants supplied with KNO 3 (KNO 3 + ).  Only two out of the top 10 most downregulated DEGs in samples exposed to 400 mM NaCl were found in both plant groups, namely RL5 and FLA12, both related to planttype secondary cell wall biogenesis. Among the top upregulated DEGs, only the sulphur deficiency induced gene SDI1 was common to both groups. At 600 mM NaCl, half of the top downregulated DEGs were common to both KNO 3 + and NOD + , namely RL5, Lectin, LAC22, LECRK91, and LRK10L-1.4. These genes were respectively involved in the following biological processes: DNA binding, iron ion homeostasis and transport, lignin catabolic process, defense response to bacterium and oomycetes, positive regulation of cell death, positive regulation of hydrogen peroxide metabolic process, cellular response to hypoxia and phosphorylation. Regarding the top upregulated DEGs, only two were common to KNO 3 + and NOD + plants, namely DTX56, associated with the cellular response to carbon dioxide and regulation of stomatal opening, and HSP70, involved in cellular response to unfolded protein and chaperone cofactor-dependent protein refolding. Overall, differential expression values, estimated by log2FC, were always higher in downregulated than in upregulated DEGs in both plant series.

Differentially Top-Expressed Genes in KNO 3 + and in NOD + Plants
In KNO 3 + plants, Clust tool showed 10 different clusters ( Figure 4, Table S9). Clusters C0 to C3, C8, and C9 consistently presented a clear downregulation of DEGs at 600 mM NaCl relative to the control, with changes in the regulation pattern at 200-and 400 mM NaCl. Cluster C0 showed a progressive downregulation with increasing salinity from 200 to 600 mM NaCl, while C1 displayed an accentuated downregulation only after 400 mM NaCl. Although clusters C2 and C3 presented a sharp downregulation between 400 and 600 mM, they also showed an increasing upregulation of DEGs from control to 400 mM NaCl, which is much more evident in C3. The downregulation pattern in clusters C8 and C9 was sharp until 400 mM NaCl, nearly maintaining the expression levels at the highest salinity condition. In contrast, clusters C4 and C5 consistently presented a strong upregulation of DEGs from 0 to 600 mM NaCl. Overall, in clusters C6 and C7, a downregulation was observed at 400 mM NaCl, more accentuated in C7, and reverted at 600 mM NaCl in both cases.  DEGs in almost all clusters were associated with DNA-, metal-, or sugar-binding, catalytic activity, transport, signaling, defense against biotic and abiotic stress and response to stimuli. In cluster C0, they were also related to chloroplast components, cell death, aging, oxidation-reduction process, and response to osmotic stress and water deprivation. Cluster C1 was more associated with membrane and mitochondrial components, and response to cold and salt stress. Although presenting similar expression patterns, while C2 was mainly associated with membrane components and defense responses, C3 was related to stomatal movement and signaling pathways involving brassinosteroids, ethylene and abscisic acid. Cluster C4 was associated with membrane and chloroplastic components, stomatal opening, glucose import, response to carbon dioxide, jasmonic acid and auxin, heat acclimation, regulation of growth. However, cluster C5, which presented a similar expression pattern, was mainly related to the response to metabolic processes and oxidative stress. Clusters C6 and C7 were both related to membrane, Golgi apparatus, mitochondria, and chloroplast components. However, while C7 was associated with the DEGs in almost all clusters were associated with DNA-, metal-, or sugar-binding, catalytic activity, transport, signaling, defense against biotic and abiotic stress and response to stimuli. In cluster C0, they were also related to chloroplast components, cell death, aging, oxidation-reduction process, and response to osmotic stress and water deprivation. Cluster C1 was more associated with membrane and mitochondrial components, and response to cold and salt stress. Although presenting similar expression patterns, while C2 was mainly associated with membrane components and defense responses, C3 was related to stomatal movement and signaling pathways involving brassinosteroids, ethylene and abscisic acid. Cluster C4 was associated with membrane and chloroplastic components, stomatal opening, glucose import, response to carbon dioxide, jasmonic acid and auxin, heat acclimation, regulation of growth. However, cluster C5, which presented a similar expression pattern, was mainly related to the response to metabolic processes and oxidative stress. Clusters C6 and C7 were both related to membrane, Golgi apparatus, mitochondria, and chloroplast components. However, while C7 was associated with the cell cycle, response to cold, root hair elongation and cell wall organization, C6 was related to cellular respiration, defense responses and the response to stress, including cold and salt stress. Clusters C8 and C9 were associated with defense responses and membrane, but while C8 was deeply related to the microtubule, auxin homeostasis, and signaling, C9 was more related to cell wall, leaf abscission, and cell death. The full list of GO terms of the three main categories associated with each cluster of DEGs from KNO 3 + plants can be seen in Table S10. In NOD + plants, the Clust tool showed 11 different clusters ( Figure 5, Table S11). Clusters C0 to C5 and C10 presented a noticeable downregulation of DEGs at 600 mM NaCl relative to the control, with different patterns of regulation at 200-and 400 mM NaCl. Clusters C0 and C1 showed a behavior similar to that of the corresponding clusters in KNO 3 + plants, where the first had a progressive downregulation with increasing salinity and the second only showed an accentuated downregulation at 600 mM NaCl. Clusters C2 and C3 showed upregulation at 200 mM NaCl, which overall was reverted at higher salinity levels. In clusters, C4 and C5 the major downregulation happened rapidly between 200 and 400 mM NaCl, slightly higher in C4 and moderately reverted in C5 at 600 mM NaCl. In cluster C10, the downregulation occurred in two marked steps, between the control and 200 mM NaCl and between 400 and 600 mM NaCl. Cluster C6 presented a strong downregulation at 400 mM NaCl, which was inverted at the highest salinity concentration to levels close to the control. Cluster C7 showed a somewhat similar behavior. However, the upregulation at 600 mM was much more evident and the transcript levels were higher than those in the control. control and 200 mM NaCl and between 400 and 600 mM NaCl. Cluster C6 presented a strong downregulation at 400 mM NaCl, which was inverted at the highest salinity concentration to levels close to the control. Cluster C7 showed a somewhat similar behavior. However, the upregulation at 600 mM was much more evident and the transcript levels were higher than those in the control. Overall, like in the case of KNO3 + plants, the majority of clusters were associated with binding, catalytic activity, transport, signaling, the defense response and responses to stimuli and stress. Cluster C0 was related to membranes and mitochondria, and highly associated with responses to stimuli, namely salicylic, jasmonic, and abscisic acids, ethylene and chitin, and stress responses, specifically to wounding, hypoxia, water depriva- Overall, like in the case of KNO 3 + plants, the majority of clusters were associated with binding, catalytic activity, transport, signaling, the defense response and responses to stimuli and stress. Cluster C0 was related to membranes and mitochondria, and highly associated with responses to stimuli, namely salicylic, jasmonic, and abscisic acids, ethylene and chitin, and stress responses, specifically to wounding, hypoxia, water deprivation and salt stress. Cluster C1 included a set of genes associated with chloroplasts, defense response, oxidation-reduction process and cell death. Cluster C2 was mainly associated with the defense response and responses to stimuli and stress, namely cold, hypoxia, wounding, and auxin, also including some genes related to cell differentiation and membrane components. Clusters C3 and C4 were enriched in genes involved in cell death and redox homeostasis. Additionally, C4 was also related to chloroplast, aging, cell wall biosynthesis and organization, response to oxidative stress, wounding, and biotic stimuli. Cluster C5, which was associated with chloroplasts and membrane components, included also genes involved in response to auxin and biotic stimuli, leaf morphogenesis and metabolic processes. Cluster C6 represented DEGs associated with regulation of growth, cell differentiation, and responses to oxidative and salt stress. Cluster C7 was mainly related to cell cycle, protein phosphorylation and the Golgi apparatus. Cluster C8 included genes involved in oxidation-reduction processes, chloroplasts, stomatal opening, and cellular responses to carbon dioxide. Cluster C9, which was related to membrane components, contained also genes related to the response to jasmonic acid and biotic stimuli, DNA protection, and defense against biotic stress. Finally, Cluster C10 included genes associated with membrane components, response to ozone, phloem development, response to hypoxia and abscisic acid. The full list of GO terms of the three main categories associated with each cluster of DEGs from NOD + plants can be seen in Table S12.
Likewise, PPI networks of DEGs were found in the STRING database only at 600 mM NaCl in both plant groups. In 600-KNO 3 + plants, three DEGs associated with enriched biological processes were mapped in the PPI network, namely SAG101, CRT3 and TAO1 ( Figure S6). The first two, which were related to the regulation of response to stress and positive regulation of response to stimuli, were directly linked in the network. Moreover, in NOD + plants, these two genes were indirectly linked to RBOHD, which was associated with cell death. Another four DEGs in NOD + plants were mapped in the PPI network, namely SAG101, CRT3, TAO1 and PUB17 ( Figure S7). Again, the first two were directly linked in the network. However, in this case, while CRT3 is associated with cell death, SAG101 is related to defense signaling, response to bacteria, regulation of response to stress and positive regulation of response to stimuli. Moreover, PUB17, which was associated with the defense response, was indirectly linked to both of these genes.

Genetic Stability and Volatility among DEGs
CodonW analysis was performed to assess the genetic stability or volatility among DEGs. It was observed that genes upregulated in NOD + plants had a GC content of 69.35% with 58.79% GC content at synonymous sites (GC3). A significant positive correlation between % GC and the frequency of optimal codons (Fop) was found (r = 0.57), whereas a significant negative correlation was obtained between GC3 and the effective number of codons (Nc) (r = −0.62). On the contrary, the upregulated genes in KNO 3 + plants showed different features. In this group, GC was 50.02% with 51.03% GC3. Moreover, statistical analyses between GC and Fop along with GC3 and Nc did note demonstrate statistical significance. This indicated distinct genetic volatility among the upregulated genes present in KNO 3 + plants. Moreover, COG analysis revealed "cellular processing and signaling" as a major functional category followed by "information storage and processing" and "metabolism" in NOD + plants ( Figure 6A), whereas "information storage and processing" was found to be the major COG family in non-nodulated KNO 3 + plants ( Figure 6B). The upregulated proteins in NOD + plants were more energy-consuming than those in nonnodulated plants.

Discussion
The impact of salt stress in Casuarina glauca has been extensively studied in depth at the morpho-physiological and biochemical levels, as well as at the proteome and metabolome scale [7]. Altogether, these studies revealed the robustness of the stress-responsive mechanisms in this species, particularly regarding protective traits at the photosynthetic and membrane level, whose activation was triggered at early stages of stress. In order to conclude this broad analysis, here, we report on the effects of increasing salt concentrations (200, 400 and 600 mM NaCl) on the transcriptome of branchlets from non-nodulated (KNO 3 + ) and nodulated (NOD + ) C. glauca plants. RNA-Seq analysis yielded an average of ca. 25 and ca. 26 million clean reads for KNO 3 + and NOD + plants, respectively, corresponding to 86,202 unigenes with a N50 size of 2792 bp and 41% GC content. These values are within the range of those reported in similar studies in Casuarina, which vary according to the species, environmental condition, and plant organ. For example, in Casuarina equisetifolia seedlings exposed to cold, ca. 21 million reads corresponding to ca. 118,000 unigenes and a N50 value of 2827 bp have been reported [38]. Developing secondary wood tissue of the same species presented ca. 43 million clean reads, comprising ca. 27,000 unigenes with a N50 contig size of 750 bp [39]. On the other hand, in C. equisetifolia roots from plants exposed to salt, these values were ca. 75-95 million clean reads comprising ca. 53,000-72,000 unigenes [37]. Finally, Wang et al. [40] reported ca. 41 million clean reads comprising ca. 60,000 unigenes with a N50 value of 2832 bp in arbuscular mycorrhiza roots of C. glauca exposed to salt. These results, together with the high mapping coverage of our data (almost 95%) indicates that a high-quality transcriptome assembly has been generated for downstream analyses.
In both plant groups, the patterns of differentially expressed genes (DEGs) clearly separate two groups, (i) control and 200 mM NaCl-treated plants, and (ii) 400 and 600 mM NaCl-treated plants (Figure 4, Figure 5, Figures S1, S4 and S5). With a few exceptions, this separation from low-to-moderate impacts (0, 200 mM) and heavy impacts (400, 600 mM) was also reported for several morpho-ecophysiological parameters, like plant growth and water use efficiency, mineral contents, leaf gas exchanges, chlorophyll a fluorescence, thylakoid electron transport rates, photosynthetic enzymes, and structural carbohydrates [22], as well for differentially expressed proteins (DEPs) of the same experimental set [30]. Notably, although the number of total transcripts was relatively high in both plant groups (9930 in the control and an average of 19,645 genes in NaCl-treated KNO 3 + plants, and 20,278 in the control and an average of 19,780 genes in NaCl-treated NOD + plants), the percentage of significant DEGs was remarkably low, ranging from 6 (200 NOD + ) to 314 (600 KNO 3 + ), mostly downregulated ( Figure 1; Table 1). The results are within the range of those reported for the invasive halophyte, Phragmites karka, exposed to 150 mM NaCl: 305 and 289 DEGs in leaves and roots, respectively [41]. Nevertheless, in this case, the fraction of positively regulated genes was higher than that of the negatively regulated ones. On the other hand, the number of DEGs reported in C. equisetifolia seedlings was considerably higher: >6000 DEGs in roots from plants exposed to 200 mM NaCl for seven days [37] and >4000 DEGs in leaves of plants exposed to cold [38]. However, in both cases the fold-change, rather than the statistical significance, has been privileged. Such observations highlight the complexity of the molecular interactions leading to salinity tolerance in C. glauca which might be related to a rapid and robust stress-responsive system [7], probably linked to long-term ecological adaptation [38].
A minor proportion of DEGs was common to both plant groups, in response to 400 (10%) and 600 mM NaCl (28%) (Figure 3), while no common DEGs were identified at the lowest salt concentration. The fact that the nodulation process triggers a set of defenseresponses [42,43] and that at 200 mM NaCl the N2 fixation by Frankia is residual [23] may explain this difference. Nevertheless, in both plant groups, the increase of the number of DEGs was gradual along the salt gradient, likely reflecting changes towards acclimation to high salt levels [30].
Considering the two key stress concentrations, i.e., 400 and 600 mM NaCl, the stressresponsive genes were in general related to regulatory processes in both plant groups. The results are rather different from those reported for the C. glauca branchlets' proteome [30], and metabolome [25][26][27], where changes were associated to major physiological changes related to photosynthesis, membrane stability and osmoprotection mechanisms [22,[24][25][26][27]. This reinforces the hypothesis that C. glauca has a strong stress-responsive system with an extensive set of constitutive defense mechanisms [7,43] complemented by rapid induction of transcriptional changes at the initial stages of salt exposure [37,40]. Despite the differences associated with the molecular changes imposed during the nodulation process (highlighted above), in both plant groups the top-ten up-regulated DEGs were related to signaling, transport, refolding, and stomatal control (Tables 1-5; Figures S4 and S5). These included auxin-(400 KNO 3 + ) and abscisic acid (ABA)-activated signaling (600 KNO 3 + ; 600 NOD + ); Myb transcription factors and jasmonic acid signaling (400 NOD + ), regulation of sulphur utilization, transporter activity, and protein refolding (400-and 600-KNO 3 + ; 400and 600 NOD + ), and response to cellular CO2/regulation of stomatal opening (600 KNO 3 + ; 600 NOD + ). It is widely known that the mechanisms used by plants to cope with salt stress are highly dependent on a coordinated set of events that regulates growth, metabolism, and homeostasis [37,44,45]. In halophytic plants, gene regulatory networks are an essential part of salt tolerance. For example, in C. equisetifolia [37], Spartina alterniflora [46], and Puccinellia nuttalliana [47], transcription factors (TFs), and ion transporters were identified as key drivers of salt stress tolerance, playing an essential role in ion homeostasis, regulation of ROS scavenging and detoxification, water balance, and water transport.
The top up-regulated genes common to both plant groups were the sulphur deficiency induced gene, SDI1, the gene encoding protein detoxification 56, DTX56, and the gene encoding the 70 KDa heat shock protein, HSP70 (Tables 2-5). The accumulation of SDI transcripts is related to sulphur deficiency, utilization of stored sulphate pools and sulphur partioning [48][49][50]. It has already been shown that the expression of the sulfate uptake transporter SULTR3 is induced in Arabidopsis and Medicago in response to various abiotic stresses [51] and that this induction depends on ABA [52]. Given that levels of cysteine and glutathione tend to increase significantly under abiotic stress conditions, which was also confirmed for glutathione in C. glauca [24,30], an upregulation of sulphate use is a plausible response to long-term salinity stress. On the other hand, DTX belongs to the group of transporters involved in detoxification in plants [53][54][55], as well as in turgor-regulation and ABA efflux in drought tolerant Arabidopsis [56]. Thus, it seems reasonable to hypothesize that DTX56 is a key multifunctional gene involved in the modulation of C. glauca responses to high salinity through the regulation of signal transduction pathways and associated mechanisms leading to stress tolerance [57]. Finally, the induction of the expression of the chaperone HSP70 is likely related to the maintenance of protein folding/unfolding as well as degradation of misfolded and denatured proteins, a common mechanism observed in plants under stress, e.g., osmotic [58], drought [59], and heat [60].
To be highlighted also is the fact that in NOD + plants, glutamate dehydrogenase (GDH) 2 is upregulated under all three salinity levels tested (Tables 5 and 6, Supplementary Table  S8). This enzyme is generally involved in stress adaptation [61] and plays an important role in replenishment of the tricarboxylic acid cycle and generally at the interface between carbon and nitrogen metabolism [62]. In fact, these authors showed that during short term salt stress, post-transcriptional regulation of GDH led to increased enzyme activity, and overexpression of the corresponding gene, improving biomass accumulation. Interestingly GDH is upregulated in NOD + , but not in KNO 3 + plants. Abscisic acid 8'-hydroxylase initiates the degradation of ABA and is thus important for maintaining ABA homeostasis [63,64]. The upregulation of its gene in NOD + plants exposed to 400 mM NaCl is interesting as during salt stress, an upregulation of ABA synthesis, not degradation, would have been expected [65], supporting a gradual stomatal closure that reached negligible values at 600 mM NaCl [22]. Another result that might seem surprising is the upregulation of WAXY, the starch synthase responsible for the synthesis of amylose, in 400 KNO 3 + plants. It was previously established for rice seeds that salt stress leads to a decrease in amylose biosynthesis [66,67]. However, another study [68] found a relative increase of amylose for seeds of triticale wheat under salt stress. In any case, despite the differences between the two plant groups, overall, these results suggest that the robustness of the stress response system in C. glauca is regulated by a limited number of genes that tightly regulate detoxification and protein/enzyme stability.
One of the mechanisms used by C. glauca to cope with salinity is the reduction of growth ensuring relative high values of hydration (as evaluated by the branchlets relative water content, RWC), ion homeostasis, and photosynthetic functioning [22]. In this context, the down-regulation of genes involved in growth and developmental processes, e.g., histones (DNA), patellins (cell polarity and patterning) [69], or cysteine proteases such as XCP1 (differentiation of xylem vessels) [70], is not surprising. In addition, previous studies on the effects of salinity stress have often shown negative regulation of genes involved in cell wall biosynthesis, a part of the growth process, or up-regulation of genes involved in cell wall degradation [40,71]. Accordingly, in this study, we found down-regulation of a set of cell-wall related genes, e.g., pectin esterase/pectin esterase inhibitor (400 NOD + ), a polygalacturonase (400 KNO 3 + ), an arabinogalactan (400 NOD + and 400 KNO 3 + ), a cellulose synthase (600 KNO 3 + ), a laccase, (600 NOD + and 600 KNO 3 + ), and a regulator of cell expansion, cobra-like protein (600 NOD + ).
As mentioned above, the expression of several genes encoding TFs putatively involved in signal transduction was differentially regulated under salt stress. The striking effect here, however, was the consistent downregulation of RADIALIS-LIKE SANT/MYB TFs (RADIALIS-LIKE 3 and RADIALIS-LIKE5) in both NOD + and KNO 3 + plants (400 and 600 mM NaCl). These TFs form a small family in Arabidopsis [72], and one of them (RSM1), which is down-regulated under salinity, was recently shown to modulate seedling development in response to ABA and salinity [73]. The function of RSM1 in the response to long-term salinity stress was confirmed in Castor Bean (Ricinus communis) by Han et al. [35], who suggested that the downregulation of RSM1 by salt was regulated via epigenetic modifications of the promoter.
Hierarchical cluster analysis identified 10 (KNO 3 + ) and 11 (NOD + ) gene clusters (Figures 4 and 5), all including genes associated with binding, catalytic activity, transport, signaling, defense against biotic and abiotic stress, and response to stimuli. Among these, five clusters (C0-C2, C8-C9) were sharply downregulated in KNO 3 + plants exposed to 600 mM NaCl, while in NOD + plants the number of negatively regulated gene clusters was seven (C0-C5, and C10). Despite the differences between the two plant groups, likely related to the nitrogen source [23,24], the overall set of results shows that C. glauca is able to balance the transcriptional activity at 200 mM NaCl to levels similar to that of the control (Figures 4 and 5). At 400 mM NaCl the expression patterns tend to be more variable, while at 600 mM NaCl transcript levels are severely affected, i.e., ca. 60 % of the DEGs are downregulated in both KNO 3 + and NOD + plants. Indeed, PPI networks of DEGs were found only at 600 mM NaCl in both plant groups, mainly related to the response to stress and stimuli, as well as to cell death. In our view, despite the positive regulation of some stress-related genes, at the highest salt concentration the plant is unable to sustain the necessary transcriptional activity to cope with salinity. Comparing this trend with the overall impact of 600 mM NaCl on plant growth, eco-physiological performance, proteome, and metabolome [22,[24][25][26][27]30], it seems clear that under the experimental conditions used, this C. glauca ecotype tolerates at least 400 mM. Such tolerance is related to a constitutive defense system, an enhanced antioxidative and osmoprotectant status, as well as a tight mechanism of transcription regulation and post-transcriptional modification.
Finally, codon and amino acid usage of DEGs among KNO 3 + and NOD + plants revealed a distinct pattern. Compositional constrain and translational efficiency were found to be the major driving force in maintaining the codon usage indices among DEGs of NOD + plants. In this group, the presence of costly proteins with higher amounts of aromatic amino acids highlights the importance of energy rich amino acids. Moreover, COG analysis made it evident that in NOD + plants, the signaling proteins were mostly functional. This may be due to a strong influence exerted by the symbiosis with Frankia casuarinae Thr. While the bacteria inside nodules are mostly dead at 200 mM NaCl [23], the effect of nodulation on the plant should still be in place. In legumes, extensive research has revealed complex local and systemic control of root architecture during nodulation [74,75]. Consistent with these data, legume studies have shown that nodulation affects the response to salt stress with regard to pathways involved in photosynthesis, respiration, anion transport, and plant defense [40]. Similar processes have evolved in actinorhizal plants [76]. However, they were not yet analyzed in Casuarina spp. In short, it is not surprising that even if the microsymbionts in nodules do not survive salt stress, the symbiosis still affects the stress response.
In conclusion, the results highlight the complexity of the molecular interactions leading to salinity tolerance in C. glauca likely linked to long-term ecological adaptation and, in this case, independent of symbiotic Frankia.

Growing Conditions and Experimental Setup
Casuarina glauca clones were grown in Broughton and Dillworth's (BD) medium under the conditions previously described by Zhong et al. [77] and Tromas et al. [78]. Six-monthold plants were either nodulated by Frankia casuarinae Thr (NOD + ) or supplemented with mineral nitrogen (KNO 3 + ) and maintained in a walk-in growth chamber (10,000 EHHF, ARALAB, Portugal) under environmental controlled conditions of temperature (26/22 • C), photoperiod (12 h), relative humidity (70%), external CO 2 levels (380 µL L −1 ), and irradiance (ca. 500 µmol m −2 s −1 ). Salt stress was gradually imposed through the addition of 50 mM NaCl per week to the nutrient solution until concentrations of 200, 400 and 600 mM were obtained, respectively. Control plants were supplemented only with mineral nutrients without the addition of NaCl. In all cases, the nutrient solution was renewed twice per week. To ensure the same age of all plants at the time of sample collection, i.e., one week after the exposure to each NaCl concentration, stress implementation was imposed sequentially from the highest to the lowest salt condition, i.e., the 600 mM NaCl group was treated first, followed by the 400 mM group (when the 600 mM group surpassed 200 mM NaCl), and the 200 mM group (when the 400 mM group surpassed 200 mM NaCl). Branchlets from four plants per treatment were harvested, frozen immediately in liquid nitrogen, and stored at -80 ºC until RNA extraction.

Total RNA Extraction and Library Preparation
Total RNA from branchlets of NOD + and KNO 3 + plants exposed to the three stress levels, as well as from control plants was extracted using the GeneJet Plant RNA Purification Kit (Thermo Scientific, Waltham, MA, USA) and digested with dsDNAse (Thermo Scientific, USA) to remove genomic DNA, as per manufacturer's instructions. The quality of the RNA samples was verified by electrophoresis in 1% agarose-Tris acetate EDTA buffer containing GelRed Nucleic Acid Gel Stain (Biotium, Fremont, CA, USA), by evaluating the integrity of the 28S and 18S ribosomal RNA bands and absence of smears, as well as with a bioanalyzer (Agilent 2100, Agilent, Santa Clara, CA, USA), by determining the RNA integrity number (RIN > 8.5). Libraries were prepared with the TruSeq RNA Sample Prep Kit v2 (Illumina, San Diego, CA, USA) and sequenced on an Illumina HiSeq 2000 platform (2 × 125 bp pair-end reads; 30 million reads per sample) at Macrogen (Seul, South Korea). Raw data is available in NCBI Sequence Read Archive (SRA) BioProject SymbSaltStress under accession PRJNA706159 (https://www.ncbi.nlm.nih.gov/sra/PRJNA706159, accessed on 1 September 2022).

Processing and Mapping of Illumina Reads
The sequenced raw reads were assessed for integrity, quality and contamination through the application of the following tools: FastQC version 0.11.8 to analyze reads quality [79] and FastQ Screen version 0.13 to survey putative contaminants, run against the genome of 14 default pre-indexed species and adaptors [80]. Trimmomatic version 0.39.1 was then used to eliminate the remaining adaptors and low-quality or small reads [81], using ILLUMINACLIP with keepBothReads option, SLIDINGWINDOW:4:15, MAXINFO:36: 0.5 and MINLEN:36. After filtering and trimming, Trinity version 0.39.1 was used to perform de novo transcriptome assembly, combining all samples to generate one single assembly [82]. This software was developed specifically for short reads and is advantageous for non-model plant sequence assemblies [83]. The assembled transcriptome was assessed for completeness through the gVolante [84] online interface with BUSCO v2.0.1 option [85]. To align the reads against the assembled transcriptome, the sequences were processed with Trinity tool Bowtie2 version 2.3.5 [86,87] and the aligned reads for each condition were quantified at gene expression level with RSEM version 1.3.2. [88]. The normalized expression of all samples was estimated using the Trimmed Mean of M values (TMM). A principal component analysis (PCA) was performed to survey the relatedness of all samples using the function plotPCA in R studio version 3.6.0 [89].

Identification of Differentially Expressed Genes
To identify differentially expressed genes (DEGs), a set of R packages were applied, namely edgeR version 3.26.8 [90], DESeq version 1.36.0 [91] and NOISeq version 2.28.0 [92]. DESeq fits a generalized linear model to estimate variance-mean dependence in count data, testing for differential expression based on the negative binomial distribution. The dispersion was estimated through blind mode [91]. On the other hand, edgeR uses weighted likelihood methods to implement a flexible empirical Bayes approach to allow gene-specific variation estimates even with very few or no replicates [93]. The dispersion value, which must be fixed manually in the absence of replicates, was set to 0.1. By contrast, NOISeq is an exploratory analysis that tests for differential expression between two experimental conditions without parametric assumptions that can simulate technical replicates. This method relies on the premise that read counts follow a multinomial distribution, where probabilities for each feature are the probability of a read to map to it [92]. Differential expression was computed using a stringent threshold of q = 0.9 along with the following parameters: pnr = 0.2, nss = 5, v = 0.02.
To study the effect of salinity, these three tools were used to identify DEGs at each salinity level (200, 400, and 600 mM NaCl) compared to the control (0 mM NaCl), for KNO 3 + and NOD + plants, respectively. The results were adjusted with the Benjamini and Hochberg's approach for controlling the false discovery rate (FDR) [94]. A filter of FDR < 0.05 and normalized log2 fold change (log2FC) > |2| were set to define DEGs. DEGs commonly detected by edgeR, DESeq and NOISeq were combined to increase the accuracy of results, and only the genes detected as differentially expressed by all the three tools were used in downstream analyses. To visualize the resulting expression profiles, volcano plots, heatmaps, barplots, and Venn diagrams were prepared for each combination, using stats and graphics core R packages and Python's Matplotlib 3.2.1 library [95].
The Clust version 1.8.10 command-line tool was applied to visualize the expression patterns of the detected DEGs and to find co-expressed genes [96]. This is a fully automated method for the identification of gene clusters that are co-expressed in heterogeneous datasets that automatically normalizes the data and determines the best number of clusters, based on a selected tightness, which in this case was five (5). Then, the results were filtered to keep only the DEGs that were found uniquely in one of the clusters and to find terms related to salt stress response.

Functional Annotation and Enrichment Analysis
The Basic Local Alignment Search Tool (BLAST) version 2.9.0 command line application from the NCBI C++ Toolkit [97] was used for functional annotation of DEGs, using homology searches with the latest references, comprising only the expertly curated component of UniProtKB. Results were filtered by maximum E-Value of 1.0E-3 and minimum Identities of 40% [98]. The resultant annotated proteins were then characterized by Gene Ontology (GO) terms: Cellular Component (CC), Molecular Function (MF), and Biological Process (BP), using the Uniprot and QuickGO APIs [99] to retrieve direct terms and GO term ancestors.
An over-representation analysis (ORA) was implemented by g:GOSt functional profiling tool from gProfiler website, which was applied using g:SCS tailored algorithm that uses a minimum hypergeometric test (Fisher's exact test). ShinyGO version 0.61 webtool, which is also based on hypergeometric distribution followed by FDR correction was used to retrieve Protein-Protein Interactions (PPIs) from the STRING database [100]. In both tools, Arabidopsis thaliana was selected as the organism of interest and separate log2FC ranked lists of DEGs for each salinity condition were used as inputs.

Downstream Bioinformatics Analysis
The downstream analysis included the codon and amino acid usage analysis of DEGs. Compositional constrain (GC and GC3), effective number of codons (ENc), and frequency of optimal codons (Fop) were determined using CodonW software. Biosyn-thetic energy costs of DEGs were obtained through DAMBE software. The Cluster of Orthologue analysis (COG) was performed using an in-house python-based program where information regarding different COG families was taken from the COG database (http://www.ncbi.nlm.nih.gov/COG, accessed on 1 September 2022). A detailed study on biological networks maintained by DEGs was performed. and NCBI gene symbol and protein names were retrieved from NCBI blastx command line tool; Table S10: Full list of Gene Ontology (GO) biological processes, molecular functions and cellular components associated with each cluster of potentially co-expressed differentially expressed genes (DEGs) in branchlets of non-nodulated Casuarina glauca plants supplied with KNO 3 (KNO 3 + ), grown at 200, 400, and 600 mM NaCl. Clusters were formed by Clust software and GO names were retrieved from UniprotKB database; Table S11: Clusters of potentially co-expressed differentially expressed genes (DEGs) in branchlets of Casuarina glauca nodulated by nitrogen-fixing Frankia casuarinae Thr (NOD + ), grown at 200, 400, and 600 mM NaCl. Clusters were formed by Clust software and NCBI gene symbol and protein names were retrieved using NCBI blastx command line tool; Table S12: Full list of Gene Ontology (GO) biological processes, molecular functions and cellular components associated with each cluster of potentially co-expressed differentially expressed genes (DEGs) in branchlets of Casuarina glauca nodulated by nitrogen-fixing Frankia casuarinae Thr (NOD + ), grown at 200, 400, and 600 mM NaCl. Clusters were formed by Clust software and GO names were retrieved from UniprotKB database; Table S13: Enriched pathways among differentially expressed genes (DEGs) in branchlets of non-nodulated Casuarina glauca plants supplied with KNO 3 (KNO 3 + ), grown at 600 mM NaCl relative to control (0 mM NaCl). Pathways were retrieved from InterPro, Pfam and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases, through ShinyGO webtool, with a False Discovery Rate (FDR) threshold < 0.05; Figure S1: Principal component analysis of gene expression counts from Casuarina glauca nodulated by nitrogen-fixing Frankia casuarinae Thr (NOD + ) and non-nodulated plants supplied with KNO 3 (KNO 3 + ), grown at control (0 mM NaCl) and three salinity conditions (200, 400, and 600 mM NaCl). PCA was plotted using R function plotPCA using normalized expression matrices; Figure S2: Gene expression in branchlets of non-nodulated Casuarina glauca plants supplied with KNO 3 (KNO 3 + ) and nodulated with nitrogen-fixing Frankia casuarinae Thr (NOD + ), grown at control (0 mM NaCl) and three salinity conditions (200, 400, and 600 mM NaCl). Black dots represent genes without differential expression and red dots represent differentially expressed genes (DEGs) between control (0mM NaCl) and salinity-exposed plants: