Differential Gene Expression Responding to Low Phosphate Stress in Leaves and Roots of Maize by cDNA-SRAP

Phosphate (Pi) deficiency in soil can have severe impacts on the growth, development, and production of maize worldwide. In this study, a cDNA-sequence-related amplified polymorphism (cDNA-SRAP) transcript profiling technique was used to evaluate the gene expression in leaves and roots of maize under Pi stress for seven days. A total of 2494 differentially expressed fragments (DEFs) were identified in response to Pi starvation with 1202 and 1292 DEFs in leaves and roots, respectively, using a total of 60 primer pairs in the cDNA-SRAP analysis. These DEFs were categorized into 13 differential gene expression patterns. Results of sequencing and functional analysis showed that 63 DEFs (33 in leaves and 30 in roots) were annotated to a total of 54 genes involved in diverse groups of biological pathways, including metabolism, photosynthesis, signal transduction, transcription, transport, cellular processes, genetic information, and organismal system. This study demonstrated that (1) the cDNA-SRAP transcriptomic profiling technique is a powerful method to analyze differential gene expression in maize showing advantageous features among several transcriptomic methods; (2) maize undergoes a complex adaptive process in response to low Pi stress; and (3) a total of seven differentially expressed genes were identified in response to low Pi stress in leaves or roots of maize and could be used in the genetic modification of maize.


Introduction
Maize is a well-known economically important crop in the world. However, the maize production is commonly affected by phosphate (Pi) deficiency in soil due to the high percentage of fixed Pi that cannot be absorbed by plants [1]. Plants have evolved many adaptive strategies to cope with low Pi supply with many genes involved in the response to Pi starvation [2]. For example, the Pi uptake is increased by directly inducing the expression of Pi transporters, while the inositol polyphosphates are shown to play important roles in the intracellular Pi signaling [3][4][5]. In order to understand the molecular mechanisms responding to low Pi stress in maize, it is becoming increasingly important to characterize the transcriptional changes and signal transduction pathways in plants under low Pi stress.
Previous studies have investigated gene expression patterns in response to plant defense against low Pi [6][7][8][9][10][11][12]. For example, microRNAs (e.g., miRNA399b and Zma-miR3) are induced to express in the roots of maize under low Pi stress [12], while studies of metabolite profiling and genome-wide association have shown that the gene encoding glucose-6-phosphate-1-epimerase is capable of correlating traits based on the Pi levels in maize seedlings [9]. Furthermore, abscisic acid (ABA) plays a dominant role in regulating candidate genes in response to Pi stress in Arabidopsis [11].
The initial response to Pi starvation in plants includes the sensation and communication of the amount of Pi at both the cellular and molecular levels in various organs (e.g., roots and leaves) with the roots establishing the initial response to low Pi stress [2]. Studies show that Pi starvation promotes growth of primary roots but reduces the numbers of lateral roots and lateral root primordia, while the Pi intake under Pi stress is enhanced by greater lateral root branching density [13][14][15]. Transcriptomic analyses of maize gene expression under Pi stress have been carried out in maize roots [6,7,11,13]. These studies showed that auxin signaling is involved in the response to Pi stress, while the lateral root development is regulated by DNA replication, transcription, protein synthesis and degradation, and cell growth [13].
Plant survival under Pi stress depends on the tolerance at the whole-plant level rather than individual tissues or organs because low Pi usually leads to detrimental effects on the entire plant. For example, studies have shown the accumulations of the flavonoid pigment anthocyanin, di-and trisaccharides, and metabolites of ammonium metabolism in leaves under low Pi stress [6,16]. A total of 78 differentially expressed fragments (DEFs) in leaves and roots of maize have been identified using the cDNA-amplified fragment length polymorphism (cDNA-AFLP) technique [7], while the induction of acid phosphatase and activities of catalase/peroxidase in both stems and roots of Pi-tolerant cultivar of maize under Pi stress is identified much earlier than that in the Pi-sensitive cultivar [6]. Recent studies have shown that the nucleic acids, organic acids, and sugars were increased, while the phosphorylated metabolites, certain amino acids, lipid metabolites, and nitrogenous compounds were decreased in maize seedlings under low Pi stress [9]. Furthermore, it has been reported that the peroxidase and superoxide dismutaserelated genes and the lipid peroxidation genes were upregulated in maize leaves [17]. Therefore, simultaneous analyses of gene expression patterns in response to Pi stress in multiple plant organs, e.g., roots and leaves, provide more accurate understanding of the molecular mechanisms of Pi tolerance in maize.
This study focused on the gene expression in both roots and leaves of a maize inbred line in response to Pi stress using the sequence-related amplified polymorphism (SRAP) technique [18]. Due to its advantages of being easy, effective, and fast, the SRAP analysis has been widely used to construct a genetic map and to investigate molecular diversity and comparative genomics in various species of crops [18][19][20][21][22]. Besides the genomic DNA, cDNA can also be used as templates to identify the fragment length polymorphism using SRAP primers. For example, the cDNA-SRAP technique has been utilized successfully to study the differential gene expression in several species of plants [22][23][24][25][26]. In our study, we used this simple and effective method with low cost to study the gene expression patterns in both roots and leaves of maize under low Pi stress. Our goals were to (1) identify the DEFs that are transcriptionally regulated in response to Pi stress in maize and (2) to further identify the candidate genes of maize involved in response to low Pi stress.

Plant Material and Growth
Conditions. The previously isolated maize inbred line 478 showing high efficiency of intaking Pi was used in this study [27]. Seeds were surfacesterilized with 75% ethanol for 30 sec, rinsed three times with sterile distilled water, and planted in boxes filled with steril-ized perlite. A total of 40 3-day-old seedlings were transferred to a continuously aerated nutrient solution with high Pi (HP) content (1 mM KH 2 PO 4 ) for 7 days. Then, 20 of the seedlings were transferred into the low Pi (LP) solution (0.01 mM KH 2 PO 4 )), while the other 20 seedlings were maintained in the HP solution (as the control group). The basal nutrient solution contained 2 mM Ca(NO 3 ) 2 , 0.65 mM MgSO 4 , 25 μM Fe-EDTA, 5 μM MnSO 4 , 50μM KCl, 2 μM ZnSO 4 , 0.5 μM CuSO 4 , 0.005 μM (NH 4 ) 6 Mo 24 , and 25 mΜ H 3 BO 4 . In addition, K + was supplied in the LP solution in the form of 0.99 mM KCl. The seedlings were kept in an artificial climate chamber under a photoperiod of 8 h of darkness (18°C) and 16 h of light (28°C; 100 μmol m −2 s −1 photon flux density). Humidity was maintained at~70%.

Sample
Collection and RNA Preparation. Samples of leaves and roots were harvested at 0 day from the HP solution and at 3, 5, and 7 days from the LP solution, respectively. These eight samples were frozen immediately using liquid nitrogen and stored at -80°C for further analysis. Leaf or root tissues were ground to fine powder in liquid nitrogen to extract the total RNA with Trizol Reagent (Sangon Biotech, Shanghai, China) according to the manufacturer's instructions.
2.3. cDNA-SRAP Analysis. Synthesis of the first cDNA strand was performed with M-MLV Reverse Transcriptase (Promega). Each PCR reaction mixture (20 μL) of the cDNA-SRAP amplification contained 140 ng of cDNA, 0.25 mmol/L of dNTP mixture, 1.8 mmol/L of Mg 2+ , 0.5 μmol/L of forward and reverse primers, respectively, and 1.0 U of Taq DNA polymerase (Takara) [28]. In order to reduce the amplifications of nonspecific fragments in the cDNA-SRAP analysis, the conditions of PCR amplification were optimized to be 35 cycles of PCR amplification with a 20-fold dilution of DNA template. Further analyses were based on a total of 60 pairs of SRAP primers showing reliable amplifications indicated by rich and strong bands ( Table 1). The SRAP thermal cycling using PCR was as follows: an initial predenaturation at 94°C for 5 min, followed by five cycles of denaturation at 94°C for 30 sec, annealing at 35°C for 30 sec, and extension at 72°C for 1 min. Then, the annealing temperature was increased to 50°C for another 35 cycles with a terminal extension at 72°C for 10 min. The PCR products were separated using 6% urea-polyacrylamide gel electrophoresis (PAGE) and examined using the silver-staining method. For each primer combination, the final PCR products from a series of days (i.e., 0, 3, 5, and 7) of treatments of Pi stress were loaded in order into lanes next to each other in the PAGE gel for the comparison of band density for bands of the same size.
2.4. Identification, Isolation, and Sequencing of Differentially Expressed Fragments (DEFs). The slices of the PAGE gel containing the target bands showing high expression, over 200 bp in length, and representing one of the identified differential expression patterns were excised and then hydrated in 50 μL of Tris-EDTA buffer (pH 8.0) and incubated overnight at 37°C. The eluted fragments were reamplified using PCR 2 BioMed Research International with the same primers and the same conditions as those for the cDNA-SRAP analysis. The reamplified products were cloned using the pMD18-T vector (Takara) and then sequenced. These sequences of the transcript-derived fragments were compared to the protein database using the BLASTx algorithm at the National Center of Biotechnology Information (https://blast.ncbi.nlm.nih.gov/Blast.cgi). These sequences were also annotated using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (https://www .genome.jp/kegg/) to reveal the biochemical and physiological pathways.

Relative Quantitative Real-Time PCR (qRT-PCR)
Analysis. To further verify the results of the cDNA-SRAP analysis, the cDNAs derived from the eight leaf and root samples of maize were randomly selected to quantitatively examine the expression levels of selected fragments using the qRT-PCR. The amplification mixture contained 1 μL of cDNA, 10 μL of 2x SYBR® Premix Ex Taq™ II (Perfect Real Time) (Takara), and 0.5 μmol/L forward and reverse primers, respectively. The thermal cycling protocol consisted of preincubation at 95°C for 3 min, followed by 40 cycles of denaturation at 95°C for 10 s, and annealing at 60°C for 20 s. A standard melting curve was generated at the end of the amplification with the measurements recorded between 60°C and 95°C used to calculate the PCR efficiency (E). The qRT-PCR amplifications were conducted in parallel in triplicate with the normalization performed using the transcript level of the constitutively expressed Tubulin gene as control in all samples. Relative expression ratios were calculated according to Advanced Relative Quantification of Roche LightCycler 480 software release 1.5.0.
In comparison to other transcriptomic techniques, the cDNA-SRAP method has been widely used to study the differential gene expression because it does not have the disadvantages of the other methods and it has several unique advantages. First, the cDNA-SRAP analysis is a valuable tool for high-throughput gene expression analysis because of its generation of reliable and repeatable gene expression data and identification of a high coverage of transcript tags. Second, in contrast to the hybridization-based techniques, the cDNA-SRAP method allows a more thorough and comprehensive gene expression analysis at a large scale because of its capability of direct, simultaneous comparison of gene expression levels in different organs or tissues at different developmental stages under the treatment of various types of stress investigations. Third, in comparison to the cDNA-AFLP method, the cDNA-SRAP analysis costs less but is easier to operate because it does not require double-enzyme digestion of cDNA. The effectiveness of the cDNA-SRAP method is demonstrated evidently in our study showing the annotations of the cDNA fragments identified in both roots and leaves of maize to the diverse classes of known genes.

Expression Patterns of DEFs Revealed by cDNA-SRAP.
The cDNA-SRAP transcript profiling in both roots and leaves of maize under low Pi stress over 0, 3, 5, and 7 days was performed using 60 primer pairs to investigate the gene expression patterns in response to low Pi stress in maize. The sizes of a total of 2571 DNA fragments that we revealed ranged from 80 to 800 bp with an average of~43 fragments obtained per primer pair. A total of 2494 identified fragments (97%) were identified as DEFs with an average of~42 DEFs obtained per primer pair. These DEFs were categorized into 13 classes of differential expression patterns, while the nondifferentially expressed fragments were grouped in Class 14 ( Table 2) In leaves, there were a total of 994 DEFs (~80%) revealed in Class 1-5 expression patterns with each class containing more than 80 DEFs. Class 13 contained 67 DEFs, while less than 50 DEFs were identified in each of the other eight classes. Similar distributions of DEFs were identified in roots as those in leaves. A total of 960 DEFs (~72%) were identified in four classes (1, 2, 4, and 5) of expression patterns with each class containing more than 90 DEFs, while less than 60 DEFs were categorized in each of the other nine classes.
Studies using other transcriptomic analyses to investigate the response of maize to low Pi stress identified varied numbers of DEFs. For example, microarray analysis of the lateral root primordium zone responding to low Pi showed that in 2 days of treatment, a total of 148 differentially expressed transcripts contained 71 upregulated and 77 downregulated transcripts, and in 8 days of low Pi treatment, a total of 549 DEGs contained 270 upregulated and 279 downregulated genes in the roots of maize [13]. Studies using strand-specific RNA-Seq transcriptomic analyses of leaves and roots of low Ptolerant and P-sensitive maize inbred lines identified a total of 5900 DEFs in the low-P-sensitive and 3389 DEFs in the low-P-tolerant maize [6]. A total of 142 DEFs with 121 upregulated and 21 downregulated in LP were detected using the cDNA-AFLP method with 136 primer pairs to screen LPtolerant maize under HP and LP stress conditions [7]. In comparison to our results, these studies suggested that cDNA-SRAP analysis is more powerful and advantageous than other transcriptomic analyses in identifying the number of DEFs in maize responding to low Pi stress, while our results provide a more complete and integrated variations at the transcriptional level between the treatment of HP and LP stress.

Sequence Analysis and Functional Annotation.
A total of 63 DEFs (33 in leaves and 30 in roots) representing 11 out of the 13 classes of differential expression patterns with high expression and longer than 200 bp in length were selected and isolated from the polyacrylamide gels, reamplified by PCR, and sequenced (GenBank accession numbers JZ983140-JZ983202). The sequences of the 63 DEFs were annotated based on the GenBank protein database at the National Center of Biotechnology Information (https://blast.ncbi.nlm.nih.gov/Blast.cgi) using the BLASTx algorithm (Table 3). Significant homology (Evalue: 0.001) was revealed for 60 DEFs (~95%) to genes with known functions, while 3 DEFs (~5%) were annotated to genes without noting that different DEFs of each one of these four genes showed compatible rather than conflicting expression patterns, probably because these genes were not detected due to either having a low expression or containing multiple copies of transcripts, suggesting that the different expression patterns existed for different organs (e.g., roots or leaves) under various treatments (i.e., the length of treatments under low Pi stress) in relation to the distance from the sites of treatments, e.g., the quick response in roots, while the delayed response in leaves.
These 63 DEFs in maize under low Pi stress were further annotated into eight metabolic pathways with a group of "unclassified" based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (https://www.genome.jp/ kegg/). The "unclassified" group contained the largest number of DEFs (~35%), followed by the group of "metabolism" (~24%) containing proteins of unknown functions and hypothetical or unclassified proteins (Figure 1). The lowest numbers of DEFs (~1.5%) were revealed in three pathways including "cellular processes," "organismal systems," and "photosynthesis." These results indicated that the identified DEFs were involved in diverse groups of biological pathways in response to low Pi stress in maize. It was noted that the varied numbers of DEFs were revealed in different pathways in leaves and roots, suggesting that there are alternative molecular mechanisms responding to low Pi stress in leaves and roots of maize. Furthermore, previous studies of the response of maize to low Pi stress have showed that the pathway of "metabolism" contained the most number of DEFs [7,8,13], suggesting that "metabolism" is likely one of the most important pathways involved in response to low Pi stress in both roots and leaves of maize. These results further indicate that among these identified pathways, "signal transduction" and "genetic information processing" may also play important roles in response to low Pi stress in leaves and roots of maize, respectively. To summarize, many biological pathways are involved in the response to low Pi stress in maize, ultimately affecting its development, growth, and production; therefore, it is important to investigate comprehensively the molecular mechanisms responding to the low Pi stress in maize.

DEFs in Response to Low Pi Stress in Maize.
Although roots are the first organs in plants to respond to low Pi, plant survival depends on whole-plant tolerance to low Pi stress. Therefore, it is important to investigate both the root and the leaf in order to achieve a comprehensive understanding of the molecular mechanisms responding to low Pi stress in plants. Previous studies detected the total number of DEFs in maize under low Pi stress without any differentiation between the numbers of DEFs identified specifically from either leaves or roots [7,8,13]. Furthermore, tissue-specific gene expression patterns under low Pi conditions were revealed. For example, genes in auxin biosynthesis and signaling and cell defense response protein degradation were upregulated and the expression of genes involved in cell proliferation and growth decreased, while the meristem region and other tissues related to genes not differentially expressed in the lateral root primordium zone of the primary root [13].    Our results showed that there were more DEFs in roots than those in leaves in seven out of the 13 classes of DEFs, including the transient-induced (both 5 and 7 days), induced, switched-off-induced, induced-switched-off-induced, down, up, and up-down-up expression patterns. More DEFs was also realized in roots than those in leaves showing the constitutive expression pattern. Comparative proteomic analyses of Pi responses in the roots of maize revealed that pathways of citrate secretion, sugar metabolism, and root-cell proliferation played important roles in enhancing tolerance to low Pi conditions [8]. In comparison to both Arabidopsis thaliana and rice, maize showed different responsive patterns to LP stress; specifically, the growth of primary roots was promoted, and the formation of lateral roots was inhibited in maize [13]. Using the cDNA-AFLP analysis, Jiang et al. [7] identified a total of 78 DEFs and 9 genes in maize highly homologous to genes in both Arabidopsis thaliana and Oryza sativa involved in phosphorus metabolism, including phosphorus circulation, transportation, or response. Our results showed that there were more DEFs showing induced than down expressions in roots, while more DEFs were identified with transient-induced, rapidly switched-off, and induced expressions in roots than those in leaves, indicating that roots may respond more strongly to low Pi stress than leaves at transcriptional level.   BioMed Research International analyses. The correlation coefficient (R 2 ) for the tublin reference gene was 0.998, while the PCR efficiency (E) calculated from the slope of the standard melting curve was 105%. The R 2 and E values for the 8 target gene transcripts varied from 0.991 to 1 and from 94% to 133%, respectively. The results of quantitative expression of these 8 DEFs were consistent with the expression patterns revealed by the cDNA-SRAP analysis (data not shown). For example, there were two bands on the PAGE gels of DEFs 120 and 181, suggesting that there might be two transcripts or homologous genes. This was confirmed by the DEF 181 annotating the vacuolar ATPase subunit H protein because there were indeed two vacuolar ATPase subunit H proteins in maize [34]. The verification of the cDNA-SRAP analysis provided by the qRT-PCR method strongly indicate that the cDNA-SRAP analysis is reliable, easy to operate, and cost efficient to study the differential gene expressions at a large scale in plants.

Candidate Genes in Response to Low Pi Stress in Maize.
Our results showed that among the 54 genes annotated by the 63 DEFs (33 in leaves and 30 in roots), there were a total of seven genes found to be expressed differentially in both roots and leaves of maize under low Pi stress (Table 3). These genes were annotated in most of the biological pathways classified by the KEGG database, except for pathway "organismal systems." Among these seven genes, the gene encoding the photosystem II protein H was induced in both roots and leaves of maize (i.e., induced expression in leaves and updown expression in roots), while the other six genes were expressed in opposite expression patterns in leaves and in roots, suggesting different molecular mechanisms responding to the low Pi stress in leaves and roots, respectively.
These seven differentially expressed genes (DEGs) in roots and leaves of maize under low Pi stress annotated with known functions were considered as the candidate genes of maize involved in response to low Pi stress (Table 3). These candidate genes encode aspartic proteinaseo oryzasin-1 (annotated by DEFs 005 and 013 in the pathway of "signal transduction"), bZIP transcription factor ABI5 (annotated by DEFs 093, 094, 096, and 101 in the pathway of "transcription"), copia-type pol polyprotein (annotated by DEFs 018 and 243 as "unclassified"), photosystem II protein H (annotated by DEFs 115 and 139 in the pathway of "photosynthesis"), rRNA intron-encoded homing endonuclease (annotated by DEFs 262, 267, and 268 in the pathway of "genetic information processing"), vacuolar ATPase subunit H protein (annotated by DEFs 176 and 181 in the pathway of "transport"), and violaxanthin de-epoxidaserelated (annotated by DEFs 231 and 235 in the pathway of "metabolism"). These DEFs showed different expression patterns in roots and leaves in response to low Pi stress. For example, the opposite expression patterns were observed for genes encoding vacuolar ATPase subunit H protein and violaxanthin de-epoxidase-related in roots and leaves. These results suggested that these candidate genes play important roles in response to low Pi stress in maize. No previous studies have reported the functions of these genes in response to low Pi stress in maize. Further function analysis of these can-didate genes is necessary to help investigate the molecular mechanism of maize responding to the low Pi stress. It was noted that other genes identified by the cDNA-SRAP analysis were also candidate genes involved in the response to low Pi stress. Further studies were needed to characterize the roles these genes play in responding to Pi starvation. Here, we discuss the potential functions these genes play in response to low Pi stress.
The aspartic proteinase oryzasin-1 was expressed during seed ripening and germination in rice [35]. Aspartic protease in guard cell 1 (ASPG1) is suggested to be involved in ABAdependent responsiveness and the overexpression of the ASPG1 gene can confer drought avoidance in Arabidopsis [36]. Therefore, we speculated that the aspartic proteinase oryzasin-1 is probably involved in ABA-dependent responsiveness to low Pi stress in maize.
Our results revealed for the first time that genes encoding photosystem II protein H showed differential expression patterns in maize under LP treatments. Previous studies identified DEGs related to photosynthesis with downregulation in maize under low Pi treatment [17], while our results indicated that the photosystem II protein H showed induced expression pattern in leaves and up-down regulation in root, suggesting that the photosystem II protein H may play an important role in photosynthesis in response to low Pi stress.
Previous studies identified the functions of the Pideficiency-induced long-noncoding RNA1 (PILNCR1) in the inhibitions of ZmmiR399-guided cleavage of ZmPHO2 [37]. We report for the first time that the rRNA intronencoded homing endonuclease is involved in the response to the low Pi stress, while further investigations are needed to reveal the exact molecular mechanisms of how it participates in the response to low Pi stress in maize.
The vacuolar ATPase plays many key roles in plant growth and development and in stress response, while subunit H is vital to the activity and stability of the vacuolar ATPase [37,38]. Two vacuolar ATPase subunit H protein genes are reported in maize (ZmVHA-H1 and ZmVHA-H2) [34]. These genes contain cis-acting elements responsive to circadian control, abscisic acid, auxin, anaerobic, MeJA, and drought (ZmVHA-H1), as well as elements responsive to low-temperature, auxin, anaerobic, MeJA, drought, and wound (ZmVHA-H2). Studies have shown that overexpression of ScVHA-H in Suaeda corniculata improves the tolerance in transgenic alfalfa to salt and saline-alkali stresses [39]. These results suggest that the vacuolar ATPase subunit H protein may respond to low Pi stress by regulating the pH values in maize.
Studies have shown that the bZIP transcription factor ABI5 is involved in drought stress tolerance in barley [40], strawberry [41], Arabidopsis [42,43], and other species of plants [44]. The main Pi response pathway (i.e., SIZ1-PHR1-miR399-PHO2) has been identified in Arabidopsis [45,46], while SIZ1 negatively regulates the ABA signaling through sumoylation of ABI5 [47]. As a small ubiquitinlike modifier E3 ligase in plants, AtSIZ1 is identified as the focal controller of Pi starvation-dependent responses [46]. Furthermore, PHOSPHATE1 (PHO1) plays important roles in Pi homeostasis in Arabidopsis [48], while binding of ABI5 to the PHO1 promoter causes the downregulated expression of PHO1, which stops the ABA-insensitive germination of the ABI5 mutant [49]. These results suggest that ABI5 is involved in the response to low Pi stress by participating in either the SIZ1 or ABA signaling processes.
Violaxanthin de-epoxidase (VDE) is a type of lipocalin, which is characterized by a conserved structural organization with an eight-strand b-barrel and often binding hydrophobic molecules [50]. Studies have shown that VDEs play a key role in the pathway of carotenoid biosynthesis, which is involved in response to high light in bamboo [51], protecting the photosynthesis apparatus from damage caused by excessive light and chilling stress in tomato and Arabidopsis [52,53], and improving tolerance to drought and salt stress in transgenic Arabidopsis [54]. Furthermore, studies showed that VDE1 in maize and teosinte share the conserved functions of VDE1 in other plants [55]. These results suggest that the VDE is probably involved in the response to low Pi stress in maize by participating in the pathway of photosynthesis.

Conclusions
We apply the cDNA-SRAP method to identify the DEFs in both roots and leaves of maize under low Pi stress. In comparison to other transcriptomic analyses, the cDNA-SRAP technique shows substantial advantages of being easy and effective with low cost in the analysis of differential gene expression in maize. A total of 13 classes of differential gene expression patterns are revealed with a total of 2494 DEFs (1202 in leaves and 1292 in roots) identified in response to Pi starvation in maize. Results of sequencing and functional analyses demonstrate that 63 DEFs are involved in diverse groups of biological pathways (i.e., metabolism, photosynthesis, signal transduction, transcription, transport, cellular processes, genetic information, and organismal system), suggesting evidently that maize undergoes a complex adaptive process in response to low Pi stress. We speculate that the eight candidate genes identified involved in response to low Pi stress could be used in the genetic modification of maize.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.