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 Reference Bagavathiannan and Davis2018).
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. Reference Gaines, Duke, Morran, Rigon, Tranel, Küpper and Dayan2020). 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 Reference Ghanizadeh and Harrington2017; Jugulam and Shyam Reference Jugulam and Shyam2019; Rigon et al. Reference Rigon, Gaines, Küpper and Dayan2020).
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 Reference Dopson2022). Recently, soybean cultivars with a dicamba resistance trait were released, increasing the usage of this specific herbicide (Behrens et al. Reference Behrens, Mutlu, Chakraborty, Dumitru, Jiang, LaVallee, Herman, Clemente and Weeks2007). Dicamba is a synthetic auxin herbicide (HRAC Group 4) commercialized since the 1960s (Egan and Mortensen Reference Egan and Mortensen2012), 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 Reference Heap and Duke2018), 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. Reference Hager, Wax, Stoller and Bollero2002; Korres et al. Reference Korres, Norsworthy, Young, Reynolds, Johnson, Conley, Smeda, Mueller, Spaunhorst, Gage, Loux, Kruger and Bagavathiannan2018; Steckel and Sprague Reference Steckel and Sprague2004). 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. Reference Evans, Strom, Riechers, Davis, Tranel and Hager2019; Strom et al. Reference Strom, Hager, Seiter, Davis and Riechers2020).
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. Reference Bobadilla, Giacomini, Hager and Tranel2022). 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. Reference Giacomini, Patterson, Küpper, Beffa, Gaines and Tranel2020). 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 F2 population derived from the CHR population.
Material and Methods
Plant Material and Sequencing
For the RNA-seq experiment, a previously characterized pseudo-F2 population was used (Bobadilla et al. Reference Bobadilla, Giacomini, Hager and Tranel2022). Briefly, CHR plants were treated with 560 g ai ha−1 of dicamba, and the most resistant was selected to cross with a universal sensitive population (WUS) to generate an F1 population, from which a pseudo-F2 line was generated. An F2 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. Reference Giacomini, Gaines, Beffa and Tranel2018). All seeds were subjected to a 50% commercial bleach treatment for 10 min and rinsed twice with water for 10 min each. Pseudo-F2 seeds were suspended in a 0.1 g L−1 agarose solution and placed in a 4 C refrigerator for 4 wk to reduce dormancy (Bell et al. Reference Bell, Hager and Tranel2013).
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-cm3 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-F2 plants were treated with 560 g dicamba ha−1 (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. Reference Bobadilla, Giacomini, Hager and Tranel2022). 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-F2 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. Reference Bobadilla, Giacomini, Hager and Tranel2022). Photos were taken 14 d after treatment with dicamba at 560 g ai ha−1. 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. Reference Simms, Cizdziel and Chomczynski1993) 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. Reference Chen, Zhou, Chen and Gu2018). Reads were further filtered to remove rRNA using the software SortMeRNA (Kopylova et al. Reference Kopylova, Noé and Touzet2012). 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. Reference Montgomery, Giacomini, Waithaka, Lanz, Murphy, Campe, Lerchl, Landes, Gatzmann, Janssen, Antonise, Patterson, Weigel and Tranel2020) using STAR (Dobin and Gingeras Reference Dobin and Gingeras2015), and all bam files were combined using the SAMtools merge tool (Li et al. Reference Li, Handsaker, Wysoker, Fennell, Ruan, Homer, Marth, Abecasis and Durbin2009). 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. Reference Bryant, Johnson, DiTommaso, Tickle, Couger, Payzin-Dogru, Lee, Leigh, Kuo and Davis2017), 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. Reference Simão, Waterhouse, Ioannidis, Kriventseva and Zdobnov2015) 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 Reference Langdon2015) 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 Reference Wu and Watanabe2005). All reads were pseudo-aligned to the assembled transcriptome using Kallisto for read count quantification (Bray et al. Reference Bray, Pimentel, Melsted and Pachter2016).
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. Reference Robinson, McCarthy and Smyth2010). Quantification files generated from Kallisto were loaded into R using the R package Tximport and summarized to the gene level (Soneson et al. Reference Soneson, Love and Robinson2015). 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 Reference Alexa and Rahnenführer2009).
Weighted Co-expression Network Analysis (WGCNA)
A WGCNA was conducted using the WGCNA R package (Langfelder and Horvath Reference Langfelder and Horvath2008). 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 log2 scale (Love et al. Reference Love, Huber and Anders2014). 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 Reference Quinlan and Hall2010). 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. Reference Tian, Yang, Meng, Jin and Gao2020).
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-F2 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 Reference Livak and Schmittgen2001). 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-F2 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−1. 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 F2 population greatly reduced background noise from the data (Giacomini et al. Reference Giacomini, Gaines, Beffa and Tranel2018; 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. Reference Gleason, Foley and Singh2011; Johnston et al. Reference Johnston, Malladi, Vencill, Grey, Culpepper, Henry, Czarnota and Randell2020; Tsuchisaka and Theologis Reference Tsuchisaka and Theologis2004).
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 −log10 false discovery rate (FDR), and the x axis refers to the log2 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. Reference Mateo-Bonmatí, Casanova-Sáez, Šimura and Ljung2021). 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. (Reference Giacomini, Patterson, Küpper, Beffa, Gaines and Tranel2020) 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. Reference Rigon, Gaines, Küpper and Dayan2020; Tani et al. Reference Tani, Chachalis and Travlos2015; Zhang and Yang Reference Zhang and Yang2021). 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. Reference Giacomini, Patterson, Küpper, Beffa, Gaines and Tranel2020).
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. Reference Huang, Zhao, Zhang, Li, Shen, Li and Hou2021; Mateo-Bonmatí et al. Reference Mateo-Bonmatí, Casanova-Sáez, Šimura and Ljung2021; Meech et al. Reference Meech, Miners, Lewis and Mackenzie2012).
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 Reference Mhamdi and Noctor2015; Pan et al. Reference Pan, Liu, Wang and Hu2020). 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 Reference Mhamdi and Noctor2015; Pan et al. Reference Pan, Liu, Wang and Hu2020; Rodríguez-Serrano et al. Reference Rodríguez-Serrano, Pazmiño, Sparkes, Rochetti, Hawes, Romero-Puertas and Sandalio2014). 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 Reference Gaines2020).
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. Reference Wang, Niu and Zheng2021). 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 Reference Martín-Trillo and Cubas2010).
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. Reference Andersen, Algreen-Petersen, Hoedl, Jurkiewicz, Cvitanich, Braunschweig, Schauser, Oh, Twell and Jensen2007; Song et al. Reference Song, Zhang, Wu and Zhang2016; Zhou et al. Reference Zhou, Hu, Ye, Jiang and Liu2018). 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. Reference Theissen, Becker, Di Rosa, Kanno, Kim, Münster, Winter and Saedler2000). 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 Reference Bobadilla and Tranel2023).
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. Reference Giacomini, Gaines, Beffa and Tranel2018, Reference Giacomini, Patterson, Küpper, Beffa, Gaines and Tranel2020). 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. Reference Giacomini, Gaines, Beffa and Tranel2018).
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. Reference Gleason, Foley and Singh2011). Previous studies also show that peroxidase overexpression can lead to resistance to herbicides, such as paraquat, that cause high oxidative stress (Murgia et al. Reference Murgia, Tarantino, Vannini, Bracale, Carravieri and Soave2004).
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. Reference Evans, Strom, Riechers, Davis, Tranel and Hager2019; Strom et al. Reference Strom, Hager, Seiter, Davis and Riechers2020). 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. Reference Christoffoleti, Figueiredo, Peres, Nissen, Gaines, Christoffoleti, Figueiredo, Peres, Nissen and Gaines2015; Gaines Reference Gaines2020).
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. Reference Christoffoleti, Figueiredo, Peres, Nissen, Gaines, Christoffoleti, Figueiredo, Peres, Nissen and Gaines2015; Gaines Reference Gaines2020). 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. Reference Cummins, Cole and Edwards1999; Edwards et al. Reference Edwards, Dixon and Walbot2000; Kawano Reference Kawano2003). 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. Reference Murgia, Tarantino, Vannini, Bracale, Carravieri and Soave2004). 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. Reference Edwards, Dixon and Walbot2000; Hatton et al. Reference Hatton, Cummins, Price, Cole and Edwards1998; Reade et al. Reference Reade, Milner and Cobb2004)
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. Reference Pan, Liu, Wang and Hu2020). Peroxisome-generated ROS are involved in diverse cellular and physiological functions (Sandalio and Romero-Puertas Reference Sandalio and Romero-Puertas2015). 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. Reference Pan, Liu, Wang and Hu2020). 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. Reference Reumann, Quan, Aung, Yang, Manandhar-Shrestha, Holbrook, Linka, Switzenberg, Wilkerson, Weber, Olsen and Hu2009). 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. Reference Ghahremani, Hassannejad, Alizadeh and Eslam2022; Radwan et al. Reference Radwan, Mohamed, Fayez and Abdelrahman2019).
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. Reference Johnston, Vencill, Grey, Culpepper, Henry and Czarnota2019, Reference Johnston, Malladi, Vencill, Grey, Culpepper, Henry, Czarnota and Randell2020). 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. Reference Gaines, Duke, Morran, Rigon, Tranel, Küpper and Dayan2020). In addition, some enzymes may show broad substrate specificity and can metabolize multiple types of molecules (Han et al. Reference Han, Yu, Beffa, González, Maiwald, Wang and Powles2021; Huang et al. Reference Huang, Zhao, Zhang, Li, Shen, Li and Hou2021; Mateo-Bonmatí et al. Reference Mateo-Bonmatí, Casanova-Sáez, Šimura and Ljung2021). 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. Reference Smith, Nourizadeh, Peer, Xu, Bandyopadhyay, Murphy and Goldsbrough2003).
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. Reference Jin, Ma, Han, Wang, Sun, Zhang, Li and Hou2013; Mateo-Bonmatí et al. Reference Mateo-Bonmatí, Casanova-Sáez, Šimura and Ljung2021). 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. Reference Huang, Zhu, Liu, Chen, Li and Hou2018; Jin et al. Reference Jin, Ma, Han, Wang, Sun, Zhang, Li and Hou2013; Mateo-Bonmatí et al. Reference Mateo-Bonmatí, Casanova-Sáez, Šimura and Ljung2021; Meech et al. Reference Meech, Miners, Lewis and Mackenzie2012; Priest et al. Reference Priest, Ambrose, Vaistij, Elias, Higgins, Ross, Abrams and Bowles2006; Šmehilová et al. Reference Šmehilová, Dobrůšková, Novák, Takáč and Galuszka2016). 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. Reference Meech, Miners, Lewis and Mackenzie2012).
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. Reference Do, Martinoia and Lee2018). Synthetic auxins are known to lead to an increase in ABA production, and through this, ROS production leading to plant death (Gaines Reference Gaines2020). 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. Reference Do, Martinoia and Lee2018). 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. Reference Giacomini, Patterson, Küpper, Beffa, Gaines and Tranel2020). 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. (Reference Giacomini, Patterson, Küpper, Beffa, Gaines and Tranel2020) 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 Reference Bobadilla and Tranel2023).
NTS resistance is known to often affect multiple herbicides (Rigon et al. Reference Rigon, Gaines, Küpper and Dayan2020). 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. Reference Bobadilla, Giacomini, Hager and Tranel2022). 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. Reference Bobadilla, Giacomini, Hager and Tranel2022; Evans et al. Reference Evans, Strom, Riechers, Davis, Tranel and Hager2019). 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.
Supplementary material
To view supplementary material for this article, please visit https://doi.org/10.1017/wsc.2023.73
Acknowledgments
This work was partially supported by the USDA National Institute of Food and Agriculture (grant no. 2020-67013-31854).
The authors declare no competing interests.
Open access
