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

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.


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 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.

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 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. 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.

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.

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.   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.

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. 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.

Class III (SPX-MFS)
Class IV (SPX-RING) 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.

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.

Genes Linked to Ca 2+ 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 Ca 2+ 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.

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.

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. 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.   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. 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.

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, 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.

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. 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.

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

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. 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.

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, 0% 5% 10% 15% 20% 25% 30% 35% 40% 45% 50% Percentage of genes with elements Up-regulated genes Down-regulated genes Randomly-selected genes ** * * * 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).

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 Ca 2+ signaling (Table 4). Ca 2+ 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 Ca 2+ concentration, which can be perceived and decoded by a series of Ca 2+ sensors containing EF-hand motifs [64,65]. However, little information is available with respect to the effect of Pi deficiency on cytosolic Ca 2+ levels. Considering the critical roles of the Ca 2+ signal in plant responses to other nutrient stress, such as potassium and nitrate deficiencies [66,67], it is conceivable that Ca 2+ signals might be an important player in Pi stress response. In addition, Ca 2+ and Pi are incompatible ions, because Ca 2+ can form insoluble compounds with phosphate derivatives at high levels. Changing the allocation patterns of Pi may also necessitate changes to the patterns for Ca 2+ . Thus, how Ca 2+ 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.

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(NO 3

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% (NH 4 ) 6 Mo 7 O 24 ·4H 2 O, 0.46 M H 2 SO 4 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).

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.

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).

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  Table S16.

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).