Gossypium barbadense inbred lines Pima S-7 (PI 560140) and Pima 3-79 susceptible to Fusarium oxysporum f. sp. vasinfectum (FOV) race 4 (FOV4), and Pima-S6 resistant to FOV4, were used in this study. Pima entries are public lines released previously (Feaster and Turcotte, 1984; Turcotte et al., 1992; Kim et al., 2005). Pima-S6, specifically, was derived from the original source Pima S-6 after being identified as an important source of resistance traits to FOV4 (Ulloa et al., 2006; Ulloa et al., 2020; Ulloa et al., 2022a). The new source was named Pima-S6, as explained in (Chávez Montes et al., 2023).

The three Pima entries were used to prepare cDNA libraries and RNA-seq analyses. The FOV4 inoculum was produced in potato dextrose liquid medium supplemented with 3 mM streptomycin to prevent contamination and incubated at 28°C and constant shaking at 100 rpm during 3–4 days. We utilized plugs (2–3 mm²) of FOV4 fresh-grown onto potato dextrose agar (PDA, streptomycin), as the starting material. Finally, a conidial suspension (1 × 10⁶ conidia/mL) was generated by filtering the culture through cheesecloth.

The infection of the roots in the greenhouse was performed according to the root-cut dip method as previously described (Ulloa et al., 2013). Briefly, uprooted one-week-old cotton seedlings were dipped in the FOV4 conidia suspension and transplanted into sterilized U.C. mix #2 substrate after 3 min. Five treated plants were planted per pot and incubated in the greenhouse at 24°C ± 2°C. None inoculated plants treated with water only were used as controls. Four replicates for all treatments were used following a randomized complete block design. At 11 days after inoculation (dai), the plant and soil ball were removed from the pots, and the soil was removed gently from the roots by washing under continuous water flow. The washed roots were placed immediately in 50 mL screw-top plastic vials, liquid nitrogen was added, and the vials were placed in a −80°C freezer with the tops off to allow vaporization of the liquid nitrogen. Then the tops were screwed on, and the frozen roots were stored at −80°C until RNA extraction.

We used a similar RNA-seq data processing pipeline as that reported in (Ojeda-Rivera et al., 2022). Briefly, quality control of Illumina reads was performed using FastQC (version 0.11.9; https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Adapter sequences were removed from libraries with Trimmomatic (version 0.39) (Bolger et al., 2014). We then quantified gene expression using kallisto (version 0.46.1) (Bray et al., 2016). Pima cotton libraries were aligned to our assembly of the Pima-S6 G. barbadense genome (Chávez Montes et al., 2023). edgeR was used to perform differential expression analysis (Robinson et al., 2010), the tximport package (version 3.130) (Soneson et al., 2015) was used to import estimated counts from kallisto to edgeR. We applied a base 2 logarithm of fold change (logFC) threshold (−1>logFC>1) and a false discovery (FDR) threshold (FDR<0.05) to detect differentially expressed genes. We used a similar pipeline for GO analysis to that reported in (Ojeda-Rivera et al., 2022). Briefly, annotations for the G. barbadense Pima-S6 genome (Chávez Montes et al., 2023) were assigned using the GO annotations of the blastp reciprocal best hits from multiple databases as reference. Databases included Araport11 (Cheng et al., 2017), UniProt (Apweiler et al., 2004), Interproscan (v5.48.83) (Blum et al., 2021), Swiss-Prot proteins from multiple plant species (Wang et al., 2021), and GO annotations from the PANNZER2 (Törönen et al., 2018) functional annotation webserver with a PPV value of at least 0.5. We combined analyses into a non-redundant gtf file and used that file for GO enrichment analyses using the clusterProfiler package for R (Wu et al., 2021). The datasets generated and analyzed for this study can be found in the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database repository under accession number GSE233711.

For non-targeted metabolomics analysis, non-infected roots of the three cotton entries, Pima-S6, Pima 3-79, and Pima S-7 were utilized. Cotton seedlings were grown vertically in Petri dishes in a Conviron growth chamber set at 28°C, 150 µM/m²s light intensity, 16/8 h light/dark cycle, and 45% relative humidity. Culture media was 0.25 X MS with slight modifications (pH 6, 5 mM CaCl₂, and 1 mM MgCl₂) and 0.35% Phytagel (Sigma-Aldrich). Root samples for metabolomic analysis were harvested 8 days after germination. Five biological replicates were performed (n = 5).

Non-targeted metabolomic analysis and quantitative analysis of esculetin were performed in non-infected cotton root samples. Non-targeted metabolomics analysis was conducted using a reverse phase chromatography method and positive electrospray ionization mode mass spectrometry. Freeze-dried samples of cotton roots were powdered and extracted with a mixture of methanol-water (60:40 v/v) with 0.1% formic acid. The supernatant of this extraction was further concentrated in a vacuum centrifuge (Labconco, United States) and filtered before chromatographic analysis. Chromatographic separation was performed in a Vanquish Flex UHPLC system (Thermo Scientific, United States) using an Acclaim VANQUISH column (150 × 2.1 mm, with 2.2 µm) and a solvent gradient method starting 5% phase B for 4 min, then increasing to 100% B from 4–17 min, and finally maintaining 100% B for 4 min prior to re-equilibration time of 4 min. Phases consisted of Water (A) and Methanol (B) supplemented with 0.1% formic acid. Detection was carried out in an Orbitrap Exploris 240 (Thermo Scientific, United States) with electrospray ionization source in positive mode scanning from 70–800 Da with a resolution of 120,000 and Data Dependent Acquisition (DDA) of MS/MS with a cycle time of 3 s. Raw files were imported into Compound Discoverer 3.3 (Thermo Fisher Scientific) and processed with an untargeted metabolomics workflow. Chromatograms were aligned, compounds were detected and grouped with a 5 ppm tolerance window. Detected metabolites were identified using MZCloud MS/MS database with a similitude factor of at least 50%. Exact masses were also searched against BioCyc and KEGG databases (Wu et al., 2021), and the final lists were curated manually. Pathway enrichment analysis was done in the MetaboAnalyst online platform (Pang et al., 2022), using the exact m/z values of detected compounds for the MS Peak-to-Pathway enrichment (Li et al., 2013), retention times, and statistical values from Compound Discoverer analysis. The KEGG pathway analysis presented in Figure 4D was generated by graphing the results of the Metaboanalyst platform (Pang et al., 2022), which uses the KEGG database as a reference (Kanehisa et al., 2023), in the ComplexHeatmap R package (Bolger et al., 2014). Figure 5B, an illustration of the phenylpropanoid pathway, was generated from our gene-by-gene analysis of Arabidopsis orthologs in cotton.

Quantitative analysis of esculetin was done using an external standard methodology. High-purity esculetin (>98%) was acquired from Cayman Chemical (item:19286) and injected into the same UHPLC-MS system described before using a modified method. The separation column was a Waters BEH C18 (2.1 × 150 mm, 2.5 µm) with water and methanol added with 0.1% formic acid as mobile phases. The retention time for esculetin was 6.36 min, when methanol concentration in the chromatographic method was 50%. Detection in the Orbitrap MS was carried out in negative mode, using H-ESI. A full MS scan from m/z = 100–300 with a resolution of 120,000 preceded a targeted MS/MS using an intensity threshold filter of 5 × 103 and targeted mass for esculetin [M-H]- = 177.0193 ± 10 ppm. Collision energy for the MS/MS fragmentation was normalized to 30% and resolution to 15,000. A calibration curve from 10–10000 ng/mL was injected in duplicate (r.s.d. < 10%) with R2 = 0.9978 and used to quantify samples (n = 3) fragmentation pattern of the standard.

Total RNA was isolated from the root tissue of non-infected cotton plants, as described in previous sections. 100 ng of total RNA was used to carry out the reverse transcription reaction with SuperScript III reverse transcriptase (Invitrogen), following the manufacturer’s protocol. The RT-qPCR was performed in a MIC thermocycler (Magnetic Induction Cycler from Biomolecular Systems) using the reagent SensiFast TM SYBR No-ROX Kit from Bioline (BIO-11060). The PCR conditions were 95°C for 2 min, followed by 40 cycles at 95°C for 5 s, 65°C for 10 s, and 72°C for 20 s. The relative expression of each gene was calculated relative to the mean of the Cq of POLYUBIQUITIN 14 (UBQ14) gene for each sample. UBQ14 qPCR was made using previously reported primers (Artico et al., 2010). Primers were designed for PHENYLALANINE AMMONIA-LYASE 2 (PAL2) and PLEIOTROPIC DRUG RESISTANCE 12 (PDR12); the forward primer 5′-CAT​TGT​GAG​CCA​AGT​TGC​TAA​G-3′ and reverse primer 5′-GTC​CAC​GGC​TTT​GAG​TAA​GT-3′ for PAL2 gene and the forward primer 5′-GGA​GAA​TGG​TGA​GGC​ATT​TAG​A-3′ and reverse primer 5′-CAC​GAG​CAG​ACA​TGG​AGA​AA-3′ for PDR12) gene. To calculate the FC values, the Pfaffl formula, which corrects the FC accounting for differences in primer PCR efficiencies, was used (Pfaffl, 2001).

To characterize the transcriptional response to Fusarium wilt disease caused by FOV4 in Pima cotton, we performed RNA-seq analysis of plant roots grown under control and FOV4-infected conditions. Three Pima cotton cultivars were studied: one FOV4-resistant cultivar, Pima-S6, and two FOV4-susceptible cultivars, Pima S-7 and Pima 3-79 (see Materials and Methods). Plants of the three cotton entries were visually inspected after FOV4 infection throughout the duration of the experiment, looking for disease symptoms. Infected susceptible plants exhibited severe leaf wilt, chlorosis, and leaf abscission, leading to plant death, as reported previously (Ulloa et al., 2013). Resistant plants showed normal growth and minimal symptoms, and survived FOV4-infection. Representative photographs of infection experiments are presented in Figure 1. Based on our observations on disease severity, it was determined that 11 day was appropriate for tissue harvesting and conducting RNA-seq studies. Otherwise, advanced disease symptoms would prevent high-quality total RNA extraction. To this end, three biological replicates per treatment per cotton entry were used to generate 18 RNA-seq libraries. Alignment statistics and library names are listed in Supplementary Table S1.

To get insights into the molecular mechanisms underlying FOV4 resistance, we first performed a comparative analysis of the transcriptional response of all FOV4-susceptible (Pima S-7 and Pima 3-79) and resistant (Pima-S6) cultivars (Figure 2). We performed sample relationship analysis to determine whether the samples clustered by the tested genotypes and treatments. RNA-seq libraries were observed to cluster together in response to both genotype and treatment conditions (Supplementary Figure S1). Susceptible Pima 3-79 and Pima S-7 cultivars present similar expression profiles in response to FOV4 infection as they cluster closely in the sample relationship analyses, while Pima-S6 sample libraries cluster more closely and separated from those of Pima 3-79 and Pima S-7 in response to FOV4 treatments (Supplementary Figure S1). This result indicates a differential response to FOV4 infection in Pima-S6 compared to Pima S-7 and Pima 3-79.

To further analyze the differences in response to FOV4 infection among the tested cultivars, we performed a comparative transcriptomic response analysis to determine which genes are commonly regulated in response to FOV4 treatment (Figure 2A). We determined that 6,657 and 7,393 genes were upregulated and downregulated in Pima 3-79 in response to FOV4 treatment (Figure 2A) and that up to 63% and 70% of these gene sets are specific to this cultivar. In the case of Pima S-7, we detected the upregulation and downregulation of 3,977 and 5,028 genes, respectively, in response to FOV4 treatment (Figure 2A), of which 30% and 55% were specific to this cultivar. Lastly, in the case of Pima-S6, we detected that 1,586 and 2,820 genes were upregulated and downregulated, respectively, in response to FOV4 treatment (Figure 2A), of which 70% and 83% of these gene sets are specific to Pima-S6. Genes commonly upregulated and downregulated in susceptible cultivars Pima S-7 and Pima 3-79 are 2,318 and 1,928, respectively, which correspond to a ∼47% of the Pima S-7 transcriptional response to FOV4 (∼30% of the Pima 3-79 response), indicating that a relatively large number of genes of the response to FOV4 is commonly regulated in susceptible cultivars. In contrast, Pima-S6 and Pima S-7 had only 352 upregulated and 212 downregulated genes in common, corresponding to ∼6% and ∼12% of the FOV-4 responsive genes in Pima S-7 and Pima-S6, respectively. Considering the three genotypes, 105 genes were upregulated and 114 downregulated in response to FOV4; this number equals less than 4% of the global transcriptional response for all the genotypes. The transcriptomic landscape of the response to FOV4 (Figure 2B) agrees with the observed differences in gene expression between susceptible and resistant cultivars. It also indicates differences in the magnitude of the fold change in gene expression. Overall, our data showed that susceptible genotypes share a more significant portion of their transcriptional response to FOV4 among themselves than with Pima-S6, and that the transcriptional response of the resistant cultivar is specific and different from that of the susceptible cultivars.

To get biological insights into the FOV4-resistance mechanism(s) that are transcriptionally activated in Pima-S6 in response to FOV4 treatment, we performed differential Gene Ontology (GO) enrichment analysis (Figure 3). The complete list of GO-enriched categories, p-values, and gene identifiers is included in Supplementary Dataset S1. First, we determined the GO categories enriched in the gene sets that are transcriptionally activated in Pima 3-79, Pima-S6, and Pima S-7, respectively, in response to FOV4 treatment (Figure 3A). We then focused our analysis on categories significantly enriched only in Pima-S6, the FOV4-resistant genotype (Figure 3B). The identified GO categories presented in Figure 3B are related to a variety of biological processes including transmembrane transport (GO:0042626, GO:0005372), glycoside hydrolysis (GO:0016161, GO:0102229, GO:0000272, GO:0004565), metabolism of small molecules and amino acids (GO:0006520, GO:0004471), phospholipid metabolism (GO:0009395), phytohormone metabolism (GO:0010102, GO:0009850, GO:0009690, GO:0051213) and phenylpropanoid metabolism (GO:0046417, GO:0016841). Because the Arabidopsis genome is more functionally characterized than the cotton genome, we searched for Arabidopsis orthologs of the genes included in these GO categories to get more insights into their possible roles in FOV4 resistance (Figure 3C, Supplementary Dataset S1). GO enrichment and KEGG pathway analysis (Kanehisa et al., 2023) of Arabidopsis ortholog genes (Figure 3C) resulted in multiple categories clustered to similar biological processes as the categories from Figure 3B. We found clusters of GO categories that include the metabolism of small molecules, phenylpropanoid biosynthesis and metabolism, ABC-type transporter activity, lipid metabolism, regulation of enzymatic activity, diterpenoid biosynthetic processes, and phytohormone signaling-related processes.

The phenylpropanoid pathway in plants is the starting point for other metabolic pathways that produce critical secondary metabolites, including flavonoids, coumarins, and lignans, which have roles in plant resistance to pathogens (Fraser and Chapple, 2011). Gene-by-gene analysis of Arabidopsis orthologs (Supplementary Dataset S1) allowed the identification of genes that encode key enzymes involved in the phenylpropanoid pathway, including PHENYLALANINE AMMONIA LYASE 1 and 2 (PAL1,2; PimaS6.A09G041328.1, PimaS6.A11G018686.1, PimaS6.D11G014299.1), which catalyze the initial step of the phenylpropanoid pathway (Huang et al., 2010), and CHORISMATE MUTASE 2 (CM2; PimaS6.A06G072050.1, PimaS6.D06G069084.1) which catalyzes the first committed step in the synthesis of Phenylalanine (Fraser and Chapple, 2011). We also found an ortholog of CHALCONE SYNTHASE (CHS; PimaS6.D05G007786.1, PimaS6.A08G034631.1, PimaS6.A12G026213.1, PimaS6.D12G022648.1) that encodes the enzyme that catalyzes the initial step in the synthesis of flavonoids (Dao et al., 2011). Root exudation of secondary metabolism compounds through ABC-type transporters has also been shown to play a role in plant-microbe interactions and resistance to fungal pathogens (Chaparro et al., 2013; Fourcroy et al., 2014; He et al., 2019). Concomitantly, we found multiple cotton genes encoding ABC-type transporters orthologs to multiple Arabidopsis ABC-type transporters. Some of these transporters have been reported to play roles in fungal pathogen resistance, like PLEIOTROPIC DRUG RESISTANCE 12 (PDR12), by facilitating root exudation of secondary metabolites (He et al., 2019). The differential transcriptional response of Pima-S6 compared with susceptible genotypes Pima S-7 and Pima 3-79 indicates that secondary metabolism might play a key role in the FOV4 resistance mechanism. Interestingly, we found that genes in GO categories related to phenylpropanoid regulation and metabolism (GO:0009698, GO:0009715), have a higher basal expression in Pima-S6, the FOV4-resistant cotton, compared to the susceptible varieties, which is enhanced upon FOV4 infection (Supplementary Figure S2). This suggests that the enhanced basal expression of these genes might contribute to the FOV4 resistance phenotype.

To further explore the differences in secondary metabolism between the resistant and the susceptible genotypes, we performed mass spectrometry-based non-targeted metabolomics analysis in the roots of non-infected Pima-S6, Pima S-7, and Pima 3-79 seedlings (Figure 4). We detected 3,115 compounds in the roots of the three studied genotypes (Figure 4A). Cluster profiling of the normalized relative abundance of the identified compounds indicates that Pima-S6 and Pima S-7 metabolomes are more closely related to each other than to Pima 3-79 (Figure 4A). We performed pairwise comparisons to determine statistically significant differences in relative compound concentration of these 3115 metabolites between FOV4-resistant Pima-S6 and the susceptible genotypes, namely, Pima-S6 vs. Pima S-7 and Pima-S6 vs. Pima 3-79 (Figures 4B, C). We found 2,251 metabolites with higher levels and 18 with lower levels in Pima-S6 relative to Pima 3-79, indicating that Pima-S6 has a more diverse root metabolome than Pima 3-79 (Figure 4B). The Pima-S6 vs. Pima S-7 comparison showed that 204 metabolites have higher and 190 lower levels in Pima-S6 (Figure 4C), which agrees with the clustering profiling analysis (Figure 4A).

To get functional insights into the metabolomic differences between Pima-S6 and FOV4-susceptible genotypes (Pima-S6 vs. Pima S-7 and Pima-S6 vs. Pima 3-79), we performed a functional analysis using MS Peak-to-Pathway analysis (Li et al., 2013; Pang et al., 2022), which harnesses metabolite changes at pathway level to calculate enrichment and combines it with metabolite set enrichment analysis to determine an overall pathway enrichment p-value (Figure 4D; Supplementary Dataset S1). In the case of the pairwise comparison between Pima-S6 and Pima 3-79, it was determined that the “cutin, suberin, and wax biosynthesis” pathway was significantly more enriched in Pima-S6 than in Pima 3-79 (Figure 4D). Along with lignin, cutin, suberin, and wax are secondary cell wall components that, amidst other functions, provide protection against pathogen attack and interact metabolically with the phenylpropanoid synthesis pathway (Miedes et al., 2014; Andersen et al., 2021). In the case of the pairwise comparison between Pima-S6 and Pima S-7, six pathways were found significantly more enriched in Pima-S6 than in Pima S-7, including “flavonoid biosynthesis”, “lysine biosynthesis”, “ubiquinone and another terpenoid-quinone biosynthesis”, “flavone and flavonol biosynthesis” and “phenylpropanoid biosynthesis” (Figure 4D). All these categories, except for lysine biosynthesis, are closely related to secondary metabolism, specifically the phenylpropanoid pathway.

We rationalized that the set of metabolites key to FOV4-resistance in Pima-S6 must have enhanced levels compared to susceptible genotypes. To test this hypothesis, we performed Venn diagram analysis to determine more abundant metabolites in Pima-S6 compared to FOV4-susceptible Pima 3-79 and Pima S-7 genotypes (Figure 4E). We identified 175 compounds more abundant in Pima-S6 than in the two susceptible genotypes (Supplementary Dataset S1, spreadsheet 9.1). Their identity was determined by matching their spectra to known databases using the Compound Discoverer software (Thermo Scientific). We found 57 hits to known molecules (∼32%) among the Pima-S6-enhanced metabolites (Supplementary Dataset S1, spreadsheet 9.2). Inspection of the list revealed multiple molecules related to secondary metabolism, specifically phenylpropanoid metabolism, as suggested by the presence of a phenyl group in their structure (Supplementary Dataset S1), which agrees with the pathways presented in Figure 4D. Interestingly, we found enhanced levels of the coumarin esculetin in the roots of Pima-S6 (Figure 4F). Root coumarin exudation has been recently shown to inhibit Fusarium oxysporum f. sp. raphani growth in vitro and promote the establishment of plant-beneficial microorganism communities in the rhizosphere (Stringlis et al., 2018). A quantitative determination was performed using an external high-purity standard to corroborate the metabolite esculetin’s identity. UHPLC-MS analysis suggests that the fragmentation pattern of the esculetin standard and our samples matches and is one of the most abundant in Pima-S6 (Figure 4F; Supplementary Figure S3).