Complex Gene Regulation Underlying Mineral Nutrient Homeostasis in Soybean Root Response to Acidity Stress

Proton toxicity is one of the major environmental stresses limiting crop production and becomes increasingly serious because of anthropogenic activities. To understand acid tolerance mechanisms, the plant growth, mineral nutrients accumulation, and global transcriptome changes in soybean (Glycine max) in response to long-term acidity stress were investigated. Results showed that acidity stress significantly inhibited soybean root growth but exhibited slight effects on the shoot growth. Moreover, concentrations of essential mineral nutrients were significantly affected by acidity stress, mainly differing among soybean organs and mineral nutrient types. Concentrations of phosphorus (P) and molybdenum (Mo) in both leaves and roots, nitrogen (N), and potassium (K) in roots and magnesium (Mg) in leaves were significantly decreased by acidity stress, respectively. Whereas, concentrations of calcium (Ca), sulfate (S), and iron (Fe) were increased in both leaves and roots. Transcriptome analyses in soybean roots resulted in identification of 419 up-regulated and 555 down-regulated genes under acid conditions. A total of 38 differentially expressed genes (DEGs) were involved in mineral nutrients transportation. Among them, all the detected five GmPTs, four GmZIPs, two GmAMTs, and GmKUPs, together with GmIRT1, GmNramp5, GmVIT2.1, GmSKOR, GmTPK5, and GmHKT1, were significantly down-regulated by acidity stress. Moreover, the transcription of genes encoding transcription factors (e.g., GmSTOP2s) and associated with pH stat metabolic pathways was significantly up-regulated by acidity stress. Taken together, it strongly suggests that maintaining pH stat and mineral nutrient homeostasis are adaptive strategies of soybean responses to acidity stress, which might be regulated by a complex signaling network.


Introduction
Acid soils with pH lower than 5.5 comprise up to approximate 30% of total areas of the planet and one-half of world's potentially arable lands [1][2][3][4]. In the tropics and subtropics, soil acidification is a natural and very slow process. However, anthropogenic activities, such as inappropriate application of fertilizers and massive use of fossil energy sources, have accelerated soil acidification in recent time [5][6][7][8][9].
In acid soils, H + rhizotoxicity is recognized as a major limiting factor for crop production [2]. Field experiments in Northern Idaho showed that yield of important crops decreased on even moderately low-pH soils without Al toxicity [10]. It has been suggested that physiological response of plant roots to low-pH stress could occur in a very short time period. For example, root cells of With the aid of genome wide transcriptomic analysis and functional characterization of sets of genes, molecular mechanisms underlying complex responses of plants to acidity stress have been elucidated [33][34][35]. A microarray assay on wheat (Triticum aestivum) identified 1057 differentially expressed genes (DEGs) after two hours of acid treatment, in which rhizosphere alkalification was suggested to be the earliest responses to acidity stress [34]. Moreover, longer acid treatment (8 h) on Arabidopsis identified a total of 881 genes exhibited at least 2-fold changes in transcripts, suggesting redundancies and interactions among the responses to pH, auxin and pathogen elicitors [33]. Interestingly, a recent RNA-seq analysis in soybean found a total of 3242 genes were up-regulated by pH3.0 medium for 9 days. Among them, 27 genes were identified to be involved in the biosynthesis of glycollins, which are isoflavonoid-derived pathogen-inducible defense metabolites [35]. Moreover, genes function in systemic acquired resistance (SAR), which is a component of the plant immune system to pathogen, were also significantly up-regulated [35]. These results indicated that acidity stress is an elicitor of plant defensive responses to pathogens. Except for transcriptomic assays, complex regulatory mechanisms underlying plant acid tolerance have been widely documented through identification and functional analysis of STOP1 (Sensitive to Proton Rhizotoxicity 1) and its down-stream genes such as GAD (Glutamate Decarboxylase), ME (Malic Enzyme) and GDH (Glutamate Dehydrogenase) [36][37][38][39][40][41][42][43]. However, most of these studies mainly focus on the mechanisms underlying plant root growth adaptation to acidity stress, few studies have tried to investigate the effects of acidity stress on mineral nutrient acquisition and translocation in plants, especially in crops.
Soybean (Glycine max) is one of the most economically important leguminous crops, which provides more than a quarter of the world's protein for food and animal feed [44]. Approximate 2000 years ago, soybean cultivation was expanded worldwide, and domesticated with wide adaptation to different soil types, especially acid soils [45,46]. For instance, in Brazil, the world's second largest soybean producer, soybean is mainly grown on savannahs regions, which comprises prevalence of acid infertile soils with low soil pH, high aluminum (Al) availability, and low mineral nutrient-holding capacity [47,48]. Although intensive studies have been conducted to investigate soybean adaptation to Al toxicity and mineral nutrient deficiencies (e.g., phosphorus and nitrogen) on acid soils [42,[49][50][51][52][53][54], few studies have been carried out to elucidate complex mechanisms underlying soybean adaptation to proton toxicity through integration of mineral nutrients accumulation and transcriptomic analyses. Therefore, the present study investigated the effects of acidity stress on soybean growth and mineral nutrients accumulation. Subsequently, transcriptome analysis was conducted to identify DEGs in soybean roots responding to acidity stress, which highlighted transcriptional regulatory mechanisms underlying its adaptive strategies, especially for mineral nutrient acquisition and translocation.

Plant Material and Growth Conditions
Soybean genotype Huachun 6 was used in the study. Soybean seeds were surface-sterilized with 10% (v/v) H 2 O 2 , then directly sown in the prepared peat soils (JIFFY, Hoek van Holland, The Netherlands) moistened with nutrient solution, of which soil pH was adjusted to 4.2 and 5.8 by addition of diluted H 2 SO 4 . After seed germination, seedlings were watered with deionized water every day and supplied with 200 mL nutrient solution with pH of 4.2 and 5.8 once a week. The nutrient solution was composed of (in mM): For transcriptome analysis, soybean seeds were surface-sterilized and germinated in paper rolls moistened with half-strength nutrient solution as mentioned above. After germination, uniform seedlings were transplanted to hydroponic culture, the component of which was the same as mentioned above with pH 5.8 and pH 4.2, respectively. The solution pH was adjusted every day with diluted KOH or H 2 SO 4 and renewed once a week. After 30 days, soybean roots were harvested for transcripteomic assay. Each replicate had two plants and was used for one cDNA library construction. Similarly, roots of other two plants in each replicate were harvested for quantitative real-time PCR (qRT-PCR) analyses. Meanwhile, dry weight, leaf number, shoot height, branch number, along with total root length, average root diameter, root surface area and mineral nutrient concentrations in roots of another two plants in each replicate were determined. All experiments had three replicates.

Determination of Mineral Nutrient Concentration
Measurement of N and P concentrations in soybean samples was conducted using the continuous flowing analyzer (Series SA1100, SKALAR, Breda, The Netherlands). Briefly, about 0.1 g dry samples were incubated with 1 mL H 2 O 2 and 5 mL H 2 SO 4 overnight. The processed samples were then heating digested and constant-volumed to 50 mL. The resulted products were analyzed following the instruction in the manuals.
The other dry samples were digested using the START D microwave digestion system (Series: 131825, Milestone, Sorisole, Italy) following the instruction in the manuals. Briefly, samples were incubated with 8 mL diluted HNO 3 and 2 mL 30% H 2 O 2 for overnight. The resulting products were applied to microwave digestion following the consecutive procedures, in which the temperature was raised to 130°C within 10 min, then hold for 5 min at 130°C. After that, the temperature was raised to 180°C within 10 min, and then hold for 30 min at 180°C and ventilation for 30 min. After microwave digestion, mixtures were moved to plastic volumetric flasks and adjusted to 100 mL, then filtered with filter papers. The products were subsequently used for determination of potassium (K), calcium (Ca), magnesium (Mg), iron (Fe), molybdenum (Mo), zinc (Zn), copper (Cu), and manganese (Mn) by atomic absorption spectrometer (Series: ZA-3300, Hitachi, Tokyo, Japan), and sulfate (S) by inductively couple plasma atomic emission spectrometer (Series: 710-ES, Agilent, Santa Clara, CA, USA).

Transcriptome Analysis of Soybean Roots
Total RNA was extracted from soybean roots under acid (pH 4.2) and normal (pH 5.8) conditions using an RNeasy Mini Kit (Qiagen, Hilden, Germany). Then, the total RNA was digested with RNase free DNAase I (Qiagen, Hilden, Germany) for the purification of mRNA using Oligo (dT) magnetic beads.
Purified mRNA was broken into short fragments and used for the cDNA synthesis. Short fragments were purified with a QiaQuick PCR Extraction Kit (Qiagen, Hilden, Germany) and connected with sequencing adapters. The proper products were selected for PCR amplification and cDNA library construction. The library was sequenced using BGISEQ500 (The Beijing Genomics Institute, Beijing, China). After sequencing, the low-quality, adaptor-polluted and high content of unknown base reads were filtered from the raw data. The clean reads were mapped to soybean reference genome (Glycine max Williams82.a2.v1) to predict novel transcripts. Subsequently, the coding novel transcripts were merged to get a complete reference, which was then used for the analysis of differentially expressed genes (DEGs). During the processes, software including HISAT, StringTie, Bowtie2 and RSEM was used to perform genome mapping, reconstruct transcripts, map clean reads to reference and calculate gene expression level, respectively [55][56][57]. Finally, DEseq2 was used to identify DEGs by different treatments with p value less than 0.05, and log 2 Fold change 1 or-1 as the threshold. The transcriptomic data was up-loaded to GEO repository (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE129320) with number of GSE129320.

Quantitative Real-Time PCR Analysis
Total RNA was extracted from soybean roots using the RNA-solve reagent (OMEGA Bio-tek, Norcross, GA, USA). After being treated with RNase-free DNase I (Invitrogen, Carlsbad, CA, USA), total RNA was used to synthesize the first-strand cDNA using MMLV-reverse transcriptase (Promega, Madison, WI, USA) following instructions in manuals. The qRT-PCR analysis was conducted using SYBR Green-monitored on a ABI StepOne Plus real-time PCR system (Thermo Fisher Scientific, Waltham, MA, USA). Primer pairs (Table S1) for qRT-PCR were designed. Housekeeping gene GmEF-1a (accession no. X56856) was used as an internal control to normalize the expression of corresponding genes in soybean.

Statistical Analysis
All data were analyzed by Student's t-test using SPSS software (IBM SPSS Statistics, Chicago, IL, USA).

Effects of Low-pH Stress on Soybean Growth in Soils
Effects of low-pH stress on soybean growth were investigated in soils after 30 days treatment. Results showed that compared to control, soybean root dry weight was decreased by 26% under acid conditions ( Figure 1C). Consistently, total root length and root surface area were also significantly decreased by low-pH stress, as reflected by 18% and 17% reduction, respectively ( Figure 1G,H). In contrast, no significant differences were observed in shoot dry weight, leaf number, shoot height, branch number and average root diameter between two pH-treatments ( Figure 1). The results strongly suggest that soybean roots are more sensitive to low-pH stress than shoots.

Concentration of Mineral Nutrients Affected by Acidity Stress
Concentrations of six mineral macronutrients (i.e., N, P, K, Ca, Mg and S) were separately determined in leaves and roots of soybean grown on soils. It was found that concentrations of the mineral macronutrients were affected by low-pH stress, mainly differing among plant organs and types of mineral nutrients ( Figure 2). For N, P, and K, their concentrations in soybean roots were separately decreased by 21%, 38% and 19% by low-pH treatment compared to the control. However, only P concentration was decreased in leaves, while neither N or K concentration in soybean leaves was found to be significantly affected by low-pH stress (Figure 2A-C). Therefore, the results indicate

Concentration of Mineral Nutrients Affected by Acidity Stress
Concentrations of six mineral macronutrients (i.e., N, P, K, Ca, Mg and S) were separately determined in leaves and roots of soybean grown on soils. It was found that concentrations of the mineral macronutrients were affected by low-pH stress, mainly differing among plant organs and types of mineral nutrients ( Figure 2). For N, P, and K, their concentrations in soybean roots were separately decreased by 21%, 38% and 19% by low-pH treatment compared to the control. However, only P concentration was decreased in leaves, while neither N or K concentration in soybean leaves was found to be significantly affected by low-pH stress (Figure 2A-C). Therefore, the results indicate that acidity stress might affect the absorption of K and N, but not their translocation from roots to shoots.
Concentrations of three other macronutrients (i.e., Ca, Mg and S) were also analyzed in both leaves and roots. Concentrations of Ca and S in both leaves and roots were significantly increased under low-pH conditions ( Figure 2D,E). Compared to the control, low-pH treatment increased Ca concentration by 20% in leaves and 35% in roots, respectively ( Figure 2D). Similarly, S concentration was increased by 3 folds in leaves and 2 folds in roots, respectively ( Figure 2E). However, low-pH led to significant decreases of Mg concertation in leaves, as reflected by 15% reductions, but not in roots ( Figure 2F). These results strongly suggest that low pH exhibits different effects of mineral macronutrient acquisition and translocation in soybean.
Concentrations of five mineral micronutrients, including Fe, Mo, Mn, Zn and Cu, were also determined in soybean leaves and roots at two pH-treatments. It showed that low-pH significantly affected concentrations of Fe, Mo, and Mn, but just slightly influenced the concentrations of Zn and Cu ( Figure 2G-K). Under acid conditions, Fe concentration was increased by about 2 folds in both leaves and roots ( Figure 2G). Similarly, Mn concentration was increased by about 5 folds in roots by low-pH treatment, whereas, no significant difference was observed in leaves ( Figure 2I). In contrast, Mo concentration was significantly reduced under acid conditions, as reflected by 71% decreases in leaves and 64% decreases in roots, respectively ( Figure 2H).

Growth and Mineral Nutrient Concentrations in Roots of Soybean Grown in Nutrient Solution
Soybean seedlings were further subjected to hydroponic culture with two pH-treatments (pH 5.8 and pH 4.2) to determine effects of acidity stress on soybean growth and the mineral nutrient concentration in roots. It showed that similar to the results in soil culture, the acidity stress significantly inhibited soybean root growth as indicated by reduced root dry weight and total root length, which were 31% and 19% decreased under acid condition, respectively (Figure S1C,G). However, no significant differences were observed for the shoot dry weight, plant height, leaf number, branch number, root surface area and average root diameter between two pH-treatments ( Figure S1). Mineral nutrient concentrations in soybean roots were also measured. Results showed that concentrations of N, P, K and Mo in acid-treated (pH 4.2) soybean roots were 15%, 23%, 28% and 34% less than those of the control (pH 5.8), respectively. In contrast, concentrations of Ca, S and Mn in acid-treated soybean roots were 1.3, 1.5 and 1.9 folds higher than those of control, respectively. These results are similar with that in soil culture. However, no significant differences were observed for concentrations of Fe, Mg, Zn and Cu between two pH-treatments ( Figure S2).

Transcriptomic Analysis of Soybean Roots Responding to Acidity Stress
In order to investigate molecular mechanism underlying soybean adaptive strategies, transcriptomic analysis was conducted in soybean roots grown in the nutrient solution at two pH-treatments. For each treatment, three libraries were constructed and produced around 53 million raw reads for each treatment (Table S2). After excluding the low-quality readings, about 44.8 and 44.5 million clean reads were obtained for libraries of control and acid-treated soybean roots, respectively (Table S2). Total number of transcribed genes in control and acid-treated roots were 66,391 and 65,400, respectively (Table S2). Moreover, transcripts resulted in a total of 41,055 and 40,476 known-genes in control and acid-treated roots, respectively (Table S2).
The DEGs between control and acid-treated roots were defined as those with at least 2-fold changes in expression, as well as p value 0.05. Using the criteria, a total of 974 genes were considered to be differently expressed in soybean roots between two pH-treatments (Table 1 and Table S3). Among them, 419 genes were up-regulated, while 555 genes were down-regulated by low-pH stress (Table 1 and  Table S3). Gene ontology (GO) category analysis showed that 974 DEGs could be divided into 18 biological processes, 14 cellular components, and 11 molecular function terms ( Figure S3). The pathway functional enrichment of DEGs showed that most DEGs were predicted to function in the metabolic pathway ( Figure S4).

Transcription Factors Involved in Soybean Responses to Acidity Stress
A total of 79 DEGs were predicted to encode transcription factors (Figure 3, Tables S4 and S5). Among them, genes encoding AP2-EREBP transcription factors were the most abundant, with twelve up-regulated and two down-regulated. Other DEGs encoding transcription factors included ten MYBs, nine bHLHs, eight WRKs, seven GRASs, four G2-likes and NACs, three C2H2s, HSFs, PLATZs, and RWP-RKs, two C2C2-Dofs and C3Hs, as well as one MADS, LOB, TCP, ARR-B, SBP, Trihelix, and ZF-HD (Figure 3 and Table S4).
Genes 2019, 9, x FOR PEER REVIEW 9 of 21 contrast, one of MYB genes, Glyma.02G177800, was predicted to encode a PHR1-LIKE (PHL) transcription factor, whose expression level was decreased by 83% under acid conditions (Table S4).

Identification of DEGs Involved in the pH Stat Pathway
A set of DEGs was predicted to be involved in the pH homeostasis process in plant cells, namely pH stat pathway [58]. The genes included one glutamate dehydrogenases (i.e., GmGDH1), two malic enzymes (i.e., GmME1-1/2), three pyruvate decarboxylases (i.e., GmPCD1 and GmPCD2-1/2) and four alcohol dehydrogenases (i.e., GmADH1-1/2/3 and GmADH4) ( Table 2). Most of these genes were up-regulated by acidity stress, except for GmME1-1 and GmADH4, whose expression was suppressed by 2.8 and 2.5 times under acid conditions. Among the up-regulated DEGs, transcript level of GmGDH1, which was involved in the GABA shunt pathway, was increased by about 3.5 folds (Table 2). Moreover, expression level of GmME1-2 and GmADH1-2 was separately increased by about 5.3 and 6.5 folds (Table 2). These results indicate that pH stat pathway is significantly influenced in soybean roots by acidity stress.

Identification of DEGs Involved in the pH Stat Pathway
A set of DEGs was predicted to be involved in the pH homeostasis process in plant cells, namely pH stat pathway [58]. The genes included one glutamate dehydrogenases (i.e., GmGDH1), two malic enzymes (i.e., GmME1-1/2), three pyruvate decarboxylases (i.e., GmPCD1 and GmPCD2-1/2) and four alcohol dehydrogenases (i.e., GmADH1-1/2/3 and GmADH4) ( Table 2). Most of these genes were up-regulated by acidity stress, except for GmME1-1 and GmADH4, whose expression was suppressed by 2.8 and 2.5 times under acid conditions. Among the up-regulated DEGs, transcript level of GmGDH1, which was involved in the GABA shunt pathway, was increased by about 3.5 folds (Table 2). Moreover, expression level of GmME1-2 and GmADH1-2 was separately increased by about 5.3 and 6.5 folds ( Table 2). These results indicate that pH stat pathway is significantly influenced in soybean roots by acidity stress.

Transcripts Analysis of DEGs Involved in Nutrients Transportation Using qRT-PCR
To confirm transcriptome results, qRT-PCR analysis was further conducted to investigate the expression of sets of DEGs involved in mineral nutrient acquisition and translocation, including twenty-four down-regulated and four up-regulated genes. Results showed a significant correlation presented between transcriptome and qRT-PCR analyses (r = 0.713, p < 0.01) (Figures 4 and S5). Among down-regulated genes, it was observed that expression of GmNRT2.4-1 and GmNRT2 was decreased by 19-and 8-folds, respectively (Figure 4). For the four up-regulated genes (i.e., GmSULTR3.1-1/2 and GmSULTR3.5-1/2), fold changes of their expression ranged from 1.7 to 7.4 ( Figure 4). Genes 2019, 9, x FOR PEER REVIEW 13 of 21

Transcripts Analysis of DEGs Involved in Nutrients Transportation Using qRT-PCR
To confirm transcriptome results, qRT-PCR analysis was further conducted to investigate the expression of sets of DEGs involved in mineral nutrient acquisition and translocation, including twenty-four down-regulated and four up-regulated genes. Results showed a significant correlation presented between transcriptome and qRT-PCR analyses (r = 0.713, p < 0.01) (Figures 4 and S5). Among down-regulated genes, it was observed that expression of GmNRT2.4-1 and GmNRT2 was decreased by 19-and 8-folds, respectively (Figure 4). For the four up-regulated genes (i.e., GmSULTR3.1-1/2 and GmSULTR3.5-1/2), fold changes of their expression ranged from 1.7 to 7.4 ( Figure 4).

Discussion
Low pH is a major constraint on crop productivity in the world's acid soil regions [2]. It has been documented that low soil pH not only limits nutrient availability, but also increases metal toxicity underground [21,59]. This poor growth condition leads to yield losses of more than 50% in grain crops [10]. So far, most studies have been focusing on plant growth limitations caused by P deficiency, or metal toxicities (i.e., Al and Mn) in acid soils [2,60]. However, the adverse effects of acidity stress alone, which is equally considered as a limiting factor for plant growth in acid soils, are often less emphasized [33,61]. In the study, effects of low-pH on soybean growth and its mineral nutrient homeostasis were investigated. Subsequently, transcriptomic analyses in soybean roots were conducted to shed light on molecular mechanisms underlying adaptation of soybean to acidity stress, especially mineral nutrients acquisition and translocation.
It showed that after the acid treatment, soybean root growth was significantly inhibited, as indicated by reduced root dry weight and total root length in both soil and hydroponic cultures (Figure 1 and Figure S1). This is consistent with previous reports in Arabidopsis, field beans (Vicia faba), corn (Zea mays) and lotus, in which significant decreases of root growth were found to be the main symptoms of proton toxicity [12][13][14][15]62]. Since mineral nutrients are mainly taken up by roots, it is generally assumed that mineral nutrient acquisition and translocation must be severely inhibited in plants under low-pH conditions. However, the present results showed that most essential mineral nutrients in soybean leaves and roots were differentially accumulated and tightly dependent on soybean organs and mineral nutrient types (Figure 2 and Figure S2). For examples, the N, P, K, and Mo concentrations in soybean roots were decreased, while concentrations of Ca, S, Fe and Mn were increased under acid conditions in soil culture (Figure 2). In the hydroponic culture, similar results were also observed in soybean roots, except for Fe concentration, which was slightly increased but not significant ( Figure S2). Other study of Erica also showed that acidity stress significantly affected the Mg and Mn concentrations in leaves, S, B and Fe concentrations in stems, and Ca, P, Fe, Cu and Zn concentrations in roots [25]. Thus, these results strongly suggested that the capability of mineral nutrients acquisition and translocation in plants could be influenced by proton toxicity.
To understand the molecular mechanisms underlying genes regulation in soybean responses to acidity stress, the transcriptome assay of soybean roots was conducted. It showed that long-term acidity stress mainly affects the expression of genes function in metabolic pathways ( Figure S4). Moreover, a proportion of low-pH up-regulated genes are categorized as pathogen response genes, including genes involved in glycollins biosynthesis, such as C4H (Glyma.20G114200), CHR (Glyma.02G307300), G4DT (Glyma.10G295300) and genes involved in SAR, such as DIR1 (Glyma.11G120400), WRKY40 (Glyma.14G103100), RLP32 (Glyma.16G169500), etc. (Table S3). This is consistent with the recent report that these genes were up-regulated in soybean variety Harosoy 63 after a 9-day pH3.0 treatment [35]. Therefore, it suggests that overlays between acidity-stress responses and pathogen responses are conserved in soybean.
Moreover, a total of thirty-eight DEGs associated with nutrient acquisition and transportation were identified, indicating that acid-responsive mineral nutrient homeostasis in soybean might be affected by the regulation of DEGs encoding mineral nutrient transporters (Table 3). For example, ten DEGs were identified to be involved in N transportation, including two AMTs, four NPFs and four NRTs (Table 3). It has been well documented that AMT1 and AMT2 proteins play critical roles in ammonium uptake and root-to-shoot translocation, respectively [30,[63][64][65][66][67]. Moreover, NPF and NRT2 are two major transporter families for NO 3 − absorption and transportation in plant species [68][69][70][71].
Our results showed that except for GmNPF3.1 and GmNPF4.6, all the other eight DEGs encoding NO 3 − transporters were significantly down-regulated, which exhibited a positive correlation with the reduced N concentration in low-pH-treated soybean roots ( Figure 2 and Table 3). Similar results were also found for P and K concentrations and expression of their corresponding transporter genes. It showed that P and K concentrations were significantly decreased in acid-treated soybean roots ( Figure 2). Meanwhile, all five Pi transporter genes (i.e., GmPT1/3/4/7/13) and five K transporter genes (i.e., GmHKT1, GmKUP3/6, GmTPK5, GmSKOR) were also down-regulated by acid-treatment (Table 3). Although functions of most Pi transporter and K transporter genes remain largely unknown in soybean, five GmPTs were suggested to encode high-affinity Pi transporters through complementary experiments in yeasts [72]. Meanwhile, the homologues of HKT, KUP, TPK, and SKOR have been demonstrated to function in K + absorption or translocation in other plant species [73][74][75][76]. Therefore, it is suggested that down-regulated expression of DEGs involved in N, P, and K transportation might play important roles in mediating N, P, and K homeostasis in soybean under acid conditions. In addition, eight GmSULTR genes were identified in the transcriptome assay, which were classified into GmSULTR2 and GmSULTR3 sub-groups [77] (Table 3). Transcripts of all three GmSULTR2s were down-regulated by acid treatments (Figure 4 and Table 3). Phylogenetic analysis showed that three GmSULTR2 members shared high similarity with SHST1, a putative sulfate transporter in Stylosanthes hamata that also could transport Mo [78]. Therefore, it was suggested that suppressed expression of GmSULTR2 genes might be responsible for the decreased Mo concentration in the acid-treated soybean roots (Figure 2). In contrast, five GmSULTR3 genes were all up-regulated by acid-treatments ( Table 3). The phylogenetic analysis showed high similarity between GmSULTR3.1/5 and AtSULTR3;1/5 in Arabidopsis ( Figure S6). It has been suggested that AtSULTR3;1 was responsible for sulfate transportation into chloroplasts [79], while AtSULTR3;5 was involved in root-to-shoot sulfate transportation and sulfate translocation within developing seeds [80,81]. Taken together with the increased sulfate accumulation in soybean plants under acid conditions (Figures 2 and  S2), it was indicated that up-regulated GmSULTR3 genes might play roles in the translocation or compartmentalization of excessive sulfate in roots.
Several metal transporter genes, including GmZIP1/2, GmZIP11-1/2, GmNramp5 and GmMTP7-1/2, were found to exhibit differential expression patterns under acid conditions. The homologues of ZIP and Nramp5 in other plant species have been suggested to mediate Mn uptake and/or translocation [82][83][84][85][86]. Meanwhile, phylogenetic analysis showed that GmMTP7-1/2 shared high similarity with Mn-CDF MTP proteins in rice and Arabidopsis ( Figure S7), which function in Mn acquisition and translocation [87,88]. Combining the higher Mn concentration in acid-treated soybean roots (Figures 2 and S2), it therefore suggests that suppressed expression of DEGs involving in Mn transportation might play roles in reducing Mn content in soybean roots under acid conditions.
Besides the DEGs involving in nutrient transportation, ten DEGs were identified to be associated with plant pH stat pathways ( Table 2). One of these pH stat pathways is the biochemical pH stat pathway, which plays important role in the H + consumption, and consists of alcohol dehydrogenase, malic enzyme and pyruvate decarboxylase [58,89,90]. The other pH stat pathway is the GABA shunt pathway, which consists of GDH, GAD, and GABA-T and also contributes greatly to the H + homeostasis in plant cells [58,91,92]. Expression of genes involving in these two pH-stat pathways was significantly suppressed in Arabidopsis stop1 mutants, strongly suggesting that they play important roles in Arabidopsis low-pH-tolerance [35]. In the present study, GmGDH1 in the GABA shunt pathway, and several other genes (i.e., GmME1-1, GmPCD1, GmPCD2-1/2, GmADH1-1/2/3) in the other biochemical pH stat pathway were significantly up-regulated under acid conditions (Table 2). Taken together, our results suggested that genes involving in pH-stat pathways could contribute greatly to the low-pH tolerance in soybean. The functions of DEGs associated with pH-stat pathways merit further characterization.
Finally, the identification of 79 transcription factors responding to low-pH treatment further shed light on the complex regulatory networks in soybean root responses to acidity stress (Table S4). It was recently reported that a NAC family transcription factor GmNAC42-1 (Glyma.02G284300), together with its down-stream genes IFS2 (Glyma.13G173500) and G4DT (Glyma.10G295300) was significantly up-regulated by 9-day low-pH treatment, indicating a regulatory link between the acidity and pathogen responses in soybean [35]. However, our transcripteomic data showed that the expression of GmNAC42-1 and IFS2 was no affected, while G4DT was up-regulated by the long-term acidity stress (30day) (Tables S3 and S4). Therefore, it indicates that there are might be other regulatory pathways involved in the pathogen responses under the long-term acid condition. Moreover, a PHR1-LIKE1 (GmPHL1) was down-regulated under acid conditions (Table S4). Although the function of GmPHL1 remains largely unknown in soybean, homologue of PHL1 in Arabidopsis is suggested to be closely related to PHR1 and act redundantly to regulate plant responses to Pi starvation [93]. It thus suggests that GmPHL1 might be involved in Pi homeostasis in response to acidity stress in soybean roots. Moreover, expression levels of two C2H2 transcription factors, GmSTOP2s, were up-regulated by acidity stress in soybean (Table S4). Since it is well known that AtSTOP2 in Arabidopsis partially confers low pH-tolerance by recovering expression of genes regulated by AtSTOP1 [94], our results indicate that GmSTOP2 genes might play important roles in soybean acid tolerance.
Overall, the present study showed that long-term acidity stress significantly inhibited soybean root growth and remarkably affected soybean mineral nutrients accumulation. Consistently, our transcriptome and qRT-PCR analyses showed that a set of DEGs might play important roles in mediating mineral nutrient homeostasis in soybean responses to acidity stress. Moreover, identification of DEGs associated with pH stat pathway strongly suggested that maintaining pH stat is one adaptive strategy of soybean responses to acidity stress. Taken together, the results herein provided new insights into physiological and molecular mechanisms underlying soybean adaptation to acidity stress.  Figure S2. Concentration of mineral nutrients in leaves and roots of soybean subjected to two pH treatments in hydroponic culture. Data in figures are means of three replicates with standard error bars. Asterisks indicate significant difference between pH 5.8 and pH 4.2 treatments in the Student's t-test (*: p < 0.05; **: 0.001 < p < 0.05). Figure S3. The GO analysis of DEGs in soybean roots. Figure S4. Pathway functional enrichment of DEGs in soybean roots. Figure S5. Correlation of gene expression levels between transcriptome data and qRT-PCR analysis. Figure S6. Phylogeny analysis of sulfate transporter proteins. Figure S7. Phylogeny analysis of metal tolerance proteins. Table S1. Primers used for qRT-PCR. Table S2. Summary of gene expression. Table S3. General information about DEGs between control (pH 5.8) and acid-treated (pH 4.2) soybean roots. Table S4. Transcription factors respond to acidity stress in soybean roots. Table S5. Abbreviations of DEGs.

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