Integrative leaf–root metabolomics and root proteomics reveal proline-associated mechanisms underlying drought tolerance in RRIT 251 versus RRIM 600 rubber trees

The rubber tree (Hevea brasiliensis) serves as the primary source of natural rubber (NR), a critical raw material for multiple industries, including the automotive, medical, and manufacturing sectors. Approximately 90% of the world's natural rubber originates from plantations in Southeast Asia, making this species economically vital for countries such as Thailand, Indonesia, and Malaysia^[1]. In addition to its industrial value, rubber cultivation provides livelihoods for millions of smallholder farmers in these regions, underscoring its significant socioeconomic importance. However, climate variability, particularly the increasing frequency of drought events, threatens the sustainability of rubber plantations. Prolonged drought can severely disrupt physiological processes such as photosynthesis and latex production, leading to reduced yields and long-term economic losses^[2,3]. In addition, rubber trees take more than 7 years to mature before harvesting can begin, so the survival and growth of seedlings and grafted plants are critical. This makes the development of water stress-tolerant rootstock clones a necessity for the industry. Understanding the biological mechanisms that enable certain rubber tree clones to maintain growth and productivity under water-limited conditions is therefore essential for improving drought resilience.

Investigations into these adaptive strategies in H. brasiliensis have largely centered on leaf-level responses, revealing significant clonal variation. Tolerant clones like RRIT 251 exhibit rapid stomatal closure to conserve water, thereby maintaining higher leaf water potential and protecting the photosynthetic systems from damage^[4]. This physiological resilience is supported by a robust molecular defense, including the accumulation of osmoprotectants like proline and the upregulation of antioxidant enzyme systems to mitigate oxidative stress^[5]. Gene expression studies have identified key transcription factors such as LEA, NAC, and MYB that control the expression of stress-responsive genes and orchestrate this defense^[6−8]. Although these mechanisms have been characterized in model plants such as Arabidopsis thaliana and rice (Oryza sativa)^[7,9], their roles in H. brasiliensis remain less defined compared with model plants, limiting the development of targeted strategies to improve drought resilience. This knowledge gap limits the development of targeted strategies to improve drought resilience in rubber trees, necessitating comprehensive studies that integrate physiological, metabolomic, and proteomic analyses.

Although previous research has documented some physiological and molecular responses of H. brasiliensis to drought, such as stomatal regulation, antioxidant enzyme activity, and osmolyte accumulation^[2,10], these studies have largely focused on leaf tissues. Few investigations have explored integrated metabolic and proteomic responses, particularly in root tissues, which are the primary sites for water sensing and uptake. Roots play a pivotal role in plants' adaptation to drought by modulating water transport, activating oxidative stress defenses, and synthesizing osmoprotectants^[11]. Moreover, root-derived signals such as pipecolic acid and azelaic acid are increasingly recognized as systemic mediators of stress tolerance, linking local root responses to whole-plant adaptations^[12,13]. Despite their central role, root-specific biochemical and proteomic responses to drought in rubber trees remain poorly characterized. This gap limits the identification of robust molecular markers that could be leveraged for the selection or engineering of drought-tolerant clones. Furthermore, the dynamic nature of drought responses, which vary over time and between different stress intensities, requires time-resolved and tissue-specific omics investigations to capture the full spectrum of adaptive mechanisms.

Multi-omics strategies that combine metabolomics and proteomics provide powerful tools for unraveling the complex biochemical networks underlying plant stress responses. Metabolomics enables the comprehensive profiling of small molecules that function as osmoprotectants, antioxidants, and signaling mediators, providing insights into the biochemical adjustments that support stress tolerance. Proteomics complements this by identifying changes in protein abundance and post-translational modifications that regulate biological processes, such as seed germination under microgravity, plant development, and metabolic pathway reprogramming^[14−16]. Together, these approaches provide a holistic view of plants' stress adaptations, linking molecular changes to physiological outcomes. In the context of H. brasiliensis, applying integrated metabolomic and proteomic analyses to compare drought responses in clones with contrasting stress tolerance profiles has the potential to reveal key biochemical and molecular determinants of resilience. Such findings could inform breeding programs by identifying candidate metabolites and proteins that serve as biomarkers for drought tolerance, ultimately contributing to the development of climate-resilient rubber plantations.

Grafted rubber, consisting of RRIM 600 or RRIT 251 scions grafted onto unselected seedling rootstocks, were utilized. This configuration allows for the investigation of how specific scion genotypes systemically modulate the biochemical and proteomic responses of the root system under drought stress. Healthy clones were obtained from a local nursery in Thung Tao Mai Subdistrict, Na San District, Surat Thani Province, operating under the supervision of the Rubber Authority of Thailand. All grafted rubber trees were six months old and derived via clonal propagation to ensure genetic uniformity. To ensure that the omics data reflected the root system's response to drought as influenced by the scion variety, root tissues were harvested directly from the rootstock portion of the grafted assembly. The plants were cultivated under greenhouse conditions with controlled environmental parameters: The temperature was maintained at 30–35 °C and relative humidity at 60%–70%, with a 12-h photoperiod (light/dark). Prior to the onset of the drought treatment, all plants were acclimated under well-watered conditions, equivalent to 85%–100% of field capacity (the amount of water retained in the soil after it has been thoroughly saturated with water and then allowed to drain freely under gravity until no further drainage occurs) for two weeks to establish their baseline physiological and biochemical status.

The experimental design included three soil moisture regimes to simulate progressive drought conditions as shown in Fig. 1. These were (1) well-watered at 85%–100% field capacity (control), (2) moderate drought stress at 50% field capacity, and (3) severe drought stress at 30% field capacity^[17]. Soil water content was continuously monitored using the weighing method and adjusted daily to maintain the target levels. Each pot was irrigated daily at 08:30 a.m. to maintain the designated conditions.

Figure 1. 

Experimental design for drought stress treatments in H. brasiliensis. Two rubber tree clones, RRIM 600 and RRIT 251, were subjected to three soil moisture regimes: Well-watered at 85%–100% field capacity (control), moderate drought stress at 50% field capacity, and severe drought stress at 30% field capacity. Soil water content was monitored continuously using the gravimetric method and adjusted daily to maintain the target water levels.

Sampling was performed at the long-term stress time point, which was identified by the occurrence of leaf abscission (leaf drop), representing prolonged water deficit. For each rubber tree clone, two biological replicates in total were collected per treatment condition. All harvested leaf and root tissues were flash-frozen in liquid nitrogen immediately upon collection and stored at –80 °C until proteomic and metabolomic analyses.

Photosystem II efficiency analysis or Fv/Fm analysis

Chlorophyll fluorescence analysis was performed to evaluate the efficiency of Photosystem II (PSII) in the leaves. The leaf samples were used for analyzing chlorophyll florescence, which was measured with a portable plant efficiency analyzer (Hansatech Instruments, UK), which provides high-resolution detection of fast fluorescence kinetics. The samples were tested for Fv/Fm during midday (11:00 a.m.−1:00 p.m. solar time) to ensure consistent light conditions and minimize diurnal variation in photosynthetic activity^[18]. Prior to the measurement of chlorophyll fluorescence, the leaves were dark-adapted for 30 min using clips to allow complete oxidation of PSII reaction centers. Dark adaptation ensured that all PSII reaction centers were open, providing a reliable baseline for the measurement of minimum fluorescence (Fo). After dark adaptation, a saturating pulse of actinic light was applied to induce maximum fluorescence (Fm). The variable fluorescence (Fv) was calculated as the difference between Fm and Fo (Fv = Fm – Fo). The maximum photochemical efficiency of PSII was then expressed as the ratio Fv/Fm^[19].

Metabolite extraction and sample preparation for liquid chromatography–tandem mass spectroscopy analysis

Leaf and root samples of RRIM 600 and RRIT 251 were harvested and immediately stored at –80 °C until metabolite extraction. During transfer and processing of the samples, tissues were kept below 4 °C to prevent metabolite degradation. Frozen tissues were ground to a fine powder in liquid nitrogen, and 100 mg of each sample was weighed for extraction. Metabolite extraction was performed following a previously described protocol with minor modifications^[20]. Briefly, 100 mg of powdered leaf or root tissue was mixed with 100% methanol and incubated at 4 °C for 12 h, followed by an additional extraction step at –20 °C for 48 h to enhance the metabolites' solubilization. After incubation, the extracts were centrifuged at 14,000 × g for 20 min at 4 °C to pellet insoluble material. The resulting supernatants were subjected to solid-phase extraction (SPE) using an extraction manifold system. The SPE cartridges were preconditioned with 25 mL of acetonitrile and equilibrated with 60 mL of water. The clarified methanol extracts were loaded onto the SPE columns and eluted with 90% acetonitrile in water. The collected eluates were evaporated under a vacuum using rotary evaporation to remove solvents. The dried extracts were reconstituted in 1 mL of methanol and subsequently diluted 1:10 (v/v) with water, yielding a final volume of 10 mL. This prepared solution was used for liquid chromatography–tandem mass spectroscopy (LC-MS/MS) analysis. All extractions were conducted in two biological replicates per treatment to ensure reproducibility and confirm phytochemical consistency across the samples.

LC-MS/MS instrumentation and acquisition parameters for untargeted metabolomics

The untargeted metabolomics analysis was performed using an UltiMate 3000 UHPLC system coupled to a Thermo Q Exactive-X Quadrupole Orbitrap Mass Spectrometer (Thermo Fisher Scientific, Waltham, MA, USA). Metabolite separation was achieved on a Hypersil GOLD™ C18 column maintained at 60 °C. A sample injection volume of 5 µL (0.3 µg/µL) was applied with a constant flow rate of 0.35 mL/min. The auto-sampler temperature was controlled at 8 °C to preserve the samples' integrity. Two mobile phase systems were used, depending on the ionization mode. For the positive mode, the mobile phases consisted of 20:80 methanol/water with 0.1% formic acid (Mobile Phase A) and acetonitrile with 0.1% formic acid (Mobile Phase B). For negative mode, 20:80 methanol/water with 0.1% acetic acid (Mobile Phase A) and acetonitrile with 0.1% acetic acid (Mobile Phase B) were used. The gradient started at 95% A and 5% B, held for 2 min, followed by a linear increase to 55% B over 6 min. The column was flushed with 99% B for 5 min and returned to the initial conditions for re-equilibration, resulting in a total run time of 20 min per sample. Blank injections (0.1% formic acid in water) were included between runs to prevent carryover. The mass spectrometer operated in both positive and negative electrospray ionization (ESI) modes with spray voltages of 4.0 and 3.8 kV, respectively. The sheath gas and auxiliary gas flow rates were set to 48 and 11 arbitrary units (AU), and the capillary temperature was maintained at 350 °C. Data acquisition alternated between full mass spectroscopy (MS) scans and data-dependent tandem MS (MS/MS) scans using dynamic exclusion. Full MS scans were acquired over a m/z range of 75–650 with a resolution of 120,000, an automatic gain control (AGC) target of 3 × 10⁶, and a maximum injection time of 10 ms. MS/MS scans were performed at a resolution of 15,000 with variable AGC targets and injection times depending on the ion intensity. The top 10 most intense precursor ions were selected for fragmentation in each cycle. Data acquisition was controlled using Xcalibur 3.1 software (Thermo Fisher Scientific, Waltham, MA, USA).

Data processing and metabolite quantification

Raw LC-MS/MS spectral data were processed using Compound Discoverer 3.1 (Thermo Fisher Scientific) with a default untargeted metabolomics workflow, modified to optimize feature detection and annotation. Retention time alignment was performed with a maximum deviation threshold of 0.4 min, and mass accuracy was maintained within a tolerance of 5 parts per million (ppm) tolerance. The signal-to-noise ratio threshold was set at 1.5 to ensure reliable peak detection. Intensity tolerance was set at 30%, and the minimum peak intensity threshold was defined as 500,000 to filter out low-abundance noise. QC-based area correction using a cubic spline regression model was applied to normalize sample-to-sample variability. Detected features were integrated to predict molecular formulas, which were subsequently matched against online databases including mzCloud (www.mzcloud.org, accessed 01/05/2025) and ChemSpider (www.chemspider.com, accessed 01/05/2025), referencing the ChEBI and ChemBank repositories for compound identification. Among the candidate metabolites, the annotations with the highest MS/MS fragmentation coverage scores with profiling were selected as the most confident identifications. Processed metabolite abundance data from leaf and root tissues were averaged across technical replicates for each experimental group, defined by cultivar (RRIM 600 or RRIT 251) and soil moisture condition (85%, 50%, or 30% field capacity). Prior to statistical analysis, abundance values were normalized using z-score transformation to standardize the data across samples. Unsupervised hierarchical clustering was conducted using average linkage and Euclidean distance metrics. The resulting clusters were visualized as heatmaps using the Seaborn package (v0.11.2) in Python. Euclidean distances between experimental conditions were converted to relative percentage differences to quantify the magnitude of metabolic shifts under severe drought stress.

Protein extraction, digestion, and sample preparation for root proteomics analysis

Protein extraction from root tissues of RRIM 600 and RRIT 251 was performed following established protocols with minor modifications^[21,22]. Briefly, 100 mg of frozen root tissue was homogenized in a lysis buffer containing 1% sodium dodecyl sulfate (SDS), 0.5 mM Tris (2-carboxyethyl) phosphine (TCEP), 10 mM NaCl, and 2 mM ethylenediaminetetraacetic acid (EDTA) in 50 mM 4-(2-Hydroxyethyl)piperazine-1-ethanesulfonic acid, N-(2-hydroxyethyl)piperazine-N′-(2-ethanesulfonic acid) (HEPES)–KOH (pH 8.0), supplemented with a protease inhibitor cocktail. Tissue homogenization was carried out using a tip-probe sonicator with 5-s bursts on ice to minimize thermal degradation. The lysate was centrifuged at 12,000 × g for 20 min at 16 °C to separate soluble proteins from cellular debris. The supernatant containing the extracted proteins was collected, and the protein concentration was quantified using the BCA Protein Assay Kit (Pierce, Thermo Fisher Scientific, USA). Protein precipitation was performed by incubating the supernatant at –20 °C for 16 h in a 1:5 (v/v) mixture of 16% trichloroacetic acid (TCA) in acetone. Precipitated proteins were collected by centrifugation at 10,000 × g for 10 min at 4 °C. The resulting protein pellets were resuspended in 0.5% RapiGest SF surfactant (Waters Corporation, Milford, MA, USA), 10 mM NaCl, and 10 mM HEPES–KOH (pH 8.0), followed by a second sonication step to ensure complete solubilization. The protein concentration was remeasured and adjusted to 50 µg per sample. Reduction and alkylation were performed by adding 1 mM TCEP and incubation at 90 °C for 15 min, followed by cooling to room temperature and the addition of 5 mM iodoacetamide (IAA). Samples were incubated in the dark at room temperature for 50 min. Desalting was performed using Zeba spin desalting columns, and the flow-through was collected and mixed with 0.5% RapiGest (1:1 v/v). Tryptic digestion was carried out by adding 100 ng/µL trypsin (at an enzyme-to-substrate ratio of 1:75 w/w) and incubation at 37 °C for 5 h. Digestion was terminated by adding 1% formic acid (1:10 v/v), and the resulting peptides were lyophilized and stored at –80 °C until LC-MS/MS analysis.

LC-MS/MS acquisition and protein quantification analysis

Tryptic peptides were analyzed using an Orbitrap HF-X mass spectrometer (Thermo Fisher Scientific) coupled to an EASY-nLC nano-LC system equipped with an EasySpray C18 column (Thermo Scientific ES903; 75 µm × 15 cm, 2.0 µm particle size). Chromatographic separation was performed using a linear gradient of 4% to 45% Mobile Phase B at a flow rate of 300 nL/min. Mobile Phase A consisted of 0.1% formic acid in water; Mobile Phase B consisted of 99% acetonitrile with 0.1% formic acid. The mass spectrometer operated in positive ionization mode using data-dependent acquisition (DDA). Full MS scans were acquired over an m/z range of 400–1,300 with a resolution of 120,000, an AGC target of 3 × 10⁶ ions, and a maximum injection time of 15 ms. The top precursor ions meeting the AGC target threshold were selected for MS/MS fragmentation using an isolation window of 1.4 m/z. MS/MS scans were acquired at a resolution of 15,000 with a maximum injection time of 25 ms. All spectra were acquired using Xcalibur software (Thermo Fisher Scientific). Raw data files were processed using Proteome Discoverer software with the SEQUEST HT search engine, referencing the UniProtKB H. brasiliensis database (taxonomy: 3,981, 44,374 sequences, retrieved 15/05/2025)^[23]. The search parameters included strict trypsin specificity, allowing up to two missed cleavages. Carbamidomethylation of cysteine residues (+57.0215 Da) was set as a fixed modification, and oxidation of methionine (+15.9949 Da) was set as a variable modification. Protein quantification was normalized using the software's default normalization algorithm to account for technical variability across samples.

Statistical analysis

For the analysis of physiological parameters, including Fv/Fm, statistical comparisons among treatments within each cultivar were performed using one-way analysis of variance (ANOVA). Significant differences were accepted at p-values less than 0.05. Data were presented as the means ± standard deviations (SD) based on two biological replicates (n = 2) per treatment condition. For the pairwise comparisons during the proteomic analysis, we performed one-way ANOVA at the protein level with two multiple testing correction methods, including Bonferroni and Benjamini–Hochberg false discovery rate (FDR) corrections. For metabolomics, significant features were defined according to an alpha level of 0.05, with subsequent validation through FISh profiling to ensure chemical and biological reliability.

Visual leaf phenotyping revealed progressive alterations in color intensity and distribution patterns in response to drought stress in both rubber tree clones. Under well-watered control conditions, the leaves of RRIM 600 and RRIT 251 maintained a uniformly green color with no visible signs of chlorosis or dehydration, indicating optimal physiological status and intact chlorophyll integrity. As water availability declined to 50% field capacity, RRIM 600 exhibited early signs of leaf yellowing, particularly along the margins and interveinal regions. This symptom suggests the partial degradation of chlorophyll and redistribution of nitrogenous compounds as a result of drought-induced metabolic shifts. In contrast, RRIT 251 exhibited a more uniform green coloration at the same stress level, suggesting a delayed onset of chlorotic symptoms. Under the most severe drought condition, both clones exhibited marked leaf yellowing; however, the extent and distribution differed. RRIM 600 displayed widespread chlorosis, curling at the leaf edges, and loss of turgor, indicative of advanced cellular dehydration and impaired pigment synthesis. RRIT 251, though also showing yellowing, retained comparatively greener central lamina regions. These visual differences provide qualitative evidence that RRIT 251 has a greater capacity to delay senescence-related leaf damage under drought stress.

Photosystem II efficiency measurement by Fv/Fm

Under drought stress conditions, a progressive decline in Fv/Fm values was observed in both rubber tree clones (Table 1), indicating impaired photosynthetic efficiency. In RRIM 600, the Fv/Fm value under well-watered conditions was 0.808 ± 0.005, which significantly decreased to 0.792 ± 0.005 under moderate drought and further dropped to 0.665 ± 0.008 under severe drought, with the latter being highly significant (p < 0.01). This pattern suggests that RRIM 600 exhibits a marked sensitivity in PSII's activity under increasing dehydration stress.

Similarly, RRIT 251 showed a control Fv/Fm value of 0.823 ± 0.005, which slightly decreased to 0.788 ± 0.035 under 50% field capacity without statistical significance, and further declined to 0.658 ± 0.014 at 30% field capacity with a significant difference (p < 0.01). These results indicate that although both clones experience reductions in Fv/Fm under water deficit, RRIT 251 maintains relatively higher photochemical efficiency under severe stress compared with RRIM 600.

Multivariate and clustering analysis reveal clone-specific metabolic reprogramming

Principal component analysis (PCA) based on replicate-level metabolite abundance in leaf tissues revealed the distinct separation of sample groups according to soil moisture levels and clonal identity. The first two principal components (PC1 and PC2) accounted for a substantial portion of the total variance, enabling effective discrimination among treatments. Samples under severe drought formed a tight cluster, clearly separated from moderate drought and well-watered control conditions (Fig. 2a). This pattern was consistent for both clones. RRIM 600 samples under 30% drought exhibited greater dispersion along PC1, whereas RRIT 251 samples showed more compact clustering, indicating a more uniform metabolic profile under stress.

Figure 2. 

Multivariate and hierarchical clustering analysis of metabolite profiles in leaf and root tissues of RRIM 600 and RRIT 251 under different soil moisture conditions (85%, 50%, and 30% field capacity). (a) PCA of leaf metabolite profiles and root metabolite profiles (b). The x-axis and y-axis represent the first and second principal components, explaining the highest proportions of total variance in the dataset. Clustering by drought severity and cultivar is observed. Hierarchical clustering heatmap of (c) leaf metabolites and (d) root metabolites. Columns represent treatment groups labeled by field capacity and cultivar; rows represent individual metabolites. Clustering was based on average linkage and Euclidean distance. Color gradients reflect standardized (z-score) abundance values across samples.

PCA of the root metabolite profiles also demonstrated condition- and genotype-dependent clustering (Fig. 2b). Samples under the 30% field capacity treatment were separated from those at 50% and 85% along PC1, though the overall dispersion was lower compared with the leaf PCA. The separation between RRIM 600 and RRIT 251 was evident, with RRIT 251 root samples clustering more tightly under drought stress. Replicates were consistent within groups, confirming the reliability of technical injections and the stability of metabolic signatures under each experimental condition. To complement the PCA, hierarchical clustering heatmaps based on Euclidean distance further clarified the global structure of metabolites' abundance profiles across treatments (Fig. 2c, d). In leaf tissues (Fig. 2c), metabolic patterns clustered primarily according to drought severity, with the 30% field capacity samples forming a distinct group, particularly in RRIT 251. The clustering in roots (Fig. 2d) was even more differentiated, where RRIT 251 under 30% drought formed a clearly separated branch from both its 50% and 85% field capacity counterparts. Quantitatively, the Euclidean distance between RRIT 251 roots under 30% and 50% field capacity was approximately 119% greater than the distance between those at 30% and 85%, whereas the distance between RRIT 251 and RRIM 600 at 30% field capacity was 76% higher than that between RRIT 251 at 30% and 85%. These relative distances underscore the uniquely altered metabolic landscape of RRIT 251 roots under severe drought, reflecting early clone-specific reprogramming in response to water deficit. Given that the proteomics analysis was performed only in the roots, the higher Fv/Fm values in RRIT 251 (Table 1) are attributed to its stronger root antioxidant/osmoprotectant role rather than unmeasured leaf proteome changes.

Differential abundance of identified metabolites in leaves and roots under drought stress

Prior to the comparative analysis, data normalization was performed to account for technical variation in sample injection and ionization efficiency, ensuring that the metabolites' intensities were comparable across all biological replicates. The effectiveness of this normalization is illustrated by the distribution of total metabolite intensities, which are presented as box plots in Supplementary Fig. S1. Comparative metabolomic analysis under 30% field capacity revealed widespread changes in the metabolites\ abundance profiles in both leaf and root tissues.

In the leaves of RRIM 600, out of 358 identified metabolites, 34 features (9.5%) were significantly upregulated and 71 features (19.8%) were significantly downregulated compared with the well-watered control (Fig. 3a). In contrast, the leaves of RRIT 251 exhibited a more balanced metabolic shift, with 54 features (15.1%) upregulated and 43 features (12.0%) downregulated (Fig. 3b). In root tissues, the RRIM 600 cultivar showed a more constrained response, with 74 out of 524 metabolites (14.1%) significantly upregulated and 38 metabolites (7.3%) downregulated (Fig. 3c). Meanwhile, the RRIT 251 roots demonstrated the highest degree of metabolic reprogramming, with 87 metabolites (16.6%) significantly upregulated and 48 metabolites (9.2%) downregulated under drought stress (Fig. 3d). These global changes highlight the differential capacity of each cultivar and tissue type to modulate metabolic pathways in response to water limitation. Notably, RRIT 251 roots showed both the highest number and percentage of upregulated metabolites, reflecting an active and coordinated drought-adaptive metabolic response.

Figure 3. 

Differential abundance of metabolites under drought stress relative to the well-watered controls (85% field capacity) in RRIM 600 and RRIT 251. (a) Leaves of RRIM 600, (b) leaves of RRIT 251, (c) roots of RRIM 600, and (d) roots of RRIT 251. Each panel displays metabolite differences, with log2 fold change plotted on the x-axis and –log10(p) on the y-axis.

Metabolite-level comparisons further revealed genotype-specific markers. In the leaves of RRIT 251, the dominant upregulated compounds included NP-005013 (log₂ fold change [FC] = +6.28, p = 0.0003; ranked first), 4-carboxy nevirapine (log₂FC = +5.46, p = 0.0003; ranked second), and (2R)-2-hydroxy-3-(phosphonooxy) propyl (9Z)-9-nonadecenoate (log₂FC = +4.31, p = 0.007; ranked third). The RRIM 600 root, in contrast, exhibited downregulation in a broader array of compounds, consistent with metabolic suppression under drought. In the root tissue of RRIT 251, the most strongly upregulated metabolite was proline (log₂FC = +4.35, p = 0.0081; ranked first), followed by 3-(sulfooxy) butanoic acid and sphondin, confirming an intense stress-responsive shift. RRIM 600 roots, by contrast, displayed relatively fewer significant changes.

Characterization of the top-ranking metabolites under drought stress revealed several genotype- and tissue-specific chemical markers. Table 2 presents the molecular details of these compounds, including their ionization mode, chemical formula, mass accuracy, and fragmentation evidence. This table supports the reliability of metabolite annotation through precise mass measurements and FISh profiling.

Table 2.  Top-ranking metabolites identified in rubber leaf and root tissues under drought stress (30% field capacity) with molecular annotation details.

MS mode	Name	Formular	Mass error (ppm)	FISh profiling
Leaf
(−)	NP-005013	C16H12O5	0.09
(−)	4-Carboxy nevirapine	C15H12N4O3	−4.25	−
(−)	(2R)-2-Hydroxy-3-(phosphonooxy)propyl (9Z)-9-nonadecenoate	C22H43O7P	−0.13
Root
(+)	proline	C5H9NO2
(−)	3-(Sulfooxy)butanoic acid	C4H8O6S	−0.25
(+)	Sphondin	C12H8O2
The table lists the ionization mode, metabolite name, chemical formula, and annotation mass error (ppm). FISh profiling is noted for compounds with confirmed fragmentation spectra supporting structural identification. In the FISh profiling '–' denotes compounds for which FISh fragmentation data could not be assigned because of database or adduct limitations.

Using root proteomics analysis to elucidate differential mechanisms in RRIM 600 and RRIT 251 under drought stress

Proteomic profiling of RRIM 600 and RRIT 251 roots under progressive water deficit established a robust dataset of drought-responsive proteins. In total, 785 proteins were identified across the well-watered control and the moderate and severe drought treatments, of which 642 proteins detected with at least two unique peptides were retained for quantitative analysis. Quality control assessments of the proteolytic digestion efficiency (Supplementary Fig. S2) and MS performance confirmed their analytical reproducibility (Supplementary Fig. S3). Differential abundance analysis between 50% and 30% field capacity in RRIM 600 yielded the most significantly up- and downregulated proteins (Table 3); analogous comparisons in RRIT 251 produced its top differentially regulated proteins (Table 4). These curated protein lists informed the subsequent comparative evaluation of biochemical mechanisms differentiating the two rubber tree clones under drought stress.

Under moderate drought, RRIT 251 exhibited a pronounced induction of antioxidative enzymes, with glutathione transferase increasing 6.27-fold (p = 0.0006) and (S)-2-hydroxy-acid oxidase rising 5.58-fold (p = 0.0036), whereas RRIM 600 primarily upregulated structural and metabolic proteins, notably a von Willebrand factor A (VWFA) domain-containing protein (100-fold; p < 0.0001) and beta-glucosidase (2.99-fold; p = 0.0461). Both clones shared a modest downregulation of peroxidase activity, suggesting a common reduction in hydrogen peroxide scavenging under water limitation. These divergent proteomic signatures under moderate stress highlight RRIT 251's early mobilization of detoxification pathways versus RRIM 600's focus on cell-wall and metabolic restructuring. These patterns under moderate drought set the stage for evaluating proteomic adjustments under more severe stress.

Under severe drought, RRIT 251 sustained a high abundance of stress-responsive proteins, with glutathione transferase further increasing to 7.85-fold (p = 0.0002) and annexin rising 4.17-fold (p = 0.0295), whereas RRIM 600 showed widespread downregulation of redox-related and transport proteins, including NADP-dependent oxidoreductase (0.22-fold; p = 0.0196), NAD(P)H dehydrogenase (0.01-fold; p < 0.0001), and the tonoplast intrinsic protein (0.01-fold; p < 0.0001). Both clones continued to repress tonoplasts' intrinsic protein levels, indicating a shared impairment of vacuolar water-handling under severe drought. Collectively, these results demonstrate that RRIT 251 mobilizes a more robust antioxidative and membrane-stabilizing proteomic network than RRIM 600 during progressive water deficit.

Drought stress imposes multiple physiological constraints on plants, particularly affecting leaf pigmentation and photosynthetic performance. In this study, the two rubber tree clones showed distinct visual and photochemical responses during progressive dehydration, reflecting different levels of drought tolerance. Under well-watered conditions, both clones maintained uniformly green leaves and optimal Fv/Fm values, indicating stable chlorophyll content and maximal PSII efficiency. These baseline measurements provided a reference for evaluating stress-induced changes.

As soil moisture declined to 50% field capacity, phenotypic differences became evident. RRIM 600 developed marginal yellowing and interveinal chlorosis, indicating chlorophyll degradation and nutrient remobilization. This was accompanied by a significant reduction in Fv/Fm, suggesting early impairment of photochemical efficiency and increased susceptibility to oxidative stress. In contrast, RRIT 251 retained greener leaves and showed only a slight, statistically nonsignificant decrease in Fv/Fm. The maintenance of pigmentation together with stable PSII efficiency under moderate stress implies more effective protective mechanisms against dehydration-induced photoinhibition.

Under severe drought, both clones exhibited extensive yellowing, although symptoms were more severe in RRIM 600. This clone also displayed curling leaf margins and visible tissue dehydration. These symptoms coincided with a sharp decline in Fv/Fm to below 0.67, indicating serious PSII impairment and potential irreversible damage to the thylakoid structures^[24]. Conversely, RRIT 251 preserved green pigmentation in its central leaf areas and maintained significantly higher Fv/Fm values, suggesting sustained energy conversion efficiency under water deficit. This resilience may arise from stronger antioxidant capacity, enhanced photoprotective processes such as nonphotochemical quenching, or improved osmotic adjustment that mitigates dehydration stress. The visual phenotype therefore provides an early indicator of stress response and complements Fv/Fm measurements by enabling a rapid assessment of physiological status. Although the root tissues analyzed here are genetically derived from seedling rootstocks, the pronounced differences in their molecular profiles, which were characterized by massive proline accumulation and GST induction in RRIT 251 assemblies, reflect a scion-driven systemic signaling mechanism. This indicates that the superior drought tolerance of RRIT 251 involves not only leaf-level resilience but also the ability to orchestrate protective metabolic reprogramming in the rootstock.

The metabolomic reprogramming observed in RRIT 251 under drought stress highlights a clone-specific adaptive mechanism that is particularly pronounced in root tissues. Among the significantly accumulated metabolites, proline emerged as the most strongly induced compound. Proline functions as a compatible osmolyte, a reactive oxygen species (ROS) scavenger, and a protein stabilizer, playing a critical role in maintaining cellular homeostasis during abiotic stress^[25]. Its accumulation under drought has been reported in several model plants, including Arabidopsis thaliana, Oryza sativa, and Glycine max, where it correlates with increased stress resilience and upregulation of the P5CS gene, a key enzyme in proline biosynthesis^[26,27]. In addition to proline, RRIT 251 root metabolomes showed higher accumulation of pipecolic acid, azelaic acid, and raffinose. Pipecolic acid is a lysine-derived metabolite implicated in systemic acquired resistance and long-distance signaling in plants, particularly through its downstream conversion to N-hydroxy-pipecolic acid, a key systemic defense signal in Arabidopsis^[13]. Azelaic acid, a lipid-derived signaling molecule, has been demonstrated to function as a systemic priming agent, amplifying stress signaling, lignin biosynthesis, and defense gene expression under abiotic stress^[28]. Meanwhile, raffinose, a known osmoprotectant and antioxidant, contributes to membrane stabilization and scavenging of hydroxyl radicals, with its drought-associated accumulation being extensively documented in Medicago truncatula and Zea mays^[29−31]. In contrast, RRIM 600 displayed a less dynamic metabolomic response. Though minor increases in some secondary metabolites were observed, the lack of any significant induction of known drought-related metabolites such as proline, raffinose, or signaling acids suggests limited activation of stress-responsive metabolic pathways. This subdued biochemical response may explain the physiological vulnerability of RRIM 600 under water-limited conditions. In addition to the top-ranking metabolites highlighted above, several other compounds with potential roles in drought adaptation were also identified (Table 5), indicating broader metabolic adjustment in root metabolism under water-deficit conditions.

Additional metabolites with antioxidant or osmoprotectant properties, including puerarin, coumarin derivatives, and quercetin 3-isobutyrate, were significantly elevated in RRIT 251 (Table 5). These compounds have been reported to enhance antioxidant enzyme activities and inhibit lipid peroxidation^[32,33]. Their consistent accumulation across drought treatments suggests coordinated metabolic reprogramming that supports RRIT 251's drought resilience. Puerarin is known for its capacity to activate antioxidant enzymes such as superoxide dismutase and catalase, enhances cellular resilience against ROS accumulation^[34]. This highlights puerarin as a drought-responsive antioxidant metabolite. Moreover, 6,7,8-trimethoxy-2H-chromen-2-one, a methoxylated coumarin derivative, also displayed strong accumulation under stress, particularly in the 50% RRIT 251 sample, with a fold change of 8.684 Although the RRIM 600 response at 30% field capacity was moderate, the overall trend aligns with coumarins' capacity to inhibit lipid peroxidation and scavenge free radicals^[35]. Such elevation supports the functional role of coumarins as secondary metabolites contributing to drought adaptation. Similarly, quercetin 3-isobutyrate, a flavonoid ester, exhibited remarkable accumulation under moderate drought in the 50% RRIT 251 replicate, with a 13.040-fold increase. Flavonoids such as quercetin are well-documented for their ability to quench singlet oxygen, regulate antioxidant pathways, and resist pathogens^[36,37]. Choline serves as a precursor for glycine betaine biosynthesis, a widely recognized osmoprotectant^[38]. Although its ratio was modest in 30% RRIT 251 (0.944), it showed notable accumulation in the 50% RRIM 600 sample (1.209). The presence of choline indicates metabolic flux toward compatible solute synthesis, which contributes to osmotic balance under drought conditions^[39]. Collectively, these findings provide strong evidence that drought stress in root tissues induces both antioxidant and osmoprotectant metabolite biosynthesis, with coordinated metabolic responses across both genotypes.

Root proteomic profiling in RRIM 600 and RRIT 251 revealed distinct drought-responsive patterns. In RRIM 600's roots, the most upregulated protein under drought (both 50% and 30% field capacity) was a VWFA domain-containing protein (100-fold increase), suggesting a role in stress-responsive protein interactions or signaling. Proteins harboring VWFA domains in plants, such as certain RING E3 ubiquitin ligases (e.g., MdRGLG3 in apple [Malus domestica]), have been implicated in abiotic stress tolerance; notably, an allelic variant of an apple RGLG with a VWFA domain conferred enhanced salt and alkali resistance^[40]. Another highly induced protein in RRIM 600 was a β-glucosidase (roughly threefold upregulation under moderate drought), which may facilitate rapid drought signaling by releasing active abscisic acid (ABA) from its stored glucose-conjugated form, thereby rapidly amplifying ABA signaling^[41]. Such activation of ABA would promote stomatal closure and other protective responses early in drought stress, aligning with the need for prompt osmotic regulation in RRIM 600 roots under water deficit. In RRIT 251, a glutathione S-transferase (GST) was among the most upregulated root proteins (a six- to eightfold increase, significant at 30% field capacity). GSTs play pivotal roles in plants' drought tolerance by detoxifying reactive electrophiles and mitigating oxidative damage, as they catalyze the conjugation of glutathione (GSH) to toxic stress-induced compounds^[42] Consistently, several GST genes are strongly induced in the roots during drought stress in other species, enhancing the scavenging of ROS and maintenance of redox balance^[42]. The high induction of a GST in RRIT 251 indicates stronger antioxidant system activation in this drought-tolerant clone. Additionally, an annexin was markedly upregulated (roughly fourfold) in RRIT 251 under severe drought. Annexins are Ca²⁺-binding proteins known to associate with membranes and stress signaling; importantly, they can exhibit peroxidase activity and thereby reduce ROS levels during stress^[43]. Overexpression of an annexin gene in wild tomato (Solanum pennellii) was shown to improve drought tolerance by modulating ABA-mediated stomatal closure and enhancing ROS scavenging capacity (including lower lipid peroxidation and higher proline accumulation)^[44]. The induction of annexin in RRIT 251's roots thus likely supports membrane stability and redox homeostasis under dehydration.

Both rubber tree clones also showed drought-induced increases in Bet v 1/Major latex protein (MLP) family proteins (pathogenesis-related Class 10 [PR-10] proteins), with RRIT 251 displaying stronger induction of specific isoforms (~5–15-fold at 50% field capacity). Moreover, expression of the HbPR-10 gene showed opposite responses in RRIM 600 and RRIT 251. PR-10/MLP proteins are commonly induced by abiotic stress and are thought to enhance tolerance by binding protective ligands or modulating hormone signaling pathways^[45]. Supporting this role, a root-specific PR-10 protein from rice is rapidly induced by drought, and its overexpression significantly improved tolerance, partly by promoting root growth and development^[46]. The accumulation of PR-10/MLP proteins in drought-stressed rubber tree roots may therefore contribute to cellular stabilization and stress signal transduction. Their conserved functions across species suggest that these proteins form part of a generalized protective mechanism that strengthens plants' defense capacity under drought.

Drought stress also elicited adjustments in proteins related to water transport and redox homeostasis. Notably, aquaporins were strongly downregulated in roots under severe drought in both clones: A tonoplast intrinsic protein (TIP) and a plasma membrane intrinsic protein (PIP) showed drastic decreases in abundance (< 0.2-fold of well-watered levels at 30% field capacity). For the aquaporin gene expression study, the HbTIPs and HbPIPs showed an increasing trend increasing trend in RRIT 251 under drought, whereas they decreased in RRIM 600. Such downregulation of aquaporins is a known plant response to water deficit, effectively reducing membranes' water permeability to conserve tissue water^[47]. By closing aquaporin channels in their roots, plants reduce hydraulic conductivity and water loss, thereby helping to maintain cellular turgor during drought. In parallel, a putative 2-hydroxy-acid oxidase (likely a glycolate oxidase; Table 4) was uniquely upregulated (5.6-fold at 30% field capacity) in RRIT 251's roots. Glycolate oxidase catalyzes the oxidation of 2-hydroxy acids to glyoxylate, producing hydrogen peroxide (H₂O₂) as a byproduct^[5]. Increased H₂O₂ production in the roots could function as a signaling mechanism to trigger stress-response pathways, since H₂O₂ acts as a secondary messenger in plant stress signaling. However, excess H₂O₂ can be damaging, so the plant must balance its production and scavenging. The concomitant induction of GST and other ROS-scavenging proteins in RRIT 251 suggests that any ROS generated (for example, via enhanced glycolate oxidase activity) is promptly detoxified, allowing H₂O₂ to play its signaling roles without accumulating to toxic levels. This coordinated control of ROS reflects an intricate drought response where RRIT 251 actively manages oxidative stress while possibly using ROS as signals to strengthen defense responses.

To extend the gene expression analysis of Bet v I and aquaporin genes, both TIP and PIP genes' expression in response to drought stress was performed through quantitative real-time polymerase chain reaction (qRT-PCR). We found that HbPR-10 was clearly different between RRIM 600 and RRIT 251 leaves. In RRIM 600, the expression of the HbPR-10 gene was downregulated under moderate drought treatment, and the lowest expression level was present under severe drought treatment. The observed downregulation of aquaporins in both rubber tree clones presents a notable discrepancy with several other drought-tolerant crops, where upregulation is often reported to facilitate water transport under stress. In H. brasiliensis, this reduction in protein abundance likely reflects a specialized, species-specific water conservation strategy. By drastically lowering its membranes' water permeability, the tree may effectively isolate its internal hydraulic system from the drying soil, preventing back-flow and minimizing water loss^[47,48]. This conservative strategy appears to be aimed at preserving cellular turgor for prolonged survival, even at the cost of reduced overall water flux, which is a critical adaptation for woody perennials in drought-prone environments. On the contrary, in RRIT 251, the expression level of HbPR-10 was upregulated under both drought stress treatments. Under severe drought, RRIM 600 presented yellow leaves, which may be a sign of water deficiency. Therefore, the gene expression of RRIM 600 under drought conditions was reduced, which is consistent with the physiological parameter and proteomic results.

Integration of the proteomic and metabolomic data reveals a cohesive stress response strategy. A key metabolite change under severe drought was the massive accumulation of proline in the roots, especially in RRIT 251. Proline is a well-known compatible solute that functions as both an osmoprotectant and a ROS quencher in plants^[25,49]. High proline levels help maintain cell turgor by osmotic adjustment and protect cellular proteins and membranes by stabilizing their structure under drought stress. Proline can also directly scavenge hydroxyl radicals, thereby mitigating oxidative damage during drought^[49]. The observed proline surge in roots aligns with the proteomic findings: Downregulation of aquaporins would work in concert with osmolyte accumulation to retain water, and the enhanced antioxidant enzymes complement proline's role in scavenging ROS. Similarly, increases in sugars like raffinose under 30% field capacity suggest the activation of another osmoprotective mechanism. Raffinose family oligosaccharides accumulate in many plants during drought and function as stabilizing osmolytes that protect membranes and proteins from dehydration-induced injury^[50]. The accumulation of raffinose in drought-exposed rubber tree roots, together with proline, indicates a concerted effort to equilibrate the cellular osmotic potential and protect intracellular structures. This osmotic adjustment strategy is a hallmark of drought tolerance and is evident in RRIT 251's metabolism, likely contributing to its better maintenance of physiological functions (photosynthetic efficiency, as shown by the higher Fv/Fm in RRIT 251) under severe stress.

Metabolomic analysis revealed increased levels of signaling-related metabolites, including pipecolic acid and azelaic acid, in roots under severe drought. This pattern suggests an overlap between abiotic stress responses and defense pathways. Pipecolic acid is a lysine-derived metabolite known as a critical regulator of systemic acquired resistance in plants, and recent evidence indicates that the pipecolic acid biosynthetic pathway also contributes to abiotic stress tolerance by enhancing antioxidant defenses during drought^[12]. Pipecolic acid has been identified as a drought-enriched metabolite in the roots of drought-stressed sorghum (Sorghum bicolor), suggesting a conserved role in stress signaling across species^[51]. Azelaic acid is another mobile signal that primes plants for systemic immunity – it has been shown to induce or 'prime' salicylic acid accumulation and defense gene activation following the initial stress exposure^[13]. This interpretation is supported by the strong induction of PR-10/MLP proteins in RRIT 251, which are typically linked to pathogen defense but are also deployed during drought. Such crosstalk between drought and immune signaling has been reported in other species, where stress hormones and signaling molecules interact to fine-tune adaptive responses. The concurrent presence of pipecolic and azelaic acids with PR-10 induction suggests that RRIT 251 uses a multifaceted defense strategy, integrating osmolyte accumulation, antioxidant protection, and stress-responsive signaling. This integration of metabolic and proteomic regulation likely underpins the superior drought tolerance of RRIT 251 by enabling both immediate cellular protection and pre-emptive activation of defense pathways during prolonged water deficit.

Collectively, the integration of metabolomic, physiological, and proteomic findings in this study reveals that RRIT 251 exhibits a more coordinated drought adaptation strategy than RRIM 600 than RRIM 600, driven by coordinated metabolic and protein-level responses in the roots. RRIT 251 roots activated an antioxidant defense network, including GST and annexin, alongside high proline and raffinose accumulation, which together stabilize cellular structures and mitigate oxidative stress. This response is further complemented by the induction of systemic signaling molecules such as pipecolic acid and azelaic acid, suggesting a primed defense status beyond osmoprotection. The concurrent downregulation of aquaporins and upregulation of stress-related proteins, including PR-10/MLP proteins and (S)-2-hydroxy-acid oxidase, highlights a controlled water conservation mechanism coupled with redox signaling. In contrast, RRIM 600 primarily relied on structural protein changes and ABA-related signaling, with limited antioxidant reinforcement, which likely contributes to its lower physiological resilience under severe drought. The preservation of PSII's efficiency and delayed leaf senescence in RRIT 251 under drought further validate these molecular findings. This integrated multi-omics evidence positions RRIT 251 as a superior model for drought tolerance in rubber trees, offering molecular markers such as proline, raffinose, GST, and annexin as promising targets for future breeding strategies aimed at enhancing stress resilience in H. brasiliensis.

The integration of metabolomic and proteomic findings in this study provides critical insights for the natural rubber industry, where drought-induced reductions in latex yield remain a primary economic threat. Because latex production and regeneration are heavily dependent on maintaining cellular turgor and metabolic homeostasis, the root-centric adaptive strategies identified in RRIT 251 offer a clear pathway for developing more resilient plantations. Specifically, the massive accumulation of proline in the rootstock—orchestrated by the tolerant scion—serves as a vital mechanism for maintaining osmotic balance and protecting the tree's energetic capacity, as evidenced by its stable PSII efficiency under stress. Furthermore, the dramatic 7.85-fold upregulation of GST in RRIT 251 provides a robust proteomic defense against oxidative damage, which otherwise leads to premature leaf senescence and reduced latex flow.

These molecular signatures, specifically proline as a metabolic marker and GST as a proteomic marker, represent high-priority candidates for marker-assisted breeding programs. Historically, the development of rubber clones has been hindered by the long maturation period of over seven years before harvestability; however, utilizing these markers allows for the early-stage screening of scion–rootstock combinations at the seedling stage. By selecting for clones that exhibit the water-saving profile of aquaporin downregulation and strong antioxidant induction seen in RRIT 251, breeders can accelerate the deployment of climate-resilient cultivars. Ultimately, framing the drought response as a scion–rootstock interaction provides a new molecular framework for the industry to mitigate the long-term economic losses associated with climate-induced water deficits.

This study presents the first integrative metabolomic and proteomic analysis of drought adaptation mechanisms in H.brasiliensis, revealing RRIT 251 as a superior drought-tolerant clone compared with RRIM 600. RRIT 251 demonstrated root-specific metabolic reprogramming characterized by high proline and raffinose accumulation, coupled with the induction of antioxidative proteins such as glutathione S-transferase and annexin. These biochemical adjustments were supported by higher photosynthetic efficiency and delayed leaf senescence under water deficit. The coordinated activation of osmoprotectants, antioxidant enzymes, and signaling molecules in RRIT 251's roots underpins its enhanced drought resilience. This multi-omics evidence establishes RRIT 251 as a promising genetic resource for developing stress-resilient rubber cultivars, providing a molecular framework for marker-assisted breeding strategies to mitigate climate-induced drought challenges in rubber plantations.