Transcriptome analysis reveals potential genes associated with glyphosate resistance in cassava

Glyphosate, widely used to manage weeds in cassava crops, simultaneously inhibits cassava growth, necessitating the development of herbicide-tolerant cassava varieties. In this study, screened 262 cassava varieties, identifying the glyphosate-resistant (GR) variety ZM8701 and the glyphosate-sensitive (GS) variety SC9. Transcriptomic analysis via Illumina sequencing revealed differentially expressed genes associated with resistance, including Cytochrome P450, GST, GT, ABC transporters, and others such as MIOX1, LHCA1, PPH, HSP26, HSP83A, and UGT73C5. Notably, the EMB3004 gene, involved in the biosynthesis of aromatic amino acids, was significantly upregulated in resistant varieties, suggesting a key role in countering glyphosate’s inhibition of the shikimic acid pathway. These genes are pivotal in enhancing cell wall biosynthesis, optimizing photosynthesis, and improving detoxification processes. This research elucidated the molecular mechanisms underlying cassava’s resistance to glyphosate, thereby laying the groundwork for breeding programs aimed at developing herbicide-resistant varieties.


Introduction
Cassava, a tropical root crop, is predominantly cultivated in regions characterized by poor, arid soils, and serves as a staple food for approximately 800 million people globally.Compared to many other crops, cassava exhibits slow initial growth.Although loose soil is favorable for cassava growth, it concurrently facilitates intense competition from weeds.Furthermore, the wide spacing between cassava plants permits the emergence and competition of weeds for sunlight, water, and nutrients [1][2][3].In East Africa, weeds frequently pose more significant threats than pests or diseases, resulting in yield reductions of approximately 50% [4].
Weeds significantly contribute to decreased crop yields, necessitating weed management as an essential component of practical agricultural production.Chemical herbicides represent the most widely used and cost-effective method for weed control [5].Glyphosate, a distinctive herbicide recognized for its broad-spectrum control, low toxicity, and minimal soil residue, has played a pivotal role in global food production [6].Since its commercial introduction in the early 1970s, glyphosate has consistently expanded its market share, establishing itself as the most extensively used herbicide worldwide [7].From 1996 to 2014, the application of glyphosate increased 15-fold [8].
Glyphosate competitively inhibits the binding of phosphoenolpyruvate (PEP) to 5-enolpyruvylshikimate-3-phosphate (EPSPS), forming a stable ternary complex that subsequently blocks the reaction of PEP with 3-phosphoshikimate (S3P).The primary mechanism of glyphosate's toxicity arises from its stringent inhibition of EPSPS activity, which leads to the obstruction of branched enzyme synthesis [9].This disruption significantly alters normal nitrogen metabolism, thereby impeding the biosynthesis of aromatic amino acids.Suppression of EPSPS function can lead to the rapid accumulation of shikimic acid within the plant, which subsequently hinders photosynthesis and leads to plant wilting and death [10].
Herbicide tolerance mechanisms are categorized into target-site resistance (TSR) and non-target-site resistance (NTSR) [11].Target-site resistance occurs via genetic alterations that render the target protein/enzyme immune to herbicides, thus preventing interference with its function and shielding plants from damage.This mechanism is comparatively well-characterized and typically involves traits governed by single-gene inheritance [12,13].Conversely, non-target-site resistance mechanisms, which do not involve the target protein, may arise when herbicides fail to sufficiently reach the site of action (SoA), thereby inadequately affecting plant physiological processes and permitting plant survival [14].Non-targetsite resistance generally entails alterations in herbicide absorption, translocation, or metabolism, frequently involving genes from extensive gene families such as cytochrome P450 (P450) and glutathione S-transferase (GST) [15,16].
To date, limited research has focused on the mechanisms of glyphosate resistance in cassava.Aaron W. Hummel et al. found that leaves expressing the enzyme under the control of the 2×35S promoter had significantly higher EPSPS activity than leaves expressing the same enzyme under the control of the native cassava EPSPS promoter.Additionally, amino acid substitution variants demonstrated enhanced glyphosate tolerance compared to the wild type.Plants expressing the enzyme with the 2×35S promoter sustained less glyphosate-induced damage than those with the native EPSPS promoter [17].However, the molecular mechanisms underlying cassava's resistance to glyphosate are still not well understood.NTSR mechanisms typically protect resistant plants from irreversible cellular damage and must remain active beyond the effect duration of glyphosate [16].Moreover, stress-responsive genes are subject to diurnal regulation; thus, transcriptional sampling during off-peak periods might not fully capture the plants' responses to herbicide exposure [18,19].Shikimic acid content may serve as a critical indicator of plant responses to glyphosate exposure.Plant damage typically initiates within 3-8 h following herbicide application [20].Analysis of RNA-seq data from glyphosate-resistant (GR) and glyphosate-sensitive (GS) cassava varieties enhances our understanding of the molecular mechanisms behind cassava's resistance to glyphosate.
This study screened 262 cassava varieties for glyphosate tolerance and identified glyphosate-resistant (GR) and glyphosate-sensitive (GS) varieties using physiological and biochemical indicators.The objective of this study was to elucidate the physiological mechanisms underlying cassava's resistance to glyphosate.Additionally, transcriptome sequencing was utilized to identify further candidate genes involved in herbicide resistance, thereby providing a reference for enhancing the efficiency of breeding for herbicide resistance in cassava.

Determination of cassava resistance levels to glyphosate in the field
To assess cassava's resistance to glyphosate, 262 varieties were cultivated at the Danzhou experimental site of Hainan University (N19° 31′ E109° 34′), with 3 biological replicates for each.Upon reaching the thirty-leaf stage, the total number of leaves (N1) and the number of withered leaves (N2) per plant were recorded before treatment.Subsequently, the plants were treated with 1000 g ai/ha of glyphosate (41% glyphosate isopropylamine salt; Roundup, Shanghai, China) using a high-pressure sprayer (KOMAX, Germany).Two weeks post-treatment, the number of withered leaves per plant (N3) was calculated as the percentage of dead leaves.Resistance levels were classified based on the proportion of withered leaves as follows: ≤ 10%, highly resistant; 10-30%, moderately resistant; 30-50%, low resistance; and ≥ 50%, sensitive.

Glyphosate dose-response experiments
To assess the response of cassava varieties to escalating doses of glyphosate, both resistant and sensitive varieties were selected for cultivation in pots.Resistant seedlings at the ten-leaf stage were treated with glyphosate at doses of 0, 125, 250, 500, 1000, and 2000 g ai/ha, whereas the sensitive variety received doses of 0, 31.25,62.5, 125, 250, and 500 g ai/ha.Each treatment was conducted in three replicated pots, and the entire experiment was independently replicated twice.The fresh weight of the aboveground parts of the plants was recorded 4 weeks post-treatment.Since no significant difference (p > 0.05) was observed between the replicates, the data were combined to calculate the 50% growth reduction (GR 50 ) using logistic regression analysis in Sigma Plot 14.0 software.
In the model, y is the fresh weight response at the glyphosate dose x, C represents the lower limit, D represents the upper limit, and b is the slope around GR 50 .

Shikimic acid accumulation experiments
Shikimic acid levels were measured following the protocols specified in the shikimic acid ELISA Kit.This assay employs a microtiter plate coated with purified shikimic acid antibodies for the quantification of shikimic acid levels.Purified plant shikimic acid (HCT) antibodies were applied to the wells of the microplate to create a solidphase antibody.Plant shikimic acid (HCT) was introduced into the wells, followed by HRP-labeled shikimic acid (HCT) antibodies, forming an antibody-antigenenzyme conjugate complex.Following thorough washing, the substrate TMB was added to develop color.TMB undergoes conversion to a blue color through the catalytic action of the HRP enzyme, subsequently changing to a yellow color upon acid addition.The color intensity directly correlates with the concentration of plant shikimic acid (HCT) in the samples.Absorbance (OD values) was measured at a 450 nm wavelength with an ELISA reader, and the concentration of plant shikimic acid (HCT) was quantified by referencing a standard curve.

Expression analysis of EPSPS gene
Two functional EPSPS genes, Manes.05G046900 and Manes.01G266800,have been identified in cassava.Leaf samples were collected from plants treated with 1000 g ai/ha of glyphosate at 0, 3, and 6 days post-treatment.RNA was extracted using the TIANGEN RNA Extraction Kit (Catalog Number: DP441).The extracted RNA was reverse transcribed into cDNA with the Thermo RevertAid ™ PreMix (Catalog Number: M1631).Primers for gene expression analysis were designed via NCBI Primer-BLAST.The forward and reverse primers for Manes.05G046900 are GGT TGT GGC GGT CAT TTT CC and ATT TCC ACC TGC TGC CGT AA, respectively.Similarly, the forward and reverse primers for Manes.01G266800 are TTT CAG CTT CAG TCG CCA CA and CGA TTG GAC AGC GAC TTG GA, respectively.The Actin gene served as the reference gene for normalization.For RT-qPCR analysis, the reaction mixture contained 1 μL of each primer, 2 μL of cDNA, 6 μL of sterile water, and 10 μL of TBGREEN.

RNA-Seq
Varieties exhibiting glyphosate resistance (R) and sensitivity (S) were cultivated to the ten-leaf stage, after which they were treated with glyphosate 1000 g ai/ha, as previously described.Six days post-treatment, leaf tissues from both treated and untreated plants were collected and immediately frozen in liquid nitrogen prior to RNA extraction.The experiment included two treatments: the untreated control (RC and SC) and the glyphosate treatment (RT and ST), with each treatment replicated three times.
RNA was extracted using the RNA prep Pure Kit (TIANGEN, China) according to the manufacturer's instructions.Leaf tissues were pulverized into a powder with liquid nitrogen in an RNA-free mortar.The powdered tissue was subsequently combined with the RNA prep Pure reagent for RNA extraction.To optimize the preservation of both non-coding and coding RNA, rRNA was depleted from the sample.The RNA was fragmented into random pieces, which served as templates for the initial DNA strand synthesis using six-base random primers.The second cDNA strand was synthesized by incorporating dNTPs, buffer, DNA polymerase, and RNase H. Subsequently, the total cassava RNA was ligated with 5′ and 3′ RNA linkers.PCR amplification was carried out, and the resultant products were purified by excising the 18-30 nt band from the agarose gel electrophoresis.Index-coded samples were clustered using the cBot system with the TruSeq PE Cluster Kit v3-cBot-HS (Illumina).Sequencing was performed on the Illumina platform to generate reads.
RNA was extracted using the RNA prep Pure Kit (TIANGEN, China) according to the manufacturer's instructions.Leaf tissues were pulverized into a powder with liquid nitrogen in an RNA-free mortar.The powdered tissue was subsequently combined with the RNA prep Pure reagent for RNA extraction.To optimize the preservation of both non-coding and coding RNA, rRNA was depleted from the sample.The RNA was fragmented into random pieces, which served as templates for the initial DNA strand synthesis using six-base random primers.The second cDNA strand was synthesized by incorporating dNTPs, buffer, DNA polymerase, and RNase H. Subsequently, the total cassava RNA was ligated with 5′ and 3′ RNA linkers.PCR amplification was carried out, and the resultant products were purified by excising the 18-30 nt band from the agarose gel electrophoresis.Index-coded samples were clustered using the cBot system with the TruSeq PE Cluster Kit v3-cBot-HS (Illumina).Sequencing was performed on the Illumina platform to generate reads.
Enrichment and refinement analyses of Gene Ontology (GO) pathways were performed using the topGo package (version 2.12.0), utilizing the elim method to precisely identify GO terms associated with differentially expressed genes (DEGs).For KEGG pathway analysis, the R package clusterProfiler was employed, using a hypergeometric test to determine the statistical significance of DEGs enriched in these pathways.Pathways were deemed significantly enriched if they exhibited a corrected q-value of less than 0.05, ensuring that the results were statistically valid and biologically meaningful.
The metabolic process of crops to herbicides typically encompasses four stages: activation, binding, transport, and degradation.This process is predominantly mediated by gene families including P450s, GSTs, GTs, and ABC transporters.In this study, we identified genes differentially expressed between glyphosate-resistant and glyphosate-sensitive cassava populations before and after treatment.We rigorously analyzed differences in gene expression levels, statistical significance, and functional implications, ultimately pinpointing candidate genes linked to glyphosate metabolism resistance in cassava.
Candidate DEGs were selected based on their statistically significant up-regulated expression (fold change > 2) in the comparisons between ST vs. RT and RC vs. RT.Subsequently, these selected DEGs were annotated as herbicide metabolism genes.To confirm the expression levels of these candidate genes, both the transcriptome sequencing samples and their corresponding parallel samples were utilized.
Primers for the candidate genes were designed with Primer Premier 5.0 (Table S2), using cassava actin as the reference gene.In this study, the first cDNA strand was synthesized from total cassava RNA using Tiangen Biochemical's FastKing One-Step RT-qPCR Kit.Subsequent qRT-PCR analysis was conducted using Takara's Tli RNaseH Plus Reagent Kit.Each sample was analyzed in triplicate to guarantee repeatability and accuracy of the results.The relative expression levels of the selected genes were quantified using the 2 −ΔΔCT method, involving three biological and two technical replicates per sample.

Identification of glyphosate resistance in 262 cassava varieties
The recommended field dosage of glyphosate, 1000 g ai/ ha, was utilized to determine the percentage of defoliated leaves in cassava after 7 days of treatment.Among the evaluated cassava varieties, ten exhibited notable resistance, characterized by withered leaf percentages of ≤ 10%.Additionally, 82 varieties showed moderate resistance, with withered leaf percentages ranging from 10 to 30%.Furthermore, 73 varieties exhibited low resistance, with withered leaf percentages between 30 and 50%.Additionally, 97 varieties were categorized as sensitive, displaying withered leaf percentages of ≥ 50% (Table 1, Table S1).Among these, ZM8701 had the lowest percentage of defoliated leaves at 3.32%, while SC9 had the highest, reaching 86.48%.

The whole-plant response
The results of the logistic regression analysis of the inhibition rates for ZM8701, C3, and SC9, 2 weeks postglyphosate treatment (Fig. 1).Additionally, the GR 50 values for ZM8701, C3, and SC9 in response to glyphosate have been calculated and are detailed in Table 2. Following glyphosate treatment, the GR 50 value for ZM8701 was 3150 g ai/ha, significantly exceeding that of C3 (low resistance) at 927.9 g ai/ha and SC9 (sensitive) at 181.3 g ai/ha.ZM8701 exhibited a 3.2-fold greater resistance to glyphosate than the standard field application rate of 1000 g ai/ha and a 7.4-fold greater resistance relative to SC9, indicating a significantly elevated level of glyphosate resistance.

Accumulation of shikimic acid
Analysis of cassava leaves treated with 1000 g ai/ha glyphosate showed a significant increase in oxalic acid content post-treatment (Fig. 2).Shikimic acid levels in the leaves of GS (SC9) and GR (ZM8701) initially increased, then declined over time.Notably, there was a significant difference in shikimic acid levels between GS and the control group (CK) during this period, whereas the difference for GR was not statistically significant.On the 6th day, GS peaked at 120.17 ng/L, exceeding CK by 3.26 times, a significant increase.Conversely, GR peaked at 45.59 ng/L, a 1.20-fold increase over CK, but this difference was not statistically significant.After 6 days, shikimic acid levels in GS continued to decline yet remained significantly higher than those in the control group.Conversely, shikimic acid levels in GR showed no statistically significant difference from the control group.

Expression levels of the EPSPS gene
Using RT-qPCR, we assessed the expression levels of the EPSPS gene in the glyphosate-resistant cassava variety ZM8701 and the glyphosate-sensitive variety SC9.Prior to treatment, no significant differences were observed in the mRNA expression levels of the EPSPS genes (Manes.01G266800,Manes.05G046900) between ZM8701 and SC9.However, following glyphosate application, the expression levels of the EPSPS genes in SC9 significantly increased by day 3, while in ZM8701, they significantly decreased.By day 6, the expression levels of Manes.01G266800 had decreased in both glyphosateresistant (GR) and glyphosate-sensitive (GS) samples relative to the control.However, the expression levels of Manes.05G046900exhibited no significant changes compared to the control group (Fig. 3).

Illumina sequencing and de novo assembly
Using Illumina sequencing technology, we obtained a total of 499,028,120 raw sequencing reads from the 12 RNA libraries (RC, RT, SC, and ST, each with three biological replicates).Following quality control and data refinement, 495,948,550 clean reads were retained, ranging from 49,270,938 to 60,784,678 per sample, suitable for subsequent de novo assembly (refer to Table 3).After this screening, the Q20 ratio exceeded 96.9%, the Q30 ratio surpassed 92%, and the GC content was over 42.41%, confirming the reliability of the transcriptome data for further analysis.

Sample-to-sample distribution of gene expression
To analyze the distribution of gene expression between the resistant variety ZM8701 and the susceptible variety SC9, we utilized RSEM software for quantifying gene expression levels.The resulting data were normalized using FPKM values and presented in distribution graphs (Fig. 4) to visualize the patterns of gene expression among biological replicates under varying conditions.This approach offers a clear and intuitive understanding of gene expression distribution.

Differential gene expression and functional analysis
Fold change (FC) represents the ratio of expression levels between two sample groups.The False Discovery Rate (FDR) is derived by adjusting the significance of p-values.
The criteria for selecting differentially expressed genes (DEGs) include an FC of ≥ 2 and an FDR of < 0.05.The differential expression of genes between glyphosateresistant and sensitive cassava varieties was analyzed.In the RC vs. RT comparison, 3162 DEGs were identified, with 1673 upregulated and 1489 downregulated in the RT group.In the SC vs. ST comparison, 1386 DEGs were found, with 699 upregulated and 687 downregulated in the ST group.In the ST vs. RT comparison, 2404 DEGs were identified, comprising 1267 upregulated and 1137 downregulated in the RT group.A total of 1058 common DEGs were identified across both the RC vs. RT and ST vs. RT comparisons (Fig. 5).
The KEGG pathways, enriched with differentially expressed genes, displayed the top 20 pathways with the lowest significant Q-values.In the comparison between RC and RT, differentially expressed genes were significantly enriched in several pathways, including metabolic pathways (425 genes), biosynthesis of secondary metabolites (268 genes), plant hormone signal transduction (60 genes), and starch and sucrose metabolism (49 genes), among others.
In the comparison between ST and RT, differentially expressed genes were predominantly enriched in several key pathways, including metabolic pathways, biosynthesis of secondary metabolites, plant-pathogen interactions, fatty acid degradation, shikimic acid biosynthesis, flavonoid biosynthesis, tryptophan metabolism, and beta-alanine metabolism.Notably, upregulated DEGs were primarily concentrated in metabolic pathways (270 genes), biosynthesis of secondary metabolites (193 genes), plant-pathogen interactions (30 genes), fatty acid degradation (20 genes), and flavonoid biosynthesis (18 genes).Additionally, 14 upregulated genes were annotated to glutathione metabolism, and another 14 to the cytochrome P450 enzymes involved in drug metabolism pathways (Fig. 7).
Based on enrichment analysis of KEGG pathways and differential gene expression analysis, we identified gene clusters with upregulated expression in the ST vs. RT and RC vs. RT comparisons, as well as in two control groups.Subsequently, we selected genes exhibiting significant differences in the RT vs. RC comparison that are involved in herbicide metabolism, including P450s, GSTs, GTs, ABC transporters, and oxidases.Considering the expression levels, statistical significance, and functional relevance of these genes, we identified 20 candidate resistance genes (Table 4).

Crop resistance screening
Glyphosate-binding sites are often highly active, which makes identifying mutations associated with glyphosate tolerance a challenging task.The resistance of crops to glyphosate is not easily observed, making the screening for glyphosate-resistant varieties exceptionally demanding [21].Tao et al. [22] conducted field treatment experiments and applied physiological and biochemical research methods to assess glyphosate resistance in 66 bean lines, identifying two lines with natural resistance.In this study, we evaluated the tolerance of 262 cassava varieties to glyphosate, identifying ZM8701 as glyphosate-resistant (GR) and SC9 as sensitive (GS) based on physiological and biochemical indicators.The aim was to elucidate the physiological mechanisms underlying cassava's resistance to glyphosate.Additionally, transcriptome sequencing was employed to identify candidate genes associated with herbicide resistance, offering valuable insights for enhancing herbicide resistance in cassava.

The impact of glyphosate on shikimic acid content in cassava
Shikimic acid content is a crucial indicator of plant response to glyphosate stress [23].After glyphosate Fig. 6 The GO function classification was performed on the annotated DEGs in the RC vs. RT and RT vs. ST comparison groups treatment, shikimic acid content in plants varies depending on the timing of sampling.Generally, maize, soybean, rice, and wheat tissues exhibit peak shikimic acid concentrations approximately 4-7 days post-treatment [24,25].Similarly, in cassava, shikimic acid content reached its peak on the 6th day after glyphosate treatment.In a study by Gao et al. [26], glyphosate spray experiments on selected soybean varieties revealed no significant difference in shikimic acid content between glyphosate-resistant (GR) and glyphosate-sensitive (GS) leaves before  treatment.After glyphosate treatment, shikimic acid levels in GS leaves gradually rose, peaking at about 9.9 times higher than in GR leaves.Furthermore, Wang et al. [27] found that the shikimic acid content in glyphosate-sensitive (GS) peanuts (C179) post-glyphosate treatment was approximately twice as high as that in glyphosate-resistant (GR) peanuts (C091).The study results show a substantial increase in shikimic acid accumulation in glyphosate-sensitive (GS) cassava (SC9) following glyphosate treatment, with levels reaching approximately 2.63 times higher than those in glyphosate-resistant (GR) cassava (ZM8701).This observation Fig. 8 The expression patterns of 20 selected genes were verified via RT-qPCR, and the correlation between RNA-seq and RT-qPCR in the leaves was verified aligns with findings in soybeans and peanuts, where glyphosate has been shown to inhibit enzyme activity in the shikimic acid pathway, resulting in increased shikimic acid content.

EPSPS gene expression analysis
In various resistant weed populations, approximately 106 mutations have been identified at different sites on the target enzyme, EPSPS.These mutations lead to the substitution of specific amino acids in EPSPS, including Ser, Thr, Ala, or Leu.Examples of resistant weeds include Eleusine indica [24], Lolium perenne L. [28], and Amaranthus tuberculatus [25].This resistance is due to a unique EPSPS gene that shows high homology with EPSPS genes in other plants.Convolvulus arvensis L. exhibits glyphosate tolerance through a unique EPSPS gene that shares significant homology with those in other plants.Introduction of the EPSPS gene from Convolvulus arvensis L. into Arabidopsis results in significant improvement in glyphosate resistance [29].Furthermore, gene duplication or promoter mutations may lead to overexpression of target genes.Initial reports indicated that EPSPS gene copy numbers in glyphosate-resistant Palmer amaranth populations from Georgia ranged from 5 to over 160 times higher than in susceptible populations [30].Subsequent studies in glyphosate-resistant multiple-flowered ryegrass populations in Arkansas showed EPSPS gene copy numbers to be 25 times higher in resistant populations compared to susceptible ones [31].This study further investigated the transcription levels of the EPSPS gene between glyphosate-resistant (GR) and glyphosatesusceptible (GS) populations using real-time fluorescent quantitative PCR.However, in the case of GR (ZM8701), no overexpression of the EPSPS gene was observed following glyphosate treatment.

Transcriptome analysis of glyphosate on cassava leaves
High-throughput omics technologies have emerged as powerful tools for unraveling plant responses to abiotic stressors [32].Recently, these techniques have been extensively used to explore plant resistance mechanisms to herbicides [33].Research on glyphosate resistance primarily focuses on genes associated with detoxification processes.Non-target site resistance (NTSR) involves the plant's detoxification process, typically including four stages: activation, binding, active transport, and degradation of herbicides [34].RNA-Seq technology is an effective tool for studying NTSR mechanisms, swiftly identifying genes and metabolic pathways linked to glyphosate resistance, thus enhancing research efficiency on cassava's resistance mechanisms to glyphosate.In this study, we conducted de novo sequencing on cassava leaf samples damaged by glyphosate over 6 days, focusing also on the gene expression differences in NTSR enzymes.
Cytochrome P450 enzymes play a crucial role in metabolizing various substances and activating genes of metabolic enzymes related to glyphosate detoxification [35][36][37].In our study, eight Cytochrome P450 genes: CYP704B1, CYP707A4, CYP71D10, CYP76A2, CYP82A1, CYP82C2, CYP93A1, and CYP94A2 were identified as key candidates for resistance.Specifically, the CYP71D subfamily exhibits distinct functional specificity, participating in the biosynthesis of indole alkaloids, flavonoids, and triterpenoids [38].The CYP71 family in plants primarily metabolizes phenylurea and sulfonylurea herbicides [39].Research has indicated that CYP94 and CYP71 genes may play a significant role in resistance to glufosinate [40].Moreover, CYP71C6v1 in wheat [41] and CYP71A10 in soybean [42] have been demonstrated to metabolize herbicides.The CYP76 subfamily is crucial in plant hormone biosynthesis, including the synthesis of secondary metabolites, hormone signal transduction, and environmental stress response [43,44].CYP76B1 from Helianthus tuberosus has been shown to metabolize herbicides within the phenylurea class [45].In wild radish, expression levels of RrCYP704C1 and RrCYP909B1 in herbicide-resistant plants were significantly higher compared to those in susceptible plants [46].CYP94A1 plays a role in defending against chemical damage in wild pea plants [47].The Cytochrome P450 genes identified in our cassava study may confer herbicide resistance, but this hypothesis requires empirical testing and validation.
Glutathione S-transferases (GSTs) may protect cells from lipid oxidation induced by oxidative stress or participate in exogenous detoxification processes.GSTs perform crucial functions by conjugating glutathione (GSH) to various molecules, including exogenous substances, during metabolic detoxification and removing these substances from the cytoplasm [48].GSTs and GTs are key enzyme families involved in herbicide detoxification [47].The biochemical response of non-target site resistance (NTSR) to herbicides is well documented [49].In our study, four GT genes: GT2 (Manes_10G121800), GT-2 (Manes_02G179200), and GT-3B (Manes_08G092900, Manes_09G106100) were identified as major candidates for resistance.In A. thaliana, enhanced GT gene activity was observed when using AHAS inhibitors to combat herbicides [50].GT genes have been shown to glycosylate the primary metabolite 6-hydroxybentazone in a herbicide-tolerant variety [51].However, the specific roles of GT2, GT-2, and GT-3B in herbicide resistance remain unclear, necessitating further research to elucidate their functions.
ATP-binding cassette (ABC) proteins are crucial to various physiological processes in plants, such as auxin polar transport, lipid catabolism, exogenous detoxification, disease resistance, and stomatal function [52].Unlike the CYP450, GST, and GT gene families that confer herbicide resistance through metabolic processes, ABC transporters neutralize herbicides and their metabolites by sequestering them [53].In this study, two instances of the gene Manes_01G234500 were annotated as ABCB11, suggesting a potential association with glyphosate resistance.Furthermore, similar findings were observed in bifidobacteria resistant to fenoxaprop-P-ethyl, where an upregulated ABCB11 gene was also identified [54].These observations suggest the possible involvement of ABCB11 in conferring metabolic resistance to glyphosate.
In the target pathway of glyphosate, the gene Manes_09G160900 (EMB3004) was significantly upregulated in the GR treatment group, while there was no significant change in the expression of this gene in the GS treatment group.The EMB3004 gene encodes two enzymes, 3-dehydroquinate dehydratase I and shikimate dehydrogenase, both crucial for the biosynthesis of aromatic amino acids.Specifically, 3-dehydroquinate dehydratase I, the third enzyme in this pathway, converts 3-dehydroquinate into 3-dehydroshikimate, while shikimate dehydrogenase, acting in the fourth step, converts shikimate into 3-dehydroshikimate.We hypothesize that this gene's upregulation represents an adaptive response to the inhibition of EPSP synthase (EPSPS) by glyphosate.By enhancing the activity of early pathway enzymes, plants may increase the production of intermediate metabolites needed for aromatic amino acid synthesis, thereby partially compensating for reduced EPSPS activity.This compensatory mechanism by upstream metabolites might help circumvent the biosynthetic blockade induced by glyphosate, thus enhancing plant tolerance to the herbicide.
Additionally, this study identified several genes, including Manes_11G030400 (MIOX1), Manes_07G104600 (LHCA1), Manes_02G095700 (PPH), Manes_01G129800 (HSP26), Manes_14G022300 (HSP83A), and Manes_08G015600 (UGT73C5), that correlate with cassava's resistance to glyphosate.Specifically, MIOX1 catalyzes the conversion of myo-inositol into d-glucuronic acid, an essential precursor for cell wall polysaccharide synthesis.This conversion may represent a defensive strategy against glyphosate-induced stress by enhancing cell wall biosynthesis.LHCA1 is part of the lightharvesting complex associated with photosystem I, essential for photosynthesis.Its primary function is to capture and transfer light energy, facilitating electron transport in the photosynthetic process.Under glyphosate stress, the elevated expression of LHCA1 in the cassava variety ZM8701 may maintain photosynthetic efficiency, supporting continuous energy production and biosynthetic processes.The PPH gene, involved in regulating chlorophyll degradation, may buffer the effects of impaired photosynthesis under glyphosate stress, preventing the accumulation of toxic intermediates.The LHCA1 and PPH genes are crucial in sustaining photosynthesis and balancing chlorophyll metabolism, thereby enhancing cassava's tolerance to glyphosate.By optimizing energy capture and utilization and regulating chlorophyll degradation, these genes contribute to cassava's adaptation to glyphosate, reducing its deleterious impacts.Additionally, the HSP26 and HSP83A genes, which encode heat shock proteins (HSPs), are critical in this context.These molecular chaperones, common in plants, assist in proper protein folding, prevent misfolding and aggregation, and support cellular recovery under stress conditions.The upregulation of these HSPs may be a key adaptive mechanism in cassava, developing resistance to glyphosate by mitigating protein misfolding and cellular stress.Furthermore, the UGT73C5 gene encodes an enzyme from the UDP-glucuronosyltransferase (UGT) family, central to detoxification.UGT73C5 enhances the solubility of xenobiotic compounds by glycosylating them, thus facilitating their excretion.This mechanism is crucial for plant survival under chemical stress.
Continued validation of these genes' roles in conferring resistance is essential.Future research should incorporate gene editing and overexpression studies to focus on the highlighted genes.Investigating how these genes interact within broader metabolic networks could reveal additional targets and strategies for developing comprehensive resistant phenotypes, thereby advancing resistance breeding in cassava.

Conclusion
This study investigates glyphosate resistance in 262 cassava varieties, identifying distinct levels of resistance based on defoliation percentages post-herbicide application.Notably, ZM8701 exhibited substantial resistance, with significantly higher glyphosate doses required to inhibit growth compared to other varieties.Through transcriptome sequencing and differential gene expression analysis, key genes involved in resistance, such as Cytochrome P450s and GSTs, were pinpointed.The study also explored the dynamics of shikimic acid accumulation, a critical response to glyphosate exposure, noting distinct patterns between resistant and sensitive varieties.Further, RT-qPCR validated the expression levels of candidate resistance genes, reinforcing the RNAseq data.Overall, this research enhances understanding of glyphosate resistance mechanisms in cassava and provides insights for developing herbicide-resistant crops through targeted genetic and breeding strategies.

Fig. 1 a
Fig. 1 a Dose-response curves of three cassava varieties to glyphosate.b Symptoms of ZM8701, C3 and SC9 12 days after treatment with 1000 g ai/ha glyphosate.The data represent the mean values, and the error bars represent the standard errors of the mean; the same applies for the following

Fig. 4 Fig. 5
Fig. 4 Gene expression density map and gene expression level violin plot of 12 groups of samples.a Gene expression abundance distribution plot.b Gene expression level violin plot

Fig. 7
Fig. 7 KEGG enrichment bar charts for the RC vs. RT and ST vs. RT comparison groups

Table 1
Susceptible level of cassava field populations to glyphosate

Table 2
Herbicide rates causing a 50% growth reduction (GR50) for the two cassava varieties

Table 3
The data table of RNA-seq of ZM8701 and SC9

Table 4
Screening of differential genes for resistance expression