Comparative Analyses Reveal Peroxidases Play Important Roles in Soybean Tolerance to Aluminum Toxicity

Comparative Analyses Reveal Peroxidases Play Important Roles in Soybean Tolerance to Aluminum Toxicity. Abstract: Aluminum (Al) toxicity is an important barrier to soybean ( Glycine max (L.) Merr.) production in acid soils. However, little is known about the genes underlying Al tolerance in soybean. We aim to ﬁnd the key candidate genes and investigate their roles in soybean tolerance to Al toxicity in this study. Comparative transcriptome analyses of the Al-tolerant (KF) and Al-sensitive (GF) soybean varieties under control and Al stress at 6, 12, and 24 h were investigated. A total of 1411 genes showed speciﬁc up-regulation in KF or more up-regulation in KF than in GF by Al stress, which were signiﬁcantly enriched in the GO terms of peroxidase (POD) activity, transporter activity (including the known Al tolerance-related ABC transporter, ALMT, and MATE), and four families of transcription factors (AP2, C3H4, MYB, WRKY). The expression levels of seven POD genes were up-regulated by Al stress for at least one time point in KF. The H 2 O 2 pretreatment signiﬁcantly improved Al tolerance of KF, which is likely due to increased POD activity induced by H 2 O 2 . Our results suggest that PODs play important roles in soybean tolerance to Al toxicity. We also propose a list of candidate genes for Al tolerance in KF, which would provide valuable insights into the Al tolerance mechanisms in soybean.


Introduction
Soybean (Glycine max (L.) Merr.) is an important legume crop, supplying high-quality oil and protein for human food, animal feed, and industrial use [1,2]. Nevertheless, Aluminum (Al) toxicity, mainly due to the phytotoxic species of Al (Al 3+ ) in acid soils (pH < 5.0), is becoming an important factor limiting crop production, including soybean, worldwide [3,4]. After exposure to Al, the inhibition of root growth would soon occur as the first visible symptom of Al toxicity [5]. Soluble Al 3+ inhibits plant root elongation rapidly by damaging the root apices [3,6,7], which will subsequently prevent the uptake and translocation of water and nutrients from roots [8]. Al disrupts calcium homeostasis [9] and induces the production of reactive oxygen species (ROS) [10], alterations in the cytoskeleton [11,12], and the formation of callose deposition in the apoplast [13], which all contribute to the Al toxicity.
Al-induced oxidative stress has been reported as a prominent symptom of Al toxicity [14,15]. Excessive ROS can give rise to oxidative damage to lipid membranes, proteins and deoxynucleic acids, which eventually leads to injured or died cells [16,17]. Antioxidant enzymes might play important roles in plant tolerance to Al stress by scavenging ROS. A study shows peroxidase (POD) expression was markedly up-regulated under Al stress, and POD activities significantly increased under Al exposure in wheat roots [18]. It has been reported that overexpression of a tobacco (Nicotiana tabacum L.) peroxidase gene (NtPox) or an Arabidopsis (Arabidopsis thaliana) POD gene AtPrx64 in Arabidopsis enhanced Al tolerance of transgenic plants [19,20]. In addition, hydrogen peroxide (H 2 O 2 ) is a key ROS that is an important signal for the activation of plant defense pathways against biotic and abiotic stresses. A previous study showed that H 2 O 2 pretreatment increased the POD activities and total antioxidant capacity, which alleviated Al-induced oxidative stress in wheat seedlings [21].
Previous studies showed that plants have developed two main mechanisms to adapt to Al stress, including Al exclusion and internal tolerance [4,[22][23][24]. Al exclusion is a well-documented and most widespread mechanism, by which Al induces the secretion of organic acid anions, such as oxalate, citrate, and/or malate from the roots to chelate Al 3+ in rhizosphere [23]. The second mechanism is internal detoxification such as the transportation of Al into the vacuoles and chelation in Al-accumulating plants [6,25].
Intraspecific variation in Al tolerance has been observed in many crop species including soybean, which is controlled by quantitative trait loci [37,[47][48][49][50][51]. Although many Al tolerance genes have been studied in plants [52], only few Al tolerance-related genes were described in soybean such as GmALMT1 [35], GmARI1 [53], and GsMATE [46]. Alresponsive genes have been identified in the Al-tolerant soybean varieties of PI416937 and Jiyu70 by transcriptome analysis [54,55]. However, a detailed comparative transcriptome analysis of soybean varieties with contrasting Al tolerance phenotypes was limited, which would be useful to identify candidate genes for Al tolerance in soybean.
In this study, the genome-wide gene expression profiling of roots from Al-tolerant (KF) and Al-sensitive (GF) soybean varieties under 0 and 25 µM AlCl 3 (pH 4.3) at 6, 12, and 24 h was investigated using Affymetrix microarray, and the gene expression patterns were verified by qRT-PCR. The metabolic pathways and candidate genes that might be involved in Al tolerance were identified by comparative transcriptome analysis. Gene Ontology (GO) and transcription factor (TF) enrichment analyses were performed to identify the important genes for Al tolerance in soybean, and the GO terms of POD activity and transporter activity were found to be significantly enriched in the up-regulated Al tolerance candidate genes by AgriGO analysis. The role of POD in Al tolerance was further investigated by comparing the changes in transcript level and enzyme activity in response to Al stress in KF and GF. This research proposed a list of candidate genes and possible mechanisms for Al tolerance in the KF soybean variety, which would provide valuable insights into the mechanisms of Al tolerance in soybean.

Plant Material, Growth Conditions, Treatments, and Tissue Harvest
The seeds of the Al-tolerant soybean cultivar KF (Kefeng 1) and Al-sensitive cultivar GF (Guangfeng Maliaodou) used in the study were obtained from the National Center for Soybean Improvement, Nanjing, China. The seeds were surface sterilized with 1% NaClO for 1 min, washed with sterile water for 5 min, then germinated in moist sterile sand and grown under a photoperiod of 14-h/10-h (day/night) and 26/24 • C (day/night) temperature circulations for three days. Then, the 3-day-old seedlings were transferred to 0.5 mM CaCl 2 (pH = 4.3) for 24 h before Al treatment. The seedlings were then exposed to 0.5 mM CaCl 2 (pH = 4.3) solution containing 0 µM AlCl 3 (control), or 25 µM AlCl 3 (Al treatment) for 6, 12, and 24 h, respectively, as described previously [56].
For H 2 O 2 pretreatment, the 3-day-old seedlings were transferred to 0.5 mM CaCl 2 (pH = 4.3) for 24 h, then pretreated by 100 µM H 2 O 2 for 6 h. After that, they were transferred to 0.5 mM CaCl 2 (pH = 4.3) solution containing 0 µM AlCl 3 (H 2 O 2 treatment), or 25 µM AlCl 3 (H 2 O 2 + Al treatment) for 6, 12, and 24 h, respectively. The lengths of the primary roots were measured with a ruler before and after the treatment. Relative root growth (RRG) is calculated following the previous study [57]. The formula is as follows: where RG t and RG c represent the root growth (RG) under treatment and control, respectively, while RL t and RL 0 represent the root length (RL) at t h (t = 6, 12, and 24 h) and 0 h, respectively.
The root tips (0-1 cm) of six seedlings were cut and pooled together in one sample and immediately frozen in liquid nitrogen and stored at −80 • C for RNA isolation. To compare the changes in POD activity after Al stress with control, uniformly grown seedlings were carefully selected and transferred into 1/2 Hoagland nutrient solution to grow for 14 days. Then 17-day-old seedlings were exposed to 0.5 mM CaCl 2 (pH = 4.3) solutions containing 0 µM (control) or 25 µM AlCl 3 . Root samples were then harvested at 6, 12, and 24 h and immediately frozen in liquid nitrogen and stored at −80 • C.

POD Activity Measurement
POD activity was measured by the previously published method [58]. Soluble proteins in soybean root samples were extracted by homogenizing the ground frozen root powder in 7 mL of 50 mM potassium phosphate buffer (pH 7.0) containing 1 mM EDTANa 2 and 1% polyvinylpyrrolidone (PVP), with the addition of 1 mM ascorbate acid (ASA) in the reaction assay. The homogenate was centrifuged at 10,000× g for 15 min at 4 • C and the supernatant was used for the following enzyme assays. The UV/visible light spectrophotometer (Qingke 722N, Shanghai, China) was used for all spectrophotometric analyses at room temperature. POD activity was determined by measuring the increase in the oxidation of guaiacol by H 2 O 2 . The reaction was performed in a final volume of 3 mL including 50 µL of the crude protein extract, 10 mM potassium phosphate (pH 6.4), 8 mM guaiacol, and 2.75 mM H 2 O 2 . The change in absorbance (A) was recorded at 470 nm after the addition of H 2 O 2 . One unit of POD activity was defined as the change in A470 nm per min (µmol min −1 mg −1 protein). Protein concentration was determined using the Bradford assay method.

Determination of Al Concentration in Soybean Roots
Roots from over 50 seedlings were pooled and collected in one sample for Al concentration measurement. The soybean roots were washed using distilled water to remove the residual Al 3+ on the surface, then dried at 85 • C for 72 h and weighed. The dried samples were digested with 5 mL nitric acid (HNO 3 ) in a microwave digestion system (ETHOS T, Milestone, Italy) following the manufacturers' recommendations. The Al 3+ concentrations were estimated by an Inductively Coupled Plasma Optical Emission Spectrometer (ICP-OES) Optima 8000 (PerkinElmer, Shelton, CT, USA).

Hematoxylin Staining of Soybean Roots
After washing by distilled water, the roots were stained with hematoxylin solution as described previously [59]. Briefly, the root tips were stained with 0.1% hematoxylin for 10 min, washed in distilled water for 30 min (10 min × 3 times), and then photographed.

Statistical Analysis
Differences between two groups or more than two groups were analyzed by Student's t-test or Duncan's new multiple range test, respectively. Figures were prepared using the Origin 9.0 software.

Total RNA Isolation, Labeling, and Array Hybridization
Total RNA was isolated from root tips (1 cm) using TRIzol (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions, and digested with DNase I (TAKARA, Shiga, Japan) to remove contaminated DNA. Then mRNAs were purified from the total RNA using Dynabeads Oligo (dT)25 (Life Technologies, Carlsbad, CA, USA). RNA concentration and purity were determined by NanoDrop 1000 (Thermo, Waltham, MA, USA) and verified by electrophoresis in 1% agarose gels.
The detailed procedure of microarray hybridization can be found in the Affymetrix GeneChip Expression Analysis Technical Manual. Soybean Gene 1.1 ST Array Strip detects 66,473 transcripts from Glycine max database and GenBank, and each transcript can be detected with 16-19 probe sets (25 bp per oligonucleotide). Briefly, mRNA (200 ng) was used to synthesize double-stranded cDNA using random primers. Then, in vitro transcriptional cDNA was generated using random primers dUTP and purified by the paramagnetic particle method. Subsequently, cDNA was hydrolyzed and purified into the first = stranded cDNA, and the first-stranded cDNA was fragmented and terminal labeled by biotin. Lastly, arrays were crossed by Hybridization Station 0381 (Affymetrix, Santa Clara, CA, USA), washed, stained in the Fluidics Station 0377 (Affymetrix, Santa Clara, CA, USA), and scanned by GeneChip ® Scanner 3000 (Affymetrix). Microarray hybridization was conducted by the CapitalBio Corporation (http://www.capitalbio.com, accessed on 26 May 2016) according to the Affymetrix standard protocols.

Microarray Data Analysis
Scanned images were transformed to digital signal by Affymetrix ® GeneChip ® Command Console ® Software, and the fluorescence signal intensity of each probe was recorded in CEL file format. Then CEL files were pre-processed by background subtraction, variance stabilization, and normalization with the RMA method [60]. Pairwise comparisons were performed between Al-treated sample and the corresponding control sample at each time-point for KF and GF, respectively, and differentially expressed genes were selected by fold change ≥2 (up-regulation) or ≤0.5 (down-regulation). The data for 12 microarrays are available in the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (accession number GSE116796). The annotation and functional categories for these transcripts were assigned based on Molecular Annotation System (MAS3.0) (http://bioinfo.capitalbio.com/mas3/, accessed on 6 August 2016) and Soybase (https://soybase.org/, accessed on 3 September 2016). The visualization of metabolism was displayed by MapMan v3.5.1 [61]. Heatmap with clustering analysis was constructed using MeV 4.9 software with the K-means method and Euclidean Distance [62]. Gene Ontology (GO) enrichment analysis was performed using AgriGO [63]. Fisher's exact test (p < 0.01) and the Benjamini and Hochberg method (FDR < 0.05) were used to determine the significance of enrichment for GO terms and the transcription factor families in the list as compared with the soybean reference genome.

Quantitative Real-Time PCR (qRT-PCR)
The first-strand cDNAs were synthesized using a PrimerScript First Strand cDNA synthesis kit (TaKaRa, Shiga, Japan) following the manufacturers' protocol, which was carried out in a total of 20 µL reaction volume including 1 µg of total RNA, 4 µL 5 × PrimeScript RT Master Mix, and RNAase-free ddH 2 O. Gene-specific primers (Table S1) were designed using primer premier 5.0 (Premier Biosoft International, Palo Alto, CA, USA) and synthesized by Invitrogen (Shanghai, China). The qRT-PCR was performed in a final volume of 15 µL containing 2 µL cDNA, 7.5 µL 2 × SYBR Premix Ex Taq (TaKaRa, Shiga, Japan), and 200 nM of forward and reverse primers, using a Roche 480 Realtime detection system (Roche Diagnostics, Basel, Switzerland) following the manufacturer's instructions. The amplification program was set as follows: initial denaturation at 95 • C for 5 min; 40 cycles of denaturation at 95 • C for 10 s, annealing at 58 • C for 20 s, and extension at 72 • C for 20 s. Each experiment was performed in triplicates. The relative expression values were calculated by the 2 −∆∆CT method [64]. The housekeeping gene GmEF-1α was used as an internal control. The relative expression level of each gene in response to Al stress (25 µM AlCl 3 ) was compared to their corresponding samples under control conditions (0 µM AlCl 3 ) at each time point. Each data point is the average ± SD of three independent replications.

Phenotypic Difference in Al Tolerance Between Soybean Varieties of KF and GF
The root lengths of two soybean varieties, KF and GF, were measured at 0, 6, 12, and 24 h after 0 (control) or 25 µM AlCl 3 treatment, and the corresponding RRG was calculated. The RRG of both varieties decreased with the duration of Al treatment, and the reduction was more obvious in GF ( Figure 1A). The RRG of KF was significantly higher than that of GF across 6, 12, and 24 h. The concentration of Al in soybean roots was detected by ICP after 6, 12, and 24 h of Al treatment. The roots of GF accumulated significantly more Al than those of KF at 6 and 24 h ( Figure 1B). The results confirmed that KF is more tolerant to Al toxicity than GF, and therefore we used T (tolerant) and S (sensitive) to represent KF and GF, respectively.

Genome-Wide Gene Expression Profiles in Soybean Roots in Response to Al Stress
Overall, 65,535 out of the 66,473 transcripts in the Affymetrix Soybean Gene 1.1 ST Array were detected in this study. The dynamic changes in gene expression profiles at 6, 12, and 24 h in the roots of T and S under control and Al treatment could be observed by the heatmap with clustering analysis ( Figure S1). The graph showed that 12 samples were classified into two clusters, including one cluster with six control samples and the other cluster with Al-treated samples. The transcripts were also classified into two groups, with one group showing lower expression levels in Al-treated samples than the control samples, while the other group showed the reverse pattern. The transcriptional changes of genomewide genes between Al stress and control at 6, 12, 24 h in two soybean varieties could be visualized using MapMan software ( Figures S2-S4). There were more transcriptional changes at 24 h than 6 and 12 h after Al stress. At 24 h, more down-regulated genes were noticed in S than T, especially in the metabolic pathways involved in glycolysis and tricarboxylic acid cycle (TCA) ( Figure S4). In the glycolysis pathway, there were three key glycolytic enzymes, including hexokinase, phosphofructokinase, and pyruvate kinase. The genes encoding phosphofructokinase and pyruvate kinase were found more downregulated in S than T by Al stress (Table S2). Most genes in TCA were down-regulated in S, while only few of them decreased in T ( Figure S4). Several key genes in TCA such as genes encoding phosphoenolpyruvate carboxylase (PEPC), isocitrate dehydrogenase, succinate dehydrogenase, and fumarate lyase were more down-regulated in S than T (Table S2).

Identification of Al-Responsive Genes in Soybean Roots
Al-responsive genes were identified by comparing Al treatment and control at 6, 12, and 24 h in T and S, respectively (Table S3). The numbers of up-or down-regulated genes in response to Al stress increased with the duration of treatment time. At 6 h, the number of up-regulated genes by Al stress was 759 and 850 in T and S, respectively. At 12 h, the number of up-regulated genes was 1233 and 1356 in T and S, respectively. At 24 h, the number of up-regulated genes was 2240 and 2908 in T and S, respectively. The number of down-regulated genes also showed a trend of increase with the duration of Al treatment. After 6 h of treatment, the number of down-regulated genes in response to Al was 451 and 544 in T and S, respectively. After 12 h treatment, the number of down-regulated genes was 613 and 1560 in T and S, respectively. After 24 h treatment, the number of down-regulated genes was 1501 and 3725 in T and S, respectively. Much more genes were suppressed in S than T at 24 h after Al stress, and a greater number of differentially expressed genes were found in S than T across all three time points.

Identification of Candidate Genes for Soybean Tolerance to Al Stress
Many known Al tolerance genes were up-or down-regulated in response to Al stress, and their fold changes in Al-tolerant plants were more than those of Al-sensitive plants.
Our results showed that 1417 probes were up-regulated by Al stress only in T but not in S, and 520 probes showed more up-regulation by Al in T than S when comparing the same time point. These 1937 (1417 + 520) probes represent 1460 genes (due to multiple probes for some genes), and 1411 of them have the correspondence in newest version of soybean reference genome (Wm82a2.v1) with gene annotations (Table S4). On the other hand, 925 probes were down-regulated only in T but not S, and 57 probes showed more downregulation by Al in T than S when comparing the same time point. These 982 (925 + 57) probes represent 671 genes (due to multiple probes for some genes), and 633 of them have the correspondence in newest version of soybean reference genome (Wm82a2.v1) with gene annotations (Table S5). These 1411 and 633 genes that showed different response to Al stress between T and S could be the potential candidate genes for soybean tolerance to Al toxicity.

Gene Ontology (GO) Analysis of Al Tolerance Candidate Genes in Soybean
Gene Ontology (GO) classification was performed using Soybase and WEGO ( Figure S5). The 1411 Al-induced genes and 633 Al-repressed genes were classified into 47 functional categories including 10 cell components, 14 molecular functions, and 23 biological processes, respectively ( Figure S5A,B).
Next, we performed AgriGO functional enrichment analysis of these Al tolerance candidate genes. The GO term of peroxidase (POD) activity (GO:0004601, FDR = 7.35 × 10 −9 , 34 genes) in molecular function was identified as the most significantly enriched group in the Al-induced genes ( Figure 2, Table S6). No GO term in molecular function was found significantly enriched among the Al-repressed genes. Al tolerance in soybean. The gene ratio on the x-axis represents the percentage of genes in each category. The Singular Enrichment Analysis (SEA) was performed using AgriGO. The GO terms and significance levels are displayed. Only the significant (p < 0.01, FDR < 0.05) GO terms are displayed. FDR: False discovery rate, which is shown as the logo on the right ranging from the lowest FDR (in red) to highest FDR (in blue).
There are 142 genes encoding transcription factors in the 1411 Al-induced genes (Table S7). According to the putative annotation in the Plant Transcription Factor Database 4.0 [65], these 142 transcription factors representing 20 different families, and four families (AP2, C3H4, MYB, WRKY), were significantly enriched in the 1411 Al tolerance candidate genes compared with the transcription factors in the soybean genome background ( Figure S6).

Validation of Microarray Data by Real-Time Quantitative PCR
Ten genes were selected from the 1411 Al-induced genes (Table S4) to validate the expression patterns of the microarray data and the significantly enriched GO terms of peroxidase (POD) activity and transporter activity, including nine POD genes (Table S6) and one gene encoding ABC transporter. The relative expression levels (log 2 fold change) of ten genes in the roots of T and S from the microarray and qRT-PCR were compared, which generally showed a similar expression pattern by the two methods ( Figure S7A,B). There was a significant correlation (r = 0.86, p < 0.01) between the microarray and qRT-PCR results ( Figure S7C), and the coefficient of determination (R 2 ) is 0.74, which is comparable to the other studies where the R 2 is~0.72 [66,67].

Role of Peroxidase in Soybean Tolerance to Al Stress
Since the GO term of POD activity was identified as the most significantly enriched group in the Al-induced candidate genes for Al tolerance in soybean, we performed qRT-PCR to confirm the expression patterns of seven POD genes (randomly chosen from Table S6) in response to Al stress (Figure 3). The relative expression of these seven POD genes were up-regulated by Al stress in T and showed a significant difference in response to Al between T and S for at least one time point (Figure 3). Next, we determined the POD activity in the roots of two soybean varieties (Figure 4). The POD activity in the roots of T increased during Al treatment time, and there was a significant difference between Al and the control condition at 6 and 24 h. In the roots of S, the POD activity only showed a significant increase at 24 h after Al stress compared with the control. The POD activity in the roots of T was generally higher than that of S under Al stress (Figure 4).  Next, we investigated if H 2 O 2 pretreatment alleviated Al toxicity in soybean ( Figure 5). The results showed that H 2 O 2 pretreatment improved Al tolerance (RRG) and reduced Al concentration in soybean primary roots of both T and S (S only showed significantly improved RRG at 6 h) compared with Al treatment alone, while H 2 O 2 treatment (without Al) did not have significant effect on RRG or Al concentration in roots ( Figure 5A,B). Moreover, hematoxylin staining of the root tips also showed that Al accumulation (reflected by a blue color) in soybean root tips was decreased by H 2 O 2 pretreatment and was more obvious in T than S ( Figure 5C). In addition, the POD activity was significantly increased by 6 h of H 2 O 2 pretreatment in the roots of T ( Figure S8). The results suggest that H 2 O 2 pretreatment increased the POD activity in soybean roots, which might alleviate Al-induced oxidative stress, leading to improved Al tolerance (RRG) of soybean roots.

Al Tolerance Related Transporter Genes in Soybean Roots
Among the 1411 Al tolerance candidate genes identified above, we found that genes with the molecular function of transporter activity (GO:0005215) were also significantly enriched (FDR = 3.09 × 10 −3 , Figure 2). Previous studies have shown that ALMT, MATE, and ABC transporters play important roles in Al resistance of many plant species. In this study, nine genes encoding ABC transporters, one ALMT gene, and four MATE genes were up-regulated by Al stress and with more fold changes in T than S for at least one time point ( Figure 6). One MATE gene, Glyma.13G203000, which has been previously designated as GmMATE75 [56], was significantly up-regulated by Al stress in both T and S and with more fold changes in T than S at all three time points including at 6, 12, and 24 h ( Figure 6). Therefore, GmMATE75 might play an important role in soybean tolerance to Al toxicity.

Discussion
Al toxicity in acid soils is one of the primary environmental challenges limiting crop productivity. To identify Al-responsive genes and characterize the molecular mechanisms underlying Al tolerance, transcriptomic analysis has been performed in Arabidopsis, buckwheat (Fagopyrum tataricum), common bean (Phaseolus vulgaris), maize, rice, and soybean response to Al [55,66,[68][69][70][71]. Previous studies showed that there was intraspecific variation in Al tolerance in many crop species [37,47,50]. In this study, two soybean varieties, KF (T) and GF (S), showed a significant difference in Al tolerance. Therefore, we compared the genome-wide expression profiles of these two soybean varieties in response to Al stress to identify the candidate genes for soybean Al tolerance and investigated the roles of PODs in soybean tolerance to Al toxicity.
The genome-wide transcriptional changes of two soybean varieties, T and S, in response to Al stress were analyzed, and the metabolic pathways involved in T and S soybean response to Al stress were compared ( Figure S2-S4). At 24 h, more down-regulated genes involved in glycolysis and tricarboxylic acid cycle (TCA) were found in S than in T ( Figure S4). It has been suggested that the operation of TCA in response to stress, coupled with the citrate pathway that enhances Al tolerance, may lead to modified metabolic homeostasis and/or improvements in productivity and overall cellular physiology in plant response to Al [72]. In this study, genes encoding PEPC in TCA (Glyma.06g229900, Glyma.06g277500, Glyma.13g270400) were more down-regulated in S than T. Plants overexpressing PEPC showed enhanced oxalate exudation [73], while oxalate is a strong chelator of Al. Therefore, the role of PEPC in soybean tolerance to Al stress should be further investigated.
GO enrichment analysis revealed that the GO term of POD activity in the molecular function was significantly enriched in the Al-induced candidate genes for Al tolerance compared with the genome-wide genes ( Figure 2). Hence, we performed qRT-PCR to determine the expression patterns of POD genes in response to Al stress. The expression levels of seven POD genes in the roots of T were up-regulated by Al and with more fold changes in T than in S for at least one time point (Figure 3). Moreover, the POD activity was also significantly increased by Al stress. These results indicate that POD is likely to play an important role in soybean tolerance to Al toxicity. A previous study has showed that Al stress induced ROS production, and higher levels of POD activity were observed in the Al-tolerant rice variety than in the Al-sensitive rice variety [74].
H 2 O 2 is considered one of the important members of ROS that is an important signal in the activation of plant defense pathways against biotic and abiotic stresses. Previous studies found that exogenous application of H 2 O 2 was able to improve the Al tolerance of wheat and tobacco [21,75]. This study showed that 6 h of H 2 O 2 pretreatment decreased the Al content in soybean root tips and enhanced Al tolerance (RRG) of soybean ( Figure 5). A previous study found that Al toxicity induced massive endogenous H 2 O 2 accumulation in soybean roots [76], which might act as signals for the activation of the expression of ROS scavenging enzymes such as PODs to alleviate the oxidative damage and root inhibition caused by Al stress.
We also found that four families of transcription factors (AP2, C3H4, MYB, WRKY) were significantly enriched in the 1411 Al-induced candidate genes for Al tolerance ( Figure S6). The AP2 family is well known for its important functions in plant growth and development, hormonal regulation, and environmental stress response [77,78], but the role of AP2 in Al tolerance is unknown. MYB transcription factors possess diverse functions in secondary metabolism, hormone signal transduction, and response to environmental stresses. A previous study found that the targets of Al-responsive miRNA include MYB transcription factor mRNA in wild soybean [79]. The role of MYB transcription factors in Al tolerance needs to be further investigated.
Plant WRKY transcription factors are known to be involved in biotic and abiotic stress responses. Previous studies found that WRKY are crucial regulators of defense response and disease resistance [80,81]. It has been reported that AtWRKY46 negatively regulated AtALMT1, and the wrky46 mutant increased the release of malate and decreased the accumulation of Al 3+ which resulted in enhanced Al tolerance in Arabidopsis [82]. In this study, 15 WRKY genes were identified in the Al-induced candidate genes for Al tolerance, and the expression of GmWRKY70 (Glyma.13g267600) was significantly induced by Al at 6, 12, and 24 h in T, but with almost no changes in S. This pattern was different from the AtWRKY46 expression pattern in response to Al stress in Arabidopsis. Thus, GmWRKY70 might play a role in soybean tolerance to Al toxicity but differed from AtWRKY46.
In addition, transporters such as the ABC transporter, ALMT, and MATE, have been reported to play important roles in plant tolerance to Al toxicity [83]. In Arabidopsis, AtALS3 encoding an ABC transporter is required for Al tolerance [84]. OsSTAR1 and OsSTAR2 encode a nucleotide binding domain and a transmembrane domain of an ABC transporter, respectively, and are required for detoxification of Al in rice [83,85]. The ALMT and MATE genes encode membrane transporters that control malate and citrate efflux in response to Al, respectively [83]. The function of ALMT genes have been studied in Arabidopsis, barley, rape, rye, and wheat [32,34,38,86,87]. The overexpression of GmALMT1 enhanced malate secretion and Al tolerance of soybean hair roots [35]. The role of MATE genes in Al tolerance has also been reported in many crops such as rice, sorghum and wild soybean [36,37,46]. In the present study, the expression levels of nine genes encoding ABC transporters, one ALMT, and four MATE genes were induced by Al and showed more up-regulation in T than in S for at least one time point (Figure 6), indicating their potential roles in soybean tolerance to Al toxicity. One MATE gene, Glyma.13G203000 (GmMATE75), was significantly up-regulated at 6, 12, and 24 h of Al stress and showed more fold changes in T than S in the present study and also in our previous study [56]. Overexpression of GsMATE (a homolog gene of GmMATE75 from wild soybean) enhanced the Al tolerance of transgenic Arabidopsis [46]. A recent study reported that GmMATE75 is a plasma-membrane-localized citrate transporter, and its overexpression led to increased citrate efflux and less Al content in soybean hairy roots [88], suggesting GmMATE75 plays important roles in soybean tolerance to Al toxicity.

Conclusions
In summary, the genome-wide gene expression profiling of roots from Al-tolerant (T) and Al-sensitive (S) soybean varieties under control and Al stress at 6, 12, and 24 h was investigated, and 1411 genes showed specific up-regulation in T or more up-regulation in T than in S by Al stress, which were significantly enriched in the GO terms of POD activity, transporter activity, and four families of transcription factors (AP2, C3H4, MYB, WRKY). The POD activities in the roots of T and S increased after Al stress. The H 2 O 2 pretreatment significantly improved Al tolerance of T, which is likely due to the increased POD activity induced by H 2 O 2 . Fourteen genes encoding transporters that were related to Al tolerance, such as the ABC transporter, ALMT, and MATE, were up-regulated by Al and with more fold changes in T than in S for at least one time point. Our results suggested that POD and transporters play important roles in soybean tolerance to Al toxicity.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/agronomy11040670/s1, Figure S1: Hierarchical cluster analysis of the genome-wide transcripts in the roots of soybean varieties T and S under control (CK) or Al stress. The hierarchical cluster analysis was performed using complete method with Euclidean distance metric by Average linkage in Cluster 3.0. The gene expression levels were transformed by log 10 . Figure S2: The transcriptional changes of genome-wide genes in response to Al stress at 6 h in two soybean varieties, KF (T) and GF (S), by MapMan overview of metabolism. The log 2 (fold change) of genome-wide genes between Al stress and control in two soybean varieties could be visualized using the MapMan software. The metabolic pathways were shown by different shapes with a group of genes in small square boxes. The red and blue square boxes indicate the genes were up-regulated or down-regulated in response to Al, respectively. Figure S3: The transcriptional changes of genome-wide genes in response to Al stress at 12 h in two soybean varieties, KF (T) and GF (S), by MapMan overview of metabolism. The log 2 (fold change) of genome-wide genes between Al stress and control in two soybean varieties could be visualized using the MapMan software. The metabolic pathways were shown by different shapes with a group of genes in small square boxes. The red and blue square boxes indicate the genes were up-regulated or down-regulated in response to Al, respectively. Figure S4: The transcriptional changes of genome-wide genes in response to Al stress at 24 h in two soybean varieties, KF (T) and GF (S), by MapMan overview of metabolism. The log 2 (fold change) of genome-wide genes between Al stress and control in two soybean varieties could be visualized using the MapMan software. The metabolic pathways were shown by different shapes with a group of genes in small square boxes. The red and blue square boxes indicate the genes were up-regulated or down-regulated in response to Al, respectively. Figure Figure S6: Transcription factor family enrichment analysis of the 1411 up-regulated candidate genes for Al tolerance in soybean. The percentage of transcription factor (TF) families in the test set and all TFs in the soybean genome reference set were shown by a green bar and grey bar, respectively. * indicates significant enrichment of the TF family (p < 0.01, FDR < 0.05). Figure S7: Comparison of the results from microarray data and qRT-PCR analysis. Log 2 fold change of ten genes in the roots of T after exposure to Al stress for 6, 12, and 24 h (A); log 2 fold change of ten genes in the roots of S after exposure to Al stress for 6, 12, and 24 h (B); and the correlation of the log 2 (fold change) obtained by microarray (x-axis) with log 2 (fold change) obtained by qRT-PCR (y-axis) (C). Figure S8: The effect of H 2 O 2 treatment on POD activity in soybean roots. POD activity in the roots from 4-d-old soybean seedlings of T and S exposed to 100 µM H 2 O 2 treatment for 6 h and control (CK). Means and standard errors were calculated from three biological replications, 12 plants per replication. ** indicates significant difference (student's t-test) in POD activity between CK and H 2 O 2 treatment at a significance level of 0.01. Table S1: Primers used for quantitative RT-PCR. Table S2: Genes in glycolysis and tricarboxylic acid cycle pathways in response to Al at 24 h in the roots of KF (T) and GF (S). Table S3: Number of differently expressed genes between Al stress and control. Table S4: Differentially up-regulated genes by Al stress between Al-tolerant soybean T and Al-sensitive soybean S. Table S5: Differentially down-regulated genes by Al stress between Al-tolerant soybean T and Al-sensitive soybean S. Table S6: The 34 genes in the enriched GO term of peroxidase activity (GO:0004601). Table S7: Genes encoding transcription factors in the 1411 Al tolerance candidate genes.

Data Availability Statement:
The data presented in this study are available in the supplementary material, and the microarray data have been deposited in the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (accession number GSE116796).

Conflicts of Interest:
The authors declare no conflict of interest.