Phosphate Deficiency Negatively Affects Early Steps of the Symbiosis between Common Bean and Rhizobia

Phosphate (Pi) deficiency reduces nodule formation and development in different legume species including common bean. Despite significant progress in the understanding of the genetic responses underlying the adaptation of nodules to Pi deficiency, it is still unclear whether this nutritional deficiency interferes with the molecular dialogue between legumes and rhizobia. If so, what part of the molecular dialogue is impaired? In this study, we provide evidence demonstrating that Pi deficiency negatively affects critical early molecular and physiological responses that are required for a successful symbiosis between common bean and rhizobia. We demonstrated that the infection thread formation and the expression of PvNSP2, PvNIN, and PvFLOT2, which are genes controlling the nodulation process were significantly reduced in Pi-deficient common bean seedlings. In addition, whole-genome transcriptional analysis revealed that the expression of hormones-related genes is compromised in Pi-deficient seedlings inoculated with rhizobia. Moreover, we showed that regardless of the presence or absence of rhizobia, the expression of PvRIC1 and PvRIC2, two genes participating in the autoregulation of nodule numbers, was higher in Pi-deficient seedlings compared to control seedlings. The data presented in this study provides a mechanistic model to better understand how Pi deficiency impacts the early steps of the symbiosis between common bean and rhizobia.


Plant Material
Common bean (P. vulgaris L. cv Negro Jamapa) seeds were produced under greenhouse conditions in the Facultad de Estudios Superiores Iztacala, Universidad Nacional Autónoma de México (UNAM), at Tlalnepantla, Estado de Mexico, Mexico. Seeds were surfaced sterilized by soaking in 70% ethanol for 1 min, followed by treatment for 10 min with 10% bleach. Seeds were subsequently washed 10 times in sterile water. Sterilized seeds were germinated for three days at 25 • C in Petri dishes containing wet and sterile germination paper (Anchor paper CO, St. Paul, MN, USA). Petri dishes containing common bean seeds were kept under dark conditions at 25 • C.

Bacterial Strains and Culture Conditions
Rhizobium tropici CIAT 899 strain was used to inoculate common bean seedlings. R. tropici cells were grown for 2 days at 30 • C on PY medium (5 mg/L peptone; 3 mg/L yeast extract) supplemented with 0.7 M CaCl 2 and 20 µg/mL nalidixic acid (Nal). After two days, R. tropici cells were harvested and resuspended in sterile water at O.D 600 = 0.3. 1 mL of this bacterial suspension was used to individually inoculate common bean seedlings.
E. coli pRK2013 strain was used as a helper for triparental mating with R. tropici. pRK2013 cells were grown in LB broth supplemented with 30 µg/mL kanamycin overnight at 37 • C.
R. tropici::uidA strain was obtained by triparental mating using the E. coli DH5α strain carrying the plasmid pFAJ1700::placZ:uidA, containing the uidA gene under the lacZ promoter as a donor, and the pRK2013 as conjugation helper. R. tropici derivatives carrying the plasmid were selected in PY medium as Nal and Tc resistant transjugants. R. tropici::uidA cells were grown in PY medium supplemented with 0.7 M CaCl 2 , 10 µg/mL Tc, and 20 µg/mL Nal for 2 days at 30 • C. After two days, cells were harvested and resuspended in sterile water at OD 600nm = 0.3. 1 mL of this bacterial suspension was used to individually inoculate common bean seedlings.

Treatment
Three-day-old common bean seedlings were transferred into a nitrogen-free Fähraeus medium [53] plates supplemented with 1 mM PO 4 (optimal-Pi conditions) or 5 µM PO 4 (low-Pi conditions). Three days after transplanting, seedlings were inoculated with autoclaved-water (Mock) or with R. tropici CIAT899. Inoculated seedlings were kept under dark conditions at room temperature. At 1, 24 and 48 h post-inoculation, the zone of the roots susceptible to rhizobium infection (hereafter referred as susceptible zone) was isolated, flash frozen in liquid nitrogen, and stored at −80 • C until used for total RNA extraction. The susceptible zone encompassed a region of~1 cm and is analogous to the distal region of the elongation zone and the entire differentiated zone of the root [54]. Four biological replicates, each one containing susceptible zones from four different seedlings, were independently generated.

Quantification of Soluble Phosphate
Soluble Pi content was determined in roots from common bean seedlings grown for three days in optimal-or low-Pi conditions using the colorimetric assay reported in [55]. Briefly, roots were harvested, rinsed with ultra-pure water, weighed, immediately homogenized in 10 N trichloroacetic acid and fractioned by centrifugation. A soluble fraction was incubated with the staining solution Genes 2018, 9, 498 5 of 26 (350 mM FeSO 4 , 16% (NH 4 ) 6 Mo 7 O 24 ) for 10 min at room temperature. Optical density was determined at 660 nm. For each experimental condition, five replicates, each one containing four roots from different seedlings, were analyzed.

Root Hair Deformation Analysis
Three-day-old common bean seedlings grown for three days under optimal-Pi or low-Pi conditions were inoculated with 1 mL of R. tropici suspension (OD 600nm = 0.3). Forty-eight hours after inoculation, roots were collected and stained with methylene blue to maximize contrast, and then observed with a bright-field microscope. A total of 20 independent biological replicates were generated, each one including 10 plants.

Quantification of Number of Infection Threads
To quantify the number of infection threads per root, common bean seedlings grown for three days under optimal-Pi or low-Pi conditions were inoculated with 1 mL of R. tropici::uidA suspension (OD 600nm = 0.3). Three days post-inoculation, susceptible zones were collected, rinsed twice with autoclaved ultra-pure water and immersed in uidA staining solution (0.05% 5-bromo-4-chloro-3-indolxyl-β-D-glucuronic acid, 100 mM sodium phosphate buffer (pH7), 0.5 mM potassium ferrocyanide, 0.5 mM potassium ferricyanide, 10 mM Na 2 EDTA and 0.1% Triton X-100) and incubated for one hour at 37 • C. After rinsing in phosphate buffer, the number of infection threads formed at the susceptible zone was quantified under a bright-field microscopy. Additionally, to evaluate the rhizobia adhesion to the common bean roots, the intensity of the blue precipitate (as a product of the activity of the R. tropici::uidA) formed in the susceptible zone was densitrometrically quantified. To this end, pictures of susceptible zones showing blue precipitate were captured using a SZX10 stereomicroscope (Olympus, Center Valley, PA, USA) equipped with an Olympus UC50 camera (Olympus). The intensity of the blue precipitate in each susceptible zone was quantified as number of pixels by using the software CellSens Standard (Olympus). For these experiments, four biological replicates, each one with 10 roots from different seedlings, were included.

Gene Expression Analysis
To analyze the expression of PvNSP2, PvNIN, PvFLOT2, PvRIC1 and PvRIC2 genes, total RNAs were extracted from 0.5 g of rhizobia-inoculated or mock-inoculated susceptible zones from common bean seedlings growing under optimal-Pi or low-Pi conditions using ZR Plant RNA Miniprep kit (Zymo Research, Irvine, CA, USA) according to manufacturer's instructions. Genomic DNA (gDNA) was removed from purified RNA using DNaseI RNase-free (Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer's instructions. 1 µg of gDNA-free total RNA was used to synthesize cDNA using Thermo Scientific RevertAid Reverse Transcriptase (Thermo Fisher Scientific, USA) according to manufacturer's instructions. cDNA samples were used to analyze the expression of the aforementioned genes by qRT-PCR in a Step-One qPCR thermocycler (Applied Biosystems, Foster City, CA, USA). The housekeeping gene PvActin was used to normalize gene expression levels [56]. The expression level of different genes was calculated according to the equation E = nP eff (−∆Ct) where P eff is the primer set efficiency calculated using LinRegPCR [57] and ∆Ct is the difference of the cycle threshold (Ct) values of the housekeeping gene and a given gene. qRT-PCR primers were designed by using the online available software Primer 3 v.0.4.0 [58] and using the following parameters: 80-120 pb as product size, 19-23 nt as primer size, 60 • C as primer temperature melting, and 40-60% of primer GC content. To design primers, the sequence from the 3'UTR region of each gene was obtained from the P. vulgaris genome Version 2.1 [59,60]. The nucleotide sequences of the qRT-PCR primers used in this study are provided in Table S1. For this experiment, four biological replicates, each one containing roots from four different seedlings, were included.

Preparation of Messenger RNA-Seq Library and High-Throughput Sequencing
Total RNA was isolated from 0.5 g of 48 h rhizobia-inoculated or mock-inoculated roots from common bean plants growing under optimal-Pi or low-Pi conditions, as described earlier in this paper. Stranded messenger RNA-seq (mRNA-seq) libraries were generated from 1.5 µg of gDNA-free total RNA from each experimental condition and prepared using TruSeq RNA Sample Prep kit (Illumina, San Diego, CA, USA) according to the manufacturer's instructions. Briefly, sample concentration was determined by Qubit flourometer (Invitrogen, Carlsbad, CA, USA) using the Qubit HS RNA assay kit, and the RNA integrity was checked using the Fragment Analyzer automated electrophoresis system (Advanced Analytical, Ankeny, IA, USA). Poly-A containing mRNAs were purified from total RNAs and subsequently fragmented. cDNAs were generated from fragmented RNAs, and the index-containing adapters were ligated to their ends. The amplified cDNAs were purified by the addition of Axyprep Mag PCR Clean-up beads (Axygen, Corning, NY, USA). The quality of each purified RNA-seq library was evaluated using the Fragment Analyzer automated electrophoresis system, and quantified with the Qubit flourometer using the Qubit HS dsDNA assay kit.
In total, 12-stranded mRNA-seq [two Pi conditions (sufficient and deficient), two treatments (mock or rhizobia), and three independent biological replicates] were high-throughput sequenced in 75 bp single-end at the University of Missouri DNA Core Facility by using an Illumina NextSeq500 sequencer (Illumina, USA) according to the manufacturer's instructions. Image analysis and base calling were performed using the Illumina pipeline (http://www.illumina.com).

Mapping and Processing mRNA-Seq Reads
mRNA-seq analysis was conducted using a custom in-house developed informatics pipeline which first performs a quality check on the raw sequencing FASTQ files using the tool FastQC [61] and trim-galore tool [62] to remove low-quality reads and trim adaptor sequences. We then indexed the P. vulgaris Version 2.1 [59,60] reference genome using Bowtie2 short read aligner [63]. The trimmed mRNAseq reads were then aligned, allowing two mismatches, to this index reference genome using TopHat, which also reports splice junction sites [64].

Quantification and Identification of Differentially Regulated Genes
The next step in our pipeline was to quantify the expression level of each transcript/gene by measuring the level of mRNA-seq alignments using the tool Cufflinks [65]. The quantified expression level was represented by Fragments per Kilobase of transcript per Million mapped reads (FPKM) values. The differential expression levels of quantified genes were then calculated using the tool Cuffdiff [65] which takes the observed log fold change of a gene's expression from control and treatment samples and then reports if the gene is significantly differentially expressed (q-values). The thresholds for Cuffdiff were set at 2-fold change and q-value ≤ 0.05 before calling genes to be differentially expressed, which are usually annotated with their unique identifiers defined by the gene models annotation and gene coordinates files. Data obtained from this differentially expression analysis was independently validated by qRT-PCR.

Gene Functional Classification
The biological relevance of the differentially regulated genes was assessed by a gene function enrichment analysis using the method Singular Enrichment Analysis (SEA) available in the web-based tool AgriGo [66,67], followed by a summarization and visualization of statistically significant non-redundant gene ontology (GO) terms by using the web-based tool REVIGO [68,69]. Briefly, GO terms enriched in each set of genes were compared to the P. vulgaris Version 2.1 gene reference background [59,60]. p-values for the GO terms were obtained through Fisher's exact test, and a q-value was computed to produce lists of significant GO terms with an estimated False discovery rate (FDR) of 5%. Enriched GO terms with q < 0.05 were further summarized and visualized on REVIGO. Additionally, MapMan gene functional classification was used [70,71]. For MapMan analysis, the common bean mapping file, available at MapMan website [72], was used.

Statistical Analyses
All the statistical analyses were conducted using R software 3.0.1. The specific tests performed are indicated in the legend of the corresponding figure.

Rhizobia-Induced Root Hair Deformation and Infection Thread Formation is Affected in Pi-Deficient Common Bean Seedlings
It has been demonstrated that the reduction in nodule number is one of the negative effects of Pi deficiency in legumes [40]. However, it is unclear whether this nutritional stress interferes with the molecular dialogue between the common bean and rhizobia; therefore, the precise stages of this molecular dialogue that are potentially affected by Pi deficiency remain unknown. To better characterize the effect of the Pi deficiency in the molecular dialogue between the common bean and rhizobia, we first developed a system which allows us to track early molecular and physiological responses to rhizobia under Pi limiting conditions. This system allows the growth of common bean seedlings under optimal-(1 mM) or low-(5 µM) Pi conditions, as well as under sterile and controlled environmental conditions ( Figure S1). By using this system, we were able to obtain Pi-deficient common bean seedlings, which led to a 50% reduction in plant Pi contents ( Figure S1). Once we confirmed that our experimental system yields Pi-deficient common bean seedlings, we proceeded to assess early responses to rhizobia in common bean seedlings growing under optimal-or low-Pi conditions.
One of the early physiological responses triggered by rhizobia is the root hair deformation, which is required for rhizobia entrapping and the formation of the infection chamber [9,10]. To test whether Pi deficiency affects the rhizobia-induced root hair deformation and infection thread formation, common bean seedlings growing for three days under optimal-or low-Pi conditions were inoculated with R. tropici CIAT899::uidA. Forty-eight hours post inoculation, 95% (n = 190) of the seedlings growing under optimal-Pi conditions showed rhizobia-induced root hair deformation (Figure 1a,c). In contrast, only 50% (n = 200) of the Pi-deficient seedlings displayed deformed root hairs (Figure 1b,c). Interestingly, most of the deformations observed in Pi-deficient common bean seedlings were aberrant ( Figure 1d). For instance, we observed root hairs with swelled tip, root hairs with more than one tip outgrowth, and root hairs with spatula-like deformation ( Figure 1e). Additionally, we observed that Pi deficiency significantly affected rhizobia attachment, which was reflected by a reduction in the blue precipitate (as a product of the activity of the R. tropici::uidA strain) at the susceptible zone (Figure 2a,b). Accordingly, we observed that Pi-deficient common bean seedlings developed 60% less infection thread than the seedlings growing under optimal-Pi condition and inoculated with rhizobia (Figure 2c-e). Similar aberrant rhizobia-induced root hair deformations have been observed in M. truncatula ern1/ern2, M. truncatula dmi1 and in L. japonicus scarn mutant plants [9,53,73,74]. This defect has been associated with the inhibition in the polar growth of the root hair cells, which affects the rhizobia infection process including the infection thread formation [9,73,74]. Additionally, an early study reported that Pi deficiency affects the rhizobia attachment to Medicago roots [75]. Similarly, our data suggest that Pi deficiency negatively affects the polar elongation of the root hairs, likely leading to inefficient attachment and entrapping of rhizobia.  seedlings growing under OP conditions or LP conditions. Box plots represent first and third quartile (horizontal box side), minimum and maximum (outside whiskers). Asterisk indicates a significant difference according to one-way analysis of variance (ANOVA) (p-value < 0.001). Data showed was obtained from 20 biological replicates, each one containing 10 roots from different common bean seedlings. (e) Types of aberrant rhizobia-induced root hair deformation observed in common bean seedlings growing under LP conditions. Scale bars in panel represent 50 μm. . Asterisk indicates a significant difference according to one-way ANOVA (p-value < 0.001). Data showed was obtained from 20 biological replicates, each one containing ten roots from different common bean seedlings.

Pi Deficiency Decreases the Expression of Genes Required for Root Hair Curling and Rhizobia Infection
The coordinated action of different transcription factors is essential for a successful symbiosis with rhizobia. Upon NF perception, the transcription factors NSP2 and NSP1 form a DNA binding complex [24,76]. This complex positively regulates the expression of the transcription factor NIN, which is required to activate the expression of different genes involved in rhizobia infection and colonization processes [77]. It has been demonstrated that mutations in the NSP2 gene significantly affect the rhizobia-induced root hair deformation and rhizobia infection process, particularly in the formation of the infection chamber and infection thread [74,78]. Likewise, early reports have Box plots represent first and third quartile (horizontal box side), minimum and maximum (outside whiskers). Asterisk indicates a significant difference according to one-way ANOVA (p-value < 0.001). Data showed was obtained from 20 biological replicates, each one containing ten roots from different common bean seedlings.

Pi Deficiency Decreases the Expression of Genes Required for Root Hair Curling and Rhizobia Infection
The coordinated action of different transcription factors is essential for a successful symbiosis with rhizobia. Upon NF perception, the transcription factors NSP2 and NSP1 form a DNA binding complex [24,76]. This complex positively regulates the expression of the transcription factor NIN, which is required to activate the expression of different genes involved in rhizobia infection and colonization processes [77]. It has been demonstrated that mutations in the NSP2 gene significantly affect the rhizobia-induced root hair deformation and rhizobia infection process, particularly in the formation of the infection chamber and infection thread [74,78]. Likewise, early reports have demonstrated that M. truncatula nin mutants undergo excessive rhizobia-induced root hair curling but are impaired in infection and nodule formation [77].
Pi deficiency reduces both rhizobia-induced root hair deformation and infection thread formation in common bean seedlings suggesting that the expression of genes participating in the root hair curling and rhizobia infection is compromised under this nutritional condition. To test this hypothesis, we evaluated the expression of PvNSP2 (Phvul.009G122700), PvNIN (Phvul009G115800), and PvFLOT2 (Phvul.009G090700) in common bean seedlings growing under optimal-and low-Pi conditions and inoculated with rhizobia. Our analysis revealed that the expression of these three symbiosis-related genes was induced upon one-, 24, and 48 h post-inoculation with R. tropici in common bean seedlings growing under optimal Pi conditions ( Figure 3). In contrast, we observed a 50% reduction in the expression of these three symbiotic genes in Pi-deficient common bean seedlings inoculated with rhizobia ( Figure 3). Interestingly, at 48 h post-inoculation with rhizobia, PvNIN was the only gene that showed a substantial induction (in average by 5 fold-change) in Pi-deficient seedlings in response to rhizobia (Figure 3c).
To explore the mechanism behind the reduction in the fold-change expression of these three symbiosis-related genes in response to Pi-deficiency, we analyzed their expression levels in response to mock and rhizobia treatment (Figure 3b,d,f). By comparing the expression levels of PvNSP2 and PvNIN from both optimal-Pi and Pi-deficient seedlings, we did not find significant differences at one and 24 h post-mock treatment (Figure 3b,d). Similarly, we did not find significant differences in the expression levels of PvFLOT2 across the entire time-course when comparing optimal-and low-Pi seedlings post-mock treatment (Figure 3f). At 48 h, we observed an increase in the expression levels of PvNSP2 and a diminution in the expression levels of PvNIN in mock-inoculated Pi-deficient seedlings. By comparing the expression levels obtained from mock-and rhizobia-inoculated Pi-deficient seedlings, we observed a significant increase in the expression levels of PvNSP2 and PvNIN upon 24 and 48 h of inoculation with rhizobia (Figure 3b,d). Altogether, these transcriptome analyses suggest that the defects in root hair deformation and infection thread formation observed in common bean seedlings grown under Pi-deficient conditions are the consequences of an inefficient activation of critical transcription factors like PvNIN and PvNSP2, and symbiotic genes like PvFLOT2. The fact that substantial induction of the expression of PvNIN in Pi-deficient seedlings was detected upon 48 h post-inoculation with rhizobia only also indicates that Pi deficiency delays the activation of critical symbiosis-related genes.
Genes 2018, 9, x FOR PEER REVIEW 10 of 25 symbiosis-related genes was induced upon one-, 24, and 48 h post-inoculation with R. tropici in common bean seedlings growing under optimal Pi conditions ( Figure 3). In contrast, we observed a 50% reduction in the expression of these three symbiotic genes in Pi-deficient common bean seedlings inoculated with rhizobia ( Figure 3). Interestingly, at 48 h post-inoculation with rhizobia, PvNIN was the only gene that showed a substantial induction (in average by 5 fold-change) in Pi-deficient seedlings in response to rhizobia (Figure 3c).

Pi Deficiency Increases the Expression of the PvRIC1 and PvRIC2 Genes in Common Bean Seedlings
It has been reported that the expression of PvRIC1 and PvRIC2 is induced one-and five-days post-inoculation with rhizobia, respectively [32]. Based on the fact that Pi deficiency triggers a reduction in the number of nodules in common bean [40] and that the peptides PvRIC1 and PvRIC2 actively participate in the AON process, we hypothesized that the expression of PvRIC1 and PvRIC2 might be affected in Pi-deficient common bean seedlings. To test this hypothesis, we quantified the expression of PvRIC1 (Phvul.005G096901) and PvRIC2 (Phvul011G135900) in common bean seedlings grown under optimal-or low-Pi conditions, and under mock and rhizobia-treated conditions. This analysis revealed that the expression of PvRIC1 is induced as early as one-hour post inoculation with rhizobia in both control and Pi-deficient common bean seedlings and have its maximum induction level at 48 h post-inoculation with rhizobia ( Figure 4a). In contrast, we were not able to detect a consistent and reproducible induction of PvRIC2 at 1 and 24 h post-inoculation with rhizobia. However, we observed a clear induction of the expression of PvRIC2 in response to rhizobia in both control and Pi-deficient common bean seedlings after 48 h post-inoculation (Figure 4c). Although the induction of the expression of these two AON-related genes was observed in both control and Pi-deficient seedlings, we observed that the expression levels of both PvRIC1 and PvRIC2 genes were always higher in Pi-deficient seedlings than in control seedlings (Figure 4b,d). Additionally, by comparing the expression levels of mock-inoculated seedlings, we observed that Pi-deficient seedlings always showed higher expression levels similar to those observed in control-seedlings inoculated with rhizobia (Figure 4b,d). Altogether, these data indicate that Pi deficiency induces the expression of PvRIC1 and PvRIC2 in the absence of rhizobia, and that the expression of these AON-related genes is highly induced in response to rhizobia in Pi-deficient common bean seedlings. Under this scenario, it is likely that: (1) Pi-deficient common bean seedlings are preconditioned to reduce the number of rhizobia infections. This hypothesis can be supported by the fact that the overexpression of RIC1 and RIC2 inhibits nodulation in soybean [79]. (2) Higher expression of PvRIC1 and PvRIC2 in response to Pi deficiency might be critical for common bean to properly adapt to Pi deprivation. This hypothesis is supported by the fact that other CLE peptides (e.g., CLE14) [80] play a critical role in the plant adaptation to Pi deficiency [81,82].  or LP (light-gray) conditions. Despite several attempts, we were not able to detect reproducible and consistent expression (non-detected: ND) of PvRIC2 at one and 24 h post-inoculation with mock or rhizobia. Data showed was obtained from four independent replicates, each one containing susceptible zones from four different plants. Box plots represent first and third quartile (horizontal box side), minimum and maximum (outside whiskers). One-way ANOVA, followed by a Tukey HSD test was performed (p-value < 0.01). Statistical classes sharing a letter are not significantly different.

mRNA-Seq Analysis
To obtain a better insight into the genetic responses of common bean plants leading to a reduced nodulation under Pi-deficient conditions, we conducted an mRNA-seq analysis. For this genomewide transcriptional analysis, based on our transcriptional data obtained by qRT-PCR, we selected 48 h post-inoculation to capture more transcriptional responses to rhizobia.
A total of 12 cDNA libraries were generated from control (optimal Pi: OP) and Pi-deficient roots conditions. Despite several attempts, we were not able to detect reproducible and consistent expression (non-detected: ND) of PvRIC2 at one and 24 h post-inoculation with mock or rhizobia. Data showed was obtained from four independent replicates, each one containing susceptible zones from four different plants. Box plots represent first and third quartile (horizontal box side), minimum and maximum (outside whiskers). One-way ANOVA, followed by a Tukey HSD test was performed (p-value < 0.01). Statistical classes sharing a letter are not significantly different.

mRNA-Seq Analysis
To obtain a better insight into the genetic responses of common bean plants leading to a reduced nodulation under Pi-deficient conditions, we conducted an mRNA-seq analysis. For this genome-wide transcriptional analysis, based on our transcriptional data obtained by qRT-PCR, we selected 48 h post-inoculation to capture more transcriptional responses to rhizobia.
A total of 12 cDNA libraries were generated from control (optimal Pi: OP) and Pi-deficient roots (low Pi: LP) inoculated 48 h with mock (OPM or LPM, respectively) or rhizobia (OPR or LPR, respectively). Three independent biological replicates were generated for each condition. These libraries were sequenced using the Illumina NextSeq500 platform. After filtering low-quality reads, a total of 176,365,724 single-end reads (75 bp) were aligned to the P. vulgaris genome reference sequence Version 2.1 [59,60] using Bowtie and Tophat software and allowing two mismatches [64]. Of these, 164,037,947 were uniquely mapped to the common bean genome and were used for further analysis. The full RNA-seq datasets generated from the 12 cDNA libraries used in this study were deposited in the Gene Expression Omnibus [83] under the accession number GSE118968.

Global Transcriptional Responses of Pi-Deficient Common Bean Seedlings Interacting with Rhizobia
In order to identify additional genes with a potential role in the regulation of the early stages of the symbiosis between common bean and rhizobia under low-Pi conditions, the mRNA-seq data sets were analyzed using Cuffdiff [84] with an additional cutoff of 2-fold in pairwise comparisons (e.g., OPR/OPM). In total, 867 (511 up-regulated and 356 down-regulated) and 383 (285 up-regulated and 98 down-regulated) differentially regulated genes were identified when comparing OPR/OPM and LPR/LPM, respectively (Tables S2 and S3). Among them, PvNSP2 (Phvul.009G122700) and PvFLOT2 (Phvul.009G090700) were also classified as differentially regulated by mRNA-seq approach (Tables S2 and S3). However, we did not detect the expression of PvNIN (Phvul009G115800) in our mRNA-seq analysis, which we explain by a difference in the sensitivity of mRNA-seq and qRT-PCR approaches.
To confirm these mRNA-seq results, the expression of 25 randomly selected genes was analyzed via qRT-PCR and using independent biological material from the one used for mRNA-seq analysis ( Figure S2). The expression pattern obtained by qRT-PCR showed the same trend observed by mRNA-seq in response to rhizobia. We found some differences in the fold-change computed by the two independent methods, which we explain by the relative sensitivity and data processing (e.g., algorithms used for data normalization) of each method, along with technical and biological aspects (e.g., efficiency/specificity of qRT-PCR and inherent variation in the responses of common bean to both Pi availability and rhizobia). In spite of these differences, the fold-changes computed for each analyzed gene by qRT-PCR were statistically different between both Pi conditions. Upon confirmation of the reliability of our mRNA-seq datasets, we performed additional pairwise comparisons to identify genes that might explain the reduction in the nodule number triggered by Pi deficiency. To this end, we compared the differentially regulated genes found in the pairwise comparisons OPR/OPM and LPR/LPM to identify commonly and uniquely regulated genes under these Pi conditions and interacting with rhizobia. This comparative analysis revealed that 234 (180 up-regulated and 54 down-regulated) genes were commonly regulated in both Pi conditions while interacting with rhizobia (Figure 5a,b and Table S4). In contrast, 633 (331 up-regulated and 302 down-regulated) and 149 (105 up-regulated and 44 down-regulated) genes were uniquely regulated in Pi-optimal and Pi-deficient seedlings, respectively (Figure 5a,b and Tables S5 and S6).

Expression of Hormone-and Signal-Transduction-Related Genes Is Compromised in Rhizobia Inoculated Pi-Deficient Common Bean Plants
To gain more insights into the molecular functions of the differentially regulated genes identified in this study, we performed a gene-enrichment analysis and a functional classification in MapMan to identify the molecular mechanisms leading to an inefficient symbiosis between common bean and rhizobia under Pi-deficient conditions. This analysis indicated that Regulation of transcription, Hormones, Regulation/Signal transduction, Transport, Biotic interactions and Protein modification/degradation were the most overrepresented functional categories in both control and Pi-deficient common bean seedlings interacting with rhizobia ( Figure 6). Further hierarchical clustering analyses on the 234 common regulated genes revealed a clear variation in the expression levels computed for these genes (Figure 5c,d). For instance, we identified two clusters showing obvious differences when comparing the expression levels computed for both Pi-optimal and Pi-deficient seedlings inoculated with rhizobia ( Figure 5c). The first identified cluster (cluster I) contains genes whose expression levels were higher in rhizobia-inoculated Pi-deficient seedlings than in rhizobia-inoculated optimal-Pi seedlings, whereas the second cluster (cluster II) contains genes whose expression levels were higher in optimal-Pi seedlings than in Pi-deficient seedlings, both of them inoculated with rhizobia (Figure 5c). To extend this analysis and identify genes differentially regulated between Pi-optimal and Pi-deficient seedlings while interacting with rhizobia, we performed a pair-wise comparison of the expression levels between rhizobia-inoculated optimal-Pi and rhizobia-inoculated Pi-deficient conditions. This comparison led us to the identification of 131 differentially expressed genes (77 showing higher expression in Pi-deficient seedlings and 54 showing higher expression in control seedlings) ( Figure S3 and Table S7). The observed variation in the expression of these 131 genes was on average five-fold. For instance, Phvul.001G226900 and Phvul.008G244400, which encode for a Wall-associated receptor kinase galacturonan-binding Phytohormones play a critical role in different stages of the symbiosis between legumes and rhizobia [85]. For instances, it has been demonstrated that both auxin and cytokinin positively regulate the expression of genes that conforms the common symbiosis pathway, as well as the progression of the symbiotic signaling events from the epidermis to the cortex [27,86]. Recently, it was demonstrated that the auxin accumulation in the root hairs upon rhizobia recognition is essential for the infection thread formation [10]. Furthermore, it has also been demonstrated that the correct perception of cytokinin is crucial for the polar auxin transport and the subsequent auxin accumulation in the rhizobia-infected root hairs [87].
Based on the fact that auxin and cytokinin are crucial regulators of both early and late (e.g., seedlings were grown under Pi-deficient condition (Figure 7). The fact that Pi deficiency affects the expression of auxin-and cytokinin-related genes in response to rhizobia, lead us to propose that the defects in early molecular (e.g., expression of PvNIN and PvFLOT2) and physiological (e.g., root hair deformation) events observed in Pi-deficient common bean seedlings might be partially explained by defects in the auxin-cytokinin balance.  In addition to auxin-and cytokinin-related genes, we also found that genes involved in the biosynthesis of the phytohormones ethylene, abscisic acid (ABA)-and jasmonate, as well as in the signal transduction related to these three phytohormones, were differentially regulated in response to rhizobia inoculation. Although we observed that common bean seedlings growing under optimal-Pi conditions and inoculated with rhizobia showed more differentially regulated genes (six ABA-related genes, seven ethylene-related genes, and six jasmonate-related genes) related to these three hormones than Pi-deficient seedlings inoculated with rhizobia, the majority of them were down-regulated ( Figure 8).
In contrast, all ethylene-, ABA-and jasmonate-related genes identified in common bean seedlings grown under Pi-deficient conditions and inoculated with rhizobia were up-regulated ( Figure 8). related genes, seven ethylene-related genes, and six jasmonate-related genes) related to these three hormones than Pi-deficient seedlings inoculated with rhizobia, the majority of them were downregulated ( Figure 8). In contrast, all ethylene-, ABA-and jasmonate-related genes identified in common bean seedlings grown under Pi-deficient conditions and inoculated with rhizobia were upregulated ( Figure 8).  Ethylene, jasmonate, and ABA are essential to activate a series of molecular and physiological responses allowing plants to cope with Pi deficiency [89]. However, growing evidence indicates that these three phytohormones negatively control the establishment of the root nodule symbiosis [90][91][92]. For instances, it has been demonstrated that the accumulation of ethylene or jasmonate inactivates the expression of early symbiosis-related genes (e.g., ENOD11), inhibits the initiation and maintenance of the calcium spiking, decreases the number of rhizobia-induced root hair deformation and infection thread, and reduces nodule formation in M. truncatula plants [90,91]. Recently, it was demonstrated that the negative impact of ethylene on the rhizobia infection process and nodule formation is under the control of the regulator EIN2 [28]. Hence, our transcriptional data suggest that the activation of genes related to the biosynthesis and signal transduction of these three phytohormones is required in order to allow common bean to properly adapt to Pi deficiency. However, the activation of the ethylene-, jasmonate-and ABA signaling pathways compromise the activation of early molecular and physiological events of the root nodule symbiosis, as we demonstrated in this study.
The proper recognition of rhizobia by the host legume is crucial for a successful symbiosis. To this end, legumes rely on the LysM-domain receptor kinases NFR5 [93], NFR1 [94] and eNFR [7], as well as in the LRR receptor-like kinase SYMRK, to properly detect and decode the NF signal [15,16]. Recently, it was demonstrated that a fourth LysM receptor kinase receptor, EPR3, is required for the perception of compatible rhizobial exopolysaccharides and promotes a successful infection and colonization [95]. Additionally, based on the fact that a gene encoding cell wall-associated receptor kinases (WAKs) is significantly up-regulated upon one-hour treatment with NF in M. truncatula roots, it has been proposed that WAKs might play a critical role in the early steps of this symbiosis, primarily in the cell wall-cytoplasm signaling [96]. In our transcriptomic analysis, we also found different types of protein receptors, including LysM receptor kinases, LRR receptor kinases, and WAKs, that were differentially regulated in both control and Pi-deficient seedlings inoculated with rhizobia. However, we observed that Pi-deficient seedlings inoculated with rhizobia showed less differentially regulated genes encoding for protein receptors in response to rhizobia ( Figure S4). These data suggest that the misregulation of the expression of additional protein receptor might contribute to an inefficient communication between common bean and rhizobia under Pi-limiting conditions.

Conclusions
The data presented in this study lead us to conclude that Pi deficiency negatively affects early events of the symbiosis between legumes and rhizobia. The negative effects in the symbiosis between the common bean and rhizobia under Pi deficiency might be triggered by an imbalance between auxin and cytokinin, as well as by the activation of the ethylene, jasmonate and ABA-related signaling pathway, which are negative regulators of early molecular events of this symbiosis ( Figure 9). Finally, the reduction in the number of nodules observed under Pi-deficient conditions can be also explained by a constitutive activation of the AON process ( Figure 9); however, further experimentation is needed to confirm this hypothesis.  Figure 9. Proposed model of the effects of Pi deficiency in the early steps of the symbiosis between common bean and rhizobia. Higher induction, lower induction and down-regulation of the expression of different rhizobia-induced genes is indicated in red, orange, and green, respectively. Upon NF perception, a variety of molecular responses occurs. For example, the accumulation of mevalonate at the cytoplasm is required for the activation of the Ca 2+ spiking [97]. These fluxes of calcium are further transduced by CCaMK, which in turn activates the transcription factor CYCLOPS. In turn, CYCLOPS is required to activate the transcription factors NSP2/NSP1 and NIN. The coordinated action of these transcription factors activates the expression of different genes (e.g., FLOT2) required for the infection of the root by rhizobia. There is evidence indicating that these molecular events are positively regulated by a delicate balance among the phytohormones auxins (IAA), cytokinins (CK) and ethylene (Eth). These molecular events occur during the first twenty-four hours of interaction between legumes and rhizobia. At this stage, the expression of RIC1 and RIC2 genes, which participate in the AON process, is not fully activated. The coordinated action of these molecular events leads to a successful symbiosis. Based on our data, Pi deficiency might affects the phytohormone balance by activating the expression of genes participating in the biosynthesis of abscisic acid (ABA), ethylene (Eth), and jasmonate (JA). These three hormones negatively affect the Ca 2+ spiking [90,91], leading to a reduction in the expression of symbiosis-related genes (e.g., NSP2, NIN, and FLOT2) required for the rhizobia infection process. Additionally, Pi deficiency induces the Figure 9. Proposed model of the effects of Pi deficiency in the early steps of the symbiosis between common bean and rhizobia. Higher induction, lower induction and down-regulation of the expression of different rhizobia-induced genes is indicated in red, orange, and green, respectively. Upon NF perception, a variety of molecular responses occurs. For example, the accumulation of mevalonate at the cytoplasm is required for the activation of the Ca 2+ spiking [97]. These fluxes of calcium are further transduced by CCaMK, which in turn activates the transcription factor CYCLOPS. In turn, CYCLOPS is required to activate the transcription factors NSP2/NSP1 and NIN. The coordinated action of these transcription factors activates the expression of different genes (e.g., FLOT2) required for the infection of the root by rhizobia. There is evidence indicating that these molecular events are positively regulated by a delicate balance among the phytohormones auxins (IAA), cytokinins (CK) and ethylene (Eth). These molecular events occur during the first twenty-four hours of interaction between legumes and rhizobia. At this stage, the expression of RIC1 and RIC2 genes, which participate in the AON process, is not fully activated. The coordinated action of these molecular events leads to a successful symbiosis. Based on our data, Pi deficiency might affects the phytohormone balance by activating the expression of genes participating in the biosynthesis of abscisic acid (ABA), ethylene (Eth), and jasmonate (JA). These three hormones negatively affect the Ca 2+ spiking [90,91], leading to a reduction in the expression of symbiosis-related genes (e.g., NSP2, NIN, and FLOT2) required for the rhizobia infection process. Additionally, Pi deficiency induces the expression of RIC1 and RIC2, even in the absence of rhizobia. All these defects in the early stages of this symbiosis lead to a reduction in the number of nodules in common bean.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2073-4425/9/10/498/s1, Figure S1: In vitro system to assess early responses to rhizobia on Pi-deficient common bean seedlings. Seedlings were grown under OP (a) or LP (b) conditions. Seedlings growing under LP conditions show a 50% decrease in their Pi content (c). Asterisk indicates a significant difference according to One-way ANOVA (p-value < 0.001), Figure S2: qRT-PCR validation of the mRNA-seq datasets. Data showed are the fold-change obtained from the pairwise comparison between the expression levels obtained from seedlings grown under optimal-Pi (OP) or low-Pi (LP) inoculated with rhizobia (R; OPR or LPR) or mock (M; OPM or LPM), Figure S3: Genes showing significantly higher expression levels in Pi-deficient common seedlings inoculated with rhizobia. Data showed are the expression levels obtained from common bean seedlings grown under optimal-Pi (OP) or low-Pi (LP) conditions and inoculated with rhizobia. Panels on the right site show the functional categories regulated under the described experimental conditions. TF: transcription factors; PM: protein modification; PD: protein degradation; IAA: auxin; ABA: abscisic acid; Eth: ethylene; CK: cytokinins; JA: jasmonate; SA: salicylic acid; GA: gibberellic acid; RK: receptor kinases; CaR: calcium regulation; GP: G-protein; Pi: phosphoinositides; C and Nutr: C and nutrients; MK: MAP kinases; U: unspecified; Thx: thioredoxin; Asc/Glut: arcorbic/gluath; Glutx: glutaredoxin; Px: periredoxin; Dism/Cat: dismutase/catalase, Figure S4: Expression pattern of protein receptors-and secondary metabolism-related genes. Data showed are the fold-change obtained from the pairwise comparison between the expression levels obtained from seedlings grown under optimal-Pi (OP) or low-Pi (LP) inoculated with rhizobia (R; OPR or LPR) or mock (M; OPM or LPM). Fold-changes were false-colored; genes showing higher expression are shown in different scales of red color, whereas genes showing lower expression were colored in green color, Table S1: Primer sequences used for quantitative real-time polymerase chain reaction (qRT-PCR) analysis, Table S2: Rhizobia-induced genes under optimal-Pi conditions, Table S3: Rhizobia-induced genes under low-Pi conditions, Table S4: Commonly regulated genes in both optimal-and low-Pi conditions, Table S5: Genes uniquely regulated under optimal-Pi conditions, Table S6: Genes uniquely regulated under low-Pi conditions, Table S7