INTRODUCTION

Water scarcity results in huge reductions in crop yield and is one of the greatest limitations on the expansion of agriculture areas. Although the total annual precipitation in Brazil is sufficient for soybean cultivation, water deficiency could be caused by a dry spell of many weeks, without any rain in the growth period in a particular location. As a result, soybean yields could be affected by a water deficit, particularly during flowering and early pod expansion [1]. Drought stress induces several changes at morphological, physiological and biochemical levels in all plant organs. On the other hand, plants have developed several strategies in dealing with drought stress, including a short life cycle, enhanced water uptake and reduced water loss, as well as osmotic adjustment, antioxidant capacity, and desiccation tolerance [2].

Metabolic adjustments in response to the environmental conditions may alter pools of metabolites that play important roles in the physiological response. The study of these metabolites may indicate which pathways have been disturbed by the stress [3–5] and could be used to indicate the cell adjustment during drought [6]. Some of the most important responses of plants against drought stress are related to the accumulation of osmoprotectant metabolites, which contribute to water retention in stressed cells [7]. These osmoprotectants maintain relatively low cellular osmotic potential and may greatly differ in their chemical features and concentrations among plant species.

Several genes respond to water stresses at pre- and post-transcriptional and translational levels [8,9]. Understanding the role of stress-induced proteins is important in explaining the plant tolerance processes. Therefore, the study of proteomics is an approach that could be used to discover proteins and pathways associated with crop physiological and stress responses [10]. Analyses of the soybean leaf proteome have been used to identify the differential expression of various abiotic stress-responsive proteins that regulate many molecular processes and signaling cascades [9,11,12]. On the other hand, the metabolomic profiling shows the metabolites involved in the tolerance to drought and heat stress [5,13]. Thus, studying metabolome and proteome of the plants could help in understanding the pathways involved in water stress tolerance.

Brazil is the second leading producer of soybean in the world and has the potential to become the largest due to increases in cultivated areas in recent years (OCDE-FAO, 2017). However, the effects of drought which is the major challenge for soybean production in Brazilian savannas must be reduced. Therefore, the development of soybean cultivars which are drought-tolerant is required in order to improve their productivity and meet the ever-increasing demand for soybean grains worldwide. Therefore, it is important to understand the mechanism of drought tolerance in soybean plants.

It has already been reported that soybean plants can adapt to these conditions by triggering protective mechanisms that enable them to survive drought conditions [14,15]. The development of two genotypes; Brazilian soybeans (Glycine max L. Merrill), Embrapa 48 and BR 16, were evaluated under drought conditions in a greenhouse [16] and in the field [17]. The results showed that the Embrapa 48 cultivar may be characterized as drought-tolerant and the BR 16 cultivar as drought-sensitive. Physiological responses and gene expression profiles were also evaluated to understand the mechanism of tolerance for these cultivars [17,18]. Until now, no studies have investigated proteome and metabolome, so this research can lead to a new understanding of the contrasting molecular mechanisms for drought tolerance between these genotypes that share a common ancestor.

Therefore, in this work, we performed a combination of biochemical, proteomics and metabolomics approaches to investigate the molecular components that have shown greater variations when compared to the responsive profiles of the Embrapa 48 and BR 16 cultivars. To achieve this goal, we focused on metabolites and proteins responsive to stress conditions that directly affect the drought tolerance between genotypes. Our results showed an agreement between the metabolic response and the physiological status, that allowed for higher rates of leaf growth, which was not caused by osmoprotection by amino acids biosynthesis or by ABA-dependent cascades. Proteomic and metabolomics data indicate that the metabolism of the Embrapa 48 tolerant genotype is not impaired, when compared to the sensitive BR 16. Thus, regulatory cascades of the metabolic activities, mediated by protein phosphorylation, may be operating for a more efficient water use.

MATERIALS AND METHODS

The soybean genotypes BR 16 and Embrapa 48 which are sensitive and tolerant to drought, respectively, were selected for a study based on previous experiments [16,17]. We prepared 6 pots containing 10 L of a mixture of soil, sand, and dung (3:1:1). 3 pots were used to sow each of the soybean cultivars BR 16 and Embrapa 48, using 3 seeds in each pot. The pots were taken to a greenhouse and kept under natural light and photoperiod, with relative humidity varying daily between 65% and 85% and a temperature varying between 15 °C and 35 °C. The seedlings of both cultivars were watered once a day with 100 mL of water, for 43 days, until the pre-flowering period (stage V4), when the water supply was interrupted to set the experimental treatments.

The experimental design was completely randomized, in a 3 × 2 factorial design. The first factor was the plants’ water potential (−1.5 MPa, −1.0 MPa and control) and the second was the two different cultivars. The water regimes were defined as irrigated (IR), related to the control treatment, and the non-irrigated treatment (NI Ψ1.0 and NI Ψ1.5). We used a Scholander pressure pump to measure the water potential of the plants and set the different treatments. In the morning (Ψam) (4:30–5:30 h) we monitored the water pressure in the fourth leaf from the apical meristem of each plant, according to DaMatta [19]. We cut off the individual leaf when the water potential of the plants reached −1 MPa and −1,5 MPa, and then stored them in liquid nitrogen at −80 °C. For each treatment 3 pots were used and each pot contained 3 plants. A trifoliate leaf from each plant from one pot was collected together (3 plants by replicate). Thus, the biochemical analyses were performed using 3 distinct pools, resulting in 3 biological replicates. These procedures were performed on both plant cultivars.

The protein extraction of the soybean leaves was performed on the plants under drought stress (−1 MPa) and on the irrigated controls, according to the method for fractionation/precipitation of abundant proteins, described by Aryal [20] with some modifications. All the washing and centrifugation steps were performed using 50 mL centrifuge tubes or 15 mL Falcon tubes. The soybean leaves were grounded with liquid nitrogen in pre-cooled porcelain mortars. The powdered tissue (approximately 3 g) was suspended in a 10 mL extraction buffer [0.5 M Tris-HCl (pH 8.0), 20 mM MgCl₂, 2% Triton X-100, 2% β-mercaptoethanol, 0.2% protease inhibitor cocktail (P-9599, Sigma Aldrich, San Luiz, USA), 1% phosphatase inhibitor cocktail (P-0044, Sigma Aldrich, San Luiz, USA), 25 mM Sodium Fluoride (NaF), 1 mM Sodium Molybdate (Na₂MoO₄) and 0.5% polyvinylpolypyrrolidone (PVPP)]. After homogenization in a vortex for 30 s, the solution was filtered using four layers of gauze and the filtrate was transferred into a new 50 mL tube. A 10 mL extraction buffer containing 30% polyethylene glycol 4000 (PEG 4000, Sigma-Aldrich, San Luiz, USA) was added to the filtrate. The solution was homogenized again using a vortex for 30 s and kept on ice for 30 min, followed by centrifugation at 3200× g for 20 min at 4 °C. After the two extraction procedures, the supernatants and the pellet (designed as PEG-pellet) were collected for further analysis.

The supernatant, containing low abundance proteins, was transferred into a new 50 mL tube. The protein precipitation was carried out by adding 20 mL of 20% trichloroacetic acid (20% TCA) to the tube, and kept overnight at −20 °C before centrifuging the solution at 6000× g for 30 min at 4 °C. The supernatant was discarded, and the pellet was then collected, dried at room temperature, followed by a solubilization in 10 mL of dense SDS-buffer (30% sucrose, 2% SDS, 0.1 M Tris-HCl pH 8.0 and 2% β-mercaptoethanol). The next step was to add 5 mL of buffered phenol to the solution, homogenized using a vortex for 30 s and transferred into a 15 mL tube. The solution was kept on ice for 30 min, while it was homogenized using a vortex 4 times, and then centrifuged at 6000× g for 30 min at 4 °C. The phenolic phase was collected and transferred into a 50 mL tube, followed by adding 5 volumes of 0.1 M ammonium acetate in methanol in order to precipitate the proteins. The tube was kept overnight at −20 °C. After that, the solution was centrifuged at 6000× g for 30 min at 4 °C. The supernatant was discarded and the pellet was washed performing 3 cycles of resuspension and centrifugation at 6000× g for 10 min at 4 °C using 80% acetone (ice cold) and 1 cycle using 70% ethanol. After being washed, the pellet was dried at room temperature. The protein extract was solubilized in a solution containing 7 M urea, 2 M thiourea, and 2.0% CHAPS. The resulting solution was sonicated using an UltraSonic Processor (Model GE 50) for cycles of 5 to 10 s until complete solubilization of the pellet. The sample was then stored at −80 ºC.

The pellet-PEG, containing high abundance of proteins, was solubilized in SDS-dense buffer (30% sucrose, 2% SDS, 0.1M Tris-HCl pH 8.0 and 2% β-mercaptoethanol) at room temperature. After the solubilization, 5 mL of buffered phenol was added, and the solution was homogenized using a vortex for 30 s and transferred to a 15 mL tube. The 15 mL tube was kept on ice for 30 min, using a vortex every 4 min and centrifuged at 6000× g for 30 min at 4 °C. The phenolic phase obtained after the centrifugation was collected and transferred into a 50 mL tube, and 5 volumes of 0.1 M ammonium acetate in methanol were added to precipitate the proteins. The tube was kept overnight at −20 °C, before centrifuging the solution at 6000× g for 30 min at 4 °C. The supernatant was then discarded and the pellet obtained was collected, washed (3 cycles of resuspension and centrifugation at 6000× g for 10 min at 4 °C) with 80% acetone (ice cold) and 70% ethanol before being dried at room temperature. The protein extracted was solubilized with 7 M urea, 2 M thiourea, and 2.0% CHAPS. The solution was sonicated using an UltraSonic Processor (Model GE 50, Cole-Parmer, Vernon Hills, USA) for cycles of 5 to 10 seconds until complete solubilization of the pellet. The sample was then stored at −80 °C.

The isoelectric focusing (IEF) was performed in a 13 cm IPG strip pH 3–10 (GE Healthcare, Chicago, USA) as described by Coutinho [21]. SDS-PAGE electrophoresis was carried out at 20 °C in the SE 600 Ruby system (GE Healthcare) at 40 mA/gel for 5 h. The gels were stained with Coomassie blue G-250 solution (8% ammonium sulfate (v/v), 0.8% phosphoric acid (v/v), 0.08% Coomassie Blue G-250 (v/v) and 20% methanol (v/v)), scanned using an ImageScanner III (GE Healthcare, Chicago, USA), and calibrated using the Labscan software (GE Healthcare, Chicago, USA). To correctly determine the protein spots with differential abundance, a comparison of the protein profiles between the genotypes submitted to irrigated and non-irrigated treatments were performed using the ImageMaster2D Platinum 7 software (GE Healthcare, Chicago, USA). For the protein spot quantitative analysis, only the spots present in three gels were considered as differentially expressed. The spots greater than 1.5, indicated differences between the protein expressions of the classes. One-way ANOVA was performed considering a p-value less than 0.05. The gels were performed in triplicates.

The abundances (%V) of the protein spots expressed differentially using the criteria above were converted to fold-change and classified as up- or down-regulated (fold-change higher than 1.3 were plotted).

The detection of the phosphorylated proteins was carried out by staining the 2DE gels using a Pro-Q^® Diamond Phosphoprotein Gel Stain kit (Invitrogen, Calrsbad, USA), according to the protocol described by Agrawal and Thelen [22]. The gels were then scanned with a FLA 5100 laser scanner and the images were analyzed as described by Vital [23].

The spots stained with Pro-Q Diamond with ratio values greater than 1.5 were considered as differentially expressed (p < 0.05) phosphorylated proteins. The three replicates of phosphoprotein spots were selected to be identified using mass spectrometry.

Protein spots were excised and destained in 50 mM ammonium bicarbonate and 50% methanol (v/v), followed by acetonitrile addition for dehydration. The proteins were then reduced by adding 200 mM DTT solution in 100 mM of ammonium bicarbonate for 30 min of incubation at 56 °C, followed by alkylation with 200 mM of iodoacetamide solution containing 100 mM of ammonium bicarbonate for 30 min at room temperature. The spots were washed with 100 mM of ammonium bicarbonate solution, dehydrated with acetonitrile and dried using vacuum centrifugation. the spots were then rehydrated using trypsin digestion solution overnight (20 h) at 37 °C. The digested peptides were extracted using an extraction buffer (50% acetonitrile, 5% formic acids) and dried using vacuum centrifugation. The peptides were dissolved in 0.1% formic acid and analyzed online using LC-MS coupled to a UPLC system (nanoAcquity, Waters, Milford, USA), equipped with a column and a capillary trap C18 BEH130 1.7 uM 100 nm × 100 mm operating at a flow rate of 0.5 µL/min. The peptides eluted were automatically injected into a mass spectrometer Ion Trap (Amazon ETD, Bruker, Billerica, USA), using a nanoESI ionization needle. The scanning of the ions from the mass spectrometer was carried out between 300 and 1500 m/z in positive mode and the data were acquired for 70 min through LC-MS/MS analysis. The mass spectrometer was operated in the auto-MSn mode. The acquisition of data from the LC-MS instrument was controlled using Hystar software (Bruker, Billerica, USA) and the spectra were processed with the Data Analysis software (Bruker, Billerica, USA) using the default settings for proteomics. A peak list was converted to the mascot generic format (.mgf). Peak lists were used for protein identification MASCOT algorithm with a local Client license connected to a remote server. The parameters used in the program were: protein list obtained from Phytozome containing all the proteins described for soybeans in the database, methionine oxidation as a variable modification of cysteine, cysteine carbamidomethylation as a fixed modification, a missed cleavage, peptide charge state 2+, 3+, 4+, trypsin-like cleavage enzyme and a mass error of 0.1 Da. SCAFFOLD 4.2.1 version was used to validate MS/MS based peptide and the protein identifications. Peptide identification was accepted with a probability greater than 95% using three unique peptides.

The FASTA sequences of the proteins stained as phosphorylated were analyzed for possible phosphorylation sites. The prediction of the number of phosphorylation sites was performed online through the MUSITE website (http://www.musite.net), designed for plant protein, allowing for variation in the probability of occurrence of phosphorylation sites. For this approach, the general phosphorylation on serine and threonine were chosen, and all the green plants were considered for group analysis. To identify the phosphorylation sites, the algorithm takes the downstream and upstream amino acid sequence into the serine and threonine residues to determine if, and in which position, these amino acids will be phosphorylated. As there are not a large number of phosphorylated tyrosine residues in plants, the site does not identify the phosphorylation sites of this amino acid residue. In order to have differentiated levels of requirement, the number of phosphorylation sites was determined to be 70% (low requirement) and 95%.

The presence of homologous proteins or phosphorylated sites were confirmed using the P3DB site (Plant Protein Phosphorylation DataBase: http://www.p3db.org). In this site, the input of the search is the name of the proteins and the output is the details of the results as phosphoproteins, such as the phosphorylated amino acid sites, the position, the protein ontology, the plant species and the description of the experimental conditions. Therefore, it was possible to verify the existence or not of the phosphorylation sites from the identified proteins in our assays. It is worth noting that the non-identification of phosphorylation sites of a protein by the P3DB site does not necessarily imply that this protein is not phosphorylated because the information on this protein is not found on the database.

The metabolite extractions and derivatization were performed as described by Coutinho [21]. The metabolic extract was analyzed using an Agilent 7890A GC System coupled with a Mass Spectrometric TruTOF^® HT TOFMS (Leco, Setor Joseph, USA) equipped with a capillary column of 30 mm (MDN-35) operating according to Lisec [24]. The samples were injected in splitless mode at 230 °C using a gas flow (continuous flow of helium) of 2 mL/min. The oven temperature was initially maintained at 80 °C and then increased to 15 °C/min, until reaching 330 °C, before being kept at this temperature for 5 min. Mass spectrum was obtained by the full-scan method with a range from 33 to 600 m/z. A series of n-alkanes was used along with samples to calculate the retention indices.

The raw GC/MS data were processed and converted to the CDF format (NetCDF) using the ChromaTOF package and analyzed using TargetSearch algorithms [25]. A script was used to identify and quantify the metabolites designed to run on the R package [25]. The processing parameters and alignment used was optimized for our GC/TOF platform. The parameters used were: mass range of 85–500 Da; threshold 50; TopMasses 10; r threshold of 0.05; quality index threshold 600. The identification of the compounds was determined by database searches based on the combination of the mass spectrum and the chromatographic retention indices according to the TAG-based method [25]. In this study, a metabolites fragmentation library was characterized by electron impact (EI) and retention index (RI) GMDB_FAMELib_TS_20110228_IS.txt, generated by the GMD Mass Spectrum Reference Library (http://gmd.mpimp-golm.mpg.de/download/). It was generated using the same configuration of our GC/MS platform. A filtered identification table containing the identified compounds and their intensities were normalized by the fresh mass from the leaves, and used as input data for processing and statistical analysis using the MetaboAnalyst platform (http://www.metaboanalyst.ca/).

Tables containing the identified metabolites by GC/MS and phytohormonal concentrations were used as input data for the characterization of metabolites with differential relative abundances. A display of the metabolic pathway maps significantly disturbed using the MetaboAnalyst 3.0 package was generated. The Quality Filters based on the standard deviation methods were used to automatically remove the low-quality data. The intensity values were then normalized. The data were analyzed using the Partial Least Squares Discriminant Analysis (PLS-DA) for grouping the treatments, and to arrange the metabolites (VIP score) in order of their importance in the grouping of the samples. The significantly different metabolites (One-way ANOVA, p < 0.05) were analyzed using the module Pathway Analysis search tool present in the MetaboAnalyst (v 3.0) platform. To do this, the intensity of each identified metabolite showing differences in relative abundance was converted to fold-change and classified as up- and down dysregulated (fold-changes higher than 2.0 were used as input). A pathway library of Arabidopsis thaliana was used.

Ten milligrams of powdered tissue (~110 mg fresh weight) was mixed with 400 μL of methanol, isopropanol and acetic acid (20:79:1). The phytohormones were extracted and analyzed as described by Coutinho [21]. The sample was scanned by MRM (multiple reaction monitoring) to detect each hormone as shown in Supplementary Table S1. Three biological replications were performed. The generated mass spectra were processed using the MassHunter software to obtain the extracted chromatograms (XIC) of each transition and the areas to quantify each hormone. A standard curve of each of the hormones varying concentrations from 0.1 to 300 ng/mL was used to convert the area values from XIC in ng/g of plant tissue. The synthetic plant hormone naphthalene acetic acid (NAA) was added into the extraction buffer and used as an internal normalizer in the quantitative analysis.

RESULTS

In the control plants of both soybean cultivars (BR 16–sensitive and Embrapa 48–tolerant) the water potential was higher than −0.25 MPa. The tolerant cultivar had lower potentials (Figure 1) throughout the experimental period in the irrigated condition. After the irrigation was interrupted, an increasing reduction of ψ_am was observed. After three days of stress conditions, the decrease of ψ_am was more noticeable in the sensitive BR 16 plants (Figure 1). The BR 16 plants reached ψ_am values of −1.0 and −1.5 MPa, on the seventh and ninth day respectively, after the irrigation suspension. The Embrapa 48 plants reached the same ψ_am levels, on the ninth and eleventh days respectively, after the irrigation suspension. These results are consistent with a more efficient water use by the cultivar Embrapa 48 and confirm this cultivar as drought-tolerant. As the severe water potential of −1.5 MPa caused considerable damage to the BR 16 plants (Supplementary Figure S1), only the irrigated control plants and those that showed water potential of −1.0 MPa were used for the proteomic analyses.

FIGURE 1

Figure 1. Time course (days) of the water potentials (ψ_am in MPa) of soybean leaves after the suspension of the irrigation. Irrigated (IR) and non-irrigated (NI) treatments, respectively, for the tolerant (Embrapa 48) and the sensitive (BR 16) plants.

A significant portion of the cellular proteins are post-translationally modified by phosphorylation. However, phosphoproteins were observed at the low stoichiometric ratio and at low levels. 2DE gel staining tests for phosphoproteins showed a few spots stained with Pro Q Diamond (data not shown). Therefore, a strategy was applied for differential precipitation of high abundance proteins based on the addition of polyethylene glycol (PEG) to the aqueous protein extract, under non-denaturing conditions. The efficiency of the fractionation by PEG to remove higher abundance protein was evaluated in one-dimensional (1-DE) and 2DE gel electrophoresis (Supplementary Figure S2) by comparing the profile of the PEG-fractionated and non-fractionated samples. In all of the non-fractionated samples, the presence of the large and small rubisco subunits was observed, as indicated in Supplementary Figure S2. In contrast, the protein profile of PEG-fractionation on the supernatant samples did not show the protein bands corresponding to the rubisco subunits (Supplementary Figure S2). This method made it possible to see a larger number of stained spots using Pro Q Diamond.

The supernatant of the PEG-fractionated protein extracts, from the sensitive and the drought-tolerant soybean plants, was separated by 2DE. As the genotypes have different genetic backgrounds, contrasts were performed for each cultivar in irrigated (IR) and non-irrigated (NI) conditions. Afterwards, the drought-responsive protein profiles were analyzed to characterize the molecular and physiological responses. For the Embrapa 48, 643 well-separated spots were used to compare plants with −1.0 MPa water potential (NI) with the corresponding irrigated control (IR). 72 out of 643 were identified as differentially expressed spots, of which 60 were identified by mass spectrometry (Supplementary Table S2). These 60 spots corresponded to 45 different proteins as putative isoforms or modifications of the same protein. Among the 45 unique spots, 38 showed differential intensities and 7 were present only in plants under drought stress. In contrast, 4 spots were detected only in the control NI and 11 spots showed a reduction in relative abundance when comparing IR × NI. For the sensitive cultivar BR 16, 574 spots were used to compare IR vs NI plants. Only 58 were classified as differentially expressed spots, of which 42 were identified by mass spectrometry (Supplementary Table S3). These 42 spots corresponded to 35 different proteins, because of some possible isoforms or modifications. Among the 42 spots identified, 25 spots had an increase in relative abundance in NI treatment compared to the IR, 11 spots were present only in the IR treatment. The abundances (%V) of the responsive proteins to drought of each genotypes were converted to fold-change (Figures 2–4), and classified as up- or down-regulated, when an increase or reduction of the protein abundance was observed in the NI treatment, respectively. The protein abundance data were also separated by the response pattern to each genotype, when expressed in both (Figure 2) and when it was genotype-specific (Figures 3 and 4). For a better visualization of the general response differences between genotypes, the proteins were classified into functional categories and frequent biological processes, as the function of the interruption of water supply. The GO classification of the responsive proteins for both genotypes was in general similar, (Supplementary Figure S3) however, the specific expression profile was distinct (Supplementary Figure S4). The protein lists were also used for functional enrichment analyses using a ClueGO plugin to detect the significantly enriched GO terms when the proteome was up or down-regulated. A distinct response pattern was observed for the genotypes. Two main clusters were seen when comparing up-regulated proteins between genotypes (Supplementary Figure S5). The up-regulated proteins of the sensitive genotype BR16 share clusters related to biological processes, such as response to oxygen radical and regulation of translational termination (Supplementary Figure S5). However, for the tolerant Embrapa 48 (Supplementary Figure S5), the up-regulated proteins showed a distinct cluster related to positive regulation of catalytic activity, photosynthesis, and photosystem II assembly. The analysis applied for down-regulated proteins also showed distinct results for the genotypes, notably a cluster related to plastid translation and transcription, and chlorophyll biosynthesis for the drought sensitive BR 16 (Supplementary Figures S4 and S5). The individual analyses of the identified proteins showed a similar up-fold change for both genotypes such as the glutathione S-transferase TAU 19, glutathione peroxidase 6, transcriptional coactivator/pterin dehydratase and chloroplast drought-induced stress protein of 32 KD (Figure 2). This is an indicator of the water potential of the leaves and the general response of these soybean genotypes. However, some proteins classified as up-regulated were only identified in the tolerant Embrapa 48, such as rubisco activase and chloroplast RNA-binding protein 31B (Figure 3). Some proteins were only identified in the sensitive BR 16, such as the eukaryotic elongation factor 5A-1, transcriptional coactivator/pterin dehydratase and ascorbate peroxidase (Figure 4). These results are in accordance with the enriched GO terms (Supplementary Figure S5). Complex expression profiles were observed in the 2DE gels for some identified proteins because different spots were detected coding isoform for the same protein. For example, ten spots were detected for carbonic anhydrase (CA): spots 872, 881 and 69 as down-regulated in the Embrapa 48; 436 and 997 as down-regulated in the BR 16. One spot (93) was up-regulated in the Embrapa 48. Two spots were observed for gamma CA 1, as up-regulated in both genotypes and one spot for CA 2 as up-regulated in the BR 16.

FIGURE 2

Figure 2. Up- and down-regulated proteins of the soybean leaves from the drought sensitive BR 16 and the tolerant Embrapa 48 genotypes. Fold change was calculated between non-irrigated (NI) and irrigated (IR) plants from the same genotype. Different values of protein fold changes are shown with the same name because spots with different PIs were identified for the same Glyma ID by LC/MS. Values of fold change for protein spots when present only in the non-irrigated (NI) or irrigated (IR) treatments are shown as 10.00 or −10.00, respectively.

FIGURE 3

Figure 3. Up- and down-regulated proteins of the soybean leaves of drought sensitive BR 16. Fold change was calculated between non-irrigated (NI) and irrigated (IR) plants. Different values of protein fold changes are shown with the same name because spots with different PIs were identified for the same Glyma ID by LC/MS. Values of fold change for protein spots present only in the non-irrigated (NI) and irrigated (IR) treatments are shown as 10.00 or −10.00, respectively.

FIGURE 4

Figure 4. Up- and down-regulated proteins of the soybean leaves of drought sensitive Embrapa 48. Fold change was calculated between non-irrigated (NI) and irrigated (IR) plants. Different values of protein fold changes are shown with the same name because spots with different PIs were identified for the same Glyma ID by LC/MS. Values of fold change for protein spot present only in the non-irrigated (NI) and irrigated (IR) treatments are shown as 10.00 or −10.00, respectively.

The spot proteins stained with Pro Q Diamond were identified by mass spectrometry (Supplementary Tables S5 and S6) and were analyzed for differential expression, classified as up- and down-regulated for both genotypes (Figure 5) and genotype-specific (Figures 6 and 7). The phosphoprotein lists were also classified for GO functional categories and functional enrichment analyses by ClueGO plugin to detect the significantly enriched GO terms when the proteome was up-regulated (Supplementary Figure S6). A distinct profile was also observed for the drought contrasting soybean genotypes when the phosphoproteome was analyzed. Significant functional enrichment results were observed for the overexpressed phosphoprotein in the tolerant Embrapa 48 (in red), in which sharing clusters related to biological processes, such as oligosaccharides biosynthesis process, S-glycoside biosynthesis process and response to disaccharide stimulus (Supplementary Figure S6). On the other hand, the sensitive BR 16 enrichments were more evident for the down-regulated phosphoproteins mainly for biological processes related to response to disaccharide stimulus and glycolysis (Supplementary Figure S6). Stress-responsive proteins, which were detected as phosphorylated, showed a similar profile for both genotypes, such as the stress-responsive alpha-beta barrel domain protein and ascorbate peroxidase 1 (Figure 5). Six spots were identified as phosphorylated isoforms of the transketolase as up-regulated in the tolerant Embrapa 48. An enzyme involved in the biosynthesis of saccharides were also identified in the Embrapa 48, while four phosphorylated isoforms were down-regulated in the sensitive BR 16 (Figure 6). Further proteins were identified including; phosphoprotein and glycine decarboxylase P which showed expression up-regulated for Embrapa 48 and down-regulated for BR 16. This is a multi-protein complex, which plays a major role in photorespiration in plants, and is involved and related to the tolerance of the oxidative stress [26]. Spots were also identified as regulated phosphoproteins involved in the photosynthetic apparatus, such as photosystem II subunit O-2 and photosystem II subunit P-1 as up-regulated for tolerant Embrapa 48 and down-regulated for sensitive BR 16 (Figures 5 and 6). Other phosphoproteins related to carbon fixation processes such as carbonic anhydrase and rubisco activase were also observed as different protein spots. Spots for isoforms of the cobalamin-independent synthase family protein were found up-regulated in the Embrapa 48 (Figure 7). Protein homologous related to methionine biosynthesis is important in oxidative stress [27]. The protein IDs, for each phosphoprotein detected in the 2DE profiles, were confirmed for the presence of phosphoprotein in the public plant. Databases were used to verify those already described as phosphorylated proteins (Supplementary Table S4). In addition, the FASTA amino acid sequences were submitted to in silico analysis to predict the presence of conserved phosphorylated sites in these proteins (Supplementary Table S4). Most of the proteins staining in the gels as phosphoproteins, were validated as phosphoproteins present in the databases, or as showing a high score for the presence of the phosphorylated sites from the prediction results.

FIGURE 5

Figure 5. Up- and down-regulated phosphoproteins of the soybean leaves drought sensitive BR 16 and tolerant Embrapa 48 genotypes. Fold change was calculated between non-irrigated (NI) and irrigated (IR) plants from the same genotype. Different values of protein fold changes are shown with the same name because spots with different PIs were identified for the same Glyma ID by LC/MS. Values of fold change for phosphoprotein spot present only in the non-irrigated non-irrigated (NI) and irrigated (IR) treatments are shown as 10.00 or −10.00, respectively.

FIGURE 6

Figure 6. Up- and down-regulated phosphoproteins of the soybean leaves drought sensitive BR 16 genotype. Fold change was calculated between non-irrigated (NI) and irrigated (IR) plants. Different values of protein fold changes are shown with the same name because spots with different PIs were identified for the same Glyma ID by LC/MS. Values of fold change for phosphoprotein spot present only in the non-irrigated (NI) and irrigated (IR) treatments are shown as 10.00 or −10.00, respectively.

FIGURE 7

Figure 7. Up- and down-regulated phosphoproteins of the soybean leaves drought sensitive Embrapa 48 genotype. Fold change was calculated between non-irrigated (NI) and irrigated (IR) plants. Different values of protein fold changes are shown with the same name because spots with different PIs were identified for the same Glyma ID by LC/MS. Values of fold change for phosphoprotein spot present only in the non-irrigated (NI) and irrigated (IR) treatments are shown as 10.00 or −10.00, respectively.

Determining the impact of differentially phosphorylated proteins under water stress deficit in molecular and physiological processes may be challenging, as phosphoproteins are shown as many isoforms, and such spots could be interconverted by the addition of a phosphate group to the protein, therefore, changing the isoelectric points (PI). In this case, a spot will show an increase and the other one will show a reduction in volume. This could be more complex when multiples phosphorylated isoforms are present in the profile. We detected 13 spot clusters for BR 16 and 14 for Embrapa 48, containing different phosphorylated isoforms in the 2DE profiles (Supplementary Tables S6–S8). This phosphorylation behavior could be important in regulating the photosynthesis in the soybean genotypes.

The analysis of the metabolome was performed by GC/MS and the metabolic profiles were compared between the contrasting genotypes to verify the fluctuation of the metabolites during drought treatments related to the tolerance mechanism. This approach allows us to identify metabolites that were quantified and used for data analysis by MetaboAnalyst platform. The PLSA-DA method was used to analyse the most significant variables and observed differential grouping of the samples, in order to compare the genotypes and water regimes (Supplementary Figure S10). The BR 16 showed a better separation of the treatments, whereas, in Embrapa 48, the treatments remained more grouped. This indicates a lower fluctuation of the metabolites in this genotype, especially under moderate stress (Supplementary Figure S9).

The PLS-DA method was also used to identify the variables that had the biggest effect in grouping the samples (Supplementary Figure S11) and indicate the responsive amino acids in both genotypes. The metabolites that showed significant abundance variations as a function of the applied treatments are more responsive to drought (Figure 8). For both genotypes, the proline was the metabolite that showed higher variations, presenting a higher value VIP score when compared to the other fourteen most important metabolites. Proline is the most responsive amino acid to drought. The main responsive metabolites belong to the classes of organic acids, amino acids, and carbohydrates. When the metabolites of these classes were evaluated individually for the response to drought between genotypes, a pattern of behavior with a trend of greater abundance and response in the sensitive genotype BR 16 (Figure 8), such as the amino acids alanine, asparagine, aspartic acid, cysteine, isoleucine, lysine, methionine, tyrosine, phenylalanine, tryptophan and ornithine was observed. The amino acid proline showed a higher abundance in the sensitive BR 16, especially in severe drought stress. Likewise, some carbohydrates, such as cellobiose, xylulose, sucrose, xylose, melibiose, and trehalose, had the same effect of higher abundance in the sensitive BR 16 (Figure 8). Other metabolites that were responsive to drought treatment, such as putrescine, spermidine and ornithine, also presented the same behavior. The identified metabolites and the fold-change were used to map the pathway affected during stress by using the MetaboAnalyst platform. The number and the impact in the pathways were higher for the sensitive BR 16 (Tables 1 and 2). In these analyses, the pathway for amino acids metabolism showed higher impact indexes.

FIGURE 8

Figure 8. Relative abundances of some of the significantly responsive metabolites to non-irrigated treatments of the soybean leaves. Black bars indicate the metabolite intensities for the tolerant genotype Embrapa 48 and grey bars for the sensitive BR 16, irrigated and non-irrigated in the leaf potentials of −0.2 MPa (IR), −1 MPa (NI) and 1.5 MPa (NI), respectively.

TABLE 1

Table 1. Metabolic pathways significantly affected by drought stress for the tolerant Embrapa 48. Fold change was calculated between non-irrigated (NI) and irrigated (IR) plants from the same genotype and was used as an input for the pathway analysis.

TABLE 2

Table 2. Metabolic pathways significantly affected by drought stress for the tolerant BR 16. Fold change was calculated between non-irrigated (NI) and irrigated (IR) plants from the same genotype and was used as an input for the pathway analysis.

FIGURE 9

Figure 9. Absolute abundance of the phytohormone Abscisic Acid (ABA) of the soybean leaves quantified by LC MS. The treatments were irrigated (T1), under moderate drought stress (T2, −1.0 MPa) and under severe drought stress (T3, −1.5 MPa).

The absolute concentrations of nine phytohormones were determined by LC/MS and were used to identify which one had an effect on a water deficit. In Embrapa 48, only abscisic acid (ABA) showed significant variations in the relative concentrations between the analyzed treatments (Figure 9 and Supplementary Figure S12). In BR 16 genotype, ABA, carboxylic-1-aminocyclopropane (ACC, Ethylene precursor) and auxin (AIA) were statistically significant (Supplementary Figure S12). For both cultivars, the phytohormones response to drought treatment showed an increase in the relative concentration with the severity of the water deficit. As expected, ABA showed big variations in the absolute concentration in response to drought treatment, it was however, more noticeable for the sensitive genotype BR 16 (Figure 9).

DISCUSSION

Soybeans have been reported to have a wide variation in drought tolerance and these genetic resources have been used to develop new genotypes. Furthermore, genetic crosses were obtained to combine agronomic and tolerance characteristics [2,15,16]. The productivity analysis aims at distinguishing the different levels of tolerance among the cultivars, this can be carried out by many different complex physiological and molecular mechanisms. Trying to determine those contributing to the phenotype requires a combination of different experimental data. The cultivars BR 16 and Embrapa 48, share a common ancestor, and were studied by Oya [16] under drought conditions. They found that at the vegetative stage, the tolerant cultivar (Embrapa 48) showed a higher growth rate and a larger leaf area when compared to the sensitive BR 16 plants.

Initially the tolerance drought phenotype of the Embrapa 48 was confirmed by monitoring water potential under a reduction of the water supply (Figure 1), which was postponed by seven days to reach to −1.0 Mpa. Afterwards, the proteome and metabolomic profiles were compared. The proteome of the plant under drought stress was extensively studied and characterized, with a list of the responsive proteins [28]. However, in this study, we focused on the protein expression patterns that were distinct between the genotypes, aiming at identifying the physiological behavior and specific candidates that explain the drought tolerance phenotype observed. The global analysis of the protein abundance of both genotypes (Figure 2) indicated that some molecular and physiological processes for the Embrapa 48 were less affected by the reduction in the water supply. This can be explained by the presence of proteins related to functional categories that are indicative of the oxidative and water stress status, such as glutathione S-transferase TAU 19, transcriptional coactivator/pterin dehydratase and chloroplast drought-induced stress protein of 32 KD (CDSP32). CDSP32 is composed of two thioredoxin modules that have a critical role in plastid defense against oxidative damage. Furthermore, this is related to its role as a physiological electron donor to the peroxiredoxin [29]. An isoform of the Mog1/PsbP/DUF1795-like photosystem II reaction center PsbP protein was up-regulated only in the sensitive BR 16, which required photosystem II core assembly and stability [30].

The proteome response profiles were different between the genotypes, and indicated that the physiological and biochemical processes, such as photosynthesis, were less affected in the tolerant genotype. Under osmotic stress, the plant cells respond to limiting cellular damages and restores a new homeostasis, which includes a coordination of biochemical and physiological changes, including the photosynthesis inhibition [8]. In this context, it was observed that the photosynthetic electron transfer C protein showed a distinct expression between the tolerant and the sensitive genotypes (Figure 2), up-regulated in Embrapa 48 and down-regulated in BR 16. This encodes the rieske FeS center of the cytochrome b6f complex and the hcf mutant, in which the biosynthesis was inhibited, reducing the electron transport at saturating light intensities. Generally, drought reduces transpiration and at the same time photosynthesis. However, we observed that in the tolerant genotype the regulation of the protein expression could maintain higher photosynthetic activities in comparison to the sensitive one. This is in accordance with the higher growth rate that was observed for Embrapa 48 [16]. We identified nine spots for carbonic anhydrase 1 (CA1) and one for carbonic anhydrase 2 (CA2). However, six protein isoforms were down-regulated (Figure 2). CA1 and CA2 are involved in photosynthesis, and are usually depressed under drought stress in plants [10], and have a great effect on photosynthesis and water use efficiency [31]. Overexpression of CA1 or CA4 decreased stomatal density and improved water use efficiency [31]. CA is a metalloenzyme located in the chloroplast stroma, very close to rubisco, that catalyzes the reversible reaction of bicarbonate to carbon dioxide and maintains the supply of CO₂ for rubisco activity [32]. These mechanisms could be important in differentiating the responses between the genotypes and improving the water-use efficiency, and the leaf growth of the Embrapa 48, as reported by [16]. Other proteins related to photosynthesis were also stained as phosphorylated in 2DE profiles (Figure 7). A phosphorylated isoform for the rubisco activase and carbonic anhydrase was observed as up-regulated only in the tolerant genotype, which might reduce the damage on rubisco by drought stress (Figures 7 and 8). Hierarchical clustering by ClueGO revealed phosphoproteins related to some biological processes that were specific for the genotypes (Supplementary Figures S4 and S6). As reported by Rocha [33], phosphorylation could regulate the behavior of the soybean’s response to tolerance. The chloroplastic drought-induced stress protein was observed in different phosphorylated isoforms stained in the gel (32 KD). (Supplementary Table S8). This protein was also detected in spots showing reduced molecular weight, which is an indication that it is also regulated by proteolytic cleavage of the polypeptide chain.

A key role in the regulation of the response to drought is related to the hormone abscisic acid (ABA). The higher levels of this phytohormone were observed in the sensitive BR 16. Additionally, it was observed that this sensitive genotype had a greater number of metabolites showing significantly altered abundances in drought conditions and therefore, more noticeable perturbations in the metabolic pathways (Supplementary Figure S11). Proline is widely regarded as the main osmoprotectant in drought stress tolerance in plants. In our study on plants at the vegetative stage, we found that drought triggers an increase in the proline synthesis or accumulation on both genotypes (Supplementary Figure S11). However, higher levels were found for the sensitive BR 16. The role of proline in stress tolerance remains controversial as some authors have reported high proline levels in the susceptible cultivars subjected to stress conditions [34,35], while others have highlighted the opposite trend [5]. It has been suggested that proline functions are related to the plant water status but not as a measure of the level of tolerance [36], as also observed in soybean by Silvente [5]. In our study, the higher levels of ABA and proline identified in the sensitive BR 16 is in agreement with these findings. In addition, the metabolic profiles indicated that in the sensitive BR 16 genotype, the leaf cellular system is suffering greater damage when exposed to the same water potential of −1.0 Mpa. Thus, the signal for drought by ABA and proline were more noticeable in the sensitive BR 16. In the same way, amino acids and sugar were more abundant during drought in the sensitive genotype, which suggests that these compounds were not involved in osmoprotection of the tolerant genotype. This indicates that BR 16 has more evident cell damage and probably presents a higher protein degradation [37]. The proteome analysis is in accordance with the hierarchical clustering analysis by ClueGO, which found clusters of proteins related to the regulation of translational termination in the sensitive BR 16. The eukaryotic translation initiation factor 5A3 isoform X1 was found up-regulated in the BR 16 (Figure 2) and could be involved as an attempt to balance protein status. Wang [38] showed that this protein is a stress-responsive gene involved in the ABA signal transduction pathway, which facilitates protein synthesis and regulates several physiological pathways to improve stress tolerance. Another up-regulated protein only found in the intolerant Embrapa 48 was the chloroplast RNA-binding protein 31B (RBP). The importance of the mRNA in stabilizing mechanisms against oxidative or enzymatic degradation was observed in plants tolerant to salt and drought stress. Recently, it was shown that the RNA-binding protein (RBP) participates in the selective mRNA translation mechanism during hypoxia and drought by sequestrating mRNAs in stress granules. RBP was found up-regulated in both cells and plants upon long-term exposure to PEG and ABA [39]. Interestingly, when exposed to drought stress, RBP-overexpressing plants also accumulate less proline, suggesting that the higher tolerance observed is independent of the proline accumulation [39]. For sensitive BR 16, a differential expression of these RNA binding proteins was not found (Figure 2). However, the phosphorylated isoform for the chloroplast stem-loop binding protein of 41 kDa, that also binds and cleaves RNA, particularly in stem-loops and participates in chloroplast ribosomal RNA metabolism, was down-regulated for both genotypes (Figure 2). This protein was also reported as down-regulated in tomato plants [40] under stress drought.

A notable distinct protein expression profile between the genotypes was also observed for the protein Transketolase (Figures 7 and 8). We found nine phosphorylated isoforms in the tolerant Embrapa 48, and eight in the sensitive BR 16. Six were up-regulated (spots 322, 323, 320, 318, 921 and 316) and one was down-regulated (spot 319) for tolerant. However, for the sensitive genotype four were down-regulated (spots 1092, 1082, 607, 1021). Transketolases are key enzymes of the reductive and oxidative pentose phosphate pathways, and are responsible for the synthesis of sugar phosphate intermediates. They are expressed in three enzymatically active isoforms and their accumulation is associated with the rehydration process [41]. This observation suggests a possible role for these transketolase isoforms in the conversion of sugars during drought stress response. The complex phosphorylation behavior of this protein was confirmed by in silico prediction analysis of the phosphorylated sites (Supplementary Tables S6 and S7) and experimentally by Rocha [33].

CONCLUDING REMARKS

Drought tolerance in plants is due to many complex mechanisms and to evaluate which is the most significant for this phenotype is challenging because an extensive gene reprogramming is activated. Thus, a physiological and molecular analysis of the contrasting genotypes, showing a similar genetic background, could be an efficient approach. In this study, two Brazilian soybean genotypes were confirmed as drought-sensitive and tolerant by dissecting curves, protein expression and metabolites profiles, showing distinct responses. Our results indicate that the metabolic response observed is in accordance with a physiological study and explains the higher leaf growth rates observed for the tolerant genotype Embrapa 48. However, the tolerance appears not to because of osmoprotection by amino acids biosynthesis or by ABA-dependent cascades, as determined by metabolic profiles. In contrast, it is indicated that the tolerant genotypes maintain a photosynthetic steady-state even during drought conditions. A lower perturbation level of the metabolic pathways for the Embrapa 48 was observed. This behavior was also found in the gene expression profiles and suggests a lack of adaptive traits, which are important for the plants resistance to drought in the sensitive genotype [18]. Another explanation could be the occurrence of damage-related responses that may not have been induced in the tolerant cultivar under the same stress level. In physiological assays [17] lower rates of chlorophyll degradation were found for Embrapa 48, which could have allowed the recovery of the photosynthesis after re-irrigation and higher stomatal conductance values after re-irrigation [17]. In addition, the transpiration rate of both cultivars decreased with the development of water stress. However, transpiration was higher for BR 16 plants and reduced the soil moisture more rapidly than the Embrapa 48 [17]. Our proteomic and metabolomic data is in accordance with physiological and expression gene results and indicate that the metabolism of the tolerant Embrapa 48 is not impaired, or relative to the sensitive BR 16. Thus, regulatory cascades of the metabolic activities, mediated by protein phosphorylation, may be operating as a result of water saving mechanisms and of a more efficient translocation of the root-shoot system.

SUPPLEMENTARY MATERIALS

Supplementary File (PDF, 2.76 MB)

DATA AVAILABILITY

All data generated from the study are available in the manuscript or supplementary files.

AUTHOR CONTRIBUTIONS

HJOR, EGB, MEL, AM and EPBF designed the study; BPB, ROM, CEV and FSC performed greenhouse experiments; BPB, HJOR performed proteomic assays; LLL, FSC, CEV and HJOR performed metabolomic assays; HJOR, LLL, JCFS, FSC, FMSC analyzed the data; LLL, FSC, and HJOR wrote the paper with input from all authors; LLL, FSC, CEV and HJOR performed data integration.

CONFLICTS OF INTEREST

The authors declare that they have no conflicts of interest.

FUNDING

This research was funded by Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), Fundação de Amparo à Pesquisa de Minas Gerais (FAPEMIG), Conselho Nacional de Desenvolvimento Cientifico e Tecnológico (CNPq) and Instituto Nacional de Ciência e Tecnologia (INCT).

ACKNOWLEDGMENTS

The authors would like to thank the Núcleo de Análise de Biomoléculas (NuBioMol), the Universidade Federal de Viçosa, MG, Brazil, the Instituto de Biotecnologia Aplicada a Agropecuária (BIOAGRO), the Instituto Nacional de Ciência e Tecnologia em Interações Planta-Praga (INCTIPP) and the Brazilian Soybean Genome Consortium (GENOSOJA). This work was funded by Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), Fundação de Amparo à Pesquisa de Minas Gerais (FAPEMIG), Conselho Nacional de Desenvolvimento Cientifico e Tecnológico (CNPq) and Instituto Nacional de Ciência e Tecnologia (INCT).

How to Cite This Article

Lima LL, Balbi BP, Mesquita RO, da Silva JCF, Coutinho FS, Carmo FMS, et al. Proteomic and Metabolomic Analysis of a Drought Tolerant Soybean Cultivar from Brazilian Savanna. Crop Breed Genet Genom. 2019;1:e190022. https://doi.org/10.20900/cbgg20190022