- 1Department of Biochemistry, Molecular Biology, Entomology & Plant Pathology, Mississippi State University, Starkville, MS, United States
- 2Faculty of Science and Engineering, Southern Cross University, East Lismore, NSW, Australia
- 3Institute for Genomics, Biocomputing and Biotechnology, Mississippi State University, Starkville, MS, United States
- 4Department of Plant and Soil Sciences, Mississippi State University, Starkville, MS, United States
The cotton chromosome substitution line, CS-B15sh, exhibits 41% lower injury from 2,4-D when applied at the field recommended rate of 1.12 kg ae ha−1 (1×) than does Texas Marker-1 (TM-1). CS-B15sh was developed in the genetic background of Gossypium hirsutum L. cv TM-1 and has chromosome introgression on the short arm of chromosome 15 from Gossypium barbadense L. cv. Pima 379. In a previous experiment, we observed reduced translocation of [14C]2,4-D outside the treated leaf tissue in CS-B15sh, which contrasted with an increased translocation of the herbicide in the tissues above and below the treated leaf in TM-1. Our results indicate a potential 2,4-D tolerance mechanism in CS-B15sh involving altered movement of 2,4-D. Here, we used RNA sequencing (RNA-seq) to determine the differential expression of genes between 2,4-D-challenged and control plants of the tolerant (CS-B15sh) and susceptible lines (TM-1 and Pima 379). Several components of the 2,4-D/auxin-response pathway—including ubiquitin E3 ligase, PB1|AUX/IAA, ARF transcription factors, and F-box proteins of the SCFTIR1/AFB complex—were upregulated with at least threefold higher expression in TM-1 compared with CS-B15sh, while both Pima 379 and TM-1 showed the same fold change expression for PB1|AUX/IAA mRNA. Some genes associated with herbicide metabolism, including flavin monooxygenase (Gohir.A01G174100) and FAD-linked oxidase (Gohir.D06G002600), exhibited at least a twofold increase in CS-B15sh than in TM-1 (the gene was not expressed in Pima 379), suggesting a potential relationship between the gene’s expression and 2,4-D tolerance. It is interesting to note that glutathione S-transferase was differentially expressed in both CS-B15sh and Pima 379 but not in TM-1, while cytochrome P450 and other genes involved in the oxidation–reduction process were significantly expressed only in CS-B15sh in response to 2,4-D. Gene set enrichment analysis on the union DEGs of the three cotton genotypes revealed the depletion of transcripts involved in photosynthesis and enrichment of transcripts involved in ABA response and signaling.
Introduction
Gossypium hirsutum and Gossypium barbadense are two cultivated species of allotetraploid cotton. While G. hirsutum accounts for more than 90% of Upland cotton production worldwide, G. barbadense is superior in fiber quality, producing extra-long fibers for superior textile products (Saha et al., 2006; Li et al., 2015; Liu et al., 2015). G. hirsutum L. is known to be sensitive to 2,4-dichlorophenoxyacetic acid (2,4-D), and understanding the genetics of herbicide tolerance and identification of specific gene(s) involved at the molecular level is crucial for the development and genetic improvement of herbicide-resistant commercial cotton. A chromosome substitution line CS-B15sh was developed in the genetic background of Texas Marker-1 or TM-1 (G. hirsutum L.) with introgressions on the short arm of chromosome 15 from G. barbadense L. cv. Pima 379 via hybridization, cytogenetic analysis of progeny, and molecular marker selection (Stelly et al., 2005; Saha et al., 2012). In previous experiments, CS-B15sh exhibited reduced herbicide injury compared to TM-1 cotton seedlings when treated with 2,4-D at 1× field rate in greenhouse and field conditions (Perez, 2021). With the complex allotetraploid genomes of the two Gossypium species in hand and recombination and segregation data for genomic regions during the development of chromosome substitution lines, complementation of alleles leading to the identification of good recombinants can be a source of a novel gene for cotton genetic improvement on traits like herbicide tolerance (Saha et al., 2006; Li et al., 2015; Zhang et al., 2015).
2,4-D is a popular synthetic auxin that kills unwanted dicot plants (Schulz and Sehobye, 2016). At very low concentrations, 2,4-D mimics natural auxin in promoting cell division and elongation, while it exhibits herbicidal activity at high concentrations (Grossmann, 2000; Song, 2014). High doses of 2,4-D applied to sensitive dicots result in abnormal growth, premature senescence, and tissue decay (Grossmann, 2010). The herbicide is believed to act at multiple sites once the compound is absorbed by the plant (Gunsolus and Curran, 1999). It was reported that the mechanism of action of 2,4-D involves the activation of the auxin receptor system resulting in upregulation of auxin responses in the plant, increased ethylene production, and upregulation of ABA biosynthesis (Song, 2014; Schulz and Sehobye, 2016). High doses of auxin have been shown to result in chloroplast damage and progressive chlorosis leading to membrane leakage, overproduction of reactive oxygen species (ROS), localized necrosis, and cell death (Grossmann, 2000; Schulz and Sehobye, 2016).
At the molecular level, tolerance to 2,4-D in wild radish has been associated with mutations in ABCB-type auxin transporters in the plasma membrane; the mutated transporters reduced the herbicide translocation rate in wild radish (Goggin et al., 2016). The biosynthesis and signaling by natural auxin, IAA, which promotes normal cell division, elongation, and normal plant growth, is tightly regulated in plants. At low auxin concentrations, Aux/IAA proteins bind auxin-response factors (ARFs), preventing the transcription and expression of auxin-inducible genes (Song, 2014). However, 2,4-D, when present at high concentrations, acts as molecular glue that binds Aux/IAA to the F-box transport inhibitor response 1 (TIR1) protein and mediates proteasome degradation of Aux/IAA proteins (Dharmasiri et al., 2005). The TIR1 protein appears to have a critical role in auxin (2,4-D) signaling (Parry et al., 2009). This leaves free ARF proteins that bind to auxin-response elements (AuxRes) and facilitates transcription of auxin-response genes leading to herbicidal responses (Song, 2014). Interestingly, TIR1 associated with ubiquitin-mediated auxin signaling was downregulated at low 2,4-D concentrations in Arabidopsis (Raghavan et al., 2006). In addition, the transcriptional gene silencing pathway was affected in response to 2,4-D herbicide treatment, which led to increased susceptibility to the herbicide in Arabidopsis (Markus et al., 2021). Other transcriptional responses as a result of herbicide treatment in other species have been associated with nutrient limitation due to the increased expression of genes involved in alternative carbon and nitrogen source metabolism (Teixeira et al., 2006). The herbicidal effects of 2,4-D on ABA biosynthesis and signaling have been reported in Arabidopsis (Raghavan et al., 2005, 2006). In citrus, ABA levels detected by high-performance liquid chromatography–tandem mass spectrometry (HPLC-MS/MS) were significantly increased with post-harvest 2,4-D treatment (Ma et al., 2014). In salt-tolerant rice cultivar, 2,4-D treatment stimulated the synthesis of the stress hormone ABA regulating specific antiporter activities in the cell associated with the efflux and influx of ions in the plasma membrane (Islam et al., 2017). The effects of 2,4-D treatment on the reduction of photosynthetic process have been reported as early as 1950s (Wedding et al., 1953; Williams and Dunn, 1966). In photosynthetic cyanobacterium Nostoc muscorum Meg 1, transcription levels of several photosynthesis-related genes were compromised with the increasing dose of 2,4-D applied (Sachu et al., 2022). In plants, a significantly reduced photosynthesis and stimulated ABA levels were detected in the weed species Erigeron canadensis following 2,4-D herbicidal application (McCauley et al., 2020). The upregulation of NCED, key enzyme in ABA biosynthesis, is said to be the principal step in the synthetic auxin herbicide mode of action (McCauley et al., 2020).
However, their findings also indicate that the increase in ABA levels were independent of the increase in ethylene which is contrary to previous reports on the pathways leading to auxin herbicidal response and symptoms in plants (Grossmann, 2000; Hansen and Grossmann, 2000). However, the general downregulation of transcript abundance of photosynthesis-related genes is proposed to be the result of ABA accumulation which then leads to loss of photosynthetic capacity, deregulation of growth, and plant death (Gaines, 2020).
Several mechanisms have been described for herbicide tolerance in weeds and cultivated crops, including upregulation or downregulation of families of enzymes, such as cytochrome P450s, glutathione S-transferases (GSTs), and glycosyltransferases; all of these families are involved in herbicide degradation and metabolism of 2,4-D into non-toxic forms (Ohkawa et al., 1999; Werck-Reichhart et al., 2000; Busi and Powles, 2017). Herbicide tolerance to 2,4-D in the weed species Raphanus raphanistrum was shown to be conferred by a single dominant locus that depicted a nuclear inheritance pattern (Busi and Powles, 2017). Similarly, a single codominant gene was responsible for 2-4-D tolerance in prickly lettuce (Lactuca serriola L.) (Riar et al., 2011). Reduced translocation and enhanced metabolism were observed in 2,4-D-resistant corn poppy (Papaver rhoeas) populations in Spain (Torra et al., 2017). Rapid metabolism of 2,4-D was also observed in common waterhemp (Amaranthus tuberculatus), showing tolerance to the herbicide (Shyam et al., 2019).
In cotton, 2,4-D tolerance was introduced by Agrobacterium-mediated transfer of a 2,4-D monooxygenase gene, tfdA, from Alcaligenes eutrophus, which degrades the active compound 2,4-D into non-toxic 2,4-dichlorophenol (Bayley et al., 1992). Dow AgroSciences released the first commercial 2,4-D-tolerant cotton varieties employing the Enlist technology which uses a transgene from bacteria encoding an aryloxyalkanoate dioxygenase that efficiently degrades 2,4-D to non-toxic derivatives (Wright et al., 2010). This technology has gained widespread acceptance and adoption by cotton farmers needing a new control technology for weeds that have developed tolerance to glyphosate while maintaining high seed and lint yield. Deployment of Enlist cotton has led to problems in areas where non-Enlist cotton is also grown. Non-Enlist cotton is highly sensitive to 2,4-D, and off-target spray drift injury incidents are increasing with concomitant damage claims (Byrd et al., 2016). To address this problem, alternative sources of 2,4-D tolerance that could be introduced into cultivated cotton without genetic engineering would be highly desirable.
High-throughput sequencing has enabled the rapid and efficient analysis of complex traits, including the genetic basis of tolerance mechanisms and the metabolic pathways that respond to the application of auxinic herbicides, like 2,4-D. Using this approach, specific-resistant allele variants of a cytochrome P450 have been associated recently with non-target site tolerance to 2,4-D in A. tuberculatus (Giacomini et al., 2020). In wild radish (R. raphanistrum), the RNA-seq analysis of herbicide-resistant plants treated with 500 g ae ha−1 2,4-D amine revealed that the tolerance mechanisms included complex and population-specific changes in auxin signaling and elevated plant defense responses (Goggin et al., 2018). RNA-seq analysis has also been used to dissect glufosinate tolerance in Amaranthus palmeri and investigate global transcriptional changes associated with tolerance to the herbicide (Salas-Perez et al., 2018).
This study aims to understand the interaction of cotton with 2,4-D herbicide at the molecular level. We expect to improve our understanding of the genetics and molecular mechanisms of 2,4-D tolerance by inference from the analysis of specific DEGs that were upregulated or downregulated in the 2,4-D-tolerant line, CS-B15sh, compared to the 2,4-D-sensitive parental lines, TM-1 and Pima 379. We anticipate that this information will inform future breeding efforts to improve the tolerance of non-engineered Upland cotton to herbicide spray drift from adjoining agricultural fields.
Materials and methods
Plant materials and 2,4-D treatment
The chromosome substitution line CS-B15sh (31-4), G. hirsutum L. cv. TM-1, and G. barbadense L. cv. Pima 379 were used in the study. Seeds were obtained from the USDA−ARS Crop Science Research Laboratory, Mississippi State, MS, United States. Seedlings were established by sowing cotton seeds in 18-cell Landmark plastic trays (53 cm × 26 cm) (Landmark Plastic, Akron, OH, United States) containing soil (BX PROMIX Growing Medium 10280, BWI Companies, Inc., Nash, TX, United States) amended with 1–2 g of the basal fertilizer Osmocote Plus (The Scotts Company, Marysville, OH, United States). Seedlings (one per cell) were maintained in the greenhouse and watered once a day prior to treatment. The cotton seedlings were treated when plants reached the 4–5 leaf growth stage by spraying 2,4-D (Weedar 64, Nufarm Americas Inc., Alsip, IL, United States) at a rate of 1 lb. acre−1 (1.12 kg ae ha−1) in a controlled spray chamber equipped with a TP8002VS Even Flat Spray Tip (TeeJet®, Spraying System Co., Wheaton, IL, United States) calibrated to deliver 20 gallons acre−1 (GPA) at 40 psi. After herbicide application, the treated plants were allowed to dry for 1–2 h before being transferred back to the greenhouse; the drying period was included to prevent the transfer of volatilized herbicide to control plants (water-sprayed only). Irrigation was resumed on the treated seedlings after 24 h. The non-treated controls (water-sprayed) were planted and grown at the same time as the treated plants and were maintained in the same manner except for the 2,4-D treatment. The experiment was conducted following a completely randomized design with two treatments (2,4-D-treated and non-treated control) and four replications per treatment for each line (TM-1, Pima 379, and CS-B15sh). Herbicide injury was evaluated 21 days after herbicide application using the 0–100% injury scale by Behrens and Lueschen (1979). Leaf tissues (2nd true leaf) were collected from the treated and non-treated cotton seedlings 12 h after treatment (12 HAT) and were immediately frozen in liquid nitrogen and stored at −80°C prior to extraction of total RNA.
RNA extraction, library preparation, and Illumina sequencing
Total RNA was extracted from the collected leaf tissues using a modified hot borate method (Wan and Wilkins, 1994). Double- and single-stranded DNA were removed from the RNA samples using deoxyribonuclease I (RQ1 RNase-Free DNase, 1000u, Promega Corporation, Madison, WI, United States). The RNA samples were further purified using the Qiagen RNeasy® Mini Kit (Catalog No. 74104, Qiagen, Germantown, MD, United States). The quality of the total RNA was assessed by agarose gel electrophoresis, and the concentration was determined using a Nanodrop spectrophotometer (NanoDropC, Thermo Scientific, Waltham, MA, United States).
The RNA samples then were used to construct cDNA libraries using the NEBNext® Ultra II Directional RNA Library Prep Kit for Illumina with Sample Purification Beads (New England Biolabs, Ipswich, MA, United States) following the manufacturer’s instructions. The concentration and quality of cDNA libraries were determined using fluorometry (Qubit™ fluorometer, Invitrogen™ by Thermo Fisher Scientific, Carlsbad, CA, United States) and capillary zone electrophoresis (Agilent 2100 Bioanalyzer, Santa Clara, CA, United States), respectively. The manufacturer’s procedure described in the Agilent DNA 1000 Reagent kit (Agilent Technologies, Santa Clara, CA, United States) was used in checking the quality of cDNA libraries. The 24 RNA-Seq libraries (treated vs. control, three genotypes, four replicates), each with unique barcodes, were pooled and sequenced on three lanes (eight libraries per lane) of an Illumina HiSeq X-Ten (paired-end, 150 bp) (Illumina, San Diego, CA, United States).
Identification and functional analysis of differentially expressed genes
Alignment of the raw reads to the reference G. hirsutum (AD1) TM-1 genome UTX_v2.1 (Chen et al., 2020) and quantification of transcript abundance were performed using Salmon v1.3.0 (Patro et al., 2017). The Salmon aligner which uses a couple forms of Bayesian inference was employed to estimate abundance with the ambiguously/multiple mapped reads (Hu et al., 2021). The R v4.0.2 package tximport (version 1.16.1) was used to convert the transcript-level data into gene-level analysis using the recommended expression normalization method (Soneson et al., 2016; R Core Team, 2020). Genes with low expression levels (average log-transformed counts per million less than 1) were filtered from further analysis (Bourgon et al., 2010). Differential gene expression analysis was conducted using the quasi-likelihood generalized log-linear model available from edgeR (version 3.30.3), a Bioconductor software package1 (Robinson et al., 2010). The response to 2,4-D was tested by pairwise comparison between the treated and non-treated controls of each line (CS-B15sh, TM-1, and Pima 379). Likewise, the differences in the responses of the three genotypes to 2,4-D were tested by pairwise comparison of the responses of each line to one another. The changes in basal gene expression were also tested by pairwise comparison of the non-treated controls of each line with one another. In a similar manner, changes in gene expression among the treated samples of each line were also identified using pairwise comparison. Genes that were significantly (FDR ≤ 0.05) different in response to 2,4-D in the CS-B15sh line vs. the two susceptible lines, TM-1 and Pima 397, and were not significantly different between TM-1 and Pima 397 were considered genes of interest. A composite DEG list (FDR < 0.001) for herbicide-treated vs. non-treated samples in all three cotton lines (Supplementary Table 1) was used in further downstream analyses such as hierarchical clustering (distance metric Pearson’s correlation, average linkage clustering), metabolic pathway/reaction analysis using CottonGen Tools2, MapMan (version 3.5.1R2) analysis (Thimm et al., 2004), and GO enrichment analysis using FGSEA package (release 3.14) with default settings of 1,000 permutations and P-value cutoff P < 0.05 in R Bioconductor3. Pathway analysis using the gene expression data for 59 selected DEGs was carried out using Pathway Tools Omics Dashboard (Paley et al., 2017). Functional annotation of the DEGs was obtained from the TM-1 genome UTX_v2.1 (Chen et al., 2020).
Results
This study determined that CS-B15sh plants treated with the recommended field rate (1.12 kg ae ha−1) of 2,4-D showed 58% herbicide injury on average, while TM-1 and Pima 379 plants exhibited 97 and 98.5% herbicide injury, respectively (Figure 1). Whereas CS-B15sh showed moderate leaf epinasty, stem curling, and necrosis 14 days after treatment, the TM-1 and Pima 379 plants showed severe injury with epinasty, leaf/stem curling, severe and widespread necrotic lesions, and plant death in some seedlings. Non-treated seedlings of all genotypes showed none of the visible injuries associated with 2,4-D treatment (data not shown).
Figure 1. 2,4-D herbicide injury and symptoms observed on CS-B15sh, TM-1, and Pima 379. Clockwise from upper left, (A) histogram of herbicide injury observed 21 days after application of 1× rate (1.12 kg ae ha–1) 2,4-D on 4–5 leaf stage cotton seedlings in the greenhouse; 2,4-D injury observed in CS-B15sh (B), TM-1 (C), and Pima 379 (D) after 2 weeks from herbicide spraying; (E) in another experiment, CS-B15sh and TM-1 sprayed with drift rate (0.05×) of 2,4-D herbicide at 5 weeks from spraying, with TM-1 showing leaf epinasty while CS-B15sh showing normal leaf growth.
RNA sequencing and differential expression
Sequencing of the short-read Illumina libraries generated approximately 60 million fragments (average) per library. Analysis of the cDNA fragment libraries by capillary zone electrophoresis (Bioanalyzer) showed the expected 250–300 bp size range. The fragment sequences mapped to the G. hirsutum (AD1) TM-1 genome UTX_v2.1 at a rate of 61–74% with 43.5 million fragments per library aligned to the reference genome for cotton. A total of 23,821 genes were found in the samples, of which 920 (3.9%), 7,448 (31.3%), and 560 (2.4%) were exclusively expressed (herbicide-treated vs. untreated) in CS-B15sh, TM-1, and Pima 379, respectively (Figure 2). Common DEGs to both CS-B15sh and TM-1 were 11,698 (49.1%), while 161 (0.7%) were shared between CS-B15sh and Pima 379. A total of 289 (1.2%) DEGs were shared between TM-1 and Pima 379 only, while 2,745 (11.5%) exhibited significant differential expression in the three cotton lines used in the study. The total number of DEGs in response to 2,4-D treatment was 15,524, 22,180, and 3,755 in CS-B15sh, TM-1, and Pima 379, respectively.
Figure 2. Distribution of differentially expressed genes (DEGs) (N = 23,821) with FDR P-value < 0.001 either in CS-B15sh, TM-1, or Pima 379 cotton lines in response to 2,4-D herbicide applied at 4–5 leaf stage cotton seedlings in the greenhouse.
The gene set enrichment analysis (2,745 DEGs, intersect) revealed the significant Gene Ontology (GO) terms involved in cotton’s response to 2,4-D, including regulation of transcription, catalytic activity, response to auxin, transferase activity, extracellular region, cellular amino acid metabolic process, plant-type cell wall organization, ion channel inhibitor activity, biosynthetic process, response to oxidative stress, and xenobiotic transmembrane transporter activity (Supplementary Table 2).
The auxin herbicide response and signaling in cotton
The molecular pathways for auxin and 2,4-D response in plants have been reported (Badescu and Napier, 2006; Hao and Yang, 2010; Song, 2014). Of the key molecular factors identified, no information is available in cotton, and knowing as well as understanding the several components involved is important in the genetic manipulation for herbicide response in Upland cotton. Since 2,4-D mimics the natural auxin hormone in plants, we found several genes reported significantly expressed and upregulated in the cotton lines treated with 2,4-D. This includes ubiquitin E3 ligase, F-box domain, PB1 domain|AUX/IAA proteins, and auxin-response factor (Table 1). The ubiquitin E3 ligase genes (RBR family and RNF126-like protein) exhibited twofold higher in TM-1 compared to CS-B15sh. For instance, Gohir.A05G142700 showed a 106-fold change in TM-1 while CS-B15sh showed 53.3. Similarly, Gohir.D04G063800 showed 9.2 in TM-1, while 3.4-fold change was observed in the 2,4-D-treated CS-B15sh plants. Both genes were not significant in Pima 379. The PB1|AUX/IAA domains, Gohir.A08G062000, and Gohir.A07G037300 exhibited similar fold change in both TM-1 and Pima 379, and the level of expression is threefold higher (FC, 40–49) compared with CS-B15sh (FC, 13–16). Gohir.D08G261100 was not significant in Pima 379 but exhibited threefold higher expression in TM-1 (FC, 72) than CS-B15sh (FC, 23). The auxin-response factor, Gohir.A01G112300, was significant in TM-1 with a sevenfold change, while the gene was not significant in CS-B15sh and Pima 379. A gene coding for F-box domain (Gohir.A12G196800) is highly expressed in TM-1 with a 182-fold change, 3× higher than CS-B15sh (FC, 53). The number of transcripts in terms of counts per million (CPM) detected on each gene corresponded with the fold change of expression described above (Figure 3). E3 ubiquitin ligase showed 8× higher transcript level after 12 h in TM-1 while 2× higher in CS-B15sh compared to the baseline transcript level of untreated plants (Figure 3A). The PB1|AUX/IAA domain exhibited the same increase of transcript level (38×) at 12 h for TM-1 and Pima 379, while only a 12× increase was shown in CS-B15sh (Figure 3B). Interestingly, the F-box domain exhibited a sharp rise in transcript level with a 180-fold increase after 12 h of 2,4-D treatment in TM-1. The gene was not significant in Pima 379 and showed 38× higher expression in CS-B15sh (Figure 3D). On the contrary, the auxin-response factor transcript amounts were 6× higher at 12 h in TM-1, while the gene was not significant in CS-B15sh and Pima 379 (Figure 3C). Other DEGs known as hormone-responsive genes were also significantly upregulated, including AP2/ERF domain protein (Gohir.A13G053300) with a 65-fold change in CS-B15sh. At the same time, oxoglutarate/iron-dependent dioxygenase (Gohir.A13G155500) exhibited a 442-fold increase in expression in treated TM-1 seedlings, 7× higher than the transcript level detected in the treated CS-B15sh seedlings (Table 1).
Table 1. Differentially expressed genes of the 2,4-D auxin-response pathway and other hormone-responsive genes in the cotton lines, CS-B15sh, TM-1, and Pima 379 in response to 2,4-Da.
Figure 3. Number of transcripts detected in DEGs associated with the auxin (2,4-D) herbicide response pathway in plants; (A) Gohir.D04G063800, E3 ubiquitin-protein ligase; (B) Gohir.A07G037300, PB1 domain|AUX/IAA; (C) Gohir.A01G112300, auxin-response factor; and (D) Gohir.A12G196800; F-box domain. A number of transcripts for Gohir.D04G063800, Gohir.A01G112300, and Gohir.A12G196800 detected in Pima 379 are not significantly different. Analysis of variance and comparison of means based on Tukey–Kramer (HSD) was applied using JMP 14 statistical package (SAS Institute, NC, United States).
Effects on photosynthesis and abscisic acid metabolism
A total of 193 genes associated with photosynthesis were significant in all three cotton lines treated with 2,4-D (Supplementary Table 3). Of these, 95% were downregulated and composed of DEGs involved in photosystems I and II, light-harvesting complex, chlorophyll biosynthesis, ATP synthesis, and electron transport chain. In Pima 379, the top five significantly downregulated DEGs were involved in PS I light-harvesting complex, chlorophyll biosynthesis, and ATP synthesis, while TM-1 also included PS I light-harvesting and PS II DEGs. CS-B15sh also exhibited the most significant downregulation of PS I light-harvesting followed by DEGs for the electron transport chain. From the major group of downregulated genes, 2Fe-2S ferredoxin-type iron-sulfur binding domain (Gohir.A13G222600), which is involved in the electro-transport chain, is upregulated in all the herbicide-treated cotton lines. At the same time, Mog1/PsbP of PS II, an extrinsic membrane component, showed upregulation of gene expression in Pima 379 but was downregulated in both CS-B15sh and TM-1. The small group of DEGs for photosynthesis were composed of upregulated genes involved in electron transfer activities, including ATP synthase (Gohir.D01G067400 and Gohir.A01G080200), 2Fe-2S ferredoxin-type iron-sulfur binding domain (Gohir.D13G227400), and Ferredoxin–NADP reductase (Gohir.A05G110400), with two DEGs of unknown function from chromosome 19 (D05) of G. hirsutum L., Gohir.D05G110800 and Gohir.D05G0837000. Photosystem I reaction center subunit IX (PsaJ) superfamily is upregulated in both Pima 379 and TM-1 while exhibiting downregulation of gene expression in CS-B15sh. An overview of metabolic pathways and gene regulatory networks observed in the three cotton lines is presented in Figure 4 and Supplementary Figures 1, 2. Generally, almost all genes found and expressed in the cotton lines and involved in the Calvin cycle, photorespiration, PS I, and PS II exhibited negative log2 fold changes (red boxes) in the three genotypes. More downregulation of photosynthesis-related genes was observed in TM-1 compared to CS-B15sh and Pima 379 (Supplementary Figure 3).
Figure 4. MapMan analysis showing the overview of metabolic pathways represented by genes differentially expressed (2,4-D-treated vs. untreated) in the cotton line CS-B15sh at 12 h after treatment with 1× rate of herbicide 2,4-D. Each data point (square) represents a DEG gene. Blue and red colors indicate the DEGs are upregulated and downregulated, respectively.
The effects of 2,4-D herbicide on ABA signaling and biosynthesis in cotton showed 11 of 18 DEGs were upregulated, exhibiting positive log2 fold change and involving genes for abscisic acid binding and receptor activity (Supplementary Table 4). The top five highly expressed ABA-associated genes encoded mainly Bet v I/Major latex protein (START-like domain superfamily) and were significantly upregulated in both TM-1 and CS-B15sh, except for Gohir.A10G169466, which showed the highest 7.7 log2 fold change only found in TM-1. Aspartic peptidase A1 family (Gohir.A09G177900) was significantly upregulated (FDR < 0.001) in the three cotton lines evaluated with log2 fold change of 1.7, 2.1, and 1.6 for Pima 379, CS-B15sh, and TM-1, respectively (Supplementary Table 4). The rest of the ABA genes found were downregulated including nodulin-related protein 1, dual-specificity phosphatase, protein SHORT HYPOCOTYL IN WHITE LIGHT 1, aspartic peptidase A1 family, and forkhead-associated (FHA) domain. The patterns of gene expressions of DEGs associated with ABA metabolism and response are similar in the cotton lines.
DEGs responsive to 2,4-D in CS-B15sh
Based on the analysis of the different responses of each cotton line to 2,4-D treatment, 59 genes passed the response comparison filter with significantly different (FDR ≤ 0.05) responses to 2,4-D treatment in CS-B15sh compared to susceptible lines but not significantly different (FDR > 0.05) between the susceptible lines (Table 2). Functional annotation revealed that the genes were involved in protein and DNA binding, transmembrane transporter activity, transferase activity, hydrolase activity, photosystem 1, iron ion binding, and nucleosome and transcription processes. Fifty-one genes were significantly (FDR ≤ 0.05) differentially expressed in the CS-B15sh response to 2,4-D, with 18 upregulated and 33 downregulated (Table 2). The top three upregulated DEGs in the CS-B15sh response to 2,4-D are expansin, cellulose-binding-like domain (Gohir.A05G184100), zinc finger, FYVE/PHD-type (Gohir.A07G147600), and a novel cotton gene Gohir.A07G188800 of unknown function with log2 fold change of 4.4, 2.4, and 2.9, respectively. The top five downregulated genes in the CS-B115sh response to 2,4-D were glycoside hydrolase family 16, photosystem I PsaJ, Cytochrome P450 E-class, NDRG|alpha/beta hydrolase fold, and glycosyl transferase family 8 with log2 fold change ranging from −5.4 to −2.6 (Table 2). Interestingly, eight genes were not significant (FDR > 0.05) in the CS-B15sh response to 2-4,D but had a significantly (FDR ≤ 0.05) different responses to 2,4-D compared with the susceptible lines, including nucleotide-diphospho-sugar transferases (Gohir.A06G026500), Myc-type, basic helix–loop–helix (bHLH) domain (Gohir.A13G035300), Ctr copper transporter (Gohir.A03G046900), Cytochrome P450, E-class, group I (Gohir.D06G184100 and Gohir.A12G201300), Heat shock factor binding 1 (Gohir.D11G120200), and two genes of unknown function (Gohir.D12G144700 and Gohir.D06G053801). Four genes, Ctr copper transporter, nucleotide-diphospho-sugar transferases, Myc-type, basic helix–loop–helix (bHLH) domain, and Heat shock factor binding 1, showed an inverted response to 2,4-D in the tolerant compared to the susceptible cotton lines. Interestingly, 16 DEGs were novel cotton transcripts with no functional annotation available (Table 2). Two of these genes, Gohir.A07G188800 and Gohir.D05G048200, were upregulated in response to 2,4-D with eight and threefold change, respectively, while the rest were downregulated in the CS-B15sh response to 2,4-D. Gene co-expression analysis for the 59 DEGs using hierarchical clustering with TM-1 and Pima 379 is presented in Figure 5. GO analysis showed enrichment in genes in pathways involving protein dimerization activity, protein binding, DNA and zinc ion binding, regulation of transcription, hydrolase activity, apoplast, cell wall biogenesis, and organization, as well as photosynthesis (Table 3). Further pathway and network analysis using ptools (Paley et al., 2017) revealed that five genes, including Gohir.A12G201300, Gohir.D08G169900, Gohir.D05G079500, Gohir.A13G049200, and Gohir.D11G154700, are involved in ABA degradation, hormone biosynthesis, L-arginine degradation, and protein modification (Table 4). Both Gohir.A12G201300 and Gohir.D08G169900 are cytochrome P450 with iron ion and heme binding functions and oxidoreductase activity associated with the degradation of ABA into phaseic acid (Figure 6). Gohir.D05G079500 has no known function yet, but is predicted to be involved in ABA biosynthesis (Supplementary Figure 4), while Gohir.A13G049200 and Gohir.D11G154700 are involved in the arginine monooxygenase pathway and protein ubiquitination, respectively (Supplementary Figures 5, 6).
Table 2. Listing of cotton genes with significantly different responses to 2,4-D in CS-B15sh compared with TM-1 and Pima 379 at 12 h after application of 1.12 kg ae ha−1 herbicide1.
Figure 5. Gene co-expression heatmap with hierarchical clustering of 59 genes differentially expressed (12 HAT treated vs. untreated) in the three cotton lines, CS-B15sh, TM-1, and Pima 379, sprayed with 1× rate 2,4-D herbicide. DEGs with “–” indicate unknown function.
Table 3. Gene ontology (GO) enrichment analysis results for 59 genes responsive to 2,4-D treatment in the herbicide-tolerant cotton line CS-B15sh.
Table 4. Ptools pathway analysis results for 59 selected DEGs responsive to 2,4-D treatment in the herbicide-tolerant cotton line CS-B15sh showing the identified pathways and reactions for five DEGs.
Figure 6. Ptools pathway analysis shows that two of the 59 DEGs responsive to 2,4-D treatment in the herbicide-tolerant cotton line CS-B15sh (Gohir.D08G169900 and Gohir.A12G201300) are involved in the 2-ci-abscisate → 8′-hydroxyabscisate conversion step in the abscisic acid degradation pathway. The pathway collage generated using the Pathway Tools Omics Dashboard (https://ptools.cottongen.org).
Of the 59 genes that passed the response comparison filter, 19 had significantly (FDR ≤ 0.05) different expressions between the non-treated samples of the susceptible cotton lines (Supplementary Table 5). Among them were two genes with protein kinase domain (Gohir.A05G288000, Gohir.A09G213644), zinc finger CCCH-type (Gohir.D07G023700), zinc finger RING-type (Gohir.D11G154700, Gohir.D09G182000), JmjC/JmjN domain|zinc finger (Gohir.A03G000035), leucine-rich repeat F-box-like domain (Gohir.A03G025000 and Gohir.A06G098214), amidase signature superfamily (Gohir.A13G049200), NRAMP family (Gohir.D03G026100), protein BRANCHLESS TRICHOME-like (Gohir.A11G093500), HEAT repeat|Armadillo-type fold (Gohir.D06G110600), Cytochrome P450 E-class group 1 (Gohir.D06G184100), homeobox|ELK domain (Gohir.A07G030300), DUF1084 (Gohir.A05G168800), COBRA (Gohir.A07G180900), and three genes of unknown function (Gohir.D12G144700, Gohir.A06G049200, and Gohir.D02G083700).
From the composite DEGs list, some genes were significant (FDR ≤ 0.001) and differentially expressed only in CS-B15sh in response to 2,4-D (treated vs. non-treated) but are not significantly expressed (FDR > 0.001) in TM-1 and Pima 379 plants. These DEGs were not found in the non-treated control pairwise comparisons between CS-B15sh vs. TM-1, CS-B15sh vs. Pima 379, and TM-1 vs. Pima 379. A total of 27 DEGs significantly upregulated in CS-B15sh with 8- to 133-fold change expression were detected (Supplementary Table 6). Functional annotations revealed genes primarily involved in oxidation–reduction and metabolic processes demonstrating that herbicide metabolism is associated with the reduced herbicide injury symptoms observed in CS-B15sh. Several genes involved in transmembrane transporter activities were highly upregulated and explicitly found in CS-B15sh, including Gohir.A05G02220 (SLC26A/SulP transporter domain) and Gohir.A11G129600 (EamA domain|WAT1-related protein) with 133- and 46-fold change, respectively. Other genes involved in plant-type cell wall organizations include the expansin, cellulose-binding-like domain, Gohir.A13G076500, which exhibited 57-fold change of gene expression. Interestingly, several genes involved in oxidation–reduction process (Gohir.D08G249300, Gohir.A06G152100, Gohir.D07G226200, and Gohir.A05G173900) are significant and differentially expressed in CS-B15sh in response to 2,4-D but not in the susceptible lines nor the non-treated control pairwise comparisons.
Herbicide degradation and metabolism
Flavin monooxygenase (FMO) (Gohir.A01G174100) associated with herbicide metabolism is differentially expressed in CS-B15sh. The gene showed 6.8 log2 FC, which is twofold higher than TM-1 but is not differentially expressed in Pima 379 (Supplementary Table 7). After 12 h, the transcript level is 5× higher in CSB15sh than TM-1 (Figure 7). Flavin adenine dinucleotide (FAD)-linked oxidase (Gohir.D06G002600) also showed a twofold higher expression in CS-B15sh with 5.7 log2 FC compared with the susceptible cotton line TM-1. However, the transcript expression profile showed a similar pattern for CS-B15sh and TM-1 at 43–49 CPM. Other genes associated with herbicide degradation in plants were also highly expressed in CS-B15sh, including amine oxidase (Gohir.A01G025200), glutathione S-transferase (GST) (Gohir.A11G195400 and Gohir.D11G232100), and cytochrome P450 superfamily (Gohir.D11G187200). Log2 fold change of amine oxidase is 2.7 and 1.9 for CS-B15sh and TM-1, respectively (Supplementary Table 7). Both GSTs were significantly expressed in CS-B15sh and Pima 379 with 1.1–1.4 log2 fold change, while the DEGs were not significant in TM-1. Cytochrome P450 showed 1.6 log2 fold change after 2,4-D treatment, while it was not significant in both TM-1 and Pima 379. The transcript levels detected for amine oxidase at 12 h after 2,4-D treatment were 333 and 271 CPM in CS-B15sh and TM-1, respectively. Both DEGs for GST exhibited twofold higher transcript expression in CS-B15sh than TM-1 (Figure 7). Gohir.A11G195400 and Gohir.D11G232100 exhibited 120 and 198 CPM, respectively, in CS-B15sh, while TM-1 exhibited 50% of the CPM detected in both genes. Other genes involved in the metabolism of endogenous substrates including xenobiotics include glycoside hydrolase (Gohir.D09G116300) whose expression was 130-fold higher in treated versus control CS-B15sh plants, which was 6.4 times the change seen in treated versus control TM-1 plants (Supplementary Table 7).
Figure 7. Transcript level expression of DEGs associated with herbicide degradation and metabolism, including flavin monooxygenase, Gohir.A01G174100 (A); FAD-linked oxidase, Gohir.D06G002600 (B); amine oxidase, Gohir.A01G025200 (C); glutathione S-transferase, Gohir.A11G195400 and Gohir.D11G232100 (D,E); and Cytochrome P450 monooxygenase, Gohir.D11G187200 (F).
Discussion
This paper highlights an attempt to understand the interactions of Upland cotton with 2,4-D herbicide at the molecular level. We selected 12 HAT time points for analysis of cotton’s response to 2,4-D based on studies of gene expression responses to abiotic stress in other plants like Arabidopsis and Brassica napus (Ishitani et al., 1998; Chen et al., 2011). The cotton chromosome substitution line CS-B15sh carries substitution on the short arm of chromosome 15 (D01) from G. barbadense L. In previous experiments, a selection of this line (CS-B15sh) showed reduced 2,4-D injury under greenhouse conditions, while TM-1, the genetic background of CS-B15sh, showed sensitivity to 2,4-D and exhibited severe herbicide injury symptoms and plant death. Apparently, Pima 379 is also sensitive to 2,4-D, and thus, the herbicide tolerance in CS-B15sh could be the result of genetic interactions, complementation of alleles of G. hirsutum and G. barbadense in the development of chromosome substitution line, or formation of transgressive segregants that contributed to phenotypes not observed in the parents such as the level of 2,4-D tolerance observed in CS-B15sh (deVicente and Tanksley, 1993; Rieseberg et al., 1999; Saha et al., 2006). The molecular pathway involved in plant response to 2,4-D has been outlined in previous reports (Grossmann, 2010; Song, 2014; Gaines et al., 2020). Since Upland cotton is sensitive to 2,4-D, it should be interesting to see these molecular key players in the auxin-response pathway differentially expressed with 2,4-D treatment. Several elements of the auxin-signaling pathway outlined in those previous studies were differentially expressed in this study but with varying fold changes observed between the 2,4-D-treated CS-B15sh and TM-1 plants. For instance, the expression of the ubiquitin E3 ligase, part of the ubiquitin-proteasome pathway that degrades AUX/IAA repressor proteins under high auxin/2,4-D conditions, was higher in TM-1 than in CS-B15sh. When AUX/IAAs are degraded, auxin-response factors (ARF) are released, activating auxin-responsive genes and elevating expressions in the herbicide-treated plants (Mockaitis and Estelle, 2008; Guilfoyle, 2015). However, both Gohir.A05G142700 and Gohir.D04G063800 were not differentially expressed in Pima 379 (herbicide 2,4-D-treated vs. untreated) signaling that these genes are not responsive to 2,4-D. The molecular intricacies surrounding how the cotton genome responds to herbicidal 2,4-D and what leads to sensitivity in the plant are still largely unknown. It is possible that other ubiquitin ligase protein families are involved in Pima 379 which is different from those expressed in CS-B15sh and TM-1. Interestingly, three DEGs (Gohir.D08G261100, Gohir.A08G062000, and Gohir.A07G037300) are characterized with PB1 domains showing threefold higher expression in TM-1 compared to CS-B15sh. Proteins with the PB1 domain have been shown to engage in protein–protein interactions with ARF transcription factors (Guilfoyle, 2015) to activate auxin-responsive genes, which leads to severe herbicide injury symptoms as those observed in treated TM-1 plants. This would be consistent with the observation that the expression of an auxin-response factor (Gohir.A01G112300) was threefold higher in treated TM-1 plants than CS-B15sh. In addition, the expression of a gene encoding an F-box (Gohir.A12G196800) is associated with auxin signaling through the promotion of interaction between AUX/IAA proteins and the SCFTIR1/AFB was also higher in TM-1. The SCFTIR1/AFB complex facilitates ubiquitination and subsequent degradation of AUX/IAA repressor proteins, leading to the activation of auxin-responsive genes and disarray of physiological responses caused by 2,4-D at herbicidal rates (Dharmasiri et al., 2005; Song, 2014; Takahashi et al., 2017). The activation and increased expression of hormone-responsive genes such as AP2/ERF proteins and an oxoglutarate/iron-dependent dioxygenase involved in ethylene signaling and biosynthesis indicate that these genes are responsive to auxin herbicide application (Yang and Hoffman, 1984; Houben and Van de Poel, 2019). Future studies in this area are suggested to focus on quantitative gene expression studies of these key elements in the auxin-response pathway to better understand how the auxin/2,4-D response pathway is differentially modulated in CS-B15sh and TM-1 plants exposed to 2,4-D. The effects of auxin herbicides on photosynthesis and ABA metabolism in cotton were similar to recent studies reported (Gaines, 2020; McCauley et al., 2020). Genes associated with photosynthesis were severely downregulated in the 2,4-D-treated cotton plants. This confirms the findings on auxin responses in the weed species Erigeron canadensis treated with dicamba and 2,4-D (McCauley et al., 2020). At 12 HAT, genes involved in PS I, PS II, light-harvesting complex, electron transport chain, and chlorophyll biosynthesis were downregulated, indicating the whole scale shut down of photosynthetic process in cotton. Interestingly, both Gossypium species (G. hirsutum L. and G. barbadense L.) showed similar patterns of downregulation of genes involved in photosynthesis, and the differential expression observed is due to genetic background and is not associated with herbicide tolerance in Upland cotton. Although a few genes were found upregulated, the effects of these genes could translate into differences in the coping mechanisms of cotton to 2,4-D and survival from herbicide injury. Similar observations were obtained on ABA and the effects of 2,4-D on cotton, as previously reported in other plants (Gaines, 2020; McCauley et al., 2020). An increase in ethylene and ABA production is well known to be associated with auxin herbicide injury symptoms, such as leaf epinasty, tissue swelling, growth inhibition, tissue decay, and senescence (Grossmann, 2010). Our transcriptomic data also suggested the increased gene expression of DEGs involved in ABA signaling and biosynthesis in Upland cotton treated with herbicidal 2,4-D. Although the ABA biosynthetic gene, 9-cis-epoxycarotenoid dioxygenase (NCED), was not detected in the analysis, several genes particularly associated with zeaxanthin synthesis are downregulated. In contrast, genes involved in the synthesis and conversion of intermediates to the formation of ABA from 9′-cis-Violaxanthin to xanthoxin and abscisic acid aldehyde exhibited positive log2 fold change. Aspartic peptidase was also reported affecting ABA accumulation in Arabidopsis, while the gene (Gohir.A09G177900) was also found significantly upregulated in all three cotton lines in this study (Kalladan et al., 2017). However, our transcriptome data do not provide evidence on the association of ABA signaling and biosynthesis on herbicide tolerance mechanisms in Upland cotton.
The gene enrichment and pathway analysis of the 59 DEGs responsive to 2,4-D revealed genes associated with photosynthesis and ABA, the five DEGs, including Gohir.A12G201300, Gohir.D08G169900, Gohir.D05G079500, Gohir.A13G049200, and Gohir.D11G154700, present possible pathways that may be involved in herbicide treatment and tolerance response and deserve further investigation. It is interesting to note that Gohir.D05G079500 predicted to be involved in ABA biosynthesis is downregulated in CS-B15sh revealing a relationship with the stress hormone accumulation and possible reduction of herbicide injury symptoms observed in the herbicide-tolerant cotton genotype. While a myriad of molecular patterns and interactions are going on due to 2,4-D’s action, it is worth noting that CS-B15sh exhibits a tolerance level to the herbicide. Comparison of the response of CS-B15sh to those of the two susceptible lines yielded genes involved in cell wall organization and cellular transport mechanisms, indicating possible roles in the diffusion and cell-to-cell movement of 2,4-D active compounds and the tolerance mechanisms exhibited by CS-B15sh to the herbicide. However, some of the genes that showed a significant difference in response to the herbicide were also found significantly expressed in the non-treated plants of TM-1 and Pima 379, indicating that these genes are not associated with defense response to 2,4-D in CS-B15sh. Among these genes are protein kinase-like domain, homeobox domain, zinc finger (CCH-, RING-, and C5HC2-type) protein BRANCHLESS TRICHOME-like, amidase signature (AS) superfamily, NRAMP family, leucine-rich repeat F-box-like domain superfamily, HEAT repeat, COBRA plant, and cytochrome P450 E-class group I (Pace and Brenner, 2001; Roudier et al., 2005; Nevo and Nelson, 2006; Guo et al., 2009). Interestingly, pathway and GO analysis of the DEGs responsive to 2,4-D in CS-B15sh revealed protein dimerization and binding activities, regulation of transcription, hydrolase activity, cell wall, and carbohydrate metabolic processes are among the most enriched pathway signaling regulation of proteins such as enzymes, co-factors, ion channels, and transcription factors (Marianayagam et al., 2004). Both protein containing assembly and hydrolase activity were also associated with the non-target site resistance phenotypes to herbicide phenylurea chlorotoluron and aryloxyphenoxypropionate fenoxaprop acid in black-grass weeds (Franco-Ortega et al., 2021). Glycoside hydrolase is primarily involved in the processing of carbohydrates and works on the metabolism of endogenous substrates, including xenobiotics (Oesch-Bartlomowicz and Oesch, 2007). Other genes involved in xenobiotic metabolism like glycosyltransferases were expressed only in CS-B15sh (Brazier-Hicks et al., 2007). The downregulation of these genes after 2,4-D treatment signals a different molecular pathways contributing to the herbicide tolerance exhibited by CS-B15sh. The three cotton lines we used exhibited a differential response to 2,4-D, and CS-B15sh showed some degree of tolerance to the herbicide compared to TM-1 and Pima 379. Likely, some physiological and metabolic processing of the herbicide compound within the plant occurs with the high transcript expression patterns of DEGs associated with herbicide metabolism, such as those involved in oxidation–reduction processes, including GSTs, amine oxidase, and FAD-linked oxidase. Although recent reports have implicated L-lectin domain-containing receptor kinases, ABC transporters, and cytochrome P-450s in the tolerance of various weeds to 2,4-D (Giacomini et al., 2020; Goggin et al., 2020). We found examples of these genes upregulated in 2,4-D-treated CS-B15sh and TM-1 plants. However, there was no significant increase in the expression of any of these genes that will contribute to herbicide tolerance in our study. Future experiments should focus on elucidating the genetic functions of DEGs highly expressed in CS-B15sh, such as glutathione S-transferase and flavin monooxygenase, because these genes are involved in herbicide detoxification. It would be interesting to do further assays to validate its functions in terms of herbicide tolerance in cotton, develop genetic markers useable as molecular markers for 2,4-D tolerance in marker-assisted breeding, and therefore contribute to the genetic improvement of modern Upland cotton varieties.
Conclusion
In summary, the interaction of cotton to 2,4-D herbicide at the molecular level was elucidated in this study. Several genes involved in the auxin (2,4-D) response pathway previously reported have been detected. Downregulation of a whole suite of genes involved in photosynthesis was observed. The response of the 2,4-D-tolerant line CS-B15sh showed upregulation of genes involved in the oxidation–reduction process. It will be interesting to conduct a follow-up investigation on the specific function of these genes in cotton’s response to 2,4-D. This information will allow the design of appropriate breeding strategy toward the development of modern Upland cotton varieties with improved tolerance to herbicides.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA824818.
Author contributions
TT, LP, ZM, JD, and DP designed the experiments. LP conducted the experiments and wrote the original version of this manuscript. RM and MA performed data analysis. All authors have read, contributed, and approved the manuscript.
Funding
The funding for this project was provided by the Cotton Incorporated and the Mississippi Agricultural and Forestry Experiment Station. It is based upon work supported by the National Institute of Food and Agriculture, United States Department of Agriculture, Hatch project under accession number: 230100.
Acknowledgments
We acknowledge Chuan-Yu Hsu for her technical expertise and assistance in RNA isolation from cotton, library prep, and inputs in design of the experiment and Sukumar Saha and Evita Lopez for providing the seeds of cotton lines used in the study.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2022.910369/full#supplementary-material
Supplementary Figure 1 | MapMan analysis showing the overview of metabolic pathways in TM-1. Each data point (square) represents a DEG gene. Blue and red colors indicate the DEGs are upregulated and downregulated, respectively.
Supplementary Figure 2 | MapMan analysis showing the overview of metabolic pathways in Pima 379. Each data point (square) represents a DEG gene. Blue and red colors indicate the DEGs are upregulated and downregulated, respectively.
Supplementary Figure 3 | Gene co-expression heatmap with hierarchical clustering of DEGs associated with photosynthesis.
Supplementary Figure 4 | Ptools pathway analysis on hormone biosynthesis (Gohir.D05G079500).
Supplementary Figure 5 | Ptools pathway analysis on L-arg degradation (Gohir.A13G049200).
Supplementary Figure 6 | Ptools pathway analysis on protein modification (Gohir.D11G154700).
Supplementary Table 1 | List of DEGs in the composite set (union of 3 genotypes).
Supplementary Table 2 | Gene ontology.
Supplementary Table 3 | List of DEGs associated with photosynthesis.
Supplementary Table 4 | List of DEGs associated with ABA signaling and biosynthesis.
Supplementary Table 5 | Gene expression in UTD (untreated) and 12 HAT contrasts of the 59 DEGs.
Supplementary Table 6 | List of DEGs unique and significantly expressed (FDR < 0.0001) in CS-B15sh (treated vs. untreated).
Supplementary Table 7 | DEGs associated with herbicide metabolism.
Supplementary Table 8 | Read counts and mapping rate.
Footnotes
- ^ http://www.bioconductor.org
- ^ http://ptools.cottongen.org/
- ^ https://bioconductor.org/packages/release/bioc/html/fgsea.html
References
Badescu, G. O., and Napier, R. M. (2006). Receptors for auxin: Will it all end in TIRs? Trends Plant Sci. 11, 217–223. doi: 10.1016/j.tplants.2006.03.001
Bayley, C., Trolinder, N., Ray, C., Morgan, M., Quisenberry, J. E., and Ow, D. W. (1992). Engineering 2,4-D resistance into cotton. Theor. Appl. Genet. 83, 645–649. doi: 10.1007/BF00226910
Behrens, R., and Lueschen, W. E. (1979). Dicamba volatility. Weed Sci. 27, 486–493. doi: 10.1017/S0043174500044453
Bourgon, R., Gentleman, R., and Huber, W. (2010). Independent filtering increases detection power for high-throughput experiments. Proc. Natl. Acad. Sci. U. S. A. 107, 9546–9551. doi: 10.1073/pnas.0914005107
Brazier-Hicks, M., Offen, W. A., Gershater, M. C., Revett, T. J., Lim, E. K., Bowles, D. J., et al. (2007). Characterization and engineering of the bifunctional N- and O-glucosyltransferase involved in xenobiotic metabolism in plants. Proc. Natl. Acad. Sci. U. S. A. 104, 20238–20243. doi: 10.1073/pnas.0706421104
Busi, R., and Powles, S. B. (2017). Inheritance of 2,4-D resistance traits in multiple herbicide-resistant Raphanus raphanistrum populations. Plant Sci. 257, 1–8. doi: 10.1016/j.plantsci.2017.01.003
Byrd, S. A., Collins, G. D., Culpepper, A. S., Dodds, D. M., Edmisten, K. L., Wright, D. L., et al. (2016). Cotton stage of growth determines sensitivity to 2,4-D. Weed Technol. 30, 601–610. doi: 10.1614/WT-D-15-00191.1
Chen, L., Zhong, H., Ren, F., Guo, Q. Q., Hu, X. P., and Li, X. B. (2011). A novel cold-regulated gene, COR25, of Brassica napus is involved in plant response and tolerance to cold stress. Plant Cell Rep. 30, 463–471. doi: 10.1007/s00299-010-0952-3
Chen, Z. J., Sreedasyam, A., Ando, A., Song, Q., De Santiago, L. M., Hulse-Kemp, A. M., et al. (2020). Genomic diversifications of five allopolyploid Gossypium species and their impact on cotton improvement. Nat. Genet. 52, 525–533. doi: 10.1038/s41588-020-0614-5
deVicente, M. C., and Tanksley, S. D. (1993). QTL analysis of transgressive segregation in an interspecific tomato cross. Genetics 134, 585–596. doi: 10.1093/genetics/134.2.585
Dharmasiri, N., Dharmasiri, S., and Estelle, M. (2005). The F-box protein TIR1 is an auxin receptor. Nature 435, 441–445. doi: 10.1038/nature03543
Franco-Ortega, S., Goldberg-Cavalleri, A., Walker, A., Brazier-Hicks, M., Onkokesung, N., and Edwards, R. (2021). Non-target site herbicide resistance is conferred by two distinct mechanisms in black-grass (Alopecurus myosuroides). Front. Plant Sci. 12:636652. doi: 10.3389/fpls.2021.636652
Gaines, T. A. (2020). The quick and the dead: A new model for the essential role of ABA accumulation in synthetic auxin herbicide mode of action. J. Exp. Bot. 71, 3383–3385. doi: 10.1093/jxb/eraa178
Gaines, T. A., Duke, S. O., Morran, S., Rigon, C. A. G., Tranel, P. J., Küpper, A., et al. (2020). Mechanisms of evolved herbicide resistance. J. Biol. Chem. 295, 10307–10330. doi: 10.1074/jbc.REV120.013572
Giacomini, D. A., Patterson, E. L., Kupper, A., Beffa, R., Gaines, T. A., and Tranel, P. J. (2020). Coexpression clusters and allele-specific expression in metabolism-based herbicide resistance. Genome Biol. Evol. 12, 2267–2278. doi: 10.1093/gbe/evaa191
Goggin, D. E., Bringans, S., Ito, J., and Powles, S. B. (2020). Plasma membrane receptor-like kinases and transporters are associated with 2,4-D resistance in wild radish. Ann. Bot. 125, 821–832. doi: 10.1093/aob/mcz173
Goggin, D. E., Cawthray, G. R., and Powles, S. B. (2016). 2,4-D resistance in wild radish: Reduced herbicide translocation via inhibition of cellular transport. J. Exp. Bot. 67, 3223–3235. doi: 10.1093/jxb/erw120
Goggin, D. E., Kaur, P., Owen, M. J., and Powles, S. B. (2018). 2,4-D and dicamba resistance mechanisms in wild radish: Subtle, complex and population specific? Ann. Bot. 122, 627–640. doi: 10.1093/aob/mcy097
Grossmann, K. (2000). Mode of action of auxin herbicides: A new ending to a long, drawn out story. Trends Plant Sci. 5, 506–508. doi: 10.1016/S1360-1385(00)01791-X
Grossmann, K. (2010). Auxin herbicides: Current status of mechanism and mode of action. Pest Manage. Sci. 66, 113–120. doi: 10.1002/ps.1860
Guilfoyle, T. J. (2015). The PB1 domain in auxin response factor and Aux/IAA proteins: A versatile protein interaction module in the auxin response. Plant Cell 27, 33–43. doi: 10.1105/tpc.114.132753
Gunsolus, J. L., and Curran, W. S. (1999). Herbicide mode of action and injury symptoms. Minnesota: University of Minnesota Extension Service.
Guo, Y. H., Yu, Y. P., Wang, D., Wu, C. A., Yang, G. D., Huang, J. G., et al. (2009). GhZFP1, a novel CCCH-type zinc finger protein from cotton, enhances salt stress tolerance and fungal disease resistance in transgenic tobacco by interacting with GZIRD21A and GZIPR5. New Phytol. 183, 62–75. doi: 10.1111/j.1469-8137.2009.02838.x
Hansen, H., and Grossmann, K. (2000). Auxin-induced ethylene triggers abscisic acid biosynthesis and growth inhibition. Plant Physiol. 124, 1437–1448. doi: 10.1104/pp.124.3.1437
Hao, G. F., and Yang, G. F. (2010). The Role of Phe82 and Phe351 in auxin-induced substrate perception by TIR1 ubiquitin ligase: A novel insight from molecular dynamics simulations. PLoS One 5:e10742. doi: 10.1371/journal.pone.0010742
Houben, M., and Van de Poel, B. (2019). 1-aminocyclopropane-1-carboxylic acid oxidase (ACO): The enzyme that makes the plant hormone ethylene. Front. Plant Sci. 10:695. doi: 10.3389/fpls.2019.00695
Hu, G., Grover, C. E., Arick, M. A. II, Liu, M., Peterson, D. G., and Wendel, J. F. (2021). Homoeologous gene expression and co-expression network analyses and evolutionary inference in allopolyploids. Brief. Bioinf. 22, 1819–1835. doi: 10.1093/bib/bbaa035
Ishitani, M., Xiong, L., Lee, H., Stevenson, B., and Zhu, J. K. (1998). HOS1, a genetic locus involved in cold-responsive gene expression in Arabidopsis. Plant Cell 10, 1151–1161. doi: 10.1105/tpc.10.7.1151
Islam, F., Farooq, M. A., Gill, R. A., Wang, J., Yang, C., Ali, B., et al. (2017). 2,4-D attenuates salinity-induced toxicity by mediating anatomical changes, antioxidant capacity and cation transporters in the roots of rice cultivars. Sci. Rep. 7:10443. doi: 10.1038/s41598-017-09708-x
Kalladan, R., Lasky, J. R., Chang, T. Z., Sharma, S., Juenger, T. E., and Verslues, P. E. (2017). Natural variation identifies genes affecting drought-induced abscisic acid accumulation in Arabidopsis thaliana. Proc. Natl. Acad. Sci. U. S. A. 114, 11536–11541. doi: 10.1073/pnas.1705884114
Li, F., Fan, G., Lu, C., Xiao, G., Zou, C., Kohel, R. J., et al. (2015). Genome sequence of cultivated Upland cotton (Gossypium hirsutum TM-1) provides insights into genome evolution. Nat. Biotechnol. 33, 524–530. doi: 10.1038/nbt.3208
Liu, X., Zhao, B., Zheng, H. J., Hu, Y., Lu, G., Yang, C. Q., et al. (2015). Gossypium barbadense genome sequence provides insight into the evolution of extra-long staple fiber and specialized metabolites. Sci. Rep. 5:14139. doi: 10.1038/srep14139
Ma, Q., Ding, Y., Chang, J., Sun, X., Zhang, L., Wei, Q., et al. (2014). Comprehensive insights on how 2,4-dichlorophenoxyacetic acid retards senescence in post-harvest citrus fruits using transcriptomic and proteomic approaches. J. Exp. Bot. 65, 61–74. doi: 10.1093/jxb/ert344
Marianayagam, N. J., Sunde, M., and Matthews, J. M. (2004). The power of two: Protein dimerization in biology. Trends Biochem. Sci. 29, 618–625. doi: 10.1016/j.tibs.2004.09.006
Markus, C., Pecinka, A., and Merotto, A. Jr. (2021). Insights into the role of transcriptional gene silencing in response to herbicide-treatments in Arabidopsis thaliana. Int. J. Mol. Sci. 22:3314. doi: 10.3390/ijms22073314
McCauley, C. L., McAdam, S. A. M., Bhide, K., Thimmapuram, J., Banks, J. A., and Young, B. G. (2020). Transcriptomics in Erigeron canadensis reveals rapid photosynthetic and hormonal responses to auxin herbicide application. J. Exp. Bot. 71, 3701–3709. doi: 10.1093/jxb/eraa124
Mockaitis, K., and Estelle, M. (2008). Auxin receptors and plant development: A new signaling paradigm. Annu. Rev. Cell Dev. Biol. 24, 55–80. doi: 10.1146/annurev.cellbio.23.090506.123214
Nevo, Y., and Nelson, N. (2006). The NRAMP family of metal-ion transporters. Biochim. Biophys. Acta 1763, 609–620. doi: 10.1016/j.bbamcr.2006.05.007
Oesch-Bartlomowicz, B., and Oesch, F. (2007). “Mechanisms of toxification and detoxification which challenge drug candidate and drugs,” in Comprehensive medicinal chemistry, eds B. Testa and H. van de Waterbeemd (Oxford: Elsevier), 193–214. doi: 10.1016/B0-08-045044-X/00124-3
Ohkawa, H., Tsuji, H., and Ohkawa, Y. (1999). The use of cytochrome P450 genes to introduce herbicide tolerance in crops: A review. Pestic. Sci. 55, 867–874. doi: 10.1002/(SICI)1096-9063(199909)55:9<867::AID-PS31>3.0.CO;2-S
Pace, H. C., and Brenner, C. (2001). The nitrilase superfamily: Classification, structure and function. Genome Biol. 2:REVIEWS0001. doi: 10.1186/gb-2001-2-1-reviews0001
Paley, S., Parker, K., Spaulding, A., Tomb, J. F., O’Maille, P., and Karp, P. D. (2017). The Omics Dashboard for interactive exploration of gene-expression data. Nucleic Acids Res. 45, 12113–12124. doi: 10.1093/nar/gkx910
Parry, G., Calderon-Villalobos, L. I., Prigge, M., Peret, B., Dharmasiri, S., Itoh, H., et al. (2009). Complex regulation of the TIR1/AFB family of auxin receptors. Proc. Natl. Acad. Sci. U. S. A. 106, 22540–22545. doi: 10.1073/pnas.0911967106
Patro, R., Duggal, G., Love, M. I., Irizarry, R. A., and Kingsford, C. (2017). Salmon: Fast and bias-aware quantification of transcript expression using dual-phase inference. Nat. Methods 14, 417–419. doi: 10.1038/nmeth.4197
Perez, L. M. (2021). The genetics and molecular mechanisms of tolerance to 2,4-dichlorophenoxyacetic acid (2,4-D) in upland cotton (Gossypium hirsutum L.). Ph.D. dissertation. Starkville, MS: Mississippi State University.
R Core Team (2020). R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing.
Raghavan, C., Ong, E. K., Dalling, M. J., and Stevenson, T. W. (2005). Effect of herbicidal application of 2,4-dichlorophenoxyacetic acid in Arabidopsis. Funct. Integr. Genomics 5, 4–17. doi: 10.1007/s10142-004-0119-9
Raghavan, C., Ong, E. K., Dalling, M. J., and Stevenson, T. W. (2006). Regulation of genes associated with auxin, ethylene and ABA pathways by 2,4-dichlorophenoxyacetic acid in Arabidopsis. Funct. Integr. Genomics 6, 60–70. doi: 10.1007/s10142-005-0012-1
Riar, D. S., Burke, I. C., Yenish, J. P., Bell, J., and Gill, K. (2011). Inheritance and physiological basis for 2,4-D resistance in prickly lettuce (Lactuca serriola L.). J. Agric. Food Chem. 59, 9417–9423. doi: 10.1021/jf2019616
Rieseberg, L. H., Archer, M. A., and Wayne, R. K. (1999). Transgressive segregation, adaptation and speciation. Heredity 83, 363–372. doi: 10.1038/sj.hdy.6886170
Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2010). EdgeR: A Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. doi: 10.1093/bioinformatics/btp616
Roudier, F., Fernandez, A. G., Fujita, M., Himmelspach, R., Borner, G. H. H., Schindelman, G., et al. (2005). COBRA, an Arabidopsis extracellular glycosyl-phosphatidyl inositol-anchored protein, specifically controls highly anisotropic expansion through its involvement in cellulose microfibril orientation. Plant Cell 17, 1749–1763. doi: 10.1105/tpc.105.031732
Sachu, M., Kynshi, B. L., and Syiem, M. M. (2022). A biochemical, physiological and molecular evaluation of how the herbicide 2, 4-dichlorophenoxyacetic acid intercedes photosynthesis and diazotrophy in the cyanobacterium Nostoc muscorum Meg 1. Environ. Sci. Pollut. Res. 29, 36684–36698. doi: 10.1007/s11356-021-18000-5
Saha, S., Jenkins, J. N., Wu, J., McCarty, J. C., Gutiérrez, O. A., Percy, R. G., et al. (2006). Effects of chromosome-specific introgression in Upland cotton on fiber and agronomic traits. Genetics 172, 1927–1398. doi: 10.1534/genetics.105.053371
Saha, S., Stelly, D. M., Raska, D. A., Wu, J., Jenkins, J. N., McCarty, J. C., et al. (2012). “Chromosome Substitution Lines: Concept, Development and Utilization in the Genetic Improvement of Upland Cotton,” in Plant Breeding, ed. I. Y. Abdurakhmonov (Croatia: InTech), 107–128. doi: 10.5772/35585
Salas-Perez, R. A., Saski, C. A., Noorai, R. E., Srivastava, S. K., Lawton-Rauh, A. L., Nichols, R. L., et al. (2018). RNA-Seq transcriptome analysis of Amaranthus palmeri with differential tolerance to glufosinate herbicide. PLoS One 13:e0195488. doi: 10.1371/journal.pone.0195488
Schulz, B., and Sehobye, K. (2016). 2,4-D transport and herbicide resistance in weeds. J. Exp. Bot. 67, 3177–3179. doi: 10.1093/jxb/erw199
Shyam, C., Jhala, A. J., Kruger, G., and Jugulam, M. (2019). Rapid metabolism increases the level of 2,4-D resistance at high temperature in common waterhemp (Amaranthus tuberculatus). Sci. Rep. 9:16695. doi: 10.1038/s41598-019-53164-8
Soneson, C., Love, M. I., and Robinson, M. D. (2016). Differential analyses for RNA-seq: Transcript-level estimates improve gene-level inferences. F1000Research 4:1521. doi: 10.12688/f1000research.7563.2
Song, Y. (2014). Insight into the mode of action of 2,4-dichorophenoxyacetic acid (2,4-D) as an herbicide. J. Integr. Plant Biol. 56, 106–113. doi: 10.1111/jipb.12131
Stelly, D., Saha, S., Raska, D., Jenkins, J., McCarty, J., and Gutierrez, O. (2005). Registration of 17 upland (Gossypium hirsutum) germplasm lines disomic or different G. barbadense chromosome or arm substitutions. Crop Sci. 45, 2663–2665. doi: 10.2135/cropsci2004.0642
Takahashi, M., Umetsu, K., Oono, Y., Higaki, T., Blancaflor, E. B., and Rahman, A. (2017). Small acidic protein 1 and SCFTIR1 ubiquitin proteasome pathway act in concert to induce 2,4-dichlorophenoxyacetic acid-mediated alteration of actin in Arabidopsis roots. Plant J. 89, 940–956. doi: 10.1111/tpj.13433
Teixeira, M. C., Fernandes, A. R., Mira, N. P., Becker, J. D., and Sá-Correia, I. (2006). Early transcriptional response of Saccharomyces cerevisiae to stress imposed by the herbicide 2,4-dichlorophenoxyacetic acid. FEMS Yeast Res. 6, 230–248. doi: 10.1111/j.1567-1364.2006.00041.x
Thimm, O., Bläsing, O., Gibon, Y., Nagel, A., Meyer, S., Kruger, P., et al. (2004). MAPMAN: A user-driven tool to display genomics data sets onto diagrams of metabolic pathways and other biological processes. Plant J. 37, 914–939. doi: 10.1111/j.1365-313X.2004.02016.x
Torra, J., Rojano-Delgado, A. M., Rey-Caballero, J., Royo-Esnal, A., Salas, M. L., and De Prado, R. (2017). Enhanced 2,4-D metabolism in two resistant Papaver rhoeas populations from Spain. Front. Plant Sci. 8:1584. doi: 10.3389/fpls.2017.01584
Wan, C. Y., and Wilkins, T. A. (1994). A modified hot borate method significantly enhances the yield of high-quality RNA from cotton (Gossypium hirsutum L.). Anal. Biochem. 223, 7–12. doi: 10.1006/abio.1994.1538
Wedding, R. T., Erickson, L. C., and Brannaman, B. L. (1953). Effect of 2,4-Dichlorophenoxyacetic Acid on Photosynthesis and Respiration, Paper No. 775. California: University of California Citrus Experiment Station. doi: 10.1104/pp.29.1.64
Werck-Reichhart, D., Hehn, A., and Didierjean, L. (2000). Cytochromes P450 for engineering herbicide tolerance. Trends Plant Sci. 5, 113–123. doi: 10.1016/S1360-1385(00)01567-3
Williams, G., and Dunn, S. (1966). Effects of 2,4-D on plant metabolism as modified by light quality. Weeds 14, 17–22. doi: 10.2307/4041115
Wright, T. R., Shan, G., Walsh, T. A., Lira, J. M., Cui, C., Song, P., et al. (2010). Robust crop resistance to broadleaf and grass herbicides provided by aryloxyalkanoate dioxygenase transgenes. Proc. Natl. Acad. Sci. U. S. A. 107, 20240–20245. doi: 10.1073/pnas.1013154107
Yang, S. F., and Hoffman, N. E. (1984). Ethylene biosynthesis and its regulation in higher plants. Annu. Rev. Plant Physiol. 35, 155–189. doi: 10.1146/annurev.pp.35.060184.001103
Keywords: herbicide resistance mechanism, plant physiology, abiotic stress tolerance, transcriptome analysis, RNA sequencing, auxin response and signaling, herbicide metabolism, Illumina sequencing
Citation: Perez LM, Mauleon R, Arick MA II, Magbanua ZV, Peterson DG, Dean JFD and Tseng TM (2022) Transcriptome analysis of the 2,4-dichlorophenoxyacetic acid (2,4-D)-tolerant cotton chromosome substitution line CS-B15sh and its susceptible parental lines G. hirsutum L. cv. Texas Marker-1 and G. barbadense L. cv. Pima 379. Front. Plant Sci. 13:910369. doi: 10.3389/fpls.2022.910369
Received: 01 April 2022; Accepted: 21 July 2022;
Published: 22 August 2022.
Edited by:
Linghe Zeng, United States Department of Agriculture (USDA), United StatesReviewed by:
Ram Kumar Basnet, Rijk Zwaan, NetherlandsRagunathan Devendran, The University of Georgia, Tifton Campus, United States
Jonathan Elias Maldonado, University of Santiago, Chile
Copyright © 2022 Perez, Mauleon, Arick, Magbanua, Peterson, Dean and Tseng. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Te Ming Tseng, t.tseng@msstate.edu
†Present address: Loida M. Perez, Philippine Rice Research Institute, Muñoz, Philippines