Transcriptomic Analysis of Banana in Response to Phosphorus Starvation Stress

Bananas are an important part of the diets of millions of people around the globe. Low P absorption and use efficiency significantly restrict banana yields. To further explore the molecular mechanisms of P regulation in banana plants, we used RNA sequencing-based transcriptomic analysis for banana plants subjected to Pi deficit stress for 60 days. We detected 1900 significantly differentially expressed genes (DEGs) in aboveground plant parts and 7398 DEGs in root parts under low P stress. Gene ontology (GO) classification analysis showed that 156,291 GO terms belonging to molecular functions, 53,114 GO terms belonging to cellular components, and 228,544 GO terms belonging to biological processes were enriched in the aboveground and root components. A number of DEGs involved in energy metabolism-related processes, signal transduction, control of rhizosphere P activation, and Pi mobilization were found, which were confirmed by quantitative reverse-transcription Polymerase Chain Reaction (qRT-PCR) analysis. At the transcriptomic level, we detected 13 DEGs from different organs and with different functions in the time-course response to phosphorus deficiency stress. These DEGs may include some key genes that regulate the phosphorus network, increasing our understanding of the molecular mechanism of Pi homeostasis in banana. These findings will also help develop biotechnologies to create a variant of banana with more effective Pi absorption and utilization.


Introduction
Phosphorus is one of the essential macronutrients required for normal plant growth and development.However, low P availability is a global problem restricting yield.The weathering of minerals changes the solubility of P; as the relative abundance of Fe and Al increases, the solubility of P becomes controlled by Fe or Al phosphates [1].As a result, Pi deficiency in the soil limits P absorption by the plant.In order to keep up with the growing global population, the demand for food and energy has led to excessive consumption of Pi fertilizer.Yet, according to the Food and Agriculture Organization of the United Nations (FAO; http://faostat.fao.org/), up to almost 70% of large-scale applications of Pi fertilizer might not be absorbed by crops due to low use efficiencies.
Banana was one of the earliest crops cultivated in the history of human agriculture; it is the most widely grown tropical fruit, spread over 130 countries throughout the tropics and subtropical areas.Most edible bananas are either derived solely from Musa acuminate Colla are a hybrid between two wild diploid species, M. acuminate Colla and Musa balbisiana Colla [2].In many countries of Africa and Asia, bananas are a vital food source.More than 110 million tons of bananas were produced across the world in 2014 (http://faostat.fao.org/).This fruit is not only becoming more popular in industrialized countries but also facilitating the economies of developing countries.However, faced with a worsening low-phosphate soil situation, arable land for bananas is becoming scarcer.Plants have adapted various strategies involving multiple gene regulation networks to survive under low P stress [3].Therefore, understanding the detailed gene regulation for maintaining Pi homeostasis is also important for understanding P use efficiency.Furthermore, developing banana varieties with effective Pi uptake and Pi homeostasis mechanisms under P-limited conditions will be significant for reducing P fertilizer use and enhancing agricultural sustainability.

Plant Growth Conditions
Banana seedlings at a uniform four-leaf stage were taken from planting bags and their roots carefully washed clean.The seedlings were then transferred to 3 L pots filled with clean sand for recovery, cultivated in a greenhouse at 28 • C to 35 • C in the July to August at Hainan, and provided with clean water during this time.After 2 weeks of growth recovery, the seedlings were supplied with 1/4-strength Hoagland's nutrient solution [4] for 2 weeks and then supplied with half-strength Hoagland's nutrient solution containing either 136 mg/L PO 4 3− , or 0 mg/L PO 4 3− , designated WT (Plant under normal conditions) and −P treatments, respectively.Every banana seedling was supplied with about 120 mL of Hoagland's nutrient solution every 2 days, and treated for a total of 60 days.

Measurement of Plant Weight and Total P Concentration
Twenty seedlings from two treatments were selected to record total leaf and symptomatic leaf number and their roots and the aboveground parts were separated to measure fresh weight using an electronic balance.The aboveground and root components were baked in an oven at 105 • C for 30 min, then at 80 • C for 5 h, before determining dry weight.For the measurement of total P concentration in plant tissues, dry samples of each part were ground to a fine powder and screened with a 0.5 mm aperture sieve, and then measured using the molybdenum antimony anti-colorimetric method [5].

Sample Preparation and RNA-Seq
Banana seedling roots and aboveground parts were harvested at 60 days after the start of experimental treatments, and these samples were then frozen immediately using liquid nitrogen and stored at −80 • C before being sent to BGI-Tech (Shenzhen, China) for RNA-seq by Illumina HiSeq2000 (Illumina Inc., San Diego, CA, USA).Raw data were filtered by removing reads including adapter, poly-N > 10%, and low quality reads in order to obtain clean reads as a foundation for further analysis.Clean reads were aligned to the banana reference genome database (http://banana-genome-hub.southgreen.fr/)by BWA (Broadband Wireless Access) [6].FPKM (Fragments Per Kilobase of exon model per Million mapped fragments) was used to determine the expression level of each gene, which could be calculated as the fold expression of a gene expressed in different samples.Our raw data and processed data have been uploaded in the Gene Expression Omnibus (GSE115792) at the NCBI (National Center for Biotechnology Information).

Identification of Differentially Expressed Genes (DEGs)
In our study, genes with a false discovery rate (FDR) ≤0.001 and an absolute value of Log 2 Ratio ≥1 were determined to be significantly DEGs that could be classified into up-or down-regulated DEGs (Tables S1 and S2).Gene ontology (GO) analysis is able to recognize the main biological functions that DEGs exercise.GO terms with a corrected p value ≤0.05 were defined as significantly enriched GO terms in significantly differential expression.

Quantitative Rreverse-Transcription PCR (qRT-PCR)
The total RNA was separately extracted from aboveground parts and roots in WT and −P treatments, respectively, using the CTAB (cetyltrimethylammonium bromide) method.The resulting RNA (1 µg) was treated with gDNA Eraser to remove contaminating genomic DNA and synthesized to the first-strand complementary DNA using a PrimeScript RT reagent kit with gDNA Eraser (TaKaRa, Shiga, Japan) according to the manufacturer's instructions.Primer pairs for quantitative RT-PCR were designed using Primer Premier 5.0 (PREMIER Biosoft International, Palo Alto, CA, USA).The cDNA reverse-transcription products were used as templates for qRT-PCR.An Applied Biosystems 7500 Real-Time PCR System (Life Technologies, Foster City, CA, USA) was used for quantitative RT-PCR, and the reaction consisted of 10 µL 2 × SYBR Premix Ex Taq II, 0.4 µL 50 × ROX Reference Dye or Dye II, 2 µL cDNA, 6 µL sterilized purified water, 0.8 µL forward primer (10 µM), and 0.8 µL reverse primer (10 µM) in a total volume of 20 µL.No-template controls were also prepared for each primer pair.The relative quantification method (Delta-Delta cycle threshold) was used to evaluate quantitative variation between the replicates examined.

Physiological Indicators of Banana Response to Deficient P Stress
At the onset of treatment, there were no phenotypic differences between plants under P deficit treatment or control treatment.Symptoms of the −P stress treatment were obvious after 48 days: Leaves became more blue, smaller, and curlier.Withered spots were dramatically increased in older leaves, and the whole plant was shorter compared with WT plants after 60 days (Figure 1A).There was no doubt that the banana seedlings exhibited reduced P usage efficiency and absorption of light due to phosphorus deprivation.We measured total leaf number, number of symptomatic leaves, fresh weight of aboveground parts, fresh weight of roots, dry weight of aboveground parts, and dry weight of roots in plants from the WT and −P groups, respectively, after 60 days (Figure 1B).The number of symptomatic leaves under low P stress was more than that in the control group.After 60 days, the shoot fresh weight under P deficiency was decreased by 30.57% compared with that of the WT group.This result was similar to that for root fresh weight.P deficiency led to an evident decrease in P concentration in roots and leaves (Figure 1C).After 60 days of treatment, the P concentration in the aboveground component was 6.33 mg/g and that in the roots was 1.78 mg/g, which were 50.66% and 65.83% lower, respectively, than the results for the WT.These observations suggest that banana growth will be restricted under phosphorus deficiency stress.

RNA Sequencing (RNA-Seq) and Alignment of Clean Reads to the Banana Reference Genome
In order to understand the underlying molecular events involved in the response to low P stress, we sent aboveground and root samples of the WT and −P treatment plants to BGI-Tech (Shenzhen, China) for RNA-seq.These RNA libraries were sequenced by Illumina HiSeq 2000 and yielded 24,243,199 and 24,073,614 clean reads in the WT and −P stress treatments, respectively, after filtering (Table 1).Clean reads were aligned to the banana reference genome database (http://banana-genome-hub.southgreen.fr/).The clean reads of the four samples mapped to the reference genome with up to 80% similarity.Mismatch and multi-position-match reads were excluded from further analysis, and unique-match reads were used for the subsequent analysis.

General Features of Differential Expression Profiling
The abundance of each gene was determined by counting fragments per kilobase of transcript per million reads mapped (FPKM) to infer the expression level.Therefore, this method could be used to directly compare the difference in gene expression among samples.The correlation of gene expression level among samples is a key criterion to test whether the experiments are reliable and whether the samples chosen are reasonable.If one sample is highly similar to another one, the correlation value between them is very close to 1.We calculated the correlation values between samples based on the FPKM results.According to the standard that the ENCODE (Encyclopedia of DNA Elements) project [7]

General Features of Differential Expression Profiling
The abundance of each gene was determined by counting fragments per kilobase of transcript per million reads mapped (FPKM) to infer the expression level.Therefore, this method could be used to directly compare the difference in gene expression among samples.The correlation of gene expression level among samples is a key criterion to test whether the experiments are reliable and whether the samples chosen are reasonable.If one sample is highly similar to another one, the correlation value between them is very close to 1.We calculated the correlation values between samples based on the FPKM results.According to the standard that the ENCODE (Encyclopedia of DNA Elements) project [7] recommends, the square of the correlation value should be ≥0.92(under ideal experimental conditions and with reasonable samples).The heatmap of correlations for these samples is shown in Figure 2. We can see that the correlation between WT root (CR) and −P root (PR) samples is not close indicating that the difference between CR and PR is significant.Put another way, there are a large number of differentially expressed genes (DEGs) between CR and PR.The same situation also appears for WT aboveground (CL) and −P aboveground (PL) samples (Figure 2).In our study, there were 26,215 and 28,921 genes expressed in the aboveground and root components, respectively, under P starvation conditions (Tables S3 and S4).Of these, 24,793 genes were commonly expressed in both aboveground samples (CL and PL), and a similar situation was detected in both root samples (CR and PR).After low P stress, a total of 9298 significantly DEGs were detected with the threshold of FDR ≤ 0.001 [8] and the absolute value of Log 2 Ratio ≥ 1 (Figure 3).There were 1064 up-regulated genes and 836 down-regulated genes in the aboveground samples (PL) under P deficit conditions.A total of 7398 DEGs were detected in the root parts under low P stress, of which 3290 genes were up-regulated and 4108 genes were down-regulated.The number of genes showing increased expression under P starvation stress was about three times that of the control, and the number of down-regulated genes was 4.9 times more than that of the control group.Gene ontology (GO) functional enrichment analysis was applied in our subsequent study in order to annotate the expression patterns of the DEGs to the given GO terms, which were classified as molecular function, cellular component, and biological process.In our RNA-seq report, 156 GO terms belonging to molecular function, 53 GO terms belonging to cellular component, and 228 GO terms belonging to biological process were enriched, respectively, in the aboveground component after P deficit treatment (Tables S5-S7).Our GO analysis identified that these genes were related to various functions.For example, Ma11_g17540 related to nucleotide binding was up-regulated, Ma03_g21400 (nucleic acid binding transcription factor activity) was down-regulated, Ma06_g03680 (catalytic activity) was up-regulated, and eight genes related to phosphatase activity (Ma06_g28160, Ma07_g02150, Ma02_g07260, Ma07_g18140, Ma04_g26760, Ma07_g22910, Ma07_g28770, and Ma04_g18100) were also up-regulated.Under the same stress, DEGs in the roots were enriched by 291,114, and 544 GO terms for molecular function, cellular component, and biological process respectively (Tables S8-S10), and these genes also revealed different functions during phosphorus starvation stress.
ideal experimental conditions and with reasonable samples).The heatmap of correlations for these samples is shown in Figure 2. We can see that the correlation between WT root (CR) and −P root (PR) samples is not close indicating that the difference between CR and PR is significant.Put another way, there are a large number of differentially expressed genes (DEGs) between CR and PR.The same situation also appears for WT aboveground (CL) and −P aboveground (PL) samples (Figure 2).In our study, there were 26,215 and 28,921 genes expressed in the aboveground and root components, respectively, under P starvation conditions (Tables S3 and S4).Of these, 24,793 genes were commonly expressed in both aboveground samples (CL and PL), and a similar situation was detected in both root samples (CR and PR).After low P stress, a total of 9298 significantly DEGs were detected with the threshold of FDR ≤ 0.001 [8] and the absolute value of Log2Ratio ≥ 1 (Figure 3).There were 1064 up-regulated genes and 836 down-regulated genes in the aboveground samples (PL) under P deficit conditions.A total of 7398 DEGs were detected in the root parts under low P stress, of which 3290 genes were up-regulated and 4108 genes were down-regulated.The number of genes showing increased expression under P starvation stress was about three times that of the control, and the number of down-regulated genes was 4.9 times more than that of the control group.Gene ontology (GO) functional enrichment analysis was applied in our subsequent study in order to annotate the expression patterns of the DEGs to the given GO terms, which were classified as molecular function, cellular component, and biological process.In our RNA-seq report, 156 GO terms belonging to molecular function, 53 GO terms belonging to cellular component, and 228 GO terms belonging to biological process were enriched, respectively, in the aboveground component after P deficit treatment (Tables S5-S7).Our GO analysis identified that these genes were related to various functions.For example, Ma11_g17540 related to nucleotide binding was up-regulated, Ma03_g21400 (nucleic acid binding transcription factor activity) was down-regulated, Ma06_g03680 (catalytic activity) was up-regulated, and eight genes related to phosphatase activity (Ma06_g28160, Ma07_g02150, Ma02_g07260, Ma07_g18140, Ma04_g26760, Ma07_g22910, Ma07_g28770, and Ma04_g18100) were also up-regulated.Under the same stress, DEGs in the roots were enriched by 291,114, and 544 GO terms for molecular function, cellular component, and biological process respectively (Tables S8-S10), and these genes also revealed different functions during phosphorus starvation stress.

Energy Metabolism-Related Genes during P Deficiency Stress
In our analysis, 38 genes involved in glycolysis were significantly expressed under P deprivation conditions in the aboveground samples (

Energy Metabolism-Related Genes during P Deficiency Stress
In our analysis, 38 genes involved in glycolysis were significantly expressed under P deprivation conditions in the aboveground samples (Table 2).Of these 38 DEGs, 23 genes were up-regulated, and the expression of the remaining genes was decreased.Glyceraldehyde-3-phosphate dehydrogenase, phosphoglycerate kinase, and pyruvate kinase are three key enzymes in glycolysis.Two glyceraldehyde-3-phosphate dehydrogenases (Ma11_g17540 and Ma06_g03680) were up-regulated; one phosphoglycerate kinase gene (Ma08_g33800) and a pyruvate kinase gene (Ma00_g00620) were also up-regulated.We identified 14 genes significantly expressed in the TCA (tricarboxylic acidcycle), of which 11 genes (Ma05_g07350, Ma10_g09450, Ma03_g09370, Ma02_g01890, Ma06_g27170, Ma06_g03480, Ma06_g18050, Ma05_g02120, Ma02_g23000, Ma05_g03680, and Ma02_g22920) were up-regulated and three genes (Ma04_g27480, Ma06_g33960, and Ma10_g29220) were down-regulated.Isocitrate dehydrogenase, which catalyzes the conversion of isocitric acid into alpha-ketoglutaric acid and reduces NAD + or NADP + , is one of the key enzymes in the TCA cycle.One isocitrate dehydrogenase gene (Ma02_g22920) increased in expression during Pi starvation and then augmented the NADH and NADPH content to maintain the antioxidant system, thus avoiding oxidative toxicity.Among the genes involved in glycolysis in the roots, three glyceraldehyde-3-phosphate dehydrogenase genes (Ma11_g08300, Ma11_g17540, and Ma08_g33830) were down-regulated during −P treatment; three phosphoglycerate kinase genes (Ma05_g00310, Ma11_g16870, and Ma06_g01900) and three pyruvate kinase genes (Ma03_g28290, Ma04_g28340 and Ma11_g09970) were also down-regulated.The expression of these genes encoding these three key enzymes in roots is the opposite of the expression in aboveground parts.For genes involved in the TCA cycle in roots, three isocitrate dehydrogenase genes (Ma04_g00870, Ma04_g24900 and Ma10_g11000) and two malate dehydrogenase genes (Ma05_g03680 and Ma04_g28520) were down-regulated under P deficiency.

Signal Transduction-Related Genes under Low P Conditions
Auxin is a vital signal molecule that acts during various abiotic stresses [9].In our study, five auxin influx carrier genes including Ma08_g02110, Ma11_g22520, Ma06_g35720, Ma06_g35630, and Ma05_g31240 were up-regulated under −P treatment.On the contrary, the expression of an auxin efflux gene (Ma01_g09070) was down-regulated.These results suggest that an increase in the auxin content of banana roots may promote the growth of the plant root and the absorption of phosphorus, which provides better resistance to P starvation stress.
The universal second messenger Ca 2+ plays an important role in resisting P deficiency stress.According to some reports, the intracellular Ca 2+ concentration increases rapidly during cold stress, followed by increases a number of signals mediated by a series of protein phosphorylation cascades [10].Interestingly, our study demonstrated that the expression of most Ca 2+ signal DEGs was significantly increased during P deficiency (Table 3).In aboveground parts, five calcium/calmodulin-dependent protein kinase genes were up-regulated by −P treatment.The expression of one Ca 2+ -transporting ATPase gene (Ma07_g27150), and nine calcium-binding protein genes was also increased.In the roots of banana, 10 calcium-dependent protein kinase genes were up-regulated during −P treatment.The expression of three calcium ion transport genes (Ma04_g05840, Ma09_g21740 and Ma05_g07830), one calcium exchanger 1-like gene (Ma06_g25970), and two calcium-binding protein genes (Ma04_g24480 and Ma03_g18040) was increased under P deficit stress.The altered expression of these DEGs indicated that the density of Ca 2+ accumulating between or inside cells and higher efficiency of Ca 2+ transportation could facilitate the response to phosphorus deprivation stress.

Control of Rhizosphere P Activation Genes
A large fraction of soil P is fixed into organic or inorganic (e.g., Ca-P, Fe-P, Al-P) forms, which cannot be directly utilized by plants [11].At this time, acid phosphatases (APases), RNases, protons (H + ) and organic acids will be excreted from the root to activate insoluble phosphorus [12].
An important component of the plant P-starvation response is the excretion of H + and organic acids by roots into the rhizosphere, which would activate the fixed Pi from soils [13].In our study, the expression of three H + -ATPase genes (Ma03_g29190, Ma03_g12710, and Ma09_g10590) was increased under −P treatment.Enhanced synthesis of organic acid by P-starved plants has been correlated with up-regulation of phosphoenolpyruvate carboxylase (PEPC) [14].Furthermore, overexpression of a maize PEPC gene increased organic acid synthesis and secretion in rice [15].During phosphorus deficiency stress, the expression of five PEPC genes (Ma04_g18750, Ma04_g08080, Ma02_g14160, Ma10_g25380, and Ma05_g31880) belonging to banana roots was distinctly increased.

Genes Involved in Pi Mobilization
To maintain the energy requirement balance in each part of a plant, P mobilization must operate efficiently.Pi translocation in different plant tissues is conducted through both PHT (phosphate transporter) and SPX domain (The SPX domain is named according to a homologous sequence shared by yeast SYGI, PHO81 and human XPRI) proteins [16].In our study, five SPX domain-containing protein genes were up-regulated (Table 4).
In our study, the expression of 12 phosphate transporter genes (PHTs) from the roots (Table 4) was significantly up-regulated.In addition, the expression of four phosphate transporter genes we selected from the aboveground plant components, including Ma03_g26260 (inorganic phosphate transporter 1-11-like), Ma02_g06130 (probable inorganic phosphate transporter 1-8), Ma01_g01890 (probable inorganic phosphate transporter 1-8) and Ma04_g36790 (probable inorganic phosphate transporter 1-8), was obviously up-regulated.Interestingly, the four up-regulated PHT genes were expressed in both aboveground and root regions.

SPX Domain-Containing Protein
Pi Transporters WRKY transcription factors are a large family of regulatory proteins in plants, which possess a WRKY domainand a type of zinc-finger motif [17].According to a previous study, a high percentage of WRKY transcription factors were responsive to biotic stress [18], but there is little information about roles of WRKY proteins in plant responses to abiotic stress in general and Pi starvation in particular.In this study, some WRKY transcription factor genes were found to be up-regulated under −Pi treatment (Table 5).The expression of one WRKY transcription factor 26 gene, three WRKY transcription factor 71 genes, and one WRKY transcription factor 31 gene was increased in banana roots under low phosphorus stress.In the aboveground compartments of banana, one WRKY transcription factor 50, one WRKY transcription factor 70, one WRKY transcription factor 26, one WRKY transcription factor 40, and three WRKY transcription factor 31 genes were up-regulated in response to phosphate deficit conditions.

Validation of the Selected Degs by Quantitative RT-PCR Analysis
We deem that the process of controlling rhizosphere P activation and P mobilization are more important than other courses in phosphorus deficit stress.One gene encoding purple acid phosphatase (PAP) and two genes encoding phosphoenolpyruvate carboxylase (PEPC) involved in rhizosphere P activation were selected to detect mRNA levels by quantitative RT-PCR.One WRKY transcription factor gene in each of the root and aboveground parts, respectively, was detected by quantitative RT-PCR.In addition, five and three phosphate transporter (PHT) genes, associated with the course of Pi mobilization, in the root and aboveground parts, respectively, were chosen for real-time quantitative PCR (Table 6).The actin gene was used as a reference gene.The quantitative RT-PCR results from the aboveground and root samples are shown in Figures 4 and 5.The expression profiles of all 13 detected genes revealed a similar trend and consistent results between the RT-PCR and the RNA-seq.Four aboveground DEGs displayed significant up-regulation in plants under Pi deficiency conditions compared with WT conditions, including Ma03_g26260 (inorganic phosphate transporter 1-11), Ma01_g01890 (inorganic phosphate transporter 1-6), Ma04_g36790 (probable inorganic phosphate transporter 1-4), and Ma08_g01650 (probable WRKY transcription factor 26). Interestingly, two genes exhibiting functional homology to inorganic phosphate transporter 1-6 in the root also showed significant enhancement in transcript levels under P deficiency.A similar situation was found for the gene encoding the inorganic phosphate transporter 1-11 protein of the root part, for which the mRNA abundance was substantially increased in the aboveground sample following −Pi treatment.These circumstances indicate that some PHT genes will be expressed in the whole plant to load and transport Pi faster in response to phosphorus deficit stress.A gene related to PAP and WRKY11 displayed a higher mRNA level during −P treatment than WT treatment, as determined by qRT-PCR.Strikingly, the expression of a PEPC gene was up-regulated; however, a PEPC kinase (PPCK) gene revealed no remarkable increase in its mRNA level.

Discussion
Bananas are vital commercial crops in Hainan, China.However, no previous research has examined phosphorus deficit stress.Low P absorption and low P use efficiency problems are becoming more and more severe, so it is necessary to study the connection between phosphate utilization and the key genes controlling the regulation network.After 48 days of P deficit stress, adverse symptoms in banana were more serious in terms of plant height, leaf color, and withered spots than under WT conditions.The phenotypic differences in the phosphorus response between the −P and WT treatments allowed us to explore differential gene expression in response to P stress.Energy metabolism consists of two processes: Glycolysis and the tricarboxylic acid (TCA) cycle.Glycolysis is the pathway that converts glucose into pyruvate in the cytoplasm.The TCA cycle is completed in the mitochondria.During the course of the process, the NADH, H + , and reduced flavin adenine dinucleotide (FADH2) generated will subsequently be oxidized in the electron transport chain on the mitochondrial inner membrane and then produce a great deal of ATP.Proverbially, ATP provides energy for all life activities.Accurately, for our study, phosphorus transporters are

Discussion
Bananas are vital commercial crops in Hainan, China.However, no previous research has examined phosphorus deficit stress.Low P absorption and low P use efficiency problems are becoming more and more severe, so it is necessary to study the connection between phosphate utilization and the key genes controlling the regulation network.After 48 days of P deficit stress, adverse symptoms in banana were more serious in terms of plant height, leaf color, and withered spots than under WT conditions.The phenotypic differences in the phosphorus response between the −P and WT treatments allowed us to explore differential gene expression in response to P stress.Energy metabolism consists of two processes: Glycolysis and the tricarboxylic acid (TCA) cycle.Glycolysis is the pathway that converts glucose into pyruvate in the cytoplasm.The TCA cycle is completed in the mitochondria.During the course of the process, the NADH, H + , and reduced flavin adenine dinucleotide (FADH 2 ) generated will subsequently be oxidized in the electron transport chain on the mitochondrial inner membrane and then produce a great deal of ATP.Proverbially, ATP provides energy for all life activities.Accurately, for our study, phosphorus transporters are energized by ATP-dependent proton efflux and actively assimilate P against a steep concentration gradient, as the soluble Pi concentration in the rhizosphere can be up to 10,000-fold lower than that in root cells [12].In our study, genes involved in glycolysis and the TCA cycle in banana roots were down-regulated, indicating that the two energy metabolism-related processes were restricted in roots during P deficiency.As a consequence of the large decline in cytoplasmic Pi levels that follows severe Pi stress, large (up to 80%) reductions in intracellular levels of ATP, ADP, and related nucleoside phosphates also occur [19].Alternative pathways of cytosolic glycolysis and mitochondrial electron transport may promote the survival of Pi-deprived plants [12].Plants also increase P use efficiency during Pi starvation via alteration of metabolic pathways.PEPC-malate dehydrogenase-malic enzyme glycolytic bypass, which reduces the requirement for Pi, is a well-documented alteration of glycolysis under Pi starvation [11].Up-regulation of PEPC and PPCK by P deficiency has been documented in Arabidopsis, rice, white lupin, orange, and melon [12,20,21].
In our study, we found some PAP (Purple acid phosphatases) genes, which play roles in intraand/or extracellular Pi scavenging and recycling during Pi starvation.PAPs comprise the largest class of plant APases and have been documented to function in Pi release from other organic P forms [22].Take for instance AtPAP10, AtPAP12, or AtPAP26, which release Pi from glycerol 3-phosphate, ADP, or DNA in Arabidopsis [23].The transcription factors PHR1, WRKY75 and ZAT6 have been implicated in the control of PAP expression in −Pi Arabidopsis, while other studies have revealed PAPs whose synthesis is controlled by posttranscriptional mechanisms [24,25].An important component of the plant PSR (phosphate starvation-response) is root excretion of organic acids into the rhizosphere.The mechanism responsible for the release of organic acids was unknown until plasma membrane H + -pumping ATPases were shown to be involved in plant adaptation to Pi starvation.A previous study showed that a tomato (Solanum lycopersicum Mill.) 14-3-3 member, TFT, functions as a root plasma membrane H + -ATPase that increases the release of H + under low P conditions [26].This has been correlated with the up-regulation of many H + -ATPase genes that encode proteins to transport anions such as citrate and malate from root cells into the rhizosphere [27].Organic acid excretion can mobilize Pi bound in humic-metal complexes, enhancing both the solubility of organic P and its amenability to dephosphorylation by secreted PAPs, and may also promote the growth of symbiotic rhizosphere microbes that facilitate root Pi acquisition [14,15,28,29].
Several proteins in yeast, including Pi sensors, sharing an SPX domain maintain Pi homeostasis [30] and some proteins bearing this domain are involved in Pi starvation signaling in plants [31].In recent years, studies have showed the SPX domain proteins in yeast and Arabidopsis could be involved in the fine tuning of Pi transport and signaling through mechanisms such as physical interactions with other proteins [32,33].Pho87, Pho90, and Pho91, part of the regulon, possessing an SPX domain, mediate Pi translocation under P-deficient conditions [16].Additionally, proteins consisting of SPX domain play a critical role in loading Pi into the xylem of roots [31].OsSPX-MFS1, identified as a Pi transporter functioning in leaf Pi reallocation, also contains an SPX domain, and the mRNA levels of OsSPX-MFS1 were induced by Pi deficiency [34].
P is transported into and within plants mainly through high-affinity Pi transporters encoded by Pi transporter genes, which is one strategy plants use to adapt to low P soil [35].Most of the high-affinity Pi transporter genes (PHT) identified in the Arabidopsis and rice genomes over the past decades were expressed either exclusively or predominately in the roots [36,37].In rice, the expression of OsPT2 and OsPT6 were showed to be dominant in roots but also increased in leaves in response to Pi starvation [38].In our analysis, some PHTs could be significantly expressed in both leaves and roots.According to phylogenetic analysis, sequence, and functional similarity, phosphorus transporters can be classified into the Pht1, Pht2, Pht3, Pho1, and Pho2 families.Most members of the Pht1 family involved in Pi acquisition from the soil solution and the redistribution of Pi within the plant are up-regulated under P-deficit stress, and our results were consistent with this tendency.Mutants with knockdown PHR1 revealed impaired Pi starvation responses, including reduced root-to-shoot ratio, decreased accumulation of anthocyanins, and low Pi uptake [39].There is no report on Pi transporters involved in maintaining Pi homeostasis in banana.However, in our study, we identified several Pht1 family members in banana roots.
The Arabidopsis phosphate starvation response regulator 1 (PHR1) is a MYB transcription factor, playing a key role in Pi starvation signaling by binding to a cis-element of GNATATNC, named the PHR1 binding sequences (P1BS) [40].Genes downstream of PHR1 include genes encoding the signal molecules SPXs [32,41], PHTs (phosphate transporters) [37], and PAPs [42,43].However, there is no report of the detection of PHR1 in banana.In our study, we found two PHR1-like genes Ma07_g15660 (protein PHR1-LIKE 1-like) and Ma10_g23180 (protein PHR1-LIKE 1-like), which may regulate SPXs and PHTs involved in this phosphorus deficit stress.
Under low phosphorus stress, banana plants activate a cascade of cellular events including a number of specific transcription factors, among of which those of the WRKY family regulate a multitude of abiotic stress responses.They are essential for modulating Pi homeostasis in Arabidopsis.WRKY42, regulating PHT1;1 and PHO1 expression under Pi-sufficient conditions, mediates Pi uptake and translocation inside the plant [44].WRKY75 acts as a positive regulator of Pi stress responses, and suppression of WRKY75 expression results in impaired Pi starvation responses in Arabidopsis [45].WRKY33 is also a member of the WRKY family.Through KEGG pathway analysis, we could detect a few WRKY33 genes in our study, including Ma06_g28830, Ma07_g26900, Ma02_g16960 and Ma06_g07330.The over-expression of WRKY33 produces effective responses to various microbial pathogens in Arabidopsis [46], and WRKY33 modulates the expression of salicylate and jasmonate [17].Interestingly, the expression of a series of WRKY33 genes was greatly increased under Pi deficiency stress, indicating that WRKY33 may induce the signaling molecules of resistance genes, especially in response to plant abiotic stress.

Conclusions
In summary, this study reveals comparative transcriptomic profiles under −P and control treatments as investigated by RNA-seq and quantitative RT-PCR.At the transcriptomic level, we detected several DEGs in different organs that carried out different functions in the time-course response to phosphorus deficiency stress.Those DEGs may interact with each other, and some key genes may regulate the phosphorus network.These findings will enhance our understanding of the molecular mechanism of Pi homeostasis in banana.These findings will enhance our understanding of the molecular mechanism of Pi homeostasis in banana and help develop biotechnologies to create a variant of banana with more effective Pi absorption and utilization through genetic modification or targeted gene mutation.Biotechnology, such as this is, needed to reduce P fertilizer use and sustain the banana industry.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2073-4395/8/8/141/s1,Table S1: DEGs of the aboveground, Table S2: DGEs of the root, Table S3: Expressed genes in the aboveground, Table S4: Expressed genes in the root, Table S5: GO terms belonging to molecular function, Table S6: GO terms belonging to cellular component, Table S7: GO terms belonging to biological process, Table S8: GO terms for molecular function, Table S9: GO terms for cellular component, Table S10: GO terms for biological process.

Figure 2 .
Figure 2. Correlations of samples.CL: Aboveground part of WT; CR: Root of WT; PL: Aboveground section of −P treatment; PR: Root of −P treatment.

Figure 2 . 17 Figure 3 .
Figure 2. Correlations of samples.CL: Aboveground part of WT; CR: Root of WT; PL: Aboveground section of −P treatment; PR: Root of −P treatment.Agronomy 2018, 8, x FOR PEER REVIEW 7 of 17

Figure 3 .
Figure 3. Differentially expressed genes (DEGs) in aboveground and root components during P starvation treatment.CL: Aboveground part of WT; CR: Root of WT; PL: Aboveground section of −P treatment; PR: Root of −P treatment.

Figure 4 .
Figure 4. Relative expression of P deficit-induced genes in aboveground plant parts in WT and −P treatments was determined by quantitative RT-PCR analyses.Expression was normalized to that of actin.The transcript levels from the WT sample were set to 1.The red columns represent the expression levels determined by qRT-PCR (left axis), while the lines represent the gene expression levels determined by RNA-seq (right axis).The transcript abundances of genes encoding inorganic phosphate transporter 1-11 protein (A); norganic phosphate transporter 1-6 (B); probable inorganic phosphate transporter 1-4 (C), and probable WRKY transcription factor 26 (D) were determined and compared across the time course of phosphorus deficit stress.Data represent means ± SD of three replicates (n = 3).Different lowercase letters above the columns indicate a significant difference at p ≤ 0.05 between the columns by T test using SPSS statistical software.Data in columns with the same letters showed no significant difference (p > 0.05).

Figure 4 .Figure 5 .
Figure 4. Relative expression of P deficit-induced genes in aboveground plant parts in WT and −P treatments was determined by quantitative RT-PCR analyses.Expression was normalized to that of actin.The transcript levels from the WT sample were set to 1.The red columns represent the expression levels determined by qRT-PCR (left axis), while the lines represent the gene expression levels determined by RNA-seq (right axis).The transcript abundances of genes encoding inorganic phosphate transporter 1-11 protein (A); norganic phosphate transporter 1-6 (B); probable inorganic phosphate transporter 1-4 (C), and probable WRKY transcription factor 26 (D) were determined and compared across the time course of phosphorus deficit stress.Data represent means ± SD of three replicates (n = 3).Different lowercase letters above the columns indicate a significant difference at p ≤ 0.05 between the columns by T test using SPSS statistical software.Data in columns with the same letters showed no significant difference (p > 0.05).

Figure 5 .
Figure 5. Relative expression of Pi deficit-induced genes in roots in WT and −Pi treatment was determined by quantitative RT-PCR analyses.Expression was normalized to that of actin.The transcript levels from the WT sample were set to 1.The red columns represent the expression levels determined by qRT-PCR (left axis), while the lines represent the gene expression levels determined by RNA-seq (right axis).The transcript abundances of genes encoding purple acid phosphatase 2 protein (A); hosphoenolpyruvate carboxylase (B); putative phosphoenolpyruvate carboxylase kinase 2 (C); inorganic phosphate transporter 1-11 (D); phosphate transporter PHO1-2 (E); putative WRKY11 (F), inorganic phosphate transporter 1-6 (G,H), and phosphate transporter PHO1 homolog 1 (I) were determined and compared across the time course of phosphorus deficit stress.Data represent means ± SD of three replicates (n = 3).Different lowercase letters above the columns indicate a significant difference at p ≤ 0.05 between the columns by T test using SPSS statistical software.Data in columns with the same letters showed no significant difference (p > 0.05).

Table 1 .
Statistical analysis of clean reads mapped to the reference genome.

Table 1 .
recommends, the square of the correlation value should be ≥0.92(under Statistical analysis of clean reads mapped to the reference genome.
CL: Aboveground part of WT; CR: Root of WT; PL: Aboveground section of −P treatment; PR: Root of −P treatment.

Table 2 .
Energy metabolism-related genes of aboveground samples.

Table 4 .
Up-regulated genes involved in Pi mobilization under −P treatment.