Next Article in Journal
Cartilage Regeneration in Humans with Adipose Tissue-Derived Stem Cells and Adipose Stromal Vascular Fraction Cells: Updated Status
Next Article in Special Issue
Ectopic Expression of a Thellungiella salsuginea Aquaporin Gene, TsPIP1;1, Increased the Salt Tolerance of Rice
Previous Article in Journal
Dynamic DNA Methylation in Plant Growth and Development
Previous Article in Special Issue
Genome-Wide Identification and Characterization of CIPK Family and Analysis Responses to Various Stresses in Apple (Malus domestica)
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Early Transcriptomic Response to Phosphate Deprivation in Soybean Leaves as Revealed by RNA-Sequencing

1
College of Life and Environmental Sciences, Hangzhou Normal University, Hangzhou 310036, China
2
College of Resources and Environmental Sciences, Nanjing Agricultural University, Nanjing 210095, China
*
Authors to whom correspondence should be addressed.
Int. J. Mol. Sci. 2018, 19(7), 2145; https://doi.org/10.3390/ijms19072145
Submission received: 29 June 2018 / Revised: 17 July 2018 / Accepted: 17 July 2018 / Published: 23 July 2018
(This article belongs to the Special Issue Ion Transporters and Abiotic Stress Tolerance in Plants)

Abstract

:
Low phosphate (Pi) availability is an important limiting factor affecting soybean production. However, the underlying molecular mechanisms responsible for low Pi stress response and tolerance remain largely unknown, especially for the early signaling events under low Pi stress. Here, a genome-wide transcriptomic analysis in soybean leaves treated with a short-term Pi-deprivation (24 h) was performed through high-throughput RNA sequencing (RNA-seq) technology. A total of 533 loci were found to be differentially expressed in response to Pi deprivation, including 36 mis-annotated loci and 32 novel loci. Among the differentially expressed genes (DEGs), 303 were induced and 230 were repressed by Pi deprivation. To validate the reliability of the RNA-seq data, 18 DEGs were randomly selected and analyzed by quantitative RT-PCR (reverse transcription polymerase chain reaction), which exhibited similar fold changes with RNA-seq. Enrichment analyses showed that 29 GO (Gene Ontology) terms and 8 KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways were significantly enriched in the up-regulated DEGs and 25 GO terms and 16 KEGG pathways were significantly enriched in the down-regulated DEGs. Some DEGs potentially involved in Pi sensing and signaling were up-regulated by short-term Pi deprivation, including five SPX-containing genes. Some DEGs possibly associated with water and nutrient uptake, hormonal and calcium signaling, protein phosphorylation and dephosphorylation and cell wall modification were affected at the early stage of Pi deprivation. The cis-elements of PHO (phosphatase) element, PHO-like element and P responsive element were present more frequently in promoter regions of up-regulated DEGs compared to that of randomly-selected genes in the soybean genome. Our transcriptomic data showed an intricate network containing transporters, transcription factors, kinases and phosphatases, hormone and calcium signaling components is involved in plant responses to early Pi deprivation.

1. Introduction

As a non-substitutable macronutrient, phosphorus (P) is essential for plant growth and development by being part of fundamental bio-molecules and participating in various cellular activities. However, plants frequently suffer from P deficiency, which is a limiting factor for crop productivity worldwide, because of the low availability of the available form of P (phosphate, Pi) in soils, especially in acid soils [1,2]. Substantial use of Pi fertilizers is the main solution for this problem but the reserves of Pi rock (the primary source of Pi fertilizers) are finite and the cost of Pi fertilizers is high. In order to improve Pi acquisition and use efficiency of crop plants, many studies have been conducted to understand the physiological, biochemical and molecular mechanisms of plant adaptation and tolerance to low Pi stress [3,4,5,6].
Plants have a series of adaptive strategies to cope with low Pi stress, such as changes in shoot-to-root biomass ratio, exudation of organic acids, secretion of acid phosphatases, induction of high-affinity Pi transporters (PHTs), association with soil microorganisms, release of Pi from vacuolar stores and the optimization of internal Pi utilization by metabolic changes [7,8,9]. A large number of Pi-responsive genes have been functionally identified to coordinate these Pi deficiency responses. For example, in Arabidopsis, the MYB-CC (MYB-coiled–coil) transcription factor PHR1 and its close homologs PHR1-LIKE 1 (PHL1) and PHL2 are considered to be the central signaling components controlling transcriptional and metabolic responses to variations in Pi supply [10,11,12]; SPX-MFS (SPX-Major Facilitator Superfamily) protein VPT1 (Vacuolar Phosphate Transporter 1, also named PHT5;1) is regarded to mediate vacuolar Pi sequestration, which is critical for plant acclimation to varying Pi availability in the environment [13,14,15]; microRNAs are important players in controlling Pi transport by regulating the expression levels of some key components in Pi signaling networks, such as PHO2 and NLA (Nitrogen Limitation Adaptation) [16,17,18].
Transcriptomic analyses by microarray and recently developed high-throughput sequencing methods provide a useful tool for holistic understanding of the transcriptional responses under Pi deficiency in the post-genome era. Significant progress has been made on model plant Arabidopsis [19,20,21,22], rice [23,24] and some other plants [25,26,27]. Numerous genes were shown to be responsive to Pi deprivation, including genes involved in Pi acquisition and recycling, hormonal signaling and transcriptional regulation. These findings also support previous physiological and biochemical studies [28,29]. However, most of these studies are based on mid-term or long-term Pi deprivation. Early transcriptional responses to Pi deprivation are still largely ambiguous.
Soybean is an important leguminous crop for providing cooking oil and protein-rich food for humans. Low Pi availability is considered to be the most limiting nutrient condition for its production [30]. Some studies have been conducted to demonstrate the mechanisms of Pi deficiency response and tolerance. For instance, GmACP1 and GmACP2 both encoding an acid phosphatase was found to contribute to low Pi tolerance by combining quantitative trait loci (QTL) linkage, comparative transcriptome and genome-wide association analyses [31,32]; GmPHT1;4 (also named GmPT5) is critical for maintaining Pi supply in nodules by transporting Pi from roots to nodules [33]; three metabolism pathways including methane metabolism, phenylalanine metabolism and phenylpropanoid biosynthesis were found to be enriched in the common Pi-responsive genes between a low-Pi tolerant soybean accession and a low-Pi sensitive soybean accession [34]. Recently, we studied the transcriptomic responses to medium-term Pi deficiency (7 days) by digital gene expression deep sequencing and identified a total of 1612 differentially expressed genes (DEGs) in soybean roots [35]. However, the early signaling events in low Pi stress responses remain elusive. Leaves are important sources of the shoot-to-root systemic signals in Pi deficiency response, such as phloem mobile RNAs, proteins, sugars, Pi and other metabolites [4,9] but the sensing and signaling of Pi deficiency in leaves is largely unclear. In order to capture early transcriptional response to Pi deprivation in soybean leaves, we measured the transcriptional response to short-term (24 h) Pi deprivation by RNA-Seq. The DEGs and signaling components identified here represent new candidates for understanding the molecular mechanisms of early Pi stress responses in leaves and the improvement of Pi stress tolerance in soybean.

2. Results

2.1. Transcriptome Profiling by RNA-Seq of Soybean Leaves in Response to Short-Term Pi Deprivation

To assess the global transcriptome profile of soybean leaves in response to short-term Pi deprivation, we performed RNA-seq analyses following Pi deprivation for 24 h, a time point where the Pi concentration of Pi-deprived tissues decreases [23,36]. After Pi deprivation treatment for 24 h, Pi concentration in the roots, the first and the second trifoliate leaves of soybean significantly decreased but there was no significant difference for the third trifoliate leaves (Figure 1). Two biological replicates of RNA-seq were included for both Pi-sufficient leaves (PSL) and Pi-deprived leaves (PDL), and, therefore, a total of four libraries were constructed. By Illumina’s deep sequencing, a total of 38.8 to 49.6 million reliable clean reads were obtained from each library after excluding the low-quality reads (Table S1). Most of the clean reads (69.5–74.1%) from each library uniquely mapped to the soybean reference genome (v2.0). The correlation of the two biological replicates for PDL and PSL was calculated to determine the variability between the replicates. The Pearson’s correlation (R value) of the two comparisons was both around 85% (Figure S1), indicating a high correlation between the biological replicates. The abundance of transcripts from each gene mapped was measured in terms of the fragments per kilobase of transcript per million mapped reads (FPKM). A total of 39,505 gene loci were detected in PDL and/or PSL (Table S2).
Differential expression analysis showed that a total of 533 loci were differentially expressed, including 36 putative mis-annotated loci and 32 putative novel loci (Tables S3–S5). Here, mis-annotated loci were considered as loci assembled by Cufflinks but spanning two or more loci annotated in the soybean genome assembly. This is often due to the incorrect annotation of the initial gene model or annotated genes being too close together to be resolved. Unannotated loci (i.e., the novel loci) are an assembled group of reads not overlapping with any known loci annotated in the soybean genome assembly. Among the DEGs that have been annotated previously, 270 were up-regulated and 195 were down-regulated by short-term Pi deprivation (Figure S2). Interestingly, transcripts from nine DEGs were found only in PDL, while transcripts from 12 DEGs were found only in PSL (Table 1). The exclusive expression patterns suggest that these genes could have a potential role in plant response to short-term Pi deprivation.

2.2. Verification of RNA-Seq Data by qRT-PCR

To validate the reliability of the RNA-seq data, 18 DEGs were randomly selected and investigated by qRT-PCR. The fold changes of these genes under short-term Pi deficiency observed by qRT-PCR were similar to that revealed by RNA-seq (Figure 2), indicating that the RNA-seq data obtained in this study is reliable. However, there were differences between RNA-seq and qRT-PCR results for some genes. For example, the induction of Gm17g100000 transcripts appeared much stronger when detected by RNA-seq than by qRT-PCR, while the repression of Gm03145600 transcripts appeared much stronger when detected by qRT-PCR than by RNAseq. These discrepancies may be due to the facts that the samples used for RNA-seq and qRT-PCR were not the same and there would be inevitable differences between different batches of samples and perhaps some primer pairs used in qRT-PCR were not optimal for detecting the target transcripts.

2.3. Enrichment Analysis of Gene Ontology (GO) Functional Annotation and KEGG Pathways of Differentially Expressed Genes (DEGs)

Half of the DEGs (266/533) can be assigned to 183 GO terms (Table S3). GO enrichment analysis revealed that 29 and 25 GO terms were significantly enriched in the up-regulated DEGs and the down-regulated DEGs, respectively (Tables S6 and S7). For the up-regulated DEGs, carbohydrate metabolic process, cellular glucan metabolic process and lipid metabolic process were the three most significantly enriched GO terms within biological processes; each term contains at least three DEGs (Table S8). For the down-regulated DEGs, dTMP (deoxy-thymidine monophosphate) biosynthetic process, glycine biosynthetic process and nucleotide biosynthetic process were the three most significantly enriched GO terms within biological processes; each term contains at least two DEGs (Table S9). In addition, 15.6% of the DEGs (83/533) can be assigned to 152 KEGG pathways (Table S3). Eight and 16 KEGG pathways were significantly enriched in the up-regulated DEGs and the down-regulated DEGs, respectively (Table S10). Ether lipid metabolism, cutin, suberine and wax biosynthesis and pentose and glucuronate interconversions were the three most significantly enriched KEGG pathways in the up-regulated DEGs; each contains two to three DEGs. Drug metabolism-cytochrome P450, chemical carcinogenesis and metabolism of xenobiotics by cytochrome P450 were the three most significantly enriched KEGG pathways in the down-regulated DEGs; each contains three to four DEGs.

2.4. Genes Potentially Involved in Pi Signaling and Utilization Are Induced by Short-Term Pi Deprivation in Soybean Leaves

In this study, at least 13 classical Pi-responsive genes that are potentially involved in Pi signaling and utilization were found to have changed expression in response to short-term Pi deprivation (Table 2). These genes included five genes encoding SPX domain-containing proteins, two purple acid phosphatase (PAP) genes, three phospholipase genes, one glycerol-3-phosphate permease (G3Pp) gene, one sulfolipid sulfoquinovosyldiacylglycerol (SQD) gene and one gene encoding subfamily IIIB acid phosphatase (Gm16g220900). All of these genes were significantly up-regulated by Pi deprivation, with the exception that Gm16g220900—which potentially encodes a class IIIB acid phosphatase—was repressed.
SPX-containing proteins can be classified into four sub-families based on the presence of additional domains in their structure, namely the SPX, SPX-EXS, SPX-MFS and SPX-RING families. They are conserved in higher plants and are essential for Pi signaling and utilization [4,37,38]. Previously, nine SPX proteins and 14 SPX-EXS proteins have been identified in the soybean genome [39,40]. Here, an additional six SPX-MFS genes and four SPX-RING genes were identified in the soybean genome based on homology searches using Arabidopsis counterparts (AtPHT5;1 and AtNLA1) (Tables S11 and S12). Alignments with the sequences of SPX-MFS proteins and SPX-RING proteins were performed by ClustalW (Figures S3 and S4), suggesting that the SPX, MFS and RING domains are conserved among these proteins. Phylogenetic analysis indicated that soybean SPX-containing proteins are closely related to their homologs in Arabidopsis and rice (Figure 3), suggesting their conserved functions in land plants.

2.5. Expression of Genes Potentially Involved in Transportation of Water, Sugars and Mineral Nutrients Is Altered by Short-Term Pi Deprivation

Twenty-seven transporter genes were found to be differentially expressed in Pi-deprived soybean leaves (Table 3). Among them, several genes potentially involved in transportation of water; sugar; sulfate; and copper were up-regulated, while three genes potentially involved in zinc/iron transport were down-regulated by short-term Pi deprivation. In addition, one nitrate transporter gene (Gm18g127200) and one malate transporter gene (Gm11g179100) were down-regulated by Pi deprivation. The responsiveness of these transporter genes suggested that transport and allocation of water, sugars and nutrients are quickly altered by Pi deprivation in order to adapt to the stressed condition.

2.6. Genes Linked to Ca2+ and Hormonal Signaling are Regulated by Short-Term Pi Stress

At least 10 genes differentially expressed in Pi-deprived soybean leaves were putatively linked to Ca2+ signaling (Table 4). In addition, at least 15 DEGs are potentially involved in the transport and signaling of diverse hormones, including auxin, cytokinin, gibberellin (GA), brassinosteroids (BRs), jasmonate and ethylene (Table 4). Among these genes, DEGs related to GA, BR and jasmonate signaling were all up-regulated by Pi deficiency. The expressional changes of these hormone-related genes suggested that short-term Pi deficiency could affect hormone synthesis, transport and sensitivity.

2.7. Diverse Transcription Factor Family Genes Are Responsive to Short-Term Pi Deprivation in Soybean Leaves

Here, expression of at least 31 transcription factor genes was affected by short-term Pi deprivation in soybean leaves; among them, 14 were up-regulated and 17 were down-regulated (Table 5). Several transcription factors are also potentially involved in hormonal signaling. For example, Gm14g127400 encoding a BES1/BZR1 homolog protein is potentially involved in BR signaling and Gm05g144500 encoding an RR protein is potentially involved in cytokinin signaling (Table 4). These differentially expressed transcription factors belong to diverse families, such as MYB, WRKY, NAC (NAM, ATAF, and CUC), ERF (ethylene response factor), bHLH (basic helix-loop-helix), TCP, bZIP (basic leucine zipper domain), HD-ZIP (homeodomain leucine zipper), YABBY and zinc finger proteins (ZFPs) of various types (C2H2, C3H, B-box, DOF and HD). MYB/MYB-like, bHLH, C2H2 ZFP and YABBY were the most abundant transcription factor families differentially expressed; each of them contained at least three DEGs. Interestingly, four DEGs encoding bHLH transcription factors were all up-regulated, whereas two genes encoding TCP transcription factors were all repressed by Pi deprivation.

2.8. Short-Term Pi Deprivation Modifies the Expression of Genes Encoding Diverse Protein Kinases and Phosphatases

In this study, at least 31 protein kinase genes and nine phosphatase genes (including three acid phosphatase genes listed in Table 2) were found to be responsive to short-term Pi deprivation (Figure 4). These DEGs potentially encode protein kinases of diverse families. In addition, most of the DEGs encoding potential phosphatases were up-regulated by Pi deprivation. However, a gene encoding phosphoinositide phosphatase (Gm10g060600) and a PP2C gene (Gm11g050900) were repressed by Pi deprivation, suggesting their differential roles in Pi-deprivation response.

2.9. The Expression of Genes Associated with Metabolism Is Affected by Short-Term Pi Deprivation in Soybean Leaves

At least 77 DEGs were found to be associated with metabolism and the majority of them (72.7%) were up-regulated by short-term Pi deprivation (Figure 5, Table S13). Of the 15 DEGs potentially involved in lipid metabolism (including glycolipid synthesis, sulfolipid synthesis, fatty acid synthesis and elongation and lipid degradation), most were up-regulated under Pi deprivation. Ten DEGs potentially involved in cell wall degradation (including 1,4-β-mannan endohydrolases, pectate lyases and polygalacturonases) were all up-regulated by Pi deprivation. Seven of the 10 DEGs that are potentially involved in secondary metabolism of metabolites like carotenoids, terpenoids, isoflavonols and phenylpropanoids, were up-regulated by Pi deprivation. Moreover, nine DEGs were found to be potentially involved in the synthesis or degradation of some amino acids, like serine, cysteine, arginine, lysine, phenylalanine, tryptophan and proline (Table S13). Four of them were up-regulated and five were down-regulated by Pi deprivation, suggesting a complex regulation of amino acid metabolism under Pi deficiency and the different roles of these amino acids in Pi deficiency acclimations. In addition, nearly all of the 10 DEGs potentially involved in cell wall modification, such as expansins, xyloglucan endotransglycosylases, endoxyloglucan transferases and pectinesterases were up-regulated by short-term Pi deprivation. However, two DEGs potentially involved in cell wall cellulose synthesis (encoding cellulose synthases) were down-regulated by short-term Pi deprivation. These results suggest that plant metabolic acclimation to Pi deprivation stress occurs rapidly after encountering the stress.

2.10. Identification of Pi-Responsive Cis-Regulatory Elements in the Promoters of DEGs

We examined the distribution of Pi-responsive cis-regulatory elements in the putative promoter regions up to 1000 bp upstream of the transcription start sites of DEGs. In total 501 promoter regions were obtained from 288 up-regulated loci and 213 down-regulated loci (Tables S14 and S15). Eight types of previously identified Pi-responsive cis-elements were found to exist in Pi-responsive genes [11,41,42,43]. The cis-elements of PHO (phosphatase) element, PHO-like element and P responsive element were present more frequently in the promoter regions of up-regulated DEGs as compared to that of randomly-selected genes in the soybean genome (Table S15, Figure 6). In addition, the Helix-loop-helix element was present less frequently in the promoters of down-regulated DEGs. However, no significant difference was found for P1BS, TATA-box like, TC element and NIT 2-like. These results suggest that the Pi-responsive cis-elements enriched in the promoters of DEGs may be involved in the transcriptional regulation events at the early stage of Pi deprivation stress.

3. Discussion

Because of the low availability of Pi in soils, understanding the molecular acclimation mechanisms under Pi stress is of critical importance for developing crops with enhanced P-use efficiency in modern agriculture. Several studies utilizing microarrays or deep sequencing have documented the genome-wide transcriptional responses of plants to Pi stress but most of these studies are based on mid-term or long-term Pi deficiency and the early responses remain poorly understood [9]. In this study, we measured the transcriptional responses of soybean leaves to short-term Pi deprivation using an RNA-Seq approach. A total of 533 loci were found to be responsive to early Pi deprivation. By comparison with the recent transcriptomic analyses in two recombinant inbred lines of soybean that have different Pi stress tolerance [32], 12 genes were found to be commonly responsive to Pi deprivation in leaves; five of them have similar responses to short-term or long-term Pi stress (Figure S5). In addition, a small portion of them (28/533) were found to be responsive to long-term Pi deprivation in soybean roots [35] (Figure S6). The majority of the 28 common genes showed a similar transcriptional response under short-term and long-term Pi deprivation in the different organs. There were 21 DEGs expressed exclusively in either PDL or PSL (Table 1) but none of them were found to be responsive to long-term Pi deprivation in soybean roots. This may suggest their different responses to short-term and long-term Pi deficiency in leaves and roots.
At least 13 genes which are commonly responsive to Pi deprivation in plants were affected by short-term Pi deprivation, including five SPX domain-containing genes, two PAP genes, three phospholipase genes, one G3Pp gene and one SQD gene (Table 2). Some of them, such as PLDZ2, PAP13, PAP31, SPX3, SPX8 and Gm09g223700 were also found to be responsive to long-term Pi deficiency [35,39,44,45]. Whether the other genes are responsive to long-term Pi deprivation or only responsive to early Pi deprivation needs further investigation. Replacement of phospholipids in internal cellular membranes with galactolipids and sulfolipids is an important adaptive mechanism under Pi deficiency, which have been found in many plant species [46,47]. Lipid remodeling is also considered to be related to mobilization of P from less-essential cellular components to be exported to growing sinks [8]. PLDZ2 plays a critical role in galactolipid biosynthesis and thus facilitates Pi recycling from phospholipids [48]. SQD1 and SQD2 are Pi-responsive enzymes essential for sulfolipid biosynthesis under Pi-deficient conditions [49,50]; they also facilitate Pi recycling. In addition, PLDζ1 and PLDζ2 are Pi-responsive phospholipases that can hydrolyze phospholipids and thus contribute to Pi supply for Pi-deficient Arabidopsis [51]. Thus, the quick responses of these genes to Pi deprivation in leaves suggest that the replacement of phospholipids could be activated in a short time in order to adapt to Pi deficiency. G3Pp family genes which are potentially involved in transporting G3P, the hydrolysis products of diacylglycerol (a product of phospholipid breakdown), were found to be increased by Pi starvation in Arabidopsis [52]. Here, the homologs of these genes were also found to be induced by short-term Pi deprivation in soybean leaves but further researches are required to characterize their exact roles in Pi stress response and tolerance.
In this study, a total of 33 genes encoding SPX-containing proteins were found in soybean genome (Figure 3). The number is much higher than that in Arabidopsis (20) and rice (15) [38], which could be related to the whole-genome duplication events which occurred about 59 and 13 million years ago [53]. The early Pi deprivation-responsive GmSPX3, GmSPX4 and GmSPX8 are closely related to SPX1 and SPX2 in Arabidopsis and rice. In Arabidopsis and rice, SPX1, SPX2, SPX4 and SPX6 have been demonstrated to be involved in Pi sensing and signaling by inhibiting the transcriptional activity of PHR1 [54,55,56,57]. It has also been shown that overexpression of soybean SPX1 and SPX3 increases and decreases total P concentration in plant tissues, respectively, suggesting their contrasting roles in the regulation of Pi distribution in the plant [39,45]. Whether the short-term Pi deprivation-responsive SPX genes in soybean would function in Pi deficiency response by repressing the activity of conserved PHR1 and/or by participating in other regulatory pathways deserve further investigation. In addition to SPX genes, two SPX-EXS genes GmPHO1.H12 and GmPHO1.H14 were also induced by short-term Pi deprivation. Their closely related homologs in Arabidopsis, such as AtPHO1 and AtPHO1.H1 are also responsive to Pi deficiency and are associated with Pi economy [37,58]. However, the functions of Pi-responsive SPX-containing proteins in soybean remain to be investigated in the future.
In addition to the genes that are potentially involved in Pi signaling and utilization, many genes possibly involved in transportation of water, sugars and other mineral nutrients were also altered by short-term Pi deprivation (Table 3). It has been known for a long time that Pi deficiency decreases leaf water potential and transpiration rate [59]. Thus, it would be interesting to determine whether the observed up-regulation of aquaporins could be a compensatory response for the lower hydraulic conductance under Pi deficiency. Enhancement of the uptake and translocation of sugars and sulfate could also promote the synthesis of galactolipids and sulfolipids to substitute for the decline of phospholipids under Pi deficiency [8]. Sugar is an important systemic signal for regulating Pi starvation responses and root system architecture [60]. In addition, the homeostasis of iron, zinc and copper could also be affected by Pi deficiency. The most abundant Cu proteins in green tissues are plastocyanin and Cu/Zn-superoxide dismutase (Cu/ZnSOD), which are associated with electron transfers in photosynthesis and the scavenging of stress-induced reactive oxygen species (ROS), respectively [61]. Therefore, the up-regulation of copper transporters may be a part of a mechanism to boost ROS scavenging and photosynthesis which are impaired by Pi deficiency. Both iron and zinc have been shown to accumulate in Pi-deficient leaves in Arabidopsis [21]. Thus, the down-regulated iron/zinc transporter genes shortly after Pi deficiency could impede the excessive accumulation of these metals in the long run.
In this study, at least 10 DEGs are potentially involved in Ca2+ signaling (Table 4). Ca2+ is a universal signal playing a critical role in plant responses to environmental stresses [62,63]. It is well-known that diverse external environmental stimuli can quickly trigger specific spatial-temporal patterns of changes in cytosolic Ca2+ concentration, which can be perceived and decoded by a series of Ca2+ sensors containing EF-hand motifs [64,65]. However, little information is available with respect to the effect of Pi deficiency on cytosolic Ca2+ levels. Considering the critical roles of the Ca2+ signal in plant responses to other nutrient stress, such as potassium and nitrate deficiencies [66,67], it is conceivable that Ca2+ signals might be an important player in Pi stress response. In addition, Ca2+ and Pi are incompatible ions, because Ca2+ can form insoluble compounds with phosphate derivatives at high levels. Changing the allocation patterns of Pi may also necessitate changes to the patterns for Ca2+. Thus, how Ca2+ signal is linked to Pi deficiency response should deserve further investigations.
Many hormones have been involved in Pi stress responses by regulating root development and architecture as well as shoot development [4]. Here 15 genes potentially involved in signaling of auxin, CK (cytokinin), GA, BRs, jasmonate and ethylene were responsive to short-term Pi deprivation (Table 4). It has been shown that bioactive GA levels are reduced by Pi deficiency, which leads to the accumulation of DELLA proteins and, therefore, modulates root system architecture and anthocyanin accumulation in leaves [68]. BRs are a class of plant polyhydroxysteroids playing a pivotal role in plant growth and development as well as a wide variety of environmental stress responses [69]. However, there is little information available on the role of BRs in Pi deficiency. Whether BRs are involved in regulating the Pi economy during leaf development under Pi stress remains to be examined. Initiation and expansion of soybean leaves were demonstrated to decline under Pi deficiency [70]. However, whether hormonal signals like auxin, BRs, CK are involved in these physiological processes remains to be answered. In addition, leaf senescence is usually accelerated under nutrient stress in order to enhance remobilization of nutrients from senescing leaves [8] but the exact roles of hormones like ethylene, CK, jasmonate, auxin in nutrient stress-induced leaf senescence remain to be investigated in the future.
Transcription factors are critical components mediating gene regulatory networks under Pi stress. In the present study, thirty-one transcription factor genes belonging to 10 diverse families including MYB/MYB-like, bHLH, WRKY and ERF were found to be responsive to Pi deprivation (Table 5). Many of these transcription factor families were previously found to be responsive to Pi starvation and mediate transcriptional regulation of Pi-responsive genes in plants [21,71,72,73,74]. The identification of these short-term Pi stress-responsive transcription factors may provide preliminary evidence for further characterization of their functions in early Pi stress signaling.
Many genes possibly involved in protein phosphorylation and dephosphorylation are transcriptionally affected by short-term Pi stress (Figure 4). This result conforms with previous reports in Arabidopsis, demonstrating that many kinds of kinase and phosphatase genes are differentially expressed upon Pi deficiency [21,22,75]. Protein phosphorylation and dephosphorylation are linked with phosphate transporter trafficking and metabolic acclimations under Pi deficiency stress [76,77]. Recently, rice kinases CK2 and PSTOL1 were illustrated to be involved in regulating Pi stress response and tolerance [78,79]. Arabidopsis plasma membrane-localized receptor-like kinase BIK1 (botrytis-induced kinase 1) and MKK9-MPK3/MPK6 (mitogen-activated protein kinase kinase9-mitogen-activated protein kinase3/mitogen-activated protein kinase6) cascade were also shown to have functions in Pi signaling [80,81]. Moreover, some kinases like CIPKs can regulate the activities of transporters of nutrients like nitrate and potassium [82,83]. It would be interesting to examine the roles of Pi-responsive kinases or phosphatases in the early Pi stress signaling.
The expressions of many genes affected by short-term Pi deprivation are associated with metabolisms, such as lipid metabolism, cell wall degradation and modification, carbohydrate metabolism and amino acid metabolism (Figure 5). Earlier studies also demonstrated these primary or secondary metabolic changes upon medium-term or long-term Pi deficiency in plants [21,22,27]. Thus the modifications in expression levels of genes involved in metabolism can exist from the early stage of Pi stress. Cell walls have crucial functions in regulating the rate and direction of growth and determining the morphology of plant cells and organs. Synthesis and remodeling of the cell wall were documented to be associated with the acclimations to many kinds of environmental stresses [84].
In conclusion, our RNA-seq analysis revealed an early transcriptomic response of soybean leaves to Pi deprivation, suggesting an intricate regulatory network of signaling components upon short-term Pi stress. Quick changes in the transcript levels of various genes allow the plant to properly and accurately acclimate to Pi limitation conditions. Although the exact roles of these early Pi stress-responsive genes remain to be investigated, our data provide a platform for further functional characterizations of these genes in Pi stress sensing, signaling and tolerance.

4. Materials and Methods

4.1. Plant Material and Growth Conditions

Williams 82 is the soybean cultivar used for the reference soybean genome [53]. Soybean seeds (Glycine max var. Williams 82) (kindly provided by Prof. Haijian Zhi from Nanjing Agricultural University) were soaked in sterilized water for 4 h and then incubated at room temperature in the dark between two layers of moistened filter paper. Four days later, seedlings were transferred and grown hydroponically in a 10 L tank filled with a half-strength modified Hoagland nutrient solution containing 2.5 mm Ca(NO3)2, 2.5 mm KNO3, 0.5 mm KH2PO4, 1.25 mm MgSO4, 10.0 μm Fe-EDTA, 3.4 μm MnSO4, 0.16 μm CuSO4, 0.38 μm ZnSO4, 23.0 μm H3BO3, 0.25 μm Na2MoO4, with pH adjusted to 5.6. Nutrient solution was changed every two days. Plants were grown in a growth chamber with a photoperiod set at 16-h-light/8-h-dark at 26/22 °C and light intensity set at 150 μL·m−2·s−1. After 18 days of cultivation in hydroponics, soybean seedlings with the first trifoliate true leaves fully expanded were transferred into Pi-sufficient (0.5 mm KH2PO4) or Pi-deprived (0 mm KH2PO4, K2SO4 was substituted for KH2PO4 to keep the concentration of K+) nutrient solutions. Roots and leaves of soybean seedlings were separately sampled after 24 h Pi deprivation, frozen in liquid nitrogen and stored at −80 °C until RNA preparation.

4.2. Pi Concentration Determination

Pi concentrations were analyzed as described previously [35]. About 0.2 g fresh tissue that were frozen in liquid nitrogen were ground into fine powder and suspended in extraction buffer (10 mm Tris, 1 mm EDTA, 100 mm NaCl and 1 mm β-mercaptoethanol, pH 8.0) at a ratio of 1 mg of fresh weight sample to 10 μL of extraction buffer. A total of 100 μL of sample suspension was mixed with 900 μL of 1% glacial acetic acid and incubated at 42 °C for 30 min. Then, the suspension was centrifuged at 13,000× g for 10 min and 500 μL of the supernatant was used for the Pi quantitation assay. The reaction mixture containing 1000 μL of Pi assay solution (0.34% (NH4)6Mo7O24·4H2O, 0.46 M H2SO4 and 1.4% ascorbic acid) and 500 μL of supernatant was incubated at 42 °C for 30 min, cooled on ice and the absorbance at 820 nm was measured using a UV-Vis spectrum meter (Thermo Scientific BioMate 3S, Chino, CA, USA).

4.3. RNA Isolation, Library Construction and RNA Sequencing

Soybean leaves (the first trifoliate true leaves) were collected after Pi-deprivation treatment for 24 h. Four samples (two biological replicates of both Pi-deprived and Pi-sufficient leaves) were used for mRNA library construction and sequencing. Each biological replicate was sampled from three different randomly-selected plants. Total RNA was extracted using Trizol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s procedure. The total RNA quantity and purity were determined by Agilent Bioanalyzer 2100 with RNA 6000 Nano LabChip Kit (Agilent, Santa Clara, CA, USA). The RIN (RNA integrity number) of all RNA samples were determined to be more than 7.0. Approximately 10 µg of total RNA was subjected to Poly (A) mRNA isolation using poly-T oligo-attached magnetic beads (Invitrogen). The mRNA was fragmented into small pieces using divalent cations under elevated temperature after purification. The cleaved RNA fragments were reverse-transcribed to create the final cDNA library in accordance with the protocol for the mRNA-Seq sample preparation kit (Illumina, San Diego, CA, USA). The average insert size for the paired-end libraries was 300 ± 50 bp. Subsequently, paired-end sequencing was performed on an Illumina Hiseq 2000 at the LC-BIO TECHNOLOGIES (Hangzhou, China) following the instructions from Illumina.

4.4. RNA-Seq Reads Mapping and Differential Counting

The initial base calling and quality filtering of the reads generated with the Illumina analysis pipeline (Fastq format) were implemented using a custom Perl script and the default parameters of the Illumina pipeline (http://www.illumina.com). Additional filtering for poor-quality bases was carried out using the FASTX-toolkit available in the FastQC software package (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). To facilitate the read mapping, the Glycine max reference genome (Gmax2.0 version) was indexed by Bowtie2 (http://www.phytozome.net) [85]. The read mapping was conducted using the TopHat software package [86]. TopHat allows multiple alignments per read (up to 40) and a maximum of two mismatches when mapping the reads to the reference genome. The reads were first mapped directly to the genome using indexing and then the unmapped reads were used to identify novel splicing events. The aligned read files were processed by Cufflinks to measure the relative abundances of the transcripts by using the normalized RNA-seq fragment counts [87]. The estimated abundance of genes was measured in terms of the fragments per kilobase of transcript per million mapped reads (FPKM). Differentially expressed genes (DEGs) between the two sets of samples were identified using Cuffdiff [87]. Only the genes with a log2 fold change ≥1 or ≤−1 and a p-value ≤ 0.05 were considered as significantly DEGs.
The datasets were deposited in NCBI’s Gene Expression Omnibus and are accessible through GEO accession number GSE104286 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE104286) (the secure token for review is mbuhcgoylbanjyx).

4.5. Quantitative RT-PCR Analysis

Total RNA was extracted from soybean tissues using RNApure Plant Kit (with DNase I) (CoWin Biotech, Beijing, China) and digested with DNase I to eliminate genomic DNA contamination according to the manufacturer’s instruction. cDNA was synthesized from 1.0 μg total RNA in a 20 μL reaction by SuperRT Reverse Transcriptase (CoWin Biotech) using oligo (dT) primers. Quantitative RT-PCR (qRT-PCR) was performed on a MyiQ Single Color Real-time PCR system (Bio-Rad, Hercules, CA, USA) as described previously [88]. Briefly, two μL of a 1/10 dilution of cDNA in water was added to 10 μL of 2 × UltraSYBR (with Rox) (CoWin Biotech); together with two gene-specific primers (200 nm each) the final volume was brought to 20 uL by adding DNase-free water. The procedures for PCR were as follows:15 °C for 10 min; 40 cycles of 95 °C for 15 s, 60 °C for 60 s. Amplifications were run in triplicate together with controls that contained no template and no reverse transcription for each of the examined genes. Relative expression levels were normalized to that of an internal control ACTIN11 (Glyma.18g290800) using the Pfaffl method (Ratio = (Etarget)Cttarget(control−sample)/(Eref)Ctref(control−sample)) [89]. The calculated efficiency (E) of all primers was between 1.7 and 2.2. The gene-specific primers are listed in Table S16.

4.6. Functional Annotation and Gene Ontology (GO) Enrichment

The DEGs were annotated for gene ontology (GO) terms [90] and categorized into molecular function, cellular component and biological process categories. A gene enrichment test was performed on each of the gene lists to acquire the terms that were significantly enriched among the DEGs. Fisher’s exact test, which is based on hyper-geometric distribution, was used to calculate the p-value. A GO category (http://geneontology.org/) or KEGG pathway (http://www.genome.jp/kegg/) with a p-value ≤ 0.05 was regarded as significantly enriched. GO and KEGG enrichment analyses were conducted with the help of LC-BIO company (Hangzhou, China).
The MultiExperiment Viewer program (http://mev.tm4.org/) was used to draw a heatmap of the significant DEGs in response to short-term Pi deprivation. The metabolic pathways were plotted using MapMan (http://mapman.gabipd.org/) [91].

4.7. Promoter Cis-Element Analysis

Promoter sequences (1000 bp upstream) of the transcription start sites of DEGs were extracted from the SoyBase database (http://www.soybase.org/). The presence of the Pi-responsive cis-elements was analyzed using Regulatory Sequence Analysis Tools (http://floresta.eead.csic.es/rsat/) [92]. The sequences of eight types of Pi-responsive cis-elements were as follows: P1BS (GNATATNC) [11]; PHO (CACGT(G/C)); PHO-like (C(G/T/A)(C/T/A)GTGG) [19]; P responsive (ATGCCAT) [43]; TATA box like (TATAAATA) [19]; TC element (TCTCTCT); NIT 2-like (AAATATCT); Helix-loop-helix (CA(T/G)(A/C)TG) [41].

Supplementary Materials

Supplementary materials can be found at https://www.mdpi.com/1422-0067/19/7/2145/s1.

Author Contributions

H.Z. and Y.Z. conceived the study and wrote the manuscript. H.Z., X.Z. (Xiajun Zhang), X.Z. (Xin Zhang) and L.X. performed the experiments. H.Z. analyzed the data. All authors read and approved the final manuscript.

Funding

This work was supported by grants from the National Key Basic Research and Development Program of China (2017YFD0200206), the National Natural Science Foundation of China (31201679), Zhejiang Provincial Natural Science Foundation of China (LY15C020006).

Acknowledgments

We thank Keke Yi (Chinese Academy of Agricultural Sciences) for critical reading of the manuscript. We thank Haijian Zhi (Nanjing Agricultural University) for kindly providing the soybean seeds.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Shen, J.; Yuan, L.; Zhang, J.; Li, H.; Bai, Z.; Chen, X.; Zhang, W.; Zhang, F. Phosphorus dynamics: From soil to plant. Plant Physiol. 2011, 156, 997–1005. [Google Scholar] [CrossRef] [PubMed]
  2. Rao, I.M.; Miles, J.W.; Beebe, S.E.; Horst, W.J. Root adaptations to soils with low fertility and aluminium toxicity. Ann. Bot. 2016, 118, 593–605. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Heuer, S.; Gaxiola, R.; Schilling, R.; Herrera-Estrella, L.; Lopez-Arredondo, D.; Wissuwa, M.; Delhaize, E.; Rouached, H. Improving phosphorus use efficiency: A complex trait with emerging opportunities. Plant J. 2017, 90, 868–885. [Google Scholar] [CrossRef] [PubMed]
  4. Chiou, T.J.; Lin, S.I. Signaling network in sensing phosphate availability in plants. Annu. Rev. Plant Biol. 2011, 62, 185–206. [Google Scholar] [CrossRef] [PubMed]
  5. Lopez-Arredondo, D.L.; Leyva-Gonzalez, M.A.; Gonzalez-Morales, S.I.; Lopez-Bucio, J.; Herrera-Estrella, L. Phosphate nutrition: Improving low-phosphate tolerance in crops. Annu. Rev. Plant Biol. 2014, 65, 95–123. [Google Scholar] [CrossRef] [PubMed]
  6. Wu, P.; Shou, H.; Xu, G.; Lian, X. Improvement of phosphorus efficiency in rice on the basis of understanding phosphate signaling and homeostasis. Curr. Opin. Plant Biol. 2013, 16, 205–212. [Google Scholar] [CrossRef] [PubMed]
  7. Gu, M.; Chen, A.; Sun, S.; Xu, G. Complex Regulation of Plant Phosphate Transporters and the Gap between Molecular Mechanisms and Practical Application: What Are Missing? Mol. Plant 2016, 9, 396–416. [Google Scholar] [CrossRef] [PubMed]
  8. Lambers, H.; Finnegan, P.M.; Jost, R.; Plaxton, W.C.; Shane, M.W.; Stitt, M. Phosphorus nutrition in Proteaceae and beyond. Nat. Plants 2015, 1, 15109. [Google Scholar] [CrossRef] [PubMed]
  9. Zhang, Z.; Liao, H.; Lucas, W.J. Molecular mechanisms underlying phosphate sensing, signaling, and adaptation in plants. J. Integr. Plant Biol. 2014, 56, 192–220. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  10. Sun, L.; Song, L.; Zhang, Y.; Zheng, Z.; Liu, D. Arabidopsis PHL2 and PHR1 Act Redundantly as the Key Components of the Central Regulatory System Controlling Transcriptional Responses to Phosphate Starvation. Plant Physiol. 2016, 170, 499–514. [Google Scholar] [CrossRef] [PubMed]
  11. Rubio, V.; Linhares, F.; Solano, R.; Martin, A.C.; Iglesias, J.; Leyva, A.; Paz-Ares, J. A conserved MYB transcription factor involved in phosphate starvation signaling both in vascular plants and in unicellular algae. Genes Dev. 2001, 15, 2122–2133. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  12. Pant, B.D.; Pant, P.; Erban, A.; Huhman, D.; Kopka, J.; Scheible, W.R. Identification of primary and secondary metabolites with phosphorus status-dependent abundance in Arabidopsis, and of the transcription factor PHR1 as a major regulator of metabolic changes during phosphorus limitation. Plant Cell Environ. 2015, 38, 172–187. [Google Scholar] [CrossRef] [PubMed]
  13. Liu, J.; Yang, L.; Luan, M.; Wang, Y.; Zhang, C.; Zhang, B.; Shi, J.; Zhao, F.G.; Lan, W.; Luan, S. A vacuolar phosphate transporter essential for phosphate homeostasis in Arabidopsis. Proc. Natl. Acad. Sci. USA 2015, 112, E6571–E6578. [Google Scholar] [CrossRef] [PubMed]
  14. Liu, T.Y.; Huang, T.K.; Yang, S.Y.; Hong, Y.T.; Huang, S.M.; Wang, F.N.; Chiang, S.F.; Tsai, S.Y.; Lu, W.C.; Chiou, T.J. Identification of plant vacuolar transporters mediating phosphate storage. Nat. Commun. 2016, 7, 11095. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  15. Wang, C.; Yue, W.; Ying, Y.; Wang, S.; Secco, D.; Liu, Y.; Whelan, J.; Tyerman, S.D.; Shou, H. Rice SPX-Major Facility Superfamily3, a Vacuolar Phosphate Efflux Transporter, Is Involved in Maintaining Phosphate Homeostasis in Rice. Plant Physiol. 2015, 169, 2822–2831. [Google Scholar] [PubMed]
  16. Zeng, H.; Wang, G.; Hu, X.; Wang, H.; Du, L.; Zhu, Y. Role of microRNAs in plant responses to nutrient stress. Plant Soil 2014, 374, 1005–1021. [Google Scholar] [CrossRef]
  17. Wang, L.; Zeng, H.Q.; Song, J.; Feng, S.J.; Yang, Z.M. miRNA778 and SUVH6 are involved in phosphate homeostasis in Arabidopsis. Plant Sci. 2015, 238, 273–285. [Google Scholar] [CrossRef] [PubMed]
  18. Lin, W.Y.; Huang, T.K.; Chiou, T.J. Nitrogen limitation adaptation, a target of microRNA827, mediates degradation of plasma membrane-localized phosphate transporters to maintain phosphate homeostasis in Arabidopsis. Plant Cell 2013, 25, 4061–4074. [Google Scholar] [CrossRef] [PubMed]
  19. Hammond, J.P.; Bennett, M.J.; Bowen, H.C.; Broadley, M.R.; Eastwood, D.C.; May, S.T.; Rahn, C.; Swarup, R.; Woolaway, K.E.; White, P.J. Changes in gene expression in Arabidopsis shoots during phosphate starvation and the potential for developing smart plants. Plant Physiol. 2003, 132, 578–596. [Google Scholar] [CrossRef] [PubMed]
  20. Wu, P.; Ma, L.; Hou, X.; Wang, M.; Wu, Y.; Liu, F.; Deng, X.W. Phosphate starvation triggers distinct alterations of genome expression in Arabidopsis roots and leaves. Plant Physiol. 2003, 132, 1260–1271. [Google Scholar] [CrossRef] [PubMed]
  21. Misson, J.; Raghothama, K.G.; Jain, A.; Jouhet, J.; Block, M.A.; Bligny, R.; Ortet, P.; Creff, A.; Somerville, S.; Rolland, N. A genome-wide transcriptional analysis using Arabidopsis thaliana Affymetrix gene chips determined plant responses to phosphate deprivation. Proc. Natl. Acad. Sci. USA 2005, 102, 11934–11939. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  22. Morcuende, R.; Bari, R.; Gibon, Y.; Zheng, W.; Pant, B.D.; Blasing, O.; Usadel, B.; Czechowski, T.; Udvardi, M.K.; Stitt, M.; et al. Genome-wide reprogramming of metabolism and regulatory networks of Arabidopsis in response to phosphorus. Plant Cell Environ. 2007, 30, 85–112. [Google Scholar] [CrossRef] [PubMed]
  23. Secco, D.; Jabnoune, M.; Walker, H.; Shou, H.; Wu, P.; Poirier, Y.; Whelan, J. Spatio-temporal transcript profiling of rice roots and shoots in response to phosphate starvation and recovery. Plant Cell 2013, 25, 4285–4304. [Google Scholar] [CrossRef] [PubMed]
  24. Cai, H.; Xie, W.; Zhu, T.; Lian, X. Transcriptome response to phosphorus starvation in rice. Acta Physiol. Plant. 2012, 34, 327–341. [Google Scholar] [CrossRef]
  25. Oono, Y.; Kobayashi, F.; Kawahara, Y.; Yazawa, T.; Handa, H.; Itoh, T.; Matsumoto, T. Characterisation of the wheat (Triticum aestivum L.) transcriptome by de novo assembly for the discovery of phosphate starvation-responsive genes: Gene expression in Pi-stressed wheat. BMC Genom. 2013, 14, 77. [Google Scholar] [CrossRef] [PubMed]
  26. Sun, Y.; Mu, C.; Chen, Y.; Kong, X.; Xu, Y.; Zheng, H.; Zhang, H.; Wang, Q.; Xue, Y.; Li, Z.; et al. Comparative transcript profiling of maize inbreds in response to long-term phosphorus deficiency stress. Plant Physiol. Biochem. 2016, 109, 467–481. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  27. O’Rourke, J.A.; Yang, S.S.; Miller, S.S.; Bucciarelli, B.; Liu, J.; Rydeen, A.; Bozsoki, Z.; Uhde-Stone, C.; Tu, Z.J.; Allan, D.; et al. An RNA-Seq transcriptome analysis of orthophosphate-deficient white lupin reveals novel insights into phosphorus acclimation in plants. Plant Physiol. 2013, 161, 705–724. [Google Scholar] [CrossRef] [PubMed]
  28. Theodorou, M.E.; Plaxton, W.C. Metabolic Adaptations of Plant Respiration to Nutritional Phosphate Deprivation. Plant Physiol. 1993, 101, 339–344. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  29. Rao, I.M.; Terry, N. Leaf Phosphate Status, Photosynthesis, and Carbon Partitioning in Sugar Beet (IV. Changes with Time Following Increased Supply of Phosphate to Low-Phosphate Plants). Plant Physiol. 1995, 107, 1313–1321. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  30. Wang, X.; Yan, X.; Liao, H. Genetic improvement for phosphorus efficiency in soybean: A radical approach. Ann. Bot. 2010, 106, 215–222. [Google Scholar] [CrossRef] [PubMed]
  31. Zhang, D.; Song, H.; Cheng, H.; Hao, D.; Wang, H.; Kan, G.; Jin, H.; Yu, D. The Acid Phosphatase-Encoding Gene GmACP1 Contributes to Soybean Tolerance to Low-Phosphorus Stress. PLoS Genet. 2014, 10, e1004061. [Google Scholar] [CrossRef] [PubMed]
  32. Zhang, D.; Zhang, H.; Chu, S.; Li, H.; Chi, Y.; Triebwasser-Freese, D.; Lv, H.; Yu, D. Integrating QTL mapping and transcriptomics identifies candidate genes underlying QTLs associated with soybean tolerance to low-phosphorus stress. Plant Mol. Biol. 2017, 93, 137–150. [Google Scholar] [CrossRef] [PubMed]
  33. Qin, L.; Zhao, J.; Tian, J.; Chen, L.; Sun, Z.; Guo, Y.; Lu, X.; Gu, M.; Xu, G.; Liao, H. The high-affinity phosphate transporter GmPT5 regulates phosphate transport to nodules and nodulation in soybean. Plant Physiol. 2012, 159, 1634–1643. [Google Scholar] [CrossRef] [PubMed]
  34. Wang, Q.; Wang, J.; Yang, Y.; Du, W.; Zhang, D.; Yu, D.; Cheng, H. A genome-wide expression profile analysis reveals active genes and pathways coping with phosphate starvation in soybean. BMC Genom. 2016, 17, 1. [Google Scholar] [CrossRef] [PubMed]
  35. Zeng, H.Q.; Wang, G.P.; Zhang, Y.Q.; Hu, X.Y.; Pi, E.X.; Zhu, Y.Y.; Wang, H.Z.; Du, L.Q. Genome-wide identification of phosphate-deficiency-responsive genes in soybean roots by high-throughput sequencing. Plant Soil 2016, 398, 207–227. [Google Scholar] [CrossRef]
  36. Zhang, Z.; Zheng, Y.; Ham, B.K.; Chen, J.; Yoshida, A.; Kochian, L.V.; Fei, Z.; Lucas, W.J. Vascular-mediated signalling involved in early phosphate stress response in plants. Nat. Plants 2016, 2, 16033. [Google Scholar] [CrossRef] [PubMed]
  37. Jung, J.-Y.; Ried, M.K.; Hothorn, M.; Poirier, Y. Control of plant phosphate homeostasis by inositol pyrophosphates and the SPX domain. Curr. Opin. Biotechnol. 2018, 49, 156–162. [Google Scholar] [CrossRef] [PubMed]
  38. Secco, D.; Wang, C.; Arpat, B.A.; Wang, Z.; Poirier, Y.; Tyerman, S.D.; Wu, P.; Shou, H.; Whelan, J. The emerging importance of the SPX domain-containing proteins in phosphate homeostasis. New Phytol. 2012, 193, 842–851. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  39. Yao, Z.; Tian, J.; Liao, H. Comparative characterization of GmSPX members reveals that GmSPX3 is involved in phosphate homeostasis in soybean. Ann. Bot. 2014, 114, 477–488. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  40. He, L.; Zhao, M.; Wang, Y.; Gai, J.; He, C. Phylogeny, structural evolution and functional diversification of the plant PHOSPHATE1 gene family: A focus on Glycine max. BMC Evol. Biol. 2013, 13, 103. [Google Scholar] [CrossRef] [PubMed]
  41. Hammond, J.P.; Broadley, M.R.; White, P.J. Genetic responses to phosphorus deficiency. Ann. Bot. 2004, 94, 323–332. [Google Scholar] [CrossRef] [PubMed]
  42. Zeng, H.Q.; Zhu, Y.Y.; Huang, S.Q.; Yang, Z.M. Analysis of phosphorus-deficient responsive miRNAs and cis-elements from soybean (Glycine max L.). J. Plant Physiol. 2010, 167, 1289–1297. [Google Scholar] [CrossRef] [PubMed]
  43. Mukatira, U.T.; Liu, C.; Varadarajan, D.K.; Raghothama, K.G. Negative regulation of phosphate starvation-induced genes. Plant Physiol. 2001, 127, 1854–1862. [Google Scholar] [CrossRef] [PubMed]
  44. Li, C.; Gui, S.; Yang, T.; Walk, T.; Wang, X.; Liao, H. Identification of soybean purple acid phosphatase genes and their expression responses to phosphorus availability and symbiosis. Ann. Bot. 2012, 109, 275–285. [Google Scholar] [CrossRef] [PubMed]
  45. Zhang, J.; Zhou, X.; Xu, Y.; Yao, M.; Xie, F.; Gai, J.; Li, Y.; Yang, S. Soybean SPX1 is an important component of the response to phosphate deficiency for phosphorus homeostasis. Plant Sci. 2016, 248, 82–91. [Google Scholar] [CrossRef] [PubMed]
  46. Andersson, M.X.; Stridh, M.H.; Larsson, K.E.; Liljenberg, C.; Sandelius, A.S. Phosphate-deficient oat replaces a major portion of the plasma membrane phospholipids with the galactolipid digalactosyldiacylglycerol. FEBS Lett. 2003, 537, 128–132. [Google Scholar] [CrossRef] [Green Version]
  47. Lambers, H.; Cawthray, G.R.; Giavalisco, P.; Kuo, J.; Laliberté, E.; Pearse, S.J.; Scheible, W.-R.; Stitt, M.; Teste, F.; Turner, B.L. Proteaceae from severely phosphorus-impoverished soils extensively replace phospholipids with galactolipids and sulfolipids during leaf development to achieve a high photosynthetic phosphorus-use-efficiency. New Phytol. 2012, 196, 1098–1108. [Google Scholar] [CrossRef] [PubMed]
  48. Cruz-Ramirez, A.; Oropeza-Aburto, A.; Razo-Hernandez, F.; Ramirez-Chavez, E.; Herrera-Estrella, L. Phospholipase DZ2 plays an important role in extraplastidic galactolipid biosynthesis and phosphate recycling in Arabidopsis roots. Proc. Natl. Acad. Sci. USA 2006, 103, 6765–6770. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  49. Essigmann, B.; Guler, S.; Narang, R.A.; Linke, D.; Benning, C. Phosphate availability affects the thylakoid lipid composition and the expression of SQD1, a gene required for sulfolipid biosynthesis in Arabidopsis thaliana. Proc. Natl. Acad. Sci. USA 1998, 95, 1950–1955. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  50. Yu, B.; Xu, C.; Benning, C. Arabidopsis disrupted in SQD2 encoding sulfolipid synthase is impaired in phosphate-limited growth. Proc. Natl. Acad. Sci. USA 2002, 99, 5732–5737. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  51. Li, M.; Welti, R.; Wang, X. Quantitative profiling of Arabidopsis polar glycerolipids in response to phosphorus starvation. Roles of phospholipases D zeta1 and D zeta2 in phosphatidylcholine hydrolysis and digalactosyldiacylglycerol accumulation in phosphorus-starved plants. Plant Physiol. 2006, 142, 750–761. [Google Scholar] [CrossRef] [PubMed]
  52. Ramaiah, M.; Jain, A.; Baldwin, J.C.; Karthikeyan, A.S.; Raghothama, K.G. Characterization of the phosphate starvation-induced glycerol-3-phosphate permease gene family in Arabidopsis. Plant Physiol. 2011, 157, 279–291. [Google Scholar] [CrossRef] [PubMed]
  53. Schmutz, J.; Cannon, S.B.; Schlueter, J.; Ma, J.; Mitros, T.; Nelson, W.; Hyten, D.L.; Song, Q.; Thelen, J.J.; Cheng, J.; et al. Genome sequence of the palaeopolyploid soybean. Nature 2010, 463, 178–183. [Google Scholar] [CrossRef] [PubMed]
  54. Puga, M.I.; Mateos, I.; Charukesi, R.; Wang, Z.; Franco-Zorrilla, J.M.; de Lorenzo, L.; Irigoyen, M.L.; Masiero, S.; Bustos, R.; Rodriguez, J.; et al. SPX1 is a phosphate-dependent inhibitor of Phosphate Starvation Response 1 in Arabidopsis. Proc. Natl. Acad. Sci. USA 2014, 111, 14947–14952. [Google Scholar] [CrossRef] [PubMed]
  55. Lv, Q.; Zhong, Y.; Wang, Y.; Wang, Z.; Zhang, L.; Shi, J.; Wu, Z.; Liu, Y.; Mao, C.; Yi, K.; et al. SPX4 Negatively Regulates Phosphate Signaling and Homeostasis through Its Interaction with PHR2 in Rice. Plant Cell 2014, 26, 1586–1597. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  56. Wang, Z.Y.; Ruan, W.Y.; Shi, J.; Zhang, L.; Xiang, D.; Yang, C.; Li, C.Y.; Wu, Z.C.; Liu, Y.; Yu, Y.A.; et al. Rice SPX1 and SPX2 inhibit phosphate starvation responses through interacting with PHR2 in a phosphate-dependent manner. Proc. Natl. Acad. Sci. USA 2014, 111, 14953–14958. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  57. Zhong, Y.; Wang, Y.; Guo, J.; Zhu, X.; Shi, J.; He, Q.; Liu, Y.; Wu, Y.; Zhang, L.; Lv, Q. Rice SPX6 negatively regulates the phosphate starvation response through suppression of the transcription factor PHR2. New Phytol. 2018, 219, 135–148. [Google Scholar] [CrossRef] [PubMed]
  58. Stefanovic, A.; Ribot, C.; Rouached, H.; Wang, Y.; Chong, J.; Belbahri, L.; Delessert, S.; Poirier, Y. Members of the PHO1 gene family show limited functional redundancy in phosphate transfer to the shoot, and are regulated by phosphate deficiency via distinct pathways. Plant J. 2007, 50, 982–994. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  59. Radin, J.W. Stomatal responses to water stress and to abscisic acid in phosphorus-deficient cotton plants. Plant Physiol. 1984, 76, 392–394. [Google Scholar] [CrossRef] [PubMed]
  60. Hammond, J.P.; White, P.J. Sugar signaling in root responses to low phosphorus availability. Plant Physiol. 2011, 156, 1033–1040. [Google Scholar] [CrossRef] [PubMed]
  61. Yruela, I. Copper in plants: Acquisition, transport and interactions. Funct. Plant Biol. 2009, 36, 409–430. [Google Scholar] [CrossRef]
  62. Yuan, F.; Yang, H.; Xue, Y.; Kong, D.; Ye, R.; Li, C.; Zhang, J.; Theprungsirikul, L.; Shrift, T.; Krichilsky, B.; et al. OSCA1 mediates osmotic-stress-evoked Ca2+ increases vital for osmosensing in Arabidopsis. Nature 2014, 514, 367–371. [Google Scholar] [CrossRef] [PubMed]
  63. Zeng, H.; Xu, L.; Singh, A.; Wang, H.; Du, L.; Poovaiah, B.W. Involvement of calmodulin and calmodulin-like proteins in plant responses to abiotic stresses. Front. Plant Sci. 2015, 6, 600. [Google Scholar] [CrossRef] [PubMed]
  64. Zeng, H.; Zhang, Y.; Zhang, X.; Pi, E.; Zhu, Y. Analysis of EF-Hand proteins in soybean genome suggests their potential roles in environmental and nutritional stress signaling. Front. Plant Sci. 2017, 8, 877. [Google Scholar] [CrossRef] [PubMed]
  65. Poovaiah, B.W.; Du, L.; Wang, H.; Yang, T. Recent advances in calcium/calmodulin-mediated signaling with an emphasis on plant-microbe interactions. Plant Physiol. 2013, 163, 531–542. [Google Scholar] [CrossRef] [PubMed]
  66. Behera, S.; Long, Y.; Schmitz-Thom, I.; Wang, X.-P.; Zhang, C.; Li, H.; Steinhorst, L.; Manishankar, P.; Ren, X.-L.; Offenborn, J.N.; et al. Two spatially and temporally distinct Ca2+ signals convey Arabidopsis thaliana responses to K+ deficiency. New Phytol. 2017, 213, 739–750. [Google Scholar] [CrossRef] [PubMed]
  67. Riveras, E.; Alvarez, J.M.; Vidal, E.A.; Oses, C.; Vega, A.; Gutierrez, R.A. The Calcium Ion Is a Second Messenger in the Nitrate Signaling Pathway of Arabidopsis. Plant Physiol. 2015, 169, 1397–1404. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  68. Jiang, C.; Gao, X.; Liao, L.; Harberd, N.P.; Fu, X. Phosphate starvation root architecture and anthocyanin accumulation responses are modulated by the gibberellin-DELLA signaling pathway in Arabidopsis. Plant Physiol. 2007, 145, 1460–1470. [Google Scholar] [CrossRef] [PubMed]
  69. Kagale, S.; Divi, U.K.; Krochko, J.E.; Keller, W.A.; Krishna, P. Brassinosteroid confers tolerance in Arabidopsis thaliana and Brassica napus to a range of abiotic stresses. Planta 2007, 225, 353–364. [Google Scholar] [CrossRef] [PubMed]
  70. Chiera, J.; Thomas, J.; Rufty, T. Leaf initiation and development in soybean under phosphorus stress. J. Exp. Bot. 2002, 53, 473–481. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  71. Yi, K.; Wu, Z.; Zhou, J.; Du, L.; Guo, L.; Wu, Y.; Wu, P. OsPTF1, a novel transcription factor involved in tolerance to phosphate starvation in rice. Plant Physiol. 2005, 138, 2087–2096. [Google Scholar] [CrossRef] [PubMed]
  72. Chen, Y.F.; Li, L.Q.; Xu, Q.; Kong, Y.H.; Wang, H.; Wu, W.H. The WRKY6 transcription factor modulates PHOSPHATE1 expression in response to low Pi stress in Arabidopsis. Plant Cell 2009, 21, 3554–3566. [Google Scholar] [CrossRef] [PubMed]
  73. Ramaiah, M.; Jain, A.; Raghothama, K.G. Ethylene Response Factor070 regulates root development and phosphate starvation-mediated responses. Plant Physiol. 2014, 164, 1484–1498. [Google Scholar] [CrossRef] [PubMed]
  74. Ruan, W.; Guo, M.; Wu, P.; Yi, K. Phosphate starvation induced OsPHR4 mediates Pi-signaling and homeostasis in rice. Plant Mol. Biol. 2017, 93, 327–340. [Google Scholar] [CrossRef] [PubMed]
  75. Lan, P.; Li, W.; Schmidt, W. Genome-wide co-expression analysis predicts protein kinases as important regulators of phosphate deficiency-induced root hair remodeling in Arabidopsis. BMC Genom. 2013, 14, 210. [Google Scholar] [CrossRef] [PubMed]
  76. Bayle, V.; Arrighi, J.F.; Creff, A.; Nespoulous, C.; Vialaret, J.; Rossignol, M.; Gonzalez, E.; Paz-Ares, J.; Nussaume, L. Arabidopsis thaliana high-affinity phosphate transporters exhibit multiple levels of posttranslational regulation. Plant Cell 2011, 23, 1523–1535. [Google Scholar] [CrossRef] [PubMed]
  77. Gregory, A.L.; Hurley, B.A.; Tran, H.T.; Valentine, A.J.; She, Y.M.; Knowles, V.L.; Plaxton, W.C. In vivo regulatory phosphorylation of the phosphoenolpyruvate carboxylase AtPPC1 in phosphate-starved Arabidopsis thaliana. Biochem. J. 2009, 420, 57–65. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  78. Chen, J.; Wang, Y.; Wang, F.; Yang, J.; Gao, M.; Li, C.; Liu, Y.; Liu, Y.; Yamaji, N.; Ma, J.F.; et al. The rice CK2 kinase regulates trafficking of phosphate transporters in response to phosphate levels. Plant Cell 2015, 27, 711–723. [Google Scholar] [CrossRef] [PubMed]
  79. Gamuyao, R.; Chin, J.H.; Pariasca-Tanaka, J.; Pesaresi, P.; Catausan, S.; Dalid, C.; Slamet-Loedin, I.; Tecson-Mendoza, E.M.; Wissuwa, M.; Heuer, S. The protein kinase Pstol1 from traditional rice confers tolerance of phosphorus deficiency. Nature 2012, 488, 535–539. [Google Scholar] [CrossRef] [PubMed]
  80. Zhang, H.; Huang, L.; Hong, Y.; Song, F. BOTRYTIS-INDUCED KINASE1, a plasma membrane-localized receptor-like protein kinase, is a negative regulator of phosphate homeostasis in Arabidopsis thaliana. BMC Plant Biol. 2016, 16, 152. [Google Scholar] [CrossRef] [PubMed]
  81. Lei, L.; Li, Y.; Wang, Q.; Xu, J.; Chen, Y.; Yang, H.; Ren, D. Activation of MKK9-MPK3/MPK6 enhances phosphate acquisition in Arabidopsis thaliana. New Phytol. 2014, 203, 1146–1160. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  82. Xu, J.; Li, H.-D.; Chen, L.-Q.; Wang, Y.; Liu, L.-L.; He, L.; Wu, W.-H. A Protein Kinase, Interacting with Two Calcineurin B-like Proteins, Regulates K+ Transporter AKT1 in Arabidopsis. Cell 2006, 125, 1347–1360. [Google Scholar] [CrossRef] [PubMed]
  83. Ho, C.H.; Lin, S.H.; Hu, H.C.; Tsay, Y.F. CHL1 functions as a nitrate sensor in plants. Cell 2009, 138, 1184–1194. [Google Scholar] [CrossRef] [PubMed]
  84. Tenhaken, R. Cell wall remodeling under abiotic stress. Front. Plant Sci. 2015, 5, 771. [Google Scholar] [CrossRef] [PubMed]
  85. Langmead, B.; Salzberg, S.L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 2012, 9, 357–359. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  86. Trapnell, C.; Pachter, L.; Salzberg, S.L. TopHat: Discovering splice junctions with RNA-Seq. Bioinformatics 2009, 25, 1105–1111. [Google Scholar] [CrossRef] [PubMed]
  87. Trapnell, C.; Roberts, A.; Goff, L.; Pertea, G.; Kim, D.; Kelley, D.R.; Pimentel, H.; Salzberg, S.L.; Rinn, J.L.; Pachter, L. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat. Protoc. 2012, 7, 562. [Google Scholar] [CrossRef] [PubMed]
  88. Wang, G.; Zeng, H.; Hu, X.; Zhu, Y.; Chen, Y.; Shen, C.; Wang, H.; Poovaiah, B.W.; Du, L. Identification and expression analyses of calmodulin-binding transcription activator genes in soybean. Plant Soil 2015, 386, 205–221. [Google Scholar] [CrossRef]
  89. Pfaffl, M.W. A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res. 2001, 29, e45. [Google Scholar] [CrossRef] [PubMed]
  90. Ashburner, M.; Ball, C.A.; Blake, J.A.; Botstein, D.; Butler, H.; Cherry, J.M.; Davis, A.P.; Dolinski, K.; Dwight, S.S.; Eppig, J.T.; et al. Gene ontology: Tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000, 25, 25–29. [Google Scholar] [CrossRef] [PubMed]
  91. Usadel, B.; Poree, F.; Nagel, A.; Lohse, M.; Czedik-Eysenberg, A.; Stitt, M. A guide to using MapMan to visualize and compare Omics data in plants: A case study in the crop species, Maize. Plant Cell Environ. 2009, 32, 1211–1229. [Google Scholar] [CrossRef] [PubMed]
  92. Medina-Rivera, A.; Defrance, M.; Sand, O.; Herrmann, C.; Castro-Mondragon, J.A.; Delerce, J.; Jaeger, S.; Blanchet, C.; Vincens, P.; Caron, C.; et al. RSAT 2015: Regulatory Sequence Analysis Tools. Nucleic Acids Res. 2015, 43, W50–W56. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Figure 1. Pi concentration in the first trifoliate leaf (1st TL), the second TL (2nd TL), the third TL (3rd TL) and the roots of soybean seedlings after being subjected to Pi-deficiency treatment for 24 h. Data presented are mean ± SD (standard deviation) of three independent experiments. Asterisk indicates the statistically significant difference between control (+P) and Pi-deficiency (−P) treatments (Student’s t-test, * p < 0.05) (Microsoft Excel 2007).
Figure 1. Pi concentration in the first trifoliate leaf (1st TL), the second TL (2nd TL), the third TL (3rd TL) and the roots of soybean seedlings after being subjected to Pi-deficiency treatment for 24 h. Data presented are mean ± SD (standard deviation) of three independent experiments. Asterisk indicates the statistically significant difference between control (+P) and Pi-deficiency (−P) treatments (Student’s t-test, * p < 0.05) (Microsoft Excel 2007).
Ijms 19 02145 g001
Figure 2. Heatmap of the expression of 18 randomly selected differentially expressed genes (DEGs) as revealed by RNA-seq and qRT-PCR. The intensities of the color represent the fold changes in log2 values according to RNA-seq data or qRT-PCR results. Red color indicates induction and green color indicates repression. Accession number “Glyma.” is named as “Gm” for short. −P refers to PDL and +P refers to PSL.
Figure 2. Heatmap of the expression of 18 randomly selected differentially expressed genes (DEGs) as revealed by RNA-seq and qRT-PCR. The intensities of the color represent the fold changes in log2 values according to RNA-seq data or qRT-PCR results. Red color indicates induction and green color indicates repression. Accession number “Glyma.” is named as “Gm” for short. −P refers to PDL and +P refers to PSL.
Ijms 19 02145 g002
Figure 3. Unrooted phylogenetic tree of the SPX domain-containing proteins in soybean (Glycine max), Arabidopsis (Arabidopsis thaliana) and rice (Oryza sativa). The alignment for the phylogenetic tree was performed with ClustalW using full-length amino acid sequences. The phylogenetic tree was created with the MEGA6 software and the neighbor-joining method with 1000 bootstrap replications. The bar indicates the relative divergence of the sequences examined. The red arrow indicates the up-regulation of SPX genes upon short-term Pi deprivation by RNA-seq in the present study. Soybean SPX domain-containing proteins are marked with blue diamond.
Figure 3. Unrooted phylogenetic tree of the SPX domain-containing proteins in soybean (Glycine max), Arabidopsis (Arabidopsis thaliana) and rice (Oryza sativa). The alignment for the phylogenetic tree was performed with ClustalW using full-length amino acid sequences. The phylogenetic tree was created with the MEGA6 software and the neighbor-joining method with 1000 bootstrap replications. The bar indicates the relative divergence of the sequences examined. The red arrow indicates the up-regulation of SPX genes upon short-term Pi deprivation by RNA-seq in the present study. Soybean SPX domain-containing proteins are marked with blue diamond.
Ijms 19 02145 g003
Figure 4. Heatmap of DEGs that are putative protein kinases and phosphatases. The intensities of the color represent the fold changes in log2 values according to RNA-seq data. Red color indicates induction and green color indicates repression. Accession number “Glyma.” is named as “Gm” for short. −P refers to PDL and +P refers to PSL.
Figure 4. Heatmap of DEGs that are putative protein kinases and phosphatases. The intensities of the color represent the fold changes in log2 values according to RNA-seq data. Red color indicates induction and green color indicates repression. Accession number “Glyma.” is named as “Gm” for short. −P refers to PDL and +P refers to PSL.
Ijms 19 02145 g004
Figure 5. Mapman representation of DEGs involved in metabolism in soybean leaves upon short-term Pi deficiency. The results are the mean of two biological replicates. All results are shown on a log2 scale for corresponding FPKM ratios. Red color represents up-regulated genes and blue color represents down-regulated genes. TCA: tricarboxylic acid; OPP: oxidative pentose phosphate.
Figure 5. Mapman representation of DEGs involved in metabolism in soybean leaves upon short-term Pi deficiency. The results are the mean of two biological replicates. All results are shown on a log2 scale for corresponding FPKM ratios. Red color represents up-regulated genes and blue color represents down-regulated genes. TCA: tricarboxylic acid; OPP: oxidative pentose phosphate.
Ijms 19 02145 g005
Figure 6. The occurrence of Pi-related cis-elements previously identified as common to Pi-responsive genes in the promoter regions (1000 bp) of up-regulated DEGs, down-regulated DEGs and randomly-selected genes. Promoter regions of 250 genes randomly selected from all the chromosomes of soybean genome were used as control. 288 and 213 promoter regions were acquired for up-regulated and down-regulated DEGs, respectively. The hypergeometric p-value was calculated online (http://systems.crump.ucla.edu/hypergeometric/index.php). Asterisk means significantly different from the genes with elements in the soybean genome that are predicted based on the results of randomly-selected genes (* p < 0.05, ** p < 0.01).
Figure 6. The occurrence of Pi-related cis-elements previously identified as common to Pi-responsive genes in the promoter regions (1000 bp) of up-regulated DEGs, down-regulated DEGs and randomly-selected genes. Promoter regions of 250 genes randomly selected from all the chromosomes of soybean genome were used as control. 288 and 213 promoter regions were acquired for up-regulated and down-regulated DEGs, respectively. The hypergeometric p-value was calculated online (http://systems.crump.ucla.edu/hypergeometric/index.php). Asterisk means significantly different from the genes with elements in the soybean genome that are predicted based on the results of randomly-selected genes (* p < 0.05, ** p < 0.01).
Ijms 19 02145 g006
Table 1. Differentially expressed genes where transcripts were found only in Pi-deprived leaves (PDL) or Pi-sufficient leaves (PSL) after 24 h Pi deprivation. Expression level of each gene was measured in terms of fragments per kilobase of transcript per million mapped reads (FPKM). Gene functional descriptions were shown according to soybean genome annotation (V2.0). Accession number “Glyma.” is named as “Gm” for short. −P refers to PDL and +P refers to PSL. NA: no annotation.
Table 1. Differentially expressed genes where transcripts were found only in Pi-deprived leaves (PDL) or Pi-sufficient leaves (PSL) after 24 h Pi deprivation. Expression level of each gene was measured in terms of fragments per kilobase of transcript per million mapped reads (FPKM). Gene functional descriptions were shown according to soybean genome annotation (V2.0). Accession number “Glyma.” is named as “Gm” for short. −P refers to PDL and +P refers to PSL. NA: no annotation.
Accession No.−P (FPKM)+P (FPKM)p ValueDescription
Gm04g2234001260.240.001.98 × 10−2NA
Gm17g10000021.710.001.65 × 10−3mini zinc finger
Gm17g0904003.820.003.95× 10−3calcium-dependent protein kinase 6
Gm16g0283002.720.004.10 × 10−2NA
Gm15g0992002.270.002.00 × 10−4NA
Gm06g2409001.250.001.50 × 10−2PLATZ transcription factor family protein
Gm09g0230001.110.005.00 × 10−5peroxidase superfamily protein
Gm06g1772000.770.002.40 × 10−3NA
Gm19g2466000.750.001.00 × 10−4SSXT family protein
Gm13g0683000.0047.119.65 × 10−3NA
Gm05g0997000.0010.713.85 × 10−3NA
Gm04g1285000.005.877.50 × 10−4late embryogenesis abundant protein, group 1 protein
Gm13g2363000.005.779.00 × 10−4NA
Gm04g0140000.003.305.00 × 10−5polyol/monosaccharide transporter 5
Gm13g1780000.002.184.19 × 10−2NA
Gm07g0229000.002.118.65 × 10−3phospholipid N-methyltransferase
Gm04g0948000.001.215.00 × 10−5plant-specific transcription factor YABBY family protein
Gm11g1791000.001.025.00 × 10−5aluminum activated malate transporter family protein
Gm03g1136000.000.825.00 × 10−5auxin efflux carrier family protein
Gm04g0443000.000.801.61 × 10−2cytochrome B5 isoform D
Gm15g1973000.000.735.00 × 10−5CBL-interacting protein kinase 23
Table 2. DEGs potentially involved in Pi signaling and utilization. Expression level of each gene was measured in terms of FPKM. Gene functional descriptions were shown according to soybean genome annotation (V2.0). Accession number “Glyma.” is named as “Gm” for short. −P refers to PDL and +P refers to PSL. SPX: domain found in Syg1, Pho81, XPR1, and related proteins; HAD: haloacid dehydrogenase.
Table 2. DEGs potentially involved in Pi signaling and utilization. Expression level of each gene was measured in terms of FPKM. Gene functional descriptions were shown according to soybean genome annotation (V2.0). Accession number “Glyma.” is named as “Gm” for short. −P refers to PDL and +P refers to PSL. SPX: domain found in Syg1, Pho81, XPR1, and related proteins; HAD: haloacid dehydrogenase.
Accession No.−P (FPKM)+P (FPKM)Log2(−P/+P)p ValueDescription
Gm03g0783001.590.094.154.42 × 10−2Sulfoquinovosyldiacylglycerol 2 (GmSQD2)
Gm04g1476002.800.412.761.03 × 10−2SPX domain-containing protein 1-related (GmSPX3)
Gm08g0564009.381.432.715.00 × 10−5Purple acid phosphatase (GmPAP13)
Gm01g0918000.900.222.024.90 × 10−3Phosphate transporter PHO1 homolog 1 (GmPHO1;H12)
Gm02g1302001.280.322.026.50 × 10−4Phosphate transporter PHO1 homolog 1 (GmPHO1;H14)
Gm20g2380001.420.381.891.00 × 10−4Phospholipase D P1 (GmPLDZ2)
Gm06g0690002.680.911.562.52 × 10−2SPX domain-containing protein 1-related (GmSPX4)
Gm17g11470019.967.531.415.00 × 10−5SPX domain-containing protein 1-related (GmSPX8)
Gm08g1941000.710.291.292.20 × 10−2Phospholipase D α 4
Gm01g00240013.806.681.051.31 × 10−2Phospholipase A2 family protein
Gm09g22370034.4316.851.035.00 × 10−5Glycerol-3-phosphate permease gene family
Gm19g0266002.781.161.271.40 × 10−3Purple acid phosphatase (GmPAP31)
Gm16g2209002.144.75−1.151.64 × 10−2HAD superfamily, subfamily IIIB acid phosphatase
Table 3. DEGs potentially involved in the transportation of water, sugars and mineral nutrients. Expression level of each gene was measured in terms of FPKM. Gene functional descriptions were shown according to soybean genome annotation (V2.0). Accession number “Glyma.” is named as “Gm” for short. −P refers to PDL and +P refers to PSL. ATOX1: antioxidant protein 1; PIP1: plasma membrane intrinsic protein 1; TIP2: tonoplast intrinsic aquaporin 2; ATP: adenosine triphospate; ABC: ATP-binding cassette.
Table 3. DEGs potentially involved in the transportation of water, sugars and mineral nutrients. Expression level of each gene was measured in terms of FPKM. Gene functional descriptions were shown according to soybean genome annotation (V2.0). Accession number “Glyma.” is named as “Gm” for short. −P refers to PDL and +P refers to PSL. ATOX1: antioxidant protein 1; PIP1: plasma membrane intrinsic protein 1; TIP2: tonoplast intrinsic aquaporin 2; ATP: adenosine triphospate; ABC: ATP-binding cassette.
SubstrateAccession No.−P (FPKM)+P (FPKM)Log2(−P/+P)p ValueDescription
SugarGm15g1059001.480.471.642.78 × 102Glucose-6-phosphate/phosphate translocator 2
Gm13g2067006.262.771.189.00 × 104Glucose-6-phosphate/phosphate translocator 2
Gm11g226200131.2963.501.055.00 × 105Nucleotide-sugar transporter family protein
Gm04g0140000.003.30+P only5.00 × 105Polyol/monosaccharide transporter 5
Gm08g0099002.8019.84−2.835.00 × 105Bidirectional sugar transporter SWEET10
Zinc/ironGm16g1682007.322.141.781.40 × 103Vacuolar iron transporter (VIT) family protein
Gm13g0044001.337.82−2.555.00 × 105Zinc transporter 3-related protein
Gm20g0631001.809.06−2.335.00 × 105Zinc transporter 3-related protein
Gm06g05200024.5658.72−1.265.00 × 105Zinc/iron transporter
SulfateGm01g1752002.160.831.373.55 × 103Sulfite exporter TauE/SafE family protein
Gm15g0140001.780.841.081.60 × 102Sulfate transporter 3;4
CopperGm08g1805004.150.812.364.03 × 102Copper transport protein ATOX1-related
Gm14g0721006.812.581.405.00 × 105Copper transport protein ATOX1-related
Gm15g05190028.0213.361.071.85 × 103Copper transport protein ATOX1-related
NitrateGm18g1272001.073.08−1.534.50 × 104Nitrate transporter 1.7
MalateGm11g1791000.001.02+P only5.00 × 105Aluminum activated malate transporter
WaterGm08g0153008.723.581.293.50 × 104Aquaporin PIP1;4-related
Gm05g20870015.226.501.231.50 × 104Aquaporin PIP1;4-related
Gm15g01810069.6533.301.065.00 × 105Aquaporin TIP2;1
UnknownGm07g08100018.374.442.055.00 × 105ABC transporter family protein
Gm11g10630054.1517.681.615.00 × 105Major facilitator superfamily protein
Gm04g0278001.280.351.862.71 × 102Protein Walls Are Thin 1
Gm07g0940000.432.97−2.775.30 × 103ATP-binding cassette A2
Gm11g0641000.872.05−1.232.15 × 103Major facilitator superfamily protein
Gm08g1015001.012.36−1.225.00 × 105Multidrug resistance-associated protein 3
Gm17g1652005.6811.51−1.025.00 × 105Metal transporter NRAMP2-related
Gm10g2767002.394.79−1.006.00 × 104Major facilitator superfamily protein
Table 4. DEGs possibly related to Ca2+ and hormonal signaling. Expression level of each gene was measured in terms of FPKM. Gene functional descriptions were shown according to soybean genome annotation (V2.0). Accession number “Glyma.” is named as “Gm” for short. −P refers to PDL and +P refers to PSL. VCX1: vacuolar Ca2+ exchanger 1; AUX: auxin; IAA: indole-3-acetic acid; GA: gibberellin; BES1/BZR1: brassinosteroid-insensitive1-ethyl methanesulfonate-suppressor 1/brassinazole-resistant 1.
Table 4. DEGs possibly related to Ca2+ and hormonal signaling. Expression level of each gene was measured in terms of FPKM. Gene functional descriptions were shown according to soybean genome annotation (V2.0). Accession number “Glyma.” is named as “Gm” for short. −P refers to PDL and +P refers to PSL. VCX1: vacuolar Ca2+ exchanger 1; AUX: auxin; IAA: indole-3-acetic acid; GA: gibberellin; BES1/BZR1: brassinosteroid-insensitive1-ethyl methanesulfonate-suppressor 1/brassinazole-resistant 1.
Ca2+/HormoneAccession No.−P (FPKM)+P (FPKM)Log2(−P/+P)p ValueDescription
CalciumGm02g2265001.100.242.192.14 × 102Ca2+/H+ antiporter VCX1 and related protein
Gm14g0001001.596.01−1.925.00 × 105Glutamate receptor 3.1-related protein
Gm16g0618000.370.83−1.151.52 × 102Glutamate receptor 2.5
Gm07g2037001.583.28−1.062.50 × 104Glutamate receptor 2.7
Gm17g0904003.820.00−P only3.95 × 103Calcium-dependent protein kinase 6
Gm17g1309004.201.381.605.00 × 104Calmodulin-binding protein IQ-domain 22
Gm09g2253001.100.371.571.30 × 102Protein IQ-domain 11
Gm01g00240013.806.681.051.31 × 102Phospholipase A2 family protein
Gm18g2308001.072.74−1.365.00 × 105Ca2+/lipid-binding phosphoribosyltransferase
Gm09g2368008.0716.28−1.011.70 × 103Ca2+-binding EF-hand family protein (GmCML116)
AuxinGm16g1155002.691.151.231.65 × 102Auxin efflux carrier family protein
Gm03g1136000.000.82+P only5.00 × 105Auxin efflux carrier family protein
Gm19g1685000.501.59−1.662.75 × 102Auxin-responsive protein IAA10-related
Gm19g1611009.8225.41−1.375.00 × 105AUX/IAA family
CytokininGm04g2478003.321.031.701.80 × 102Response regulator of two-component system
Gm05g1445000.461.09−1.221.58 × 102Response regulator of two-component system
GibberellinGm06g044400439.6189.832.295.00 × 105Gibberellin-regulated family protein
Gm06g18530071.4215.602.191.00 × 104Gibberellin-regulated family protein
Gm05g03450019.806.631.581.89 × 102Gibberellin-regulated family protein
Gm06g193800174.5371.571.292.00 × 104Gibberellin-regulated family protein
Gm17g09280017.657.921.164.94 × 102Gibberellin-regulated family protein
Gm15g0022005.272.621.011.35 × 103GA requiring 3
BrassinosteroidGm14g1274001.360.242.521.28 × 102BES1/BZR1 homolog protein 1
JasmonateGm01g2044004.221.711.304.99 × 102JASMONATE-ZIM-domain protein 2-related
EthyleneGm10g0085002.204.84−1.142.61 × 102Ethylene response sensor 2-related
Table 5. Transcription factor genes differentially expressed in soybean leaves under short-term Pi deprivation. Expression level of each gene was measured in terms of FPKM. Gene functional descriptions were shown according to soybean genome annotation (V2.0). Accession number “Glyma.” is named as “Gm” for short. −P refers to PDL and +P refers to PSL. MYB: (myeloblastosis); DOF: DNA-binding one zinc finger; YABBY: YABBY domain;CCT: CONSTANS, CO-like, and TOC1; NAC: NAM, ATAF, and CUC; TCP: a family of transcription factors named after: teosinte branched 1 (tb1, Zea mays (Maize)), cycloidea (cyc) (Antirrhinum majus) (Garden snapdragon) and PCF in rice (Oryza sativa); ERF: ethylene response factor; WRKY: a protein domain composed of a conserved WRKYGQK motif; RING: really interesting new gene.
Table 5. Transcription factor genes differentially expressed in soybean leaves under short-term Pi deprivation. Expression level of each gene was measured in terms of FPKM. Gene functional descriptions were shown according to soybean genome annotation (V2.0). Accession number “Glyma.” is named as “Gm” for short. −P refers to PDL and +P refers to PSL. MYB: (myeloblastosis); DOF: DNA-binding one zinc finger; YABBY: YABBY domain;CCT: CONSTANS, CO-like, and TOC1; NAC: NAM, ATAF, and CUC; TCP: a family of transcription factors named after: teosinte branched 1 (tb1, Zea mays (Maize)), cycloidea (cyc) (Antirrhinum majus) (Garden snapdragon) and PCF in rice (Oryza sativa); ERF: ethylene response factor; WRKY: a protein domain composed of a conserved WRKYGQK motif; RING: really interesting new gene.
Accession No.−P (FPKM)+P (FPKM)Log2(−P/+P)p ValueDescription
Gm14g1274001.360.242.521.28 × 10−2BES1/BZR1 homolog protein 1-related
Gm17g2577001.250.312.004.86 × 10−2C2H2-like zinc finger protein
Gm08g2154001.860.471.993.36 × 10−2Basic helix-loop-helix (bHLH) DNA-binding superfamily protein
Gm06g2849001.930.541.823.15 × 10−3Basic-leucine zipper (bZIP) transcription factor family protein
Gm16g0177004.931.561.662.50 × 10−4Basic helix-loop-helix (bHLH) DNA-binding superfamily protein
Gm20g18650012.124.341.485.00 × 10−5MYB-like transcription factor family protein
Gm12g0372007.122.781.361.59 × 10−2B-box type zinc finger family protein
Gm07g0121008.873.591.302.10 × 10−3DOF zinc finger protein 1
Gm18g0423002.280.941.295.35 × 10−3Zinc finger protein 4
Gm13g0619001.190.511.224.94 × 10−2MYB domain protein 16
Gm05g0560005.892.621.171.16 × 10−2Plant-specific transcription factor YABBY family protein
Gm01g1335007.163.431.064.33 × 10−2Basic helix-loop-helix (bHLH) DNA-binding superfamily protein
Gm11g1171006.843.351.031.05 × 10−3Basic helix-loop-helix (bHLH) DNA-binding superfamily protein
Gm04g0948000.001.21+P only5.00 × 10−5Plant-specific transcription factor YABBY family protein
Gm13g2037000.131.03−3.024.61 × 10−2C2H2-type zinc finger family protein
Gm01g0293001.376.15−2.163.77 × 10−2Plant-specific transcription factor YABBY family protein
Gm10g0214000.351.44−2.062.51 × 10−2B-box type zinc finger protein with CCT domain
Gm11g1820000.923.34−1.871.33 × 10−2NAC domain containing protein 90
Gm13g0474000.311.10−1.851.38 × 10−2TCP family transcription factor
Gm12g1627000.622.11−1.762.43 × 10−2ERF domain protein 9
Gm02g1778000.852.71−1.673.78 × 10−2MYB-like HTH transcriptional regulator family protein
Gm05g0027000.290.89−1.624.98 × 10−2NAC (No Apical Meristem) domain transcriptional regulator
Gm01g0417001.092.98−1.467.15 × 10−3Homeobox 1
Gm02g2178004.2711.31−1.401.00 × 10−4Transcription factor bHLH61-related
Gm03g2217001.403.66−1.392.28 × 10−2MYB domain protein 2
Gm08g2564001.433.54−1.317.20 × 10−3TCP family transcription factor
Gm02g1529000.942.25−1.252.90 × 10−2B-box type zinc finger protein with CCT domain
Gm05g1445000.461.09−1.221.58 × 10−2Response regulator 11
Gm07g2272006.3412.74−1.011.50 × 10−4WRKY DNA-binding protein 3
Gm08g2860001.412.83−1.001.25 × 10−2Zinc finger (CCCH-type/C3HC4-type RING finger) family protein

Share and Cite

MDPI and ACS Style

Zeng, H.; Zhang, X.; Zhang, X.; Pi, E.; Xiao, L.; Zhu, Y. Early Transcriptomic Response to Phosphate Deprivation in Soybean Leaves as Revealed by RNA-Sequencing. Int. J. Mol. Sci. 2018, 19, 2145. https://doi.org/10.3390/ijms19072145

AMA Style

Zeng H, Zhang X, Zhang X, Pi E, Xiao L, Zhu Y. Early Transcriptomic Response to Phosphate Deprivation in Soybean Leaves as Revealed by RNA-Sequencing. International Journal of Molecular Sciences. 2018; 19(7):2145. https://doi.org/10.3390/ijms19072145

Chicago/Turabian Style

Zeng, Houqing, Xiajun Zhang, Xin Zhang, Erxu Pi, Liang Xiao, and Yiyong Zhu. 2018. "Early Transcriptomic Response to Phosphate Deprivation in Soybean Leaves as Revealed by RNA-Sequencing" International Journal of Molecular Sciences 19, no. 7: 2145. https://doi.org/10.3390/ijms19072145

APA Style

Zeng, H., Zhang, X., Zhang, X., Pi, E., Xiao, L., & Zhu, Y. (2018). Early Transcriptomic Response to Phosphate Deprivation in Soybean Leaves as Revealed by RNA-Sequencing. International Journal of Molecular Sciences, 19(7), 2145. https://doi.org/10.3390/ijms19072145

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop