The collection and storage of seeds from the resistant (R, AH18-origin) and susceptible (S, SD01-origin) A. aequalis populations were described in Guo et al. (2016). The sensitivity of these two populations was tested; a purified R population (AH18) homozygous for the ACCase and AHAS gene mutations and a confirmed S population (SD01) were both acquired (Guo et al., 2016). Prior to their planting, the seeds were germinated as described by Guo et al. (2015a). The pre-germinated seeds were sown in 12-cm-diameter, 11-cm-deep plastic pots (20 seeds per pot) containing moist loam soils. Next, the pots were randomly placed in a controlled greenhouse (natural light, 25/15°C, ~75% relative humidity), where they were watered every other day. For each pot, the weed seedlings that emerged were thinned to 15 that were evenly sized at the two- to three-leaf stage.

Whole-plant dose-response experiments were carried out to determine the GR₅₀ values of the purified R (AH18) and confirmed S (SD01) populations to mesosulfuron-methyl in the absence and presence of malathion. The malathion used in this experiment was formulated in a mixture of a commercial solvent and an emulsifier. The mesosulfuron-methyl was applied by using a commercial formulation—30 g L⁻¹ oil-miscible flowable concentrate (Bayer, Hangzhou, China)—diluted with distilled water.

When the A. aequalis seedlings reached the three- to four-leaf stage they were ready for use in the whole-plant dose-response experiments. These experiments were done twice in a completely randomized design that had three replicates. In each experimental repeat there was a total of 84 pots used (2 malathion rates × 2 populations × 3 replicates × 7 treatments). Malathion was applied as an active ingredient (ai), at a rate of 0 or 1,000 g ai ha⁻¹, at 1 h before the herbicides were applied (Preston et al., 1996). Mesosulfuron-methyl application rates of 0, 3, 9, 27, 81, 243, and 729 g ai ha⁻¹ were used to treat the R plants, while rates of 0, 0.04, 0.11, 0.33, 1, 3, and 9 g ai ha⁻¹ were used on the S plants. These herbicide application rates were adopted to represent the GR₅₀ value for each population in the absence of CytP450 inhibitors. The control plants were treated with the mixture (commercial solvent and emulsifier) used to formulate the malathion. All the herbicides were applied by using an air-compressed, moving-nozzle cabinet sprayer which delivered 450 L ha⁻¹ of water at a pressure of 0.275 MPa, via a flat-fan nozzle positioned 50 cm above the foliage. After being treated, the experimental plants were returned to the controlled greenhouse. Three weeks later, their aboveground shoots were harvested and oven-dried at 80°C for 72 h and their dry weight data were recorded.

Plant data obtained from the mesosulfuron-methyl dose-responses were expressed in relative terms, as a percentage of the control. The datasets from the replicate experiments were analyzed by ANOVA that used the general linear model procedure in SPSS v.19.0.0 (IBM Corporation, Armonk, USA). Because the variance among the replicated experiments was not significant (p > 0.05), the data were pooled and fitted to a non-linear regression analysis in SigmaPlot v.12.0 (Systat Software, Inc., San Jose, USA). The GR₅₀ was determined by using a four-parameter log-logistic equation (Seefeldt et al., 1995). The fitted model was as follows:

where C is the lower limit of the response, D is the upper limit of the response, X is the herbicide application rate, and b is the slope of the curve at the GR₅₀.

To estimate the resistance levels, the resistance index (RI) was calculated as follows:

Ten seeds were randomly selected from the purified R (AH18) and confirmed S (SD01) populations, respectively, for germination and onward cultivation to maturity. Mature seeds (AH18-F1/SD01-F1) were collected from each individual plant. To study metabolic resistance, one seed from each individual-plant progeny was randomly selected to make up a 10-seed combination. The experimental design included three biological replicates of R and S for the untreated control (CK), water-only control (WCK), and mesosulfuron-methyl (T; 1 × labeled field rate, 9 g ai ha⁻¹) treatments, for a total of eighteen 10-seed combinations of the R and S biotype.

The selected 180 seeds were germinated and cultivated to the three- to four-leaf stage under the above-mentioned experimental conditions. Leaf samples for the untreated control were collected at time point 0 (R_CK/S_CK), while those from the water-only control (R_WCK/S_WCK) and the mesosulfuron-methyl treatments (R_T/S_T) were collected 24 h after the treatments were applied. Each sample consisted of a two-leaf mixture: a 4-cm-long, newest emerging leaf, and a 3-cm-long, first fully expanded leaf (both taken from each of the 10 selected plants). A total of 18 samples were thus collected (3 biological replicates × 3 treatments × 2 populations). All of these samples were frozen immediately in liquid nitrogen to avoid RNA degradation and/or induction of plant responses to wounding; they were then pulverized for the total RNA extraction. The remaining seeds (AH18-F1/SD01-F1) from each individual-plant progeny were used in parallel to those for RNA-Seq—they were likewise planted, cultivated to the three- to four-leaf stage, and used as plant materials in the subsequent validation experiment by quantitative real-time reverse transcription polymerase chain reaction (qRT-PCR).

The total RNAs from the 18 samples were extracted with Transzol Up (TransGen Biotech, Beijing, China) according to the manufacturer's protocol and they were treated with DNase I (Takara, Beijing, China). Total RNA degradation and contamination were monitored on 1% agarose gels. RNA purity was checked using a NanoPhotometer® spectrophotometer (IMPLEN, West Lake Village, USA), and the RNA concentration was measured using a Qubit RNA Assay Kit and a Qubit 2.0 Fluorometer (Life Technologies, Gaithersburg, USA). The RNA integrity was assessed with a RNA Nano 6000 Assay Kit as part of the Agilent Bio analyzer 2,100 system (Agilent Technologies, Palo Alto, USA). The highest-quality RNA samples were selected for preparing the cDNA library.

The cDNA library construction and the Illumina sequencing were both performed by Annoroad Gene Technology Co., Ltd (Beijing, China). Briefly, the sequencing libraries were generated via a NEBNext Ultra™ RNA Library Prep Kit for Illumina (NEB, Ipswich, USA) following the manufacturer's instructions. A total of 3 μg of RNA per sample was used as the input material for the RNA-sample preparations. Poly (A) mRNA was purified from total RNA by using poly-T oligo-attached magnetic beads, and then fragmented into short sequences. These cleaved RNA fragments were transcribed into first-strand cDNA by a Moloney Murine Leukemia Virus (M-MuLV) Reverse Transcriptase (RNase H-) and by random hexamer primers; this step was followed by a second-strand cDNA synthesis in a reaction buffer which included dNTPs, DNA polymerase I, and RNase H. After end-repair and adapter-ligation, the products were amplified through PCR and purified by using a QIAquick PCR Purification Kit (QIAGEN, Valencia, USA) to construct the cDNA library. The cDNA libraries were sequenced on the Illumina Hiseq™ 4000 platform (Illumina Inc., San Diego, USA) and 150-bp paired-end reads were generated.

Raw image data files from the Illumina Hiseq™ 4000 were transformed to raw reads by the Consensus Assessment of Sequence and Variation (CASAVA, v.1.8.2) base recognition program, and then stored in FASTQ-formatted files. Raw reads containing the adapter poly-N (>5% of the unknown sequences designated as “N”) and the low-quality reads (>15% bases with a quality value <20) were filtered out by in-house Perl scripts to obtain high-quality, clean reads. All downstream analyses were based only on these clean reads of high quality.

A de novo transcriptome assembly was accomplished in Trinity v.20140717 (Grabherr et al., 2011) under its default parameter values, by which the transcripts and unigenes—the longest transcript of a set of transcripts that appears to originate from the same transcription locus—were obtained. To assess the accuracy of assembly, all of the unigenes were matched to the sequences of an annual model plant, false-brome (Brachypodium distachyon [L.] Beauv.), using the BLAST-Like Alignment Tool (BLAST) with E ≤ 1e-20. The assembled transcripts were translated into all six possible open reading frames (ORFs) through TransDecoder v.20140717, and their proper translation was defined as the one that gave the longest amino acid sequence (Saha et al., 2002).

To obtain comprehensive, functional information regarding the sequences, gene function was annotated through publicly accessible databases—NCBI non-redundant protein sequences (Nr), NCBI non-redundant nucleotide sequences (Nt), Universal Protein Resource (UniProt), protein family (PFAM), and orthologous groups and functional annotation (eggNOG)—by means of local BLAST programs (NCBI, NIH, Bethesda, USA) which used a significance threshold of E < 1e-5. Meanwhile, the software programs HMMER v.3.1b1 (Finn et al., 2011), signalIP v.4.1 (Petersen et al., 2011), and TMHMM v.2.0 (Krogh et al., 2001) in Trinotate were run to identify the protein domains, to predict the signal peptides, and to locate the trans membrane regions, respectively.

Next, the gene ontology (GO, http://www.geneontology.org/) annotation for the unigenes was implemented with Blast2go (Conesa et al., 2005), set to an E-value cutoff = 1e-5. The GO functional classifications were performed using WeGO software (http://wego.genomics.org.cn/cgi-bin/wego/index.pl) (Ye et al., 2006). The pathway assignments were carried out by sequence searches against the Kyoto Encyclopedia of Genes and Genomes database (KEGG, http://www.kegg.jp/), which used the KEGG Automatic Annotation Server (KAAS, http://www.genome.jp/kegg/kaas/) set to an E-value threshold < 1e-10.

The assembled transcriptome was used as the reference database, and the clean reads were mapped back onto the reference transcriptome by Bowtie 2 v.2.2.3 (Langmead and Salzberg, 2012). The read count for each gene in each sample was obtained from the mapping results by RSEM (RNA-Seq by Expectation Maximization) (Li and Dewey, 2011); then, the mapped read counts for each transcript were normalized, to eliminate the effects of the sequencing depth and gene length on the gene expression levels, by RPKM (Reads Per Kilobase Millon Mapped Reads) (Mortazavi et al., 2008). In this way, the gene expression levels were determined for each sample and could be directly compared. Next, the differences in the expression abundance of each gene per transcript between sample pairs were calculated by using DESeq 2 v.1.4.5 (Love et al., 2014). Each resulting p-value was adjusted to a q-value, following the Benjamini-Hochberg procedure for controlling the false discovery rate (Storey and Tibshirani, 2003). Those genes with a q < 0.05 and |log2(fold change)| ≥ 1 were identified as DEGs.

Expression differences were compared between the untreated R and untreated S at time point 0 (R_CK and S_CK), between R and S at 24 h after the water control treatment (R_WCK and S_WCK), and between R and S at 24 h after the mesosulfuron-methyl treatment (R_T and S_T). Expression differences were also compared within R and S among the untreated (CK), water-only control (WCK), and mesosulfuron-methyl treatments (T). The DEGs were later analyzed by GO and a KEGG enrichment analysis. Both analyses were implemented with a hyper geometric test (Young et al., 2010), in which the p-value was adjusted for multiple comparisons as the q-value, and the background data are the genes in the whole genome. Those GO or KEGG terms with a q < 0.05 were considered as significantly enriched.

To select candidate contigs for metabolic resistance, three factors were simultaneously considered: their statistical significance, their magnitude of expression differences, and their annotations related to known herbicide metabolism genes and signaling functions. The expression level of each contig was first measured by qRT-PCR in the original RNA samples that were used in the RNA-Seq experiment. For the normalization of gene expression and to ensure the reliability of the qRT-PCR data, the ubiquitin (UBQ, JN599096.1) and glyceraldehyde-3-phosphate dehydrogenase genes (GAPDH, JN599100.1) were used as the internal control genes (Petit et al., 2012). Gene-specific primers were designed according to the sequences of the selected unigenes in Primer Premier v.5.0 (Premier Biosoft International, Palo Alto, USA). The sequences of these primers are listed in Supplementary Table S4.

The qRT-PCR was performed in 96-well plates on the Bio-Rad CFX96 Touch™ real-time PCR system (Bio-Rad, Richmond, USA), which used the TransStart® Top Green qPCR SuperMix (TransGen Biotech, Beijing, China). Each reaction was conducted in a 20-μL mixture—consisting of 10 μL of 2 × TransStart® Top Green qPCR SuperMix, 0.4 μL of forward primer, 0.4 μL of reverse primer, 0.4 μL of cDNA, and 8.8 μL of RNase-free ddH₂O—with four replicates per cDNA sample. The qRT-PCR programs consisted of a 2-min incubation phase at 95°C, followed by 45 cycles, with each cycle consisting of a 10-s incubation at 95°C, a 15-s incubation at 60°C, and a 20-s incubation at 72°C. At the end of the amplification cycle, melting curve analyses were performed to confirm the specificity of qRT-PCR primers. The amplification efficiencies were derived from a standard curve generated by a 10-fold serial dilution points of cDNA (Pfaffl, 2001). A similar amplification efficacy (87.6–99.7%) for the target and internal control genes was observed. The fold-change in gene expression (as 2^−ΔCt) was calculated by using the comparative Ct method (Schmittgen and Livak, 2008), relative to the S samples, where ΔCt = [Ct target gene-Ct mean of the two internal control genes]; three biological replicates were obtained. Statistical analysis of the qRT-PCR data was conducted using the Student's t-test (p < 0.05) procedure of SPSS software.

Furthermore, the expression patterns of those contigs, whose relative expression data provided by qRT-PCR agreed with the RNA-Seq-detected profiles (i.e., the contigs whose expression levels were significantly higher in the R than in the S samples in both the RNA-Seq and qRT-PCR experiments), were also examined in the parallel plant materials—for which the seed germination, plant cultivation, and total RNA extraction followed the above-mentioned methods. Their cDNA was later obtained, and qRT-PCR was carried out using the same primers and procedure as before.

The whole-plant dose-response experiments demonstrated that the purified R population (AH18) had evolved a high level resistance (RI = 31.15-fold) to mesosulfuron-methyl (Table 1); this was consistent with our previous report (Guo et al., 2016). When malathion alone was applied at 1,000 g ai ha⁻¹, there were no visual effects on A. aequalis seedling growth, in either the R or S population (Figure 1). Importantly, malathion greatly reduced the resistance level of the R population when it was used before the mesosulfuron-methyl treatment. Under a malathion pretreatment of 1,000 g ai ha⁻¹, the GR₅₀ value of mesosulfuron-methyl decreased by 41% for the R population (Table 1, Figure 2). By contrast, malathion did not lead to any obvious increase in the susceptibility of the S population to mesosulfuron-methyl, as its GR₅₀ decreased little (Table 1). Malathion, the P450 inhibitor, has long been used as an indicator of CytP450 involvement in metabolic resistance to AHAS-inhibiting herbicides. The observation that, in our present study, malathion caused a significant reduction in the resistance level when it was applied with the herbicide indicated the involvement of at least one CytP450 gene in the resistance mechanisms operating in the R biotype.

Using the Illumina sequencing technology, a total of 975,586,352 raw sequencing reads were generated from the 18 RNA libraries (S_CK, R_CK, S_WCK, R_WCK, S_T, and W_T; each had three biological replicates). After the quality control and data cleaning, there were 831,846,736 clean reads left over, ranging from 49,270,938 to 60,784,678 per sample, for use in the de novo assembly (Table 2, Supplementary Table S1). The clean reads assembled to 179,989 transcripts, with an average length of 1,233 bp. Figure 3 shows the length distribution of these transcripts. Of the 831,846,736 clean reads sequenced, we were able to map the great majority of them (93.02%) back onto the transcripts using Bowtie 2. Among these transcripts, 43,290 unigenes >200 bp and 14,352 unigenes >500 bp, with an overall mean length of 944 bp and an N50 length of 1,857 bp, were obtained by using the longest transcript at each locus of each gene. The length distribution of these unigenes is also shown in Figure 3. B. distachyon is now intensively utilized as a model grass species in various biological studies (Draper et al., 2001). By using the BLAST search, 34.40% of the unigenes were successfully matched to the sequences of B. distachyon, at an average accuracy of 84.18%.

The numbers of unigenes annotated by the databases are summarized in Table 3. Of the 95,479 assembled unigenes, over half were annotated in at least one of the seven databases, whereas <5% had annotations in all seven databases, and the greatest sequence similarity was obtained from the Nr database (Table 3). Based on the Nr database annotations, the unigene sequences of our study species, A. aequalis, were most similar to the gene sequences from B. distachyon, Hordeum vulgare subsp. vulgare, Aegilops tauschii Coss., Triticum Urartu L., Oryza sativa L. Japonica, and Triticum aestivum L. (Figure 4).

To characterize the functional classifications of the annotated unigenes, GO and KEGG analyses were performed to access the distributions of functional categories. A total of 29,889 (29.29%) unigenes were annotated by the GO database; these were classified into 67 functional subgroups consisting of 23 subgroups as “biological process” (BP), 22 as “cellular component” (CC), and 22 as “molecular function” (MF) (Figure 5). In the BP category, cellular process (18,218, 60.95%) and metabolic process (17,088, 57.17%) represented the two major contributors. In the molecular function category, binding (18,469, 61.79%) and catalytic activity (16,732, 55.98%) represented the major contributors. In addition, 14,246 unigenes were assigned to 32 KEGG pathways; these unigenes were primarily involved in “carbohydrate metabolism” (1,722, 12.09%), “translation” (1,244, 8.73%), “signal transduction” (1,228, 8.62%), and “amino acid metabolism” (1,200, 8.42%)” (Figure 6).

Constitutive differential expression (i.e., q < 0.05 and |log₂(fold change)|≥1) was apparent between the untreated R and S samples, such that 5,840 contigs were constitutively up-regulated in R, and 5,648 contigs were up-regulated in S (Figure 7A). Differential expression between R and S in all three treatments was evident for 5,654 contigs (Figure 7B). When comparing the mesosulfuron-methyl treatment with the water-only control, a total of 5,549 contigs were up-regulated in R, and 5,649 contigs were up-regulated in S. However, considering only those contigs up-regulated by the mesosulfuron-methyl treatment revealed an overlap of 409 contigs between the R and S populations, including contigs putatively annotated as 10 CYPs, 1 GST, 4 GTs, and 2 ABC transporters (Supplementary Table S2). DEGs between the mesosulfuron-methyl treatments and the untreated control were also identified: there were 5,837 contigs up-regulated and 4,492 contigs down-regulated in R_T relative to R_CK, as well as 8,561 contigs up-regulated and 6,859 contigs down-regulated in S_T relative to S_CK (Figure 7A).

To further characterize the function of the DEGs, we conducted GO enrichment analyses. In total, there were 11,915 DEGs found between the R_T and S_T samples; these were classified into 57 functional groups, namely 22 groups as BP, 17 groups as CC, and 18 groups as MF (Supplementary Figure S1). Among these DEGs, the GO terms of “cellular process” (1587, 13.32%), “metabolic process” (1200, 12.59%), and “single-organism process” (1202, 10.09%)” in the BP category, “cell part” (1912, 16.05%)” in the CC category, and those of “binding” (1691, 14.19%)” and “catalytic” (1691, 13.00%) in the MF category, were all significantly enriched between the R and S population samples.

We also conducted a detailed KEGG enrichment analysis of the DEGs between the R and S populations for each of the three treatments (Supplementary Figures S2–S4). The DEGs between the mesosulfuron-methyl-treated R and S samples were significantly enriched in the 12 pathways (Table 4). Specifically, 17 up-regulated genes were enriched in the “metabolism of xenobiotics by CytP450” pathway, and 16 up-regulated genes were enriched in the “drug metabolizing-cytochrome CytP450” pathway. Notably, the genes involved in the “metabolism by cytochrome CytP450” changed significantly in the KEGG pathway analysis; this result suggested that the CytP450 genes could play a key role in the metabolic resistance to mesosulfuron-methyl in A. aequalis.

Considering the contigs differentially expressed between the R and S populations, those having predicted annotations related to metabolism and signaling pathways are more likely to be involved in conferring metabolic mesosulfuron-methyl resistance (Gaines et al., 2014). Therefore, we selected those overlapping genes up-regulated in all three comparative groups—i.e., R_T relative to S_T, R_T relative to R_WCK, and R_T relative to R_CK—which had such predicted annotations (Supplementary Table S2). We also selected those up-regulated genes annotated as related to metabolism that had a significantly different expression between R and S in all three treatments (Supplementary Table S3). In particular, the overlapping genes having |log2(fold change)| values in R_T relative to S_T exceeding those of the other two comparable groups (i.e., R_WCK vs. S_WCK and R_CK vs. S_CK) were identified (Supplementary Table S3). In addition, contigs were also chosen that had predicted annotations related to oxidation (e.g., oxidase and peroxidase), hydrolysis (e.g., esterase and hydrolase), and plant stress physiology (e.g., response to plant hormone stimulus) in the above-mentioned comparative groups, along with those contigs with the largest fold-change expression differences but of an unknown function (i.e., no assigned annotation). Based on these criteria, a total of 31 contigs were selected as candidate genes that might confer mesosulfuron-methyl resistance to A. aequalis (Table 5). Of these contigs, eight were annotated to the CytP450 family, five to the GST family, five to GT family, and two to the ABC transporter family, with the remainder (11) annotated to other families.

The expression levels of the 31 candidate contigs were first validated by qRT-PCR using the original RNA samples (i.e., R_T and S_T). The relative expression data provided by the qRT-PCR showed that 24 of the 31 candidate contigs exhibited significantly higher expression levels in the R (AH18) than in the S (SD01) population samples (Table 5).

The expression patterns of the 24 candidate contigs were evaluated further using the parallel A. aequalis plant materials. These parallel samples had a genetic background identical to the RNA-Seq samples. Based on the results of the RNA-Seq and the two rounds of the qRT-PCR validations, a total of 17 candidate contigs were expressed significantly higher in the R (AH18) than in the S (SD01) population samples (Table 5, Figure 8, Supplementary Figure S5). These 17 contigs possibly encoded proteins with homology to four CytP450s (CYP94A1, CYP94A2, CYP71A4, and CYP734A6), two GSTs (GST-T3 and GST-F1), two GTs (GT75D1 and GT83A1), two ABC transporters (ABCC8 and ABCB11), three oxidases, two hydrolases, and two kinases. Although some variation among individual plants was observed, the average expression levels of the candidate contigs provided by qRT-PCR were consistent with the RNA-Seq expression data, thus revealing little deviation in the fold change (Figure 8).

In our study, four CytP450 genes (CYP94A1, CYP94A2, CYP71A4, and CYP734A6) were identified as major candidates for metabolic resistance. The biochemical role of CytP450-mediated herbicide metabolism has been well-established in resistant L. rigidum populations, and several distinct CytP450 isoforms were shown to be involved in metabolic resistance to ACCase- and AHAS-inhibiting herbicides (Christopher et al., 1994; Preston et al., 1996; Preston and Powles, 1998; Busi et al., 2011; Duhoux and Délye, 2013; Gaines et al., 2014). Several CytP450 genes, such as CYP71AK2 and CYP72A254, were reported to show a high expression level in a multiple-herbicide-resistant E. phyllopogon population (Iwakami et al., 2014b), while CYP94A1 was implicated in the process of plant defense against chemical injury in the common vetch (Vicia sativa L.) (Benveniste et al., 2005). The CYP734A family typically functions in the brassinosteroid catabolism of plants, and brassinosteroids are known to protect plant crops against toxicity from herbicides (Ohnishi et al., 2006; Xia et al., 2009). Transgenic Arabidopsis (Arabidopsis thaliana) expressing either of the two CytP450 genes, CYP81A12 or CYP81A21, confers resistance to AHAS inhibitors bensulfuron and penoxsulam (Iwakami et al., 2014a). Moreover, many CytP450s conferring herbicide resistance have been identified in several major crop species, such as CYP81A6 in rice (O. sativa) (Pan et al., 2006), CYP71C6v1 in wheat (T. aestivum) (Wen-Sheng et al., 2006), CYP71A10 in soybean (Glycine max [Linn.] Merr.) (Siminszky et al., 1999), and CYP76B1 in Jerusalem artichoke (Helianthus tuberosus L.) (Didierjean et al., 2002). Importantly, it has been demonstrated that the purified recombinant protein expressed by the T. aestivum CYP71C6v1 gene could specifically metabolize three herbicides, namely chlorsulfuron, triasulfuron, and metsulfuron-methyl, via the phenyl ring hydroxylase (Wen-Sheng et al., 2006; Xiang et al., 2006). The CytP450 genes that we identified in the R A. aequalis population may also confer potential resistance to multiple herbicides, but this hypothesis will require empirical testing for its validation.

The GSTs and GTs are two other key families of enzymes involved in herbicide detoxification. The biochemical role of GST and GT in the NTSR to herbicides is now well-established (Brazier-Hicks and Edwards, 2005; Cummins et al., 2011, 2013). Prior studies indicated that GSTs and GTs were associated with NTSR to the ACCase- and AHAS-inhibiting herbicides in L. rigidum (Gaines et al., 2014; Duhoux et al., 2015). In the present study, two GST genes (GST-T3 and GST-F1) and two GT genes (GT83A1 and GT75D1) were identified as associated with the metabolic resistance of A. aequalis. The predicted GST-F1 gene shows a high sequence homology with A. myosuroides (86%, UniProt accession No. Q9ZS17). Another study showed GST-F1 plays a central role in the multiple-herbicide resistance of A. myosuroides, and that the expression of its orthologous GST-F1 in a multiple-resistant L. rigidum population also showed enhancement (Cummins et al., 2013). Both the Tau (U) and Phi (F) classes of GSTs play roles in the herbicide detoxification performed by plants for several compounds (Cummins et al., 2011). Some GST genes conferring herbicide tolerance were also identified in several crop species, such as TaGSTU4 in wheat (T. aestivum) (Thom et al., 2002), GST I, and GST-27 in maize (Zea mays L.) (Karavangeli et al., 2005), and GmGSTU21 and GmhGS in soybean (G. max) (Skipsey et al., 2005). Additionally, the GT75D1 found in the present study shares a 73% sequence identity with the predicted 5-O-glucosyltransferase gene of B. distachyon (GenBank accession No. XP_003571494.1). Coincidentally, an A. myosuroides population, one resistant to multiple classes of herbicides, due in part to the elevated expression of CytP450s and GSTs active in herbicide detoxification, displayed higher O-glucosyltransferase activity than did the herbicide-sensitive populations of A. myosuroides (Brazier et al., 2002). Similarly, GT gene activities, for example GTA triggered by the application of AHAS-inhibiting herbicides, were observed in A. thaliana (Manabe et al., 2007), and two GT genes were shown to glycosylate the primary major bentazone metabolite, 6-hydroxybentazone, in the herbicide-tolerant G. max (Leah et al., 1992). Together, these findings suggest to us that the GST-F1 and GT75D1 genes up-regulated in the R A. aequali population may be associated with its mesosulfuron-methyl resistance. However, the respective roles of GST-T3 and GT83A1 remain unknown, and more investigative work is needed to elucidate their putative roles in herbicide resistance.

Plant ABC transporters, one of the most diverse gene families, are implicated in the detoxification of xenobiotics, which include herbicides (Rea, 2007). However, in contrast to the CytP450, GST, and GT gene families which confer herbicide resistance through metabolism, the ABC transporters detoxify herbicides by sequestering the herbicides and their metabolites (Yuan et al., 2007). In the present study, the two contigs c34002_g1 and c39076_g1 were annotated to ABCC8 and ABCB11, respectively, and so they may be associated with the mesosulfuron-methyl resistance. Similar to our result, an ABCC1, which could play an important role in the metabolic resistance to tribenuron-methyl, was identified in the resistant D. sophia population (Yang et al., 2016); an ABCB11 was also identified and found up-regulated in a resistant B. syzigachne population having NTSR to fenoxaprop-P-ethyl (Pan et al., 2016). Indeed, there are considerable lines of evidence for the involvement of ABC transporters in protecting plant tissues from toxin-induced damage. Two ABCC transporters, AtABCC1 and AtABCC2, were reported to mediate the tolerance of arsenic and arsenic-based herbicides in A. thaliana (Song et al., 2010). An ABC transporter, AtOPT6, showed an association with herbicide resistance in A. thaliana (Cagnac et al., 2004). Finally, the overexpression of AtPgp1, a multidrug resistance family member, and that of its garden pea homolog psNTP9, was able to confer resistance to multiple herbicides in A. thaliana (Windsor et al., 2003). Therefore, it is likely that ABCC8 and ABCB11 participate in driving metabolic resistance to mesosulfuron-methyl.

In this study, we found five further genes possibly associated with mesosulfuron-methyl metabolic resistance in A. aequalis. These genes are annotated in three major enzyme families functioning in the herbicide degradation process: hydrolases, oxidases, and peroxidases. Hydrolases can cleave the herbicide molecule, after which the cleaved molecule is transformed into a more hydrophilic metabolite by the oxidases and peroxidases (Délye, 2013). In addition, the increased expression of oxidases and peroxidases could effectively protect plant cells against the oxidative damage caused by a herbicide's action, which may also contribute to herbicide resistance (Délye, 2013).

Additionally, we selected two up-regulated contigs potentially encoding proteins with homology to kinase involved in response to the stimulus of salicylic acid (a plant hormone). Salicylic acid is a signaling molecule involved with local plant defense reactions at infection sites and the induction of systemic resistance (Durner et al., 1997). Plants can withstand various biotic and abiotic stresses in variable environments by having evolved complicated networks of stress detection, signaling, and response pathways which trigger both general and specific responses, often with an adjustment of the response over time (Délye, 2013). As a powerful abiotic stress, herbicides are known to trigger some of these pathways depending on the specific type of study, the herbicide, and the plant species involved (Das et al., 2010; Unver et al., 2010). NTSR is clearly part of an abiotic stress response, but the extent of this involvement varies among genotypes. Importantly, hormones participate in optimizing plant responses to abiotic stresses (Peleg and Blumwald, 2011), and yet there are few reports of such responses involving NTSR. The up-regulated genes annotated as hormone-related in this study may well be related to the mesosulfuron-methyl-induced stress response pathways. Considering the complexity of most plant defense mechanisms, further rigorous and comprehensive studies are required.

In this study, the RNA-Seq experimental design consisted of three treatments: an untreated control, a water-only control, and a mesosulfuron-methyl treatment. The objectives of using two controls were to measure the effects on transcription due specifically to the mesosulfuron-methyl treatment, and to control for any transcriptional effects from the spraying of plants with water (the solvent for mesosulfuron-methyl). According to our results, more DEGs were evidently observed due to mesosulfuron-methyl treatment relative to the untreated control than due to mesosulfuron-methyl treatment relative to the water-only control (Figure 7A). Hence, the DEGs triggered specially by mesosulfuron-methyl may play a more important role in the NTSR of A. aequalis. An earlier RNA-Seq experiment likewise established two controls—an untreated control and adjuvant-only control—to eliminate the potential effects on transcription due to the adjuvant (Gaines et al., 2014). By contrast, many other published reports pay more attention to the comparison between the untreated control and the herbicide treatment, which neglects the possible effects of the solvent or adjuvant on the RNA-Seq results (e.g., Yang et al., 2013, 2016; An et al., 2014; Chen et al., 2016).

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.