Comparative characterization of aluminum responsive transcriptome in Arabidopsis roots: comparison with other rhizotoxic ions at different stress intensities

ABSTRACT

Under abiotic stress conditions, plants actively modify the transcriptome either for activation of stress tolerance mechanisms or adjustment of biological processes for adaptation to the stress. Comparative analysis of transcriptome data under different stressors and stress intensities is useful to understand how plants balance growth and survival under stress conditions. In this study, the transcriptome of Arabidopsis roots was investigated in response to four rhizotoxic stressors, namely aluminum (Al³⁺), cadmium (Cd²⁺), and copper (Cu²⁺) ions and sodium chloride (NaCl), under mild stress (50% inhibition of root elongation compared with the control) and severe stress (>90% inhibition of root elongation). Only Al treatment showed similar transcriptome responses between the stress intensities. The shared Al-specific up-regulated genes included known Al-tolerance genes and downstream genes of SENSITIVE TO PROTON1 (STOP1) such as ALUMINUM ACTIVATED MALATE TRANSPORTER1 and malic enzymes, suggesting that Al-defense mechanisms regulated by STOP1 might be activated over a broader range of Al stress intensities. The sulfur metabolism pathway was activated only under severe Al stress, suggesting that rapid activation of sulfur metabolism is crucial for survival of root cells in the early stage of severe Al stress. A large number of genes associated with the response mechanism to reactive oxygen species production were up-regulated by Cd and Cu treatment under both severe and mild stress, but the members of the genes sets differed between the stress severities. The results indicated that Cd and Cu stress trigger common transcriptional responses, but the responses are likely associated with the stress severity. The transcriptional response to mild NaCl stress was relatively smaller compared with those to the other stressors, suggesting that in Arabidopsis different mechanisms other than rapid transcriptomic responses are involved in the response to NaCl stress. Comparison between Al and the other stressors indicated that the STOP1 regulation system and central metabolic pathways contribute to the similarity of the transcriptional response to different intensities of Al stress. The results suggest that robustness of transcriptionally regulated Al tolerance systems to different stress intensities is important for both of growth and survival in plant roots under Al stress.

KEY WORDS:

• Aluminum stress
• Arabidopsis thaliana
• cadmium stress
• co-expression
• copper stress
• microarray
• rhizotoxic ions
• root growth
• salt stress

1. Introduction

As sedentary organisms, plants must adapt to environmental stresses, but are required to maintain a balance between growth and survival in agricultural lands and natural areas (Claeys and Inze 2013). The root system is important for maintaining growth because it is essential for acquisition of water and nutrients, and thus it sometimes determines growth of the whole plant. Rhizotoxicity of cadmium (Cd) and copper (Cu) ions may be encountered and enhanced by man-made activities, such as soil contamination resulting from mining activities, whereas aluminum (Al) and NaCl rhizotoxicities occur naturally (Wu et al. 1996; Kobayashi et al. 2016). In particular, Al rhizotoxicity, which disrupts root growth of diverse crops, is common in acid soils of pH < 5.5 and affects 40% of arable lands worldwide. Identification of molecular mechanisms for tolerance to rhizotoxicity in plants, especially for Al tolerance, would contribute to improvement of crop productivity globally in stressed soil conditions.

Although many Al tolerance genes have been identified by forward genetic approaches, such as fine mapping of mutants (see reviews by Liu et al. 2014; Daspute et al. 2017), additional studies using genome-wide approaches, including transcriptome analyses, are required to understand the complex molecular physiological mechanisms of Al tolerance (Yokosho et al. 2014; Li et al. 2017; Xu et al. 2017). Comparative transcriptomic analysis that includes several types of rhizotoxic stressors is one experimental approach to investigate complex stress tolerance mechanisms, such as that of Al tolerance, because the responsive genes and biological pathways can be classified into common and specific groups. In a previous study, we analyzed the root transcriptome response to relatively strong rhizotoxic stresses, namely Al, Cd, and Cu ions and NaCl, in Arabidopsis (Zhao et al. 2009). The transcriptome response to rhizotoxic treatments using concentrations that caused almost complete cessation of growth (i.e. > 90% growth inhibition) indicated that the transcriptional response to reactive oxygen species (ROS) and calcium (Ca) homeostasis were shared with all tested ions for survival rather than growth. By contrast, several known Al tolerance genes, such as ALUMINUM ACTIVATED MALATE TRANSPORTER1 (AtALMT1; Hoekenga et al. 2006), are relatively specifically inducible by Al rather than other rhizotoxic stressors (Zhao et al. 2009). Additional comparative analyses of responses to different stress conditions, such as milder stress conditions (e.g. 50% growth inhibition), may provide molecular information on Al tolerance mechanisms that operate with maintenance of root growth. In addition, comparative transcriptomic analyses between Al and other rhizotoxic stressors may help to reveal differential tolerance strategies of plants to other rhizotoxic ions and highlight the growth/survival strategies against Al toxicity by considering commonality and specificity among the stressors.

Previous transcriptome analyses of responses to Al stress usually compared Al treatment and the control condition and dose-response transcriptional information is limited at the transcriptome level. These analyses identified inducible expression of several major Al tolerance genes, which individually contribute to Al tolerance, such as AtALMT1 and ALUMINUM SENSITIVE 3 (ALS3) in Arabidopsis (Sawaki et al. 2009). However, the pH and Al concentration varies with soil depth and fluctuates with fertilizer application (You et al. 2015). It is well known that different defense mechanisms are employed for different intensities of drought stress. Severe drought reduces maximum Rubisco carboxylation and ribulose 1,5-bisphosphate regeneration results in decreased photosynthesis (Flexas and Medrano 2002; Chaves et al. 2009), whereas stomatal closure limits photosynthesis under milder drought stress (Flexas et al. 2014). The increased isoprene emission and zeaxanthin concentration observed under mild drought stress may protect the phytosynthetic apparatus against stresses (Tattini et al. 2015). These differences imply that plants have multiple physiological and transcriptional mechanisms that are employed to overcome a particular stress with adjustment of the balance between growth and survival. Clarification of common and different transcriptional mechanisms will contribute to the establishment of advanced breeding strategies. In addition, it remains unclear whether Al alters transcription in other biological processes in accordance with the balance between maintenance of growth and survival, such as metabolic pathways, hormone metabolism, signaling pathways that regulate development, and cell division. It also unknown whether the overall pattern of transcriptomic adaptation to Al stress is similar to that for other ions.

In this study, we adopted a comparative transcriptomic approach for (1) determination of similarities and differences in the transcriptional response to mild stress compared with that for severe stress, (2) characterization of growth/survival strategies against Al stress by comparison with multiple ion stresses, and (3) improvement in understanding how plants adapt to different rhizotoxic stress conditions through maintaining a balance between growth and survival. The root transcriptome of Arabidopsis thaliana plants grown under mild stress (50% inhibition of root elongation compared with root growth under the control condition) induced by the rhizotoxic stressors Al, Cd, and Cu ions and NaCl was compared with previously published root transcriptome data for the same rhizotoxic stressors under severe stress (>90% root growth inhibition). The results show that distinct pattern of transcriptomic adaptation by each rhizotoxic ion, including greater similarity of inducible genes group by Al at different stress intensities. This greater similarity of transcriptomic pattern of Al than other ions, in other words robustness of transcriptomic adaptation pattern at different stress intensities of Al, would be explained by the involvement of STOP1 regulated Al-inducible tolerance system, including upregulation of AtALMT1 with changes in gene expression for malate metabolism.

2. Methods

2.1. Rhizotoxic stressor treatments

Rhizotoxic stressor treatments and root sampling were performed as described by Zhao et al. (2009). Arabidopsis thaliana accession Col-4 (N933, NASC; Nottingham Arabidopsis Stock Center, Nottingham, UK) seedlings were hydroponically pre-grown in modified 1/50 MGRL medium (Fujiwara et al. 1992) with 200 μM of CaCl₂ and inorganic phosphorus at pH 5.0 for 10 days and transferred to the 1/50 MGRL solution with 200 μM of CaCl₂ except inorganic phosphorus containing either 6 μM AlCl₃, 3 μM CdCl₂, 1.3 μM CuSO₄, or 10 mM NaCl for 24 h. Concentrations of each ion and NaCl were adjusted to achieve 50% root growth inhibition compared with root growth under the control condition during a 1-week growth test. As a control, pre-grown seedlings were transferred to modified 1/50 MGRL medium lacking rhizotoxic stressors. Al free activity was calculated by GEOCHEM-EZ software version 1.0 (Shaff et al. 2010). Free Al ion at the plasma membrane (PM) surface was calculated by a SGCS electrostatic model (Kinraide and Wang 2010).

2.2. RNA isolation

Total RNA was isolated using the method described by Suzuki et al. (2004), then quantified at A₂₆₀ using a NanoDrop™ spectrophotometer (ND-1000, NanoDrop Technologies, Wilmington, DE, USA). The quality of RNA used for microarray analysis was measured using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) in accordance with the manufacturer’s instructions.

2.3. Microarray experiment

Microarray analyses were performed using the Agilent microarray system (Agilent Technologies, Palo Alto, CA, USA) as described by Zhao et al. (2009). All procedures were carried out in accordance with the manufacturer’s protocols. Briefly, 1 μg total RNA from each sample was used to synthesize cRNA and was labeled with cyanine-5 (Cy5)- or cyanine-3 (Cy3)-labeled CTP (Perkin Elmer/NEN Life Sciences, Tokyo, Japan). The labeled cRNAs (including a treatment and a control sample) were competitively hybridized to the Arabidopsis 3 Oligo Microarray and then washed. The hybridized slides were scanned using Agilent DNA Microarray Scanner version 6.1 software and data points were extracted using Agilent Feature Extraction version 8.1 software. Three comparisons, including one dye-flip, were made between biologically independent samples. All microarray data have been deposited in the Gene Expression Omnibus of the National Center for Biotechnology Information under the accession number GSE108751.

2.4. Data analysis

All statistical analyses of microarray data were performed using R version 3.4.3 available at CRAN (http://cran.r-project.org/). Previously studied microarray data for the severe rhizotoxic ion stress condition (E-MEXP-1907) and the stop1 mutant under the Al stress condition (E-MEXP-1908) were downloaded from the ARRAYEXPRESS database to compare transcriptional similarities and differences between different stress intensities with consideration of STOP1 regulation. Quality control of the microarray data was performed using the default settings of GeneSpring GX version 12.6.1 software (Agilent, Tokyo, Japan). Probes for the mild and severe stress conditions were connected by their AGI codes. A total of 21,399 genes that passed the quality control in all 24 microarray experiments were used for statistical analyses and gene selection. Up- and down-regulated genes in response to an individual stressor were defined if their fold change values were in the upper or lower 2.5% (i.e. upper or lower 534 genes) in all three replications. Genes up-regulated by all mild stressors and unique gene groups were identified by Venn diagrams using the R package VennDiagram (Chen and Boutros 2011). Functional categorization of differentially expressed genes was performed using PANTHER (Thomas et al. 2003). Pathway enrichment analysis for up-regulated genes was performed using the BioAnalytic Resource Classification SuperViewer (Provart et al. 2003) and MapMan version 3.5.1 software (Thimm et al. 2004). Co-expression gene networks of up-regulated genes were prepared using NetworkDrawer implemented in ATTED-II (Aoki et al. 2016) with the ‘add a few genes’ option. Network density of individual gene modules was calculated using the NetWorkAnalyzer function implemented in Cytoscape version 3.2.0 software (Shannon et al. 2003). Enrichment analysis of octamers in the AtALMT1 promoter was performed using the method described by Yamamoto et al. (2011).

3. Results

3.1. Global analysis of transcriptomic response to mild and severe rhizotoxic stress treatments

Rhizotoxic stress treatments were applied using a hydroponic culture system to study the transcriptomic response of four rhizotoxic stressors at concentrations that induce 50% root growth inhibition (see Kobayashi et al. 2007). Either of 6 μM AlCl₃, 3 μM CdCl₂, 1.3 μM CuCl₂, and 10 μM NaCl was applied in the treatments. The transcriptome was obtained after incubation of seedlings for 24 h in each rhizotoxic solution. Transcriptomic responses to the severe stress condition (i.e. medium containing the rhizotoxic stressor at the concentration that induced > 90% root growth inhibition; Zhao et al. 2009), which were obtained in the previous study, were reanalyzed with the current data to compare transcriptomic responses at different levels of rhizotoxic stressors. It was predicted that Al-free activity and PM surface Al ion activity were 1.4 μM and 13.2 μM under the mild Al stress treatment, respectively. Al toxicity at the root plasma membrane was occurred even at moderate acidic condition when relatively lower Ca concentration used in this study (Kobayashi et al. 2013).

A total of 21,399 genes that passed the quality filters of all 24 microarray chips (two stress intensities, four stressors, three replications) were used for analyses. Scatter plots showed that NaCl altered the transcription of fewer genes than the other tested stressors (Fig. 1A; fold change > 2 or < 0.5). The correlation matrix (Fig. 1B) showed that similar transcriptomic responses between the mild and severe stress conditions were observed for Al ions only (r = 0.72), in contrast to the lower Pearson correlation coefficients observed for the other stressors (r = 0.37–0.52). This result suggests that the transcriptomic response to Al stress is broadly occurred against different stress intensities and degrees of root growth inhibition. Similarities in transcriptomic responses between Cd and Cu ions were observed both under mild and severe stress. Pearson correlation coefficients between Cd and Cu were 0.73 and 0.77 under mild and severe stress, respectively, whereas the correlation coefficients between mild and severe stress were 0.52 and 0.42 for Cd and Cu, respectively. This result suggests that the heavy metals Cd and Cu induce a similar transcriptional response to that of Al and NaCl. In addition, different transcriptional systems might operate for the different stress intensities under Cd and Cu stress.

Figure 1. Whole transcriptome response to rhizotoxic stressor treatment. (A) Scatter plot of competitive microarray data from roots of Arabidopsis subjected to mild stress of rhizotoxic stressors. Roots of hydroponically grown seedlings were treated by control (pH 5.0, no toxicant) and rhizotoxic solutions containing AlCl₃ (6 μM), CdCl₂ (3 μM), CuCl₂ (1.3 μM), or NaCl (10 mM) at pH 5.0 for 24 h. X and Y axes indicate signal intensities in the control and rhizotoxic treatments, respectively. The mean signal intensities from three biologically independent replications are plotted. The slope of the lines in each panel represent 2, 1, and 1/2 fold changes, respectively. (B) Pearson correlation matrix for the whole transcriptome in response to four rhizotoxic stressors. Correlation coefficients of 21,399 genes are shown for different stressors (Al, Cd, and Cu ions and NaCl) compared with the mild (upper matrix) or severe (lower matrix) stress conditions.

4. Stressor-specific and common responsive genes

To accentuate the responsive genes for each stressor, the most up-regulated or down-regulated genes were grouped by ranking the fold change of the genes on the basis of whether the fold change value was in the upper and lower 2.5% in each microarray experiment (Zhao et al. 2009). The current gene selection methodology using ranking-based method has an advantage to compare transcriptome profiles with different dynamic range of fold change (i.e. NaCl induced much less genes than the other ions; Fig. 1 and Supplementary Figure S1) and to extract the outstanding biological processes activated by individual stresses Additionally, our ranking-based criterion placed importance on reproducibility of gene lists in the three replicates (Guo et al. 2006; MAQC Consortium 2006) rather than collecting all of possible induced/suppressed genes. Of the 21,399 genes, 403, 382, 380, and 110 genes were classified as up-regulated by Al, Cd, Cu, or NaCl, respectively (Fig. 2). On the other hand, 370, 239, 248, and 28 genes were classified as down-regulated by Al, Cd, Cu, or NaCl, respectively. Those genes were overlapped for the most part with up- or down-regulated genes selected by a conventional gene selection method using statistical test based on t-statistic and fold change (Supplementary Figure S1). Salt stress induced and suppressed fewer genes under mild stress intensity than those under severe stress intensity, although a similar number of genes were up- or down-regulated in response to the other stressors. These results indicated that the transcriptional response to NaCl was relatively mild or slow compared with that to the other stressors, despite inducing the same degree of root growth inhibition.

Figure 2. Number of differentially expressed (up- and down-regulated) genes in Arabidopsis roots after treatment with rhizotoxic stressors for 24 h. Up- and down-regulated genes under the mild and severe stress conditions for Al, Cd, and Cu ions and NaCl. The genes were defined as those with fold change (treatment/control) values in the upper and lower 2.5% in all three replicates.

Venn diagrams identified gene groups most specifically responsive to each stressor, and the overlap between the most up- and down-regulated genes of each stressor at the different stress intensities (Fig. 3A). These gene groups reflect the specificity of the transcriptional response to the four tested rhizotoxic stressors. The highest ratio was observed in the overlapped genes between Cd and Cu treatment (Fig. 3B), which corresponded to the high correlation of the whole transcriptome response between Cd and Cu. Thirty-one genes were up-regulated in response to all four stressors under mild stress, whereas 58 genes were up-regulated in common in response to severe stress. No genes were suppressed by all four rhizotoxic stressors, although 13 genes were down-regulated in common under the severe stress intensity. Although fewer genes were up-regulated or down-regulated by all four stressors under the mild stress intensity because significantly fewer genes were selected in the mild NaCl stress treatment, many genes associated with ROS scavenging, such as GLUTATHIONE S-TRANSFERASE F2, GLUTATHIONE S-TRANSFERASE F3, PEROXIDASE 34, PEROXIDASE 52, PEROXIDASE 58, and a Ca-homeostasis associated gene that encodes a putative calcium-binding protein CML47, were in the common up-regulated gene group (Table 1). Although most of the up-regulated genes in common under mild stress were up-regulated in at least one severe stress condition, only three genes, namely NODULIN MTN21-LIKE TRANSPORTER UMAMIT36, ankyrin (At3g13950), and LACCASE 12, responded only to milder stress conditions and thus might be associated with stress tolerance without growth retardation.

Table 1. Upregulated genes by all four stressors (Al, Cd, Cu, and NaCl) under milder stress intensity.

		Average log2FC value under mild stress				Average log2FC value under severe stress
Gene ID	Description	Al	Cd	Cu	Na	Al	Cd	Cu	Na
AT1G01680	U-box domain-containing protein 54	4.4	3.8	2.6	1.9	4.0	3.7	3.8	3.0
AT1G08430	Aluminum-activated malate transporter 1	6.5	3.8	3.4	1.3	6.3	−0.3	0.3	1.1
AT1G26420	FAD-binding and BBE domain-containing protein	2.5	3.2	2.4	2.0	2.3	4.5	3.6	3.0
AT1G47890	Receptor-like protein 7	2.8	3.1	2.6	1.1	1.8	4.4	2.1	0.1
AT1G55230	Hypothetical protein	6.0	5.3	5.5	2.3	3.8	3.5	1.2	0.9
AT1G57630	Toll-Interleukin-Resistance domain-containing protein	3.7	4.0	4.2	1.8	2.4	4.5	5.1	2.9
AT1G64160	Dirigent-like protein DIR5	4.3	5.4	6.6	3.3	2.6	6.0	5.1	4.1
AT1G70260	Nodulin MtN21-like transporter UMAMIT36	2.5	3.9	3.5	2.2	0.1	−0.3	−0.3	0.3
AT1G79680	Wall-associated receptor kinase-like 10	3.7	2.9	5.9	2.1	0.4	3.0	3.4	0.8
AT2G02930	Glutathione S-transferase F3	2.4	3.6	4.6	1.3	2.6	4.4	4.4	3.1
AT2G26560	Phospholipase A 2A	3.5	4.4	5.7	1.6	2.1	6.0	7.1	1.7
AT2G43510	Trypsin inhibitor protein 1	5.8	4.0	4.3	2.1	3.7	3.1	5.4	1.8
AT3G12230	Serine carboxypeptidase-like 14	4.1	5.4	5.4	1.9	3.1	1.9	1.6	−0.2
AT3G13950	Hypothetical protein	3.8	3.6	2.7	1.3	0.3	0.5	0.3	0.1
AT3G22910	Putative calcium-transporting ATPase 13	4.8	3.9	3.2	2.3	2.9	3.7	3.1	3.2
AT3G23250	MYB domain protein 15	2.6	3.2	3.2	1.2	1.7	4.1	4.0	2.7
AT3G28510	AAA-type ATPase family protein	4.9	5.7	5.1	3.6	4.3	4.2	0.9	1.5
AT3G44870	S-adenosyl-L-methionine-dependent methyltransferase	3.5	2.8	4.0	1.2	1.9	4.2	4.9	3.1
AT3G47340	Asparagine synthetase [glutamine-hydrolyzing]	4.9	4.0	4.2	3.1	3.8	4.0	5.0	1.6
AT3G47480	Putative calcium-binding protein CML47	4.3	3.8	2.8	1.1	3.8	4.7	4.0	2.5
AT3G48850	Phosphate transporter 3;2	2.8	3.0	2.6	1.7	0.7	3.0	2.3	1.1
AT3G49120	Peroxidase 34	3.2	2.2	2.1	1.2	2.2	2.5	3.1	1.4
AT4G02520	Glutathione S-transferase F2	2.4	2.8	3.3	1.1	2.3	3.3	3.4	2.6
AT4G35770	Senescence-associated protein DIN1	3.9	3.0	4.1	3.1	1.6	1.1	4.7	0.8
AT4G36880	Cysteine proteinase1	4.8	2.8	4.4	1.8	2.4	0.5	1.2	1.0
AT4G37290	Hypothetical protein	3.9	2.8	4.3	1.8	1.0	3.6	4.4	0.8
AT5G05340	Peroxidase 52	4.9	4.1	4.9	3.6	4.3	5.2	5.5	3.0
AT5G05390	Laccase 12	2.4	3.2	2.8	1.2	0.2	0.0	0.9	0.4
AT5G19880	Peroxidase 58	2.9	4.2	3.9	3.3	3.1	4.9	5.7	3.1
AT5G22530	Hypothetical protein	2.6	4.2	4.3	2.1	1.9	4.1	1.9	2.6
AT5G22570	Putative WRKY transcription factor 38	2.9	3.3	3.3	1.6	1.8	2.3	1.0	1.5

Averaged values of fold change in three microarray replicates were shown.

Figure 3. Venn diagram of differentially expressed genes in response to individual rhizotoxic stressors. (A) The top 2.5% genes whose expression was most increased or decreased by each stressor drawn using the VennDiagram R package. The number of genes that was specific to each stressor is indicated in the small circle for comparison between the severe and mild stress intensities. (B) The proportion of genes overlapped between each two upregulated gene sets. Different letters indicate a significant difference in the ratio (Fisher’s exact test; p < 0.05). (C) Evaluation of specificity in individual stressor-dependent gene expression. Specifically expressed genes were defined as individual stressor-specific expressed genes in the Venn diagram without up-regulation by the other stressors (FC <2). Asterisks indicate significant difference in the ratio of specifically expressed genes (Fisher’s exact test; p < 0.05).

Venn diagrams also indicated that 182, 104, 114, and 33 genes showed stress-specific up-regulation under mild Al, Cd, Cu, and NaCl stress, respectively (Fig. 3A; Supplementary Tables S1–S4). A higher proportion of genes were specifically responsive under mild stress than those under severe stress (Fig. 3C). Typical Al tolerance genes, such as MULTIDRUG AND TOXIC COMPOUND EXTRUSION (MATE; Magalhaes et al. 2007) and ALUMINIUM-SENSITIVE 3 (ALS3; Larsen et al. 2005), were among the up-regulated genes under mild Al stress (Supplementary Table S1), whereas AtALMT1 was identified in the severe stress intensity. Mild Al stress specifically induced five genes, which were previously identified as expression level polymorphic genes in a RNA-Seq analysis employing six Arabidopsis thaliana accessions with contrasting degrees of Al tolerance: MYB112, LACCASE 7, NAC74, CHITINASE C (ChiC), and PLASMODESMATA-LOCATED PROTEIN 1 (PDLP1) (Kusunoki et al. 2017). These results suggested that several Al tolerance genes are highly specifically inducible by Al compared with the other tested stressors. Up-regulated genes under mild Cd stress included NAD(P)-binding Rossmann-fold superfamily protein (At1g07440) and MITOCHONDRION-LOCALIZED SMALL HEAT SHOCK PROTEIN 23.6 (ATHSP23.6-MITO), which were detected by proteome analysis in response to Cd stress in A. thaliana (Brahim et al. 2010; Sarry et al. 2006). Genes specifically responsive to NaCl included C2-DOMAIN ABA-RELATED 5 (CAR5), which is a positive regulator of abscisic acid-activated signaling pathways (Rodriguez et al. 2014). Enriched GO-BP terms were identified for the genes up-regulated in response to mild Cu stress, including ethylene-activated signaling pathways consisting of six ethylene response factors (ERF11, ERF016, ERF094, ERF096, ERF114, and ERF115) (Supplementary Table S5).

5. Possible mechanisms for the induction of the genes groups in response to both mild and severe Al stress

To investigate the possible mechanism for the up-regulated genes in response to different intensities of Al stress, we examined Al-inducible genes that overlapped between the mild and severe stress intensities and performed a co-expression network analysis. When the specific genes were compared between the stress intensities, a relatively larger overlap was identified for Al (29 genes) than for the other stressors (Fig. 3A; Table 2). Genes regulated by the STOP1 transcription factor (i.e. fold change (stop1/wild type) < 0.5) were significantly enriched in the 29 genes shared in both stress conditions (p = 5.257E-05 in Fisher’s Exact Test), including six genes such as an NADP-dependent malic enzyme (ME2) and a sulfate transporter (SURTR3;1). In addition, we performed a co-expression network analysis of Al-specific up-regulated genes under mild and severe stress intensities. Eleven small co-expression genes networks were identified by the ATTED-II co-expression database for the severe Al stress condition. Two (of which one was large) co-expression networks were identified for mild Al stress. Three of the co-expression networks identified for the severe Al stress condition were completely or partially included in the large network identified for mild Al stress. Two of the shared networks for severe Al stress contained four STOP1-regulated genes, whereas 11 STOP1-regulated genes were included in the large network derived for mild Al stress. These results suggested that a larger number of STOP1-regulated genes is one factor in the transcriptomic similarity to different stress intensities. Pearson’s correlation coefficient value of whole transcriptome between severe and mild Al stress was slightly decreased when it was calculated excluding STOP1-regulated genes, while it was not observed in the other stressors (Supplementary Figure S2). This hypothesis is further supported by higher relative appearance of part of STOP1-binidng sequence in the AtALMT1 promoter under both Al stress intensities (Supplementary Figure S3; Tokizawa et al. 2015). Pathway enrichment analysis of the Al-specific up-regulated genes identified three enriched MapMan-pathways (‘TCA/org transformation’, ‘development’, and ‘transport’) shared by the two stress conditions (Supplementary Table S6). When the up-regulated genes were mapped on the MapMan metabolic map, malate-related pathways, such as malic acid and malate dehydrogenase, were enriched (Supplementary Figure S4), which suggested that malate metabolic processes operated to support malate release by ALMT1.

Table 2. Al-specific upregulatd genes under both severe and mild stress intensities.

		Average log2FC value under mild stress				Average log2FC value under severe stress
Gene ID	Description	Al	Cd	Cu	Na	Al	Cd	Cu	Na	Regulated by STOP1
AT1G11190	Endonuclease 1	2.8	0.3	0.6	0.0	2.0	0.7	1.7	0.3
AT1G35140	Protein EXORDIUM LIKE 1	2.4	1.5	0.0	0.1	2.1	2.0	1.9	1.4
AT1G66180	Aspartyl protease family protein	2.6	1.6	0.8	−0.3	3.2	1.4	1.2	0.5
AT1G70170	Matrix metalloproteinase	2.6	1.4	1.3	0.4	2.0	1.8	2.6	1.2
AT1G72890	TIR-NBS class of disease resistance protein	3.8	1.0	1.0	0.2	2.6	1.9	2.1	0.1
AT1G80110	Phloem protein 2-B11	2.4	0.2	0.1	−0.1	2.4	1.7	1.1	1.3
AT2G01300	Hypothetical protein	3.1	0.9	0.6	−0.3	1.6	1.5	3.0	1.5
AT2G17500	Auxin efflux carrier family protein	2.2	0.0	−0.1	0.0	2.3	0.6	1.3	1.0
AT2G18480	Putative polyol transporter 3	3.5	1.1	0.1	0.3	2.8	−0.2	−0.1	0.9	Regulated
AT2G19900	NADP-dependent malic enzyme 1	2.5	1.4	1.8	−0.1	2.0	1.5	2.4	0.3
AT2G26690	Nitrate transporter 1.4	2.7	−0.4	0.2	−0.1	1.9	0.2	1.1	−0.1
AT2G38790	Hypothetical protein	3.9	1.0	1.4	0.2	1.8	1.5	1.9	1.5
AT2G40880	Cysteine proteinase inhibitor 3	3.1	0.6	−0.3	−0.2	2.2	1.5	1.2	1.2
AT3G06210	Hypothetical protein	2.9	0.2	−0.3	−0.1	1.8	0.5	0.3	−0.1
AT3G22240	Hypothetical protein	2.8	1.5	1.1	0.0	2.2	0.6	2.5	0.1	Regulated
AT3G48340	KDEL-tailed cysteine endopeptidase CEP2	3.7	0.7	0.8	−0.7	2.7	−1.0	−0.3	−0.4
AT3G51895	Sulfate transporter 3;1	3.8	1.3	1.9	−0.4	2.9	1.0	1.4	0.3	Regulated
AT3G52820	Purple acid phosphatase 22	3.6	0.2	0.9	0.0	1.9	0.3	1.0	0.3
AT3G60130	Beta glucosidase 16	2.2	0.1	0.7	0.0	1.9	1.8	2.4	0.4
AT4G22070	WRKY DNA-binding protein 31	2.6	1.3	−0.1	0.8	2.0	1.9	1.9	1.5
AT4G24690	Autophagy substrate NBR1	2.8	0.5	0.3	−0.1	2.1	1.8	2.0	0.3
AT4G34330	Hypothetical protein	3.0	0.1	0.8	0.1	1.8	0.2	0.5	0.2
AT4G34710	Arginine decarboxylase 2	3.2	1.4	0.5	0.1	3.4	1.3	2.7	1.0	Regulated
AT5G01640	Prenylated RAB acceptor 1.B5	2.6	1.3	0.4	0.3	1.8	2.1	1.3	0.7
AT5G11670	NADP-dependent malic enzyme 2	2.3	1.5	1.9	−0.1	2.1	1.8	2.6	0.3	Regulated
AT5G39520	Hypothetical protein	7.6	1.4	1.4	0.2	2.4	1.1	2.6	−0.2	Regulated
AT5G48180	Nitrile specifier protein 5	2.8	0.0	−0.2	0.2	1.9	1.2	2.5	0.6
AT5G58070	Temperature-induced lipocalin TIL1	2.2	0.0	−0.2	−0.2	1.5	0.2	0.9	0.6
AT5G59820	High light responsive zinc finger protein ZAT12	3.5	1.6	1.7	0.4	1.6	2.4	2.4	1.7

Averaged values of fold change in three microarray replicates were shown.

Enrichment of the GO-BP terms ‘glutathione metabolic process’ and ‘response to salicylic acid’ was observed in both Al stress treatments (Fig. 5; Supplementary Table S7), which suggested that a similar transcriptional response to ROS is an additional factor in the transcriptomic similarity to different stress intensities. Five and eight ROS-responsive genes up-regulated by REDOX RESPONSIVE TRANSCRIPTION FACTOR1 (RRTF1) overexpression (Matsuo et al. 2015) also were included in the co-expression gene networks constructed for the Al-specific up-regulated genes under severe and mild stress intensity, respectively (Fig. 4). Several different GO-BP terms were observed between the mild and severe Al stress conditions; ‘sulfate transmembrane transport’, ‘suberin biosynthetic process’, and ‘lateral root formation’ were only observed for severe Al stress, whereas ‘response to aluminum ion’ and ‘calcium ion transmembrane transport’ were only observed for mild Al stress. The KaPPA metabolic process viewer showed up-regulation of genes associated with the sulfur metabolism pathway for severe Al stress (Supplementary Figure S5). These results suggested that specific metabolic pathways were responsive to severe Al stress, although most of pathways were similar between the mild and severe Al stress conditions.

Figure 4. Co-expression gene network of Al-specific up-regulated genes under severe (A) and mild (B) stress. A co-expression gene network was constructed, using Al-specific up-regulated genes under severe and mild stress intensity as query genes, with the ATTED-II database. Closed and opened circles indicate query and co-expressed genes added according to co-expression relationships. Triangles and rectangles indicate genes regulated by SENSITIVE TO PROTON RHIZOTOXICITY1 (STOP1) and REDOX RESPONSIVE TRANSCRIPTION FACTOR1 (RRTF1), respectively. Dashed circles and ellipses indicate overlapped co-expression gene networks between the severe and mild stress intensities.

Figure 5. Common and specific biological pathways responsive to severe and mild Al stress intensity. Representative genes are listed in parentheses for each gene ontology (GO) term. Representative enriched GO terms in the biological process category for genes up-regulated under severe or mild Al stress. GO analysis was performed using the PANTHER classification system. All enriched GO-BP terms are listed in Supplementary Table S7.

6. Comparative analysis of genes up-regulated in response to both Cd and Cu

Under severe stress intensity, 164 genes were up-regulated in response to both Cd and Cu, whereas 218 up-regulated genes overlapped under mild Cd or Cu stress (Fig. 6). Responsive genes that overlapped between Cd and Cu stresses were compared between the severe and mild stress conditions to clarify the similarities and differences in the transcriptional response to severe and mild stress intensity (Fig. 3B). Many similarities were observed in the transcriptional responses to Cd and Cu stress, that is transcription factors (C2H2, WRKY), peroxidase, ethylene responsive transcription factors, auxin responsive proteins, laccases, germin-like proteins, gluthathione S-transferases, and ABC transporters were up-regulated in response to both Cd and Cu under severe or mild stress intensity. Many GO groups overlapped between Cd and Cu for each stress intensity, and contained ROS response-related GO terms such as ‘toxin catabolic process’, ‘hydrogen peroxide catabolic process’, and ‘response to reactive oxygen species’ including PEROXIDASE 15/49/54/58, GSTF2/3, and GSTU4. It is suggested that both Cd and Cu induce similar ROS-mediated signaling pathways. Up-regulated gene groups contained several hydrogen peroxide (H₂O₂) inducible genes, which were up-regulated in H₂O₂-accumulating transgenic Arabidopsis that overexpresses RRTF1 (Fig. 7).

Figure 6. Comparison of enriched GO-BP terms between Cd and Cu stresses under severe and mild stress intensity. Venn diagram showing overlap in genes up-regulated by Cd and Cu stress. The number in parentheses indicates genes regulated by REDOX RESPONSIVE TRANSCRIPTION FACTOR1 (RRTF1; Matsuo et al. 2015). Representative enriched GO terms in the biological process category were classified into three groups. GO analysis was performed using the PANTHER classification system.

Figure 7. Proportion of genes regulated by RRTF1 among the up-regulated genes by individual stressors and transcriptome wide. Proportions are out of 222 genes upregulated in transgenic Arabidopsis overexpressing REDOX RESPONSIVE TRANSCRIPTION FACTOR1 (RRTF1; Matsuo et al. 2015). Different letters indicate a significant difference in the ratio (Fisher’s exact test; p < 0.05).

By contrast, only a few GO terms were enriched only under severe or mild stress conditions. For example, genes for ‘salicylic acid mediated signaling pathway’ were enriched in response to severe stress treatment, whereas genes for ‘jasmonic acid mediated signaling pathway’ were enriched under the mild stress condition. Only one GO term (systemin acquired resistance in Cd-specific genes) was shared by Cd- or Cu-specific gene groups under the mild and severe stress conditions (Fig. 6). The GO terms ‘Indole glucosinolate biosynthetic process’ and ‘callose deposition’ were enriched only in the severe stress-specific gene group, including NSL1, MYB51, and CYP81F2. In the mild stress-specific gene group, ‘oxidation-reduction process’ included LACCASE 12, L-ASPARTATE OXIDASE, PEROXIDASE 28/34/39, CYP71B3, and CYP81F1 (Fig. 8; Supplementary Tables S8–S9).

Figure 8. Common and specific biological pathways between severe and mild stress intensities in the overlapped genes up-regulated by Cd and Cu stress. (A) Venn diagram of genes up-regulated in response to Cd and Cu stress. (B) Representative enriched GO terms in the biological process category were classified into three groups. GO analysis was performed using the PANTHER classification system. All enriched GO-BP terms are listed in Supplementary Table S9.

7. Discussion

Transcriptional regulation plays a major role in stress tolerance of plants during adaptation to stress conditions. In a previous comparative transcriptomic analysis of rhizotoxic ions in Arabidopsis, we identified both ion-specific responses and common transcriptional responses, including activation of several stress tolerance genes. In the previous study, we applied severe stress of Al, Cd, Cu, and NaCl using concentrations of the rhizotoxic stressors that almost completely inhibited root growth. In the present study, we obtained the transcriptome under milder stress conditions that maintained 50% root growth relative to the control condition. The present analysis further identified different transcriptional responses to the same four rhizotoxic stressors under the mild stress condition, and examined differences in the transcriptomic responses between the two stress intensities.

Our previous study demonstrated that ROS-scavenging enzymes and Ca-binding proteins play an important role in the response to exposure to Al, Cd, Cu, and NaCl under severe stress (Zhao et al. 2009). It is notable that the ROS response to the four rhizotoxic stressors under severe stress was also observed under the mild stress intensity, although fewer ROS-responsive genes were up-regulated (Table 1). This finding suggested that protection against ROS and maintenance of Ca-homeostasis is a basic defense response to any rhizotoxic stressor, even under the mild stress intensity. On the other hand, stressor-specific gene responses were expanded under the mild stress intensity for all of the tested rhizotoxic stressors (Fig. 3B). These results suggested that the stressor-specific gene response is dominant for stress tolerance while maintaining growth under the mild stress condition. Thus, utilization of mild stress intensities might be suitable for candidate gene selection in transcriptome experiments.

Comparative analysis between the severe and mild stress conditions showed that the transcriptomic response patterns differed between stress intensities. The four tested rhizotoxic stressors were classified into three groups on the basis of the transcriptomic similarity/difference between the severe and mild stress intensities (Fig. 9). The transcriptional response to Al stress was relatively similar between the severe and mild stress intensities in both the correlation analysis at the transcriptome level and in the overlap of Al-specific up-regulated genes compared with the other stressors (Figs 2 and 3). These results suggested that a similar transcriptional tolerance mechanism is triggered by different Al concentrations. This similarity can be partly explained by the Al-specific activation of genes regulated by the STOP1 transcription factor (Iuchi et al. 2007; Fig. 4; Table 2; Supplementary Figure S2–3). The mild Al stress intensity induced transcription of not only ALMT1 but also MATE and ALS3, which are representative Al-tolerance genes. Given that STOP1 is likely to be a ‘master switch’ for plant stress responses, recent studies indicate that it may regulate different gene sets, such as activation of genes responsive to phosphate deficiency (Balzergue et al. 2017; Mora-Macías et al. 2017). It suggests that several STOP1-regulated genes may require Al-specific activation, including other mechanisms such as up-regulation of the co-activator of target genes, or an Al-specific mechanism of post-translational activation of STOP1. Additionally, part of the Al-activated co-expression gene networks did not contain STOP1-regulated genes (Fig. 4). These gene groups may be regulated by other critical transcription factors, which show greater specificity to Al than to other ions, such as WRKY46 in Arabidopsis (Ding et al. 2013) and ABSCISIC ACID, STRESS AND RIPENING (ASR5; Arenhart et al. 2013) in rice. Up-regulation of ALMT1 and enrichment of tricarboxylic acid cycle-related genes, including malic enzymes, were observed under both severe and mild Al stress (Fig. 5; Table 2; Supplementary Figure S3). Rapid and robust activation of malate release mechanism against different Al stress intensities also likely contributed to the similar transcriptome between the severe and mild Al stress conditions, because Al detoxification by malate chelation is a foremost defense system against Al ions. In addition, up-regulation of transcripts for malate metabolism under different Al stress intensities suggested activation of operation of biochemical pH-stat pathway under both Al stress conditions can also contribute to the transcriptomic similarity (Sakano 1998). By contrast, GO terms associated with ROS response, such as ‘glutathione metabolic process’ and ‘response to salicylic acid’, overlapped between the severe and mild Al stress conditions, which showed that ROS response is another important response regardless of stress intensity. Conversely, genes associated with sulfate transport and metabolism were activated only under severe Al stress (Fig. 5; Supplementary Figure S4). Sulfur fertilizer alleviates Al toxicity in barley, citrus, wheat, and oilseed rape via improvement in antioxidant protection and decrease in ROS concentrations (Chen et al. 2013; Guo et al. 2017; Qian et al. 2014). Activation of sulfur transport and metabolism is likely involved in scavenging the high quantities of ROS produced under severe Al stress. The commonality of shared tolerance mechanisms (i.e. organic acid excretion, pH regulation, and sulfur-mediated metabolisms) at different stress intensities would explain the similar transcriptomic adaptation pattern between the different Al stress intensities. Greater transcriptomic similarity between different stress intensities of Al could be also explained by the limited molecular target of Al than that of Cd/Cu, then the simplicity of phytotoxicity of Al compared to that of the other stressors. For example, Cd and Cu competitively bind to the metal binding domains in the metal binding proteins regulating metal homeostasis such as heavy metal transporting ATPases and metal chaperons (Dudev and Lim 2014).

Figure 9. Schematic illustration of the proposed model for the difference in transcriptional response between mild and severe stress intensities.

In contrast to the other tested rhizotoxic stressors, under mild NaCl stress transcriptomic change was much smaller than that observed under severe NaCl stress and in response to the other stressors. Significantly fewer genes were up-regulated under mild NaCl stress (Fig. 1), which suggested that NaCl tolerance is regulated by non-transcriptional responses under mild NaCl stress. Several NaCl tolerance genes, such as SALT OVERLAY SENSITIVE 1, which encodes a protein kinase (CIPK24), were not induced in response to NaCl treatment. By contrast, the short-term (24 h) exposure to severe NaCl stress rapidly induced the drought/salt responsive pathway regulated by DEHYDRATION RESPONSIVE ELEMENT BINDING PROTEIN 1A (DREB1A), DREB subfamily A-5 (At1g22810 and At5g21960), and ZAT10, and drastically altered the root transcriptome for defense (Zhao et al. 2009). Overexpression of DREB1A results in severe growth retardation under normal growing conditions (Kasuga et al. 1999). These results suggested that the non-transcriptional defense system might be employed for root growth under moderate and weak salinity stress conditions.

The current analysis also detected a high correlation between the responses to Cd and Cu ions, both in the entire transcriptome and among up-regulated genes (Figs 2 and 3B; Supplementary Tables S8–S9). Both Cd and Cu are known to induce oxidative stress at the cellular level (Smeets et al. 2009; Cuypers et al. 2011). Overlapped up-regulated genes between the Cd and Cu stress treatments were mainly associated with ROS response/scavenging (Figs 6–8). Interestingly, glucosinolate biosynthesis and callose deposition genes were notable only under severe stress, which suggested that activation of these biological pathways is involved with survival under cessation of root growth. Mild stress also induced ROS response but in different gene members from those induced under severe stress.

In summary, we characterized differences in the transcriptional response to different stress intensities of four representative rhizotoxic stressors in relation to maintenance of the balance between growth and survival. The results suggested that the STOP1 regulation system and central metabolic pathways contributed to the similarity of the transcriptional response to different intensities of Al stress, and that robustness of Al tolerance systems to different stress intensities is more important for both of growth and survival in plant roots under Al stress than under the other tested rhizotoxic stressors. Commonality and specificity of the transcriptional response among the individual stressors clarified the biological processes involved in the tolerance to common and specific stressors. Stressor-specific up-regulated genes identified in this study can be investigated in future studies as candidate genes and physiological processes in relation to survival under severe stress and maintenance of plant growth under moderate stress, which are frequently encountered in agriculture. Combination of a genome-wide approach with the current approach may be required for efficient isolation of stress tolerance genes. For example, we recently identified several Al tolerance genes by studying expression level polymorphism between accession groups with contrasting Al tolerance (Kusunoki et al. 2017). In addition, integration of co-expression gene network analysis and genome-wide association studies for NaCl tolerance revealed previously unidentified NaCl tolerance mechanisms (e.g., involvement of the cell wall in NaCl tolerance; Kobayashi et al. 2016). A similar integrated approach may accelerate identification of genes involved in regulation of root growth under stress conditions.

Acknowledgments

We thank Robert McKenzie, PhD, from Edanz Group (www.edanzediting.com/ac), for editing a draft of this manuscript.

Disclosure statement

The authors have no conflicts of interest to declare.

Supplemental data

Supplemental data for this article can be accessed here.

Additional information

Funding

Part of this research was supported by JSPS KAKENHI grant number 22248006.