Analysis of the common bean (Phaseolus vulgaris L.) transcriptome regarding efficiency of phosphorus use

Common bean is a highly important food in tropical regions, where most production occurs on small farms with limited use of technology and, consequently, greater vulnerability to abiotic stresses such as nutritional stress. Usually phosphorus (P) is the most limiting nutrient for crop growth in these regions. The aim of this study was to characterize the gene expression profiles of the genotypes of common bean IAC Imperador (P-responsive) and DOR 364 (P-unresponsive) under different P concentrations using RNA-seq transcriptome sequencing technology. Plants were grown hydroponically, with application of two P concentrations (4.00 mg L-1 restrictive level and 8.00 mg L-1 control level). Differential expression analyses, annotation, and functional classification were performed comparing genotypes within each P rate administered and comparing each genotype response to the different P levels. Considering differential expression analyses within genotypes, IAC Imperador exhibited 1538 up-regulated genes under P restriction and 1679 up-regulated genes in the control, while DOR 364 exhibited 13 up-regulated genes in the control and only 2 up-regulated genes under P restriction, strongly corroborating P-unresponsiveness of this genotype. Genes related to phosphorus restriction were identified among the differentially expressed genes, including transcription factors such as WRKY, ERF, and MYB families, phosphatase related genes such as pyrophosphatase, acid phosphatase, and purple acid phosphatase, and phosphate transporters. The enrichment test for the P restriction treatment showed 123 enriched gene ontologies (GO) for IAC Imperador, while DOR 364 enriched only 24. Also, the enriched GO correlated with P metabolism, compound metabolic processes containing phosphate, nucleoside phosphate binding, phosphorylation, and also response to stresses. Thus, this study proved to be informative to phosphorus limitation in common bean showing global changes at transcript level.


Introduction
Most common bean (Phaseolus vulgaris L.) production in developing countries occurs on small farms with limited technology, resulting in vulnerability to pest attack and abiotic stress, including water deficit and low soil fertility. Thus, development of new, higher-yielding cultivars bred for resistance to biotic and abiotic stresses is the main objective of plant breeding programs around the world [1] Crop yield depends on the nutrients available in the soil solution and on nutrient uptake by the plant. Different nutrients are structurally or functionally involved in the plant life cycle and are essential for plant growth and reproduction [2] Phosphorus (P) is considered the most limiting nutrient for growth of legume crops in tropical and subtropical regions and, since it is a non-renewable resource, successive crops result in continual degradation of soil P in the absence of additional fertilization [3] (Due to the interaction of inorganic phosphorus (Pi) with other elements, up to 80% of Pi applied can be fixed in the soil, forcing farmers to apply up to four times more than necessary for production.
Thus, increasing the efficiency of P use by crops has been a challenge, largely due to the complex interactions among the range of acquisition mechanisms and their effectiveness in different environments. MOURICE & TRYPHONE (2012) [4] also emphasize the importance of the use of bean varieties that are able to acquire phosphorus in environments with limiting soils, and this ability is attributed to genetic variability. Therefore, the study of gene relationships that determine greater efficiency of specific genotypes to acquisition and use of the nutrient is important.
According to GEORGE & RICHARDSON (2008) [5], plants have a wide range of physiological mechanisms and characteristics that facilitate an increase in the availability and acquisition of P from the soil. It is known that P restriction rapidly reduces rates of plant growth. Phosphorus deficiency alters cell biochemistry, biomass allocation, and root morphology to meet plant P demands. Typical responses of plants to P restriction include remobilization, reduction, or substitution of P in non-essential cell compounds, exudation of metabolites and enzymes to the rhizosphere, and changes in root morphology. There may be associations with microorganisms to increase the efficiency of P acquisition from the soil [6]Nevertheless, P demand in plant tissues and responses in regard to availability of the element vary among genotypes. SILVA et al. (2014) [7] assessed 20 common bean genotype responses to P nutritional deficiency in hydroponics through evaluation of morphophysiological and agronomic traits such as leaf chlorophyll content, shoot and root dry matter, leaf area, root area, root diameter, root length, and root volume, as well grain yield components Five genotypes stood out and were classified as efficient and responsive to P use: IAPAR 81, Carioca Comum, IAC Carioca Tybatã, IAC Imperador, and G 2333, whereas DOR 364 was one of the most inefficient and unresponsive genotypes. According to FAGERIA et al. (2015) [8], ER characterizes plants that have high grain yield at low P availability and are responsive to phosphate fertilization, whereas NENR indicates plants that have low yield at low P levels and low response to phosphate fertilization.
HERNÁ NDEZ et al. (2007) [9], studying the functional genomics of common bean plants under reduction in P supply, obtained the expression profile of the root system of the genotype Negro Jamapa 81 through macroarray analyses. They identified 126 differentially expressed genes with diverse functions in plant adaptation to P deficiency, contributing information on regulation and on signaling pathways during deficiency of the nutrient.
The aim of this study was to evaluate the agronomic traits and the gene expression profiles of the root system regarding P use efficiency in two contrasting common bean genotypes, IAC Imperador and DOR 364, grown hydroponically, and compare their response to application of a restrictive P concentration, 4.00 mg L -1 , and a control P concentration, 8.00 mg L -1 . RNAseq analysis was used to detect differences in gene expression between genotypes and P application rates. Correlations between the function of differentially regulated genes and the physiological traits of both bean genotypes are discussed in light of potential applications for P management in bean crops.

Analysis of agronomic traits
At 41 days after transplanting, the flowering stage (R6) was reached, and plants of both genotypes were collected for biometric analyses. All the shoot parameters (LA, LDM, SDM, PH, NNP) exhibited highly significant statistical differences (P>0.01) in relation to the P application rate factor ( Table 1). The greater availability of P provided by the control treatment resulted in an increase in all the parameters, except for plant height ( Table 2). In relation to the genotype factor, a statistical difference could be seen only for leaf dry matter (Table 1)-IAC Imperador exhibited 54.22% more leaf dry matter production than DOR 364, namely, mean values of 11.61 g for IAC Imperador and 6.29 g for DOR 364 (Table 2). In addition, a significant difference was also observed for the P application rate vs genotype interaction, likewise for the LDM variable (Table 1).
In relation to the root system traits, a significant difference was not observed for the P application rate factor. In relation to the genotype factor, significant differences (P>0.01) were observed (Table 1). IAC Imperador had 27.97% greater root length (RL), 27.55% greater root surface area (RSA), 26.64% greater root volume (RV), and 40.15% greater root dry matter (RDM) than the DOR 364 genotype ( Table 2).
In relation to yield components, analyses of variances showed statistical difference for the phosphorus application rate factor only for the GY. For the genotype factor, highly significant differences (P>0.01) could be seen for all the traits evaluated (NP, NS, 100SW, and GY), except for NSP. There was a significant effect for the P application rate vs genotype interaction only for the 100SW trait (Table 1).
Although the P application rate supplied did not result in statistical differences for the yield component traits (NP, NS, NSP, and 100SW), it resulted in a 41.37% greater increase for NP, 45.96% for NS, 8.48% for NSP, and 2.42% for 100SW in the control P rate than in the restricted P rate. The superiority of IAC Imperador was also clear as it had 62.8, 61.97, 2.34, 11.84, and 66.19% better performances than DOR 364 for the NP, NS, NSP, 100SW, and GY traits, respectively (Table 2).

Uptake efficiency, translocation, and P use indexes
The mean concentrations of P in the leaves, branches, roots, and grain in both treatments were assessed in order to calculate the indexes of nutrient uptake and use (Table 3). Mean concentrations in each tissue were higher in the treatment that received the control application rate.
The P uptake and use efficiency indexes were calculated and subjected to analyses of variance, showing significant effects for all the traits in relation to the P application rate treatment. In relation to the genotype factor, only the P translocation efficiency index did not show a significant difference. The coefficients of experimental variation had satisfactory values, ranging from 26.52% for the uptake efficiency characteristic to 0.66% for the P translocation efficiency characteristic, which had better experimental effectiveness (Table 4).
In regard to the P uptake efficiency index (PUE), significant differences were observed between the P treatments applied, showing that the greater the supply of the nutrient, the greater its uptake; thus, the control treatment (index of 32.34) showed greater uptake of the element than the P restriction treatment (index of 11.66). DOR 364 showed greater efficiency in uptake of the element, with an index of 26.46, while IAC Imperador had an index of 17.54. In regard to the translocation efficiency index (PTE), there was a statistical difference at 1% probability only for the P application rate factor, with indexes of 0.98 and 0.96 for the control and restrictive P level, respectively. The genotypes showed very similar values for this index, and it was not possible to detect a difference in translocation between them. This index showed the lowest coefficient of environmental variation (0.66%), which shows high experimental effectiveness for this trait.
The efficiency index of P use in the shoots (PUES) showed greater use of the element in the restricted P application rate, which was 31.45% greater than in the control rate, a highly significant difference. In relation to the genotypes, IAC Imperador showed a better P use index for the shoots, which was 0.48; whereas DOR 364 exhibited an index of 0.39.
A similar result was observed for the total P use efficiency index (PUET); there was greater total use of the element by the plants that received the restricted application rate, 34.48% more than the control rate. IAC Imperador proved to be more efficient in the use of total P, with an index of 0.53, while DOR 364 exhibited an index of 0.43.
Analysis of variance in regard to the efficiency index for P use in the shoots for grain formation (PUEGP) showed statistical difference for P application rates and for genotype. The treatment that received the restricted P rate showed a greater PUEGP index, 1.55, whereas the control treatment had an index of 1.01. IAC Imperador, higher yielding in both treatments, had an index of 1.89, almost three times greater than the DOR 364 index, which was 0.66.
For the P harvest index (PHI), a highly significant difference was observed for the P rates applied, and the treatments with the restricted and control P rates had indexes of 9.59 and 5.01, respectively. This result implies greater accumulation of P in the grains of the stressed treatments, that is, in the treatments in which lower grain yield was observed. There was also a significant difference for the genotypes; DOR 364 had a higher P harvest index, 8.51, while the P harvest index of the more productive genotype, IAC Imperador, was 6.09.

RNA-seq data
Sequencing of the libraries generated approximately 450 million paired reads. After filtering the low quality reads, the total number dropped to approximately 400 million sequences, from which 207 million reads were from the libraries that received P restriction (Samples 1-6), and 193 million were obtained from the control treatment (Samples 7-12). The 12 libraries data set were published in a public repository (NCBI) in SRA accession: PRJNA498535 (https://www. ncbi.nlm.nih.gov/sra/PRJNA498535).
Approximately 396 million reads were mapped on the common bean genome used as a reference (Pvulgaris 218 v1.0), from the Phytozome v.9.1 database, with an average of 99.2% of mapped reads. From a total of 373 million reads that were mapped in only one region, 193 million were observed in the treatment with P restriction and 180 million in the control condition.

Differential expression analysis
Comparison1. IAC Imperador vs. DOR 364 at the restrictive P level. Considering the six libraries that received the P restriction treatment, 30,060 expressed genes were identified, of which 27,191 are annotated according to the common bean unigenes. After FDR correction through the Benjamini-Hochberg procedure and baseMean filtering, 24,632 genes with unique sequences were identified, of which 4,123 showed significant differential expression levels (16.72%). In this condition, IAC Imperador showed 2,159 overexpressed genes, whereas DOR 364 had 1,964 (Fig 1).
Among the genes differentially expressed under P restriction, 48 transcription factor related genes were up-regulated in IAC Imperador and another 53 were up-regulated in DOR 364 (Table 5). Notably, the majority of these transcriptional regulator genes were identified as the WRKY family (12 up-regulated in IAC Imperador and 3 up-regulated in DOR 364). Ten upregulated ethylene-responsive transcription factors were also found in IAC Imperador and four in DOR 364 (Table 5). Another five up-regulated transcription factors were identified from the MYB family, one in IAC Imperador and 4 in DOR 364 (Table 5).
In our analysis, 53 up-regulated phosphatase related genes were found; 32 of them were over expressed in IAC Imperador and 19 in DOR 364. Among these genes, four of them code for pyrophosphatase; six of them code for acid phosphatase, which is known to catalyze the cleavage of inorganic phosphate from a broad array of phosphomonoesters and anhydrides; and two of them code for purple acid phosphatase (Table 6).
Six phosphate transporters genes were over expressed, three in IAC Imperador and three in DOR 364. Among them, two genes code for MSF transporters (Phvul.007G127400, Phvul.006G064900), the first up-expressed in IAC Imperador and the second in DOR 364. One up-expressed phosphate ABC transporter ATP-binding protein was found in IAC Imperador (Phvul.002G330700), which is responsible for coupling the energy of ATP hydrolysis to the import of phosphate across cellular membranes.
Finally, three genes with the SXP and EXS family domain were identified, Phvul.002G169700, Phvul.003G024600, and Phvul.007G275300, the first one up-expressed in IAC Imperador and the others in DOR 364.
In the enrichment test, 147 enriched gene ontologies (GO) were identified. From them, 123 were enriched for the IAC Imperador genotype under P restriction and 24 enriched for DOR 364, also under P restriction.
The genotype DOR 364 under P restriction exhibited 24 GOs terms, one classified as a cell component (CC), 15 as molecular functions (MF), and eight as biological processes (BP). The    cell component subcategory exhibited nine gene sequences attributed to the codifying domains of the protein histidine kinase complex. The molecular function category exhibited 15 subcategories of gene ontology, with the most frequent codifying domains related to oxidoreductase activities (28%), ADP binding (13%), oxidoreductase activity acting on paired donors with incorporation or reduction of molecular oxygen and binding of iron ions (9%), iron ion binding (9%), tetrapyrrole binding (8%), heme binding (8%), monooxygenase activity (7%), molecular transducer activity (6%), and signal transducer activity (5%). In regard to the biological processes category, the unigenes were distributed in eight functional categories, in which the domains most frequently codify for the oxidation-reduction process (51%), defense response (34%), and signal transduction by phosphorylation system (7%). Inside the BP category, our analysis showed the participation of 324 phosphorus metabolic process related genes. These genes are involved in chemical reactions and pathways involving P or P compounds, usually in the form of a phosphate group (PO4), anion, or salt of some phosphoric acid. All these genes were up-expressed in the IAC Imperador genotype under P restriction. Among these genes related to the P metabolic process, nine sequences were observed for the first time without homology with another gene already annotated.
Comparison 2. IAC Imperador vs. DOR 364 at the control P level. Considering libraries 7-12, IAC Imperador and DOR 364 in the control P level, 2161 differentially expressed genes were identified, in which 1069 genes were overexpressed in IAC Imperador and 1093 overexpressed in DOR 364 (Fig 2). In addition, 533 new genes were identified. It is also noteworthy that the number of differentially expressed genes is approximately half of the genes found in P restriction.
At the control P level, 49 transcription factors were identified, 21 being up-regulated in IAC Imperador and 27 in DOR 364. Ten up-regulated WRKY transcription factor families were found, eight in IAC Imperador and two in DOR 364. Five ethylene-responsive transcription factors were up-regulated in IAC Imperador, and two MYB transcriptional regulators were up-expressed in DOR 364. The response of the genotypes under the control P level in the upregulation of transcription factors was about half of that expressed under P restriction.
Twenty-three up-expressed phosphatase genes were found, 11 in IAC Imperador and 12 in DOR 364. Among them, only three acid phosphatase up-regulated genes were found in DOR 364.
Only one up-regulated phosphate transporter gene was found in DOR 364 (Phvul.006G064900), which was identified as a sodium-dependent phosphate transporter.
Thirty-seven terms of gene ontology were assigned to the DOR 364 genotype under the control, and the enrichment test performed showed roles related to cell components, molecular functions, and biological processes, with higher participation of genes related to binding (11.54%), heterocyclic compound binding (8.84%), organic cyclic compound binding (8.84%), nucleoside phosphate binding (4.93%), and nucleotide binding (4.93%).
Comparison 3. IAC Imperador at the restrictive P level vs. IAC Imperador at the control P level. The libraries of IAC Imperador under P restriction (Libraries 1-3) and IAC Imperador under control (Libraries 7-9) treatments were compared, with identification of 3,217 differentially expressed genes, of which 1,537 were overexpressed in the nutrient restriction condition and 1680 under the control (Fig 3).
Among the 1538 genes differentially expressed under P restriction, the 25 most up-regulated and the 25 most down-regulated genes ( Table 7) of greatest statistical significance were chosen as candidate genes in response to P deficit. Relative expression of the 25 most induced genes of the IAC Imperador genotype was from 10.6 to 33.2 times higher than under the control condition. In relation to the top-25 under-regulated genes, relative expression was from 9.9 to 44.5 times lower. In addition, 13 differentially expressed genes were not aligned, revealing new transcripts sequences related to stress from the P nutrient that can potentially be used for greater understanding of the functional responses of the genotype to P deficiency.
Comparison 4. DOR 364 at the restrictive P level vs. DOR 364 at the control P level. The libraries of DOR 364 under P restriction (Libraries 4-6) and under the control (Libraries 10-12) were compared, and, in contrast with the efficient genotype IAC Imperador, which exhibited 3217 differentially expressed genes, DOR 364 had only 15; two of them were up-regulated under P restriction and 13 up-regulated in the control ( Table 8). The profile of genotype expression is shown in Fig 4. Relative expression of the two up-regulated genes observed in the DOR 364 genotype was 4.8 and 3.7 times greater under the P restriction condition than under the control condition, similar to relative expressions of the thirteen repressed genes (Table 8).  Table 7. Selection of the 25 candidate unigenes most up-regulated and the 25 candidate unigenes most down-regulated of IAC Imperador (decreasing order of log2-FoldChange). From left to right: gene id (Phytozome); annotation of the gene based on blastx search using the database of P. vulgaris (Phytozome) with a cutoff of 1e-10; regulation; relative expression comparing restrictive P vs. control P level; baseMean; log2FoldChange; lfcSE (standard error); padj.

Discussions
Analyses of variance show that the level of P applied influences development of all the shoot traits and seed yield. In addition, the root system traits did not show significant differences in accordance with P level. Thus, it can be inferred that even at low P concentrations, roots had the same level of development as the control rate. A common response of plants to P deficiency is an increase in the size of the root system in relation to the shoots. The rate of formation of the shoots normally decreases with an increase in root growth, thus increasing the root/shoot ratio [10]. In this experiment, it was observed that plants grown under phosphorus restriction had a root/shoot ratio of 0.186, compared to 0.094 of plants grown under the control condition. IAC Imperador, the efficient and responsive genotype, had a ratio of 0.160, greater than the ratio of DOR 364, at 0.121.
In regard to the uptake, translocation, and P use indexes, DOR 364 had greater P concentrations in all plant tissues in comparison to IAC Imperador. In spite of exhibiting lower concentrations of P in its tissues, IAC Imperador had greater dry matter production, resulting in greater content of the element in the plant but consequent dilution of the nutrient. Considering the translocation efficiency index, although DOR 364 had a greater uptake index, IAC Imperador proved to translocate the element in the same manner, even with a smaller amount of P available in the plant tissue. In addition, IAC Imperador had a better efficiency index of P use in the shoots, total P use, and P use in the shoots for grain formation.
This greater efficiency and responsiveness of IAC Imperator can probably be explained by greater activation of the stress response genes and by greater dry matter production in relation to lower P concentration in the shoots under stress conditions. According to FAGERIA (1992)  [11], the values of P use generally increase with a decrease in nutrient levels, and P use efficiency is maximized at the lowest nutrient level and minimized at the highest nutrient level. According to SCHEIBLE AND ROJAS-TRIANA (2015) [12] plants respond to P limitation by changes in gene transcript expression, including genes involved in Pi transport, breakdown of P-containing molecules, photosynthesis, respiration, and amino acid and lipid metabolism. There are also changes in regulatory genes involved in signaling events, genes encoding transcriptional regulators (TRs) and microRNAs (miRs), components of the ubiquitin-proteasome system (UPS) for targeted protein degradation, and genes with functions that are yet unknown. As P-starvation develops gradually upon P withdrawal, P-starvation-responsive genes can be classified into early and late responsive groups. Early-responsive genes include many related to signal transduction events, general stress-related genes, genes encoding Pi transporters, SPXdomain proteins, purple acid phosphatases, ribonucleases, and also some cell wall-related genes. Late-responsive genes tend to be related to primary carbon metabolism, secondary metabolism, photosynthesis, protein synthesis, and hormone synthesis/signaling.
According to HIZ et al. (2014) [13], recent approaches to studying the profile of large-scale gene expression have proven to be an important tool for understanding how plants respond to biotic and abiotic stresses. The authors report that in spite of accumulation of information on transcriptome sequencing data of model plants and of agronomical important crops, few studies analyze the transcriptome under abiotic stress conditions in these species.
Under differential expression analysis, we found 4,123 differentially expressed genes under P restriction in comparison 1, almost double the 2,161 differentially expressed genes observed  (increasing order of log2FoldChange). From left to right: gene id (Phytozome); annotation of the gene based on blastx search using the database of P. vulgaris (Phytozome) with a cutoff of 1e-10; regulation; relative expression comparing restrictive P vs. control P level; baseMean; log2FoldChange; lfcSE (standard error); padj. Gene expression regarding phosphorus use efficiency under the control P level in comparison 2. Regarding comparisons 3 and 4, within genotypes (restrictive P vs. control P), IAC Imperador, the P-responsive genotype, exhibited 3217 differentially expressed genes, whereas the P-unresponsive genotype, DOR 364, had only 15. Among the differentially expressed genes, some were identified above with P-starvation as transcription factors, phosphatase related genes, phosphate transporters, and sequences that were not aligned, i.e., those that revealed new transcripts related to P stress that can potentially be used for greater understanding of functional responses to P-deprivation. Regarding the P restriction treatment, genes encoding proteins related to phosphorus restriction were identified, including P metabolism, phosphate-containing compound metabolic processes, nucleoside phosphate binding, phosphorylation, and response to stresses.  [9] to evaluate gene expression from P-deficient roots of bean plants compared to control P-sufficient roots. They found 126 differentially expressed cDNAs, which were then compared by the BLASTX tool to the Uniprot protein database to assign putative function. Functional categories were as follows: 23% of genes were regulation/signal transduction; 23% were those coding for genes that participate in secondary metabolism pathways and/or were related to several stress/defense plant responses; 13% were classified as membrane proteins or proteins that participate in transport; 8% were classified in cell structure, cell cycle, or developmental functions; 24% were classified in different metabolic pathways, such as Pi cycling, C and N metabolism, amino acid/protein synthesis or degradation, and lipid metabolism; and 9% had no known function.

N Unigene
O'Rourke et al. (2012) [14] studied the transcriptome of White Lupin in regard to phosphorus deficiency and observed 904 differentially expressed genes in the roots, of which 396 were up-regulated in the P-sufficient treatment and 535 were up-regulated in the P-deficient treatment, i.e., greater up-regulation under P-restriction. In addition, the authors found 23 transcription factor families, with bHLH and AP2_EREB as the two largest families responding to the P deficiency stress identified.
OONO et al. (2013) [15] analyzed the differential expression profile of four rice genotypes with variation in growth response to P starvation, as indicated by the shoot/root dry weight ratio, also through RNA-seq. They identified differences in expression of transcripts under nutrient deficiency, observing 3,080 differentially expressed transcripts, in which most of the up-regulated transcripts were more strongly expressed in the tolerant genotypes, showing specificity in genotype response to P restriction. These differentially expressed transcripts include the NAC transcription factor, and it has been reported that NAC family mediated auxin and ethylene signaling promote lateral root development.
Considering the transcription factors, 101 differentially expressed genes were identified in the restrictive P level and 48 in the control P level. Among these genes, greater participation of the WRKY, ethylene-responsive element binding (ERF), and MYB families was observed.
According to SCHEIBLE AND ROJAS-TRIANA (2015) [12], reduction in the WRKY transcript halts gene expression responses, causing reduction in Pi uptake, early anthocyanin production, and increase in lateral roots and root-hair formation, regardless of P-status. WRKY regulators are also involved in negative regulation of phosphate by binding to the PHO1 promoter and by activating the Pi transporter in response to P-starvation.
We found 15 WRKY differentially expressed genes under P-starvation; 12 of them were upregulated in IAC Imperador, i.e., 12 down-regulated and 3 up-regulated in DOR 364. According to CHEN et al. (2012) [16], several studies have indicated that WRKY transcription factors also participate in the nutrient deficiency response signaling pathways. AtWRKY75 was the first WRKY member reported to be involved in regulating phosphate starvation. It is strongly induced in the plant during Pi-deficiency and its suppression by RNAi silencing led to plants more susceptible to Pi stress, causing a decrease in Pi uptake and early accumulation of anthocyanin. In addition, the expression of several genes involved in Pi starvation responses, including phosphatases, Mt4/TPS1-like genes, and high-affinity Pi transporters, decreased when WRKY75 was suppressed.
The differential expression of 14 ethylene-responsive element binding factor (ERF) genes were found and 10 of them were up-expressed in IAC Imperador. The ERFs are known to be members of a novel family of transcription factors that are specific to plants [17]. ERF expression is up-regulated by various forms of biotic and abiotic stress, resulting in ethylene biosynthesis, and it is believed that ethylene might be a mediator of some stress responses.
Another type of differentially expressed transcription factor observed was MYB. The MYB transcription factor is arguably the most influential transcriptional regulator in the P-starvation response described in plants [12]. MYB also participates in metabolism, response to stress, and cell wall synthesis [18] (DALA VIA et al. 2015).
Five MYB differentially expressed transcripts were found. Four of them, Phvul.003G222900, Phvul.004G057800, Phvul.007G139300, and Phvul.007G211800, had repressed in IAC Imperador compared to DOR 364 under P deficiency. RAMIREZ et al. (2013) [19] performed a comparative gene expression analysis through qRT-PCR of the regulatory genes that code for the MYB family TF: PvRHR1, Pv4, PvPHO2, and three isoforms of PvmiR399 from BAT477 and DOR364, classified as P-deficiency tolerant and P-deficiency sensitive common bean genotypes. According to the authors, PvRHR1 is a positive regulator of PvmiR399 and P-responsive genes, whereas Pv4 is a negative regulator of PvmiR399 by the target mimicry mechanism. PvPHO2 is a negative regulator of P-responsive proteins, and it is the target gene of PvmiR399. Their analysis, performed on roots, showed significant increases in the transcript levels of PvPHR1, Pv4, and PvmiR399 (isoforms a, b, and e) under P deprivation compared to the control in both genotypes. Furthermore, as expected, the transcript levels of the negative regulator PvPHO2 decreased 3.7-fold from BAT477 under P deficiency, though not in the P-sensitive genotype DOR364. However, DOR364 had a PvPHO2 transcript level that was 2.7 more in BAT477 under P deficiency.
According to GEORGE & RICHARDSON (2008) [5], the hydrolysis of organic P, a process necessary for Pi uptake by plant roots, is mediated by the action of phosphatase enzymes in the extracellular environment. This process of phosphatase activity is induced under P-deficiency conditions. LAN et al. (2015) also point out that acquisition of Pi is supported by the up-regulation of several intracellular and secreted Pi-releasing enzymes such as, pyrophosphorylases and other phosphatases of various types within the Pi starvation-inducible (PSI) core genes. In our analysis, 53 differentially expressed genes were found at the restrictive P level and 23 at the control P level, of which nine code for acid phosphatase, four for pyrophosphatase, and two for purple acid phosphatase.
In genome-wide associations (GWAs), ZHANG et al (2014) [20] identified a candidate gene GmACP1 that encoded an acid phosphatase. Its over-expression in soybean hairy roots increased P efficiency by 11-20% compared to the control. A candidate-gene association analysis indicated that six natural GmACP1 polymorphisms explained 33% of the phenotypic variation. The authors also conclude that the favorable alleles and haplotypes of GmACP1 associated with increased transcript expression correlated with higher enzyme activity.
According to HERNÁ NDEZ-DOMÍGUEZ (2012) [21], there are still a number of cellular events poorly characterized under low Pi conditions. This is the case of the processes involving pyrophosphate (PPi), a byproduct of the cellular metabolism of biosynthesis of nucleic acids, carbohydrates, proteins, and fatty acids. It is known that the soluble inorganic pyrophosphatases (PPase) recycle pyrophosphate, and may play a role in plant adaptation to phosphorus deficiency. Considering that, the authors measured PPase mRNA expression in leaves, stems, and roots of bean plants by qRT-PCR; their results reveal changes in the expression and activity of PPases under long-term Pi starvation, suggesting a possible role for PPi during plant adaptation to this condition. Such variations in expression and activity indicate the existence of individual regulatory mechanisms for each PPase isoform.
LIANG et al. (2012) [22] observed that throughout the period of P deficiency, the growth of both genotypes studied, G19833 (P-efficient genotype) and DOR364 (P-inefficient genotype), was inhibited after Pi starvation; whereas the internal acid phosphatase activities of both genotypes increased. The authors emphasize that the group of purple acid phosphatases (PvPAPs) plays a vital role in plant adaptive strategies to phosphorus deficiency, and it is known that PvPAPs can be classified into two groups, including a small group with low molecular weight, and a large group with high molecular weight. Analyses of qPCR expression showed that the large group, PvPAP1, was constitutively expressed in the roots, whereas PvPAP2 was specifically expressed in roots, but PvPAP2 expression levels in G19833 were close to those in DOR364. It is noteworthy that the transcripts of small PvPAPs dramatically increased during Pi starvation in both genotypes, and G19833 showed significantly more than those in DOR364.
Phosphate acquisition and distribution across plant tissues is also dependent on an array of transporters, which include proton-phosphate cotransporters belonging to the family of PHT proteins. Seven differentially expressed phosphate transporters were found; six were up-regulated at the restrictive P level and one at the control level. Among these transporters, two code for MSF transporters, one for the phosphate ABC transporter ATP-binding protein, three for SXP with EXS family domain, and one at the control P level for sodium-dependent phosphate transporters.
The major facilitator superfamily (MFS) transporter is known to facilitate transport across cytoplasmic or internal membranes from a variety of substrates, including ions, sugar phosphates, nucleosides, amino acids, and peptides. The SPX domain is found at the amino terminus of a variety of proteins. The PHO1 gene family conserved in plants is involved in a variety of processes, most notably the transport of inorganic phosphate from the root to the shoot of the plant and mediation of the response to low levels of inorganic phosphate, since the EXS family is involved in signaling transduction mechanisms [23].
According to GLASS and GIRVAN (2014) [24], enrichment tests for evaluating the functional properties of gene sets are a routine step in understanding high-throughput biological data and are used to verify that the genes implicated in a biological experiment are functionally relevant and to discover functions between those genes. Gene ontology (GO) is one of the most widely used functional enrichment tools in annotation of genes in different species.
In our enrichment test, 147 enriched GOs were identified under P restriction, while 71 GOs were found at the control P level. The GO terms associated with the unigenes were categorized in three main functional categories: cell components, molecular functions, and biological processes.
In relation to the treatment with P restriction, 21725 GO terms were associated and distributed at 0.04%, 59.97%, and 39.98% in the functional categories of cell components, molecular functions, and biological processes, respectively. In the control treatment, 7876 GO terms were associated; the molecular function category had the highest representation, at 89.02%, followed by the biological process category, at 10.74%, and, finally, the cell component category with a small share, at 0.22%. Under P restriction, the attribution of GO terms was 63.75% greater than under the control, with a bigger effective participation of the molecular function category in response to stress. PATEL et al. (2014) [25] evaluated 20.4 million common bean unique reads from the NCBI database by the RNA-seq technique. Among these reads, 6,999 were mapped mapped on the common bean genome, and 1,679 known genes were identified, of which only 629 unigenes were annotated. In functional characterization, 3,724 terms of GO were attributed to the 629 annotated genes, which were distributed in the following manner: 46.37% molecular functions, 31.36% biological processes, and 22.26% cell components.
Differential expression analyses and enrichment test results corroborate the analyses and results observed regarding the efficiency indexes for P use, root/shoot ratio, grain yield, and genotype p-responsiveness under P restriction.

Conclusions
Sequencing and analysis of the transcripts by the RNA-seq technique revealed differences in the gene expression profile between genotypes in both treatments, and confirmed the superiority of the IAC Imperador genotype in relation to DOR 364. The candidate genes related to P deprivation responses were able to be identified as transcriptor factors, phosphate transporters, and phosphatases. In addition, enrichment analysis showed that IAC Imperador enriched more subcategories than DOR 364, including P metabolic processes.

Stress induction, analysis of agronomic traits, and estimation of efficiency in P uptake and use
The P-challenge experiment was conducted in a hydroponic greenhouse at the Instituto Agronômico de Campinas (Fazenda Santa Elisa, Campinas, SP, Brazil). The contrasting genotypes of common bean, IAC Imperador, previously classified as efficient and responsive, and DOR 364, as inefficient and non-responsive by SILVA et al. (2014) [7], were evaluated in regard to efficiency of P use.
A 2 × 2 factorial experimental design was applied, with the first factor constituted by two phosphorus application rates and the second factor constituted by the two genotypes. Nine replicates were used-three plants were collected in the full flowering stage (R6) to obtain gene expression profiles, three plants were collected to perform shoot and root system analysis, and three other replicates were grown until grain production.
The plants were grown in a closed hydroponic system in 3.5 L pots filled with medium particle size sand. The nutrient solution used in the two treatments was prepared with deionized water and was composed of 143.0 mg L -1 N, 132.5 mg L -1 K, 121.0 mg L -1 Ca, 25.5 mg L -1 Mg, 33.0 mg L -1 S, 1.81 mg L -1 Fe, 0.45 mg L -1 Cu, 0.18 mg L -1 Zn, 0.45 mg L -1 Mn, 0.45 mg L -1 B, 0.09 mg L -1 Mo, and 0.09 mg L -1 Ni. The treatments with P consisted of 4.00 and 8.00 mg L -1 P, constituting the restrictive P and control P levels, respectively. The plants were irrigated with a "half strength" solution, that is, with half the total concentration in the first two weeks after setting up the experiment. Electrical conductivity (1.8-2.0 mS cm -1 ) and pH of the nutrient solution (5.8) were checked daily and adjusted if necessary.
When the flowering stage (R6) was reached, three replicates were collected for evaluation of the following morphophysiological and agronomic traits: plant height (PH); number of nodes per plant (NNP); leaf area (LA); leaf, stem, and root dry matter (LDM, SDM, and RDM); root surface area (RSA), root diameter (RD), root length (RL), and root volume (RV). Chemical analyses of shoots, roots, and seeds were also performed. Three other replicates were collected to obtain the gene expression profile of the root system. Roots were softly washed with tap water to remove the sand and then immediately frozen in liquid nitrogen and stored in an ultra-freezer (-80˚C).
For evaluation of morphophysiological traits, the three replicates were separated into shoots and roots. Total leaf area (cm 2 ) was determined using an area meter (LI-COR, LI-3100C). These same plants were used for evaluation of shoot dry matter and were separated into leaves and stems and weighed on a precision balance. To obtain dry matter, the samples were dried in a forced air circulation laboratory oven at a temperature of 60˚C.
After washing the root system for removal of the substrate, it was evaluated in regard to total length (cm), surface area (cm 2 ), volume (cm 3 ), diameter (mm), and dry matter (g). Images of the roots of each plant were made in a scanner (LA2400-EPSON), and then root characteristics were calculated through use of the WinRHIZO software (Regent Instruments Inc., Quebec, Canada). Samples (1.0 g) of dry and ground plant tissue underwent nitric perchloric digestion to determine P concentration in the shoots, roots, and grains, which was determined by atomic absorption spectrophotometry [26].
The yield components evaluated were number of pods per plant (NPP), number of seeds per pod (NSP), 100 seed weight (100SW) in grams, and grain yield (GY) in g plant -1 .
In accordance with the methodology proposed by SIDDIQI & GLASS (1981) [27] and modified by MOURA et al. (1999) [28], the following phosphorous-related indexes were estimated: P uptake efficiency (PUE), P translocation efficiency (PTE), P use efficiency in shoots (PUES), total P use efficiency (PUET), P use efficiency in shoots for grain dry matter production (PUEGP), and P harvest index (PHI).
The results of each variable were subjected to ANOVA and the F test in a 2 × 2 factorial design, and the mean values were compared by the Tukey test at 5% probability for quantification of the effects of the treatments.

RNA-seq library preparation
The root systems of 12 samples, three biological replicates each of IAC Imperador and DOR 364 in the restrictive P and control P levels, were macerated, and approximately 100 mg were used for extraction of the total RNAs by the TRIzol Reagent extraction method [29]. The samples were sent to the Animal Biotechnology Laboratory of ESALQ/USP, where they were quantified in a Bioanalyser (Agilent) apparatus and adjusted to a concentration of 500 ng/μL. The libraries were prepared, using the TruSeq RNA Sample Prepv2 Low Throughput (LT) kit from Illumina according to manufacturer's instructions, and sequenced with the Hiseq 2500 apparatus from Illumina.

Pre-processing and data analysis
Demultiplexing was carried out using the CASAVA 1.8.2 (Illumina) software, which makes the base call of the raw data and transforms them into fastq format reads together with the phred quality scores. Reads quality was visualized using the FastQC program (www.bioinformatics. bbsrc.ac.uk/projects/). Filtering of low quality reads, sequences of adaptors, and vectors was performed by the Seqyclean v1.9.7 program (https://bitbucket.org/izhbannikov/seqyclean), using the 24 QScore. Contaminants were removed using the Univec database (http://www.ncbi.nlm. nih.gov/VecScreen/UniVec.html) as a reference. After filtering, reads with less than 65pb were removed.

Mapping
The samples were mapped against the Phaseolus vulgaris genome (Pvulgaris 218 v1.0) on the Phytozome v.9.1 database. The Bowtie2 v2.1.0 [30] program was used for this step, selecting the 'very-sensitive-local' option. After that, mapping quality was checked for each sample using the Samtools v.0.1.18 package (FlagStat tool, LI et al., 2009). With the same program, mapping was organized (sort tool), and the file was used for obtaining isoforms through the Cufflinks 2.2.1 package [31]. In this step, assembly was obtained using the RABT option, along with the gtf annotation file available for P. vulgaris on Phytozome and the reference genome.

Differential expression analysis
The files of the transcripts obtained for each sample were then compared using the Cuffmerge tool of the Cufflinks 2.2.1 package [31], in which the transcripts and isoforms are grouped into only one .gtf file. This file was used as the annotation offered to the script of the HTSeq-count v.05.4.p2 for extraction of raw read counts (http://www-huber.embl.de/users/anders/HTSeq/ doc/index.html). After obtaining the counts, the groups were analyzed using the DESeq2 package [32] (LOVE et al., 2014) from R/Bioconductor [33]. Within the DESeq2 program, the data was normalized by library size. In order to avoid artifacts caused by low expression profiles and high expression variance, only transcripts that had an average baseMean > 5 and the mean greater than the standard variation were analyzed. After that, the normalized counts from transcripts that passed in filtering were analyzed using the generalized linear model (GLM). The Benjamini-Hochberg [34] (BENJAMINI & HOCHBERG, 1995) False Discovery Rate (FDR) correction for multiple tests was applied.
To better understand the response mechanisms among the P levels applied and the genotypes, four different comparisons of differential expression tests were performed. First, genotypes were compared within each P rate applied as follows: Comparison 1. IAC Imperador vs. DOR 364 at the restrictive P level; and Comparison 2. IAC Imperador vs. DOR 364 at the control P level. Second, each genotype response was compared in accordance with the different P levels applied as follows: Comparison 3. IAC Imperador at the restrictive P level vs. IAC Imperador at the control P level; and 4. DOR 364 at the restrictive P level vs DOR 364 at the control P level.
The transcripts obtained from the Cufflinks merged .gtf file were aligned against a selection from the non-redundant NCBI database (nr, date: Nov. 12, 2013) using the Viridiplantae based protein subgroup (tax id 33090, August 2014). The blastx program with the 1e-10 cutoff was used [35,36] (ALTSCHUL et al., 1990; CAMACHO et al., 2009). Functional annotation of the terms of gene ontology (GO) and its derivatives was carried out using the Blast2GO-b2g4pipe v2.5 program [37] (CONESA et al., 2005), with 20 hits of Blast for each contig. Enrichment analysis was performed by B2GO using the Fisher exact test and considering only the categories with FDR <0.05. The genes that passed through the filtering of the baseMean >5 of each analysis of differential expression were used as background.