Two cowpea cultivars with contrasting responses to drought but with similar flowering times and grain yield reductions were used. Fradel seeds (accession Cp5051) were obtained from INIAV (Instituto Nacional de Investigação Agrária e Veterinária); and Pinhel seeds (accession Vg50) were kindly supplied by the Departamento de Genética e Biotecnologia of UTAD (Universidade de Trás-os-Montes e Alto Douro). The drought stress response in the cultivars is manifest in two discrete phenotypes, namely increased root length in Pinhel after exposure to stress, and increased stem diameter in Fradel upon rewatering (Figure 12A).

The trials were conducted using 11-L plastic pots (25.5 cm in height) with drainage holes. Each pot was filled with 12 Kg of an 8:3 (w:w) mixture of soil and sand (Supplementary Table 1). Prior to sowing, the pots were thoroughly soaked with water to ensure 100% field capacity, which was the weight of the pots with saturated soil. To determine 0% field capacity, saturated soil was oven-dried to a constant weight. Plants were grown in a greenhouse at UTAD (N 41°17′7.28″, W 7°44′36.83″) with mean 43.55 ± 5.47/19.07 ± 2.88 °C day/night air temperatures, and a 16-h photoperiod, with photosynthetic photon flux density (PPFD) of up to 1701.42 ± 50.14 μmol m⁻² s⁻¹ (Supplementary Figure 1).

The layout of the experiment was a factorial design consisted of three adjacent blocks, one for each sampling event. Within each block, the pots were arranged in a randomized complete block design with four or six replicate pots per cultivar and treatment. The pots were sown on 18 May with four seeds at 80% field capacity, and thinned to two plants 15 days later.

Based on results from a pilot study to establish drought regimes, the drought experiment started 40 days after sowing, when plants began branching (Figure 1). The experiment focused on the progressive development of stress, as well as on two steady-state situations (e.g., Sicher et al., 2012), the first in which the plants were significantly stressed during the vegetative stage, but not lethally damaged (first block—first harvest), and the second in which stressed plants recovered during the reproductive stage upon rewatering (second and third blocks—second and third harvests). The first block included well-watered and drought-stressed treatments with six replicate pots. The well-watered control was maintained at 80% field capacity (by weighing the pots two times per day and adding the amount of water equal to the loss in weight), while in the other treatment plants were progressively stressed by stopping watering until the net photosynthetic rate (A) reached 0 on day 6 (D6), when plants were harvested (vegetative stage). The second block included well-watered, drought-stressed, and rewatered treatments, also with six replicate pots. When A reached 0 on D6, half of the plants were rewatered to return their soil to 80% field capacity and allowed to grow until A matched the control on day 12 (D12), when plants were harvested (reproductive stage). The third block mirrored the second one, but with four replicates and pots weighed every 2 days after day 12; plants were grown until grain maturity (D40), which provided a more detailed assessment of the treatment effects on yield and yield-related parameters.

At harvest, plants were uprooted before midday, washed gently with water to remove soil from roots, which were rinsed three times with 10 mM MES + KOH buffer (pH 6.5) to remove biological contaminants, and then blotted dry. Homogeneous healthy leaves of the same sizes and ages (three to five youngest fully expanded leaves) were removed from stems, and roots were cut off. Samples from both plants in each pot were pooled, flash frozen in liquid nitrogen, and stored at −80°C until analysis.

Leaf gas exchange measurements were performed on the uppermost fully expanded leaf of each plant in the second block every other day, using a portable infrared gas analyser (LCpro+, ADC, Hoddesdon, England) under the same conditions of growth. Net photosynthetic (CO₂ assimilation) rate (A), stomatal conductance (g_s), ratio between intercellular and atmospheric CO₂ concentration (Ci/Ca), and transpiration rate (E) were determined. Concurrently, chlorophyll a fluorescence was measured in situ using a pulse amplitude modulation system (FMS 2, Hansatech Instruments, Norfolk, England). Minimum fluorescence (F₀) was measured in 30 min dark-adapted leaves by applying a weak modulated light, and maximum fluorescence (F_m) was measured after applying a 0.7 s saturating pulse of white light (>1,500 μmol m⁻² s⁻¹). The maximum quantum yield of photosystem II was calculated as F_v/F_m = (F_m− F₀)/ F_m (Roháček, 2002). Following F_v/F_m estimation, after a 20-s exposure to actinic light (1,500 μmol m⁻² s⁻¹), light-adapted steady-state fluorescence yield (F_s) was averaged over 2.5 s, followed by exposure to saturating light (15,000 μmol m⁻² s⁻¹) for 0.7 s to establish F'_m. The sample was then shaded for 5 s with a far-red light source to determine F'₀. From these measurements the electron transport rate (ETR), photochemical quenching (qP), non-photochemical quenching (NPQ), effective quantum yield of PSII (Fv′/Fm′), quantum efficiency of PSII (Φ_PSII = ΔF/F'_m = (F'_m−F_s)/F'_m), and the fraction of PPFD neither utilized in photochemistry nor dissipated thermally (P_E) were calculated as described by Maxwell and Johnson (2000).

Primary metabolites were extracted and derivatized following the well-established method of Lisec et al. (2006). Frozen tissue was ground in liquid nitrogen to a fine powder with a mortar and pestle, transferred to a tube (100 mg) with 1.4 mL methanol, vortexed for 2 s, and 60 mL ribitol (0.2 mg/mL) added as internal standard; the samples and consumables were kept on ice to prevent thawing. The mixtures were shaken in a ThermoMixer at 950 rpm for 10 min at 70°C; tubes were opened occasionally to release pressure, and then centrifuged for 10 min at 11,000 × g at 4°C. The supernatant was collected, and the polar fraction was partitioned in 0.75 mL chloroform with 1,400 mL water. After centrifugation at 2,200 × g for 15 min, 150 mL of the upper layer was dried in a centrifugal concentrator overnight. In a derivatization reaction, the dry residue was dissolved in 40 μL of 20 mg/mL methoxyamine hydrochloride in pyridine at 37°C for 2 h and then ultrasonicated at 37°C for 30 min in 70 μL N-methyl-N-(trimethylsilyl)trifluoroacetamide that contained 20 μL/mL of a retention time index standard mixture composed of fatty acid methylesters in chloroform (0.4 mL/mL).

An aliquot (1 μL) of derivatized sample was injected into an Agilent 7890A gas chromatography (GC) (San Joseph, MI) in both splitless and split (1:30) modes at 230°C. Helium was used as the carrier gas with a flow rate of 2 mL/min, and separation was accomplished using a DB-35MS column (30 m × 0.32 mm, 0.25 μm; Agilent Technologies). The gradient program was 2 min at 85°C, and then 15°C/min to 360°C. Column eluates were introduced into a Pegasus HT time-of-flight mass spectrometer (TOF-MS, LECO, San Joseph, MI) through an electron ionization source operating in negative ion mode and ionized using a filament bias current of 70 eV. The transfer line and ion source were all held at 250°C. Mass spectra were collected from m/z 50 to 500.

Acquired chromatograms and mass spectra were evaluated using the built-in deconvolution algorithm of ChromaTOF®, and TagFinder (Luedemann et al., 2008). Compounds were identified by comparing mass spectral and retention time indices to those in the Golm metabolome database (www.mpimp-Golm.mpg.de/csbdb/gmd/gmd.html). Metabolite levels were determined as relative abundances normalized to the signal intensity of the internal standard ribitol and sample fresh weight (Supplementary Table 2).

Secondary metabolites, primarily phenolics, were extracted and characterized using liquid chromatography/photodiode array detector (LC-DAD), based on the method of Hachibamba et al. (2013). Naringin (20 μL) was added to 40 mg of sample as a peak reference, and phenolics were extracted with 1 mL of 70% methanol and continuously shaken for 30 min at ~25°C. Following centrifugation at 10,000 × g for 15 min, the resulting supernatants were passed through a 0.2 μm Spartan 13/0.2 RC Whatman filter.

Subsequently, 10 μL of solution was analyzed and quantified using a Gilson LC (Villers-le-bel, France) equipped with a Finnigan/Surveyor 81401 DAD (Thermo Electron, San Jose, CA). Samples were separated using a C18 column (150 × 4.6 mm, 2.6 μm; Phenomenex, Torrance, CA) at 25°C. The mobile phase consisted of 0.1% aqueous trifluoroacetic acid as solvent A and 0.1% trifluoroacetic acid in acetonitrile as solvent B, with a linear gradient elution at a flow rate of 1 mL/min as follows: 8–10% B (0–10 min), 10–20% B (10–35 min), 20–20% B (35–50 min), 20–45% B (50–60 min), and 45–8% B (60–80 min). Proanthocyanidins, phenolic acids, flavonoids, and anthocyanins were monitored at 280, 320, 360, and 520 nm, respectively.

Data acquisition and processing were performed with Xcalibur 2.0 and compounds were identified by comparing their retention times and spectroscopic data with those of authentic standards. In the absence of a standard, identification was made by matching DAD signals and comparing elution profiles to those reported in the literature. Quantification using peak areas was based on response factors (Supplementary Table 3), after normalization with the internal standard naringin and sample fresh weight.

The water content of shoots, leaves, and roots was estimated after oven drying. Immediately after completing the gas exchange analyses, the measured leaves were removed for leaf water potential (ψ_w) determination. A different set of leaves was used for the relative water content (RWC) and leaf solute (osmotic) potential (ψ_s), since the ψ_w measurements were followed by leaf destruction. RWC was calculated as the weight ratio between the leaf water content in situ and the completely turgid water content after rehydration (Kwasniewski et al., 2016). ψ_w was determined using a pressure bomb Model 1000 (PMS Instruments Co., Albany, OR) (Mendes et al., 2001; Warren et al., 2011). Immediately after the leaves were excised from the plant by using an oblique razor from 2 cm from the petiole, they were enveloped in the instrument chamber. The rate of increase of the inlet of gas into the cylinder was maintained constant at 0.02 MPa s⁻¹, until the xylem sap appeared at cut surface of the petiole; the readings were recorded in Bars and expressed in MPa. ψ_s was determined with a thermocouple psychrometer (SC-10, Decagon Devices Inc., Pullman, WA) in punched leaf discs (ca. 0.70 cm² each) frozen in liquid nitrogen, and kept at −80°C (Mendes et al., 2001). Frozen leaf discs were placed in psychrometer steel cups and left to equilibrate for 2 h before measurements and readings taken after 2 min. To convert the μV output to corresponding pressure values (MPa), a standard regression curve was obtained with eight known KCl concentrations in the expected ψ_s range. The pressure potential (ψ_p) was simply calculated by the equation: water potential—solute potential. In order to estimate osmotic adjustment, the obtained ψ_s values were corrected to full turgor with the RWC value of the corresponding sample (ψπ100 = ψ_s × RWC/100). The degree of adjustment was calculated as the absolute difference between ψπ100 values of the different treatments.

C and N were determined from 0.2 g of freeze-dried sample using the Dumas' combustion principle with Primacs SNC-100 (Skalar Analytical B.V, Breda, Netherlands). P and K were extracted from 0.3 g sample by digestion with H₂SO₄, followed by Skalar San⁺⁺ segmented flow analysis and Jenway PFP7 flame photometry (Essex, UK), respectively.

Shapiro-Wilk's and Levene's tests conducted using SPSS 16.0 (SPSS Inc., Chicago, IL) indicated that the assumption of normal distribution and variance homogeneity was violated in most cases. Therefore, data were arcsine or log transformed before analyses. Data are mainly presented as means ± standard deviation (SD, n = 4–6) or as response ratios calculated as log fold-changes in contents between drought-stressed and well-watered, rewatered and drought-stressed, or rewatered and well-watered plants. ANOVA, MANOVA, and MANCOVA were performed using SPSS, and differences between means determined by Tukey's and Bonferroni's tests. Probabilistic principal component analyses (PCA) and Student's t-tests were performed using the statistical tool of MS/data-independent-acquisition (Tsugawa et al., 2015). Correlations were established by calculating Pearson coefficients, using the cor.prob function of R (www.R-project.org/). Parameters for which comparisons met a threshold of significance were chosen for further analysis using general linear regression modeling in EXCEL. Bonferroni correction was applied for multiple t-tests or Pearson's correlation. Differences between treatments were highlighted in heatmaps with hierarchical clustering (HCA, R package “pheatmap”) and in pathway maps according to the KEGG database (www.genome.jp/kegg/pathway.html).

Under drought (Figure 1A), A decreased with stress duration and reached 0 at D6 (Figure 1B), at a 34% field capacity (Figure 1C). Photosynthesis clearly depended on soil water status, with a linear relationship between A and field capacity (Figure 1D). Similar relationships were found for transpiration rate (E) and stomatal conductance (g_s), but not for the other gas exchange and chlorophyll fluorescence parameters (not shown). A recovered to control levels (P < 0.05) after 4 days of rewatering (D10); however, plants were grown for another 2 days (D12) for equilibration. Pinhel had constitutively lower (P < 0.05) average g_s values than Fradel, 242.53 and 468.63 mmol m⁻²·s⁻¹, respectively (Figure 2A). Under water stress, Pinhel exhibited a slower rate and lower reduction (P < 0.05) of Fv′/Fm′ (effective quantum yield of PSII) and qP (photochemical quenching coefficient) (Figure 2B), which suggested that Pinhel maintained more efficient photochemical reactions. To further test that inference, hierarchical clustering was performed, and demonstrated that the primary differences in responses to stress and recovery were in Fv/Fm (maximum quantum yield of PSII) and P_E (fraction of PPFD neither utilized in photochemistry nor dissipated thermally), which were correlated with gas exchange parameters in Fradel but not in Pinhel (Supplementary Figure 2; Supplementary Table 4). More specifically, Fv/Fm decreased in Fradel from D8 to D12 but not in Pinhel; and P_E increased in Fradel but not in Pinhel by the end of the experiment (P < 0.05; Figure 2B), suggesting that under more severe stress conditions, the capacity of protective mechanisms in Fradel could be exceeded.

GC-MS-based metabolite profiling identified 41 primary metabolites, including sugars (5), polyols (4), amino acids (18), amino acid derivatives (6), and organic acids (8) (Supplementary Table 2). Within the LC-DAD profiling dataset, 35 peaks with spectral characteristic of phenolic acids (15), flavonoids (17), and proanthocyanidins (3) were annotated (Supplementary Table 3; Supplementary Figure 3). Peaks matching no known structure (12) were labeled “unidentified” or “unknown,” for those absorbing at 290 and 435 nm, respectively.

MANOVA data (not shown) indicated that cultivar had no significant effect on metabolite content (P = 0.367). Instead, most differences could be attributed to spatial distribution, with sugars (0.39 vs. 0.28 on average), amino acids (0.58 vs. 0.36), and proanthocyanidins (31.96 vs. 21.15) predominantly observed in roots vs. leaves (P < 0.05), most likely owing to their higher energy demand for the assimilation of soil resources. Changes induced by the treatments were also larger than differences between cultivars, indicating little metabolic divergence between Fradel and Pinhel.

Drought responses during the vegetative (first harvest at D6) and reproductive (second harvest at D12) stages were compared using MANCOVA, keeping sampling time as a within-subject factor. Time and all it interacting factors were significant predictors of metabolite changes, F_{(1, 2776)} = 21.723, P = 0.000. The time × treatment interaction indicated that the influence of stress on metabolites was more prolonged over time, with overall increases and decreases in the responses of 42.11 and 19.74% of metabolites, respectively, between D6 and D12 (data not shown). The MANCOVA model used here in no way assumed that metabolites had identical response directions, but simply indicated with good confidence a cumulative effect of stress on the metabolome. In Fradel leaves, for example, levels of phenylalanine increased 0.22-fold at D6 but 1.12-fold at D12 (Figure 3). Figure 4 showed that not only qualitative, but also quantitative changes occurred in the metabolome over time. After 6 days of stress, 26, 31, and 10 metabolites increased (P < 0.05) in Fradel leaves, Pinhel leaves, and Fradel roots, respectively; after 12 days of stress, 41, 35, and 35 metabolites were increased. Conversely, 25, 23, and 43 metabolites decreased (P < 0.05) in the same organs at first harvest and were reduced to 14, 13, and 26 at second harvest. Pinhel roots were the exception in that fewer metabolites increased with stress progression (34–18), whereas more metabolites decreased (5–15). Interestingly, in most cases, the rate at which a metabolite changed in response to stress closely matched the rate at which it reverted upon rewatering. For example, after 12 days, proline levels in Pinhel stressed roots decreased 0.20-fold, whereas those in rewatered roots increased 0.23-fold (Figure 5B).

First, t-tests were used to identify metabolites with the greatest responses to drought stress. Only metabolites with a log fold-change of at least 0.5 relative to the control at D6 or D12 were plotted: 18 for Fradel leaves, 22 for Pinhel leaves, 13 for Fradel roots, and 10 for Pinhel roots (Figure 3). In the leaves, the largest changes were in quercetin 3-O-6″-malonylglycoside (tentative identification; 3.09-fold at D12 in Pinhel), and quercetin (1.37-fold at D12 in Pinhel) for secondary metabolites, and phenylalanine (1.33-fold at D12 in Pinhel), and ornithine (1.18-fold at D12 in Fradel) for primary metabolites. The levels of these compounds generally increased, with the exception of glutamine, glycerate, γ-aminobutyrate, kaempferol derivatives, and unknown7, which decreased at both harvests. Pronounced increases in roots included kaempferol 3-O-diglycoside (2.55-fold at D6 in Pinhel) and quercetin 3-O-6″-malonylglycoside (1.67-fold at D12 in Pinhel).

Second, a snapshot of drought-induced changes in the metabolome across harvesting times, cultivars, and organs was obtained using MANOVA. In stressed plants, levels of 63.15% of the 88 assessed compounds were different (P < 0.05) from those in control plants (Supplementary Table 5). Metabolites whose fold-change was greater than 0.5 were listed as discriminatory in decreasing order by absolute value: quercetin 3-O-6″-malonylglycoside (2.14-fold), quercetin (1.12-fold), erythritol (0.80-fold), unknown7 (−0.75-fold), galactinol (0.65-fold), proline (0.63-fold), phenylalanine (0.62-fold), and glycerate (−0.53-fold). A 0.57-fold-change was observed in isoleucine, but with P = 0.082.

Third, an integrated PCA was used to generate an overview of the variance in the whole data matrix, including stress and recovery responses (Figure 6A). MANOVA results were consistent with PCAs. Organ-specific differences in metabolic products dominated the variances, with roots clustering apart from leaves without overlap on PC1, which represented 68.53% of the variation between groups. Metabolites differed more between treatments than between cultivars, although for secondary metabolites, Fradel samples clustered apart from Pinhel samples. PC2 (13.44% of total variance) was associated with metabolic responses to stress/rewatering, with isoleucine, proline, quercetin 3-O-6″-malonylglycoside, quercetin, galactinol, quercetin 3-O-diglucoside-4-O-glucoside, and phenylalanine the main contributors with positive coordinates, and unknown7, glycerate, glucose, and kaempferol 3-O-diglycoside the main contributors with negative coordinates (Figure 6B), as indicated by loading values (Supplementary Table 6). The loading plots of each class of compound separately unequivocally revealed galactinol, proline, glycerate (PC2−), and quercetin 3-O-6″-malonylglycoside, quercetin (PC2+), kaempferol 3-O-diglycoside, and unknown7 (PC2−) as the primary and secondary metabolites, respectively, contributing to the separation between watering treatments (Figure 6B).

Fourth, relationships between metabolites in roots and leaves were explored separately for each cultivar in additional PCAs (Figure 7). The PCAs indicated that nearly all control replicates were well-separated from drought-stressed replicates. In agreement with MANCOVA, progressive water deficit led to increased separation between treatments in the PCAs. With the exception of primary metabolites in Fradel roots, the drought treatment (PC1) was the main source of variance in the data. The top positive and negative loadings of PCs, distinguishing watering treatments for each cultivar and organ, corresponded to compounds that were previously selected using t-tests (Figure 3), with the addition of alanine, valine, catechin 7-O-glucoside and kaempferol 3-O-glucoside in Pinhel roots, and γ-aminobutyrate, galactinol, and unidentified4 in Fradel roots (Supplementary Table 6). The main compounds contributing to the dispersion of samples on the PCs and common to leaves and roots were proline, galactinol, quercetin 3-O-6″-malonylglycoside, and kaempferol 3-O-diglycoside. Rewatered replicates usually formed a distinct cluster near control replicates (except those corresponding to secondary metabolites in Fradel roots), indicating that plants reverted to their initial metabolic configuration (Figure 7). Indeed, Figure 4 showed that only 57% of Fradel root metabolites returned to normal status, whereas the percentage was over 75% for the other organs. Primary metabolites recovered more effectively, with only mannitol, glycerate (Fradel leaves), glucose, fumarate (Fradel roots), serine, isoleucine, ornithine, phenylalanine, and glycerate (Pinhel leaves) lagging behind (Figure 8). Besides quercetin and quercetin 3-O-6″-malonylglycoside, which were only partially restored after drought relief in all organs (suggesting they may be important during recovery), other compounds showed no biologically relevant pattern.

Fifth, treatments were grouped by similarities among variables using HCA (Figure 5A). The groups clustered according to watering regime, with a clear separation between stress and rewatering treatments, but also between the responses of leaves and roots to drought. However, in Fradel roots, recovery performance differed. For cultivars, the response was more variable and no clear separation was observed. Overall, quercetin 3-O-6″-malonylglycoside, kaempferol 3-O-diglycoside, quercetin, unknown7, galactinol, and proline (Figure 5B) responded the most to drought, as indicated by all statistical analyses. Five of these compounds grouped in the upper cluster of the hierarchy in Figure 5A, suggesting a coordinated response to drought.

Interestingly, compounds that increased or decreased in response to drought in the leaves were similar for Fradel and Pinhel. The following sustained trends were observed: (i) half of the organic acids increased while the others decreased; (ii) kaempferol and amino acid derivatives decreased; (iii) polyols, phenolic acids, unidentified phenolics, myricetin, quercetin, and catechin derivatives all increased; (iv) alanine, serine, asparagine, aspartate, glutamate, glutamine, glycine, and threonine decreased while the remaining amino acids increased; (v) raffinose and trehalose increased while sucrose, glucose, and fructose decreased; and (vi) half of the unknown compounds increased while the others decreased (Figures 9, 10).

In the roots, however, metabolites of Fradel and Pinhel generally shifted in opposite directions (Figure 11; Supplementary Figure 4). Pinhel roots initially accumulated both primary and secondary metabolites, but as drought stress continued, primary metabolites returned to control levels or decreased; although several secondary metabolites continued to increase, others were reduced. Meanwhile, in Fradel roots, primary and secondary metabolites initially decreased in response to drought stress; but as stress was prolonged, primary metabolites increased and several secondary metabolites recovered somewhat to levels that were still lower than the control. In general, metabolites belonging to same biochemical group changed in concert in roots but not in leaves. Another marked difference was with kaempferol 3-O-diglycoside, which decreased in the leaves in response to water deficit, but increased in the roots.

Water stress negatively affected all measured yield parameters, even after rewatering (Figure 12B). Correlations between metabolite levels and grain yield under all growth conditions showed that 18 primary and 32 secondary metabolites were correlated with grain yield with |r| > 0.800 in at least one of the four organ-by-cultivar comparisons (Figure 13A; Supplementary Table 7). Compounds (n = 18) present in at least two of the comparisons were further analyzed using HCA, which showed strong correlations (|r| > 0.990) between yield and quercetin 3-O-6″-malonylglycoside, proline, hydroxybenzoate, and threonate (Figure 13B).

The leaf and root water content of drought-stressed plants were similar or higher (P < 0.05) than those of control plants at every harvest. The shoot water content was also similar or higher at D6 and decreased by D12 (Figure 14A). Relative water content correlated positively with water potential; for the two parameters, differences between treatments were more evident at D12, with decreased levels in stressed plants. At D6, decline of osmotic potential at full turgor occurred in both cultivars, indicating the development of a foliar osmotic adjustment process in response to water stress. The degree of adjustment was 0.20 MPa in Fradel, and 0.14 MPa in Pinhel. However, no significant differences in adjustment were apparent between well-watered and drought-stressed leaves at D12 (Figure 14B).

In general, drought decreased the C/N ratio except in Pinhel leaves and Fradel roots at D6 (Figure 15A). Drought-stressed plants usually had higher K content (resulting in a low C/K ratio), whereas P content decreased and increased in leaves and roots, respectively (Figure 15B). A PCA that included data from the elemental analyses clearly identified K and C/K as important factors in drought response (Supplementary Figure 5). With few exceptions, values returned to control levels after rewatering.

Harvest dates were scheduled based on the loss and recovery of A, which occurred during the vegetative and reproductive stages, respectively. Cowpea responded to water deficit by showing a strong relationship between A and field capacity. Decreased photosynthesis is usually linked to limitations to CO₂ assimilation imposed by stomatal closure, but can also result from impaired mesophyll conductance and/or biochemical and photochemical constraints (Pinheiro and Chaves, 2011). Increased Ci/Ca (intercellular CO₂/ambient CO₂) and reduced A/g_s (intrinsic water use efficiency) were observed by D6 (Figure 2A), which indicated non-stomatal photosynthesis limitation. These changes coincided with an increase in NPQ (non-photochemical quenching coefficient; Figure 2B), suggesting prevention of photo-oxidative damage of the photosynthetic apparatus. Moreover, the fast recovery after irrigation suggests that other non-assimilatory processes had a role; the Mehler reaction and/or photorespiration, for example, may be enhanced to partially dissipate excess excitation energy or scavenge any reactive oxygen species produced (Pinheiro and Chaves, 2011; Sánchez-Martín et al., 2015). Therefore, photoinhibition was unlikely a major driver of drought stress in the cultivars. Pinhel, however, exhibited a finer modulation of the photosynthetic process than Fradel.

An important objective of this work was to identify metabolite markers for drought adaptation in cowpea. The approach used centered on the hypothesis that comparatively better survival and grain yield under drought and recovery must be reflected in developmental-, organ-, and cultivar-specific similarities in metabolic responses. Levels of 88 metabolites were measured during the vegetative and reproductive stages to assess metabolic acclimation to drought stress. Data showed that changes in the metabolome were more severe with extended periods of stress, irrespective of developmental stage (Figure 4). Therefore, flowering did not seem to alter metabolic responses to drought as seen with rice (Li et al., 2015; Raorane et al., 2015), and a progressive metabolic acclimation was indicated for the two cultivars. A unified response was evident from all statistical analyses (Student t-test, MANOVA, PCA), and quercetin 3-O-6″-malonylglycoside, kaempferol 3-O-diglycoside, quercetin, unknown7, galactinol, and proline were identified as having the most significant responses to drought stress (Figure 5B), assuming that fold-changes under stress, not basal levels, indicate resistance markers (Johnson et al., 2015). These conserved metabolic responses likely reflect the basic lynchpin metabolic acclimation of cowpea to drought stress.

The present study also examined whether roots or leaves would yield more metabolic information regarding drought adaptation. Metabolic responses to drought in leaves were similar during vegetative and reproductive growth (Figure 9). In roots, however, cultivar-specific metabolic responses were observed and were modulated by drought intensity, duration, and/or rate of progression (Figure 11). The changes in metabolic responses observed in roots could not be explained by differential growth inhibition, as suggested by Sanchez et al. (2012); genetic variation between cultivars seems more plausible. As cowpea responds to drought, metabolic pathways may be regulated in dissimilar orders before an optimal balance is established. However, harvest dates in this study were too distant to establish the timing of these changes. Moreover, changes in root metabolites could be part of a signaling network. The hormone abscisic acid acts as a dominant root-to-shoot signal under water deficit (Costa et al., 2011; Kwasniewski et al., 2016). Other signals have not been fully identified, and the interplay of hormones with sugars and redox signals is the focus of current research (Pinheiro and Chaves, 2011; Blum, 2017). The root metabolome as a whole, rather than a specific compound, could be acting as a signal transducer, thus modulating photosynthate investment in different parts of the plant. In summary, the data showed that both leaves and roots could be used to record changes in water conditions; however, leaves are probably more suitable for identifying water stress, while roots may exhibit distinct metabolic changes that indicate specific drought resistance mechanisms.

Potential stress markers are often identified as amenable to biotechnological exploitation based on the magnitude of their change. However, high levels of specific compounds that do not provide concordant yield increases are irrelevant (Hill et al., 2013). Therefore, the metabolic responses obtained in this study were placed in an agricultural perspective by combining them with yield data (Figure 13). Owing to their responsiveness, as highlighted by fold-changes and metabolite-yield correlations, three compounds were chosen as promising markers for yield performance under drought stress, namely proline, galactinol, and quercetin 3-O-6″-malonylglycoside. Although quercetin and unknown7 were significantly correlated with grain yield, they were not very responsive to drought stress in the roots. Curiously, kaempferol 3-O-diglycoside, which was consistent between cultivars and organs, failed to correlate with grain yield; it may respond regardless of the plant's ultimate ability to survive under stress.

To determine whether the observed metabolic responses indicated tolerance mechanisms, such as osmotic adjustment or cell wall extensibility, changes in osmosis with progressing drought stress were evaluated using osmotic potential, C/N ratio, and K data (Lugan et al., 2010; Gargallo-Garriga et al., 2015; Blum, 2017).

An “osmotic adjustment” indicates that solutes in plant tissues (osmotic potential) have accumulated at least as much as water potential or relative water content has been reduced (Blum, 2017); Fradel and Pinhel demonstrated such an adjustment, but only during the first days of stress (Figure 14B). The small and transient adjustment observed was unlikely the result of organic solute accumulation in the cells, given the cumulative effect of stress on the metabolome. Under water stress, N generally accumulated to comparable levels in both cultivars at the expense of C (Figure 15), in line with previous studies in grapevine (Hochberg et al., 2013) and rice (Raorane et al., 2015). The reduction in photosynthetic capacity during drought is thought to contribute to the shrinking of the total C pool (Jorge et al., 2015), which could explain the negative C balance observed. Enzymatic and molecular evidence argue for a negative effect of water deficit on N assimilation (Raorane et al., 2015). Thus, increased N can better be explained by increased catabolism of stored amino acids and protein pools. Cowpea evidently regulated its metabolism to balance the production of N- and C-containing metabolites, as indicated by increases and decreases in an equal number of primary metabolites in the leaves (Figure 9), which suggests that the total solute content remained constant.

Osmotic effects are controlled not only by organic compounds but also by inorganic solutes. K⁺ has a particularly well-known role in plant water balance (Gargallo-Garriga et al., 2015). As shown in Figure 15, K levels were generally much higher in drought-stressed plants, which could better account for the increase in cell osmolality of drought-stressed leaves. Therefore, K-based osmotic adjustment during the early stages of stress may constitute an important and energetically efficient adaptive feature in cowpea, but is insufficient in sustaining turgor under more severe stress conditions.

Overall, the data indicate that metabolic changes in cowpea are closely regulated and do not occur in isolation. Although these changes might not make a significant contribution to osmotic adjustment, they could still ensure appropriate homeostatic maintenance via non-osmotic roles, such as redox buffering, membrane stabilization, protein hydrotroping, radical scavenging, signaling, or N and C repository building (Sanchez et al., 2012; Johnson et al., 2015; Jorge et al., 2015; Blum, 2017). In the present study, this is indicated by the accumulation of several flavonoids with antioxidant activities (Nakabayashi et al., 2014; Simova-Stoilova et al., 2015) such as quercetin derivatives and proanthocyanidins with cell wall-stiffening functions (Hachibamba et al., 2013; Hochberg et al., 2013) such as catechin derivatives. The general emerging picture is that the cowpea metabolic configuration under water deficit is compatible with osmoprotective strategies, at least in the leaves, indicating that changes in proline, galactinol, and quercetin 3-O-6″-malonylglycoside levels are an integral part of an adaptive response rather than stress indicators. In the particular case of proline, the observed changes were unlikely a symptom of damage, given the high level of congruence between the responses to stress and rewatering and the differing responses in roots. Moreover, the accumulation of several other known compatible solutes (Hochberg et al., 2013; Meyer et al., 2014; Shelden et al., 2016) along with proline, such as phenylalanine, raffinose, and isoleucine, is unlikely to be coincidental. Thus, it is highly likely that the metabolite changes were functionally connected with an increase in resistance, and that their correlations with yield indicated beneficial effects.

Finally, it was important to determine whether these compounds were under common genetic control via a specific pathway, which could facilitate breeding for drought tolerance. Regulation of most pathways was consistent in the leaves of both cultivars (Figure 10) but differed in the roots (Supplementary Figure 4). The leaf metabolic map indicated a major shift in C partitioning away from glycolysis/gluconeogenesis into the raffinose/stachyose (via galactinol) and the shikimate (via phenylalanine) pathways. There is an increasing body of evidence for the implication of raffinose family oligosaccharides (RFOs) in drought tolerance (Lugan et al., 2010; Hochberg et al., 2013; Li et al., 2015). Myo-inositol, for example, has recently been identified as the most promising metabolite marker for increased maize yield under drought stress (Obata et al., 2015). However, RFOs may primarily serve as transient C storage sources, with advantages in mobility over sucrose and starch (Obata et al., 2015). The accumulation of branched-chain amino acids (valine/leucine/isoleucine pathway) may also be associated with storage of substrates for key metabolic pathways. For all combinations of cultivars and organs, differences in the amounts of phenylalanine were strongly related to its downstream metabolism, i.e., accumulation of phenolics. Interestingly, enhanced drought tolerance in plants via overaccumulation of phenolics and other secondary metabolites has been documented (Nakabayashi et al., 2014; Corso et al., 2015). In the roots (Supplementary Figure 4), opposing regulation of glycolysis/gluconeogenesis was found, with a significant increase in sucrose levels, which suggested a possible shift in the localization of energy-intensive processes, such as biomass production, from leaves to roots (Raorane et al., 2015; Simova-Stoilova et al., 2015). In contrast, the arginine/proline pathway seemed to receive additional N inputs (Figure 10); the pathway was far more stimulated in Fradel than in Pinhel by D12, which could be related to the greater damage sustained during water stress in the former. The map also highlighted two other pathways with important roles, with decreased levels of metabolites of the glycine/serine/threonine pathway, suggesting faster consumption of these metabolites through enhanced photorespiration (Sánchez-Martín et al., 2015) as indicated by chlorophyll fluorescence data, and the alanine/aspartate/glutamate pathway, suggesting increased energy consumption from reserves (Hill et al., 2013).

Overall, the mechanisms by which cowpea modifies its metabolism to meet the demands of diverse resistance functions when exposed to water deficit appear to be determined by the interplay between the shikimate and arginine/proline pathways, giving rise to three drought-responsive metabolites, namely proline, galactinol, and quercetin 3-O-6″-malonylglycoside. Mapping chromosomal regions jointly associated with these pathways and investigating their co-localization with quantitative trait loci in a larger population may indicate promising candidate genes for breeding.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.