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 Labelfree quantitative proteomics and GC-MS-based metabolomics. We chose sugarcane (Saccharum spp.) as our model plant because this perennial grass, cultivated in tropical and subtropical 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 twodimensional 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.

Plant Material and experimental Conditions
To study changes in protein and metabolite profiles in response to drought, two sugarcane varieties were used: CTC15, 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,. 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 (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-MS e Analysis
The peptides mixture was analyzed by reverse-phase ultraperformance 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 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.

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 matched 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 MS E 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 Retch 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 (CHCl 3 ): water (H 2 0) (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 ( 13 C 3 )-myristic acid, ( 13 C 4 )-palmitic acid and ( 2 H 4 )-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  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  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.

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 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 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).

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. CTC15 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 (Figure 1). 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 (Figure 2). Leaf transpiration and 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. 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 (Supplementary Figure 1), providing an overview of whether the changes represented common or specific responses. In CTC15 and SP90-3414 plants were found 1218 and 1336 common proteins, respectively. In both varieties, 12DI plants showed the highest number of unique proteins ( Figure 3A) and 31.5% (Figure 3B) of the total variation, respectively. 4DI and 12DI plants were separated by PC2, representing 15.7% ( Figure 3A) and 18.8% (Figure 3B) 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 (Supplementary Figure 2). 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 (Supplementary Figure 3), most of them belongs to the category Metabolic process-GO:0008152 (25%) and Cellular process-GO:0009987 (17%).

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 (Supplementary Tables 4 and 5). 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 (Supplementary Tables 4 and 5). Venn diagram was performed to identify the common and unique metabolites in each variety (Supplementary Figures 4A, B), 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 (Supplementary Figures 4A, B). 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 (Supplementary Figures 4A, B). 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 (Figures 6A, B). 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 T 2 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 ( Figure 6A). 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 ( Figure 6B). 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 (Supplementary Tables 3 and 4). 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 (Figure 7). 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.

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 (Figures 8 and 9). Due to the high stringency of the analysis. 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 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 Supplementary Tables 6 and 7. The unit representation and the correlation circle plot of the rCCA can be visualized in the Supplementary Figures 5  and 6, respectively. The unit representation plot of the tolerant variety in PI-4DI analysis (Supplementary Figure 5A) 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. Supplementary  Figures 5B-D 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 Supplementary Figure  5A, 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 (Supplementary Figure 6) 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 (Supplementary Table 6). CTC15-PI-12DI cluster 1 (Figure 8) is the biggest one, refers to Hexadecanoic acid (ID 1391, positively correlated to 15 proteins). Cluster 2 refers to three metabolites: 3-Methyl-2oxopentanoic-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 (Supplementary Table 6), 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-score = 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 (Supplementary Table 7; Figures 9A, B). 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-3acetaldehyde 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 (Supplementary Table 7). 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 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 . Therefore, photosynthesis is particularly impacted to drought once stomata closure reduces CO 2 diffusion to the fixation sites, in leaf mesophylls (Ghannoum, 2009). In accordance to this, all parameters measured (A, E and g s ) were reduced at 4DI and 12DI in SP90-3414 (Figure 2). 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  , 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 (Figure 1)

Proteome Dynamics of Sugarcane Leaves Under Drought
It is known that in C 4 plants photosynthetic inhibition under drought is mainly due to stomatal closure (Ghannoum, 2009), resulting in lower carbon fixation, followed by the overreduction of the electron transport system and production of ROS (Nouri et al., 2015). Although physiological data showed that tolerant variety, at 4DI (Figure 1), 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 (Supplementary  Table 3). 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, 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 . 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 accordance 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-1thiogalactopyranoside 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 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 CTC15-PI-4DI, the tolerant variety did not show major different with the control plants at 4 days of drought (Suplementary Figure 4A). 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 ATPdependent 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  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.
Galactaric acid is an 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 . 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 . Authors suggested that drought could lead leaves to oxidative stress, thus GalA increased to averting oxidative damage under severe water stress . 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 ( Figure  9A). 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 droughtsusceptible 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 cinnamylalcohol 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-3hydroxy-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 ( Figure 9B). PAP is an acid phosphatase (APase) and acts in phosphorus (P) scavenging and recycling under conditions of P deficiency . 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 sugarcane leaves from two droughttolerant 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 CTC15, 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/metaboliteprotein 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. 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.