Identification of candidate genes involved with dicamba resistance in waterhemp (Amaranthus tuberculatus) via transcriptomics analyses

Abstract

Waterhemp [Amaranthus tuberculatus (Moq.) Sauer] is one of the most troublesome weeds in the United States. An A. tuberculatus population (CHR) was identified in Illinois, USA, as resistant to herbicides from six different site-of-action groups. Recently, the same population was also recognized as dicamba resistant. This study aimed to identify key resistance genes and the putative dicamba resistance mechanism in A. tuberculatus via transcriptomics analysis. Multiple differentially expressed (DE) genes and co-expression gene modules were identified as associated with dicamba resistance. Specifically, genes encoding glutathione S-transferases (GSTs), ATP-binding cassette transporters, peroxidases, and uridine diphosphate (UDP)-glycosyltransferases (UGTs) were identified. Results indicated enhanced oxidative stress tolerance as the primary mechanism for reducing dicamba toxicity. Results also point to potential glycosylation via UGTs and conjugation via GSTs of dicamba and its by-products. This is the first transcriptomics characterization of dicamba resistance in A. tuberculatus. Multiple non-target-site resistance genes were identified, indicating a cross-resistance pattern in the CHR population leading to a putative-enhanced oxidative stress response. Regions of multiple DE genes (i.e., genomic hot spots) across the A. tuberculatus genome corroborate previous results and potentially add to the complexity of non-target-site resistance traits.

Introduction

Proper weed management is a crucial step for achieving high yields in agriculture and is constantly facing many challenges, including the evolutionary phenomenon of herbicide resistance. Chemical control is a heavily relied upon weed management tool, creating intense selection pressure for the evolution of herbicide-resistant weeds (Bagavathiannan and Davis 2018).

Herbicide resistance can be classified into two main types: target-site (TS) and non-target-site (NTS) resistance. TS resistance mainly refers to when one or more nucleotide mutations in the gene encoding a protein bound by a herbicide disrupts the toxicity capabilities of the herbicide but can also be due to increased expression or copy number variation of the gene (Gaines et al. 2020). NTS resistance refers to any resistance mechanism not involving a change in the herbicide target site, including metabolic resistance, reduced translocation, and sequestration (Ghanizadeh and Harrington 2017; Jugulam and Shyam 2019; Rigon et al. 2020).

With the advent of genetically modified crops, herbicide-resistance traits became overwhelmingly present in several crops, with an 80% to nearly 100% adoption rate in the U.S. corn (Zea mays L.) and soybean [Glycine max (L.) Merr.] system (Dopson 2022). Recently, soybean cultivars with a dicamba resistance trait were released, increasing the usage of this specific herbicide (Behrens et al. 2007). Dicamba is a synthetic auxin herbicide (HRAC Group 4) commercialized since the 1960s (Egan and Mortensen 2012), and it has recently become one of the most used herbicides in the United States with the release of dicamba-tolerant cultivars. For instance, for soybean and cotton (Gossypium hirsutum L.) production, dicamba is among the top five most-used herbicides, with 3.76 and 1.08 million kg of dicamba applied in the last 5 yr, respectively (USDA-NASS 2022). As observed after the release of glyphosate-resistant cultivars (Heap and Duke 2018), it is predicted that the overuse of dicamba could lead to a substantial selection of dicamba-resistant weed populations.

One of the main targets for dicamba application is one of the most troublesome weed species in the midwestern United States: waterhemp [Amaranthus tuberculatus (Moq.) Sauer]. This weed is a dioecious, herbaceous species, highly prolific and with persistent seed in the soil seedbank, and it can reduce yields in soybean and corn by 40% and 70%, respectively (Hager et al. 2002; Korres et al. 2018; Steckel and Sprague 2004). An A. tuberculatus population from Champaign County, Illinois, USA (designated CHR) was identified as resistant to herbicides from multiple site-of-action groups, including synthetic auxins (2,4-D) and inhibitors of acetolactate synthase, protoporphyrinogen oxidase, 4-hydroxyphenylpyruvate dioxygenase (HPPD), photosystem II, and very-long-chain fatty-acid elongases (Evans et al. 2019; Strom et al. 2020).

Subsequently, the CHR population was characterized as dicamba resistant, despite no history of repeated dicamba selection, with the trait identified with moderate heritability and likely being multigenic (Bobadilla et al. 2022). Previous studies investigated gene expression patterns of CHR related to 2,4-D and HPPD-inhibitor resistance, indicating an NTS resistance case involving cytochrome P450 (CYP450), ATP-binding cassette (ABC) transporter, and uridine diphosphate (UDP)-glycosyltransferase (UGT) as the primary candidates for resistance (Giacomini et al. 2020). Interestingly, genomic hot spots, that is, clusters of differentially expressed (DE) genes, were also observed. This scenario raises the question of whether other resistance traits in CHR share similar patterns and candidate resistance genes. CHR also was the first reported case of dicamba resistance in A. tuberculatus, and no prior study investigated the expression patterns of dicamba resistance in this species.

The goal of this study was to characterize the gene expression patterns and identify candidate genes for dicamba resistance in the CHR population. An RNA-seq study was conducted comparing dicamba-resistant and dicamba-susceptible individuals within a segregating F₂ population derived from the CHR population.

Material and Methods

Plant Material and Sequencing

For the RNA-seq experiment, a previously characterized pseudo-F₂ population was used (Bobadilla et al. 2022). Briefly, CHR plants were treated with 560 g ai ha⁻¹ of dicamba, and the most resistant was selected to cross with a universal sensitive population (WUS) to generate an F₁ population, from which a pseudo-F₂ line was generated. An F₂ line was used for RNA-seq to reduce background noise and genetic differences between resistant and sensitive individuals, focusing only on the genes segregating based on the selected trait (Giacomini et al. 2018). All seeds were subjected to a 50% commercial bleach treatment for 10 min and rinsed twice with water for 10 min each. Pseudo-F₂ seeds were suspended in a 0.1 g L⁻¹ agarose solution and placed in a 4 C refrigerator for 4 wk to reduce dormancy (Bell et al. 2013).

After stratification, seeds were germinated in petri dishes containing blotting paper with 2.0 ml of water. Petri dishes were closed with sealing film to avoid water loss and placed in a growth chamber for 48 h set for 12/12 h day/night and 32/15 C temperature regimes. After germination, seedlings were transplanted into 164-cm³ Cone-tainers (Ray Leach SC10 “Cone-tainer,” Tangent, OR). Plants were kept under a mist bench programmed to water plants twice a day. All Cone-tainers were filled with a custom sandy loam growth medium containing 1:1:1 (soil:peat:sand) and 3.3% organic matter, 6.8 pH. About 0.45 kg of slow-release complete fertilizer (Osmocote^® 13–13–13 slow-release fertilizer, Scotts, Marysville, OH) was mixed into 200 kg of the medium before planting. One week after transplant, a supplement of about 80 mg of additional Osmocote^® fertilizer was added to the top of the growth medium in each Cone-tainer. The greenhouse was kept in a temperature and light regime of 28/22 C and 16/8 h, respectively, with sunlight supplemented with metal-halide lamps.

Four-hundred pseudo-F₂ plants were treated with 560 g dicamba ha⁻¹ (XtendiMax^®, Bayer CropScience, St Louis, MO) when plants had 6 to 7 leaves and were 7 to 10 cm in height. Before treatment, a single, young, fully formed leaf from each plant was collected, promptly frozen in liquid nitrogen, and stored at −80 C for later extraction after phenotyping and plant selection. Delimiting rate definition and phenotyping approach were conducted as previously described (Bobadilla et al. 2022). Briefly, plants were phenotyped 16 d after treatment based on visual damage estimation, biomass, and plant area reduction (Figure 1 A and B). After phenotyping, 16 individual pseudo-F₂ plants (8 resistant and 8 sensitive) were chosen for further RNA extraction.

Figure 1. Plant selection and phenotype classification for RNA-seq. Photos show the differences in phenotypes of some of the individuals selected for sequencing: (A) resistant plants and (B) sensitive plants. Selection was done based on visual damage estimation, biomass, and plant area measured via image analysis (Bobadilla et al. 2022). Photos were taken 14 d after treatment with dicamba at 560 g ai ha⁻¹. The graph shows the relationship between biomass and plant area across resistant and sensitive individuals.

RNA from the 16 selected individuals was extracted using a Trizol-based method (Simms et al. 1993) followed by a DNAse I treatment. Samples were checked for quality and quantity via Qubit analyzer (Thermo Fisher Scientific, Waltham, MA), gel electrophoresis, and Bioanalyzer (Agilent Technologies, Santa Clara, CA). Samples were sent to the Roy J. Carver Biotechnology Center at the University of Illinois, Urbana-Champaign, for Illumina library construction and sequencing. Libraries were prepared using an Illumina TruSeq Stranded mRNAseq Sample Prep kit (Illumina, San Diego, CA). Single reads (100 bp) were sequenced using a NovaSeq 6000 system with one NovaSeq SP flow cell.

Data Preprocessing

All fastq files were filtered for low-quality reads, and adapters were trimmed using Fastp (Chen et al. 2018). Reads were further filtered to remove rRNA using the software SortMeRNA (Kopylova et al. 2012). A genome-guided transcriptome assembly was conducted with the trimmed and filtered data. Reads were mapped to the A. tuberculatus draft genome (Montgomery et al. 2020) using STAR (Dobin and Gingeras 2015), and all bam files were combined using the SAMtools merge tool (Li et al. 2009). A merged bam file was input into Trinity for the transcriptome assembly (Table 1). Functional transcriptome annotation was done using the Trinotate pipeline (Bryant et al. 2017), and the final transcriptome was generated using the annotated transcripts. Transcriptome quality was assessed using BUSCO (Benchmarking Universal Single-Copy Orthologs; Simão et al. 2015) completeness score based on the viridiplantae_odb10 database. Overall quality was accessed by mapping a set of samples to the assembled transcriptome with Bowtie2 (Langdon 2015) and obtaining the N50 parameter using a Trinity custom script (https://github.com/trinityrnaseq/trinityrnaseq/blob/master/util/TrinityStats.pl). The transcriptome was later mapped back to the genome to identify genomic locations using GMAP (Wu and Watanabe 2005). All reads were pseudo-aligned to the assembled transcriptome using Kallisto for read count quantification (Bray et al. 2016).

Table 1. Genome-guided transcriptome assembly statistics.

^a BUSCO, Benchmarking Universal Single-Copy Orthologs.

^b N50, sequence mean length.

Differential Expression

A differential expression analysis was conducted using the R package EdgeR (Robinson et al. 2010). Quantification files generated from Kallisto were loaded into R using the R package Tximport and summarized to the gene level (Soneson et al. 2015). Genes with low expression in more than 80% of samples were filtered, and the remaining genes were subjected to normalization via the trimmed mean of M-values method. Common, trended, and tagwise dispersions were estimated using a quantile-adjusted conditional maximum likelihood. A negative binomial model was fit to the data, and the exact test was conducted based on the fitted dispersion to identify DE genes. A 5% false discovery rate (FDR) was used to identify DE genes. To identify significant terms within the list of DE genes, overrepresented and conditional GO-term enrichment analysis was conducted using the TopGO R package (Alexa and Rahnenführer 2009).

Weighted Co-expression Network Analysis (WGCNA)

A WGCNA was conducted using the WGCNA R package (Langfelder and Horvath 2008). Such analysis can provide valuable information regarding which genes are co-expressed, creating a venue to better interpret the tested trait. Gene expression counts were normalized using the DESeq2 method and converted to a log₂ scale (Love et al. 2014). Hierarchical clustering was used to identify outlier samples, and the constant-height tree-cut function was used to remove the outliers. Soft-threshold power was estimated to approximate network topology to a free-scale model. A signed adjacency matrix was calculated via bi-weight midcorrelation and a signed topological overlap matrix by dissimilarity. Genes were clustered via hierarchical clustering, and the dynamic tree-cutting algorithm was used to separate genes into modules. Module eigengenes were calculated to merge similar modules and identify modules associated with dicamba resistance. An intra-modular analysis was conducted to identify hub genes within modules associated with dicamba resistance. Genes were considered hubs of a module based on their module membership (0 to 1 scale according to the overall connectivity) and their gene-trait significance (0 to 1 Pearson correlation between expression and the trait). GO-term enrichment analysis was conducted for each significant module using TopGO.

Promoter Analysis

A regulatory prediction analysis was conducted to explore the regulatory nature of the identified DE genes. The promoter regions from all DE genes were extracted using the A. tuberculatus reference genome. The promoter region was defined as 2,000 bp before the transcription starting site. Promoter regions were extracted using BEDtools software (Quinlan and Hall 2010). The regulatory prediction was conducted via a transcription factor enrichment analysis using the PlantRegMap tool in the Plant Transcription Factor Database (PlantTFDB) utilizing the transcription factor motifs from Arabidopsis thaliana (Tian et al. 2020).

Quantitative PCR Validation

A quantitative PCR (qPCR) assay was conducted to validate results from the RNA-seq experiment. Five genes (PEROX12, GST-nt, GST-ct, UGT85A, and ABC10) were selected as our validation genes. Primers were designed aiming at an annealing temperature of 60 C, amplicon size of 80 to 180 bp, and targeting an exon/exon junction (Supplementary Table 1). Initial primer efficiency was tested by conducting a qPCR assay using a five-step cDNA dilution. The genes GADPH and EF-alpha were used as the housekeeping genes to estimate efficiency and relative expression. A subset of pseudo-F₂ plants was selected (six resistant and six susceptible plants), including individuals that were and were not used in the RNA-seq experiment. Only primer sets ranging from 95% to 100% efficiency were selected (Supplementary Figure 1). RNA was converted to cDNA using a ProtoScript II First Strand cDNA Synthesis Kit (New England Biolab, Ipswich, MA). qPCR was performed using a QuantStudio 7 (Thermo Fisher Scientific, Waltham, MA) in triplicate for each sample for each primer set by combining 5 μl of iTaq Universal SYBR Green Supermix (Bio-Rad Laboratories, Hercules, CA), 0.5 µl of forward primer (10 µM), 0.5 µl of reverse primer (10 µM), 3 µl of nuclease-free water, and 1 µl of cDNA. Relative expression was calculated using the 2^−ΔΔCt method (Livak and Schmittgen 2001). The entire experiment was repeated to ensure consistency. Expression results were subjected to a two-sided t-test.

Temporal qPCR Expression Assay

Temporal qPCR analysis was conducted to quantify the differential expression after dicamba application of genes confirmed as DE by qPCR validation. A total of 200 pseudo-F₂ plants were grown under the same conditions as previously described, and 48 plants (24 resistant and 24 susceptible) were selected. Tissue from the youngest fully formed leaf was collected at four periods: 0, 0.5, 4, and 48 h after application, with a different set of 12 plants (6 resistant and 6 susceptible) sampled during each period. Tissue collection, RNA extraction, and cDNA conversion were performed as previously described. All plants were treated with dicamba at a rate of 560 g ha⁻¹. qPCR for relative expression assay was conducted utilizing the gene GADPH as the housekeeping gene. Relative expression using the 2^−ΔΔCt method and significance via t-test analysis were estimated as previously described.

Results and Discussion

Differential Expression Analysis

Differential expression analysis yielded a relatively small number of DE genes, 67, indicating that utilizing an F₂ population greatly reduced background noise from the data (Giacomini et al. 2018; Figure 2). The fact that this comparison was made in the absence of any herbicide application also accounts for the small number of DE genes. From those 67 genes, 49 were functionally annotated (Supplementary Table 2), including gene family members previously characterized with involvement in NTS resistance, such as glutathione S-transferases (GSTs), UGT, ABC transporters, and CYP450. GO-term enrichment analysis indicated that DE genes were mainly enriched for terms related to response to oxidative stress and multiple other metabolism-related terms such as glutathione and UDP−glucose metabolic processes (Figure 3). Major candidates identified included a GST from the phi family with a C-terminal domain (GST-ct), a GST from the tau family with an N-terminal domain (GST-nt), a peroxidase (PEROX12), and multiple UGTs, including UGT85A and UGT2. Other important putative genes for dicamba resistance were also identified, including CYP76A and the gene ACS10, which were previously characterized to respond to synthetic auxins in other species (Gleason et al. 2011; Johnston et al. 2020; Tsuchisaka and Theologis 2004).

Figure 2. Volcano plot of all genes with key differentially expressed genes highlighted. Major genes with potential involvement in dicamba resistance are labeled according to their homologous UniprotKB ID. Genes in red and blue were significantly up- and downregulated, respectively, in dicamba-resistant relative to sensitive plants. The y axis refers to −log₁₀ false discovery rate (FDR), and the x axis refers to the log₂ expression fold change (FC).

Table 2. Quantitative PCR (qPCR) validation assay results.

^a Average relative expression calculated in relation to two housekeeping genes (GADPH and EF1-α).

^b P-value refers to two-sided t-test significance. NS, nonsignificant P-value.

Figure 3. Biological process GO-term enrichment analysis. Circle size represents the significance of overrepresented enrichment, and color gradient represents the significance of conditional enrichment. The x axis represents the number of genes annotated with each GO-term in the y axis.

Differential expression results indicated a scenario where the dicamba resistance is potentially caused by an enhanced response to oxidative stress via GSTs and peroxidases and by metabolism of dicamba and abscisic acid (ABA) via UGT glycosylation (Mateo-Bonmatí et al. 2021). The previously characterized involvement of glycosylation in regulation of indole-3-acetic acid (IAA) and results pointing to UGT genes being overexpressed after dicamba treatment could allow one to extrapolate a hypothesis that glycosylation could also be involved in dicamba detoxification. DE genes were mapped to the A. tuberculatus genome to identify their genomic position and look for clusters (Figure 4). Interestingly, some similarities were noticed with what was observed by Giacomini et al. (2020) in a previous study with a subpopulation also derived from CHR but segregating for 2,4-D and HPPD-inhibitor resistance. Multiple DE genes from the dicamba RNA-seq study were mapped to regions in chromosomes 1, 2, 4, and 16, indicating that those regions potentially contain quantitative trait loci (QTLs) important for NTS resistance in the A. tuberculatus CHR population. Genes including ABC transporter 10 (ABC10) and other transporters were found as DE and within those genomic regions in both studies, indicating that compartmentalization of the herbicides, secondary metabolites, and reactive oxygen species (ROS) may play a key role in NTS resistance in the CHR population (Rigon et al. 2020; Tani et al. 2015; Zhang and Yang 2021). From genes selected for validation, GST-nt, GST-ct, PEROX12, and ABC10 showed significant elevated expression in dicamba-resistant plants relative to the sensitive plants (Table 2), validating observed expression differences obtained from the RNA-seq experiment.

Figure 4. Genomic distribution in sliding 50-kb window plots of differentially expressed genes (DEs). The y-axis units refer to windows in mega base pairs (Mbp) Each plot represents one of the 16 pseudo-chromosomes of Amaranthus tuberculatus. Peaks represent clusters of DEs. Blue dashed lines represent previously identified hot-spot locations for 2,4-D resistance (Giacomini et al. 2020).

WGCNA

The WGCNA yielded 33 modules identifying genes with similar expression patterns within the transcriptomics data (Figure 5; Supplementary Figure 2). Module sizes range from 41 to 3,730 genes. Module eigengenes were used to identify modules associated with dicamba resistance. A GO-term enrichment analysis was conducted for each module to identify the important biological terms enriched within each. Black, Tan, and Brown modules were found enriched for response to oxidative stress. Module Grey60 was enriched for terms involved with glutathione metabolic processes and peroxisome organization, implying that genes within this module are also potentially associated with oxidative stress regulation. The Midnightblue module was enriched for multiple UGT-related terms and hormone metabolic processes, implying a potential involvement of phytohormone glycosylation (Huang et al. 2021; Mateo-Bonmatí et al. 2021; Meech et al. 2012).

Figure 5. Weighted co-expression network analysis results. (A) Gene expression dendrogram for module assignment where a total of 33 modules were identified. (B) Trait-module correlation plot with values outside parentheses representing Pearson correlation and values inside parentheses representing the significance correlation P-values. Correlation values range from −1 to 1, with red values indicating a positive association and blue values indicating a negative association with dicamba resistance. ME refers to modules followed by their color code.

WGCNA results identified multiple modules associated with dicamba resistance that were enriched for oxidative stress–related genes, furthering the hypothesis that dicamba resistance in this multiple-resistance population is due to an enhanced response to oxidative stress caused by the herbicide. The gene UGT85A was one of the main hubs within the Black module, indicating that this gene could play an important role in dicamba resistance. For the Midnightblue module, the only annotated hub gene was a gene encoding an isocitrate dehydrogenase, which has a putative role in oxidative stress response (Mhamdi and Noctor 2015; Pan et al. 2020). The Grey60 module’s only annotated hub gene encoded a peroxisomal membrane protein. Peroxisomes are organelles that sequester diverse oxidative reactions and play essential roles in metabolism, detoxification of ROS, and signaling (Mhamdi and Noctor 2015; Pan et al. 2020; Rodríguez-Serrano et al. 2014). In summary, results from WGCNA corroborate results observed in the differential expression analysis, suggesting that dicamba resistance is primarily due to detoxification of ROS generated after dicamba molecules bind to auxin receptors (Gaines 2020).

Regulatory Prediction

Results from regulatory prediction analysis showed 4,100 potential regulatory relationships among 510 transcription factors with the identified DE genes. From these regulatory relationships, 65 transcription factors showed overrepresented binding targets in the DE genes using a 0.05 q-value threshold. From these significant interactions, 12 transcription factor families were identified, with most transcription factors as members of the MYB family, followed by the TCP family (Figure 6). MYB transcription factors represent one of the largest protein families in plants and are known to be involved in plant development and responses to stresses (Wang et al. 2021). TCP is a plant-specific transcription factor family with a bHLH motif and is also known to play critical roles in abiotic stress responses (Martín-Trillo and Cubas 2010).

Figure 6. Differentially expressed (DE) genes’ regulatory prediction analysis results: (A) total number of transcription factors (TFs) enriched per family; (B) number of DE genes with motifs specific for each enriched TF family. TF family names in the x axis are according to the nomenclature in the Plant Transcription Factor Database (PlantTFDB).

Interestingly, when comparing the number of genes with motifs for the binding of these enriched transcription factors, a unique transcription factor from the CPP family shows motifs in many DE genes. Cysteine-rich polycomb-like protein (CPP) transcription factors are a small gene family that regulates plant growth, development, and stress response (Andersen et al. 2007; Song et al. 2016; Zhou et al. 2018). Many genes also contain motifs for three enriched MADS-box transcription factors. MADS-box transcription factors are involved in a wide range of developmental processes in plants, including the formation of flowers, leaves, and roots, as well as the determination of cell fate and the regulation of stress-response pathways (Theissen et al. 2000). Transcription factors are the major regulators of gene expression. Identifying how DE genes are regulated could aid in identifying the actual expression modulator. Understanding the regulatory nature of the involved genes is key to building further predictions for cross-resistance and creating strategies to minimize it (Bobadilla and Tranel 2023).

Temporal qPCR Expression

Previous studies suggested that NTS-related genes are constitutively expressed, allowing for a quick plant response after herbicide exposure (Giacomini et al. 2018, 2020). We conducted a temporal qPCR assay to validate this hypothesis and test the expression variation after dicamba application. Results indicated that, in the CHR population, all tested genes were upregulated in resistant plants relative to sensitive plants before herbicide application (Figure 7) but had varying expression patterns after dicamba application. PEROX12 and GST-nt expression increased in resistant plants 30 min after treatment, while expression of GST-ct decreased in resistant plants after dicamba exposure but remained elevated relative to sensitive plants. At 4 and 48 h after treatment, PEROX12 expression was not significantly different between resistant and sensitive plants. For ABC10, expression decreased with time in resistant plants, but continued to be elevated relative to sensitive plants. In summary, results indicated that GSTs and ABC10 were constantly upregulated in resistant relative to sensitive plants, while PEROX12 was only upregulated before and shortly after dicamba application. These results support the hypothesis that genes for NTS resistance are constitutively overexpressed in resistant plants before herbicide application (Giacomini et al. 2018).

Figure 7. Temporal quantitative PCR results for (A) GST-nt, (B) ABC10, (C) PEROX12, and (D) GST-ct before and at multiple time points after dicamba treatment. Asterisks (*) indicate comparisons that were significant (t-test P-value < 0.05), with error bars indicating variability across replicates.

Interestingly, expression of the peroxidase gene substantially increased in sensitive plants by 48 h after dicamba treatment, suggesting a native response to dicamba exposure. Dicamba was previously characterized to induce a reduction in peroxidase expression (Gleason et al. 2011). Previous studies also show that peroxidase overexpression can lead to resistance to herbicides, such as paraquat, that cause high oxidative stress (Murgia et al. 2004).

Dicamba Resistance in CHR

The A. tuberculatus CHR population was previously characterized as resistant to herbicides spanning six sites of action and containing a combination of TS and NTS mechanisms (Evans et al. 2019; Strom et al. 2020). Based on the dicamba RNA-seq results, we hypothesize the main mechanism of dicamba resistance is an enhanced ability to deal with oxidative stress, with potentially additional contribution by conjugation and glycosylation of dicamba molecules (Figure 8). Multiple DE genes potentially involved in reducing the effect of ROS in the plant were identified, including peroxidases and GSTs. ROS overproduction was previously characterized as one of the main venues whereby synthetic auxins cause damage and plant death (Christoffoleti et al. 2015; Gaines 2020).

Figure 8. Proposed dicamba resistance mechanisms in the CHR population. Currently, knowledge about the synthetic auxin effect on plants indicates an overproduction of abscisic acid (ABA), leading to a large production of reactive oxygen species (ROS) and plant death (Christoffoleti et al. 2015; Gaines 2020). The proposed resistance mechanism is that enhanced response to oxidative stress via peroxidases and glutathione S-transferases alleviates dicamba toxicity. Other putative resistance mechanisms, such as glycosylation of dicamba and ABA, are also proposed with transport via ATP-binding cassette (ABC) transporters for further degradation. Overproduction of salicylic acid is also proposed as a potential tool for alleviating oxidative stress. Created with BioRender.com.

GSTs and peroxidases are known to be agents for dealing with ROS and reducing the effect of oxidative stress (Cummins et al. 1999; Edwards et al. 2000; Kawano 2003). Peroxidases play an important role in the plant’s defense against stress, including herbicide exposure. These enzymes are involved in the detoxification of ROS, which can damage cellular components and act as signaling molecules in response to stress (Murgia et al. 2004). GSTs catalyze the conjugation of the tripeptide glutathione (GSH) to various hydrophobic, electrophilic, and often cytotoxic substrates. Plant and animal GSTs conjugate GSH with such endogenously produced electrophiles, which results in their detoxification. Some GSTs also function as glutathione peroxidases to detoxify such products directly (Edwards et al. 2000; Hatton et al. 1998; Reade et al. 2004)

Peroxisome-related genes were also identified (Supplementary Tables 2 and 3). In plant cells, peroxisomes are highly dynamic compartments with subcellular movement and distribution that generally depend upon the actin cytoskeleton rather than the microtubules. Peroxisomes are also involved in hormone production, including auxin. These subcellular organelles have an oxidative type of metabolism and are considered one of the major ROS formation sites in plant cells (Pan et al. 2020). Peroxisome-generated ROS are involved in diverse cellular and physiological functions (Sandalio and Romero-Puertas 2015). For instance, the overproduction of ROS in the peroxisomes in response to metabolic or environmental changes was reportedly involved in stress-induced oxidative damage in the plant (Pan et al. 2020). Other studies have shown that peroxisomes’ ROS metabolism acts as an oxidant signaling source that can participate in plant cell metabolism under both physiological and stress conditions (Reumann et al. 2009). Another association with oxidative stress–enhanced response is based on DE genes related to salicylic acid production. Previous studies showed that salicylic acid could alleviate herbicide-induced oxidative stress (Ghahremani et al. 2022; Radwan et al. 2019).

A DE gene encoding an aminocyclopropane-1-carboxylic acid synthase (ACS) also could be related to dicamba resistance. Auxinic herbicides mimic endogenous auxins in plants and are known to upregulate the synthesis of ACS, which carries out the rate-limiting step in the biosynthesis of ethylene via the production of the intermediate 1-aminocyclopropane-1-carboxylic acid. Previous studies in Palmer amaranth (Amaranthus palmeri S. Watson) indicated that dicamba increased ACS expression, leading to ethylene biosynthesis (Johnston et al. 2019, 2020). The constitutive upregulation of this gene in resistant plants may indicate that such plants are predisposed to deal with upregulation of the ethylene biosynthesis pathway.

The potential for detoxification of ROS and oxidative stress alleviation is a clear pattern observed; however, detoxification via the metabolism of dicamba molecules could also be possible. The metabolism of different herbicides, including synthetic auxins, can be specific, as different molecules may require specific enzymes for detoxification or modification (Gaines et al. 2020). In addition, some enzymes may show broad substrate specificity and can metabolize multiple types of molecules (Han et al. 2021; Huang et al. 2021; Mateo-Bonmatí et al. 2021). Direct conjugation of dicamba by GST was not tested. However, studies previously showed that some GSTs in Arabidopsis thaliana show interaction with natural IAA and some synthetic auxins (Smith et al. 2003).

Many DE UGTs, which are known to play a role in the metabolism of both endogenous and synthetic auxins in plants, were also identified (Jin et al. 2013; Mateo-Bonmatí et al. 2021). Even though UDP85A qPCR validation did not corroborate RNA-seq results (Table 2), other UGTs were also identified, indicating that glycosylation of dicamba molecules could be happening. UGTs are a large and diverse family of enzymes that catalyze the transfer of a glycosyl group from a nucleotide diphosphate sugar donor to an acceptor molecule, such as a plant hormone or xenobiotic compound. UGTs are involved in the biosynthesis, storage, and detoxification of various plant metabolites, including flavonoids, alkaloids, and phytohormones (Huang et al. 2018; Jin et al. 2013; Mateo-Bonmatí et al. 2021; Meech et al. 2012; Priest et al. 2006; Šmehilová et al. 2016). Glycosylation plays an important role in the changes in solubility and activity of acceptor molecules, regulation of metabolic homeostasis, and detoxification of xenobiotics. Several UGTs associated with metabolites of herbicides/pesticides have been isolated and characterized for their enzyme activity (Meech et al. 2012).

ABC10 was identified as a potentially key gene both from differential expression and co-expression analyses. This gene could play a role in the compartmentalization of dicamba or its metabolites into vacuoles or the apoplast. ABC transporters are also known to have a function in ABA transport (Do et al. 2018). Synthetic auxins are known to lead to an increase in ABA production, and through this, ROS production leading to plant death (Gaines 2020). ABC transporters also play a role in transport of auxin and other phytohormones, indicating a potential direct effect on the transport of synthetic auxin molecules in the plant (Do et al. 2018). It is important to mention that the same ABC transporter from Family 10 (ABC10) was identified as DE in a subpopulation derived from CHR when conducting a differential expression analysis for 2,4-D and HPPD-inhibitor resistance (Giacomini et al. 2020). Interestingly, some genomic regions where DE genes were identified overlap with some previously identified genomics hot spots in the CHR population. This scenario supports the previously proposed hypothesis by Giacomini et al. (2020) that selection pressure over genomic regions can lead to the selection of cross-resistance. There is still a need to better understand this phenomenon whereby a large region can have multiple genes with altered expression. However, the presence of cis-regulatory elements can play an important role in altering the expression of multiple, proximal genes (Bobadilla and Tranel 2023).

NTS resistance is known to often affect multiple herbicides (Rigon et al. 2020). A commonly shared effect of herbicides is oxidative stress, which will eventually lead to plant death. Evolving a system to deal with ROS overproduction quickly is fundamental for dealing with herbicide stress. Because there was no history of dicamba application on the CHR field population, this scenario indicates a case of cross-resistance wherein the selection pressure of one herbicide led to dicamba resistance. Results suggest that oxidative stress regulation is a primary dicamba resistance mechanism in the CHR population. The observed phenotype correlates well with this hypothesis. Classic synthetic auxin symptoms (e.g., epinasty, leaf curling) are still noticeable after dicamba application, indicating that dicamba molecules can still bind with their respective receptors. However, even with dicamba binding, the further oxidative stress effect that would lead to plant death is not going forward. Many herbicides induce oxidative stress; the hypothesized cross-resistance mechanism selected in CHR could result in reduced sensitivity across all of them.

The putative multiple fronts to reduce dicamba toxicity illustrate the complexity of NTS resistance and its multi-genetic nature, as previously predicted (Bobadilla et al. 2022). Dose–response data from previous work showed a 5-fold dicamba resistance level and a 9-fold 2,4-D resistance in CHR, indicating higher resistance to 2,4-D (Bobadilla et al. 2022; Evans et al. 2019). These previous results, in conjunction with RNA-seq results, are consistent with dicamba resistance in the CHR population resulting from cross-resistance rather than direct selection.

In summary, our results provide the first glimpse into transcriptomics changes on the first reported case of dicamba resistance in A. tuberculatus. Although further validation is needed, our results suggest that an enhanced response to oxidative stress caused by dicamba is reducing dicamba toxicity. Two GST-encoding genes and one peroxidase-encoding gene were identified as the putative major contributors to reducing herbicide oxidative stress. Dicamba metabolism could also happen, as genes such as UGTs were also identified as DE, indicating putative glycosylation of herbicide molecules or their by-products. Genes encoding ABC transporters were also implicated as potential players in transporting metabolites and other molecules to vacuoles for further degradation. Future studies should focus on functional validation of the putative resistance genes, as well as understanding their regulatory nature, and on the contribution of major QTLs to multiple herbicide resistance in the CHR population.