Network Analyses and Data Integration of Proteomics and Metabolomics From Leaves of Two Contrasting Varieties of Sugarcane in Response to Drought.

Uncovering the molecular mechanisms involved in the responses of crops to drought is crucial to understand and enhance drought tolerance mechanisms. Sugarcane (Saccharum spp.) is an important commercial crop cultivated mainly in tropical and subtropical areas for sucrose and ethanol production. Usually, drought tolerance has been investigated by single omics analysis (e.g. global transcripts identification). Here we combine label-free quantitative proteomics and metabolomics data (GC-TOF-MS), using a network-based approach, to understand how two contrasting commercial varieties of sugarcane, CTC15 (tolerant) and SP90-3414 (susceptible), adjust their leaf metabolism in response to drought. To this aim, we propose the utilization of regularized canonical correlation analysis (rCCA), which is a modification of classical CCA, and explores the linear relationships between two datasets of quantitative variables from the same experimental units, with a threshold set to 0.99. Light curves revealed that after 4 days of drought, the susceptible variety had its photosynthetic capacity already significantly reduced, while the tolerant variety did not show major reduction. Upon 12 days of drought, photosynthesis in the susceptible plants was completely reduced, while the tolerant variety was at a third of its rate under control conditions. Network analysis of proteins and metabolites revealed that different biological process had a stronger impact in each variety (e.g. translation in CTC15, generation of precursor metabolites, response to stress and energy in SP90-3414). Our results provide a reference data set and demonstrate that rCCA can be a powerful tool to infer experimentally metabolite-protein or protein-metabolite associations to understand plant biology.

Introduction

Understanding how plants adapt to limiting stress conditions, such as drought, is key to improve global food security. Too much or too little water induces different responses in plant gene expression and metabolism, which we still poorly understand. In Brazil, around 9 million ha are used for sugarcane plantations, producing mainly sugar and ethanol. In recent years, considerable drought events during vegetative growth has led to major losses in productivity, prompting the breeding programs to search for drought tolerant genotypes. Such plant materials are particularly useful for the identification of genetic and metabolic mechanisms involved in sugarcane tolerance to drought, by exploring the omics technologies. Single omics analysis often does not provide enough information to understand the biology of an organism. While the integration of multiple omics may give a better understanding of the biological system, as a whole (Montastier et al., 2015). Global network-based approaches authorize multiple datasets analyses and carry the advantage of highlighting functionally related pathways and biological entities of potential relevance, as hubs (Naylor et al., 2008). rCCA can be applied as an exploratory approach that aims exploiting the linear relationships between two datasets of quantitative variables from the same experimental units (Gonzalez et al., 2008). These assumptions rely on the fact that similar expression patterns across a set of samples are hypothesized to have a functional relationship (Lee et al., 2004). In the present study, we explored the changes in proteome and metabolome responses of two contrasting sugarcane varieties, submitted to a progressive drought. CTC15 (tolerant) and SP90-3414 (susceptible) were submitted to 12 days of drought regime, under greenhouse conditions. The aim of this work was to identify how the leaf proteome and metabolome in these two contrasting varieties interact, leading to different responses to drought. In parallel, we propose the use of rCCA to explore the correlation between metabolites and proteins identified by Label-free quantitative proteomics and GC-MS-based metabolomics. We chose sugarcane (Saccharum spp.) as our model plant because this perennial grass, cultivated in tropical and sub-tropical countries, is the world’s most important sugar crop with a high capacity to store sucrose in its culms, as a primary energy source (D’Hont et al., 2008). It has become an important bioenergy source and is classified among the most important crop for first and second-generation ethanol production (Heringer et al., 2015; Hoang et al., 2015). During sugarcane vegetative growth, drought may lead to significant decreases in productivity and yield (Zingaretti et al., 2012). Besides, drought impacts also in membrane integrity, pigment content, osmotic adjustment, water relations, photosynthetic activity (Benjamin and Nielsen, 2006), limiting the areas suitable to sugarcane cultivation. Understanding the molecular mechanisms involved in drought tolerance is a key step to develop sugarcane varieties capable to use water efficiently. Such approach becomes highly significant when one considers the fact that sugarcane has one of the most complex plant genomes, with a high chromosome number and high degree of aneuploidy (D’Hont et al., 1996). This chromosomal complexity makes it difficult to obtain improvements of sugarcane through conventional breeding. In the last few years some publications reported advances in the study of sugarcane responses to drought, using different strategies such as differential gene expression (Vargas et al., 2014; Vantini et al., 2015; Li et al., 2016), sRNA regulation (Thiebaut et al., 2014), morphological and physiological analysis (Cia et al., 2012; Zhao et al., 2013; Jain et al., 2015) and proteomics by two-dimensional gel electrophoresis (Jangpromma et al., 2010; Zhou et al., 2012; Almeida et al., 2013; Khueychai et al., 2015). To our knowledge this is the first exploratory analysis applying rCCA to associate metabolite and protein data in sugarcane response to drought. These omics approaches provide the opportunity to evaluate cellular behaviors from a multi-level perspective and enhance our understanding of sugarcane biology.

Methods

Plant Material and Experimental Conditions

To study changes in protein and metabolite profiles in response to drought, two sugarcane varieties were used: Cpan class="Chemical">TC15, previously reported as tolerant to drought (Graça et al., 2010; Kido et al., 2012; Silva et al., 2013; Thiebaut et al., 2014), and SP90-3414 susceptible to drought (Cia et al., 2012; Kido et al., 2012; Silva et al., 2013). Stalks from both varieties were kindly provided by CTC (Centro de Tecnologia Canavieira, Piracicaba, Brazil). The experiment was conducted in a greenhouse with a temperature range of 28 °C ± 2 during the day and 22 °C ± 2 during the night. Single stalks were planted, germinated and grown in 50 L pots with a mixture (2:1) of substrate (BasaPlant®, Brazil) and vermiculit. Plants were arranged in a completely randomized experimental design, under daily irrigation and the soil moisture was maintained near to 100% of field capacity (FC, equivalent to -234 KPa). After 5 months of growth, leaf gas exchange measurements and samples were obtained from plants at three levels of water availability in the soil: permanently irrigated (PI), 4 and 12 days without irrigation (4 DI and 12 DI, respectively). The soil moisture was maintained in the range of 100%, 40% (-286 KPa), and 25% (-850 KPa) of FC, respectively. Gas exchange rates were determined from five plants of each variety, chosen randomly in the greenhouse. After gas exchange measurements, the top visible dewlap leaves (+1) were detached from nine sugarcane plants for each treatment, frozen immediately in liquid nitrogen and stored at -80°C. For metabolomics each leaf separately was used for metabolite extraction while for proteomics, each three leaves were used to form one biological replicate. In total we used nine and three biological replicates for metabolomics and proteomics, respectively. In order to build the network, we needed to have the same number of biological samples for metabolomics and proteomics. The drought treatment was applied at the vegetative stage of sugarcane as it has been identified to be critical of water demand (Ramesh, 2000); mainly because 70-80% of sugarcane yield is achieved during this developmental stage (Singh and Rao, 1987).

Gas Exchange Measurements

Photosynthetic rate (A), transpiration rate (E) and stomatal conductance (g s) were measured using the portable photosynthesis system Li-6400 (LI-COR Biosciences, Inc.). Light curves were obtained with photosynthetic active radiation (pan class="Chemical">PAR) values of 2000, 1500, 1000, 500, and 0 µmol m-2 s-1. The measurements were recorded between 8:00 to 11:00 am in the medium portion of leaf +1 completely expanded. Leaf temperature was kept at 25 °C during the IRGA analysis.

Protein Extraction and Trypsin Digestion

Total protein was extracted from each biological sample, using the phenolic method according to Hurkman and Tanaka (1986), with minor modifications. Leaf tissues were ground into a fine powder in liquid nitrogen, and 200 mg was homogenized in 1 ml of extraction buffer (0.7 M sucrose, 0.5 M Tris-HCl, pH 7.5, 50 mM EDTA, 0.1 M KCl, 1% w/v polyvinylpolypirrolidone (PVPP), 2% v/v 2-mercaptoethanol, and 2 mM PMSF), by shaking for 30 min at 4°C. After, 800 µl of Tris-HCl saturated phenol pH 8.5 was added to the protein suspension and samples were shaken for 30 min at 4°C, the phases were separated by centrifugation (10,000 g for 30 min at 4°C). The supernatant was recovered and re-extracted with an equal volume of extraction buffer. This step was repeated two times. Proteins were precipitated by adding 5 vol of 0.1 M ammonium acetate in methanol and incubated overnight at 22°C. The samples were then centrifuged (9,000 g, 30 min at 4°C) and the resulting pellets were washed three times with 0.1 M ammonium acetate in methanol, followed by a wash with acetone. After complete evaporation of acetone at room temperature, the protein pellet was resuspended in 400 µl of solubilization buffer [7 M urea, 2 M thiourea, 0.4% v/v Triton X-100 and 10 mM dithiothreitol (DTT)]. Complete protein solubilization was achieved by vigorous shaking using a vortex for 2 min. Protein extracts were desalinized using Amicon Ultra-0.5 ml 3K-NMWL filter devices (Millipore Corporation). Total proteins were quantified using the Bradford method (Bradford, 1976). Fifty micrograms of proteins were denatured with 25 µl of 0.2% RapiGest (Waters, USA), reduced with 2.5 µl of 100 mM DTT and alkylated with 2.5 µl of 300 mM iodoacetamide. Trypsin digestion was performed with sequencing Grade Modified Trypsin (Promega) at a 1:100 (w/w) enzyme: protein ratio and proteins were incubated at 37°C overnight. After, 10 µl of 5% (v/v) trifluoroacetic acid (TFA) was added to the digested mixture to hydrolyze the RapiGest (Waters, USA). The peptide mixture was then desalted using ZipTip C18-columns (Millipore Corporation). The final volume of 40 µl was obtained by the addition of 20 mM ammonium formate (pH 10) solution containing 200 fmol/µl of rabbit phosphorylase (internal standard to data normalization and label free protein quantification- P00489) to the lyophilized, desalted peptide sample.

LC-MSE Analysis

The peptides mixture was analyzed by reverse-phase ultra-performance liquid chromatography (ACQUITY UPLC M-Class System with 2D Technology—Waters, USA) using a Synapt G2 HDMS (Waters, Manchester, UK). First dimension chromatographic separation was achieved using an AQUITY UPLC M-Class pan class="Chemical">peptide BECH C18 columns (5 µm, 300 µm x 50 mm). Elution was performed using five different binary gradients with 20 mM pH 10 ammonium formate in acetonitrile, at a flow rate of 2 µl min-1. Eluted peptides from the first-dimension column were trapped in a Symmetry 2D C18 column (5 µm, 180 µm x 20 mm) and diluted, online, with acetonitrile containing 0.1% formic acid. Second dimension separation was performed in an AQUITY UPLC M-Class CSH C18 column (1.7 µm, 75 µm x 150 mm), using a binary gradient from 7% to 85% of acetonitrile with 0.1% formic acid, during 75 min, at a flow rate of 400 µl min-1.

Mass spectrometry acquisition was achieved in a Synapt G2 HDMS (multiplexed DIA—data-independent acquisition) (Waters, Manchester, UK) mass spectrometer equipped with ion mobility cell and a nanolockspray source in the positive ion and “V” mode.

The doubly-charged ion [(M + 2H)2+] was used for initial single-point and MS/MS fragment ions of GFP (Glu 1)-Fibrinopeptide B m/z 785,84,206 [(M + 2H)2+] (Waters, pan class="Chemical">Corp., Milford, USA) were used as lock masses and instrument calibration, respectively. Data-independent scanning (MSE) experiments were performed by switching between low (3 eV) and elevated collision energies HDMSE (19–45 eV), applied to the trap “T-wave” CID (collision-induced dissociation) cell filled with argon gas. Scan time of 0.8 s were used for low and high energy scans from m/z 50 to 2000 (Silva et al., 2006).

Processing Parameters and Database Search

The raw data processing, protein identification and relative quantitative analyses were all performed using ProteinLynx Global Server v2.5.1 (PLGS, Waters). Proteins identifications were performed with PLGS searching into the SUCEST translated EST database (http://sucest-fun.org/), containing specific annotation of sugarcane ESTs (237.954 sequences) (Vettore et al., 2003). The processing parameters included: automatic tolerance for precursor and product ions, minimum of three fragment ions mapan class="Chemical">tched per peptide, minimum of seven fragment ions matched per protein, minimum of two peptides matched per protein, one possible trypsin missed cleavage, carbamidomethylation of cysteine as fixed modification and oxidation of methionine as variable modification, and a maximum false positive discovery rate (FDR) lower than 1%, determined based on the search of a reversed database, which was generated automatically using PLGS by reversing the sequence of each entry. To identify and quantify the proteins, the intensities of the spectra were calculated by the stoichiometric method, according to the internal standard, the sequence of rabbit phosphorylase (Uniprot entry: P00489), by MSE analysis (Silva et al., 2006) and normalized using the PLGS auto normalization function. The amount and sequence of the matched peptides and protein fmols were obtained, based on the ratio of its three most abundant peptides (“High Top 3” method) determined in each individual experiment (Silva et al., 2006), considering the data from three biological replicates for each treatment. All protein hits were identified with a confidence of > 95%. Only those proteins identified in at least two out of three biological replicates were considered for analysis. Proteins identified but not quantified were removed from the analysis. The mass spectrometry proteomics data was deposited to the ProteomeXchange Consortium via the PRIDE (Vizcaíno et al., 2016) partner repository with the dataset identifier PXD014920 (username: reviewer11150@ebi.ac.uk; password: bs65Tv52)

Functional Annotation and Classification

Gene Ontology (GO) annotation, to assign biological process terms, was performed based on Blast2GO software using the blastp algorithm whith the nr database (11/2018) and E-value ≤ 1e-05, the other parameters were setted as default (Conesa et al., 2005).

Metabolites Extraction and Derivatization

Prior to metabolite extraction frozen leaf tissue (50 mg) was grounded under liquid nitrogen using a vibration mill (MM 301 Repan class="Chemical">tch GmbH & Co, Haan, Germany) set to a frequency of 30 Hz s-1 for 45 s, with pre-chilled holders and 3 mm tungsten beads. Leaf metabolites were extracted by adding 500 µl of pre-chilled methanol (MeOH): chloroform (CHCl3): water (H20) (6:2:2 v/v/v). The extract was vortexed vigorously, sonicated at 40 Hz s-1 for 15 min and centrifuged at 14.000 x g (Eppendorf Centrifuge 5415R) for 10 min at 4°C. The supernatant was filtered (Millipore filter PVDF 0.22 µm) and 100 µl of each sample was transferred to vials and evaporated until dryness.

Samples were derivatized according to Gullberg et al. (2004) with 30 µl of methoxyamine hydrochloride (15 mg ml-1) in pyridine for 16 h at room temperature. The samples were trimethylsilylated by adding 30 µl of N-methyl-N-(trimethylsilyl) trifluoroacetamide (MSTFA) containing 1% trimethylchlorosilane (TMCS), the resulting mixture stand at room temperature for 1 h. After silylation, 30 µl of heptane was added. Stable isotope reference compounds [1 mg ml-1 each of (13C3)-myristic acid, (13C4)-palmitic acid and (2H4)-succinic acid] were added in samples prior to derivatization and used as external standard for quality control. Derivatized samples were analyzed according to Gullberg et al. (2004). Blank control samples and a series of n-alkanes (C12–C40), which allowed retention indices to be calculated (Schauer et al., 2005) were also used.

Analysis and Metabolites Identification by GC-TOF-MS

One microliter of each derivatized sample was injected splitless into a gas chromatograph 7890A (Agilent Technologies, Santa Clara, USA) coupled with a Comb-xt Autosampler (Leap Technologies, Carrboro, USA). A 20m x 0.18 mm i.d. x 0.18 µm film thickness, DB5 GC capillary column (Agilent Technologies) was used as the primary column and the secondary GC columns was a 0.69 m x 0.1 mm i.d. x 0.1 µm film thickness (Rxi-17 Restek, Bellefonte, USA). The injector temperature was 280°C, the septum purge flow rate was 20 ml min-1 and the purge was turned on after 60 s. The gas helium flow rate through the column was 1 ml min-1, the column temperature was held at 80°C for 2 min, then increased by 15°C min-1 to 305°C, and held there for 10 min. The column effluent was introduced into the ion source of a GC×GC/TOF-MS (Pegasus 4D, Leco Corp., St. Joseph, USA). The transfer line and the ion source temperatures were 280 and 250°C, respectively. Ions were generated by a 70-eV electron beam at an ionization current of 2.0 mA, and 10 spectra s-1 were recorded in the mass range m/z 45–800.

The ChromaTOF software v. 4.51 (Leco Corp., St. pan class="Disease">Joseph, USA) was used to perform baseline correction and export all MS files in NetCDF format. Peak detection, retention time alignment and library matching were performed using the TargetSearch package (Cuadros-Inostroza et al., 2009). Metabolites were identified by comparing their retention indexes (± - 2 s) and spectra (similarity > 600) against the compounds stored in the Golm-Metabolome-Database (http://csbdb.mpimp-golm.mpg.de/csbdb/gmd/gmd.html) (Kopka et al., 2005). Metabolites intensities were normalized by dry weight and total ion chromatogram (TIC).

Statistical Analysis

Multivariate (Principal Component Analysis—PCA) univariate (ANOVA, FDR adjusted p ≤ 0.05) analysis and Heatmaps were done in MetaboAnalyst 4.0 (Xia et al., 2015). To reduce systematic variance and to improve the performance for downstream statistical analysis data were log-transformed and pan class="Chemical">Pareto scaled prior to data analysis.

Network Analysis

The network analysis was used to assess the relationship between metabolome and proteome data from PI versus 4DI treatment and PI versus 12DI, for both varieties. It was performed using the rCCA, which is a modification of the classical canonical correlation analysis (CCA), a multivariate statistical method used to assess pan class="Chemical">correlations between two multivariate datasets (Gonzalez et al., 2008). This was performed using the mixOMICS v5.2 (Rohart et al., 2017) in R (http://www.R-project.org). Regularization parameters (λ1 and λ2) were estimated using the tune.rcc function, where a set of positive values chosen to evaluate the cross-validation score (cv-score) for each point in the network, achieving an optimal value for λ1 and λ2 that offered the highest cv-score. After, the rcc function was used to define the canonical correlations and the canonical variates between the two datasets. The unit representation plots were created using the plotIndiv function, in this plot each biological sample is projected into the XY-variate space. The correlation circle plots were created using the plotVar function, where all variables are represented through their projections onto the planes, highlighting the direction and correlation values of each variable. The network function was used to produce relevant networks from the similarity matrix derived from rCCA. The threshold was set to 0.99 and only the associations exceeding this threshold were shown. This threshold values were stringent to obtain biologically interpretable networks (Gonzalez et al., 2012). In order to analyze the created networks they were exported to Cytoscape 3.3 network visualization software (Shannon et al., 2003). The used layout was the organic, with the diameter of the node relative to the number of edges and number of connected edges. The correlation values is highlighted by the intensity of edge color ranging from dark green (negative correlation) to dark red (positive correlation).

Results

Changes in Photosynthetic Capacity Under Drought

Photosynthetic rate (A), transpiration (E) and stomatal conductance (g s) were determined in leaves (+1) to verify the physiological effects of drought upon each variety, under 4DI and 12DI, corresponding to 40% and 25% FC, respectively. As expected, it was possible to observe differences in performance between tolerant and susceptible varieties, under drought. Cpan class="Chemical">TC15 showed no major changes in gas exchange at the initial drought exposure, with no significant reduction in photosynthesis, transpiration and stomatal conductance between PI and 4DI ( ). On the other hand, photosynthesis diminished at 12DI by around one third of the values observed in PI. These results indicated the ability of the tolerant variety to maintain key physiological process, during a short-term drought. In SP90-3414, photosynthesis reduced significantly at 4DI and mainly in 12DI plants ( ). Leaf transpiration and stomatal conductance also diminished after 4DI, compared to the PI ( ). In comparison to the tolerant variety, SP90-3414 was already suffering a limitation in the photosynthetic capacity at 4DI and was mainly limited at 12DI, probably due to a prolonged stomatal closure, leading to a reduced carbon fixation. During the drought treatment, SP90-3414 variety showed signs of stress, mainly at 12DI such as wilting leaves and progressive leaf rolling.

Figure 1

Effects of drought on photosynthetic rate—A (A), stomatal conductance—gs (B), and Transpiration rate—E (C) in the variety CTC15 (drought-tolerant). Comparison among permanently irrigated (PI), 4 days without irrigation (4DI) and 12 days without irrigation (12DI). Letters indicate differences significant by Tukey’s test at 5% probability level.

Figure 2

Effects of drought on photosynthetic rate—A (A), stomatal conductance—gs (B), and Transpiration rate—E (C) in the variety SP90-3414 (drought-susceptible) cultivar. Comparison between permanently irrigated (PI) and 4 days without irrigation (4DI) plants. Letters indicate differences significant by Tukey’s test at 5% probability level.

Effects of drought on photosynthetic rate—A (A), stomatal conductance—gs (B), and Transpiration rate—E (C) in the variety Cpan class="Chemical">TC15 (drought-tolerant). Comparison among permanently irrigated (PI), 4 days without irrigation (4DI) and 12 days without irrigation (12DI). Letters indicate differences significant by Tukey’s test at 5% probability level.

Proteome Changes in Sugarcane Leaves in Response to Drought

By using spectral counting-based label free quantification and searching the SUCEST project database (http://sucest-fun.org/), in total, 1650 and 1750 proteins were identified and quantified at 95% confidence level, in CTC15 and SP90-3414 varieties, respectively. Venn diagram was performed to identify the common and unique proteins in each variety ( ), providing an overview of whether the changes represented common or specific responses. In Cpan class="Chemical">TC15 and SP90-3414 plants were found 1218 and 1336 common proteins, respectively. In both varieties, 12DI plants showed the highest number of unique proteins (200—CTC15; 182—SP90-3414) followed by 4DI (176—CTC15; 142—SP90-3414) and PI (56—PI; 90—SP90-3414) plants. The unique proteins ( ) were removed from the matrix submitted to statistical analysis. Proteomics changes in response to drought were first examined to identify the protein distribution within the treatments of each variety. To reduce the dimensionality of the data and visualize samples grouping tendency and outliers, a PCA was performed ( ). All data were lying inside the 95% confidence region (Hotelling’s T2 ellipse). The score plot ( ) showed a clear trend of treatments separation, indicating a metabolic adjustment to drought at protein level, in both varieties. The first principal component (PC1) separates the PI plants, of each variety, from the 4DI and 12DI plants explaining 36.8% ( ) and 31.5% ( ) of the total variation, respectively. 4DI and 12DI plants were separated by PC2, representing 15.7% ( ) and 18.8% ( ) of the total variance in CTC15 and SP90-3414 varieties, respectively. To identify if both varieties would present similar responses to water deficit, two Venn diagrams were constructed, based on the differentially abundant (p ≤ 0.05) and unique proteins ( ). Considering the differentially abundant ones, only two proteins (1.9%) [fumarylacetoacetate hydrolase domain-containing protein 1 (SCJFRZ2015E05) and prli-interacting factor l-like (SCRULB1060G07)] were common between CTC15 and SP90-3414 varieties. However, 100 unique proteins (13.4%) were common to both varieties ( ), most of them belongs to the category Metabolic process—GO:0008152 (25%) and Cellular process—GO:0009987 (17%).

Figure 3

Principal Component Analysis (PCA) score plots of proteins identified in the three treatments (permanently irrigated (PI), 4 days without irrigation (4DI) and 12 days without irrigation (12DI) in sugarcane leaves. (A) Drought-tolerant variety CTC15. (B) Drought-susceptible variety SP90-3414.

Figure 4

Heatmap based on differentially abundant proteins (p ≤ 0.05, FDR adjusted) from the drought-tolerant variety CTC15. PI (permanently irrigated), 4DI 4 days without irrigation and 12DI 12 days without irrigation.

Heatmap based on differentially abundant proteins (p ≤ 0.05, FDR adjusted) from the drought-tolerant variety Cpan class="Chemical">TC15. PI (permanently irrigated), 4DI 4 days without irrigation and 12DI 12 days without irrigation.

Protein Changes in CTC15 and SP90-3414 in Response to Drought

By performing Anova (FDR adjusted p ≤ 0.05) 49 differentially abundant proteins were identified in the tolerant variety ( ), their distribution can be observed in the heatmap ( ). Aiming to identify which biological processes were most strongly related to drought stress, these proteins were classified, using Blast2GO. The five most represented process were: Metabolic process—GO:0008152 (30%), Translation—GO:0006412 (26%), Cellular process—GO:0009987 (20%), Catabolic process—GO:0009056 (20%), and Biosynthetic process—GO:0009058 (12%). From the 49 differentially abundant proteins ( ), 35 (71%) increased in abundance upon drought, while 6 and 7 proteins had its content reduced upon drought at 4DI and 12DI, respectively. For instance, at 4DI a zeaxantin epoxidase (SCBFRZ2017C12), a polypan class="Chemical">galacturonanse precursor (SCJFRZ3C01H09), ndph-dependent reductase (SCRUFL1114B10.b), and two ribosomal proteins, had their abundance significantly reduced. At 12DI we observed the reduction in content of the following proteins: glycine cleavage h-protein (SCEZHR1084D12), plasma membrane atpase (SCCCRT1002F07), ribosomal protein (SCAGLR1043F11), elongation factor (SCCCCL3003D04.b), and a fumarylacetoacetate hydrolase domain (SCVPCL6043F04). Most of these highly abundant proteins at 4DI and 12 DI were related to Translation process (GO:0006412). Some of the proteins found are known to be drought responsive (i.e. two glutathione-s-transferase (SCBGLR1003E05 and SCAGFL1089C03), ferredoxin-sulfite reductase (SCQGLR1019B02), lipoxygenase (SCJLRZ1024A04), nadh-ubiquinone oxidoreductase (SCAGLR1021H11), vacuolar atpase b subunit (SCQGRT1039F06), endoplasmin (SCSGFL1078A06.b), 60s ribossomal protein (SCRFLR1012C07, SCBFLR1026E01, SCCCCL3140B12).

In SP90-3414 plants, 59 proteins were differentially abundant (FDR adjusted p ≤ 0.05) ( ). Their distribution pattern can be observed by the heatmap in . The five most representative biological process found for these proteins were: Metabolic process—GO:0008152 (29%), Cellular process—GO:0009987 (22%), Generation of precursor metabolites and energy—GO:0006091 (10%), Carbohydrate metabolic process—GO:0005975 (10%), and Response to stress—GO:0006950 (8%) ( ). From the 59 proteins, 15 (25.4%) were differentially abundant in all treatments and most of them were related to the Metabolic process (GO:0008152) category. From these proteins, 60%, 27%, and 13% were highly abundant at 12DI, 4DI, and PI, respectively. Some of the proteins found were pan class="Chemical">abscisic acid stress ripening protein (SCCCST3002H03), glutathione s-transferase (SCJFRT1008A09), glutaredoxin (SCCCLR1C03E08), 60s ribosomal protein (SCBGLR1120F10), light-induced protein (SCCCRT2002C02), and mRNA binding precursor (SCRFSB1021A11). As observed for CTC15, most of the proteins (52.5%) accumulated at 4DI and 12DI, in comparison to PI. Glutathione reductase (SCRLAD1040A04), ferredoxin-chloroplastic-like (SCJLFL4100F04), glutathione transferase (SCCCCL4003D01), heat shock (SCCCCL3005G07.b), and phospholipase (SCAGLR1043E06) are some of the proteins found. Eight proteins (13.5%) were highly abundant at PI and 4DI, in comparison to 12DI and five proteins (8.4%) were highly abundant at PI and 12DI ( ). In this variety, photosynthesis related proteins such as thylakoid membrane phosphoprotein (SCMCSD1061A04), chlorophyll a-b binding protein (SCEZLR1031G02), photosystem i reaction center subunit (SCQGLR2025B12) were found.

Figure 5

Heatmap based on differentially abundant proteins (p ≤ 0.05, FDR adjusted) from the drought-susceptible variety SP90-3414. PI (permanently irrigated), 4DI 4 days without irrigation and 12DI 12 days without irrigation.

GC-TOF-MS Based Metabolic Changes in Sugarcane Leaves in Response Drought

Since metabolic regulation is a major mechanism to overcome drought in plants, metabolomics analysis was performed aiming to obtain an overview about the metabolic changes occurring in tolerant and susceptible sugarcane varieties, exposed to drought. After manual inspection, we could annotate a total of 104 and 97 metabolites, for CTC15 and SP90-3414, respectively ( and ). Few metabolites from CTC15 (probable a Glycosyl and Amino acid derivative) and SP90-3414 (probable Amino acid derivative, probable Fatty Acyls, probable Carboxylic acid, Probable Steroid/steroid derivative and probable Oxazanine) were annotated based on the best mass spectral match from the database, however, their identity refered to metabolites that were, probably, not previously reported in plants. They were kept in our analysis and annotated based on their chemical category. Organic acids, sugars, lipids and others were the main classes found in both varieties ( and ). Venn diagram was performed to identify the common and unique metabolites in each variety ( ), providing an overview of whether the changes represented common or specific responses. In CTC15 and SP90-3414 plants were found 97 and 87 common metabolites, in at least two treatments, respectively ( ). The number of unique metabolites found were 2 (PI), 1 (4DI), and 4 (12DI), in CTC15 and 5 (PI), 3 (4DI), and 2 (12DI) in SP90-3414 ( ). To obtain the metabolite distribution pattern within the treatments, of each variety, the unique metabolites (found in a single treatment) were removed from the matrix submitted to statistical analysis. To investigate how each variety responded to drought, a PCA was performed for each ( ). The PCA score plot revealed differences in the metabolome of each treatment, suggesting a metabolic reconfiguration in response to drought. All data were lying inside the 95% confidence region (Hotelling’s T2 ellipse). The sum of the first two principal components explained 49.8% and 44% of the variation found in CTC15 and SP90-3414 varieties, respectively. In CTC15, despite the dispersion inside PI and 4DI samples, components 1 and 2 separated PI from 4DI and 12DI ( ). The same dispersion was observed for SP90-3414-PI samples. For this variety, component 2 has a tendency to separate PI from samples submitted to drought 12DI ( ). To identify the metabolites significantly affected by drought, one-way Anova (FDR adjusted p ≤ 0.05) was performed for each variety and revealed 10 and 28 metabolites differentially abundant in CTC15 and SP90-3414, respectively ( and ). Sugars, organic acids and lipids were the main classes of metabolites identified. Prior to Anova, samples, from both varieties, showing high variation were removed. The distribution of differentially abundant metabolites, in each variety, was represented by heatmaps ( ). Three differentially abundant metabolites [2-(4-hydroxyphenyl)-Ethanol], hexadecane and taurine were found in both varieties. 2-(4-hydroxyphenyl)-Ethanol is a phenolic compound, reduced in CTC15 at 12DI, in comparison to PI and 4DI. In SP90-3414, it accumulated at PI. Hexadecane, an alkane, accumulated at PI in both varieties. Taurine, and amino acid reduced under drought. In both varieties, most of the differentially metabolites accumulated in PI and reduced in response to drought. In the CTC15 variety, the metabolites O-methyl-L-Threonine, 2-piperidinecarboxylic acid, and pyroglutamic acid were indicated as differentially abudant because they were not found at 12DI, however, any of them changed between 4DI and 12DI and will not be discussed.

Figure 6

Principal Component Analysis (PCA) score plots of metabolites identified in the three treatments permanently irrigated (PI), 4 days without irrigation (4DI) and 12 days without irrigation (12DI) in the leaves. (A) Drought-tolerant variety CTC15. (B) Drought-susceptible variety SP90-3414.

Figure 7

Heatmap based on differentially abundant metabolites (p ≤ 0.05, FDR adjusted). (A) drought-tolerant variety CTC15 and (B) drought-susceptible variety SP90-3414. PI (permanently irrigated), 4DI 4 days without irrigation and 12DI 12 days without irrigation.

Principal Component Analysis (PCA) span class="Chemical">core plots of metabolites identified in the three treatments permanently irrigated (PI), 4 days without irrigation (4DI) and 12 days without irrigation (12DI) in the leaves. (A) Drought-tolerant variety Cpan class="Chemical">TC15. (B) Drought-susceptible variety SP90-3414.

Heatmap based on differentially abundant metabolites (p ≤ 0.05, FDR adjusted). (A) drought-tolerant variety Cpan class="Chemical">TC15 and (B) drought-susceptible variety SP90-3414. PI (permanently irrigated), 4DI 4 days without irrigation and 12DI 12 days without irrigation.

Acessing Metabolite-Protein Networks in Sugarcane Leaves Under Drought by rCCA

It is expected that under stress conditions plants adjust their metabolic responses, in order to maintain its physiological homeostasis. Aiming to access part of these responses, rCCA was performed and metabolite–protein networks were identified, at a highly stringent threshold of 0.99 ( and ). Due to the high stringency of the analysis.

Figure 8

Correlation network depicting correlations derived from rCCA between metabolites and proteins from drought-tolerant variety CTC15. CTC-PI-12DI networks. PI (permanently irrigated) and 12DI (12 days without irrigation 12DI).

Figure 9

Correlation network depicting correlations derived from rCCA between metabolites and proteins from drought-susceptible variety SP90-3414. (A) SP-PI-4DI networks. (B) SP-PI-12DI networks. PI (permanently irrigated, 4DI (4 days without irrigation) and 12DI (12 days without irrigation 12DI).

Correlation network depicting pan class="Chemical">correlations derived from rCCA between metabolites and proteins from drought-tolerant variety CTC15. CTC-PI-12DI networks. PI (permanently irrigated) and 12DI (12 days without irrigation 12DI).

Correlation network depicting pan class="Chemical">correlations derived from rCCA between metabolites and proteins from drought-susceptible variety SP90-3414. (A) SP-PI-4DI networks. (B) SP-PI-12DI networks. PI (permanently irrigated, 4DI (4 days without irrigation) and 12DI (12 days without irrigation 12DI).

Networks were performed for each variety separately, combining PI-4DI and PI-12DI data, aiming to access the metabolites and proteins correlations at each level of stress. Given the linear pan class="Chemical">correlation of variables used in rCCA, a correlation in the same direction is called positive and in the opposite direction is called a negative correlation. The number of positive and negative correlations and the complete description of each network are reported in the and . The unit representation and the correlation circle plot of the rCCA can be visualized in the and , respectively. The unit representation plot of the tolerant variety in PI-4DI analysis ( ) did not show a clear separation between PI and 4DI. This might be due the fact that this tolerant variety, after 4 days of drought, shows no major differences in relation to control plants. In this sense, the network correlation is not represented. show a clear separation between treatments of tolerant and susceptible varieties, CTC PI-12DI, SP PI-12DI, and SP PI-4DI, respectively. As shown in , the plot did not revealed differences between CTC-15-PI and CTC-15-4DI. For this reason, the network correlation analysis was not done for CTC15-PI-4DI. All correlation circle plots ( ) showed a large number of metabolites and proteins with a high correlation value. Due to the large number of variables in our datasets, a clearer view of these correlations can be highlighted in the correlation networks. In the tolerant variety, three clusters were found in PI-12DI (lambda 1 = 1, lambda 2 = 0.99, cv-score = 0.97). In terms of biological process Metabolic process (GO:0008152) and Response to stress (GO:0006950) represented most of the proteins found in the CTC15-PI-12DI network ( ). CTC15-PI-12DI cluster 1 ( ) is the biggest one, refers to Hexadecanoic acid (ID 1391, positively correlated to 15 proteins). Cluster 2 refers to three metabolites: 3-Methyl-2-oxopentanoic-acid (114), O-Benzoyl-beta-D-glucoside (1847), and similar to Inositol (1350), positively correlated to five proteins. Cluster 3, the smaller one, contains three metabolites: Galacturonic acid (1268), Galactaric acid (1377), and 17-beta-Estradiol (1955), correlated to three proteins. In this cluster, the only positive correlation was between Galacturonic acid and protoporphyrinogen chloroplastic (SCMCRZ3064A10). The enzymatic role of the proteins identified in the networks were accessed and grouped according to metabolic pathways ( ), based on EC number.

In SP90-3414 variety, four and five clusters were found for PI-4DI (lambda 1 = 1e-07, lambda 2 = 0.049, cv-score = 0.93) and PI-12DI (lambda 1 = 0.018, lambda 2 = 1e-04, cv-span class="Chemical">core = 0.99), respectively. Based on the biological process Cellular (GO:0009987) and Metabolic process (GO:0008152) were the most represented in SP90-3414 PI-4DI and PI-12DI networks ( ; ). PI-4DI cluster 1 is the biggest, containing four metabolites (Adenine—ID 1179, Indole-3-acetaldehyde enol—ID 1311, Octadecanal—ID 1359, Octadecanoic acid—ID 1626) and 12 proteins. Octadecanoic acid was negatively correlated with six proteins that were positively correlated with the other three metabolites. Cluster 2 consist of allo-Threonine (ID 297), positively correlated to two proteins. Cluster 3 is based on Phosphoric acid monomethyl ester (ID 108), positively correlated to one protein and, cluster 4 is based on trans-3-hydroxy-Proline (417), positively correlated with two proteins. In SP90-3414 PI-12 five clusters were found. Cluster 1, the biggest, is composed of six metabolites (Norvaline—ID 153, allo-Threonine—ID 297, trans-3-hydroxy-Proline—id 417, Digitoxose—ID 616, 3-hydroxy-Indole—ID 793, Lithocholic acid—ID 2412) correlated (positive and negative) with 27 proteins. Cluster 2 is based on three metabolites (N-(3-aminopropyl) Morpholine – ID 688, alpha-1-O-methyl-Galactopyranoside—ID 981, Normetadrenaline—1469) correlated (positive and negative) with 13 proteins. Cluster 3 is based on three metabolites (Adenine—1179, Indole-3-acetaldehyde enol—ID 1311, Octadecanal—ID 1359) correlated with 10 proteins. One negative correlation was found in this cluster (adenine correlation with one protein). The metabolites in cluster 3 are the same found in SP90-3414 PI-4DI cluster 1. Cluster 4 is based on two metabolites [2-(methylaminomethyl)-Tartronic acid—ID 453, N-acetyl-Mannosamine—ID 1431] correlated (positive and negative) with two proteins. Cluster 5 is based on Threonine (308) positively correlated with a single protein. Regarding the enzymatic role of the proteins found in SP90-3414 networks, a completely different set of elements were observed. The common elements are mostly related to photosynthesis and oxidation, reflecting the major impacts caused by drought. In PI-4DI, 3 elements were unique in this condition and 15 were common to 12DI. This suggests that photosynthesis and oxidation have a higher impact at 4DI. However, at 12DI most of the elements were unique to this condition. These elements belong to important primary and secondary metabolic pathways, indicating that at 12DI, the susceptible variety could be requiring a different metabolic adjustment, to face drought. The comparison between varieties, at the same condition, revealed clear differences in response to drought ( ). At 4DI, while CTC15 required amino acid and photosynthetic ajdustments while SP90-3414 needed to adjust metabolism mainly in reponse to oxidative stress. Twelve elements were common to both varieties, other twelve were found only in CTC15 and a lower number, six were found in SP90-3414. This could also indicate that both varieties start to perceived drought at 4DI. The difference is that, CTC15 sucssesfully adjust its metabolism in order to keep stomatic conductance and photosynthesis, being able to maintain homestasis. On the other hand, to maintain cellular homeostasis, SP90-3414 needs to adjust its metabolic and physiological responses by closing stomatal conductance and reducing photosynthesis, already at 4DI. The changes between both varieties, at 12DI, showed that metabolic activity is still ongoing but, different pathways were required. CTC15 is dealing with amino acid and lipids degradation and oxidation process, while SP90-3414 faces alterations in carbon fixation and eletron transport chain. These results indicated that each variety demands a specific adjustment to metabolism maintanence and survival, under drought.

Discussion

Network biology is a rapidly developing area of research, which states that biological processes are not chiefly controlled by individual molecules or by discrete, unconnected linear pathways (Charitou et al., 2016). Here, we propose the application of rCCA to explore correlation structures between metabolites and proteins, as an exploratory strategy to improve our knowledge about leaf metabolism in pan class="Species">sugarcane under drought, which is one of the main factors negatively impacting sugarcane productivity worldwide.

Physiological Responses of CTC15 and SP90-3414 to Drought

Plants have developed physiological mechanisms to save water during drought periods. As an example, stomatal closure is the first line of defense against drought (Silva et al., 2013). Therefore, photosynthesis is particularly impacted to drought once stomata closure reduces pan class="Chemical">CO2 diffusion to the fixation sites, in leaf mesophylls (Ghannoum, 2009). In accordance to this, all parameters measured (A, E and g ) were reduced at 4DI and 12DI in SP90-3414 ( ). For the susceptible variety, stomata closure occurs earlier than observed in the tolerant, in order to reduce water loss and maintain plant metabolism. The reduction in photosynthetic capacity was already reported for several sugarcane varieties susceptible to drought (Rodrigues et al., 2009; Rodrigues et al., 2011; Kido et al., 2012; Silva et al., 2013), supporting our observations with SP90-3414. For the tolerant variety, the observed differences were not statistically significant between PI and 4DI ( ), only at 12DI. Similar observations were reported for the tolerant sugarcane variety SP83-5073. While analyzing A, E, and g the authors found that tolerant plants behaved exactly as the control plants under 4 days of drought.

Proteome Dynamics of Sugarcane Leaves Under Drought

It is known that in C4 plants photosynthetic inhibition under drought is mainly due to stomatal closure (Ghannoum, 2009), resulting in lower carbon fixation, followed by the over-reduction of the electron transport system and production of pan class="Chemical">ROS (Nouri et al., 2015). Although physiological data showed that tolerant variety, at 4DI ( ), did not differ significantly from control plants, the proteome exhibited differences, indicating that protein changes may precede physiological responses to drought, in this variety. Besides, no protein directly involved in carbon fixation was differentially abundant in CTC15 ( ), while, in SP90-3414 four proteins [chlorophyll a-b binding protein (SCEZLR1031G02), photosystem I reaction center subunit (SCQGLR2025B12), ferredoxin- chloroplastic-like (SCJLFL4100F04), ribulose-phosphate 3-epimerase (SCAGHR1017F03)] accumulated under drought ( ). The proteomic adjustment of each variety to drought seems to require different strategies. Considering the GO terms, most of the differentially abundant proteins from CTC15 and SP90-3414 belong to Metabolic process (GO:0008152). Considering that this term is a general one, attributed to a variety of proteins, this observation was not a surprise. Instead, the most represented terms, such as Translation (GO:0006412) in CTC15, Generation of precursor metabolites and energy (GO:000609) and Response to stress (GO:0006950), found in SP90-3414, contributed to reveal variety differences in response to drought.

In general, elevated levels of proteins involved in translation can be beneficial for maintenance of protein synthesis under drought (Merewitz et al., 2011). Accordingly, different proteoforms of them were found in our work. 40s (SCCCCL3080D07), 60s (SCRFLR1012C07, pan class="CellLine">SCBFLR1026E01, SCCCCL3140B12), ribosomal s8 (SCEQRT1031E01), and 26s proteasome regulatory subunit (SCVPLR1049F01, SCCCCL4005D12) proteins accumulated under drought. On the other hand, two ribosomal proteins [30s (SCCCSD1095A07), 60s (SCVPLR2019F01)] accumulated at 12DI and PI, and two proteins (eukaryotic translation initiation factor 3, ribosomal protein rps20) accumulated at 4DI and PI. Elongation factor is a core translational protein that catalyzes the initiation and elongation of newly growing peptide chains (Wan and Liu, 2008). Two proteoforms of elongation factors were found, SCVPLB1020E07 accumulating under drought and SCCCCL3003D04.b accumulating at PI and 4DI. This differential change may indicate that the form SCVPLB1020E07 is directly involved maintaining the synthesis of proteins directly involved with drought. It is important to highlight that four proteins, differentially abundant, related to translation [60s (SCBGLR1120F10) and 40s (SCCCLR1001H01, SCCCCL3080E05, SCCCLR2C02A08)] were found in SP90-3414. As mentioned above, the terms Response to stress (GO:0006950) and Generation of precursor metabolites and energy (GO:000609) contribute to characterize the protein adjustment of SP90-3414 variety, to drought. Among the Response to stress (GO:0006950), one dehydrin (dhn2-like protein—SCCCLR1001F04) was found, accumulating during drought. Dehydrins are a group of late embryogenesis abundant (LEA) proteins, which protects cellular membranes and organelles during dehydration (Wahid and Close, 2007) and it is known to be induced by drought (Wahid and Close, 2007; Hand et al., 2011; Ferreira et al., 2017). This protein was differentially abundant only in the susceptible variety, suggesting that membrane protection is earlier required in SP90-3414 than in tolerant variety. Heat shock protein that belongs to response to stress is another common protein responsive to drought (Wang et al., 2003). In our work, one heat shock protein was found in each variety. In CTC15, the abundance of SCCCRZ1003A03 increased according to drought duration. In SP90-3414, SCCCCL3005G07.b accumulated in response to drought. The Generation of precursor metabolites and energy (GO:000609) category include proteins related to the primary metabolism, such as carbon fixation, glycolysis, citric acid cycle (TCA), and electron transport chain. These proteins were found in SP90-3414, accumulating in response to drought. Since primary metabolism can be generally defined as essential for survival (Aharoni and Galili, 2011), our results suggest that the susceptible variety priorize the metabolism of energy generation as a strategy to face drought. Two differentially abundant proteins, prli-interacting factor (SCRULB1060G0) and fumarylacetoacetate hydrolase domain (SCJFRZ2015E05), were differentially abundant in both varieties. In CTC15 it accumulated in response to drought while, in SP90-3414, was gradually accumulating in response to drought. This protein is related to ubiquitin-dependent protein catabolic process and was found in sugarcane callus (Heringer et al., 2015), however, its function is still poorly understood. Fumarylacetoacetate hydrolase domain (FH) exhibited the same pattern in CTC15 and SP90-3414, accumulating in response to drought. FH is a main enzyme acting in the final step of tyrosine (Tyr) degradation pathway (Shao et al., 2015). FH accumulated in leaves of maize seedlings (Shao et al., 2015) and it was reported to be involved with cell death in Arabidopsis (Han et al., 2013). It is possible that in sugarcane leaves, FH has a role acting in cell death, under drought.

Metabolomic Responses of Sugarcane Leaves Under Drought

The metabolomic data also suggests a high plasticity of the sugarcane leaf metabolome, as a strategy to adjust its metabolism under drought. In terms of differentially abundant metabolites, a higher number was found in SP90-3414, this result may suggest that drought impact a different set of metabolites, in each variety and, the susceptible genotype shows the participation of a higher number of metabolites to face drought. This could also be a metabolic evidence of specific variety response, to drought. In acpan class="Chemical">cordance with this, only three differentially abundant metabolites [2-(4-hydroxyphenyl)-Ethanol, Hexadecane and taurine] were found in both varieties. The role of 2-(4-hydroxyphenyl)-Ethanol in plants still unknown. It was isolated from Papaver rhoeas leaves (Hillenbrand et al., 2004). Since this metabolite accumulated at 12DI, in both varieties, this phenolic compound could be involved in degradation or being the product of degradation, due to drought. Hexadecane was reported in pepper species (Kevrešan et al., 2013; Bogusz Junior et al., 2015). In our work, this alkane accumulated at PI and decrease under drought, indicating that either this metabolite or its precursors are degraded in response to drought. Taurine is an amino acid, known to act as osmolyte, accumulating in response to drought, salt and extreme temperatures (Liang et al., 2013). However, in our work, taurine reduced under drought. About the specific responses of each variety to drought, the sugars threose and rhamnose were not found at 12DI, both accumulated at PI and reduced the abundance at 4DI, in the CTC15 variety.Threose is one of the products from the oxidative degradation ascorbate, an important oxidant (Shumilina et al., 2019). Contrary to our findings, this sugar was reported to accumulate in roots of alfafa plants submitted to alkali stress (Song et al., 2017) and accumulate under drought in leaves, stem, root, and root colar in Caragana korshinskii, which is described as a highly drought tolerant plant (Zhang et al., 2017). Rhanose is a pectic component that constitutes the primary cell walls (Tenhaken, 2015), its reduction under drought might indicate the loosening of cell wall polysaccharides aiming to maintain cell expansion. The metabolite similar to Isopropyl B- d-1-thiogalactopyranoside accumulated at 4DI and was not found at 12DI, this metabolite was identified based on its probable class. In SP90-3414, the metabolites threonine, digitoxose, similar to N-(3-aminopropyl)-Morpholine, 3-ureido-Propionic acid, 4-hydroxy-Benzoic acid, indole-3-acetaldehyde enol MP, and 2-(4-hydroxyphenyl)-Ethanol, were differentially abundant in all conditions. All of them accumulated at PI and decreased in response to drought. The role of most of these metabolites is poorly understood in plants. However, they belong to different classes and might be related to different biochemical process, indicating the complexity of sugarcane responses to drought. Similar to icosapentaenoic acid, beta-D-Fructofuranosyl-(2,1)-beta-D-Fructofuranose and galactose changed their abbundace only between PI and 12DI and did not changed in comparison to 4DI. The first two metabolites accumulated at 12DI while galatose reduced at 12DI, being strongly influenced by drought.

rCCA Network Correlation

Here, we applied rCCA to explore correlation structures between metabolites and proteins, in two pan class="Species">sugarcane varieties, contrasting in drought tolerance. We propose that rCCA can be used as a powerful approach to track the main changes in response to abiotic and biotic stress. It is important to mention that as well as for metabolomic correlation network and gene co-expression, metabolite-protein correlation is not always in agreement with known biochemical pathways (Fukushima et al., 2011).

The novelty in our network results is the fact that most proteins were correlated with metabolites whose function are unknown, meaning that based on the protein function in the plant metabolism, we can have some light in deciphering the biological role of such metabolites in response to drought. By this exploratory approach, networks revealed variety specific requirements to face drought. Besides, it was observed that leaf metabolism was adjusted in a different manner, in response to 4DI and 12DI. In the analysis Cpan class="Chemical">TC15-PI-4DI, the tolerant variety did not show major different with the control plants at 4 days of drought ( ). For CTC15 PI-12DI, the biggest cluster is related to Hexadecanoic acid, also known as palmitic acid (PA). This lipid was unique at 12DI and is one of the major fatty acids forming virtually all, natural lipids, which serve either as lipid background for storage fats and oils, or hydrophobic matrix of cell membranes, or components of cuticle waxes and polymers (Sidorov et al., 2014). PA was reported to accumulate in wheat (Marček et al., 2019) and soybean (Mohamed and Latif, 2017) genotypes in response to drought. PA is positively correlated with proteins involved in different biological roles. NadH-quinone oxidoreductase (SCAGAM2126G09) and cytochrome c oxidase (SCSFLR2031H06) are involved in the electron transport chain (Melo et al., 2004; Mansilla et al., 2018). Two proteins found were related to DNA and RNA metabolism: Dead-box ATP-dependent RNA helicase (SCRUHR1077B02), that is related to replication, repair, recombination, transcription, translation, ribosome biogenesis and splicing (Tujeta and Tujeta, 2004) and polypyrimidine tract-binding protein (PTB) (SCCCCL4012D10) related to alternative splicing regulation (Rühl et al., 2012). The overexpression of dead-box ATP-dependent RNA helicase, in rice, was reported to regulate transcripts involved in the abiotic and oxidative stress responses and, therefore, it functions as a molecular switch in different signaling pathways leading to stress tolerance (Wang et al., 2016; Nawaz and Kang, 2019). The ml domain (MD-2-related lipid-recognition) protein, which is implicated in lipid recognition, is also part of this cluster. Cluster 2 is formed by three metabolites from different chemical classes. Two of them, 3-Methyl-2-oxopentanoic-acid (organic acid) and O-Benzoyl-beta-D-glucoside (glucoside) are positively correlated with all proteins in the cluster. However, the compound similar to inositol (alcohol) was correlated only with peroxidase 5-like (SCBFRT1072H11) and aldehyde dehydrogenase (SCJFLR1073D12), proteins involved in cell antioxidation (Pandey et al., 2017) and detoxification (Brocker et al., 2013). Despite the lack of information about these metabolites in plant metabolism, they correlated positively with all proteins in the cluster, which seems to be involved with cell protection (antioxidant, detoxification, and protein degradation) and photosynthesis efficiency, suggesting this cluster might have a protection role. Light-harvesting chlorophyll a/b-binding proteins (LHCB) (SCJFSB1014F03.b) increased in tolerant genotypes of Zea mays (Benesova et al., 2012), and Malus domestica (Zhou et al., 2015) when exposed to drought. The LHCBs have been predicted to modulatate ROS homeostasis (Xu et al., 2012) and this role probably requires other proteins that can cope with ROS, to protect cell metabolism. On the other hand, the proteins might suggest that the three metabolites may also be involved in cell protection.

Negative correlation was observed in cluster 3. pan class="Chemical">Galactaric acid (unique at PI) and 17-beta-streadiol (not found at 12DI) were negatively correlated to wd40-like beta propeller repeat (WDR) (SCCCRZ3004G04) and dead-box ATP-dependent RNA helicase (SCUTRZ3071H01), while Galacturonic acid (not found at 12DI) positively correlated with dead-box ATP-dependent RNA helicase and negatively correlated to protoporphyrinogen chloroplastic (SCMCRZ3064A10).

Galactaric acid is an pan class="Chemical">organic acid derived from galacturonic acid (GalA). Galactaric acid was reported to decrease by salinity in Arabidopsis (Sanchez et al., 2008) and to increase in plants exposed to elevated CO2, under heat stress, in grass (Yu et al., 2012). 17-beta-stradiol is a steroid, reported to be involved in plant growth and flowering (Janeczko and Skoczowski, 2005; Lindemann, 2015). Both metabolites correlated to proteins that have a role in DNA/RNA metabolism and scaffolding proteins. These proteins might be required to regulate important cellular processes by interacting with multiple proteins, under drought. One possibility is that, due to their role in ceullar metabolism, the metabolites levels are reduced. Another possibility is that the metabolites abundance decreased in response to drought and these two proteins are required be regulate important cellular processes by interacting with multiple proteins. GalA is an important constituent of cell wall (Bethke et al., 2014) and is an intermediate in the biosynthesis of the antioxidant ascorbic acid (AsA) (Agius et al., 2003). GalA accumulated under drought in wheat genotypes (Marček et al., 2019) and bermudagrass leaves (Du et al., 2012). Authors suggested that drought could lead leaves to oxidative stress, thus GalA increased to averting oxidative damage under severe water stress (Du et al., 2012). This metabolite has not changed significantly in CTC15; however, it is probably required to the same purpose in sugarcane leaves. GalA correlated positively with protoporphyrinogen, the abundance of this protein enhanced in salt-treated maize chloroplasts, the author suggested this increment would alleviate oxidative stress in salt-stressed maize (Zhao et al., 2013). Drought leads to oxidative stress and this could be the reason for this correlation, in our study. On the other hand, GalA correlates negatively with the dead-box ATP protein.

In SP90-3414, rCCA revelad a more complex regulatory process undergoing in the drought-susceptible variety, in comparison to the tolerant one. The PI-4DI cluster 1 represents the biggest cluster and involves positive and negative correlations. Adenine, indole-3-acetaldehyde (lAAld), and octadecanal correlated with the same eleven proteins ( ). The three metabolites act in different process; adenine is one of the major purines and play important role in plant metabolism, as it represents the major energy currency of the cell (Haferkamp et al., 2011). Sukrong et al. (2012), suggested that adenine might play a part in the signals that modulate responses to abiotic stress and plant growth. lAAld can be oxidezed to the auxin indole-3-acetic acid (IAA) and is related to the tryptophan metabolim (Tsurusaki et al., 1997). Octadecanal is a lipid, found in the cuticular wax of maize (Perera and Nicolau, 2007) and epicular wax of pear fruits (Wu et al., 2018). All metabolites were negatively correlated with two chloroplastic proteins: translocase and a protease. Besides, these metabolites correlated positively with proteins related to a variety of molecular processs, such as oxidative stress (glutathione-s-transferase), photosynthetic regulation (peroxisomal-2-hydroxy-acid-oxidase), proteins degradation (26s proteasome and proteosomal subunit), signaling (14-3-3), chromating organization (histone h1), proteins chaperoning (dnaj protein), and monolignol biosynthesis (cinnamyl alcohol dehydrogenase). Octadecanoic acid, also known as stearic acid, was reported to decrease in sunflower seeds (Petcu et al., 2001) and soybean (Mohamed and Latif, 2017) submitted to drought. But it was reported to increase in the leaves of a drought-susceptible oat variety, under drought (Sánchez-Martín et al., 2018). This metabolite was positively correlated with histone h2 and negatively correlated with histone h1, peroxisomal oxidase, proteasome, 14-3-3, glutathione-s-trasnferase, and cinnamyl-alcohol dehydrogenase. Besides the complexity of this cluster, it is clear that the same set of proteins regulates or is regulated, positively and negatively, by two different sets of metabolites, suggesting the metabolites might have completely different roles under drought. The analysis for SP90-3414 PI-12DI resulted in five clusters. Cluster 1 is the most complex, involving positive and negative correlation and six metabolites: probable steroid, 3-hydroxy-Indole, trans-3-hydroxy-Proline, probable amino acid derivative, allo-Threonine, and Digitoxose. As reported to the other cluster, this one refers to a different set of metabolites, correlated with a variety of proteins, some of them related to cellular protein modification, transport, translation, response to stress, and cellular process. Two purple acid phosphatase precursor (PAP) protoforms (SCEQLR1091F02 and SCCCLR1075F03) were found in cluster 1. Both were positively correlated to 3-hydroxy-Indole. However, SCEQLR1091F02 was also correlated to the probable steroid metabolite while, SCCCLR1075F03 was correlated to trans-3-hydroxy-proline ( ). PAP is an acid phosphatase (APase) and acts in phosphorus (P) scavenging and recycling under conditions of P deficiency (Liu et al., 2016). It is known that, drought decreases the concentration of nutrient, such as P in plant tissues (reviewed in Bista et al., 2018). The identification of both proteins in SP90-3414 PI-12DI network, suggest that, at this level of stress, the susceptible variety might be dealing with P deficiency and require PAPs. In wheat exposed to drought, the authors (Sharma and Kaur, 2008) found a significant increase in APase activity in leaves and grains of drought tolerant plant, however, no change was observed in the drought sensitive. Curiosly, the SP90-3414-PI-12DI cluster 3 is composed by the same metabolites and almost the same proteins, found in the SP90-3414-PI-4DI network, suggesting that adenine, lAAld, octadecanal, as well as the proteins present in the SP90-3414-PI-12DI cluster 3 might be strongly influenced or impacted by drought, in the sugarcane susceptible variety.

Conclusion

In summary, the analysis of pan class="Species">sugarcane leaves from two drought-tolerant contrasting varieties revealed the dynamics of proteome and metabolome in response to drought. Considering the physiologycal responses to drought, this was observed earlier in the drought-susceptible variety (SP90-3414), in comparisson to the tolerant one. Proteomics revealed that different biological process had a stronger impact in each variety (e.g. translation in Cpan class="Chemical">TC15, generation of precursor metabolites and energy in SP90-3414). The differentially abundant metabolites found, in response to drought, were in agreement with previous observations, and reveals the plasticity of metabolism under drought. rCCA network revealed the complexity of protein-metabolite/metabolite-protein regulation. In addition, we exemplify that rCCA, based on protein and metabolite data, can be a powerfull strategy to identify potential networks; assign biological conclusion on metabolite-protein/protein-metabolite regulation and elucidate biological function of predicted metabolites and proteins.

Author Contributions

IB conceived of the study, carried out the experiments, analyzed the data and wrote the manuscript. pan class="Chemical">TC assisted with metabolomics and proteomics studies. FM performed the rCCA and LF assisted with proteomics data processing. CL conceived of the study, participated in its design and coordination, and reviewed the manuscript. All authors approved the final manuscript.

Conflict of Interest

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.