An aeroponic growth chamber trial was performed in 2021. The pots were arranged in a completely randomized design with two treatments, progressive drought stress (drought) and well-watered control (control). Each treatment consisted of 10 pots and each pot had three plants. Treatments were applied for 11 days (water-deficit phase) once plants reached the match-head square stage (Elsner et al., 1979). We specifically targeted this growth stage for cotton because it is particularly vulnerable to drought, which can result in a severe reduction in yield and fiber quality (Zonta et al., 2017). The control treatment received full irrigation (130 mL hr^-1)of full-strength Hoagland solution (Hoagland and Arnon, 1950) throughout the entire experiment, while the drought treatment received a 50% reduction in irrigation every 2 to 3 days that ended with 10% full irrigation (water-deficit phase). After 11 days of treatment, irrigation was returned to 100% for all pots for 7 days (recovery phase). Root exudates were non-destructively collected at 0, 2, 4, 7, 9, 11, 14, 16, and 18 days after the experiment was initiated, which includes a baseline (day 0), water-deficit phase (day 2, 4, 7, 9, and 11), and recovery phase (day 14, 16, and 18). Five out of ten pots were randomly selected for sample collection on each sampling day. A total of 90 samples were collected, which consisted of 2 treatments and 5 replicates for 9 sampling days (Figure 1).

A moderately drought susceptible upland cotton (Gossypium hirsutum, cv. TAM 421) was chosen for this study. Cotton seeds were surface sterilized with 1 min of 70% EtOH, 3 min of 2% NaOCl, and rinsed with sterile distilled water 3 times. After the seeds germinated on the germination papers, three seedlings were placed in individual mesh net cups and supported with hydroponic sponges in each pot (~3.8 L). The hydroponic system was set up with an air pump equipped with two 10-way splitters. Each valve on splitter was connected to an air stone with airline tubing for each pot. Nutrients were supplied with a full-strength Hoagland solution (Hoagland and Arnon, 1950), where the nutrient solution was replaced twice a week, maintaining a pH range of 6.0-6.5. Plants were grown in a hydroponic system with constant aeration with dissolved oxygen levels ranging between 7.7 and 8.0 mg L^-1 until 3 to 4 true leaves had developed. Then, the plants were transferred to an aeroponic system with 3 plants per pot. The set up of aeroponic system and the misting cycle was followed by Lin et al. (2022, preprint). In brief, the nutrient solution was pressurized using a diaphragm pump with a built-in pressure shut-off switch set to 827.3 kpa (120 psi). The pressurized solution was then temporarily stored in a 0.75 L pressure accumulator tank before it was delivered to each pot. The misting cycle was controlled by digital timers. For each pot (~7.6 L) in the aeroponic system, a spray ring was installed at the bottom (Figure 1B), composed of 3 nozzles (0.1 mm diameter). Misting was set for 10 sec, and time between misting events ranged between 5 min (well-watered) and 50 min (extreme drought). Environmental conditions in the walk-in growth chamber (EGC, USA) were set as 12-hour light/12-hour dark photoperiod. The air temperature and relative humidity settings were based on the local weather station to mimic the field environment (Supplementary Table 1). A complete description of plant growth conditions, aeroponic system design, and plant performance were reported in Lin et al. (2022, preprint). In brief, the drought symptoms (less than 5% of leaf affected) were first visually observed on day 4 of the drought treatment. The symptoms reached the greatest severity (over 66% of leaf affected) on day 9 and were approaching the wilting point on day 11. The plants in the drought treatment were visually fully recovered by the end of recovery phase.

Root exudates were collected from the whole root system in the aeroponic system at 14:00-16:00 on each collection day (Lin et al., 2022, preprint). The choice of the sampling time was based on our preliminary studies, which revealed a trend of higher concentrations of a targeted compound, abscisic acid, during the afternoon collections. In brief, roots were rinsed with a sampling solution (0.05 mM CaCl₂, pH 6.0-6.5) for 2 min using the aeroponic system with a clean mist collection container. The choice of a dilute CaCl₂ solution was designed to reduce the drastic change in ionic strength and mimic the soil solution, which is typically dominated by Ca²⁺ ions. The collected exudates were then filtered through a 22-μm polyethersulfone syringe filter and aliquoted into 1.5 mL centrifuge tubes. Samples were then stored at -80°C until further analysis.

The samples were desalted and concentrated by solid phase extraction (SPE) using the method by Dittmar et al. (2008). The SPE cartridge (Bond Elut PPL, 100 mg, Agilent, CA, USA) was rinsed with 3 mL of methanol (one cartridge volume) to activate the column. A 1 mL aliquot of root exudates was first diluted with 14 mL of sterile double-distilled water and acidified with 1 M HCl until the pH was between 2 to 3. The samples were then vacuumed through the SPE cartridge under 170 mbar. Prior to elution, a total of 30 mL of 0.01 M HCl was rinsed through the cartridge to remove possible salt molecules and interferents. After allowing the sorbents to air dry, samples were eluted with 1.5 mL of methanol. Extracted samples were stored at -20°C until analysis.

Extracted samples were analyzed with a 9.4 Tesla Bruker FT-ICR spectrometer located at the University of Arizona. A standard Bruker electrospray ionization (ESI) source was used to generate negatively charged molecular ions, where samples were introduced via direct infusion to the ESI source. 144 scans were averaged for each sample and internally calibrated using an organic matter homologous series separated by 14 Da (CH2 groups). Data Analysis software (BrukerDaltonik version 4.2) was used to convert raw spectra to a list of m/z values applying FTMS peak picker module with a signal-to-noise ratio (S/N) threshold set to 7 and absolute intensity threshold to the default value of 100. Putative chemical formulae were then assigned using Formularity (Tolić et al., 2017) software as previously described in Tfaily et al. (2018). Chemical formulae were assigned based on the following criteria: S/N > 7 and mass measurement error < 1 ppm, taking into consideration the presence of C, H, O, N, S and P and excluding other elements. The data produced by FT-ICR MS (peak masses, peak intensities, and metabolic molecular formula) were then processed through Metabodirect pipeline (Ayala-Ortiz et al., 2023). Peak masses (m/z) < 200 and > 900 were filtered using -m option to target secondary metabolites (Hsu et al., 2011; Park et al., 2013). Data was z-score normalized using –norm_method option. The quality control steps, including 13C isotope filtering and error filtering (0.5 ppm), followed the default setting. Masses must be identified in at least two samples to be included in the analysis. Elemental types (CHO, CHON, CHONP, CHONS, CHONSP, CHOP, CHOS, and CHOSP) were assigned for the filtered masses. KEGG analysis and mass difference analysis were performed using -k and -t options, respectively. Mass differences were used to infer potential biochemical transformations. Potential functional characteristics of the assigned molecular formulas were obtained by mapping to the KEGG database. Despite the molecular formula being the only criterion used for mapping to KEGG, it is important to acknowledge that these assignments are tentative. As noted, there is the possibility of multiple KEGG metabolites, particularly isomers, having the same molecular formula. Therefore, the molecular formula assignments should be interpreted with caution. The drought-induced unique metabolites in water-deficit phase (day 2 to 11) were characterized by comparing with control-treated samples excluding the drought-unique peaks identified at baseline (day 0). The KEGG pathways that were associated with human metabolisms were manually removed from the annotation list to generate the count table. The number of drought-unique metabolites were further visualized with UpSet plot (Lex et al., 2014) to identify the intersections between days (day 2 to 11). The figures were generated with R 4.2.0 (R core team, 2021) or GraphPad Prism v.6.0 (GraphPad Software, La Jolla, CA).

The abundance matrix of elemental compositions was used to compute Euclidean distance matrices using vegdist function in R package vegan (Oksanen et al., 2020). The data were subset by experimental phases (i.e., baseline, water-deficit, and recovery) to generate each distance matrix. Permutational multivariate analysis of variance (PERMANOVA) analysis was performed for each experimental phase to assess the experiment effect using the adonis2 function in R package vegan (Oksanen et al., 2020). In a PERMANOVA model, the response variable was the distance matrix, and the explanatory variables were treatment (drought and control), days of treatment, and the interaction of treatment × days. Principal component analysis (PCA) was performed to detect the patterns in the elemental compositions using prcomp function in R package stats (R core team, 2021).

A total number of 33,870 m/z were identified by FT-ICR MS. After the data was filtered using quality control steps (Ayala-Ortiz et al., 2023), a total of 13,033 m/z remained. After formula assignment, there were 3,985 metabolites with an assigned molecular formula. Overall, the total number of metabolites increased by 78% as drought stress progressed from day 2 to 11 (Figure 2). The highest number of total metabolites and drought-unique metabolites were identified under severe drought (day 9), the time point prior to wilting point (day 11) (Figure 2). Up to 70% of metabolites were shared between drought and control treated plant during water-deficit phase (Figure 2D). During the water-deficit phase, there were two discrete sampling days when the total number of metabolites increased (Figures 2A, C). The first increase occurred on day 4 of the drought treatment when the initial drought symptoms visually appeared. The second increase was observed on day 9, when the plants reached severe drought. Under drought, there was a reduction of the total number of metabolites at day 11 (i.e., near permanent wilting point) and a decrease throughout the recovery phase (day 14-18) (Figures 2A, B). The described trend in number of metabolites corresponded to the drought-unique metabolites. When grouping the data by phase, a greater total number of metabolites and drought-unique metabolites were observed in the water-deficit phase compared to day 0 and recovery phase (Figures 2B, D). Unique metabolites at baseline were identified between drought and control samples and those metabolites were excluded in some of the following analysis (Figure 2).

PERMANOVA showed no differences between treatment, days, and treatment × days interaction in elemental composition (i.e., CHO, CHON, CHONP, CHONS, CHONSP, CHOP, CHOS, and CHOSP) when all compounds were used. Treatment had a nearly significant effect on elemental composition (p-value=0.07) during the water-deficit phase (Table 1). PCA analysis of elemental composition showed a minor effect of treatment during water-deficit phase (Figure 3) but not at baseline or during the recovery phase (Supplementary Figure 1). Although high percentages of explained variance in PC1 (72.4%) has been observed, the treatment effect on elemental compositions during water-deficit phase was observed at PC2 with 11.8% explained variance and the variation was mostly driven by CHO and CHON (Figure 3).

The unique peaks identified from the drought-treated samples were further visualized using an UpSet plot (Figure 4). There were 345 metabolites consistently identified throughout the water-deficit phase. The number of distinct metabolites identified at individual days ranged from 118 to 481 with the highest number identified at day 9 in the drought treatment (Figure 4). While looking at the different intersections across the water-deficit phase, there were > 500 unique metabolites identified under severe drought stress (day 9-11) and < 100 metabolites uniquely identified at early drought stress (day 2-4) (Figure 4).

KEGG annotation analysis identified a total of 68 pathways and over 17 modules that were induced under water-deficit phase (Tables 2, 3; Supplementary Tables 2, 3). In the drought treatment, 30 out of 68 pathways were consistently observed from day 2 to 11 including flavonoid biosynthesis, carotenoid biosynthesis, phenylpropanoid biosynthesis, and plant hormones biosynthesis (Table 2). A trend of increased counts was observed among those consistently induced pathways and modules (Tables 2, 3). The highest counts for drought-induced unique metabolites across most pathways were primarily observed on day 9 (severe drought) followed by a reduction at day 11 (near wilting point) (Table 2). Several pathways that might be associated with microorganisms, such as phosphotransferase system (PTS), biofilm formation and quorum sensing, and biosynthesis of several antibiotics were also identified. It is important to note that many assigned molecular formulas were mapped to higher KEGG hierarchies “metabolic pathway” and “biosynthesis of secondary metabolites” without further classification (Table 2).

A total of 99 potential biochemical transformations were identified in this study (Table 4; Supplementary Table 4). Over 97% of the identified transformations had the highest number of counts in the drought-treated exudates at day 9, including methylation (-H), oxidation/hydroxylation (-H), hydrogenation/dehydrogenation, ethyl addition (-H₂O), etc. (Table 4; Supplementary Table 4). The increase of transformation along the progressive drought stress showed a similar trend that was identified in the mass distribution and KEGG annotation results.

During water-deficit phase, > 100 assigned molecular formula in drought-induced root exudates were mapped to pathways associated with the biosynthesis of flavonoids, including isoflavonoid, flavone, and flavonol (Table 2; Supplementary Table 2). Flavonoids are a group of plant-derived secondary metabolites that consist of > 10,000 compounds (Winkel-shirley, 2001; Sugiyama and Yazaki, 2014). The biosynthesis of flavonoids starts with phenylalanine and malonyl-CoA as the direct precursors. Then it undergoes different chemical reactions, such as hydroxylation, acylation, methylation, malonylation, and prenylation, which results in nine major subgroups (e.g., isoflavonoids, flavones, flavonols, flavandiols, and condensed tannins) having diverse structures and functions (Winkel-shirley, 2001; Sugiyama and Yazaki, 2014). Therefore, it was not surprising that we also observed the pathways associated with the biosynthesis of phenylalanine, tyrosine, and tryptophan; biosynthesis of phenylpropanoids; and KEGG modules involved with biosynthesis of malonyl-CoA under drought conditions (Table 2; Supplementary Table 3). In addition, a large number of methylation and malonylation transformations were also identified in drought-treated samples (Table 4; Supplementary Table 4).

Under drought conditions, a greater concentration of flavonoids have been identified in root exudates in holm oak (Gargallo-Garriga et al., 2018) and pearl millet (Pennisetum glaucum; Ghatak et al., 2022). Flavonoids have been reported as an important player in drought tolerance by directly functioning as antioxidants or serving as signaling molecules for the symbiosis of plants and microbes (Hassan and Mathesius, 2012). In our study, several masses were mapped to the KEGG module related to “daidzein → medicarpin” reactions (Table 3; Supplementary Table 3). Daidzein is an isoflavone compound reported to be secreted by soybean roots, as signaling molecules mediating communication between plants and nitrogen-fixing bacteria (Sugiyama et al., 2007; Okutani et al., 2020). It is unclear if cotton also uses similar strategies to attract nitrogen-fixing bacteria. However, several nitrogen-uptake-related bacteria, such as Mesorhizobium, Sinorhizobium, and Rhizobium, have been identified in the cotton rhizosphere (Qiao et al., 2017), which provide potential directions for future studies. Further investigations are required to characterize and quantify the flavonoid compounds in cotton root exudates and the associated microbes.

It has been well documented that multiple ubiquitous plant hormones, including abscisic acid (ABA) and jasmonic acid (JA), are coordinated in response to drought stress (Wang et al., 2021; El Sabagh et al., 2022). In this study, drought treatment induced the pathways and modules involved in the biosynthesis of ABA and JA, as well as their associated precursors, carotenoids (Wasilewska et al., 2008) and alpha-linolenic acid (Wasternack and Hause, 2013) (Tables 2, 3). In a previous study, we were able to detect ABA in collected root exudates using GC-MS, although the differences between drought and control treatment were only detected at certain days during water deficit phase (Lin et al., 2022, preprint). We also identified the pathway associated with ATP-binding cassette (ABC) transporters induced under drought, which have been reported to transport phytohormones and root exudation (Badri and Vivanco, 2009; Kang et al., 2011). Because our root exudation collection method was non-destructive compared to other sampling approaches (Oburger and Jones, 2018), we assume that cotton undergoes active transport of those plant hormones into the rhizosphere in order to regulate hormone concentrations inside the root tissues in response to drought, as previously reported in upland rice (Oryzae sativa) (Shi et al., 2014). However, it would be necessary to identify the associated transporters to support this assumption and further understand the mechanisms behind this finding.

Interestingly, pathways associated with the biosynthesis of antibiotics (e.g., novobiocin, enediyne antibiotics, and phenazine) were identified under drought, even in soil-less environments. Those antibiotics have been reported to be produced by a variety of saprophytic or endophytic microorganisms (Steffensky et al., 1998; Laursen and Nielsen, 2004; Chen et al., 2010; Li et al., 2015). Streptomyces, which includes both endophytic and saprophytic species, has been known to produce novobiocin and enediyne antibiotics (Steffensky et al., 1998; Chen et al., 2010). Phenazine is also a common antibiotic that can be isolated from Streptomyces and Pseudomonas (Laursen and Nielsen, 2004; Li et al., 2015) that is known to be involved in biofilm formation, which is important process in drought tolerance in plant roots (Mahmoudi et al., 2019). Our growth chamber experiment was not performed in a sterile environment, so by using a soilless aeroponic setup, is it likely that those antibiotics might be produced by host endophytes. Indeed, studies have shown that Streptomyces was enriched in the root under drought and expected to be associated with improving plant fitness (Xu and Coleman-Derr, 2019). Pseudomonas also has been reported as one of the dominant endophytic organisms in cotton roots (Shi et al., 2020). Our current experiment design is not able to differentiate the metabolites that were released by plants or endophytes; however, incorporating isotope labeling techniques could further reveal the mechanisms between plants and microbes under drought.

A major shift in the metabolomic profile occurred in our study under severe drought prior to the plant reaching wilting point. Our findings provide temporal responses to progressive drought with activation of various metabolic pathways, which provide important information on relevant sampling time points to further characterize and validate the drought-response mechanisms in cotton root exudates. The transition from severe drought stress to wilting point (day 9 to 11) may be related to cell death, giving insight into the complicated senescence of plants and how root exudates respond to plant death, which has limited attention in recent studies (Prudence et al., 2021). One pathway, ferroptosis, associated with cell death, was identified exclusively on day 9 (Table 2). Ferroptosis is an iron-dependent and highly regulated cell death process (Riyazuddin and Gupta, 2021; Distéfano et al., 2022). In Arabidopsis thaliana, this programmed cell death process was induced by heat stress and was mainly dependent on reactive oxygen species (ROS)-mediated lipid peroxidation (Distéfano et al., 2017). Currently, there are few studies on this newly discovered pathway and no reports indicating whether drought stress will induce this pathway. However, many abiotic stresses, including drought, also induce the accumulation of ROS in plant tissues (Lee and Park, 2012). Shared exudation patterns among heat, drought, and the combination of both stresses have been reported in maize (Tiziani et al., 2022), indicating there are shared responses when plants encounter these types of stresses. The small amount (only 17) of assigned molecular formula that were mapped to the ferroptosis cell death pathway, which were only identified during severe drought (day 9) in our study, suggests that the plant might release “final warning stress signals” prior to wilting point, but it also suggests there might be minor confounding responses between drought and cell death at this collection timepoint.