IntroductionCassava (Manihot esculenta) is an important crop grown for its starchy roots and is a staple food for millions of people in tropical regions. Cassava can not only be used as food but also can serve as a raw material for bioethanol production. However, cassava production can be affected by several pathogens, including Xanthomonas axonopodis pv. manihotis (Xam), the pathogen causing cassava bacterial blight (CBB). Xam is a gram-negative bacteria that infects plants via stomata and wounds1 and CBB is one of the most devastating diseases affecting cassava production, leading to significant yield losses2. This disease is characterized by leaf wilting, chlorosis, and necrosis, resulting in plant death. Inside plant tissues, colonization of Xanthomonas spp. occurs in the intercellular space or xylem vessels before spreading systemically within the plant3. The dispersal of Xam is stimulated by wind-driven rain splashes4causing widespread distribution in cassava fields. Moreover, Xam has an epiphytic stage that causes persistence of the pathogen between cropping cycles5. Host-pathogen interactions exist between Xanthomonas spp. and host plants, as bacterial effectors can induce the expression of some host genes6. These effectors can specifically bind to and induce a gene in a host in a so-called gene-for-gene interaction. The induction of bacterial effectors can affect the degree of resistance in various cultivars of the host plant species7. At present, the differential expression profiles of microRNA (miRNAs), small RNA (sRNA) and genes in cassava in response to Xam exposure have been reported8,9,10,11 but protein expression profiles remain largely unexplored, particularly with respect to the evaluation of proteome wide changes. Proteomic studies have become essential for investigating complex biological processes involved in plant-pathogen interactions. In the context of CBB, proteomic analysis can help elucidate the changes in protein expression that occur in response to Xam infection, providing insights into the molecular mechanisms underlying the plant defense response. Additionally, intracellular oxidative stress occurs after plants are infected with pathogens, producing reactive oxygen species (ROS)12. Polyunsaturated fatty acids (PUFAs) can be oxidized by ROS to generate malondialdehyde (MDA), an indicator of oxidative stress13. This study compared the proteome dynamics of two cassava cultivars, the tolerant industrial cultivar Rayong 72 (R72) and the susceptible edible cultivar Hanatee (HN), in response to Xam infection. Specifically, this study investigated changes in protein expression related to MDA levels in infected plants due to the involvement of ROS in pathogen defense response. Through this comparison, the study aimed to identify key proteins and pathways involved in cassava’s resistance to CBB, with the goal of developing more effective disease management strategies and enhancing cassava production.Materials and methodsPlant materialsCassava cultivars approved by the Department of Agriculture, Thailand including R72, a CBB-tolerant cultivar, and HN, a CBB-susceptible cultivar were used in this study. The R72 and HN stems were provided by the Nakhon Ratchasima Agricultural Research and Development Center, Department of Agriculture, Thailand. The stems were cut into lengths of 10–15 cm and planted in a 6 × 4-inch soil bag in the greenhouse with temperatures ranging from 36 to 24 °C.Inoculum preparationXam isolate C (XamC) was cultured on YPGA solid medium containing 0.5% yeast extract, 0.5% peptone, 0.5% glucose, and 1.5% agar for single-colony selection at 28 °C for two days. The selected single colony was grown further on YPGA medium again for two days at 28 °C. These colonies on the plate were resuspended in YPG liquid medium and calibrated to OD600 nm = 0.2 or 108 colony forming units (CFU) per mL using a spectrophotometer. Stem puncture inoculation was performed on eight-week-old cassava plants. First, an incision point was made on the stem between the third (L3) and fourth (L4) fully expanded leaves using a sterile needle. Five microliters of calibrated inoculum of XamC were placed at the puncture point. Leaves L3, L4, and stems with incision points were collected at 0, 3, 7, 14, and 28 days after inoculation (DAI), and five plants per DAI were collected. Sample leaves and stems were immediately placed in liquid nitrogen and stored at −80 °C until use.Symptom evaluationThe cassava response to XamC infection was determined by evaluating the area under the disease progression curve (AUDPC). Tolerant and susceptible cassava cultivars were compared at 0, 3, 7, 14, and 28 DAI on a severity scale of 0–5 as established by Lozano and Laberry, 198214. Infected plants were scored as follows: 0 as healthy plants, 1 as a dark area or necrosis around the inoculation point, 2 as gum exudates on the stem, 3 as wilting of one or two leaves and exudates, 4 as > 2 wilted leaves, and 5 as complete wilting and dieback. Five biological replicates were analyzed for each cultivar.Shotgun proteomic analysisThe fourth leaves (L4) from infected R72 and HN cassava cultivars were collected, immediately snap-frozen in liquid nitrogen, and stored at − 80 °C until used. Leaf samples were subjected to protein extraction by the trichloroacetic acid (TCA)/acetone precipitation method15. Briefly, about 100 mg of L4 was ground to a fine powder in liquid nitrogen and homogenized with extraction buffer containing 1% sodium dodecyl sulfate, 0.1 M tris hydrochloride (pH 6.8), 2 mM EDTA, 20 mM dithiothreitol (DTT), and 1x protease inhibitor cocktail (Merck, Rahway, NJ). The mixture was centrifuged at 12,000 g for 5 min. The supernatant was transferred to a fresh tube, and the proteins were precipitated with a 1:1 ratio of 20% TCA in acetone and samples were incubated on ice for 30 min. The protein pellets were washed with acetone to obtain a white color and finally dissolved in a buffer containing 7 M urea, 2 M thiourea, and 20 mM DTT. The protein concentration was measured by the Bradford assay using a Bio-Rad protein assay kit (Bio-Rad Laboratories, Hercules, CA) which included bovine serum albumin (BSA) as a protein standard. Protein samples (5 µg total) for liquid chromatography with tandem mass spectrometry (LC-MS/MS) were prepared by incubation with 10 mM DTT in 10 mM NH4HCO3 at 56 °C for 1 h to reduce disulfide bonds. The protein samples were added to 10 mM iodoacetamide in 10 mM NH4HCO3 for alkylation in the dark at 28 °C for 45 min. Proteins were dehydrated using 100% acetonitrile before protease digestion. Protein samples were digested by sequencing grade trypsin at 37 °C overnight. Trypsin-digested peptides were extracted from the samples by incubating with 50% acetonitrile in 0.1% formic acid at 28 °C for 10 min, then dried at 45 °C for 4 h, and kept at −80 °C until the LC-MS/MS analysis. LC-MS/MS was performed using an Ultimate 3000/Capillary LC system (Thermo Fisher Scientific, Waltham, MA) coupled to a quadrupole Q-TOF Impact II (Bruker Daltonics Ltd., Hamburg, Germany) equipped with a Nano-captive spray ion source. The peptide enrichments were conducted on a p-Precolumn 300 pm I.D. × 5 mm C18 Pepmap 100, 5 pm, 100 A and separated on a 75 μm I.D. × 15 cm. Samples were then packed with Acclaim PepMap RSLC C18, 2 μm, 100 A, nanoViper (Thermo Fisher Scientific, Waltham, MA). Subsequently, solvents A (0.1% formic acid in water) and B (0.1% formic acid in 80% acetonitrile) were supplied to the analytical column. A gradient was maintained using a 5–55% solvent B linear concentration to elute the peptides at a constant flow rate of 0.30 µL/min for 30 min. Electrospray ionization was performed at 1.6 kV using a CaptiveSpray. Mass spectra (MS) and MS/MS spectra were obtained in positive-ion mode over the range of 150–2200 m/z using the Compass 1.9 software (Bruker, Billerica, MA) with BSA as a standard.Protein identification and data analysisMS/MS spectra were identified using MaxQuant software version 1.6.6.016. Protein identification was performed by searching the Manihot esculenta reference sequence database of UNIPROT (Manihot esculenta, UP000091857). Mass spectra data in this study were deposited to ProteomeXchange repository: PXD062377 (https://proteomecentral.proteomexchange.org/cgi/GetDataset?ID=PXD062377) and jPOST repository: JPST003734 (https://repository.jpostdb.org/entry/JPST003734)17. Proteins with a q-value greater than 0.95 were selected for further analysis. Missing intensity values were imputed using the MissForest package in R, which employs a random forest-based approach. Protein intensities for individual samples were normalized against the total intensity of each respective sample and expressed as percentage expression ratios.Differentially abundant proteins between the tolerant (R72) and susceptible (HN) cultivars were calculated as the log2 fold change of R72 relative to HN. Statistical significance was assessed using the false discovery rate (FDR), and the results were visualized as volcano plots generated with the ggplot2 package19 in R. K-means clustering (KMC) was performed to explore expression trends of significantly differentially abundant proteins in R72 over time. Clustering was supervised by specifying 7 clusters and 100 iterations, with the optimal number of clusters determined using Gap-statistic analysis20,21 via the factoextra package21 in R. Correlations between the mean expression of each cluster and observed symptoms were analyzed using the corrplot package in R22.The subcellular localization of proteins was predicted by submitting protein sequences to the Plant-mSubP website, which identifies target localization within plant cells23. Protein sequences were analyzed using the PseAAC-NCC-DIPEP feature set, which incorporates pseudo-amino acid composition, N-terminal and C-terminal amino acid composition, and dipeptide composition.Functional annotation and enrichment were calculated using the Mercator4 website24. The identified proteins were subjected to functional prediction and enrichment analysis in cassava using the MapMan BIN ontology with 28 main biological categories.Determination of leaf malondialdehydeMalondialdehyde (MDA) content was determined using 2-thiobarbituric acid (TBA) following the method described by Sarker et al., 201825. Fresh weights of ground leaves (L3 and L4) and stems (including the inoculation point) were recorded prior to homogenization with 800 µL of 0.6% TBA in 10% trichloroacetic acid (TCA). The mixture was heated at 95 °C for 15 min and then cooled on ice. After cooling, it was centrifuged at 12,000 × g for 5 min, and the absorbance of the supernatant was measured at 450, 532, and 600 nm. MDA content was calculated based on fresh weight (FW) and expressed as nmol/g FW using the following equation: MDA = 6.45 × (A₅₃₂ − A₆₀₀) − 0.56 × A₄₅₀.Quantitative reverse transcription polymerase chain reaction (qRT-PCR) analysisTo validate the candidate proteins involved in resistance to CBB, qRT-PCR was performed to assess the expression levels of the corresponding RNA. Fourth leaves from infected two cassava cultivars were immediately snap frozen in liquid nitrogen and stored at −80 °C until used. Total RNA was extracted from the fourth leaves of R72 and HN, with 5 samples collected from each time point (0, 3, 7, 14, and 28 DAI) using Fruit-mate™ (Takara, Kusatsu, Japan) and TriReagent (Molecular Research Center, Cincinnati, OH). Contaminating DNA was removed by treatment with RQ1 RNase-Free DNase (Promega, USA). Then cDNA was synthesized using a Viva cDNA synthesis kit (Vivantis, Shah Alam, Malaysia) primed with oligo (dT) as recommended by the manufacturer. For gene expression analysis, cDNA templates were amplified using KAPA SYBR® FAST qPCR kit (Thermo Fisher Scientific, Waltham, MA) by Mastercycle ep Realplex (Eppendorf) with selected primers (Table S2) that showed efficiency, slope, and R2 values according to acceptable ranges as the values of slope between − 3.6 and − 3.1, %efficiency between 90 and 105, and correlation coefficient (R2) between 0.98 and 1.00. The expression level of each gene in each sample was measured in triplicate and calculated as 2-ΔCT where CT is the cycle threshold of the gene of interest normalized against actin-7 gene (XM_021777064.2). Statistical analysis was conducted by ANOVA pair test and Duncan’s test for multiple tests.ResultsSymptom evaluationTo assess the resistance levels between the tolerant (R72) and susceptible (HN) cassava cultivars, disease symptoms were evaluated at the inoculation site, located between the third and fourth fully expanded leaves of the infected plants. Visible symptoms were scored 7 days after inoculation (DAI) in both cultivars, with HN showing a significantly higher symptom score (Fig. 1a). Disease progression in HN was further evidenced by leaf wilting observed at 28 DAI. Symptom scores between the two cassava cultivars began to diverge significantly at 7 days after inoculation (DAI). R72 did not exhibit necrosis at the inoculation site with a score exceeding 1. Disease progression over time was assessed by calculating AUDPC values, as shown in Fig. 1b. A significant difference (p 0.05). Asterisks indicate statistical significance based on ANOVA between cultivars with the same DAI: *p ≤ 0.05, **p ≤ 0.01, ***p ≤ 0.001. (b) The area under the disease progression curve (AUDPC) is the mean value ± standard deviation of five biological replications. Statistical analysis was performed using ANOVA to compare the differences between cultivars, and values indicated by the same letter are not significantly different (p > 0.05).Full size imageProteomic data analysisTotal protein was extracted from the fourth fully expanded leaf under the inoculation point at each DAI, with five biological replicates per DAI. The mass spectra were matched to Manihot esculenta in UNIPROT database. The datasets analyzed during the current study are available in the jPOST repository and ProteomeXchange repository, JPST003734 (https://repository.jpostdb.org/entry/JPST003734) and PXD062377 (https://proteomecentral.proteomexchange.org/cgi/GetDataset?ID=PXD062377). The intensity of each protein was normalized to the total intensity of each individual sample to obtain the percentage expression ratio explaining the protein expression. Sample variation was visualized using principal component analysis (PCA) (Fig. 2) and showed a distinct separation between R72 and HN, with R72 samples located on the negative x-axis and most HN samples on the positive x-axis. Additionally, samples with the same DAI were clustered together in the PCA plot. PCA provided an overview of the protein expression in our samples related to DAI and cultivars, highlighting the different protein expression profiles among the two cultivars that were induced by bacterial infection. To identify proteins responsive to XamC infection, the log2 fold changes in protein expression between the R72 and HN cultivars were calculated at each day after inoculation (DAI) and visualized using volcano plots (Fig. 3). Proteins with an FDR 2-fold change (log2 fold change ≥ 1 or ≤ −1) were considered as proteins with significantly changed expression. High protein counts were found at the start of the sample collection date (0 DAI), indicating a background difference between the two cassava cultivars, which should be subtracted from other DAIs. R72 exhibited the highest number of upregulated proteins (1,425) at 7 DAI. This time point, which corresponds to the initial appearance of disease symptoms, may represent a critical phase in the protein response to XamC in the tolerant cultivar. In contrast, HN showed the highest number of upregulated proteins (2,979) at 28 DAI. Differentially abundant proteins in R72 were further analyzed to identify protein groups potentially involved in the XamC response.Fig. 2Principal component analysis plot of every sample. Squares indicate R72, and circles indicate HN. Color labels are distinguished by DAI: blue (0 DAI), green (3 DAI), yellow (7 DAI), red (14 DAI), and purple (28 DAI).Full size imageFig. 3Volcano plot of differential protein expression of every DAI protein, with log2 fold changes as the x-axis and -log10 (FDR) of significat value as the y-axis.Full size imageDifferentially abundant proteins clusteringDifferentially abundant proteins in the tolerant cultivar were studied to understand the defense mechanisms against infection of Xam in cassava. First, significantly differentially abundant proteins at 0 DAI were subtracted from the list of proteins after 3 DAI, which could clarify the dynamics of protein changes in response to Xam infection in the tolerant cultivar. The trend of protein expression for significantly differentially abundant proteins at each DAI was classified into seven clusters using KMC, based on Gap-statistics analysis. Seven clusters were represented, and the mean expression values are shown in Fig. 4. The expression trends differed, with some peaking after infection at 3 DAI (clusters 1, 2, and 6), at 7 DAI (clusters 4 and 7), and at 14 DAI (cluster 5). However, clusters 4 and 7 were particularly interesting, as disease symptoms were first observed at 7 DAI in R72. The subcellular localization of each protein in the seven clusters is presented (as a percentage) in Figure S1. Most subcellular localization predictions for each cluster indicated a high proportion of nucleus-expressed proteins, with the highest proportion in cluster 4. However, the distinct results of cluster 7 showed that the majority of proteins were predicted to be in the plastid. Based on this information, it can be inferred that cluster 7 showed distinct behavior compared to the other clusters and this might be the most responsive cluster in the tolerant cultivar related to symptom progression.Positive values in x-axis indicate upregulated protein in R72 (tolerant) and negative value in x-axis indicate upregulated protein in susceptible HN. Colored dots represent proteins that are upregulated proteins in a particular cultivar (fold change > 2) and significant value (FDR