Comparative Analysis of Soybean Root Proteome Reveals Molecular Basis of Differential Carboxylate Efflux under Low Phosphorus Stress

Carboxylate efflux from roots is a crucial and differential response of soybean genotypes to low phosphorus (P) stress. Exudation of carboxylic acids including oxalate, citrate, succinate and fumarate was induced under low P stress, particularly in P-efficient soybean genotypes. Enhancement of root length, surface area and volume further improved P acquisition under low P stress. To understand the molecular basis of carboxylate efflux under low P stress, the root proteome of contrasting genotypes (P-efficient: EC-232019 and P-inefficient: EC-113396) was compared. Among a total of 325 spots, 105 (32%) were differentially abundant proteins (DAPs) between sufficient (250 µM) and low P (4 µM) levels. Abundance of 44 (14%) proteins decreased by more than two-fold under low P stress, while 61 (19%) proteins increased by more than two-fold. Protein identification and annotation revealed that the DAPs were involved in a myriad of functions including carboxylic acid synthesis, carbohydrate, protein and lipid metabolism. Proteins with significant abundance included malate dehydrogenase, isocitrate dehydrogenase, phosphoglucomutase, phosphoglycerate mutase, fructokinase, enolase, phosphoglycerate kinase, triosephosphate isomerase, alcohol dehydrogenase, glucan water dikinase, glutamine synthetase and argininosuccinate lyase. Inferences from proteomic analysis suggests the crosstalk between various metabolic pathways implicated in conferring superior P acquisition efficiency under stress.


Introduction
Phosphorus (P) is an essential element for plant growth and development, with structural (nucleic acids, phospholipids), metabolic (energy transfer) and regulatory functions. P nutrition positively affects crop growth, with significant influence on above-(leaf area, dry matter accumulation, leaf P content, photosynthesis) and below-ground (root morphological traits, root exudation, symbiotic association) processes [1]. In soybean (Glycine max (L.) Merr.), P nutrition is important owing to its direct effect on growth and morphological responses and indirectly influencing nodulation and N 2 -fixation, and ultimately yield. Cultivating P-efficient soybean genotypes that can utilize the unavailable forms of soil P might be a sustainable option to increase crop productivity in the face of dwindling P reserves.
Altered root architecture and morphology, along with increased synthesis and excretion of carboxylic acids is crucial to enhance P acquisition [2][3][4]. Carboxylate ions release P from Upon emergence of cotyledonary leaf, seedlings were transferred to nutrient solution with two P levels: sufficient (250 µM) and low P (4 µM). Seedlings were supported on a 5 cm thick styrofoam sheet at a spacing of 3 cm × 3 cm. The cotyledonary leaves were removed on third day of transfer to nutrient solution to minimize genotypic variation due to seed P content. Thirty-six such seedlings were accommodated in individual containers. Three replications with four seedlings each were maintained for all treatment combinations. Styrofoam sheet was placed in a plastic container with 10 L of basal nutrient solution. Composition of the nutrient solution was as described in Vengavasi et al. [27], and pH was maintained at 6.5 using either 1.0 N HCl or 1.0 N KOH. The solution was aerated continuously and renewed every third day. The experiment was conducted in greenhouse at the National Phytotron Facility, ICAR-IARI, New Delhi, India with day/night temperature of 30/26 • C, photoperiod of 12 h at a photon flux density of 850 µmol·m −2 ·s −1 and relative humidity of 85%.

Growth Traits and Tissue Phosphorus Status
Roots of 20-d-old plants were scanned using root scanner (Regent Instruments Inc., Québec, QC, Canada) ( Figure S1A) and the images were analyzed in WinRhizo Pro software (Regent Instruments, Ville de Québec, QC, Canada) to obtain total root length (cm plant −1 ), surface area (cm 2 plant −1 ), volume (cm 3 plant −1 ) and number of root tips. Root and shoot P concentration (µg g −1 dry weight (DW)) was determined by ascorbic acid method [28] after digestion of dried tissue with di-acid mixture (HNO 3 :HClO 3 ::9:4). Tissue P content was calculated from P concentration in root or shoot, the sum of which was the total P uptake, expressed as µg P plant −1 .

Collection and Quantification of Root Carboxylate Efflux
Twenty-d-old plants were placed in 100 mL Erlenmeyer flasks with their roots immersed well in 50 mL of trap solution (0.5 mM CaCl 2 , pH 4.5) ( Figure S1B). The flasks were covered with black paper and exudates were collected from 08:00 to 12:00 h. Roots were washed in deionized water and blotted dry to determine fresh weight. The exudates were passed through Whatman no. 1 filter paper and frozen immediately to avoid microbial degradation. Root exudate was processed [26] and quantified through high-performance liquid chromatography (HPLC) (Agilent 1200 Infinity, Agilent Technologies, Palo Alto, CA, USA) with Hi-Plex H column as the stationary phase. The column temperature was set at 70 • C. Mobile phase (0.005 M H 2 SO 4 ) was used at a flow rate of 0.6 mL min −1 . Individual samples were run for 25 min and peaks captured by a refractive index detector with an optical temperature of 55 • C. Concentration of acids was quantified from the calibration curves of standards and expressed in µmol g −1 root fresh weight.

Protein Isolation and Quantification
Two g of root tissue was finely ground using liquid nitrogen in a pre-chilled pestle and mortar. Fifteen mL of extraction buffer containing 50 mM 2-[4-(2-hydroethyl)piperazin-1-yl] ethanesulfonic acid (HEPES) pH 7.5, 40% (w/v) sucrose and 0.1% (v/v) β-mercaptoethanol was added to the root powder and thoroughly homogenized. The homogenate was transferred to an Oak-ridge tube and 15 mL of phenol (equilibrated with 10 mM Tris pH 8.0 and 1 mM Ethylene diamine tetraacetic acid (EDTA) was added to it. After thorough mixing for 30 min, the homogenate was centrifuged at 5000 rpm for 10 min at 4 • C. Phenol phase was transferred to a fresh Oak-ridge tube to which 40 mL of precipitation buffer (0.1 M ammonium acetate in 100% (v/v) methanol) was added and incubated overnight at −20 • C. The protein pellet was obtained by centrifugation at 15,000 rpm for 30 min at 4 • C. The pellet was washed with 40 mL of 80% (v/v) acetone by vortexing and centrifuged at 15,000 rpm for 30 min at 4 • C. Washing step was repeated thrice after which the pellet was evaporated to dryness in a vacuum evaporator and stored at −80 • C until further analysis.

Two-Dimensional Gel Electrophoresis and Image Analysis
Immobilized pH gradient (IPG) strips (ReadyStrip TM , BioRad, Hercules, CA, USA) of 11.0 cm with pH 4 to 7 were passively rehydrated overnight with 400 µg of protein. Isoelectric focusing was performed in the Protean i12 IEF cell (BioRad, Hercules, CA, USA) with the following program: 200 volts for 180 min, 500 volts for 30 min, 1000 volts for 30 min, 6000 volts for 60 min, 6000 volts until a total of 3000 volt-hours followed by 500 volts for 20 h. After completion of isoelectric focusing, the strips were equilibrated for 15 min each in reducing buffer containing 50 mM Tris-HCl (pH 8.8), 6 M urea, 30% (v/v) glycerine, 2% (w/v) SDS and 20 mM DTT followed by alkylation buffer composed of 50 mM Tris-HCl (pH 8.8), 6 M urea, 30% (v/v) glycerine, 2% (w/v) SDS and 135 mM Iodoacetamide. Separating gel containing 12% (w/v) acrylamide, 375 mM Tris-HCl (pH 8.8), 0.1% (w/v) ammonium persulfate and 0.04% Tetramethylethylenediamine (TEMED) was poured into gel slabs. After polymerization, the equilibrated IPG strips were loaded on to the separating gel along with a protein ladder (11 to 245 kDa) placed at the corner well. The gels were overlaid with 1% (w/v) agarose (containing 0.1% (w/v) Bromophenol Blue tracking dye). Separation was performed at 13 • C with a constant current of 25 mA per gel in a Protean II xi cell (BioRad, Hercules, CA, USA) electrophoresis system. When tracking dye reached the base of the gel, the current was terminated. Gels were washed with deionized water and incubated overnight in staining solution containing 10% (v/v) acetic acid, 40% (v/v) methanol and 0.1% (w/v) Coomassie Brilliant Blue G-250. After staining, the gels were placed in destaining solution containing 7% (v/v) acetic acid and 25% (v/v) methanol until the background was clear and protein spots visible.
Image of gels were digitized in a densitometer (ImageScanner III, GE Healthcare Bio-Sciences, Uppsala, Sweden) for analysis based on spot density and location ( Figure S3). Image analysis was performed with PDQuest software version 8.0.1 (BioRad, Hercules, CA, USA). Spots were quantified on the basis of their relative volume which was determined by the ratio of the volume of a single spot to the whole set of spots under low P stress. Spots with significant (more than two-fold differential expression, α = 0.05 by Student's t-test) and reproducible changes in three replicates were used for further analysis.

Trypsin Digestion of Proteins
The protein spots of interest were picked from the gels, washed twice with MilliQ water, followed by 50% (v/v) acetonitrile. To this, 100 µL of 10 mM DTT dissolved in 25 mM ammonium bicarbonate (reducing buffer) was added and incubated at 45 • C for 60 min. Reducing buffer was discarded and protein spots were incubated in dark for 15 min with 100 µL of 10 mM iodoacetamide dissolved in 25 mM ammonium bicarbonate (alkylating buffer). This was followed by two sequential washing steps with 50% (v/v) and 100% (v/v) acetonitrile. The pellet was allowed to dry thoroughly and resuspended in 50 µL of 10 µg mL −1 trypsin dissolved in 25 mM ammonium bicarbonate and incubated overnight at 35 • C. The supernatant was transferred to a fresh Eppendorf tube and lyophilized at −80 • C. Lyophilized peptides were suspended in 1% (w/v) tetrafluoroacetic acid in 50% (v/v) acetonitrile prior to mass spectrometry.

Mass Spectrometry for Protein Identification
The peptide mass fingerprints of DAPs were collected on an AB Sciex TOF/TOF™ 5800 system with Series Explorer™ 7000 (AB Sciex, Concord, ON, Canada). The parameters were set as follows: Digestion enzyme: trypsin (specificity: C-terminal to Arginine and Lysine) with one missed cleavage; fixed modification: Carbamidomethyl (C); Mass Spectrometry (precursor-ion) peak filtering: 800-4000 m/z interval, monoisotopic, minimum signal-to-noise ratio (S/N) 10, mass tolerance 250 ppm; database used: Viridiplantae taxonomic sub-database (52,74,071 sequences) of the National Centre for Biotechnology Information protein database (NCBIprot, release date 08-11-2017; 12,86,24,863 sequences). Peptide mass fingerprints were searched against the database using the online Mascot server (version 2.6.0, Matrix Science Limited, London, UK). Proteins with Mascot score greater than 80 were considered as significant (P < 0.05) hits.

In-Silico Analysis for Protein Annotation
Protein sequences identified from the Mascot search were functionally annotated using the bioinformatics platform Blast2GO [30] by assigning their associated generic Gene Ontology (GO) terms and Enzyme codes. This was based on homology to proteins from other species as determined by BLAST and the occurrence of InterPro functional domains identified by InterProScan. Annotations were further expanded using ANNEX [31]. BLAST searches were conducted for each protein (TBLASTX, nr database, report 20 hits, maximum e-value 1E −10 ), followed by mapping and annotation. The resulting GO terms were mapped onto the corresponding Plant GO Slim terms. The annotated enzymes were also assigned to their respective KEGG pathway maps [32].

Validation of DAPs at Transcript Level by Reverse Transcription-qPCR
Total RNA was isolated from root tissue with TRIzol reagent (Invitrogen) and quantified in Nanodrop1000 Spectrophotometer (Thermo Scientific, Waltham, MA, USA). Contamination by genomic DNA was removed by treating 10 µg of RNA with DNase I (Promega). RNA integrity was verified in 1% (w/v) agarose gel stained with ethidium bromide and visualized under UV light in a gel documentation system (AlphaImager, Cell Biosciences, Heidelberg, Germany). cDNA was synthesized using SuperScript III reverse transcriptase (Invitrogen) according to the manufacturer's instructions. All samples were normalized to contain 50 ng·µL −1 of cDNA confirmed by the amplification of the reference gene "elongation factor 1α" (Gm_EF1α) and visualized on a 3% (w/v) agarose gel. Gene-specific primers (Table S1A,B) of proteins with significant Mascot scores were designed using the RealTime quantitative polymerase chain reaction (qPCR) tool of Integrated DNA Technologies [33]. Reverse transcription quantitative polymerase chain reaction (RT-qPCR) was carried out using the DyNAmoColorFlash SYBR Green I qPCR kit (Thermo Scientific, Waltham, MA, USA) on a Stratagene Mx3005P qPCR System (Agilent Technologies, Santa Clara, CA, USA). Reverse transcription-qPCR was performed in triplicate of 20 µL reaction containing 10 µL of 2X SYBR Green I master mix, 0.12 µL of 50X ROX passive reference dye, 1 µL of forward primer (0.5 pmol µL −1 ), 1 µL of reverse primer (0.5 pmol µL −1 ), 1 µL of cDNA template and 6.88 µL of nuclease-free water. The cycling parameters were as follows: 10 min of pre-denaturation at 95 • C, 40 cycles of 30 s at 95 • C and 30 s at 60 • C, followed by dissociation curve analysis (1 min at 95 • C, 30 s at 55 • C and ramp up to 95 • C). The MxPro software version 4.10 (Agilent Technologies, Santa Clara, CA, USA) was used for data collection. Melt curves were examined to detect inadvertent multiple amplicons. Primer specificity was also ensured by presence of single PCR product visualized on a 3% (w/v) agarose gel. The comparative cycle threshold method [34] was used to calculate relative transcript levels at experimental condition (low P) with three reference genes: 60S ribosomal subunit (Gm_60S), F-box (Gm_F-box) and GmEF1α. Negative controls were incorporated for each primer pair and individual PCR reactions were performed in triplicates.

Experimental Design and Statistical Rationale
All experiments were completely randomized with two factors: phosphorus level (P) and genotype (G). Physiological experiments to measure growth, root traits and root exudation were technically replicated thrice (n = 9), while protein isolation and 2D separation was carried out once in each biological replicate (n = 3). For validation of gene expression, RNA was isolated from three biological replicates, and RT-qPCR was carried out in three technical replicates (n = 9). Procedures for descriptive statistics and analysis of variance (α = 0.001 for physiological data; α = 0.05 for proteomics and expression validation) were carried out in the statistical software R version 3.1.2 (R Foundation for Statistical Computing, Vienna, Austria). Graphs and figures were plotted using GraphPad Prism version 6.00 (GraphPad Software, La Jolla, CA, USA).

Results
Soybean genotypes EC-232019 and EC-113396 responded differentially to low P stress, in terms of biomass accumulation, root system traits, uptake of P and carboxylate efflux as well as the root proteome profile. The salient results of the investigation are highlighted below.

Biomass Accumulation, Root System Traits and Tissue PosphorusStatus
P level, genotype and interactive effect of P × genotype on total biomass and root-to-shoot ratio was significant (P < 0.001) ( Table S2). Reduction in biomass was more prominent in EC-113396 compared to EC-232019 ( Figure 1A). Biomass partitioning to the roots was higher in EC-232019 compared to EC-113396 under both sufficient and low P level ( Figure 1B). P level had no significant effect on root-to-shoot ratio EC-232019, while there was a slight increase under low P in EC-113396. Significant genotypic variation was observed in root length, surface area and volume under low P stress (Table S2, Figure 1C-E). Low P stress increased root surface area by 106% and root volume by 85% in the EC-232019 whereas both traits were reduced by 50% in EC-113396 at low P compared to sufficient P. P level significantly (P < 0.001) influenced tissue P status, with an overall reduction in shoot (63%) and root (73%) P concentrations at low P stress compared to sufficient P (Table S2, Figure 1F,G) in both genotypes. Irrespective of genotype, total P uptake was 70% lower at low P as compared to sufficient P ( Figure 1H), while EC-232019 exhibited least reduction (19%) in total P uptake at low P compared to sufficient P. Foundation for Statistical Computing, Vienna, Austria). Graphs and figures were plotted using GraphPad Prism version 6.00 (GraphPad Software, La Jolla, CA, USA).

Results
Soybean genotypes EC-232019 and EC-113396 responded differentially to low P stress, in terms of biomass accumulation, root system traits, uptake of P and carboxylate efflux as well as the root proteome profile. The salient results of the investigation are highlighted below.

Biomass Accumulation, Root System Traits and Tissue PosphorusStatus
P level, genotype and interactive effect of P × genotype on total biomass and root-to-shoot ratio was significant (P<0.001) ( Table S2). Reduction in biomass was more prominent in EC-113396 compared to EC-232019 ( Figure 1A). Biomass partitioning to the roots was higher in EC-232019 compared to EC-113396 under both sufficient and low P level ( Figure 1B). P level had no significant effect on root-to-shoot ratio EC-232019, while there was a slight increase under low P in EC-113396. Significant genotypic variation was observed in root length, surface area and volume under low P stress (Table S2, Figure 1C-E). Low P stress increased root surface area by 106% and root volume by 85% in the EC-232019 whereas both traits were reduced by 50% in EC-113396 at low P compared to sufficient P. P level significantly (P < 0.001) influenced tissue P status, with an overall reduction in shoot (63%) and root (73%) P concentrations at low P stress compared to sufficient P (Table S2, Figure  1F,G) in both genotypes. Irrespective of genotype, total P uptake was 70% lower at low P as compared to sufficient P ( Figure 1H), while EC-232019 exhibited least reduction (19%) in total P uptake at low P compared to sufficient P.

Carboxylate Efflux in Response to Low Posphorus Stress
Low P stress-induced carboxylate efflux differed both in type and quantity, and exhibited a significant (P < 0.001) variability among contrasting soybean genotypes (Table S2). Averaged over P levels, EC-232019 exhibited higher rate of total carboxylate efflux compared to EC-113396. Low P induced-total carboxylate exudation increased by 58% in EC-232019, while it was reduced by 35% in EC-113396. In EC-232019, the efflux of oxalate, citrate, fumarate, succinate and lactate was induced under low P ( Figure 2). Further, succinate and fumarate concentration in the root exudate increased by two-fold at low P in comparison to sufficient P in EC-232019.

Carboxylate Efflux in Response to Low Posphorus Stress
Low P stress-induced carboxylate efflux differed both in type and quantity, and exhibited a significant (P < 0.001) variability among contrasting soybean genotypes (Table S2). Averaged over P levels, EC-232019 exhibited higher rate of total carboxylate efflux compared to EC-113396. Low P induced-total carboxylate exudation increased by 58% in EC-232019, while it was reduced by 35% in EC-113396. In EC-232019, the efflux of oxalate, citrate, fumarate, succinate and lactate was induced under low P ( Figure 2). Further, succinate and fumarate concentration in the root exudate increased by two-fold at low P in comparison to sufficient P in EC-232019. and low (4 µM) P. Data correspond to mean±SE (n = 9). "nd" denotes peak "not detectable". RFW: Root fresh weight. *, ** and *** denote significance at 0.05, 0.01 and 0.001 probability levels, respectively.

Comparative Analysis of Soybean Root Proteome
Staining of two-dimensional electrophoretic gels using Coomassie Brilliant Blue dye revealed a total of 325 protein spots in the root tissue of EC-232019 and EC-113396 at low and sufficient P ( Figure S3). Out of these, 105 (32%) DAPs were observed between sufficient and low P levels ( Figure 3A). A total of 44 (14%) proteins decreased by more than two-fold under low P stress, while 61 (19%) proteins increased by more than two-fold at low P. These 105 DAP spots were picked from the gels, digested with trypsin and sequenced by matrix-assisted laser desorption/ionization (MALDI) to obtain the peptide mass fingerprints. Annotated spectra of the DAP spots are provided in Datasets S1 and S2.

Differentially Abundant Proteins in Response to Low Posphorus Stress
Functional annotation of the DAPs based on gene ontology revealed their involvement in a myriad of biological processes including biosynthetic pathways, generation of precursor metabolites and energy, carbohydrate, protein and lipid metabolism ( Figure 3B), localized to different cellular components ( Figure 3C). Out of 61 increased proteins at low P, 27 were specific to EC-232019, 16 increased only in EC-113396 and 18 were common to both genotypes (Table 1). Eighteen peptide sequences with significant Mascot scores were derived from taxonomic databases of soybean (G. max) or its wild relative (G. soja). Some of the proteins significantly increased at low P condition were phosphoglycerate mutase, malate dehydrogenase, fructokinase, phosphoglucomutase, cysteine synthase, actin, 70 kDa heat shock-related protein, heat shock cognate protein, proteasome, ATP synthase, isoflavone reductase and monodehydroascorbate reductase ( Figure 4). Among the proteins with decreased at low P, six were common to both genotypes, while 24 and 14 were specific to EC-232019 and EC-113396, respectively (Table 2). Twenty-three peptide sequences with significant Mascot scores were soybean (G. max) homologs. Twenty-two proteins were predicted from plant genome sequences, six were hypothetical and six uncharacterized or unknown proteins. Proteins decreased at low P included triosephosphate isomerase, phosphogluconate dehydrogenase, enolase, methionine synthase, chalcone isomerase, isocitrate dehydrogenase, glutathione-S-transferase, heat shock protein 70 kDa and alcohol dehydrogenase ( Figure 5).       Table 1 for entire list of increased proteins.  Table 2 for entire list of decreased proteins.

Validation of Expression of Genes Encoding DAPs in Response to Low P Stress
To validate the expression levels of genes encoding DAPs, RT-qPCR was performed on 44 out of the 105 proteins that showed differential expression under low P stress. In EC-232019, the genes encoding isoflavone reductase, actin-7-like, malate dehydrogenase and monodehydroascorbate reductase were up-regulated by more than two-fold at low P compared to sufficient P ( Figure 6A, Table S3A). Transcripts of phosphoglucomutase, phosphoglycerate mutase, heat shock cognate protein 80-like, trypsin inhibitor, cysteine synthase and an uncharacterized protein showed more than two-fold expression at low P in EC-113396. The gene encoding ATP synthase β subunit was up-regulated by more than two-fold in both the genotypes. Contrary to protein expression, transcripts of the genes encoding fructokinase and cytosolic glutamine synthetase β 2 were down-regulated at low P compared to sufficient P in both genotypes. Similarly, 19 genes were down-regulated at low P compared to sufficient P in agreement with protein expression pattern, whereas contrary to the proteome profile, transcripts of the genes encoding methionine synthase and glutathione S-transferase were up-regulated at low P compared to sufficient P in EC-232019 ( Figure 6B, Table S3B). In summary, 36 genes displayed trends at the transcript level that were consistent with protein expression, while eight genes exhibited transcription trends that were opposite to that observed in the proteome profile. Figure 5. Root proteins decreased by more than two-fold at low (4 µM) P in contrasting soybean genotypes. Arrows denote the protein spots zoomed for better visualization. Bars in the graph denote the mean (n=3) normalized intensity of protein spots on gels ofEC-232019 Sufficient P (orange, SP ER), EC-232019 Low P (red, LP ER), EC-113396 Sufficient P (violet, SP IENR) and EC-113396 Low P (green, LP IENR). Proteins with significant Mascot score are presented here. Refer Table 2 for entire list of decreased proteins.

Validation of Expression of Genes Encoding DAPs in Response to Low P Stress
To validate the expression levels of genes encoding DAPs, RT-qPCR was performed on 44 out of the 105 proteins that showed differential expression under low P stress. In EC-232019, the genes encoding isoflavone reductase, actin-7-like, malate dehydrogenase and monodehydroascorbate reductase were up-regulated by more than two-fold at low P compared to sufficient P ( Figure 6A, Table S3A). Transcripts of phosphoglucomutase, phosphoglycerate mutase, heat shock cognate protein 80-like, trypsin inhibitor, cysteine synthase and an uncharacterized protein showed more than two-fold expression at low P in EC-113396. The gene encoding ATP synthase β subunit was upregulated by more than two-fold in both the genotypes. Contrary to protein expression, transcripts of the genes encoding fructokinase and cytosolic glutamine synthetase β2 were down-regulated at low P compared to sufficient P in both genotypes. Similarly, 19 genes were down-regulated at low P compared to sufficient P in agreement with protein expression pattern, whereas contrary to the proteome profile, transcripts of the genes encoding methionine synthase and glutathione S-transferase were up-regulated at low P compared to sufficient P in EC-232019 ( Figure 6B, Table  S3B). In summary, 36 genes displayed trends at the transcript level that were consistent with protein expression, while eight genes exhibited transcription trends that were opposite to that observed in the proteome profile.    Figure S3B,D. Superscripts a and b denote increased proteins at low P compared to sufficient P in EC-232019 and EC-113396, respectively. E denotes the expectation value, M r and pI are monoisotopic mass and calculated isoelectric point, respectively. #P and C denote the number of matched/unmatched peptides and protein sequence coverage, respectively. Proteins with Mascot score > 80 are significant (P < 0.05).  Spot ID corresponds to the protein spots in Supplemental Figure S3A,C. Superscripts a and b denote decreased proteins at low P compared to sufficient P in EC-232019 and EC-113396, respectively. E denotes the expectation value, M r and pI are monoisotopic mass and calculated isoelectric point, respectively. #P and C denote the number of matched/unmatched peptides and protein sequence coverage, respectively. Proteins with Mascot score > 80 are significant (P < 0.05).

Soybean Genotypes Exhibit Improved Root System Traitsand Carboxylate Efflux under Low P Stress
Phosphorus nutrition had significant effect on the root system traits of both soybean genotypes, albeit the increase in root surface area and volume was augmented in the P-efficient EC-232019. Similar alterations to root morphology in response to low P stress was reported in lentil (Lens culinaris) [35], rapeseed [36], maize [37], green gram [10] and wheat [38]. P-efficient soybean genotype with improved root traits exhibited higher total P uptake under low P stress. The improved root system traits exhibited by the P-efficient soybean genotype EC-232019 attributed to the least reduction in total P uptake at low P compared to sufficient P. Higher root surface area improved P acquisition by increasing the amount of root exudates such as phosphatases, RNases, nucleases, apyrases and carboxylic acids [39]. Our results conform to responses of barley (Hordeum vulgare) to low P stress, wherein improved P uptake of efficient cultivars was attributed to enhanced root exudation that increased its ability to acquire more P [40]. P-efficient maize genotypes also exhibited higher carboxylic acid efflux, P uptake and biomass under low soil P availability in comparison to inefficient ones [41]. P-efficient soybean genotype EC-232019 exuded greater amounts of carboxylic acids under low P stress. Irrespective of P level or genotype, carboxylic acids in the root exudate comprised of fumarate >oxalate >lactate >pyruvate >succinate >malate. Other workers reported root exudates comprising fumarate, citrate and malate in P-stressed alfalfa [6] and soybean [7], which supported our findings. Oxalate and malate were the major carboxylic acids detected in P-efficient soybean [8,42]. Carboxylates mobilize P bound to metal ligands in the soil, hence their functionality depends on the number and arrangement of carboxyl and hydroxyl moieties. The complexing capacity for Al/Fe/Ca follows the decreasing order of tri-(citrate 3− )> di-(malate 2− , oxalate 2− , fumarate 2− , succinate 2 )> mono-carboxylic acids (lactate 1− ) [5]. Thus, soybean genotypes with higher root exudation potential efficiently maintain their tissue P status to sustain growth under low P stress condition.

Soybean Genotypes Exhibit Differential Molecular Regulation under Low Posphorus Stress
Results of comparative proteome analysis and validation by RT-qPCR have been consolidated in Figure 7, which illustrates the differential molecular regulation in soybean genotypes with contrasting carboxylic acid synthesis and efflux under low P stress; suggesting the crosstalk between various metabolic pathways implicated in conferring superior P acquisition efficiency under stress.

Tricarboxylic Acid Cycle and Glycolysis
Low P stress induced the expression of malate dehydrogenase in EC-232019 but not in EC-113396. Such a response at the molecular level might be attributed to increased synthesis and efflux of malate ions as observed in EC-232019 at low P. Further, abundance of isocitrate dehydrogenase decreased at low P (relative decline higher in EC-232019 compared to EC-113396), possibly leading to higher accumulation of malate in the root tissues [43]. Citrate synthase activity was higher at low P in EC-232019 (data published in Vengavasi et al. [26]), which might have increased the synthesis of citrate required for induction of efflux in the efficient genotype. Similarly, increase in the activity of phosphoenolpyruvate carboxylase, one of the key regulatory enzymes for replenishment of carbon was higher in EC-232019 compared to EC-113396. Other enzymes in the glycolytic cycle including pyruvate kinase, enolase, phosphoglycerate kinase and triosephosphate isomerase decreased in EC-232019 under low P stress. Such a response in the P-efficient soybean genotype might indicate inorganic P recycling through alternative glycolytic by-pass reactions (dotted lines in Figure 7) catalyzed by nucleoside diphosphate kinase, non-phosphorylating NADP-glyceraldehyde-3-phosphate dehydrogenase, pyruvate phosphate dikinase, malic enzyme and phosphoenolpyruvate carboxylase [44]. In concurrence, abundance of nucleoside diphosphate kinase increased at low P in EC-232019.

Starch Hydrolysis
Glucan water dikinase, an important enzyme involved in starch hydrolysis was increased by low P stress in EC-232019 but not in EC-113396, thereby increasing the flux of carbon through glycolysis and tricarboxylic acid cycle (TCA) cycle.

Anaerobic Respiration
Alcohol dehydrogenase, the major enzyme producing ethanol from pyruvate decreased under low P stress in both genotypes, indicating the possible production of lactate ions, as evident from its higher efflux in soybean ( Figure 2G).

Other Anaplerotic Reactions Replenishing the TCA Cycle Intermediates
The increased expression of root cytosolic glutamine synthetase under low P in EC-232019 and not in EC-113396 suggests its role in stress tolerance, which is in agreement with results obtained in creeping bentgrass (Agrostis stolonifera) [45]. Reactions occurring in the glutamine synthetase-glutamine:2-oxoglutarate amidotransferase GS-GOGAT cycle and the glutamate dehydrogenase (GDH) shunt possibly siphoned carbon in amino acids back into the TCA cycle [46]. Argininosuccinate lyase that catalyzes the irreversible breakdown of argininosuccinate to arginine and fumarate increased under low P stress in EC-232019 and not in EC-113396. Fumarate synthesized in this reaction might be available for higher efflux at low P (evident from Figure 2D) in addition to that produced in the TCA cycle. Arginiosuccinate lyase is also involved in regulating root elongation as reported in rice [47]. Similarly, abundance of fructokinase, a major contributor to glycolytic carbon flux during root growth and development [48] also increased under low P in EC-113396.

Synthesis of Sulfur Containing Amino Acids
In addition to regulation of carboxylic acid synthesis, low P stress in soybean influences several other metabolic pathways. Sustenance of sulfur assimilatory pathway under low P stress is important for recycling of P i from plastidic phospholipids [49,50]. Increased abundance of the enzymes cysteine synthase, ATP sulfurylase and S-adenosyl homocysteine synthase is augmented in EC-232019, indicating sustained sulfur assimilation reactions even under low P stress. Methionine synthase, which is also involved in increasing the carbon source under stress as reported in creeping bentgrass [45], increased at low P in EC-113396 but not in EC-232019.

Other Pathways Affected by Low P Stress
Abundance of proteins involved in maintenance of cellular homeostasis such as monodehydroascorbate reductase, 80 kDa heat shock cognate protein, 70 kDa heat shock related protein and proteasome subunit alpha and isoflavone reductase increased in EC-232019 under low P stress. Similar responses were observed in isopentenyl transferase (ipt) overexpressed creeping bentgrass tolerant to drought stress [45]. Isoflavone reductase-like gene (OsIRL) was found to be involved in scavenging reactive oxygen species in rice [51]. The abundance of actin increased at low P only in EC-232019, while tubulin increased in both genotypes under low P stress. Abundance of actin and tubulin proteins is reported under drought stress, with important roles in determining root cell structure [45]. Subunits of the energy producing protein ATP synthase showed variable response, being increased and decreased under low P in both genotypes.

•
Differentially abundant proteins with known physiological function may be tested rigorously for imparting P acquisition efficiency by creating overexpression/knockout mutant lines in soybean or other model systems.

•
Hypothetical proteins with putative or unknown function identified in the proteomic analysis may be functionally characterized to ascertain their role(s) under low P stress.

•
The identified genotypes may be potential donors in crop improvement programs to develop high-yielding P-efficient cultivars, an asset to low-input sustainable agriculture.

Conclusions
The proportion of carboxylates (including oxalate, citrate, succinate and fumarate) was highest among root exudates of P-efficient soybean EC-232019. Sustained growth of EC-232019 under low P stress may be attributed to improved root morphological traits and efflux of carboxylates. Key enzymes in the tricarboxylic acid cycle and glycolytic pathways that were induced or suppressed at low P stress showed clear-cut differential regulation among the contrasting genotypes. Enzymes including but not limited to malate dehydrogenase, isocitrate dehydrogenase, phosphoenolpyruvate carboxylase, citrate synthase, glutamine synthetase, argininosuccinate lyase and alcohol dehydrogenase might be implicated in contributing to the additional carbon required for the carboxylate synthesis and efflux in EC-232019 under low P stress. Further, alterations at the transcript and protein level suggest the interdependence of various metabolic pathways conferring higher P acquisition efficiency to plants under stress.
Supplementary Materials: The following are available online at www.mdpi.com/2073-4425/8/12/341/s1. Figure S1: (A) Digitized root images in two contrasting soybean genotypes grown at sufficient (250 µM) and low (4 µM) P, (B) soybean plants immersed in trap solution for root exudate collection, Figure S2: SDS-PAGE gels of total root protein in contrasting soybean genotypes grown at sufficient (250 µM) and low (4 µM) P, Figure S3: Two-dimensional gels of total root protein in contrasting soybean genotypes grown at sufficient (250 µM) and low (4 µM) P, Table S1A: Primers used for RT-qPCR to validate transcript levels of increased proteins at low P (4 µM), Table S1B:Primers used for RT-qPCR to validate transcript levels of decreased proteins at low P (4 µM), Table S2: Significance of F values derived from analysis of variance for various physiological traits in contrasting soybean genotypes, Table S3A: Significance of F values derived from analysis of variance of transcript levels of proteins increased at low P (4 µM), Table S3B: Significance of F values derived from analysis of variance of transcript levels of proteins decreased at low P (4 µM), Dataset S1: Annotated spectra of increased proteins at low P (4 µM) in soybean, Dataset S2: Annotated spectra of decreased proteins at low P (4 µM) in soybean.