International Journal of Molecular Sciences Article Transcriptomics Reveals Fast Changes in Salicylate and Jasmonate Signaling Pathways in Shoots of Carbonate-Tolerant Arabidopsis thaliana under Bicarbonate Exposure Laura Pérez-Martín , Silvia Busoms , Roser Tolrà and Charlotte Poschenrieder * Plant Physiology Laboratory, Bioscience Faculty, Universitat Autònoma de Barcelona, C/de la Vall Moronta s/n, E-08193 Bellaterra, Spain; laura.perez.martin@uab.cat (L.P.-M.); silvia.busoms@uab.cat (S.B.); roser.tolra@uab.cat (R.T.) * Correspondence: charlotte.poschenrieder@uab.es; Tel.: +34-93-581-2163 Abstract: High bicarbonate concentrations of calcareous soils with high pH can affect crop perfor- mance due to different constraints. Among these, Fe deficiency has mostly been studied. The ability to mobilize sparingly soluble Fe is a key factor for tolerance. Here, a comparative transcriptomic analysis was performed with two naturally selected Arabidopsis thaliana demes, the carbonate-tolerant A1(c+) and the sensitive T6(c−) . Analyses of plants exposed to either pH stress alone (pH 5.9 vs. pH 8.3) or to alkalinity caused by 10 mM NaHCO3 (pH 8.3) confirmed better growth and nutri- ent homeostasis of A1(c+) under alkaline conditions. RNA-sequencing (RNA-seq) revealed that bicarbonate quickly (3 h) induced Fe deficiency-related genes in T6(c−) leaves. Contrastingly, in A1(c+) , initial changes concerned receptor-like proteins (RLP), jasmonate (JA) and salicylate (SA)   pathways, methionine-derived glucosinolates (GS), sulfur starvation, starch degradation, and cell cycle. Our results suggest that leaves of carbonate-tolerant plants do not sense iron deficiency as fast Citation: Pérez-Martín, L.; Busoms, as sensitive ones. This is in line with a more efficient Fe translocation to aerial parts. In A1(c+) leaves, S.; Tolrà, R.; Poschenrieder, C. Transcriptomics Reveals Fast the activation of other genes related to stress perception, signal transduction, GS, sulfur acquisition, Changes in Salicylate and Jasmonate and cell cycle precedes the induction of iron homeostasis mechanisms yielding an efficient response Signaling Pathways in Shoots of to bicarbonate stress. Carbonate-Tolerant Arabidopsis thaliana under Bicarbonate Exposure. Keywords: transcriptomic; Arabidopsis thaliana; alkaline stress; glucosinolate; salicylate Int. J. Mol. Sci. 2021, 22, 1226. https://doi.org/10.3390/ijms22031226 Academic Editor: Giampiero Cai 1. Introduction Received: 11 January 2021 Several constraints can affect the growth, physiology, and metabolism of agricultural Accepted: 22 January 2021 crops cultivated on calcareous/alkaline soils of arid and semi-arid regions over the world. Published: 27 January 2021 The main anions present excessively in alkaline soils are HCO3 − and CO3 2− [1]. The pH values of most calcareous soils are within the range of 7.5 to 8.5 [2], and bicarbonate Publisher’s Note: MDPI stays neutral concentrations are found between 5–35 mmol L−1 . Soil carbonates act as pH buffers and with regard to jurisdictional claims in play an important role in rhizosphere processes such as nutrient availability for plants. The published maps and institutional affil- iations. high pH surrounding the plant roots can disrupt the membrane potential and inhibit the absorption of essential ions [3]. Moreover, high pH lowers the availability of nitrogen (N); phosphorous (P); and micronutrients such as iron (Fe), zinc (Zn), and manganese (Mn), producing nutrient deficiencies in sensitive plants [4]. Among these nutrient deficiencies induced under alkaline soil conditions, Fe has Copyright: © 2021 by the authors. received most attention. Although Fe is one of the most abundant elements in the earth’s Licensee MDPI, Basel, Switzerland. crust, under basic pH conditions, it is sparingly available to plants. Iron, among other This article is an open access article essential micronutrients, is required for multiple metabolic processes, and most crops distributed under the terms and conditions of the Creative Commons are highly sensitive to low iron availability [5]. This fact has prompted many efforts on Attribution (CC BY) license (https:// unrevealing the mechanisms of Fe uptake, transport, and regulation, and the tolerance to creativecommons.org/licenses/by/ Fe deficiency [6–8]. Iron scarcity is a relevant consequence of high pH soils, but not the 4.0/). Int. J. Mol. Sci. 2021, 22, 1226. https://doi.org/10.3390/ijms22031226 https://www.mdpi.com/journal/ijms Int. J. Mol. Sci. 2021, 22, 1226 2 of 24 only one. Deficiency of Zn, P, N, and/or Mg are common problems for crop plants on calcareous soils. Several studies have shown that plants perceive and respond in different ways to high pH, high HCO3 − , or iron deficiency [9], and the mechanisms of carbonate tolerance remain largely unknown. Alkaline pH per se can affect root growth by inhibiting cell division and elongation [10]; ethylene modulates this effect increasing the expression of AUX1 and of genes responsible for auxin synthesis [11]. Crosstalk between jasmonate and ethylene in response to alkalinity stress has been reported [12,13]. One way for plants to adapt to high pH is to accumulate small metabolites with buffering function, especially organic acids, for adjusting the internal pH value. However, the adjusting of external and internal pH is a process that consumes energy and reduces plant growth simultaneously [10]. Bicarbonate imposes further stress on plant performance. Common features are morpho-anatomical changes in roots; disturbance of water relations; reduction of total photosynthetic pigments and photosynthetic activity; accumulation of osmo-protectants, soluble sugars, and organic acids; and activation of the biosynthesis of antioxidant en- zymes [12,14]. At the molecular level, several anion transporters and gene families such as ALMT, NRT/POT, and SLAHs are related to bicarbonate response in Glycine max [15,16]. The overexpression of the anion transporters GsBOR1 (boron efflux transporter) and GsS- LAH3 (efflux of uncharacterized inorganic or organic anions) are also involved in HCO3 − tolerance. Transcriptomic experiments under alkaline stress in Tamarix hispida, Ziziphus acidojujuba, and Glycine max identified gene families and transcription factors (TFs) in- volved in NaHCO3 stress responses. These genes can be clustered in sensing and signal transduction pathways as calmodulin-like proteins, serine/threonine protein kinases, and cysteine-rich receptor-like protein kinases. P2/ERF, HD-ZIP, bHLH, MYB, WRKY, NAC, C2H2, HB, and TIFY are the main TFs [17–19]. Proteomic studies performed in Solanum lycopersicum comparing saline and alkaline stress further revealed downregulated proteins accounting for metabolism, energy conversion, and novel upregulated proteins involved in signaling or transport [20]. Nonetheless, the current understanding of plant response to bicarbonate stress is limited. Most studies have been performed on roots of crop species such as Glycine max or Oryza sativa [21,22]. Less information is available regarding the molecular-genetic mechanisms underlying naturally selected tolerance and fast bicarbonate responses at the gene expression level in shoots. The large genotypic and phenotypic diversity of natural populations of the genetic model plant Arabidopsis thaliana provides an ideal experimental scenario for exploring such mechanisms. In previous studies, we identified two contrasting demes of A. thaliana with differential tolerance to moderate soil alkalinity: the sensitive deme T6(c−) and the tolerant A1(c+) . The better performance of A1(c+) under iron deficiency conditions was mainly attributed to enhanced root release of coumarin-type phenolics [23]. However, there is increasing evidence for both quick shoot responses and shoot to root signaling under iron deficiency or bicarbonate stress [24,25]. The aim of this study was to characterize the differences in initial shoot responses of naturally selected A. thaliana demes differing in tolerance to carbonate soils. To clearly dis- criminate between the initial responses to high pH per se from those caused by bicarbonate, one of the most growth-limiting factors on alkaline soils, we submitted seedlings to control conditions (pH 5.9), high pH provided by 1,3-bis[tris(hydroxymethyl)methylamino]propane (BTP) buffer (pH 8.3), and alkaline conditions achieved by supplying bicarbonate (10 mM NaHCO3 , pH 8.3). For abbreviation, pH 5.9, pH 8.3, and bic are respectively used to describe these treatments. Primary gene responses were analyzed after 3 h of exposure, while a second analysis was performed after a medium-term treatment of 48 h. Differential expression of genes and their gene ontologies were compared between demes, treatments, and time points, giving a detailed picture of the short- and longer-term responses to pH 8.3 and bic stress of A. thaliana plants differing in adaptation to carbonated soil. Int. J. Mol. Sci. 2021, 22, 1226 3 of 24 2. Results and Discussion 2.1. Fitness Characterization of A. thaliana Plants under High pH and Bicarbonate Stress In previous studies, we identified two contrasting demes of A. thaliana with differential tolerance to moderate soil alkalinity: the sensitive population T6(c−) and the tolerant A1(c+) [23]. A1(c+) is native to areas near to calcareous soils, and soil CaCO3 content and pH in its habitat are higher than in native soils of T6(c−) , which are derived from silica parental Int. J. Mol. Sci. 2021, 22, x FOR PEERrock REVIEW(Figure 1A,B). 4 of 26 Figure 1. Location, soil parameters, physiological characterization, and plant nutrition. (A) Loca- Figure 1. Location, tion of Arabidopsis soilpopulations thaliana natural parameters, physiological A1(c+) and T6(c−) in a soil map characterization, of Catalonia. Orange and plant nutrition. (A) Location areas represent calcareous soils and white areas non-calcareous soils. (B) Soil carbonate and pH of Arabidopsis thaliana natural populations A1 (c+) and T6 levels (means ± Standard Error; n = 9) in the natural habitat of A1(c+) and T6(c−) populations. (C)(c − ) in a soil map of Catalonia. Orange areas Picture of 30-day-old represent calcareous , and T6 Col-0, A1(c+)soils (c−) plants under different alkaline treatments. Germina- and white areas non-calcareous soils. (B) Soil carbonate and pH levels tion rate (D), rosette diameter (E), and silique number (F) of A1(c+), T6(c−), and Col-0 plants submit- (means ± (pH ted to control Standard Error; 5.9, grey bars), high npH=(pH 9)8.3, in light thegreen natural habitat bars), or ofwith bic (pH 8.3 A110(c+)mMandof T6(c−) populations. (C) Picture of NaHCO3, dark green bars) in plates or hydroponic culture. Values are means ± SE; n = 45, n = 12, n 30-day-old = 9, respectively;Col-0, A1(c+) letters indicate , and differences significant T6(c−) plants(p < 0.05, under different Tukey’s honestly alkaline significant dif- treatments. Germination rate (D), ference HSD) per deme. (G,H,I) Standardized leaf mineral content of four plants per deme grown rosette diameter (E), and silique number (F) of A1 (c+) under pH 5.9 (grey lines), pH 8.3 (light green lines), or bic (dark green lines) for 2 weeks. , T6 (c − ) , and Col-0 plants submitted to control (pH 5.9, grey bars), high pH (pH 8.3, light green bars), or bic (pH 8.3 with 10 mM of NaHCO3 , dark 2.2. Nutrient Content Profile green bars) in plates or hydroponic culture. Values are means ± SE; n = 45, n = 12, n = 9, respectively; Differences in tolerance to pH 8.3 and bic treatments were also reflected in the plants’ letters indicate ionomic profiles. significant Nutrient concentrations differences in Col-0 were (p more < 0.05, Tukey’s affected honestly by bic than by pH significant difference HSD) per 8.3. In general, the tolerant deme A1(c+) maintained ion homeostasis better than the sensi- deme. (G–I) Standardized leaf mineral content of four plants per deme grown under pH 5.9 (grey tive plants (Figure 1G–I). Moreover, maintenance of ion homeostasis in A1(c+) was better lines), under the pH bic 8.3 (lightthan treatment greenunderlines), or bic pH 8.3. Only (dark P, Ca, and Mngreen lines) for concentrations were2 weeks. sig- nificantly reduced by bic in all demes. Statistical significance of treatment vs. genotype interactions (p-value < 0.05) are shown in Dataset S1. However, all tissue concentrations in A1(c+)In order exposed toremained to bic ascertain abovethe bic concentration the deficiency that efficiently thresholds for Brassicaceae species discriminates both demes under hydroponic conditions, we performed preliminary studies using a concentration [30]. Especially relevant is the ability of A1 (c+) to maintain leaf Fe concentrations under bic range from 0 to 20 mM NaHCO3 . Rosette diameter and root length growth of both demes was monitored for 5 weeks (Figure S1A,B, Dataset S1). Concentrations higher than 10 mM caused high mortality before ending the life cycle. This contrasts with more carbonate-tolerant crop plants such as Glycine max, where 50 mM NaHCO3 was established as a mild stress situation avoiding high mortalities [17]. According to our results, a concentration of 10 mM NaHCO3 was selected as an adequate discriminatory treatment in all subsequent experiments. The ability to germinate, i.e., root emergence, is a prerequisite for efficient plant estab- lishment on problem soils, and seed germination rate under stress conditions has frequently been used to assess plant tolerance to soil constraints [26]. Here, seed germination tests were performed, comparing the contrasting natural demes with the model accession Col-0. Int. J. Mol. Sci. 2021, 22, 1226 4 of 24 Seed germination tests revealed, for all accessions, a delay in germination in the high-pH treatments (pH 8.3 and bic) in comparison to pH 5.9. However, 10 days after sowing the germination rate was hardly affected in Col 0. A1(c+) seeds showed better germination rates under the high pH treatments than under pH 5.9. In contrast, both alkaline stresses reduced the germination rate of T6(c−) (Figure 1D). Accordingly, A1(c+) and Col-0 are more tolerant to alkaline conditions than T6(c−) . However, seed germination does not guarantee further growth, survival, and reproductive fitness under alkaline stress [27] Hydroponic studies using solutions with control pH (pH 5.9), high pH (pH 8.3), and bic (10 mM NaHCO3 , pH 8.3) confirmed the higher tolerance of A1(c+) to both pH 8.3 and bic. Col-0 exhibited an intermediate response, while T6(c−) was most sensitive. Both alkaline treatments, pH 8.3 and bic, induced severe chlorosis in leaves of T6(c−) (Figure 1C). Bic treatment also caused stronger chlorosis symptoms in Col-0 than in A1(c+) , while pH 8.3 apparently had little effect in both accessions. Rosette diameter data measured 15 days after the start of treatments and number of siliques further revealed alkaline sensitivity of both T6(c− ) and Col-0 and confirmed tolerance of A1(c+ ) (Figure 1E,F, Dataset S1). Interestingly, while in the sensitive accessions bic was more stressful than pH 8.3, the alkaline-tolerant A1(c+) performed better under bic than under pH 8.3 (Figure 1F). The observation that in sensitive plants bic is more stressful than pH 8.3 is in line with the inhibitory effects of bic on essential metabolic processes in plants [28]. Especially affected by bicarbonate are mechanisms activated under Fe deficiency conditions such as riboflavin biosynthesis genes and NRAMP1 in roots or bHLH38 and IRT1 in shoots, as shown in melon plants exposed to either Fe deficiency alone or in combination with bicarbonate [29]. 2.2. Nutrient Content Profile Differences in tolerance to pH 8.3 and bic treatments were also reflected in the plants’ ionomic profiles. Nutrient concentrations in Col-0 were more affected by bic than by pH 8.3. In general, the tolerant deme A1(c+) maintained ion homeostasis better than the sensitive plants (Figure 1G–I). Moreover, maintenance of ion homeostasis in A1(c+) was better under the bic treatment than under pH 8.3. Only P, Ca, and Mn concentrations were significantly reduced by bic in all demes. Statistical significance of treatment vs. genotype interactions (p-value < 0.05) are shown in Dataset S1. However, all tissue concentrations in A1(c+) exposed to bic remained above the deficiency thresholds for Brassicaceae species [30]. Especially relevant is the ability of A1(c+) to maintain leaf Fe concentrations under bic stress. According to previous investigations, this is due to the ability of this deme to mobilize Fe by the exudation of coumarins [23] acting as Fe reductants and chelators [31]. The analysis of fitness parameters and ionomic data revealed clear differences in the tolerance of A. thaliana plants to alkaline pH per se and to bic. In the sensitive accessions, Col-0 and T6(c−) , bic was more toxic than pH 8.3. T6(c−) evolved on siliceous soil with extremely low carbonate content. Contrastingly, A1(c+) , which has evolved on a moder- ately carbonated soil (12.13% carbonate content) with pH 7.36, was more tolerant to bic than to pH 8.3. Soil conditions can act as driving force for local adaptation [23] and the carbonate/bicarbonate concentration in the original soil of A1(c+) was strong enough to allow this accession to withstand the bicarbonate concentration provided in the nutrient solution (10 mM), yielding a pH of 8.3. However, this pH value achieved by organic buffer was harmful. It has previously been shown that buffers used to simulate alkaline stress may injure plants due to high salinity levels [32]. This was not the case here. The electric conductivity of both solutions was similar. We hypothesize that in the bicarbonate-tolerant deme, the buffer may not be as efficient as bic in triggering the expression of key genes required for nutrient acquisition and growth under bicarbonate stress. 2.3. Analysis of RNA-Seq Data RNA-seq analyses were performed to characterize, at the gene expression level, the mechanisms underlying the differential response to bicarbonate of the two naturally se- Int. J. Mol. Sci. 2021, 22, 1226 5 of 24 lected accessions. Early responses were assessed by using leaf samples after 3 and 48 h exposure when the plants still did not exhibit any foliar symptoms of the stress treatments. After Illumina sequencing, 72 library reads were obtained. Approximately 24,441 genes were detected to be expressed in each library (see Supplementary Datasets S2 and S3 for more details). To verify the RNA-seq transcriptome data, we selected eight genes with highly different expression levels in both demes after 3 h. The expression levels of the selected genes were quantified using q-PCR. Results obtained from both techniques were compared and we obtain high correlation between the expression level of the genes tested (Figure S2). Annotated genes were filter by log fold change (LFC > 1 and LFC < −1) and adjusted p-value < 0.05. To understand the differences between A1(c+) -tolerant and T6(c−) -sensitive lines in response to different treatments (pH 8.3 vs. pH 5.9 and bic vs. control pH 5.9) at two time points (3 h and 48 h), we performed pairwise comparisons. After filtering, a total of 6163 differentially expressed genes (DEGs) were identified, considering accession, time, and treatment (Figure 2). In general, the bic treatment caused a higher number of DEGs (5205 genes) than the high pH treatment (2420 genes), indicating that bicarbonate stress involves more complex processes than simply the specific responses to alkaline pH (Figure 2A). Exposure for 3h to bic induced differential expression of a higher number of genes in A1(c+) than in T6(c−) . Contrastingly, the highest number of DEGs was observed in the carbonate-sensitive accession T6(c−) when exposed for 48 h to bic (Figure 2A). Gene expression profiles can also be compared in the heatmap plots shown in Figure 2B. Many of2021, Int. J. Mol. Sci. these DEGs 22, x FOR that either were upregulated or downregulated by bic PEER REVIEW were not, or were 6 of 26 only slightly, affected by pH 8.3. Figure 2. Compared transcriptomic profiles under different alkalinity treatments at two time points. Differential Figure expressed 2. Compared transcriptomic genes profiles under(DEGs) bar plot different alkalinity (A), atheatmap treatments profile two time points. (B),expressed Differential and Venn diagram at 3 h (C) genes (DEGs) bar plot (A), heatmap profile (B), and Venn diagram at 3h (C) and 48 h (D) of the pairwise comparations pH and 48 h (D) of the pairwise comparations pH 8.3 vs. pH 5.9 and bic vs. pH 5.9 in A1(c+) and T6 8.3 vs. pH 5.9 and bic vs. pH 5.9 in A1(c+) and T6 (c−) demes. DEGs were filtered at log fold change (LFC > 1, LFC < −1), and (c−) demes. DEGs were filtered at log fold change (LFC > 1, LFC < −1), and adjusted p-value < 0.05. adjusted p-value < 0.05. Yellow areas indicate upregulated genes and blue areas indicate downregulated genes. Yellow areas indicate 2.4. Gene upregulated Ontology (GO) andgenes and blueof areas Kioto Encyclopedia Genes andindicate downregulated Genomes(KEGG) Pathways genes. under Bicarbonate Stress The general overview of differential gene expression comparing bic versus pH 8.3 reveals that a greater number of genes and pathways were specifically implicated in the bic responses. Moreover, A1(c+) and T6(c−) showed contrasted reactions toward bicarbonate stress. For this reason, we further focused our analysis on the DEGs and pathways found in the plants exposed to the bic treatment. Venn diagrams drawn for DEGs after 3 h (Figure 2C) and 48 h (Figure 2D) visualized DEGs that are shared between lines and time points under the different treatments. Within each deme, some of the DEGs were common to both bic and pH 8.3. This further supports the view that high pH alone cannot simulate bicarbonate stress. After 3 h exposure to bic, Int. J. Mol. Sci. 2021, 22, 1226 6 of 24 2.4. Gene Ontology (GO) and Kioto Encyclopedia of Genes and Genomes (KEGG) Pathways under Bicarbonate Stress The general overview of differential gene expression comparing bic versus pH 8.3 reveals that a greater number of genes and pathways were specifically implicated in the bic responses. Moreover, A1(c+) and T6(c−) showed contrasted reactions toward bicarbonate stress. For this reason, we further focused our analysis on the DEGs and pathways found in the plants exposed to the bic treatment. Venn diagrams drawn for DEGs after 3 h (Figure 2C) and 48 h (Figure 2D) visualized DEGs that are shared between lines and time points under the different treatments. Within each deme, some of the DEGs were common to both bic and pH 8.3. This further supports the view that high pH alone cannot simulate bicarbonate stress. After 3 h exposure to bic, more genes were differentially expressed in A1(c+) than in T6(c−) . Venn diagrams for DEGs comparing exposure times and demes for either pH 8.3 or bic are shown in Figure S3A,B, respectively. DEGs from each time point, treatment, and deme were mapped into the KEGG pathway and Gene Ontology (GO) terms to analyze their biological functions, molecular functions, and cell localization (Figure S3C,D, Datasets S4–S6). The highest log fold discovery rate (FDR) and number of DEGs corresponded to 48 h bic exposure of T6(c−) . Most affected genes were involved in biological regulation, cellular and metabolic processes, response to stimulus, cellular components, and catalytic activity. After 3 h exposure, genes related to fatty acid elongation and cutin, suberin, and wax biosynthesis were distinctive in the sensitive deme. Contrastingly, in the tolerant A1(c+) , genes related to 2-oxocarboxylic acid, alpha-linolenic acid, flavonoid biosynthesis, and glyoxylate and dicarboxylate metabolism were differentially regulated after 3 h exposure to bic. After 48 h, in this deme, genes related to carotenoid biosynthesis, photosynthesis, nitrogen metabolism, and carbohydrate metabolism were mostly affected. Venn diagrams comparing upregulated and downregulated DEGs between demes after 3 h bic exposure revealed that both demes only shared a small number of genes regulated in the same direction—41 up and 41 downregulated (Figure 3A). Differentially up- and downregulated genes in A1(c+) were mainly related to response to stimulus and catalytic activity (Figure 3C). Interestingly, in T6(c−) , most affected genes were related to biological processes (up and down), cellular components (down), and electron carrier activity (down). The 3 h exposure to bic caused upregulation of many metabolism-related genes in both demes (Figure 3C). Most notably, in the tolerant A1(c+) , a 30-fold enrichment in the upregulation of genes involved in glucosinolate biosynthesis was observed. Enrich- ments between 5- to 15-fold in upregulated genes related to 2-oxocarboxylic acids; alpha linolenic acids; fatty acid elongation; and valine, leucine and isoleucine were found. Some enrichment was further seen in upregulated genes of amino acid biosynthesis, biosynthe- sis of secondary metabolites cysteine and methionine metabolism, and other metabolic pathways. Downregulation in A1(c+) mostly affected genes related to base excision re- pair. In the sensitive T6(c−) , the most conspicuous enrichment concerned upregulation of genes related to fatty acid elongation, fatty acid degradation, and galactose metabolism (Figure 3C). Both upregulation and downregulation of genes related to secondary metabolism was observed. Downregulation also affected genes related to phenylpropanoid biosynthesis and plant hormone signal transduction (Datasets S7–S9). The large amount of DEGs in T6(c−) after 48 h exposure to bic (Figure 3B,D,F) indicated that this sensitive deme suffered from severe stress, leading to multiple secondary stress responses affecting almost all cellular and metabolic processes, even before visible foliar symptoms were detectable. Therefore, our further analytical efforts were focused on the differential responses between demes after 3 h exposure to bic (Datasets S7–S9). Int. J. Mol. Sci. 2021, 22, 1226 7 of 24 Int. J. Mol. Sci. 2021, 22, x FOR PEER REVIEW 8 of 26 Figure 3. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway of bic vs. pH 5.9 treat- Figure 3. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway of bic vs. pH 5.9 treatments ments DEGs at two time points. Venn diagram of up- and downregulated DEGs in A1(c+) and T6(c−) comparing bic vs. pH DEGs at two time points. Venn diagram of up- and downregulated DEGs in A1(c+) andofT6differentially 5.9 treatments after 3 h (A) and 48 h (B). Bubble plots indicating significant GO analysis (c−) comparing bic vs. pH 5.9 expressed genes treatments in bic vs. pH 5.9 comparison between A1(c+) and T6(c−) at 3 h (C) and 48 h (D). GO were filtered to adjusted p-value genes after 3 h (A) and 48 h (B). Bubble plots indicating significant GO analysis of differentially expressed < 0.05. in bic vs. ScalepHcolors 5.9 comparison between indicate number A1(c+) of DEGs andbubble while T6(c−) size at 3indicates h (C) and-log 48 of h (D). GO were adjusted filtered p-value. to adjusted GO terms p-value < were separated 0.05. into biological Scale function, colors indicate cellularofcomponent, number DEGs while and molecular bubble function. Arrows size indicates indicate up -log of adjusted or downregulated p-value. GO terms were genes. Heatmaps separated into of KEGG pathway analysis from DEGs in bic 8.3 vs. pH 5.9 comparison between A1(c+) and T6(c−) at 3 h (E) and 48 h (F). biological function, cellular component, and molecular function. Arrows indicate up or downregulated genes. Heatmaps of KEGG pathway analysis from DEGs in bic 8.3 vs. pH 5.9 comparison between A1(c+) and T6(c−) at 3 h (E) and 48 h (F). KEGG pathway terms were filtered by p-value < 0.05. Scale colors indicate pathway fold enrichment. BP, biological process; CC, cellular component; MF, molecular function. Int. J. Mol. Sci. 2021, 22, 1226 8 of 24 2.5. Signal Perception and Signal Transduction 2.5.1. Receptor-Like Kinases and Receptor-Like Proteins (RLK/RLP) The 3 h exposure to bic caused a considerable up- and downregulation of DEGs related to response to stimuli in both A1(c+) and T6(c−) . According to that, extreme DEGs among demes were selected (Dataset S10). Many of these genes are related to signal perception Int. J. Mol. Sci. 2021, 22, x FOR PEER REVIEW 10 of 26 and signal transduction. Therefore, heat maps for these genes were drawn (Figure 4A,B). Figure 4. Differential expression of signal perception, transduction, and transcription factor (TF) genes involved in bicar- Figure bonate 4. Differential short responses. Heatmaps expression of DEGs fromof the signal perception, receptor-like kinase transduction, andperoxidase, gene family (A), and transcription factor (TF) glutathione, genes involved and calmodulin in bicarbonate gene families short (B) in A1(c+) and T6(c+)responses. , comparing bicHeatmaps of DEGs vs. pH 5.9 exposure forfrom the receptor-like 3 h. Yellow color representskinase upregulated genes while blue color represents downregulated genes. (C) Bar plot indicating the total number of DEGs gene family (A), and peroxidase, glutathione, and calmodulin gene families (B) in A1(c+) and T6(c+) , upregulated (right) or downregulated (left) within each transcription factor family. Tolerant deme A1(c+) is marked in grey, whilecomparing sensitive deme bic vs.is pH T6(c−) 5.9 exposure visualized for 3 h. Yellow color represents upregulated genes while blue in black color. color represents downregulated genes. (C) Bar plot indicating the total number of DEGs upregulated (right) or downregulated (left) within each transcription factor family. Tolerant deme A1(c+) is marked in grey, while sensitive deme T6(c−) is visualized in black color. Int. J. Mol. Sci. 2021, 22, 1226 9 of 24 Specific analyses of genes mediating Fe deficiency-induced reversible protein phos- phorylation have revealed that in A. thaliana roots, the predominant gene groups code for proteins of the leucine-rich repeat kinase subfamily (RLK) [33]. Here, we found that 21 RLK/RLP genes were differentially expressed in shoots after 3 h exposure to bic (Figure 4A). In A1(c+) , bic distinctively upregulated genes coding for receptor-like proteins RLP 2, 32, 34, 37, and 39, and for receptor-like serine/threonine protein kinases SD1-6, also designed as ARK2, and SD18 (Figure 4A). Among these different RLPs, only RLP32 belongs to the validated defense RLPs [34]. This protein recognizes Inf-1, the bacterial translocation initiation factor of proteobacteria. Both RLP32 and RLP39 have been proposed to interact with SOBIR1; the interaction with SOBIR1 may be mediated via a stretch of negatively charged amino acids, such as aspartate and glutamate, in the apoplast close to the plasma membrane. Contrastingly, RLP2, also known as Clavata2-related gene, and the serine/threonine kinases SD1-6 or ARK2 are involved in developmental processes [35]. In both demes, bic upregulated RKL1, RLK902, and CRR9. RKL1 and RKL902 are two highly homologous receptor kinases apparently involved in mechanical and biotic stress signaling and salicylic acid-induced responses [36]. CRR9, chlororespiratory reduction 9, is a protein required for the formation of a NADH dehydrogenase-like complex [37]. The protein is involved in cyclic electron transport around photosystem I, a process favored under Fe deficiency [38]. In T6(c−) , bic specifically upregulated CRRSP7, a putative cysteine- rich repeat secreted protein. All other RLK/RLP genes were downregulated by bic in both demes. 2.5.2. Peroxidases Most of the 22 genes related to peroxidases, which were differentially expressed after 3 h bic treatment, were downregulated in both demes (Figure 4B). This is in line with reduced activity of cell wall peroxidase observed by others in Fe-deficient plants [39]. Downregulated only in T6(c−) were PER4, 5, 38, 54, and 58, coding for extracellular per- oxidases involved in oxidative stress; PER52 involved in lignin biosynthesis; and PER37 located in vascular bundles acting as a negative regulator of growth [40]. However, in the sensitive T6(c−) , PER39, 63, and 66 were specifically upregulated. The corresponding peroxidases are also located in the apoplast and are related to oxidative stress. PER66 has been found to be involved in xylem formation and lignification [41]. Specifically down- regulated in A1(c+) were PER11, 39, and 66, which are related to detoxification of reactive oxygen species (ROS), and PER17, involved in lignin biosynthesis. Interestingly, PER35 was specifically upregulated in A1(c+) . PER35 has been associated with transcription-associated mark trimethylation of H3 lysine 4 (H3K4me3) [42] and is essential for gene regulation during cell transition processes [43]. Taken together, these results reinforce the crucial role of differential expression of peroxidase genes in the early leaf responses to bicarbonate. Recent research by Santos et al. (2019) [44] hypothesized that Fe-efficient plants do not activate the antioxidant machinery in leaves because of better Fe translocation to the shoots, while inefficient plants suffering from oxidative stress induced by Fe deficiency show reduced ascorbate peroxidase (APX) activity. 2.5.3. Glutathione-Related In both A1(c+) and T6(c−) , exposure to bic upregulated GPX1 (Figure 4A). This chloro- plastic glutathione peroxidase plays an important role in degradation of H2 O2 under stress conditions [45]. Under bic exposure, main DEGs related to glutathione were glutathione-S- transferases (GSTs) of the tau type (GSTU). Most of these GSTU genes were downregulated in both demes (GSTU3, 4, 6, 10, 11, 12, and 15) (Figure 4A). Specifically upregulated in A1(c+) were GSTU17 and GSTU26, while GSTU25 and GSTU28 were specifically upregu- lated in T6(c−) . GSTs of the phi type, namely, GSTF6, 10 and 12, were downregulated in both demes. The functions of GSTUs are poorly explored. They are generally described as being associated with toxin catabolic processes. Different GSTUs (GSTU1, 3, 4, 8, 10, 11, 12, and 24) have been found strongly upregulated after fungal or bacterial attack [46]. Interest- Int. J. Mol. Sci. 2021, 22, 1226 10 of 24 ingly, certain GSTUs and GSTFs have been involved in the glucosinolate metabolism of Brassicaceae [47–49]. GSTF6 and GSTU4 seem to catalyze the GSH conjugation of indole-3- acetonitrile, a step in camalexin biosynthesis, while GSTF10 has been proposed to act in the biosynthesis of indolyl-glucosinolates. GSTU20 may act as a GSH donor in the biosynthesis of aliphatic glucosinolates [50]. Bic specifically upregulated GSTU25 and 28 in T6(c−) (Figure 4A). GSTU strongly binds to hydroxy C6 conjugates of fatty acids and it has been proposed that this enzyme may catalyze the formation of GSH conjugates of fatty acids and promote the transport to other compartments [51]. The differential expression of genes coding for glutathione- S-transferases under bic exposure suggests a general downregulation of glucosinolate biosynthesis in T6(c−) and a shift from indolyl-GS to aliphatic GS in the tolerant A1(c+) deme. 2.5.4. Calmodulin Bicarbonate exposure negatively regulated the IQ motif-containing protein IQM2 and the calmodulin-like protein coded by At5g10660 in both demes (Figure 4A). Specifically downregulated in T6(c−) were CML23, CML37, CML41, and SARD1 (SAR deficient 1). Instead, in A1(c+) , both SARD1 and CML41 (calmodulin-like 41) were specifically upregu- lated. The calcium-dependent protein kinase 2 (CPK2) was specifically upregulated only in T6(c−) and downregulated in A1(c+) . These results indicate clear differences in early Ca2+ - dependent signaling pathways between the bic-tolerant and the sensitive deme. Specific upregulation of CPK2 in T6(c−) suggests enhanced phosphorylation of nicotinamide ade- nine dinucleotide phosphate (NADPH) oxidases inducing ROS signaling and programmed cell death [52]. However, in this bic-sensitive deme, both jasmonate (JA) and salicylate (SA) signaling pathways seemed downregulated. This is supported by the observed downregu- lation of CML37, required for jasmonic acid (JA) production [53], and of SARD1, coding for a key transcription factor in salicylic acid (SA) biosynthesis [54] Contrastingly, in A1(c+) , bic-induced enhancement of both SARD1 and CML41 indicates activation of biotic stress signaling. CML41 can be induced by flagellin and promotes the deposition of callose in plasmodesmata, thus conferring resistance to Pseudomonas syringae [55]. Interestingly, studies with SA-deficient mutants have shown that accumulation of endogenous SA levels are required for activation of the Fe deficiency response. However, under Fe deficiency, pad4, a mutant with deficient SA biosynthesis, maintained higher shoot Fe levels and higher chlorophyll concentrations than Col-0 wildtype [48]. 2.5.5. Transcription Factors Exposure to bic had a strong influence on DEGs coding for transcription factors (TFs) (Figure 4C). Surprisingly, most were downregulated. In A1(c+) , 22 TFs were upregulated and 105 were downregulated; in T6(c−) , upregulation and downregulation were 37 and 70, respectively. Among the TFs exclusively upregulated in the tolerant A1(c+) , we found the Auxin Response Factor, coding for ARF20, which acts as repressor of gene transcription in response to auxin in the shoot apical meristem [56]. Upregulated only in A1(c+) were also two Heat Shock Factors—At4g19630 coding for an unknown protein, and At2g26150 coding for HSFA2, which responds not only to heat stress, but also to salinity and osmotic imbalance [57]. Among C3H type TFs, one coding for ZDW4 with unknown function was upregulated, while those coding for C3H14 and CDM1 were downregulated. Interest- ingly, C3H14 and CDM1 are regulators of the biosynthesis of hemicellulose and callose, respectively [58,59]. Specifically downregulated by bic in A1(c+) , but unaffected in T6(c−) , were the fol- lowing TFs: three Apetala2 (AP2s), coding for AIL7, ANT, and WRK14; two Cysteine-rich polycomb-like (CPPs), coding for TCX2 and TCX6; one E2F/DF, coding for DEL3; and one Teosinte branched1/Cincinnata/proliferating cell factor TCP, coding for TCP1 (Figure 4C). Except for WRK14, which is involved in the control of cuticular wax biosynthesis under stress, all others are related to cell proliferation and shoot growth regulation. AIL7 is required for high PIN1 (auxin efflux carrier component1) responsible for leaf shape, ANT controls Int. J. Mol. Sci. 2021, 22, 1226 11 of 24 cell proliferation, TCX2 has been involved in fate transition and cell division in stomatal development, and DEL3 plays a role in cell elongation and differentiation [60]. This spe- cific downregulation of TFs involved in plant growth promotion indicated that under bic-exposure, the shoot growth is quickly adjusted to the stress situation. However, this downregulation after 3 h did not result in reduced shoot growth after longer term exposure. Most probably, the downregulation is a sign of transitory changes leading to an adaptive growth response, allowing better performance under bic than in the sensitive T6(c−). In both demes, five TFs of the basic helix-loop-helix family (bHLH) were upregulated and five were downregulated. Several TFs of this family play a central role in Fe homeostasis under Fe deficiency conditions [61–63]. However, under bic exposure only PYE was upregulated in both demes. PYE is a key regulator of plant Fe deficiency responses. PYE is a strong repressor of FER coding for ferritin, which tightly binds Fe. FER repression thus reduces Fe sequestration and enhances the Fe fraction available for essential cell processes [62]. Contrastingly, only in T6(c−) were bHLHs of the Ib clade, namely, bHLH39, 38, 100, and 101, were upregulated by bic. These TFs activate FIT, the key regulator for Fe uptake and transport. Furthermore, Ib clade bHLHs promote leaf expansion in the shoot [63]. FIT repressor bHLH11 was downregulated in T6(c−) . In A1(c+) , upregulated bHLHs were bHLH23, 114, and 136 (PRE1), while bHLH10 and 45 (MUTE), SPCH, and HEC1 were downregulated. These TFs seemed to not be directly related to Fe deficiency response pathways but to growth and developmental processes. PRE1 mediates brassinoesteroid regulation of cell elongation [64]. HEC1 targets WUS and promotes cell proliferation in the shoot apical meristem [65]. SPCH and MUTE are key regulators for stomatal size and density in response to environmental changes [66]. A highly contrasting influence of bic on MYB TFs in both demes was observed. Only HOS10, a R2R3-type MYB TF, was upregulated in A1(c+) . Nine other MYBs were downregulated (MYB1, 10, 23, 46, 53, 86, 120, and 122, and MYB3-R4). MYB3-R4 is a transcriptional activator of cytokinesis [67]. MYB122 is required for indolyl-glucosinolates acting downstream of JA, SA, ethylene, and abscisic acid (ABA) [68]. These results further confirm that A1(c+) exposure to bic has a quick influence not only on shoot growth and developmental processes, but also on biotic stress signaling pathways. A1(c+) exposure to bic enhanced the transcription of WRKY70, 54, and 18, while WRK23, 44 (TTG2), and 45 were downregulated. Inversely, in T6(c−) , only downregulation of WRKYs was observed, namely, WRKY12, 14, 30, 50, 51, 62, and 71. WRKYs participate in regulatory networks of biotic and abiotic stress responses, either activating or repressing transcription [69]. In biotic stress signaling, SA-induced WRKYs can repress JA signaling, either directly or indirectly [70]. Overexpression of SA-induced WRKY70 represses the expression of PDF1.2 induced by MeJa [71,72]. WRKY23 regulates auxin transport by a local regulation of flavonol biosynthesis [73]. Among the Ethylene Responsive Factors (ERFs), 4 were upregulated and 11 were down- regulated in A1(c+) (Figure 4C). In T6(c−) , upregulation and downregulation were 4 and 7, respectively. Among the upregulated ERFs in A1(c+) , we found ERF13, 98, and 109. ERF13 responds to exogenous ABA and JA and increases expression of the defensin gene PDF1.2 [74]. ERF98 has been related to ascorbic acid biosynthesis [75], while ERF109 acts downstream of auxin, promoting cell growth under stress conditions [76]. Further influ- ence of bic on growth regulation in A1(c+) is indicated by the downregulation of Cytokinin Response Factors (CRF)1, 2, 5, and 6. Disruption of these TFs yields larger leaf rosettes and delay of senescence [77]. Interestingly, ERF72 was also downregulated in the tolerant deme. ERF72 has been described as quickly upregulated under Fe deficiency; the enhanced expression has been related to the regulation of chlorophyll degradation and Fe deficiency- induced chlorosis [78]. TFs related to drought resistance were also downregulated in A1(c+) , such as Translucent Green (TG), an activator of aquaporins [79]; DREB19 [80]; and SHINE2, involved in epicuticular wax production. Differently, in T6(c−) , bic induced upreg- ulation of stress-related ERFs—ERF54 and ESE2—and ERF39 related to primary cell wall Int. J. Mol. Sci. 2021, 22, 1226 12 of 24 biosynthesis [81]. Contrastingly, downregulation of ESR1 (Enhancer of Shoot Regeneration 1); stress-related ERF14; and ERF96, a positive regulator of ABA response [82], was observed. Bic treatment caused downregulation of eleven NAC TFs in A1(c+) ; only NAC44 and NAC95 were upregulated. In T6(c−) , only NAC71, 90, and 95 were downregulated. Most NACs are involved in plant developmental processes, phytohormone signaling, and biotic and abiotic stress responses [83]. NAC90 is one of the “NAC troika genes” that act as regulators of the shift to leaf senescence. NAC42 and 87 are targets of this NAC troika [84]. The considerable differences between A1(c+) and T6(c−) in the bic-induced changes in TF expression reveal quick regulation of growth and developmental processes in the tolerant deme in response to bic. These adaptive responses are mainly orchestrated by SA, JA, auxin, and ABA and induced by biotic and abiotic stress signaling pathways. Upregulation of PYE helping to mobilize Fe from ferritin was the only typical Fe deficiency response in A1(c+) . Contrastingly, in the sensitive T6(c−) , bic exposure caused fast FIT- regulated Fe deficiency responses in line with the impaired ability to translocate Fe to the developing shoots. 2.6. Protein–Protein Interaction Network Functional Enrichment Analysis with Specific Bic-Induced Genes The analysis of the expression of genes coding for stress signal perception/transduction and transcription factors revealed that shoots of both demes perceived and responded to carbonate stress in different ways. To further explore these differential mechanisms, we performed protein–protein interaction network functional enrichment analysis using STRING database (Figures 5 and 6). After 3 h exposure to bic, 273 genes were differentially expressed in A1(c+) and 123 genes in T6(c−) . In both demes, part of the genes was not linked to others due to lack of information on these protein interactions (Datasets S11 and S13). Nonetheless, we could observe clear differences between the tolerant A1(c+) (Figure 5) and the sensitive T6(c−) (Figure 6). In A1(c+) , many of the DEGs could be grouped into seven interconnected, functional categories related to SA, JA, GS, cell cycle, carbo- hydrate metabolism, sulfur deficit, antioxidants, and multidrug and toxic compound extrusion (MATE) efflux (Figure 5, Dataset S12). Inversely, in the sensitive T6(c−) , the main DEGs were grouped into five categories: Fe homeostasis (up); abiotic and biotic stress response; cuticular wax, nutrient, and transport activity; lipid transfer; and oxidative stress response. Only nutrient transport activity, lipid transfer, and oxidative stress response were connected by STRING (Figure 6, Dataset S14). Shen et al. (2016) [48] reported that Fe deficiency induces SA accumulation, which enhances auxin and ethylene signaling and activates Fe transport by bHLH38/39 regulation of downstream iron genes. In fact, in the sensitive deme T6(c−) , both the FIT-dependent [85] and the FIT-independent pathways [86] for iron homeostasis were activated by bic, as indicated by the upregulation of bHLH39 and bHLH100, respectively. Upregulation of NAS4 coding for nicotianamine synthase and NRAMP4 coding for Fe transporter exporting Fe from vacuoles further confirms the quick activation of Fe mobilizing mechanisms in the bic-sensitive deme. In the tolerant A1(c+) , bic enhanced SA, as shown by the upregulation of SARD1 [54] and several SA-induced genes such as WAK1 and At2g25510. SA, in turn, induced WRKY70, a key regulator of SA-induced genes and a repressor of JA-responsive genes [70]. In our tolerant plants, bic exposure caused a downregulation of ORA59, which is essential for the integration of JA/ethylene transduction pathways. ORA59 is required for the expression of defensin PDF 1.2, a typical JA signaling gene [87]. Contrastingly, VSP2, another JA-responsive marker gene, as well as AOS, AOC1, and OPR3 involved in JA biosynthesis, were upregulated (Figure 5). STRING connected AOS and OPR3 to GGP1, a γ-glutamyl peptidase that catalyzes the hydrolysis of the γ-glutamyl residue of GSH- conjugated glucosinolates. GSH is the sulfur donor for GS biosynthesis and GGP1 is a key enzyme in this process [88]. Upregulation of BCAT4, FMOGS-OX1, AKHSDH1, IMD1, and GSTU20 is in line with a bic-induced enhancement of aliphatic GS. There is a close connection between mineral nutrient supply and GS biosynthesis [89]. Sulfur deficit Int. J. Mol. Sci. 2021, 22, 1226 13 of 24 specifically reduces the production of aliphatic GS. Here, we observed an upregulation of SDI1 (Sulfur Deficit Induced 1) and LSU2 (Response to Low Sulfur). SDI1 has been reported to promote the release of sulfur from internal storage sites [90]. LSU2 is upregulated not only under S-deficiency, but also by Fe deficiency, high pH, and Cu toxicity [91]. Interestingly, BGLU28 was downregulated. This myrosinase gene usually is upregulated under sulfur deficiency, while bglu28 mutants display high GS levels [92]. Accordingly, our data indicate that in the tolerant deme A1(c+) bic does not directly cause S deficiency but specifically induces GS production and, consequently, the upregulation of S-deficiency related genes Int. J. Mol. Sci. 2021, 22, x FOR PEER REVIEW 14 of to 26 cope with the higher S-demand. Figure 5. Protein–protein Figure 5. interactionnetwork Protein–protein interaction networkfunctional functionalenrichment enrichmentanalysis analysisofofspecific specificgenes genes derived derived from from thethe response response to to bicarbonate of A1 (c+). Gene protein interaction network of A1(c+) exclusive DEGs from bic vs. pH 5.9 comparison after 3 bicarbonate of A1 (c+) . Gene protein interaction network of A1(c+) exclusive DEGs from bic vs. pH 5.9 comparison after 3 h h exposure. Each sphere corresponds to one gene and nodes represent protein interactions. Gene pathways are shown in exposure. Each sphere corresponds to one gene and nodes represent protein interactions. Gene pathways are shown in different colors. Arrows indicate up- or downregulation of the genes. different colors. Arrows indicate up- or downregulation of the genes. Shen et al. (2016) [48] reported that Fe deficiency induces SA accumulation, which enhances auxin and ethylene signaling and activates Fe transport by bHLH38/39 regula- tion of downstream iron genes. In fact, in the sensitive deme T6(C-), both the FIT-dependent [85] and the FIT-independent pathways [86] for iron homeostasis were activated by bic, as indicated by the upregulation of bHLH39 and bHLH100, respectively. Upregulation of NAS4 coding for nicotianamine synthase and NRAMP4 coding for Fe transporter export- ing Fe from vacuoles further confirms the quick activation of Fe mobilizing mechanisms in the bic-sensitive deme. In the tolerant A1(c+), bic enhanced SA, as shown by the upregulation of SARD1 [54] Int. J. Mol. Sci. 2021, 22, 1226 14 of 24 Int. J. Mol. Sci. 2021, 22, x FOR PEER REVIEW 16 of 26 FigureFigure 6. Protein–protein interaction network functional enrichment analysis of specific genes derived from the response 6. Protein–protein interaction network functional enrichment analysis of specific genes derived from the response to to bicarbonate of T6 (c−). Gene protein interaction network of T6(c−) exclusive DEGs from bic vs. pH 5.9 comparison after 3 bicarbonate of T6 (c−) sphere protein h exposure. Each . Gene interaction corresponds network to one gene of T6represent and nodes (c−) exclusive DEGs protein from bicGene interactions. vs. pH 5.9 comparison pathways are shown after in 3h exposure. Each sphere corresponds to one gene and nodes represent different colors. Arrows indicate up or downregulation of the genes. protein interactions. Gene pathways are shown in different colors. Arrows indicate up or downregulation of the genes. 2.7. Mechanisms Underlying Differential bic Toxicity and Tolerance in A. thaliana Demes STRING connected Our results GS genes demonstrate that to thecell cycleofand effects DNA-repair bicarbonate on A.related thalianagenes differsby AK3 and from WEE1 those(Figure caused 5). AK3 pH by high codes alone for(Figure an aspartate kinase 1). At equal pH, catalyzing the first bic is more toxic step inT6 to sensitive aspartate- (c−) derived amino and Col-0 thanacids [93]. WEE1 an alkaline solution is awith kinase involved organic buffer.inIncell thecycle inhibition tolerant A1(c+), therelated opposite to DNA was found. damage Induction and repair [94].ofSurprisingly, Fe deficiencyexposure is amongtothebicbest-known consequences downregulated many genesof bicar- involved bonaterepair, in DNA toxicitysuch in sensitive as XR1,plants. BUB32, In BRAC1, fact, the quick upregulation and MND1. DNAof Fe homeostasis-re- damage repair is usually lated genes in leaves of T6 (c−) exposed to bic (Figure 6) is a clear sign of the difficulty of upregulated in plants exposed to different abiotic and biotic stresses [95]. Downregulation these sensitive of WEE1 togetherplants to maintain sufficient with enhancement Fe transport of FIB coding to the leaves. for tryptophan Bicarbonate not a key aminotransferase, only interferes with the Fe-reducing capacity of the roots [101], but also specifically inhib- enzyme in auxin biosynthesis [96], supports the view that leaf cell growth was activated its the translocation from roots to shoots [102]. This inhibition of Fe translocation has clas- after 3 h exposure to bic. These results further indicate that in the tolerant A1(c+) adaptive sically been related to enhanced root production of organic acids, especially citrate. Citrate shoot growth events are activated in response to bic. However, RNA-seq only provides a snapshot of gene expression at a given time point; here, 3 h. It is possible that this growth promotion is a consequence of an earlier cell cycle arrest and upregulation of DNA repair in response to stress, followed by a downregulation after successful correction. To prove this attractive hypothesis, one must perform expression analysis after shorter exposure Int. J. Mol. Sci. 2021, 22, 1226 15 of 24 time. In addition, STRING connected GS biosynthesis to carbohydrate metabolism via adenosine 50 -phosphosulfate kinase (APK) and SEX1, a α-glucan water dikinase required for starch degradation. Enhanced starch degradation is also supported by the upregulation of BAM6 coding for ß-amylase, while BMY3 coding for α-amylase was downregulated. Glucose has been identified as a positive regulator of GS biosynthesis [97]. Furthermore, GS biosynthesis was related to organic carbon from 2-oxocarboxylic acid metabolism via AK3. DTX/MATE efflux carrier expression was downregulated by bic exposure. STRING connected downregulated DTX1 and DTX3 (At2g04050) to upregulated DHAR1 (deshy- droxyascorbate reductase) by means of the glucosyltransferase gene UGT74E2. Ascorbic acid not only plays a key role in the protection against oxidative stress, but also improves internal Fe availability [98]. The functions of DTX efflux carriers are still poorly established. DTX1 has been involved in cadmium tolerance [99]. DTX3 may facilitate membrane trans- port of xenobiotics; both DTX3 and DTX4 (At2g04070) seem to have an efflux protein with broad substrate specificity [100]. Downregulation of these DTXs could be related to the reduced efflux of some organic compound acting as a strong ligand for Fe, thus facilitating Fe availability. 2.7. Mechanisms Underlying Differential bic Toxicity and Tolerance in A. thaliana Demes Our results demonstrate that the effects of bicarbonate on A. thaliana differs from those caused by high pH alone (Figure 1). At equal pH, bic is more toxic to sensitive T6(c−) and Col-0 than an alkaline solution with organic buffer. In the tolerant A1(c+) , the opposite was found. Induction of Fe deficiency is among the best-known consequences of bicarbonate toxicity in sensitive plants. In fact, the quick upregulation of Fe homeostasis-related genes in leaves of T6(c−) exposed to bic (Figure 6) is a clear sign of the difficulty of these sensitive plants to maintain sufficient Fe transport to the leaves. Bicarbonate not only interferes with the Fe-reducing capacity of the roots [101], but also specifically inhibits the translocation from roots to shoots [102]. This inhibition of Fe translocation has classically been related to enhanced root production of organic acids, especially citrate. Citrate is a strong ligand for Fe and may favor Fe sequestration into root vacuoles [103]. More recently, studies with bic-exposed kiwi plantlets supplemented with 57 Fe related bic-induced inhibition of Fe translocation to the leaves with apoplastic accumulation of water-soluble Fe in the cell wall apoplast and the imbalance of nitrogen and carbon metabolism [98]. In fact, due to carbonic anhydrase activity, HCO3− supply favors dark fixation of inorganic C [104], leading to enhanced organic acid biosynthesis. In T6(c−) , bic exposure led to DEGs involved in nutrient balance and transport activities that were not directly related to Fe deficiency DEGs (Figure 6). Upregulation of BGAL2, coding for β-galactosidase, and downregulation of EXT3 and EXT21 (AT2G43150) suggests fast bic influence on leaf cell walls. Upregulation of AT5G13400 coding for the nitrate transporter NRT1, as well as downregulation of AAP4, responsible for phloem export of amino acids, is in line with enhanced nitrogen requirement in bic-stressed leaves. In the bicarbonate-tolerant A1(c+) , we have previously shown enhanced root exudation of flavonol-type phenolics enabling roots to mobilize sparingly soluble Fe [23]. Here, we further observed bic-induced upregulation of PYE and Ascorbate Reductase in leaves, both being factors that can contribute to increased availability of biologically active Fe in the leaves. To avoid imbalance of organic carbon caused by dark fixation of inorganic carbon derived from bic, A1(c+) may shuttle organic C into the pathway of aliphatic glucosino- lates, as indicated by the upregulation of aspartate kinase 3 (AK3) and 3-isopropylmalate dehydrogenase (IMD1) (Figure 6). Sulfur requirement for GS biosynthesis is provided by the upregulation of sulfur deficiency-induced genes such as LSU2 and SDI1 and the downregulation of BGLU28 involved in GS degradation under S deficiency conditions. A key question concerns the trigger mechanism and the regulation of these quick differential leaf responses in A1(c+) . According to our RNA-seq data, SA and JA may play a central role. Bic-induced upregulation of WRKY70 connected the enhanced SA pathway to differential regulation of the JA pathway and GS biosynthesis (Figure 5). Interestingly, in Int. J. Mol. Sci. 2021, 22, 1226 16 of 24 the tolerant A1(c+) , bic exposure specifically enhanced the expression of certain defense- related RLKs and RLPs (Figure 4). A primary objective of future research is to explore whether apoplastic bic, less Fe translocation from root to shoot, or any other root-derived signal is responsible for the activation of these primary signal perception molecules that trigger defense mechanisms in the leaves of bic-tolerant A1(c+) . 3. Materials and Methods 3.1. Plant Material In previous studies, natural variation of Arabidopsis thaliana populations from Cat- alonia [105] were tested in a multi-year small scale common garden under carbonated condition [23]. Seeds of extreme phenotype lines A1(c+) , a moderately alkaline-tolerant deme, and T6(c−) , an alkaline-sensitive deme, were collected from the last reciprocal trans- plant experiment performed in 2015 and stored under fresh (4 ◦ C) and dry conditions until the beginning of the experiments. Col-0 seeds purchased from The Nottingham Arabidop- sis Stock Centre (NASC, Nottingham, U.K.) [106] were included as a reference genome. 3.2. Growth Conditions Seed sterilization: Selected seeds were surface sterilized by soaking in 70% (v/v) ethanol for 1 min, suspended in 30% (v/v) commercial Clorox bleach and 1 drop of Tween- 20 for 5 min, and rinsed 5 times in sterile 18 MΩmilli-Q water. Hydroponic culture: Seeds were sown in 0.2 mL tubes containing 0.6% agar prepared in nutrient solution 1/2 Hoagland (pH 5.9). Seeds were kept at 4 ◦ C for 7 days in the dark to synchronize germination. Tubes containing seed were placed in the growth chamber (12 h light/12 h dark, 150 µmol cm−2 ·s−1 , 40% humidity and 25 ◦ C). After root emergence, the bottom of the tubes containing seedlings was cut off and the tubes were placed in 150 mL hydroponic containers with aerated nutrient solution 1/2 Hoagland (pH 5.9). When 15 days old, the seedlings were separated in different sets and the following treatments were applied: control ( 12 Hoagland solution at pH 5.9), high pH ( 12 Hoagland solution at pH 8.3), and bicarbonate ( 12 Hoagland solution at pH 8.3 with 10 mM NaHCO3 ). Solutions were buffered with different proportions of MES (2-(N-morpholino) ethanesul- fonic acid hydrate, 4-morpholineethanesulfonic acid) and BTP, depending on final pH in continuous aeration. The selected bicarbonate concentration has previously been shown to discriminate between the carbonate-tolerant and -sensitive deme [23]. To check whether this concentration allows reproductive growth in the tolerant accession, we performed A1(c+) preliminary essays with bicarbonate concentrations ranging from 0 to 20 mM (pH 8.3). Rosette diameter and length of the longest root were monitored every week by image analysis (Figure S1). Plate culture: For germination assays, 45 seeds from each deme were sown in plates under a flow cabinet with sterile material. Plates contained 3 treatments: control ( 12 MS at pH 5.9), high pH ( 12 MS at pH 8.3), and bicarbonate ( 12 MS at pH 8.3 and 10 mM NaHCO3 ). All plates contained Phyto-agar 0.6% (Duchefa, Haarlem, The Netherlands), and solutions were buffered using different proportions of MES and BTP depending on final pH. Plates with seeds were kept at 4 ◦ C for synchronizing germination. After 7 days under stratification treatment, plates were moved to a growth chamber (12 h light/12 h dark, 150 µmol cm−2 ·s−1 , 40% humidity and 25 ◦ C). Germination and radicle emergence were daily checked during the following 10 days. 3.3. Physiological Trait Measurements Plants from tolerant A1(c+) , sensitive T6(c−) , and Col-0 were used to test fitness in hydroponic conditions. Every week, pictures of the entire plants were taken. The pho- tographs were used to measure the length of the largest root and the rosette diameter (n = 9–12) using ImageJ software [107]. We did not use a dye prior to image acquisition. Siliques production was counted at maturity (n = 9). The leaf nutritional state of plants was assessed in plants submitted for 15 days to treatment conditions. Leaves were washed Int. J. Mol. Sci. 2021, 22, 1226 17 of 24 with 18 MΩ water before placing into Pyrex digestion tubes. Sampled plant material was dried for 4 days at 60 ◦ C and weighed before open-air digestion in Pyrex tubes using 0.7 mL concentrated HNO3 at 110 ◦ C for 5 h in a hot-block digestion system (SC154-54-Well Hot Block, Environmental Express, Charleston, SC, USA). Four plant replicates from each deme and treatment were used to determine the concentrations of the selected elements (Ca, K, Mg, Na, P, S, B, Mo, Cu, Fe, Mn, Zn) by inductively coupled plasma optical emission spectroscopy (ICP-OES; Thermo Jarrell-Ash, model 61E Polyscan, Waltham, MA, USA). 3.4. RNA Isolation, Library Construction, and Sequencing Leaf material from plants cultivated under hydroponic conditions were recollected 3 and 48 h after starting the treatments. Twelve plants per line and treatment were pooled to perform 3 biological replicates. Leaves were immersed in liquid nitrogen, homogenized to a fine powder, and stored at –80 ◦ C. The total RNA of 100 mg of leaf powder for each biological replicate was extracted using the Maxwell plant RNA kit (Promega Corporation, Madison, WI, USA) following the manufacturer’s instructions. For each sample, approxi- mately 2 µg of total RNA was used for quality evaluation. Total RNA was quantified using a Qubit 2.0 Q32866 (Life Technologies, Carlsbad, CA, USA) to prepare the complementary deoxyribonucleic acid (cDNA) library. cDNA library was composed by 72 samples, collected from A1(c+ ) and T6(c− ) plants exposed for 3 h (quick response) or 48 h (long response) to the different treatments (2 bi- ological replicates and 3 technical replicates). Library preparation was performed using Novogene protocol. The cDNA library was sequenced at the Illumina NovaSeq 6000 using standard procedures to generate paired end reads of 150 bp and number of analyzed read 30-47 million pairs (Novogene, Sacramento, CA, USA). Clean reads were obtained by removing reads containing adapters, reads containing poly-N, and low-quality reads from raw data. Simultaneously, Q20, Q30, and GC content of the clean data were calculated. After filtering, clean reads on average had 35.7 million reads, which occupied 98% of the total reads of all libraries. In all samples, the Q30 value was over 89%, and the GC content was 44% (Datasets S2 and S3). All downstream analyses were based on the clean data with high-quality reads. Reference genome and gene model annotation files were downloaded from the genome website directly. Index of the reference genome was built using Bowtie v2.2.3 (John Hopkins University, Baltimore, MD, USA) and paired-end clean reads were aligned to the reference genome using TopHat v2.0.12 (John Hopkins University, Baltimore, MD, USA). Differential expression analysis of samples was performed using a model based on the negative binomial distribution [108] and the DESeq R package (1.18.0) (www.Bioconductor. org). The resulting p-values were adjusted using the Benjamini and Hochberg’s ap- proach [109] for controlling the false discovery rate (FDR). To perform the final list of DEGs, we filtered genes by adjusted p-value < 0.05 and LFC > 1 and LFC < −1. Raw reads and normalized gene expression of RNA-seq analysis were deposited in the Gene Expres- sion Omnibus (GEO) database and in the National Center for Biotechnology Information’s (NCBI) Sequence Read Archive (GSE164502). 3.5. Gene Ontology, KEGG Pathway, and Functional Protein Association Network Analysis of DEGs Gene Ontology (GO) enrichment analysis of DEGs was implemented by AgriGO V2 (GO Analysis Toolkit and Database for Agricultural Community (cau.edu.cn) [110]. Signifi- cant GO terms (p < 0.05) were classified into 3 categories: biological function, molecular process, and cellular component. KEGG pathway and STRING version 11.0 were used to understand high-level function and gene interaction network of differential expressed genes [111]. 3.6. Relative Expression Analysis To produce cDNAs from RNAs, we used the iScriptTM cDNA Synthesis Kit (Bio-Rad, USA) with 1 µL iScript Reverse Transcriptase + 4 µL 5x iScript Reaction Mix + Sample + Int. J. Mol. Sci. 2021, 22, 1226 18 of 24 Molecular Water to obtain 20 µL volume. Samples were run in a thermocycler (48-well MJ MiniTM, Bio-rad, Hercules, CA, USA) at 5 min 25 ◦ C, 30 s 42 ◦ C, and 5 s 85 ◦ C. Dilution of the cDNAs was performed 1/50 with water (Molecular Biology Reagent, Sigma-Aldrich, St. Louis, MO, USA). Diluted cDNA (1:50) was used as a template for quantitative PCRs using iTaq Universal SYBR Green Supermix (Bio-Rad, Hercules, CA, USA). Real-time detection of fluorescence emission was performed on a CFX384 Real-Time System (Bio-Rad, Hercules, CA, USA) using the following conditions: denaturalization step 1000 95 ◦ C followed by annealing and extension 3000 60 ◦ C. A total of 40 cycles were run. Melt curve was performed, increasing from 65.0 ◦ C to 95.0 ◦ C by 0.5 ◦ C each 5 seconds Plates were edited using the CFX manager version 3.1 software. Extremely expressed genes across lines at 3 h were chosen as candidate genes to validate transcriptomic experiment using qPCR. Primers from selected genes were designed using NCBI primer blast tool [112] (Biole- gio, Nijmegen, The Netherlands). The sequences of primers used are detailed in Figure S2. The expression of target genes was normalized to the expression level of the Actin and Tubulin genes of A. thaliana [113]. The relative expression (RE) of each gene was calculated in comparison to the control treatments (pH 5.9) at each time point. The expression of the target gene relative to the expression of the reference gene was calculated using the 2−∆∆Ct method [114]. 3.7. Statistical Analysis For the physiological measurements, normal distribution of the data was confirmed by Levene’s test. Differences between treatments and different demes were analyzed by two-way ANOVA. Post hoc analyses were realized using Tukey’s test. Nutrient mineral data were standardized to avoid bias due to different magnitude orders. Radial plots were used to visualize differences among treatments in each deme. Statistical analyses and plots were performed using JMP13 software (JMP, Version 13. SAS Institute Inc., Cary, NC, USA, 1989–2019). Supplementary Materials: The following are available online at https://www.mdpi.com/1422-006 7/22/3/1226/s1. Figure S1: Growth responses of moderately bicarbonate-tolerant A. thaliana deme to different NaHCO3 concentrations. (A) Rosette diameter (cm) and (B) root length (cm) of A1(c+) plants grown in hydroponics for 5 weeks. Letters indicate significant differences between treatment (n = 9; p < 0.05 according to Tukey’s HSD). Figure S2: RNA-seq validation. (A) Primer sequence from eight selected genes used for validation of the RNA-seq expression results. (B) Bar plot from RNA-seq (dark grey) and qPCR (light grey) expression results from the selected genes after 3 h of exposure to bic and pH 5.9. Results are expressed in log2 fold change; bars indicate standard error. Figure S3: GO and KEGG pathways of demes comparing alkaline treatment at two time points. Venn diagram of up- and downregulated DEGs comparing bic vs. pH 5.9 and pH 8.3 vs. pH 5.9 treatments after 3 h and 48 h between A1(c+) (A) and T6(c−) (B) demes. (C) Bubble plot indicating significant GO analysis of DEGs in bic vs. pH 5.9 and pH 8.3 vs. pH 5.9 treatments after 3 h and 48 h in A1(c+) and T6(c−) demes. GO were filtered to adjusted p-value < 0.05. Scale colors indicate number of DEGs, while bubble size indicates -log of adjusted p-value. GO terms were separated into biological function, cellular component, and molecular function. (D) Heatmap plots of KEGG pathway analysis of DEGs in bic vs. pH 5.9 and pH 8.3 vs. pH 5.9 treatments after 3 h and 48 h between A1(c+) and T6(c−) demes. KEGG pathway terms were filtered by p-value < 0.05. Scale colors indicate pathway fold enrichment. Supplementary Dataset: Dataset S1: Statistical analysis of physiologic experiments (Figure 1, Figure S1). Dataset S2: Sequencing quality control. Dataset S3: Clean reads mapped into 24,400 genes. Pairwise comparation of alkaline treatments vs. control without filtering. Dataset S4: Slim GO from DEGs (LFC> |1| and adjusted p-value < 0.05) in pairwise comparison treatments bic vs. pH 5.9 and pH 8.3 vs. pH 5.9 at different time points (3 h and 48 h) from A1(c+) and T6(c−) demes. Dataset S5: Complete GO from DEGs (LCF> |1| and adjusted p-value < 0.05) in pairwise comparison treatments bic vs. pH 5.9 and pH 8.3 vs. pH 5.9 at different time points (3 h and 48 h) from A1(c+) and T6(c−) demes. Dataset S6: KEGG pathway of DEGs (LFC> |1| and p-value < 0.05) in pairwise comparison treatments bic vs. pH 5.9 and pH 8.3 vs. pH 5.9 at different time points (3 h and 48 h) from A1(c+) and T6(c−) demes. Dataset S7: Slim GO of upregulated and downregulated Int. J. Mol. Sci. 2021, 22, 1226 19 of 24 DEGs (LFC > |1| and adjusted p-value < 0.05) in pairwise comparison bic vs. pH 5.9, from A1(c+) and T6(c−) at 3 h and 48 h. Dataset S8: Complete GO of upregulated and downregulated DEGs (LFC > |1| and adjusted p-value < 0.05) in pairwise comparison bic vs. pH 5.9, from A1(c+) and T6(c−) at 3 h and 48 h. Dataset S9: KEGG pathway of upregulated and downregulated DEGs (LFC> |1| and p-value < 0.05) in pairwise comparison bic vs. pH 5.9, from A1(c+) and T6(c−) at 3h and 48h. Dataset S10: Top differentially expressed genes at 3 h in A1(c+) and T6(c−) . Dataset S11: Protein–protein interaction network functional enrichment analysis output from exclusive genes involved in bic response at 3 h in A1(c+) deme. DEGs filtered (LFC > |1| and adjusted p-value < 0.05). Dataset S12: GO and KEGG pathway terms of exclusive genes involved in bic response at 3 h in A1(c+) deme. DEGs filtered (LFC > |1| and adjusted p-value < 0.05). Dataset S13: Protein–protein interaction network functional enrichment analysis output of exclusive genes involved in bic 8.3 response at 3 h in T6(c−) deme. DEGs filtered (LFC > |1| and adjusted p-value < 0.05). Dataset S14: GO and KEGG pathway terms of exclusive genes involved in bic 8.3 response at 3 h in T6(c−) deme. DEGs filtered (LFC > |1| and adjusted p-value < 0.05). Author Contributions: Conceptualization, L.P.-M., S.B., C.P. and R.T.; methodology, L.P.-M., S.B., C.P. and R.T.; investigation, L.P.-M.; writing—original draft preparation, reviewing, and editing, L.P.-M., S.B. and C.P.; supervision, C.P. and R.T.; project funding acquisition, C.P. and R.T. All authors have read and agreed to the published version of the manuscript. Funding: This research was funded by Spanish MICINN, grant number Sal-CAL-MED (BFU2016- 75176-R) and IONPLAMIC (PID2019-104000RB-100) participated in by L.P., S.B., R.T. and C.P. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Not applicable. Conflicts of Interest: The authors declare no conflict of interest. Abbreviations ABA Abscisic acid ANOVA Analysis of variance A. thaliana Arabidopsis thaliana cDNA Complementary deoxyribonucleic acid BTP 1,3-Bis[tris(hydroxymethyl)methylamino]propane DEG Differentially expressed genes FDR Fold discovery rate GO Gene Ontology GS Glucosinolate GST Gluthatione-S-transferase ICP-OES Inductively coupled plasma atomic emission spectroscopy JA Jasmonic acid KEGG Kyoto Encyclopedia of Genes and Genomes LFC Log2 fold change MATE Multidrug and Toxic Compound Extrusion MES 2-(N-Morpholino)ethanesulfonic acid hydrate, 4-morpholineethanesulfonic acid NADPH Nicotinamide adenine dinucleotide phosphate NCBI National Center for Biotechnology Information Adj p-value Adjusted p-value RLK Receptor-like kinase RLP Receptor-like protein RNA Ribonucleic acid ROS Reactive oxygen species qPCR Quantitative Polymerase Chain Reaction TF Transcription factor SA Salicylic acid Int. J. Mol. Sci. 2021, 22, 1226 20 of 24 References 1. Al-Busaidi, A.S.; Cookson, P. Salinity–pH Relationships in Calcareous Soils. J. Agric. Mar. Sci. 2003, 8, 41. [CrossRef] 2. Loeppert, R.H.; Suarez, D.L. Carbonate and Gypsum. In Methods of Soil Analysis. Part 3. Chemical Methods; Sparks, D.L., Page, A.L., Hemke, P.A., Loeppert, R.H., Soltanpour, M.A., Tabatabai, M.A., Johnston, C.T., Summer, M.E., Eds.; Soil Science Society of America: Madison, WI, USA; American Society of Agronomy: Lincoln, NE, USA, 1996; pp. 437–474. ISBN 9780891188254. 3. Hu, H.; Rappel, W.-J.; Occhipinti, R.; Ries, A.; Böhmer, M.; You, L.; Xiao, C.; Engineer, C.B.; Boron, W.F.; Schroeder, J.I. Distinct Cellular Locations of Carbonic Anhydrases Mediate Carbon Dioxide Control of Stomatal Movements. Plant Physiol. 2015, 169, 1168–1178. [CrossRef] [PubMed] 4. Aboukila, E.F.; Nassar, I.N.; Rashad, M.; Hafez, M.; Norton, J.B. Reclamation of calcareous soil and improvement of squash growth using brewers’ spent grain and compost. J. Saudi Soc. Agric. Sci. 2018, 17, 390–397. [CrossRef] 5. Morrissey, J.; Guerinot, M.L. Iron Uptake and Transport in Plants: The Good, the Bad, and the Ionome. Chem. Rev. 2009, 109, 4553–4567. [CrossRef] 6. Kobayashi, T. Understanding the Complexity of Iron Sensing and Signaling Cascades in Plants. Plant Cell Physiol. 2019, 60, 1440–1446. [CrossRef] 7. Zhang, X.; Zhang, D.; Sun, W.; Wang, T. The Adaptive Mechanism of Plants to Iron Deficiency via Iron Uptake, Transport, and Homeostasis. Int. J. Mol. Sci. 2019, 20, 2424. [CrossRef] 8. López-Millán, A.-F.; Grusak, M.A.; Abadía, A.; Abadía, J. Iron deficiency in plants: An insight from proteomic approaches. Front. Plant Sci. 2013, 4, 254. [CrossRef] 9. Donnini, S.; De Nisi, P.; Gabotti, D.; Tato, L.; Zocchi, G. Adaptive strategies of Parietaria diffusa (M.&K.) to calcareous habitat with limited iron availability. Plant Cell Environ. 2012, 35, 1171–1184. [CrossRef] 10. Higuchi, K.; Ono, K.; Araki, S.; Nakamura, S.; Uesugi, T.; Makishima, T.; Ikari, A.; Hanaoka, T.; Sue, M. Elongation of barley roots in high-pH nutrient solution is supported by both cell proliferation and differentiation in the root apex. Plant Cell Environ. 2017, 40, 1609–1617. [CrossRef] 11. Li, J.; Xu, H.-H.; Liu, W.-C.; Zhang, X.-W.; Lu, Y.-T. Ethylene Inhibits Root Elongation during Alkaline Stress through AUXIN1 and Associated Changes in Auxin Accumulation. Plant Physiol. 2015, 168, 1777–1791. [CrossRef] 12. Zhu, D.; Bai, X.; Chen, C.; Chen, Q.; Cai, H.; Li, Y.; Ji, W.; Zhai, H.; Lv, D.; Luo, X.; et al. GsTIFY10, a novel positive regulator of plant tolerance to bicarbonate stress and a repressor of jasmonate signaling. Plant Mol. Biol. 2011, 77, 285–297. [CrossRef] 13. Lucena, C.; Romera, F.J.; García, M.J.; Alcántara, E.; Pérez-Vicente, R. Ethylene participates in the regulation of Fe deficiency responses in strategy I plants and in rice. Front. Plant Sci. 2015, 6, 1–16. [CrossRef] 14. Ahmad, P.; Öztürk, M.; Sharma, S.; Gucel, S. Effect of sodium carbonate-induced salinity–alkalinity on some key osmoprotectants, protein profile, antioxidant enzymes, and lipid peroxidation in two mulberry (Morus alba L.) cultivars. J. Plant Interact. 2014, 9, 460–467. [CrossRef] 15. Zhang, J.; Wang, J.; Jiang, W.; Liu, J.; Yang, S.; Gai, J.; Li, Y. Identification and Analysis of NaHCO3 Stress Responsive Genes in Wild Soybean (Glycine soja) Roots by RNA-seq. Front. Plant Sci. 2016, 7, 1842. [CrossRef] 16. Duan, X.; Yu, Y.; Zhang, Y.; Chen, C.; Duanmu, H.; Cao, L.; Sun, M.; Sun, X.; Zhu, Y. A potential efflux boron transporter gene GsBOR2, positively regulates Arabidopsis bicarbonate tolerance. Plant Sci. 2018, 274, 284–292. [CrossRef] 17. Ge, Y.; Li, Y.; Zhu, Y.-M.; Bai, X.; Lv, D.; Guo, D.; Ji, W.; Cai, H. Global transcriptome profiling of wild soybean (Glycine soja) roots under NaHCO3 treatment. BMC Plant Biol. 2010, 10, 153. [CrossRef] 18. Guo, M.; Li, S.; Tian, S.; Wang, B.; Zhao, X. Transcriptome analysis of genes involved in defense against alkaline stress in roots of wild jujube (Ziziphus acidojujuba). PLoS ONE. 2017, 12, 0185732. [CrossRef] 19. Wang, C.; Gao, C.; Wang, L.; Zheng, L.; Yang, C.; Wang, Y. Comprehensive transcriptional profiling of NaHCO3 -stressed Tamarix hispida roots reveals networks of responsive genes. Plant Mol. Biol. 2004, 84, 145–157. [CrossRef] 20. Gong, B.; Zhang, C.; Li, X.; Wen, D.; Wang, S.; Shi, Q.; Wang, X. Identification of NaCl and NaHCO3 stress responsive proteins in tomato roots using iTRAQ-based analysis. Biochem. Biophys Res. Commun. 2014, 446, 417–422. [CrossRef] 21. Waters, B.M.; Amundsen, K.; Graef, G. Gene Expression Profiling of Iron Deficiency Chlorosis Sensitive and Tolerant Soybean Indicates Key Roles for Phenylpropanoids under Alkalinity Stress. Front. Plant Sci. 2018, 9. [CrossRef] 22. Zhang, H.; Liu, X.L.; Zhang, R.X.; Yuan, H.Y.; Wang, M.M.; Yang, H.Y.; Ma, H.Y.; Liu, D.; Jieng, C.J.; Liang, Z.W. Root damage under alkaline stress is associated with reactive oxygen species accumulation in rice (Oryza sativa L.). Front. Plant Sci. 2017, 8, 1580. [CrossRef] 23. Terés, J.; Busoms, S.; Pérez-Martín, L.; Luís-Villaroya, A.; Flis, P.; Álvarez-Fernández, A.; Tolrà, R.; Salt, D.E.; Poschenrieder, C. Soil carbonate drives local adaptation in Arabidopsis thaliana. Plant Cell Environ. 2019, 42, 2384–2398. [CrossRef] 24. Chen, L.; Wang, G.; Chen, P.; Zhu, H.; Wang, S.; Ding, Y. Shoot-Root Communication Plays a Key Role in Physiological Alterations of Rice (Oryza sativa) Under Iron Deficiency. Front. Plant Sci. 2018, 9, 757. [CrossRef] 25. García, M.J.; Corpas, F.J.; Lucena, C.; Alcántara, E.; Pérez-Vicente, R.; Zamarreño, A.M.; Bacaicoa, E.; García-Mina, J.M.; Bauer, P.; Romera, F.J. A Shoot Fe Signaling Pathway Requiring the OPT3 Transporter Controls GSNO Reductase and Ethylene in Arabidopsis thaliana Roots. Front. Plant Sci. 2018, 9. [CrossRef] 26. Tudela-Isanta, M.; Ladouceur, E.; Wijayasinghe, M.; Pritchard, H.W.; Mondoni, A. The seed germination niche limits the distribution of some plant species in calcareous or siliceous alpine bedrocks. Alp. Bot. 2018, 128, 83–95. [CrossRef] Int. J. Mol. Sci. 2021, 22, 1226 21 of 24 27. Lin, J.; Li, X.; Zhang, Z.; Yu, X.; Gao, Z.; Wang, Y.; Wang, J.; Li, Z.; Mu, C. Salinity-alkalinity tolerance in wheat: Seed germination, early seedling growth, ion relations and solute accumulation. African J. Agric. Res. 2012, 7, 467–474. [CrossRef] 28. Meng, C.; Quan, T.-Y.; Li, Z.-Y.; Cui, K.L.; Liang, Y.; Dai, J.L.; Xia, G.M.; Liu, S.W. Transcriptome profiling reveals the genetic basis of alkalinity tolerance in wheat. BMC Genomics. 2017, 18, 24. [CrossRef] 29. Hsieh, E.J.; Waters, B.M. Alkaline stress and iron deficiency regulate iron uptake and riboflavin synthesis gene expression differently in root and leaf tissue: Implications for iron deficiency chlorosis. J. Exp. Bot. 2016, 67, 5671–5685. [CrossRef] 30. Bergmann, W. Ernährungsstörungen bei Kulturpflanzen, 3rd ed.; Gustav Fischer Verlag: Jena, Germany, 1993; ISBN 13: 978-3-334- 60414-4. 31. Sisó-Terraza, P.; Luis-Villarroya, A.; Fourcroy, P.; Briat, J.F.; Abadía, A.; Gaymard, F.; Abadía, J.; Álvarez-Fernández, A. Accumula- tion and Secretion of Coumarinolignans and other Coumarins in Arabidopsis thaliana Roots in Response to Iron Deficiency at High pH. Front. Plant Sci. 2016, 7. [CrossRef] 32. Ma, H.; Yang, H.; Lü, X.; Wu, P.H.; Liang, Z.; Ooi, M.K.J. Does high pH give a reliable assessment of the effect of alkaline soil on seed germination? A case study with Leymus chinensis (Poaceae). Plant Soil. 2015, 394, 35–43. [CrossRef] 33. Lan, P.; Li, W.; Schmidt, W. A Digital Compendium of Genes Mediating the Reversible Phosphorylation of Proteins in Fe-Deficient Arabidopsis Roots. Front. Plant Sci. 2013, 4, 173. [CrossRef] 34. Steidele, C.E.; Stam, R. Multi-omics approach highlights differences between functional RLP classes in Arabidopsis thaliana. Biorxiv 2020. [CrossRef] 35. Sakai, T.; Van Der Honing, H.; Nishioka, M.; Uehara, Y.; Takahashi, M.; Fujisawa, N.; Saji, N.; Seki, M.; Shinozaki, K.; Jones, M.A.; et al. Armadillo repeat-containing kinesins and a NIMA-related kinase are required for epidermal-cell morphogenesis in Arabidopsis. Plant J. 2008, 53, 157–171. [CrossRef] [PubMed] 36. Tarutani, Y.; Sasaki, A.; Yasuda, M.; Nakashita, H.; Yoshida, S.; Yamaguchi, I.; Suzuki, Y. Identification of Three Clones Which Commonly Interact with the Kinase Domains of Highly Homologous Two Receptor-Like Kinases, RLK902 and RKL1. Biosci. Biotechnol. Biochem. 2004, 68, 2581–2587. [CrossRef] [PubMed] 37. Yamamoto, H.; Fan, X.; Sugimoto, K.; Fukao, Y.; Peng, L.; Shikanai, T. CHLORORESPIRATORY REDUCTION 9 is a Novel Factor Required for Formation of Subcomplex A of the Chloroplast NADH Dehydrogenase-Like Complex. Plant Cell Physiol. 2016, 57, 2122–2132. [CrossRef] [PubMed] 38. Msilini, N.; Oueslati, S.; Amdouni Chebbi, M.; Ksouri, R.; Lachaâl, M.; Ouerghi, Z. Variability of phenolic content and antioxidant activity of two lettuce varieties under Fe deficiency. J. Sci. Food Agric. 2013, 93, 2016–2021. [CrossRef] 39. Ranieri, A.; Castagna, A.; Baldan, B.; Soldatini, G.F. Iron deficiency differently affects peroxidase isoforms in sunflower. J. Exp. Bot. 2001, 52, 25–35. [CrossRef] 40. Tognolli, M.; Penel, C.; Greppin, H.; Simon, P. Analysis and expression of the class III peroxidase large gene family in Arabidopsis thaliana. Gene 2002, 288, 129–138. [CrossRef] 41. Tokunaga, N.; Kaneta, T.; Sato, S.; Sato, Y. Analysis of expression profiles of three peroxidase genes associated with lignification in Arabidopsis thaliana. Physiol. Plant 2009, 136, 237–249. [CrossRef] 42. Xu, C.; Cao, H.; Xu, E.; Zhang, S.; Hu, Y. Genome-Wide Identification of Arabidopsis LBD29 Target Genes Reveals the Molecular Events behind Auxin-Induced Cell Reprogramming during Callus Formation. Plant Cell Physiol. 2018, 59, 749–760. [CrossRef] 43. Fiorucci, A.-S.; Bourbousse, C.; Concia, L.; Rougée, M.; Deton-Cabanillas, A.-F.; Zabulon, G.; Layat, E.; Latrasse, D.; Kim, S.K.; Chaumont, N.; et al. Arabidopsis S2Lb links AtCOMPASS-like and SDG2 activity in H3K4me3 independently from histone H2B monoubiquitination. Genome Biol. 2019, 20, 100. [CrossRef] 44. Santos, C.S.; Ozgur, R.; Uzilday, B.; Turkan, I.; Roriz, M.; Rangel, A.; Carvalho, S.M.P.; Vasconcelos, M.W. Understanding the Role of the Antioxidant System and the Tetrapyrrole Cycle in Iron Deficiency Chlorosis. Plants 2019, 8, 348. [CrossRef] 45. Rodríguez-Milla, M.A.; Maurer, A.; Huete, A.R.; Gustafson, J.P. Glutathione peroxidase genes in Arabidopsis are ubiquitous and regulated by abiotic stresses through diverse signaling pathways. Plant J. 2003, 36602–36615. [CrossRef] 46. Gullner, G.; Komives, T.; Király, L.; Schröder, P. Glutathione S-transferase enzymes in plant-pathogen interactions. Front. Plant Sci. 2018, 871, 1836. [CrossRef] 47. Piślewska-Bednarek, M.; Nakano, R.T.; Hiruma, K.; Pastorczyk, M.; Sanchez-Vallet, A.; Singkravanit-Ogawa, S.; Ciesiolka, D.; Takano, Y.; Molina, P.; Schulze-Lefert, P.; et al. Glutathione Transferase U13 Functions in Pathogen-Triggered Glucosinolate Metabolism. Plant Physiol. 2018, 176, 538–551. [CrossRef] 48. Shen, C.; Yang, Y.; Liu, K.; Zhang, L.; Gou, H.; Sun, T.; Wang, H. Involvement of endogenous salicylic acid in iron-deficiency responses in Arabidopsis. J. Exp. Bot. 2016, 67, 4179–4193. [CrossRef] 49. Liu, T.; Zhang, X.; Yang, H.; Agerbirk, N.; Qiu, Y.; Wang, H.; Shen, D.; Song, J.; Li, X. Aromatic Glucosinolate Biosynthesis Pathway in Barbarea vulgaris and its Response to Plutella xylostella Infestation. Front. Plant Sci. 2016, 7, 1–15. [CrossRef] 50. Czerniawski, P.; Bednarek, P. Glutathione S-Transferases in the Biosynthesis of Sulfur-Containing Secondary Metabolites in Brassicaceae Plants. Front. Plant Sci. 2018, 9, 1639. [CrossRef] 51. Dixon, D.P.; Edwards, R. Selective Binding of Glutathione Conjugates of Fatty Acid Derivatives by Plant Glutathione Transferases. J. Biol. Chem. 2009, 284, 21249–21256. [CrossRef] 52. Gao, X.; Chen, X.; Lin, W.; Chen, S.; Lu, D.; Niu, Y.; Li, L.; Cheng, C.; McCormack, M.; Sheen, J.; et al. Bifurcation of Arabidopsis, N.L.R. Immune Signaling via Ca2+ -Dependent Protein Kinases. PLoS Pathog. 2013, 9, 1003127. [CrossRef] Int. J. Mol. Sci. 2021, 22, 1226 22 of 24 53. Scholz, S.S.; Vadassery, J.; Heyer, M.; Reichelt, M.; Snedden, W.A.; Boland, W.; Mithöfer, A. Mutation of the Arabidopsis Calmodulin-Like Protein CML37 Deregulates the Jasmonate Pathway and Enhances Susceptibility to Herbivory. Mol. Plant 2014, 7, 1712–1726. [CrossRef] 54. Zhang, Y.; Xu, S.; Ding, P.; Wang, D.; Cheng, Y.T.; He, J.; Gao, M.; Xu, F.; Li, Y.; Zhu, Z.; et al. Control of salicylic acid synthesis and systemic acquired resistance by two members of a plant-specific family of transcription factors. Proc. Natl. Acad. Sci. USA 2010, 107, 18220–18225. [CrossRef] 55. Xu, B.; Cheval, C.; Laohavisit, A.; Hocking, B.; Chiasson, D.; Olson, T.; Shirasu, K.; Faulkner, C.; Gilliham, M. A calmodulin-like protein regulates plasmodesmal closure during bacterial immune responses. New Phytol. 2017, 215, 77–84. [CrossRef] 56. Vernoux, T.; Brunoud, G.; Farcot, E.; Morin, V.; Van den Daele, H.; Legrand, J.; Oliva, M.; Das, P.; Larrieu, A.; Wells, D.; et al. The auxin signalling network translates dynamic input into robust patterning at the shoot apex. Mol. Syst. Biol. 2011, 7, 508. [CrossRef] 57. Ogawa, D.; Yamaguchi, K.; Nishiuchi, T. High-level overexpression of the Arabidopsis HsfA2 gene confers not only increased themotolerance but also salt/osmotic stress tolerance and enhanced callus growth. J. Exp. Bot. 2007, 58, 3373–3383. [CrossRef] 58. Lu, P.; Chai, M.; Yang, J.; Ning, G.; Wang, G.; Ma, H. The Arabidopsis CALLOSE DEFECTIVE MICROSPORE1 Gene Is Required for Male Fertility through Regulating Callose Metabolism during Microsporogenesis. Plant Physiol. 2014, 164, 1893–1904. [CrossRef] 59. Nakano, Y.; Yamaguchi, M.; Endo, H.; Rejab, N.A.; Ohtani, M. NAC-MYB-based transcriptional regulation of secondary cell wall biosynthesis in land plants. Front. Plant Sci. 2015, 6, 288. [CrossRef] 60. Ramirez-Parra, E.; López-Matas, M.A.; Fründt, C.; Gutierrez, C. Role of an atypical E2F transcription factor in the control of Arabidopsis cell growth and differentiation. Plant Cell. 2004, 16, 2350–2363. [CrossRef] 61. Kim, S.A.; LaCroix, I.S.; Gerber, S.A.; Guerinot, M.L. The iron deficiency response in Arabidopsis thaliana requires the phosphory- lated transcription factor URI. Proc. Natl. Acad. Sci. USA 2019, 116, 24933–24942. [CrossRef] 62. Kroh, G.E.; Pilon, M. Connecting the negatives and positives of plant iron homeostasis. New Phytol. 2019, 223, 1052–1055. [CrossRef] 63. Schwarz, B.; Bauer, P. FIT, a regulatory hub for iron deficiency and stress signaling in roots, and FIT-dependent and -independent gene signatures. J. Exp. Bot. 2020, 71, 1694–1705. [CrossRef] [PubMed] 64. Ikeda, M.; Fujiwara, S.; Mitsuda, N.; Ohme-Takagi, M. A triantagonistic basic helix-loop-helix system regulates cell elongation in Arabidopsis. Plant Cell. 2012, 24, 4483–4497. [CrossRef] [PubMed] 65. Schuster, C.; Gaillochet, C.; Medzihradszky, A.; Busch, W.; Daum, G.; Krebs, M.; Kehle, A.; Lohmann, J.U. A regulatory framework for shoot stem cell control integrating metabolic, transcriptional, and phytohormone signals. Dev. Cell 2014, 28, 438–449. [CrossRef] [PubMed] 66. Lau, O.S.; Davies, K.A.; Chang, J.; Adrian, J.; Rowe, M.H.; Ballenger, C.E.; Bergmann, D.C. Direct roles of SPEECHLESS in the specification of stomatal self-renewing cells. Science 2014, 345, 1605–1609. [CrossRef] [PubMed] 67. Haga, N.; Kobayashi, K.; Suzuki, T.; Maeo, K.; Kubo, M.; Ontani, M.; Mitsuda, N.; Demura, T.; Nakamura, K.; Jürgens, G.; et al. Mutations in MYB3R1 and MYB3R4 cause pleiotropic developmental defects and preferential down-regulation of multiple G2/M-specific genes in Arabidopsis. Plant Physiol. 2011, 157, 706–717. [CrossRef] 68. Frerigmann, H.; Gigolashvili, T. MYB34, MYB51, and MYB122 distinctly regulate indolic glucosinolate biosynthesis in Arabidopsis thaliana. Mol. Plant 2014, 7, 814–828. [CrossRef] 69. Finatto, T.; Viana, V.E.; Woyann, L.G.; Busanello, C.; da Maia, L.C.; de Oliveira, A.C. Can WRKY transcription factors help plants to overcome environmental challenges? Genet. Mol Biol. 2018, 41, 533–544. [CrossRef] 70. Caarls, L.; Pieterse, C.M.J.; Van Wees, S.C.M. How salicylic acid takes transcriptional control over jasmonic acid signaling. Front. Plant Sci. 2015, 6, 170. [CrossRef] 71. Li, J.; Brader, G.; Kariola, T.; Palva, E.T. WRKY70 modulates the selection of signaling pathways in plant defense. Plant J. 2006, 46, 477–491. [CrossRef] 72. Li, J.; Brader, G.; Palva, E.T. The WRKY70 transcription factor: A node of convergence for jasmonate-mediated and salicylate- mediated signals in plant defense. Plant Cell. 2004, 16, 319–331. [CrossRef] 73. Grunewald, W.; De Smet, I.; Lewis, D.R.; Löfke, C.; Jansen, L.; Goeminne, G.; Vanden-Bossche, R.; Karimi, M.; De Rybel, B.; Vanholme, B.; et al. Transcription factor WRKY23 assists auxin distribution patterns during Arabidopsis root development through local control on flavonol biosynthesis. Proc. Natl. Acad. Sci. USA 2012, 109, 1554–1559. [CrossRef] [PubMed] 74. Miyamoto, T.; Uemura, T.; Nemoto, K.; Daito, M.; Nozawa, A.; Sawasaki, T.; Arimura, G. Tyrosine kinase-dependent defense responses against herbivory in Arabidopsis. Front. Plant Sci. 2019, 10, 1–9. [CrossRef] [PubMed] 75. Zhang, Z.; Wang, J.; Zhang, R.; Huang, R. The ethylene response factor AtERF98 enhances tolerance to salt through the transcriptional activation of ascorbic acid synthesis in Arabidopsis. Plant J. 2012, 71, 273–287. [CrossRef] 76. Bahieldin, A.; Atef, A.; Edris, S.; Gadalla, N.O.; Ramadan, A.M.; Hassan, S.M.; Al Attas, S.G.; Al-Kordy, M.; Al-Hajar, A.S.M.; Sabir, J.; et al. Multifunctional activities of ERF109 as affected by salt stress in Arabidopsis. Sci. Rep. 2018, 8, 1–10. [CrossRef] [PubMed] 77. Raines, T.; Shanks, C.; Cheng, C.Y.; McPherson, D.; Argueso, C.T.; Kim, H.J.; Franco-Zorrilla, J.M.; López-Vidriero, I.; Solano, R.; Vanková, R.; et al. The cytokinin response factors modulate root and shoot growth and promote leaf senescence in Arabidopsis. Plant J. 2016, 85, 134–147. [CrossRef] Int. J. Mol. Sci. 2021, 22, 1226 23 of 24 78. Liu, W.; Li, Q.; Wang, Y.; Wu, T.; Yang, Y.; Zhang, X.; Han, Z.; Xu, X. Ethylene response factor AtERF72 negatively regulates Arabidopsis thaliana response to iron deficiency. Biochem. Biophys. Res. Commun. 2017, 491, 862–868. [CrossRef] 79. Zhu, D.; Wu, Z.; Cao, G.; Li, J.; Wei, J.; Tsuge, T.; Gu, H.; Aoyama, T.; Qu, L.J. TRANSLUCENT GREEN, an ERF family transcription factor, controls water balance. Mol. Plant 2014, 7, 601–615. [CrossRef] 80. Krishnaswamy, S.; Verma, S.; Rahman, M.H.; Kav, N.N.V. Functional characterization of four APETALA2-family genes (RAP2.6, RAP2.6L, DREB19 and DREB26) in Arabidopsis. Plant Mol. Biol. 2011, 75, 107–127. [CrossRef] 81. Saelim, L.; Akiyoshi, N.; Tan, T.T.; Ihara, A.; Yamaguchi, M.; Hirano, K.; Matsuoka, M.; Demura, T.; Ohtani, M. Arabidopsis Group IIId ERF proteins positively regulate primary cell wall-type CESA genes. J. Plant Res. 2019, 132, 117–129. [CrossRef] 82. Wang, X.; Liu, S.; Tian, H.; Wang, S.; Chen, J.G. The small ethylene response factor ERF96 is involved in the regulation of the abscisic acid response in arabidopsis. Front. Plant Sci. 2015, 6, 1064. [CrossRef] 83. Mohanta, T.K.; Yadav, D.; Khan, A.; Hashem, A. Genomics, molecular and evolutionary perspective of NAC transcription factors. PLoS ONE 2020, 15, e0231425. [CrossRef] [PubMed] 84. Kim, H.J.; Park, J.H.; Kim, J.; Kim, J.J.; Sunghyun, H.; Kim, J.; Kim, J.H.; Woo, H.R.; Hyeron, C.; Lim, P.O.; et al. Time-evolving genetic networks reveal a nac troika that negatively regulates leaf senescence in arabidopsis. Proc. Natl. Acad. Sci. USA 2018, 115, 4930–4939. [CrossRef] [PubMed] 85. Naranjo-Arcos, M.A.; Maurer, F.; Meiser, J.; Pateyron, S.; Fink-Straube, C.; Bauer, P. Dissection of iron signaling and iron accumulation by overexpression of subgroup Ib bHLH039 protein. Sci. Rep. 2017, 7. [CrossRef] [PubMed] 86. Sivitz, A.B.; Hermand, V.; Curie, C.; Vert, G. Arabidopsis bHLH100 and bHLH101 Control Iron Homeostasis via a FIT-Independent Pathway. PLoS ONE 2012, 7. [CrossRef] [PubMed] 87. Pré, M.; Atallah, M.; Champion, A.; De Vos, M.; Pieterse, C.M.; Memelink, J. The AP2/ERF domain transcription factor ORA59 integrates jasmonic acid and ethylene signals in plant defense. Plant Physiol. 2008, 147, 1347–1357. [CrossRef] 88. Geu-Flores, F.; Moldrup, M.E.; Böttcher, C.; Olsen, C.E.; Scheel, D.; Halkier, B.A. Cytosolic γ-glutamyl peptidases process glutathione conjugates in the biosynthesis of glucosinolates and camalexin in arabidopsis. Plant Cell. 2011, 23, 2456–2469. [CrossRef] 89. Jeschke, V.; Weber, K.; Moore, S.S.; Burow, M. Coordination of Glucosinolate Biosynthesis and Turnover Under Different Nutrient Conditions. Front. Plant Sci. 2019, 10, 1–17. [CrossRef] 90. Howarth, J.R.; Parmar, S.; Barraclough, P.B.; Hawkesford, M.J. A sulphur deficiency-induced gene, sdi1, involved in the utilization of stored sulphate pools under sulphur-limiting conditions has potential as a diagnostic indicator of sulphur nutritional status. Plant Biotechnol. J. 2009, 7, 200–209. [CrossRef] 91. Garcia-Molina, A.; Altmann, M.; Alkofer, A.; Epple, P.M.; Dangl, J.L.; Falter-Braun, P. LSU network hubs integrate abiotic and biotic stress responses via interaction with the superoxide dismutase FSD2. J. Exp. Bot. 2017, 68, 1185–1197. [CrossRef] 92. Zhang, L.; Kawaguchi, R.; Morikawa-Ichinose, T.; Allahham, A.; Kim, S.J.; Maruyama-Nakashita, A. Sulfur deficiency-induced glucosinolate catabolism attributed to two β-glucosidases, bglu28 and bglu30, is required for plant growth maintenance under sulfur deficiency. Plant Cell Physiol. 2020, 61, 803–813. [CrossRef] 93. Clark, T.J.; Lu, Y. Analysis of loss-of-function mutants in aspartate kinase and homoserine dehydrogenase genes points to complexity in the regulation of aspartate-derived amino acid contents. Plant Physiol. 2015, 168, 1512–1526. [CrossRef] [PubMed] 94. De Schutter, K.; Joubès, J.; Cools, T.; Verkest, A.; Corellou, F.; Babiychuck, E.; Van der Schueren, E.; Beeckman, T.; Kushnir, S.; Inzé, D.; et al. Arabidopsis WEE1 kinase controls cell cycle arrest in response to activation of the DNA integrity checkpoint. Plant Cell. 2007, 19, 211–225. [CrossRef] [PubMed] 95. Nisa, M.; Huang, Y.; Benhamed, M.; Raynaud, C. The plant DNA damage response: Signaling pathways leading to growth inhibition and putative role in response to stress conditions. Front. Plant Sci. 2019, 10, 653. [CrossRef] [PubMed] 96. Yoshikawa, T.; Ito, M.; Sumikura, T.; Nakayama, A. The rice FISH BONE gene encodes a tryptophan aminotransferase, which affects pleiotropic auxin-related processes. Plant J. 2014, 78, 927–936. [CrossRef] [PubMed] 97. Miao, H.; Wei, J.; Zhao, Y.; Yan, H.; Sun, B.; Huang, J.; Wang, Q. Glucose signaling positively regulates aliphatic glucosinolate biosynthesis. J. Exp. Bot. 2013, 64, 1097–1109. [CrossRef] 98. Ramírez, L.; Bartoli, C.G.; Lamattina, L. Glutathione and ascorbic acid protect Arabidopsis plants against detrimental effects of iron deficiency. J. Exp. Bot. 2013, 64, 3169–3178. [CrossRef] [PubMed] 99. Li, L.; He, Z.; Pandey, G.K.; Tsuchiya, T.; Luan, S. Functional cloning and characterization of a plant efflux carrier for multidrug and heavy metal detoxification. J. Biol. Chem. 2002, 277, 5360–5368. [CrossRef] 100. Fukudome, A.; Aksoy, E.; Wu, X.; Kumar, K.; Jeong, I.S.; May, K.; Russell, W.K.; Koiwa, H. Arabidopsis CPL4 is an essential C-terminal domain phosphatase that suppresses xenobiotic stress responses. Plant J. 2014, 80, 27–39. [CrossRef] 101. Lucena, C.; Romera, F.J.; Rojas, C.L.; García, M.J.; Alcántara, E.; Pérez-Vicente, R. Bicarbonate blocks the expression of several genes involved in the physiological responses to Fe deficiency of Strategy I plants. Funct Plant Biol. 2007, 34, 1002–1009. [CrossRef] 102. Wang, N.; Dong, X.; Chen, Y.; Ma, B.; Yao, C.; Ma, F.; Liu, Z. Direct and bicarbonate-induced iron deficiency differently affect iron translocation in Kiwifruit roots. Plants 2020, 9, 1578. [CrossRef] 103. Marschner, H. Mineral Nutrition of Higher Plants, 2nd ed.; Academic Press: London, UK, 1995; ISBN 9780124735439. 104. Poschenrieder, C.; Fernández, J.A.; Rubio, L.; Pérez-Martín, L.; Terés, J.; Barceló, J. Transport and use of bicarbonate in plants: Current knowledge and challenges ahead. Int. J. Mol. Sci. 2018, 19, 1352. [CrossRef] [PubMed] Int. J. Mol. Sci. 2021, 22, 1226 24 of 24 105. Busoms, S.; Teres, J.; Huang, X.-Y.; Bomblies, K.; Danku, J.; Douglas, A.; Weigel, D.; Poschenrieder, C.; Salt, D.E. Salinity Is an Agent of Divergent Selection Driving Local Adaptation of Arabidopsis to Coastal Habitats. Plant Physiol. 2015, 168, 915–929. [CrossRef] [PubMed] 106. Scholl, R.L.; May, S.T.; Ware, D.H. Seed and Molecular Resources for Arabidopsis. Plant Physiol. 2000, 124, 1477–1480. [CrossRef] [PubMed] 107. Schneider, C.A.; Rasband, W.S.; Eliceiri, K.W. NIH Image to ImageJ: 25 years of image analysis. Nat. Methods 2012, 9, 671–675. [CrossRef] [PubMed] 108. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014, 15, 550. [CrossRef] 109. Benjamini, Y.; Yekutieli, D. The Control of the False Discovery Rate in Multiple Testing under Dependency. Ann. Stat. 2011, 29, 1165–1188. 110. Tian, T.; Liu, Y.; Yan, H.; You, Q.; Yi, X.; Du, Z.; Xu, W.; Su, Z. AgriGO v2.0: A GO analysis toolkit for the agricultural community, 2017 update. Nucleic Acids Res. 2017, 45, 122–129. [CrossRef] 111. Szklarczyk, D.; Gable, A.L.; Lyon, D.; Junge, A.; Wyder, S.; Huertas-Cepas, J.; Simonovic, M.; Doncheva, N.T.; Morris, J.H.; Bork, P.; et al. STRING v11: Protein–protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019, 47, 607–613. [CrossRef] 112. Ye, J.; Coulouris, G.; Zaretskaya, I.; Cutcutache, I.; Rozen, S.; Madden, T.L. Primer-BLAST: A tool to design target-specific primers for polymerase chain reaction. BMC Bioinform. 2012, 13, 134. [CrossRef] 113. Dekkers, B.J.W.; Willems, L.; Bassel, G.W.; van Bolderern-Veldkamp, R.P.M.; Ligterink, W.; Hilhorst, H.W.M.; Bentsink, L. Identification of Reference Genes for RT–qPCR Expression Analysis in Arabidopsis and Tomato Seeds. Plant Cell Physiol. 2012, 53, 28–37. [CrossRef] 114. Livak, K.J.; Schmittgen, T.D. Analysis of Relative Gene Expression Data Using Real-Time Quantitative PCR and the 2−∆∆CT Method. Methods 2001, 25, 402–408. [CrossRef] [PubMed]