Transcriptomic Analysis Identifies New Non-Target Site Glyphosate-Resistance Genes in Conyza bonariensis

Conyza bonariensis (hairy fleabane) is one of the most problematic and widespread glyphosate-resistant weeds in the world. This highly competitive weed species significantly interferes with crop growth and substantially decreases crop yield. Despite its agricultural importance, the molecular mechanisms of glyphosate resistance are still unknown. The present RNA-Seq study was performed with the goal of identifying differentially expressed candidate transcripts (genes) related to metabolism-based non-target site glyphosate resistance in C. bonariensis. The whole-transcriptome was de novo assembled from glyphosate-resistant and -sensitive biotypes of C. bonariensis from Southern Brazil. The RNA was extracted from untreated and glyphosate-treated plants at several timepoints up to 288 h after treatment in both biotypes. The transcriptome assembly produced 90,124 contigs with an average length of 777 bp and N50 of 1118 bp. In response to glyphosate treatment, differential gene expression analysis was performed on glyphosate-resistant and -sensitive biotypes. A total of 9622 genes were differentially expressed as a response to glyphosate treatment in both biotypes, 4297 (44.6%) being up- and 5325 (55.4%) down-regulated. The resistant biotype presented 1770 up- and 2333 down-regulated genes while the sensitive biotype had 2335 and 2800 up- and down-regulated genes, respectively. Among them, 974 up- and 1290 down-regulated genes were co-expressed in both biotypes. In the present work, we identified 41 new candidate target genes from five families related to herbicide transport and metabolism: 19 ABC transporters, 10 CYP450s, one glutathione S-transferase (GST), five glycosyltransferases (GT), and six genes related to antioxidant enzyme catalase (CAT), peroxidase (POD), and superoxide dismutase (SOD). The candidate genes may participate in metabolic-based glyphosate resistance via oxidation, conjugation, transport, and degradation, plus antioxidation. One or more of these genes might ‘rescue’ resistant plants from irreversible damage after glyphosate treatment. The 41 target genes we report in the present study may inform further functional genomics studies, including gene editing approaches to elucidate glyphosate-resistance mechanisms in C. bonariensis.


Results
The present study aimed to study the molecular changes of C. bonariensis in response to glyphosate treatment in the long term after herbicide treatment. The timepoints of RNA extraction to RNA sequencing were determined based on shikimic acid content in the GR biotype from 24 h to 288 h after treatment. In the last timepoint, the shikimic acid content on treated GR plants did not differ regarding to untreated plants. Thus, the gene expression changes were evaluated in GR glyphosate treated and untreated plants in contrast to GS with the same treatments.

Illumina Sequencing and De Novo Assembly
The twelve cDNA libraries sequenced from C. bonariensis leaves produced a high-quality assembled dataset (Table S1). The Phred scores of the de novo assembled data were ≥Q30 level (error probability 0.1%) for more than 90% of all raw reads. Transcriptome assembling produced 203,054 (>200 bp) transcripts and 90,124 contigs at the "gene" level (N50 of 1118 bp and the average length of 777 bp) ( Table S1). The sequence distribution of transcripts length threshold varied from 0 to 3000 bp ( Figure S1). In this way, the large scale and high-quality results provided in the present study will serve as molecular reference for further C. bonariensis studies.

Functional Annotation of Assembled Contigs
The BLASTx results of sequences indicated that 68.1% of the annotated contigs (90,124) had hits with nucleotide sequences from the NCBI database ( Figure S2). Most of the annotated nucleotide sequences of C. bonariensis transcriptome corresponded to Helianthus annuus (36.9%) and Cynara cardunculus var. scolymus (18.8%). Further, the contig identity to other plant species was 12.4%, and the number of novel genes without interspecies identity was 31.9% ( Figure S3).
The functions of C. bonariensis contigs were classified according to gene ontology assignments. A total of 11,877 contigs (13.2%) were attributed to at least one gene ontology term and classified in 30 functional categories using the complete set of GO terms in three main categories: molecular function (n = 3188, 26.9%), biological process (n = 7491, 63.1%), and cellular component (n = 1198, 10.1%) ( Figure S4). The largest proportion of genes was represented by a cellular (3.7%) and metabolic process (3.2%) in the biological process; molecular function (4.7%) and binding (3.5%) in molecular function; and a cellular component (4.7%) and cell organelle (4.2%) in cellular component ( Figure S4). The high number of novel genes (31.9%) provided in the present RNA sequencing results shows the importance of characterizing genomes of weeds and also serves to enable follow-on molecular studies in C. bonariensis and other weedy species. Further genome sequencing studies in Conyza spp. will be helpful to completely characterize the novel genes presented here. Further, the gene ontology indicates that the biological process represented the highest proportion of annotated genes regarding molecular function and cellular component.

Differentially Expressed Genes
The total of differentially expressed genes (DEGs) (9622 "genes"; 39.9% of total) was found in response to glyphosate treatment in both biotypes: 4297 (44.6%) genes were up-regulated and 5325 (55.4%) were down-regulated ( Figure 1 and Figure S5). Among the DEGs, 1770 and 2333 genes were up-and down-regulated in GR plants, while 2335 and 2800 genes were up-and down-regulated in GS plants, respectively. A total of 974 up-and 1290 down-regulated genes were co-expressed in both biotypes as a response to glyphosate treatment ( Figure 1, Figure S5). In GR and GS biotypes without glyphosate treatment (t0), there were around 100 genes each in up-and down-regulation, with no apparent co-expression ( Figure 1).

Differentially Expressed Genes
The total of differentially expressed genes (DEGs) (9,622 "genes"; 39.9% of total) was found in response to glyphosate treatment in both biotypes: 4,297 (44.6%) genes were up-regulated and 5,325 (55.4%) were down-regulated ( Figure 1 and Figure S5). Among the DEGs, 1,770 and 2,333 genes were up-and down-regulated in GR plants, while 2,335 and 2,800 genes were up-and down-regulated in GS plants, respectively. A total of 974 up-and 1,290 down-regulated genes were co-expressed in both biotypes as a response to glyphosate treatment ( Figure 1, Figure S5). In GR and GS biotypes without glyphosate treatment (t0), there were around 100 genes each in up-and down-regulation, with no apparent co-expression ( Figure 1). The most abundant GO terms for up-regulated genes were biological process -protein metabolism (87 DEGs from GR plants and 120 DEGs from GS plants) (Figure 2A Figure 2E,F). It is interesting to note that GR plants have 45.9% more upregulated genes in transmembrane transport than GS (37 vs. 17). However, in the GR biotype, this process was classified as the fifth-most prominent and was not included among the top ten upregulated genes in GS (Figure 2A,B). On the other hand, the GS biotype had 126 (50%) more upregulated genes than GR (376 vs. 250) in processes located in the cytoplasm (cellular component) ( Figure 2E,F). Further, the resistant biotype had 53 (36%) more up-regulated genes that are chloroplast localized (200 vs. 147) than in GS plants ( Figure 2E,F). The most abundant GO terms for up-regulated genes were biological process-protein metabolism (87 DEGs from GR plants and 120 DEGs from GS plants) (Figure 2A Figure 2E,F). It is interesting to note that GR plants have 45.9% more up-regulated genes in transmembrane transport than GS (37 vs. 17). However, in the GR biotype, this process was classified as the fifth-most prominent and was not included among the top ten up-regulated genes in GS (Figure 2A,B). On the other hand, the GS biotype had 126 (50%) more up-regulated genes than GR (376 vs. 250) in processes located in the cytoplasm (cellular component) ( Figure 2E,F). Further, the resistant biotype had 53 (36%) more up-regulated genes that are chloroplast localized (200 vs. 147) than in GS plants ( Figure 2E,F). The attributed GO functions for down-regulated genes show the highest proportions of characterized genes respectively in GR and GS related to biological process -cell wall metabolism (3.8%) and protein metabolism (3.5%) ( Figure 3A,B); molecular function -ATP binding for both biotypes (8.1% and 10.1%) ( Figure 3C,D); and cellular component -integral component of membrane (9,6%) and chloroplast (13.4%) ( Figure 3E,F). In the biological process, 63 genes related to translation were down-regulated in GR while just 16 in GS, which do not figure in the top ten in the sensitive biotype. Further, the transmembrane transport-related contigs were not among the top ten down-regulated DEGs in GR but were so in GS ( Figure 3A,B). In molecular function, there were 52 down-regulated genes related to ribosome metabolism in GR while just one in GS ( Figure 3C,D). Further, GR presented 40% and 31% less down-regulated genes related to chloroplast and integral component of the membrane than GS ( Figure 3E,F).
The high number of DEGs represents the high molecular changes that glyphosate action causes in the C. bonariensis transcriptome. Further, the differences on expressed genes according to each category show differences between GR and GS biotypes and might indicate the processes related with resistance to glyphosate. The attributed GO functions for down-regulated genes show the highest proportions of characterized genes respectively in GR and GS related to biological process-cell wall metabolism (3.8%) and protein metabolism (3.5%) ( Figure 3A,B); molecular function-ATP binding for both biotypes (8.1% and 10.1%) ( Figure 3C,D); and cellular component-integral component of membrane (9,6%) and chloroplast (13.4%) ( Figure 3E,F). In the biological process, 63 genes related to translation were down-regulated in GR while just 16 in GS, which do not figure in the top ten in the sensitive biotype. Further, the transmembrane transport-related contigs were not among the top ten down-regulated DEGs in GR but were so in GS ( Figure 3A,B). In molecular function, there were 52 down-regulated genes related to ribosome metabolism in GR while just one in GS ( Figure 3C,D). Further, GR presented 40% and 31% less down-regulated genes related to chloroplast and integral component of the membrane than GS ( Figure 3E,F).
The high number of DEGs represents the high molecular changes that glyphosate action causes in the C. bonariensis transcriptome. Further, the differences on expressed genes according to each category show differences between GR and GS biotypes and might indicate the processes related with resistance to glyphosate.

RNA-Seq Validation by qRT-PCR Analysis
The average expression for 19 contigs with relative expression performed in a time course indicates an amplitude of (log fold change-logFC) from 0.01 to 1990, while the transcriptome results for the same contigs ranged from logFC 0.1 to 3625 (Figure 4 and S6). The transcriptome dataset and qRT-PCR presented a significant correlation (r) of 0.9 for transcriptional results of GR and GS biotypes ( Figure S6). On the other hand, the results of the time course qRT-PCR experiments show that gene expression response to glyphosate treatment varies according to gene-related function, to biotype, and to time after herbicide exposure ( Figure 4). In general, the results of the evaluated genes in a time course experiment suggest that the highest expression levels occur between 96 and 192 hours after treatment (HAT) (Figure 4). The differences of gene expression in time might be a very good reference for further studies involving glyphosate resistance; between 96 and 192 HAT is the best time-window to capture molecular responses to glyphosate treatment in C. bonariensis.

RNA-Seq Validation by qRT-PCR Analysis
The average expression for 19 contigs with relative expression performed in a time course indicates an amplitude of (log fold change-logFC) from 0.01 to 1990, while the transcriptome results for the same contigs ranged from logFC 0.1 to 3625 ( Figure 4 and Figure S6). The transcriptome dataset and qRT-PCR presented a significant correlation (r) of 0.9 for transcriptional results of GR and GS biotypes ( Figure S6). On the other hand, the results of the time course qRT-PCR experiments show that gene expression response to glyphosate treatment varies according to gene-related function, to biotype, and to time after herbicide exposure ( Figure 4). In general, the results of the evaluated genes in a time course experiment suggest that the highest expression levels occur between 96 and 192 h after treatment (HAT) (Figure 4). The differences of gene expression in time might be a very good reference for further studies involving glyphosate resistance; between 96 and 192 HAT is the best time-window to capture molecular responses to glyphosate treatment in C. bonariensis.

EPSPS Sequence Analysis
The BLASTn analysis identified, with an E-value threshold of zero and score bits >600, one single copy of each the three EPSPS sequences assembled in each one of the individual C. bonariensis transcriptomes of GR and GS. There was no amino-acid substitution at Thr 102 and Pro 106 codons from EPSPS sequences alignment ( Figure S7A-C). On the other hand, the basal EPSPS transcription varied among the three copies but did not between biotypes. In general, EPSPS2 presented the lowest level of expression among the three copies. In response to glyphosate treatment, the EPSPS2 and EPSPS3 had increased transcription; however, there were no differences between biotypes. EPSPS1 had low variation in expression ( Figure S7D). These results indicate that the glyphosate-resistance process in GR plants is due to a NTSR.

Candidate NTSR Genes Related to Putative Functions of Glyphosate Transport and Metabolism
We identified a total of 41 candidates differentially expressed genes associated to the herbicide metabolism phases of oxidation (CYP450), conjugation (GST and GT), transport (ABC transporters), and protection/compensation (antioxidant enzyme-CAT, POD, and SOD) [20,30]. Two ABC transporters and one GT target gene were chosen for validation by qRT-PCR. From these results, the ABC transporters' transcript abundance varied between biotypes in response to glyphosate treatment and in time after herbicide exposure. The expression of ABC_12 showed the highest expression in GR rather than GS at 24 h after treatment (HAT) (fold change 29.6), whereas ABC_16 at 192 HAT (fold change 3.4). The GT_5 highest expression in GR about GS occurred at 192 HAT with a fold change of 2.8 ( Figure 4). Further, eight genes had their highest expression in GS plants, with a more significant increase in expression between 96 and 192 HAT. Among these eight genes, we highlighted the steroid 5-beta-reductase (fold change 11.1-96 HAT) and aspartyl protease (fold change 281.5-192 HAT) ( Figure 4).
We identified two copies of each ABC transporters M10 and M11 gene (M10_g1 and M10_g2; M11_g1 and M11_g2) in the C. bonariensis transcriptome. Except for M10_g2, all other genes had increased transcript abundance from glyphosate treated plants. The GS biotype presented higher M10 and M11 expression levels than the resistant biotype ( Figure S8). Because of that, M10 and M11 genes were not considered glyphosate-resistance candidate genes in the present work.
Here, we report 19 new ABC transporter candidate genes associated with glyphosate resistance in C. bonariensis. In the differential expression analysis, these genes had a higher response to glyphosate treatment in GR than in GS with different levels of transcripts abundance (Table 1, Figure 5). Fourteen ABC transporter genes (ABC_2-4, ABC_7-9, ABC_12-19) responded to treatment in both biotypes, but with higher transcription in the resistant biotype. However, two genes (ABC_5 and ABC_10) appeared to have reduced expression in GS plants but increased expression in the GR biotype; and three genes (ABC_1, ABC_6, and ABC_11) had increased transcription only in GR plants with a fold change of 3.1, 3.7, and 8, respectively, by comparing GR-and GS treated biotypes ( Figure 5, Table 1). In general, the ABC_2, ABC_5, ABC_7, and ABC_10-12 genes presented the highest expression ratio in the resistant biotype, which varied from 6-to 8-fold (Table 1).
Among the 10 CYP450 genes we report in the present work in response to glyphosate treatment, four (CYP450_1, CYP450_3, CYP450_4, and CYP450_7) had expression suppressed in the GS biotype and increased in GR; three (CYP450_8-10) increased the transcription in both biotypes, although with higher levels in GR; and three (CYP450_2, CYP450_5, and CYP450_6) showed higher expression only in the resistant biotype with a fold change of 3.3, 5.9, and 4.5, respectively, in comparison between biotypes with treatment (Table 2; Figure 6). The CYP450_1 and CYP450_4 presented the highest differences in GR rather than GS in response to treatment, which was 11.5-and 27-fold, respectively ( Table 2). Table 1. Differential expression of ABC transporter genes via RNA-Seq analysis of glyphosate-resistant (GR) and -sensitive (GS) C. bonariensis biotypes in response to glyphosate treatment up to 288 h after treatment. The transcriptional fold change data for the GR and GS biotypes are each expressed as the ratio of the ABC transporter gene expression from glyphosate treated to untreated individuals. Therefore, a value of 1 means the gene was neither induced nor repressed. The rightmost column represents the ratio of ratios; i.e., the magnitude of effect of being glyphosate-resistant.   . GR: glyphosate-resistant biotype; GS: glyphosatesensitive biotype. t0: without glyphosate treatment; t1: with glyphosate treatment -RNA was obtained from plants collected at several timepoints up to 288 h after treatment and pooled. Intervals indicate the standard error. In ABC_12 (ABC transporter C family member 14) and ABC_14 (ABC transporter C family member 9), the qRT-PCR analysis was performed (see Figure 4). The phylogenetic tree of these genes is presented in Figure S10.   14) and ABC_14 (ABC transporter C family member 9), the qRT-PCR analysis was performed (see Figure 4). The phylogenetic tree of these genes is presented in Figure S10.  Figure S10.
The transcription level of glutathione-related gene increased only in the GR biotype in response to treatment with a ratio between biotypes of 4.6 (Table 3, Figure 7). The basal GST transcription level in the GR biotype was around 50% less than in GS (Figure 7). The expression levels for the five glycosyltransferase genes increased in both biotypes as a response to glyphosate treatment, although with a higher response in the GR biotype. In this way, the ratio between biotypes was almost doubled in GR rather than GS, except GT_4, which was of 4.8 ( Figure 7, Table 3). GT_5 was evaluated in qRT-PCR in a time-course experiment and presented the highest response to treatment at 192 HAT in the GR biotype, while in GS, the expression levels had low variations ( Figure 4). Table 2. Differential expression of cytochrome P450 (CYP450) genes via RNA-Seq analysis of glyphosate-resistant (GR) and -sensitive (GS) C. bonariensis biotypes in response to glyphosate treatment up to 288 h after treatment. The transcriptional fold change data for the GR and GS biotypes are each expressed as the ratio of the ABC transporter gene expression from glyphosate treated to untreated individuals. Therefore, a value of 1 means the gene was neither induced nor repressed. The rightmost column represents the ratio of ratios; i.e., the magnitude of effect of being glyphosate-resistant.  Differentially expressed genes were selected using a p-value and FDR threshold set at ≤0.001; t0 = without glyphosate treatment; t1 = with glyphosate treatment-RNA was obtained from plants collected at several timepoints up to 288 h after treatment and pooled; GR: glyphosate-resistant biotype; GS: glyphosate-sensitive biotype. The phylogenetic tree of these genes is presented in Figure S10. Table 3. Differential expression of glutathione (GST) and glycosyltransferase (GT) genes via RNA-Seq analysis of glyphosate-resistant (GR) and -sensitive (GS) C. bonariensis biotypes in response to glyphosate treatment up to 288 h after treatment. The transcriptional fold change data for the GR and GS biotypes are each expressed as the ratio of the ABC transporter gene expression from glyphosate treated to untreated individuals. Therefore, a value of 1 means the gene was neither induced nor repressed. The rightmost column represents the ratio of ratios; i.e., the magnitude of effect of being glyphosate-resistant. Differentially expressed genes were selected using a p-value and FDR threshold set at ≤0.001; t0 = without glyphosate treatment; t1 = with glyphosate treatment-RNA was obtained from plants collected at several timepoints up to 288 h after treatment and pooled; GR: glyphosate-resistant biotype; GS: glyphosate-sensitive biotype. ** indicates that the gene was subjected to qRT-PCR analysis (See Figure 4).   Figure 4). Intervals indicate the standard error.
The two catalases had the highest expression after glyphosate treatment in both biotypes among all 41 candidate genes. The CAT_1 gene had >60,000 TPM, while CAT_2 was >7000. Glyphosate-treated GR plants had CAT_1, and CAT_2 increased transcription by 4.8-and 8.5-fold, whereas the GS biotype had increases of 1.2-and 2.3-fold, respectively ( Figure 8A,B, Table 4). The comparison between treated biotypes indicates the higher expression levels of CAT_1 and CAT_2 in the resistant biotype at a ratio of 4-and 3.7-fold, respectively (Table 4). With regard to the two differentially expressed peroxidase genes, POD_1 presented a high response to glyphosate treatment in both biotypes, whereas POD_2 had increased expression in the GR biotype and was suppressed in GS plants ( Figure 8C,D, Table 4). After glyphosate treatment, POD_1 increased the expression in GR biotype at 11.9-fold, while 4.6-fold in GS (ratio of 2.6). POD_2 expression increased 4.5-fold in GR and reduced at 0.2-fold in GS, presenting a ratio of 22.5 ( Table 4). The gene expression related to SOD was higher in the GR biotype than in GS ( Figure 8E,F, Table 4). SOD_2 presented high basal expression in the GR biotype, which reduced 28% after glyphosate treatment; however, there was 2.7 times higher expression than in GS plants ( Figure 8F). Table 4. Differential expression of catalase (CAT), peroxidase (POD), and superoxide dismutase (SOD) genes via RNA-Seq analysis of glyphosate-resistant (GR) and -sensitive (GS) C. bonariensis biotypes in response to glyphosate treatment up to 288 h after treatment. The transcriptional fold change data for the GR and GS biotypes are each expressed as the ratio of the ABC transporter gene expression from glyphosate treated to untreated individuals. Therefore, a value of 1 means the gene was neither induced nor repressed. The rightmost column represents the ratio of ratios; i.e., the magnitude of effect of being glyphosate-resistant.   In summary, the candidate gene results indicate that glyphosate resistance is a very complex process of plant defense against glyphosate action involving several molecular changes that allow GR plant survival. The new 41 candidate genes presented in the present study comprise a four-phase of a well-known process of herbicide resistance (oxidation; conjugation; transport; and degradation, detoxification, and protection).

The C. Bonariensis Transcriptome
The transcriptomic profiles performed up to 288 hours after treatment (HAT) will assist in untangling the molecular processes underlying the glyphosate NTSR mechanisms in C. bonariensis. The current transcriptome assembly produced a large dataset, which likely includes the major candidate genes involved in glyphosate resistance in Conyza. The number of assembled contigs at the gene level in our study is similar in scope to those obtained in another C. bonariensis transcriptome study (~81,000) [16].
The gene ontology assignments of the differentially expressed genes (DEGs) showed interesting differences between biotypes. The GR biotype had a higher number of up-regulated genes related to transmembrane transport (biological process) and chloroplast (cellular component) than GS. On the other hand, the GS biotype presented a higher number of down-regulated genes related to these In summary, the candidate gene results indicate that glyphosate resistance is a very complex process of plant defense against glyphosate action involving several molecular changes that allow GR plant survival. The new 41 candidate genes presented in the present study comprise a four-phase of a well-known process of herbicide resistance (oxidation; conjugation; transport; and degradation, detoxification, and protection).

The C. Bonariensis Transcriptome
The transcriptomic profiles performed up to 288 h after treatment (HAT) will assist in untangling the molecular processes underlying the glyphosate NTSR mechanisms in C. bonariensis. The current transcriptome assembly produced a large dataset, which likely includes the major candidate genes involved in glyphosate resistance in Conyza. The number of assembled contigs at the gene level in our study is similar in scope to those obtained in another C. bonariensis transcriptome study (~81,000) [16].
The gene ontology assignments of the differentially expressed genes (DEGs) showed interesting differences between biotypes. The GR biotype had a higher number of up-regulated genes related to transmembrane transport (biological process) and chloroplast (cellular component) than GS. On the other hand, the GS biotype presented a higher number of down-regulated genes related to these processes than GR (Figure 2; Figure 3). In our study, we focus on up-regulated transcripts in HR hairy fleabane, which indicate the putative importance of the transmembrane transport process mediated by ABC transporters in conferring glyphosate resistance. Upon application, glyphosate first impacts chloroplast biology [31]; therefore, the increase in transcription of genes related to chloroplast function in the GR biotype suggests the plant may be mounting a potential defense to protect plastids against glyphosate's action. These results are in accordance with a study with proteomic analysis that identified chloroplast proteins differentially expressed as the primary sites of GR in C. canadensis [32]. Thus, the results of the present study suggest that transmembrane transport and chloroplast proteins might act in association to reduce damages caused by glyphosate action and allow the GR biotype to survive.
The substantially higher number of down-regulated genes related to translation (63 DEGs in GR vs. 16 in GS) and ribosome metabolism (52 DEGs in GR vs. 1 in GS) in GR suggest that in response to glyphosate action, the protein production was reduced in GR ( Figure 3). Thus, one coping mechanism evolved in the GR biotype may be induced metabolic slowdown, which ameliorates the plant´defense against glyphosate movement and damage.
The EPSPS sequence analysis revealed no nucleotide substitutions at Thr 102 and Pro 106 codons and no increase in expression and number of copies in GR plants ( Figure S7). Therefore, we conclude glyphosate resistance in the GR biotype is NTSR. Currently, few mutations that confer resistance to glyphosate have been observed in weeds [19]; ergo, glyphosate resistance across weedy plants is also NTSR. Few peptide changes in EPSPS are tolerated without loss of function because of most mutations are lethal [14,16,19]. Our results are in accordance with those related by Hereward et al. (2018) [16].

Candidate Genes Involved in Glyphosate Resistance in C. Bonariensis
A group of well-established NTSR gene families is the cytochrome (CYP450) monooxygenases, glutathione S-transferases (GSTs), and glycosyltransferases (GTs) [20,25,31,33]. These enzyme classes have versatile and multifunctional activities and might catalyze the oxidation and herbicide conjugation to variable substrates [20]. Further, a piece of emerging knowledge about the antioxidant system complementing the glyphosate resistance process has been described [34]. Metabolic-based herbicide resistance in general follows a four-phase process: I-oxidation; II-conjugation; III-transport; and IV-degradation, detoxification, and protection [14,20,30,35]. In the present work, we report transcriptional increases of the genes belonging to these four processes in a GR biotype. This increase in transcription of the enzyme-coding gene indicates a potential enhancement of activity, and subsequently glyphosate metabolization to allow GR plants to survive after herbicide exposure. To our knowledge, this is the first report that associated the CYP450, ABC transporters, GT, GST, and antioxidant enzyme differential gene expression to glyphosate resistance in C. bonariensis.
In phase I, CYP450 acts to oxidize the herbicide molecules and expose certain functional groups to phase II enzymes [20,25,36]. The P450s catalyze several reactions in plant metabolism, and their role in herbicide conversion is through hydroxylation or dealkylation [14,37,38]. P450s use electrons from NADPH (NADPH-P450 reductase) to insert an oxygen atom in the herbicide molecule, producing a more suitable product to be conjugated to glucose and transported to the vacuole [36]. The CYP450 enzyme has been reported to be involved in herbicide resistance to herbicides belonging to groups of acetanilide, aryloxyphenoxy-propanoate, imidazolinone, phenoxyalkanoic-acid, phenylurea, sulfonamide, sulfonylurea, bentazon, and clomazone herbicides [14,33,37,38]. However, to our knowledge, there is no report associating CYP450 to glyphosate resistance. In the present work, we report the up-regulation of 10 CYP450 genes in the GR biotype followed by glyphosate treatment showing the potential to be involved with GR in C. bonariensis ( Figure 6, Table 2). Further, the annotated families of the 10 CYP450 contigs in the present work (Table 2) were not described by other authors conferring resistance to herbicide [33,36,37]. According to Powles and Yu (2010) [14], identifying the P450s conferring herbicide resistance in weeds is a vast research frontier, and this group of enzymes may indeed play significant roles in resistance owing to their ability to affect several herbicides modes of action.
In phase II, GSTs and GTs detoxify herbicides through direct conjugation or conjugate it to a wide range of substrates [20,35,39]. GSTs are found in the cytoplasm at high concentrations and work by catalyzing the glutathione conjugation to a several substrates, producing a polar product [20,39,40]. In this way, the up-regulation of GST and GT genes in the GR biotype followed by glyphosate treatment in the C. bonariensis might indicate the increase in activity of these enzymes to cope with glyphosate ( Figure 7). GTs are a huge gene family in which proteins conjugate a sugar to a large variety of lipophilic molecules and herbicides [20,41]. GTs are in the cytoplasm and act by transferring sugar to lipophilic molecules enabling access to membrane transporters, including ABC transporters [41].
In phase III, ABC transporters may direct herbicides to be compartmentalized in the vacuoles or extracellular spaces after the action of CYP450, GSTs, and GTs and thus contribute to herbicide resistance process. ABC transporters are membrane-targeted proteins that require ATP for active transport of substrates across membranes [20,30]. Studies have been reported that ABC transporters are involved in glyphosate resistance in horseweed [23,24,42]. After analysis of the M10 and M11 ABC transporters genes previously reported in the literature as being involved in GR in C. canadensis [23], the results of the present study indicate that M10 and M11 genes do not appear to be differentially regulated in hairy fleabane and do not play a role in glyphosate resistance, which is consistent with other studies [16,21,32].
On the other hand, we present 19 new candidate target genes that might be involved in NTSR GR ( Figure 5, Table 1). The relatively high number of ABC transporter genes with higher expression in the resistant biotype relative to the sensitive biotype indicate that this gene family may play an important role in the glyphosate-resistance process. Some of them were responsive to treatment in both biotypes, although with higher transcription in the GR biotype. The most interesting ABC genes are ABC_1, ABC_6, and ABC_11, which were only responsive to glyphosate treatment in the resistant biotype with a fold change ratio between biotypes of 3.1, 3.7, and 8, respectively ( Figure 5, Table 1). Further, the time-course qRT-PCR results for two ABC transporter contigs indicates that the stimulus for an increase in expression levels might vary among members of this family of the gene at different times after herbicide exposure. For example, the ABC_12 had the highest expression at 24 HAT, while that for ABC_16 occurred at 192 HAT ( Figure 4). These results clearly show that the gene transcription responses to glyphosate treatment do not occur simultaneously. A recent study involving ABC transporters genes in horseweed showed the importance of the timing of gene overexpression initiation for the resistance mechanism itself [42]. The same authors described that early overexpression of ABC transporters genes in a short time followed by glyphosate treatment could be more useful to sequester the herbicide to vacuoles or extracellular space. In this case, the delayed overexpression could not offer enough protection. Another transcriptional study with ABC transporter genes showed the influence of the environmental conditions such as glyphosate dosage and time after treatment, light, and temperature on gene induction [43]. Thus, several ABC transporters contigs seem to be gene stress-responsive because of the increase in the gene expression in both GR and GS biotypes. However, some of them presented higher regulation in the GR biotype (ABC_1, ABC_6, and ABC_11) and can be considered as related to the glyphosate-resistance mechanism in C. bonariensis.
If HR-associated ABC transporters sequester glyphosate into vacuoles or apoplasts [44], glyphosate would have mitigated early damage to sprayed tissues than in non-sequestered cells, and even more important, this would prevent effective transport to non-sprayed tissues. Indeed, this is likely the major basis of NTSR to glyphosate. In addition, if the transporters can also act in different secondary metabolites produced after glyphosate treatment, their early and late activity might help to ameliorate secondary stress.
In phase IV, the oxidized herbicides by the CYP450 action and subsequently conjugated molecule are degraded in the vacuole or extracellular space, resulting in less phytotoxic molecules [20]. Further, in this stage, the antioxidant system works to scavenge reactive oxygen species (ROS) produced by the glyphosate action.
It is well known that glyphosate inhibits EPSPS to short-circuit the shikimic acid pathway, which results in decreased aromatic amino acid (phenylalanine, tyrosine, and tryptophan) biosynthesis. In so doing, shikimic acid accumulates, which ultimately leads to oxidative stress via ROS production [45]. ROS are highly reactive toxic molecules causing damage to cellular structures and cell death [46]. Further, the lack of tyrosine inhibits the synthesis of plastoquinone, which is an electron acceptor in the photosynthetic electron transport chain in photosystem two (PSII). The non-regeneration of plastoquinone interrupts the electron transport in PSII. Therefore, GS plants die from numerous factors beyond amino acid shortages [31].
Plants have an efficient enzymatic and non-enzymatic antioxidant system in response to ROS production [46,47]. The antioxidant enzyme-coding genes with differential expression in the present study were SOD, CAT, and POD ( Figure 8). SOD removes the superoxide (O 2 •− ) through catalyzing its dismutation and providing the first line of defense against ROS [46]. CAT can directly dismutase H 2 O 2 into H 2 O and O 2 and is indispensable for ROS detoxification during stressed conditions [47,48]. PODs play an essential role in scavenging ROS and protecting cells in higher plants through scavenging H 2 O 2 in water-water and ascorbate-reduced glutathione cycles and utilize ascorbate as the electron donor [47]. The increase in the transcription of the antioxidant enzyme-coding genes SOD, CAT, and POD in resistant plants are in accordance with the previous study we reported for these enzymes' activities; that study was performed using the same C. bonariensis as that used in this study [28].
In that study, we demonstrated the higher antioxidant enzyme activity in GR than in GS and less cellular damage after glyphosate treatment. Further, other studies of glyphosate-resistance mechanisms performed in A. trifida and A. palmeri concluded that antioxidant enzymes have the potential to play a role in the resistance processes [34,49]. Indeed, in the described antioxidant enzyme action, GSTs can also cope with ROS. GSTs are isozymes known to protect the cells against chemical-induced toxicity. These enzymes catalyze the conjugation of GSH to a variety of electrophilic and hydrophobic substrate [50]. Accordingly, the action of the cellular antioxidant machinery is essential to control excess ROS to protect plant cells from oxidative damage and to restore the redox homeostasis.

Plant Accessions, Experimental Treatments, and RNA Isolation
Experiments were performed on GR and GS C. bonariensis biotypes from Pelotas, the Rio Grande do Sul State, Brazil, 32 • 04 05.91 S, 52 • 52 59.14 W. The resistance factor was twice determined to be 18.4 [28]. Seeds of each biotype were germinated in plastic trays containing sterilized soil and substrate (Mac plant-Mec Prec, Brazil) 3:1 and watered daily in a greenhouse at 30 • C (day) /20 • C (night) (±4 • C) with a 12-h photoperiod. Thirty days after emergence (30 DAE), seedlings of each biotype were transplanted to plastic pots of 3 L containing soil-substrate mix and kept in the same greenhouse and conditions. Thirty days after transplant (60 DAE; rosette stage-plants 6 to 8 cm in diameter), plants were treated with glyphosate (1480 g ae ha −1 -Roundup Original DI 370 g ae L −1 ; Monsanto) with a CO 2 sprayer and 150 L ha −1 of spray volume.
The second and third leaves (from apex) from untreated and glyphosate-treated plants were collected in a range of timepoints: GR-untreated, 24,96,192, and 288 h after treatment; GS-untreated, 24, 96, and 192 h after treatment (Figure 9). Leaves from three biological replicates were collected at each time from a total of 27 plants (Figure 9). After collection, leaves were immediately frozen in liquid nitrogen and stored at −80 • C. The times to leaf collection and RNA extraction were determined based on a shikimic acid transient accumulation in the GR biotype, which increased substantially at 24 h after treatment, peaked at 96 h, and at 288 h presented a similar content to that of the untreated plants [28] ( Figure S9). By contrast, the GS biotype had similar initial shikimic acid levels as GR plants at 24 h and 96 h, but it peaked at 192 h after treatment and died after that timepoint.
additional treatment with RNase-free DNase I (Invitrogen) was performed to remove residual genomic DNA. The RNA samples from each respective biotype and treatment (untreated and glyphosate treated) were pooled to form three technical replicates per treatment totaling 12 samples, i.e., GR untreated (GR t0), GS untreated (GS t0), GR treated (GR t1), and GS treated (GS t1) ( Figure  9).  Total RNA was extracted from leaves from 27 plant samples using Trizol reagent (Invitrogen, Carlsbad, Calif, USA) in accordance with the manufacturer's protocol (Figure 9). In the total RNA, an additional treatment with RNase-free DNase I (Invitrogen) was performed to remove residual genomic DNA. The RNA samples from each respective biotype and treatment (untreated and glyphosate treated) were pooled to form three technical replicates per treatment totaling 12 samples, i.e., GR untreated (GR t0), GS untreated (GS t0), GR treated (GR t1), and GS treated (GS t1) (Figure 9).

cDNA Library Construction and Illumina Sequencing
Before cDNA construction, the RNA integrity number (RIN) was measured in Bioanalyzer (Agilent

De Novo Assembly and Functional Annotation
Raw data were filtered using FastQC (https://www.bioinformatics.babraham.ac.uk/) and trimmed adapter by Trimmomatic [51]. A single de novo assembly was generated using the reads from all treatments (12 libraries) with Trinity for further differential expression analysis [52]. For differential expression analysis, it is essential that all treatments be included in the transcriptome assembly [53,54]. Further, a single de novo assembly for each biotype (GR and GS) was generated for EPSPS sequence alignment. Trinity was executed using default settings in all assemblies and additional assembly analysis performed using the BioPython package [55]. The assembled transcripts generated by Trinity were aligned to the UniProt-trEMBL [56] database using Diamond [57], and only those with hits on plants (E-value threshold = 1 × 10 −10 ) were selected for further analysis.
The protein prediction was performed based on amino acid sequence using the Trinotate pipeline (https://trinotate.github.io/), which included the identification of protein families from the Pfam [57] and UniProt-SwissProt database, the identification of signal peptides and transmembrane proteins by SignalP and TMHMM, respectively [58], and rRNA prediction by RNAmmer [59]. The gene ontology functional categorization was performed by the BLASTx hits [60,61] with an E-value threshold of <10 −5 from the non-redundant database according to molecular function, biological process, and cellular component ontologies.

Differential Expression Analysis
For differential expression analysis sections, we termed the assembled contig as "gene". After data normalization, gene expression levels were calculated using transcript reads per million mapped reads (TPM) ≥1. For each treatment, gene-level expression was estimated using the Kallisto method [62] implemented in the Trinity accessory, which was used to generate an expression matrix. The matrix was processed in the edgeR mode [63] to perform the differential gene expression analysis and GO binning and enrichment was performed [63].
Differentially expressed genes (DEGs) were analyzed within each biotype with (t1) and without (t0) glyphosate treatment. DEGs were also filtered based on the false discovery rate (FDR) and p-value threshold at ≤0.001. After that, lists of all filtered DEGs were exported for each comparison and MA plots produced. In that case, DEGs with a log fold change ≥2 was considered up-regulated and ≤-2 as down-regulated. The up-and down-regulated genes of both biotypes for t1 and t0 treatments were categorized according to their GO functions using the same methods as those described above.

Quantitative Reverse Transcriptase Polymerase Chain Reaction (qRT-PCR) Analysis to Validate RNA-Seq Results
The same RNA samples from above were as converted into cDNA using RevertAid First Strand cDNA Synthesis (Thermo Fisher Scientific™), following the manufacturer's protocol. The cDNA was pooled according to each biotype and time of collection after glyphosate treatment, resulting in 9 cDNA samples, and the qRT-PCR experiments performed using three technical replicates in a QuantStudio6 TM (Applied Biosystems TM ). The reaction mixture of 15 µL contained 7.5 µL SYBR Green/ROX qPCR Master Mix (2×) (Thermo Fisher Scientific™), 1 µL of 1:3 diluted cDNA, 0.3 µL of each primer (forward and reverse) (10 µM), and the final volume was adjusted with appropriate amounts of RNAse free water. The QuantStudio6 TM thermal method was a multi-step of 10 min at 95 • C for pre-denaturation, 40 cycles of 15 s at 95 • C, and 1 min at 60 • C (2×). Data were collected during the extension step, and threshold-cycle (CT) values were calculated for each reaction using the Second Derivative Maximum method in the QuantStudio Real-Time PCR software (Applied Biosystems TM ).
For control reactions, no cDNA sample was added to the reaction mix of each studied primer. The sequence annotated as hexosyltransferase was used as an internal control (Contig: Trinity_DN_28781-F: CAATGGCCAGTCAAAACCAT; R: CCAGGCTCCATCCTATCGT A). The hexosyltransferase presented the best stability results among eleven candidate internal control genes tested, which included heat shock-protein and actin previously reported for C. bonariensis and C. canadensis, respectively [21,23].
Nineteen contigs with relative expression were selected from a transcriptome dataset based on the results of TPM of treated rather than untreated biotypes (GR-t1/t0; GS-t1/t0). In this case, in each respective replicate, TPM values >2 were considered as relative expression and ≤2 as absolute. Contigs with absolute expression were not considered for further analysis. Contigs with a low, medium, and high expression based on log fold change (logFC-t1/t0) results were selected from differential expression analysis to further analysis in qRT-PCR. The relative expression was calculated using the equation QR = 2 −(∆∆CT) [64]. The fold change results from transcriptome dataset were correlated to the time-course average of those from qRT-PCR using the Pearson model (r). The real-time results were subjected to ANOVA and averages compared by orthogonal contrasts at p ≤ 0.05 [64].

EPSPS Transcript Sequence Analysis
The sequences of the three EPSPS gene copies for C. bonariensis were obtained from GenBank (EPSPS1-EF200070; EPSPS2-EF200069; EPSPS3-EF200074) and mapped into the individual assembled transcriptomes of GR and GS using BLASTn command-line to identify their respective contigs. Trinity-assembled and GenBank sequences of EPSPS copies were translated to amino acids and aligned into BioEdit (http://www.mbio.ncsu.edu/BioEdit/) using the ClustalW multiple alignment functions with default settings. The alignment was performed to verify the occurrence of amino acid substitution at Threonine 102 and Proline 106, commonly altered reported regions to confer resistance to glyphosate [16,19]. Hence, the assembled EPSPS contigs were analyzed for transcription levels (TPM) to verify the overexpression occurrence in GR and GS in response to glyphosate treatment and copy number.

Selection of Similar Contig to Transport and Metabolic Glyphosate-Resistance Coding Gene
DEGs were selected by their UniProt/SwissProt assignment to gene families related to known roles in transportation and metabolic herbicide resistance [20,30]. We selected gene families of ABC transporters that were hypothesized to be involved in NTSR in C. canadensis [23], and cytochrome P450 (CYP450), glutathione (GST), glycosyltransferase (GT), related to be involved in herbicide metabolism [20,30,33]. Further, we selected DEGs annotated as antioxidant enzymes catalase (CAT), peroxidase (POD) and superoxide dismutase (SOD) reported with differential activities and to be involved in the reactive oxygen species (ROS) scavenging after glyphosate exposure in the same biotypes of C. bonariensis of the present study [28]. The antioxidant enzymes were also described as being involved in glyphosate resistance in Amaranthus palmeri and Ambrosia trifida [34,49].
The candidate genes were chosen according to the Venn diagram results from up-regulated genes in GR t0 (81 genes) and t1 (763 genes) and co-expressed in GR and GS plants (974 genes) in response to glyphosate treatment (Figure 1). The query DEGs were selected based on expression differences between GR (t1/t0) and GS (t1/t0) biotypes in response to glyphosate treatment (fold change ratio GR/GS). In this case, DEGs with a GR (t1/t0)/GS (t1/t0) ratio ≥2 were selected, except for the two SOD genes because they had a much higher basal expression in GR relative to GS plants. Finally, each gene sequence was searched against the UniProt/SwissProt database (www.uniprot.org/blast/) to assign a putative function.
We obtained from Peng et al. (2010) [23] the sequences of ABC transporters M10 and M11 that were hypothesized to be implicated in NTSR to glyphosate in horseweed. These sequences were mapped in our whole C. bonariensis transcriptome using BLASTn, and after, their expression levels were evaluated using the TPM results for GR and GS plants in response to glyphosate treatment.

Conclusions
The present transcriptomic study revealed 41 new candidate NTSR genes that are annotated to be related to transport and metabolism in herbicide resistance. Among these candidates, there were 19 ABC transporters, 10 CYP450, one glutathione, and five glycosyltransferases. In addition, we also report the transcription results of two genes coding for antioxidant enzyme catalase, peroxidase, and superoxide dismutase. Thus, all target genes from different groups with a transcriptional increase in the glyphosate-resistant biotype might be acting in association to confer resistance to glyphosate in C. bonariensis. Further, these results indicate that gene expression in C. bonariensis varies among gene groups and within the same group, between biotypes, in response to glyphosate treatment and is dependent on the time after herbicide exposure. The present transcriptome study is the first report that associates various CYP450, ABC transporters, GT, GST, and antioxidant enzyme differential gene expression to glyphosate resistance in C. bonariensis. The results of the present work will serve as a data resource for further studies on the molecular mechanisms of glyphosate resistance in C. bonariensis. Further studies will involve functional genomic analysis using the protoplasts and gene editing approaches.
Supplementary Materials: The following are available online at http://www.mdpi.com/2223-7747/8/6/157/s1, Figure S1. Length distribution of transcripts assembled from transcriptome libraries of C. bonariensis. Figure  S2. Sequence comparison to other organisms from the distribution of BLASTx hits (e-value < 1 × 10 −10 ) against the non-redundant protein database of the National Center for Biotechnology Information. Figure S3. Sequence comparison to other plants (hit ≥1%) from the distribution of BLASTx hits (e-value < 1 × 10 −10 ) against the non-redundant protein database of the National Center for Biotechnology Information (NCBI). Figure S4. Top ten Gene Ontology (GO) terms identified in the C. bonariensis transcriptome assembly summarized in three main categories: (A) biological process, (B) molecular function, and (C) cellular component. Figure S5. MA plot of differential expression analysis generated by EdgeR from transcriptome study performed in C. bonariensis glyphosate-resistant (GR) and -sensitive (GS) biotypes in response to glyphosate treatment. Dots above zero are up-regulated, and dots below are down-regulated. Red dots indicate significant expression at an adjusted p-value and false discovery rate (FDR) threshold set at ≤0.001, and log2FC ≥ 2 (up-regulated) or ≤log2FC (down-regulated). Plots for each contig its log2FC (fold change) (A, Y-axis) vs. its counts (mean of normalized counts) (M, X-axis). t0 = without glyphosate treatment; t1 = with glyphosate treatment. GR t1: RNA sampled at 24,96,192, and 288 h after treatment and pooled; GS t1: RNA sampled at 24, 96, and 192 h after treatment and pooled. Figure  S6. Correlation of transcriptomic and qRT-PCR (average of time-course) expression levels results of 19 genes of glyphosate-resistant (GR) and -sensitive (GS) C. bonariensis biotypes. Figure S7. Partial sequence alignment of the EPSPS transcripts and amino acid sequence assembled of glyphosate-resistant (GR) and -sensitive (GS), and a sensitive C. bonariensis sequence from GenBank. (A) EPSPS1 (GenBank-accession number EF200070); (B) EPSPS2 (GenBank-accession number EF200069); (C) EPSPS3 (GenBank-accession number EF200074). The red boxed amino acids show no substitution at positions Threonine 102 and Proline 106. (D) Transcriptome expression levels (transcript reads per million mapped reads-TPM) of the three EPSPS copies in GR and GS in response to glyphosate treatment and expression difference (Fold change). GR: glyphosate-resistant biotype; GS: glyphosate-sensitive biotype. t0: without glyphosate treatment; t1: with glyphosate treatment. Figure S8. C. bonariensis transcriptome expression analysis (transcript reads per million mapped reads-TPM) of the M10