Whole Transcriptome Analysis Resulted in the Identi cation of Chinese Sprangletop (Leptochloa Chinensis) Genes Involved in Cyhalofop-Butyl Resistance

Ke Chen Longping Branch, Graduate School of Hunan University Yajun Peng Plant Protection Institute, Hunan Academy of Agricultural Sciences Liang Zhang Longping Branch, Graduate School of Hunan University Long Wang Hunan University Donghai Mao Institute of Subtropical Agriculture Zhenghong Zhao Longping Branch, Graduate School of Hunan University Lianyang Bai Agricultural Biotechnology Research Institute, Hunan Academy of Agricultural Sciences Lifeng Wang (  ifwang@hunaas.cn ) Agricultural Biotechnology Research Institute, Hunan Academy of Agricultural Sciences


Results
In this study, we utilized a resistant (LC18002) and a sensitive (LC17041) L. chinensis populations previously identi ed in our laboratory, which were divided into nine different treatment groups. We then employed whole transcriptome analysis to identify candidate genes which may be involved in cyhalofopbutyl tolerance. This analysis resulted in the identi cation of eight possible candidate genes, including six cytochrome P450 monooxygenase genes and two ATP-binding cassette transporter genes. We then carried out a phylogenetic analysis to identify homologs of the differentially expressed P450 genes. This phylogenetic analysis indicated that every genes have close homologs in pattern species, some of which have been implicated in non-target site resistance (NTSR).

Conclusions
This study is the rst to use whole transcriptome analysis to identify herbicide non-target resistance genes in L. chinensis. The differentially expressed genes represent promising targets for better understanding herbicide tolerance in L. chinensis. The eight genes belonging to classes already associated in herbicide tolerance may play important roles in the metabolic resistance of L. chinensis to cyhalofop-butyl, although the exact mechanisms require further study. Cyhalofop-butyl is an acetyl-coenzyme A carboxylase (ACCase) inhibitor herbicide used in paddy elds to control L. chinensis. It is the most commonly utilized herbicide and possesses high activity against L. chinensis when deployed just after rice emergence. Due to extensive and continuous use of cyhalofop-butyl, tolerance has developed in several L. chinensis populations [1][2][3]. Due to the high level of cyhalofop-butyl tolerance in recent years, it is highly important to get a better understanding of the mechanism of this tolerance in order to formulate feasible prevention and control measures.
Herbicide resistance refers to the ability of weeds to survive a typically lethal dose of herbicide with stable heredity, rather than a temporary phenotypic response to environmental conditions [4]. The mechanisms of weed resistance can be broadly divided into target site resistance (TSR) and non-target site resistance (NTSR) [5]. At present, research on the mechanisms of weed resistance has mainly focused on TSR. This mechanism is characterized by a change in the conformation of herbicide target enzyme in weeds that prevents herbicide binding, or an increase in the expression of the target gene which overcomes the inhibitory effects of the herbicide [6]. In addition to TSR, other resistance mechanisms are generally classi ed as NTSR. NTSR is comprised of a diverse set of mechanisms and can be related to generic plant stress response or detoxi cation of herbicides [7].
Cytochrome P450 monooxygenases (P450s) represent the largest superfamily of proteases and are involved in a host of different biological mechanisms. Hofer et al. found that the Arabidopsis thaliana genes CYP76C1, CYP76C2, and CYP76C4 were involved in the metabolism of benzoyl urea herbicides and played a role in the detoxi cation of phenylurea herbicides [8]. Saika et al. found that CYP72A31 and CYP81A6 were involved in the detoxi cation of bensulfuron-methyl herbicides in rice and A. thaliana [9][10][11]. In addition, P450 enzyme activity levels have been shown to correlate with NTSR for several different resistant weeds. During the herbicide detoxi cation process, detoxi cation products begin to accumulate in cells, eventually leading to a decrease in the activity of detoxi cation enzymes. Therefore, these detoxi cation products must be transported out of the cell in order for detoxi cation to continue. This is typically carried out by ATP-binding cassette (ABC) transporters. ABC transporters family is also one of the largest and most versatile protein families found in living organisms. It has been proposed by Hart et al. that the mechanism of plant resistance to paraquat toxicity might be either the transfer of paraquat to plant cell vacuoles by ABC transporters or the enhancement of antioxidant enzyme activity [12] Yang et al. found that the NTSR of ixweed [Descurainia sophia L.] to tribenuron-methyl is likely the result of metabolic resistance mediated by P450s and the movement of metabolites mediated by ABC transporters [13].
It has been reported that the population of L. chinensis has developed a high level of resistance to cyhalofop-butyl in Shanghai, Zhejiang, Jiangsu, Hunan, and other regions of China [1][2][3]. Despite the importance of understanding this resistance, research has thus far mainly focused by TSR mechanisms, with very little work aimed at understanding potential NTSR.
In previous studies, our laboratory identi ed several populations of L. chinensis which were resistant to cyhalofop-butyl, in addition to a susceptible population called LC17041. Additionally, we have previously shown that application of the P450 inhibitor malathion resulted in increased sensitivity of a previously resistant population (LC18002) to cyhalofop-butyl and no mutation of target ezyme amino acid was found in this populatio. A better understanding of the L. chinensis herbicide resistance mechanism will likely result in the development of improved weed control measures, which in turn will increase crop yields. In this study, transcriptomic data from LC18002 and LC17041 were generated to identify candidate genes related to cyhalofop-butyl resistance in LC18002 [1]. This work represents the rst whole transcriptomic study aimed at identifying the genes responsible for NTSR in L. chinensis.

Results
Previous studies in our laboratory have shown that LC18002 does not contain mutations in the direct target of cyhalofop-butyl, implying that its resistance is most likely operating via an NTSR mechanism. Metabolism-based resistance to herbicides that inhibit ACCases has been reported extensively, with instances increasing more in recent years [14]. RNA-seq has been successfully used to identify genes involved in metabolic resistance to acetolactate synthase (ALS) herbicides in two grass weed species, including Wimmera ryegrass (Lolium rigidum Gaud.) [15] and slender foxtail (Alopecurus myosuroides) [16]. However, no study has yet been conducted at the whole transcriptome level to understand potential NTSR mechanisms involved in L. chinensis resistance. Identifying genes involved in NTSR is crucial to better understand the evolution of metabolic resistance, which may lead to improved weed management strategies.
Transcriptome sequencing and screening RNA obtained from the different L. chinensis samples was paired-end sequenced using an Illumina HiSeq2500. The raw reads were then trimmed to remove adaptors, repetitive sequences and low-quality reads ( Table 1). The average GC content was about 53%, Q20 value was above 96%, Q30 value was above 89%, and the Q30 value of 95% of all samples was above 91%. Clean reads were then mapped against the reference genome using Hisat2. Mapping rate, number of clean reads, clean bases, and additional statistics are shown in Table 1. In addition, we also analyzed the FPKM distribution of each sample, and the results showed that the distribution of FPKM was uniform in each sample, which was suitable for subsequent analysis (Additional le 6).
We next employed correlation analysis to determine whether replicates had similar expression values, which is a well-established method for assessing data quality. Our analysis revealed that the majority of replicates had correlation coe cients of 0.97 or higher (Additional le 1). Additionally, the correlation coe cient between B4332T1h1 and B2T1h1 was only 0.83, indicating a large difference between the two materials.
Identi cation of transcription factors and differential gene expression Transcription factors are a group of protein that can bind speci cally to a speci c sequence at the upstream of a gene, so as to ensure the expression of the target gene at a speci c time and space with a speci c intensity. Transcription factor prediction was analyzed by iTAK software and HMMSCAN was used to identify TF, Transcription factor was compared to plant speci c transcription factor database transcriptome were identi ed ( Table 2). TFs form complex regulatory networks that control the expression of their target genes. These regulatory circuits are pivotal for coordinating transcriptional and posttranscriptional control of target genes [17]. SNPs in TFs have previously been shown to impact agriculturally important traits. For example, a negative regulator of DELLA genes has been shown to cause dwar sm in oil plants [18]. In Prunus, TFs control many agriculturally important traits such as the owering, fruit quality, and biotic and abiotic stress resistance [19].
EdgeR software was then used to analyze the differential expression of genes in each sample, and the corrected p-value and PADj value of the differential expression tests were calculated. The input data for differential expression analysis consisted of read counts obtained via quantitative analysis. Detailed label names and their treatment condition are shown in Additional le 8. The number of differentially expressed genes in each treatment is shown in Table 3.
As shown in Figure 1, there were 150 differentially expressed genes (DEGs) in A2 versus A3, among which 96 genes were up-regulated and 54 genes were down-regulated. There were 749 DEGs in A1 versus A2, among which 137 genes were up-regulated and 612 genes were down-regulated. There were 1753 DEGs in A1 versus A3, among which 678 genes were up-regulated and 1075 genes were down-regulated. There were a total of 14709 DEGs in B1 versus A1, among which 6147 genes were up-regulated and 8562 genes were down-regulated. There were 1649 DEGs in B2 versus B3, among which 658 genes were up-regulated and 991 genes were down-regulated. There were 17286 DEGs in B1 versus B2, among which 7275 genes were up-regulated and 10008 genes were down-regulated. There were 15194 DEGs in B1 versus B3, among which 6357 genes were up-regulated and 8837 genes were down-regulated. There were a total of 683 DEGs in B2 versus A2, among which 237 genes were up-regulated and 446 genes were downregulated. There were a total of 2704 DEGs in B3 versus A3, among which 2065 genes were up-regulated and 639 genes were down-regulated.
In order to further screen candidate herbicide genes, we believe that the following 6 contrast groups are more suitable for screening. B1 versus A1 was used to compare the difference between the sensitive and resistant material, while B3 versus A3 was used to compare the different results of cyhalofop-butyl treatment in the two L. chinensis populations. B2 versus B3 was used to measure the gene expression changes in sensitive L. chinensis genes caused by cyhalofop-butyl treatment. A2 versus A3 was used to determine which genes were induced by cyhalofop-butyl treatment in resistant L. chinensis genes. The four sets of DEGs were plotted on a Venn diagram ( Figure 2). A total of 12829 DEGs were found exclusively in the A1 versus B1 group, indicating that the two materials were signi cantly different even without cyhalofop-butyl treatment. Additionally, 1426 DEGs were found uniquely in the B3 versus A3 comparison, meaning these genes were only different during treatment with cyhalofop-butyl. These genes therefore represent possible candidates for the NTSR present in the resistant L. Chinensis line. Another 506 genes were found to be differentially expressed only in the B2 versus B3 comparison group, indicating that they represent a response present only in L. chinensis that is sensitive to cyhalofop-butyl and only 28 DEGs existed exclusively in the A2 versus A3 comparison, indicating that they represent a response present only in L. chinensis that is risistant to cyhalofop-butyl In addition, the B1 versus B2 comparison was used to determine the effect of water treatment on the expression of rice genes in the sensitive population, while the A1 versus A2 comparison was used to detect the in uence of water treatment on resistant rice. In the following analysis, we will compare the P450 family genes and ABC transporter family genes selected in the two groups with those selected in the above four groups, so as to exclude the interference of water solvent.
Functional annotation of DEGs and selected candidate genes We primarily focused on genes which were differently expressed between treated sensitive and treated resistant plants, which may represent targets for understanding the NTSR mechanism. CYP450 genes have previously been implicated in herbicide metabolism [20], while CYP76 family genes are involved in the metabolism of benzoyl urea herbicides and have been shown to play a detoxi cation role for benzoyl urea herbicides [8]. Researchers have found that CYP72A31 and CYP81A6 are involved in the detoxi cation of bensulfuron-methyl herbicides in rice and A. thaliana [9][10][11]. However, no studies on the role of P450 family genes in the regulation of cyhalofop-butyl resistance have been conducted in L. chinensis using whole transcriptomic methods. In higher plants, ABC transporters have been implicated in detoxi cation of xenobiotics, including herbicides. Unlike P450, which detoxi es herbicides through metabolic mechanisms, ABC transporters detoxify herbicides and confer herbicide resistance by transporting herbicides and the metabolites that result from their breakdown. The four comparison groups (B1vsA1, B3vsA3, A2vsA3, B2vsB3) which were most likely to be relevant to cyhalofop-butyl resistance were selected for subsequent correlation analyses ( Figure 2).
Expression analysis of sensitive versus resistant genotypes without exposure: A1vsB1 A total of 14709 genes were identi ed in this group, among which 6147 genes were up-regulated and 8562 were down-regulated ( Figure 3A). The up-regulated genes were primarily associated with biological processes such as single organism metabolic process, single organism process, and transmembrane transport. These genes were also enriched in molecular functions, such as transmembrane transport activity, transporter activity, and substrate-speci c transporter activity and enriched in Cell components such as plasmid and chloroplast ( Figure 3B). KEGG pathway analysis showed that some genes associated with metabolic pathways and secondary metabolite biosynthesis pathways were highly expressed in resistant L. chinensis ( Figure 3C). Analysis of all induced genes showed that there were 25 ABC transporter family genes with higher expression in the resistant population. Additionally, 39 CYP450 family genes were identi ed to be more highly expressed in the resistant population. Among them, 32 genes were differentially expressed only in this subset. Two CYP72A family genes (Chr4.g12261 and Chr4.g12260) and two CYP76 family genes (Chr17.g44850 and Chr17.g44847) were identi ed in this subset, all of which were highly expressed in LC1802. Members of the CYP76 family and CYP72A31 have been identi ed as herbicide metabolic resistance genes [8][9]. Therefore, we hypothesized that these genes might be impacting the NTSR of LC1802. Additionally, 48 genes related to fatty acid biosynthesis were identi ed, of which 20 were more highly expressed in the resistant population and 28 were more highly expressed in the sensitive population. Many genes associated with plant hormone signal transduction and photosynthesis were found to be more highly expressed in the resistant population. The DELLA protein has been shown to play a role in the regulation of gibberellin (GA) signaling by suppressing stress response pathways [21]. This gene is also associated with growth inhibition of weeds [22]. Interestingly, three genes encoding DELLA proteins were identi ed only in this subset of DEGs. Among them, Chr3.g10116 was found to have no expression in the resistant population, despite being present in the sensitive population.
Expression analysis of sensitive vs. resistant genotypes with herbicide exposure: A3vsB3 In this group, 2704 DEGs were identi ed, among which 2065 genes were up-regulated and 639 genes were down-regulated ( Figure 4A). A total of 1426 DEGs only existed in this subset. These genes showed no difference in expression when the two materials were not treated but showed differences in expression during cyhalofop-butyl treatment. Additionally, 1091 genes were differentially expressed in this set as well as the B1 versus A1 comparison, indicating that these genes were differentially expressed under both normal conditions and cyhalofop-butyl treatment. GO term analysis of biological processes resulted in the identi cation of up-regulated genes related to carbohydrate metabolic process, cell wall organization or biogenesis, and reactive oxygen species metabolism. These genes were also enriched in molecular functions, such as hydrolase activity, hydrolyzing O−glycosyl compounds, hydrolase activity, acting on glycosyl bonds and enriched in Cell components such as cell periphery and extracellular region ( Figure  4B). KEGG pathway analysis results were similar to the results seen in the B1 versus A1 comparison, and a total of 180 up-regulated genes were identi ed as being involved in different metabolic pathways. 130 were involved in secondary metabolite biosynthesis pathways ( Figure 4C). These ndings further support our hypothesis that metabolic genes and transport genes mediate the high cyhalofop-butyl resistance in the LC18002 line. Additionally, 11 ABC transporter family genes were identi ed among the DEGs, 9 of which were differentially expressed only in this subset. Among the up-regulated genes, we found that ABCA7 (Chr2.g06140) and ABCB2 (Chr10.g32569) were highly induced and may represent candidates for further study of cyhalofop-butyl resistance in LC18002. CYP450 family genes have been shown to play a crucial role in herbicide metabolism [20], and we found 16 CYP450 family genes which had higher expression levels in LC18002. In particular, we found that a CYP450 gene (Chr8.g25254) and CYP84A1 (Chr10.g30537) were up-regulated not only in this comparison, but also in B1 versus A1, despite showing no differential expression in any other comparisons. Given that other members of the CYP450 family have been implicated in herbicide tolerance, higher expression of these genes may play a role in the resistance of LC18002 to herbicides. Additionally, we found 12 possible xyloglucan endotransglucosylase/hydrolase encoding genes, and other genes of this class have previously been reported to participate in abiotic stress in pepper [23] and tea [24]. In addition to the above genes, some other key genes involved in herbicide resistance pathways were also found in this subset. Pectate lyases have been shown to play an important role in herbicide resistance by regulating the composition of polysaccharides in plants to affect plant stress tolerance [25]. Two pectate lyase genes were identi ed in this subset, both of which were highly expressed in LC18002. A protein strubbelig-receptor family 7 gene was also identi ed, and its homolog has been shown to be a sensor for herbicides, mediating intercellular and intracellular abscisic acid (ABA) pathways to regulate herbicide stress signaling [26]. Two aquaporin nip5-1 genes were also found to be more highly expressed in LC18002, and members of this class are known to control pore size with signi cant impacts on herbicide tolerance. Peroxidase 12-like genes have also been shown to play a part in abiotic stress tolerance and herbicide resistance by regulating the jasmonic acid biosynthesis pathway [27][28]. Pectinesterase inhibitor 34 is involved in ne-tuning cell wall remodeling processes under abiotic stress [29] and three up-regulated pectinesterase inhibitor genes were identi ed in this group, two of which were differentially expressed only in this subset.
Effects of herbicide exposure on the sensitive genotype: B3vsB2 A total of 1649 DEGs were identi ed in this group, among which 658 genes were up-regulated and 991 genes were down-regulated after herbicide treatment (Additional le 2A). This group can be used to investigate which genes were differentially expressed following cyhalofop-butyl treatment in the sensitive population. GO term analysis of biological processes indicated that many up-regulated genes were related to protein phosphorylation, phosphorylation, phosphorus metabolic process, and phosphate containing compound metabolic process. These genes were also enriched in molecular functions, such as small molecule binding, nucleotide binding and nucleoside phosphate binding and enriched in Cell components such as membrane, intrinsic component of membrane and integral component of membrane. (Additional le 2B). KEGG pathway analysis showed that thirteen genes with higher expression after cyhalofop-butyl treatment were involved in plant hormone signal transduction (Additional le 2C), including indole-3-acetic acid amido synthetase, TF PIF5, and regulatory protein NPR5. Seven genes were involved in glycerophospholipid metabolism process, while ve genes were involved in fatty acid degradation such as alcohol dehydrogenase and amino acid permease. Five other genes were involved in valine, leucine, and isoleucine degradation pathways. Eight ABC transporter family genes with differential expression were identi ed in this group. Among them, seven were induced and one was repressed. The induction of these different genes in the sensitive population likely represents a response to stress caused by cyhalofop-butyl treatment, including ABC transporters being used to transport active herbicide molecules.
Effects of herbicide exposure on the resistant genotype: A3vsA2 A total of 150 DEGs were identi ed in this group, among which 96 genes were up-regulated and 54 genes were down-regulated (Additional le 3A). This group can be used to investigate gene expression changes caused by cyhalofop-butyl treatment in the resistant line LC18002. As shown in Additional le 3B, among the up-regulated genes, twelve genes were annotated as being involved in stress response, nine genes were enriched in defense response, and nine other genes were enriched in response to acid chemical process. These genes were also enriched in cellar components, such as membrane part and intrinsic component of membrane. KEGG analysis showed that four up-regulated genes were annotated as being involved in metabolic pathways and three genes were involved in the biosynthesis of secondary metabolites (Additional le 3C). Out of the 150 genes which were differentially expressed during cyhalofop-butyl treatment, 28 showed expression change after herbicide exposures exclusively in LC18002. Additionally, three CYP450 genes related to xenobiotic metabolism were highly expressed in the herbicide-exposed resistant line. The expression levels of these genes were relatively constant between the two materials prior to treatment but their expression was signi cantly induced in the resistant line after treatment with cyhalofop-butyl. It therefore is likely that these genes play some roles in the herbicide resistance of LC18002, although further research is necessary to con rm this.
Effects of water exposure on the sensitive genotype: B2vsB1 This group can be used to identify the effects of water treatment on gene expression in sensitive populations. A total of 17283 DEGs were identi ed in this comparison, among which 7275 genes were upregulated and 10008 genes were down-regulated 24 hours after treatment with water (Additional le 4A). GO biological process enrichment analysis showed that among the up-regulated genes, 2004 genes were enriched in single organism process and 1198 genes were enriched in single organism metabolic process.
Molecular function analysis showed that 355 up-regulated genes were enriched in transmembrane transporter activity, 278 genes were enriched in substrate-speci c transmembrane transporter activity, and 145 genes were enriched in phosphoric ester hydrolase activity (Additional le 4B). Since GO enrichment analysis showed that many up-regulated genes had transmembrane transport activity, ABC transporter family genes were screened in this subset. A total of 29 ABC transporter family genes were found to be up-regulated 24 hours after water treatment, 27 of which were only present in this comparison. Although both water treatment and cyhalofop-butyl treatment affected the expression of ABC transporters, the actual genes affected differed. This result implies that the water present in the cyhalofop-butyl solution is not the sole reason behind the induction of ABC transporters seen under its application.
Effects of water exposure on the resistant genotype: A2vsA1 This comparison group can be used to identify the effects of water treatment on gene expression in LC18002. A total of 749 DEGs were identi ed in this comparison, among which 137 genes were upregulated and 612 genes were down-regulated 24 hours after treatment with water (Additional le 5A). Among the up-regulated genes, 7 genes were annotated as participating in the protein folding process, 3 genes were associated with response to light intensity, and 3 other genes were implicated in response to heat. Cellular component analysis indicated that 51 genes were associated with intracellular part, 45 genes were associated with membrane bounded organelle, 45 genes were associated with intracellular membrane bounded or ganelle and 23 genes were associated with the plastid (Additional le 5B). Unlike the A2 versus A3 comparison, only 3 CYP450 family genes were identi ed in up-regulated genes of this subset which had no differential expression in group A2 versus A3, indicating that the type of CYP family genes involved in water application and cyhalofop-butyl metabolism was different. In addition, an ABC transporter family gene was found to be up-regulated 24 hours after water treatment but was not found to be differentially expressed in the A2 versus A3 comparison. These results again indicate that water treatment alone results in a signi cantly different response compared to cyhalofop-butyl treatment.

Phylogenetic analysis of candidate genes
Phylogenetic analysis was performed to determine whether the candidate genes identi ed had close homologs in other species or belonged to any known gene families. The protein sequences of all differentially expressed P450 family genes were compared against P450 genes from A. thaliana, rice, soybean, sorghum, and maize. This clustering analysis revealed that L. chinensis P450 genes clustered closely with P450 genes of other species ( Figure 5). Chr8.g25254 was found to cluster closely with Cytochrome P450 71T4-like genes. Chr16.g43989 was annotated as a member of the CYP734A6 subfamily, but phylogenetic analysis shows that it is clustered closely with the CYP72A subfamily of rice. Additionally, we also included herbicide metabolism genes previously studied, such as OsCYP81A6 and OsCYP72A31P. This analysis revealed ve genes that were closely related to OsCYP81A6, and seven genes that were closely related to OsCYP72A31P. In addition, we selected one DELLA protein and two pectate lyase proteins as the outgroup and found that they produced two independent branches.
qRT-PCR validation of RNA-seq expression patterns In order to con rm the magnitude of differential gene expression obtained by transcriptomic approach, the expression levels of 8 genes were measured in two L. chinensis accessions using quantitative realtime PCR (qRT-PCR) ( Figure 6). The primers utilized to amplify candidate genes were designed based on the sequence of the identi ed genes and are listed in Additional le 8. There is a high level of correlation was found between the RNA-seq and qRT-PCR data.

Discussion
In this study, we show that resistance to cyhalofop-butyl in LC18002 is likely due to NTSR mechanisms.
Metabolism-based resistance to ACCase inhibitor herbicides has long been reported and is increasing in weed species [14]. However, there are few studies on herbicide metabolism and resistance genes at the whole transcriptome level, which has slowed the development of new control mechanisms. RNA-seq has been successfully used to identify genes involved in metabolic resistance to ALS herbicides in two grass weed species: L. rigidum [15] and A. myosuroides [16]. However, there has been no study focused on the NTSR mechanisms of L. chinensis at the whole transcriptome level. Identifying genes involved in NTSR is important for understanding the evolution of metabolic resistance and improving weed management strategies in the eld.
In this study, 266750156 clean reads were generated from L. chinensis by Illumina Hiseq 2500 technology. After identifying DEGs, we primarily focused on genes which were differentially expressed between sensitive and resistant plants. In this group, 2704 DEGs were identi ed, among which 2065 genes were up-regulated and 639 genes were down-regulated. A total of 1426 of these DEGs only existed in this subset. Additionally, 1091 genes were differentially expressed in this set and B1 versus A1, implying that these genes had different expression patterns both with and without cyhalofop-butyl treatment. In this set, we identi ed two P450 family genes, a CYP450 gene (Chr8.g25254) and CYP84A1 (Chr10.g30537), which may be involved in the regulation of the resistance phenotype of LC18002. In higher plants, ABC transporters have been implicated as playing a role in detoxi cation of herbicides and other xenobiotics. Unlike P450, which detoxi es herbicides through metabolic pathways, ABC transporters detoxify herbicides and confer herbicide resistance by mobilizing the herbicide and its metabolites. In our study, 11 ABC transporter family genes were identi ed among the DEGs in B3 versus A3. Among the upregulated ABC family genes, we found that the ABAC7 (Chr2.g06140) and ABC transporter b family member 2-like (Chr10.g32569) were not differentially expressed during water treatment but did change under cyhalofop-butyl treatment. This gene may therefore play a role in the cyhalofop-butyl tolerance of L. chinensis.
The CYP450 family of genes are known to play a signi cant role in herbicide metabolism [20], while CYP76 family genes have been shown to detoxify benzoyl urea herbicides [8]. Researchers have found CYP72A31 and CYP81A6 are involved in the detoxi cation of bensulfuron-methyl herbicides in rice and Arabidopsis [9][10][11]. In our study, many DEGs were found to be highly homologous to these two gene families, indicating that they may also play a role in herbicide tolerance.

Conclusions
This is the rst report utilizing whole transcriptome analysis to identify putative herbicide tolerance genes in L. chinensis. An L. chinensis population highly resistant to the ACCase inhibitor herbicide cyhalofopbutyl, which also lacks any mutations in the target proteins for this herbicide, was selected for this analysis. The NTSR observed for this line is likely due to P450-mediated metabolic resistance and ABC transporter mediated sequestration of the metabolites. CYP450 (Chr8.g25254), CYP84A1 (Chr10.g30537), two CYP72A family genes (Chr4.g12261 and Chr4.g12260), two CYP76 family genes (Chr17.g44850 and Chr17.g44847), ABCA7 (Chr2.g06140), and ABCB2 (Chr10.g32569) were all differentially expressed in resistant compared to susceptible L. chinensis lines, making them possible markers or new target genes. In addition, we identi ed several P450 family members in the DEGs, some of which were revealed by phylogenetic analysis to be similar to known herbicide tolerance genes. We found several other DEGs belonging to families associated with herbicide tolerance, including ABCs, DELLA genes, xyloglucan endotransglucosylase/hydrolases, glutamate dehydrogenases, methyl crotonoyl carboxylases, aquaporin genes, and thaumatins, all of which represent targets for additional investigations. Overall, this study represents a signi cant step forward in understanding NTSR in L. chinensis and provides a number of new avenues to explore in future studies.

Plant Material
The seeds of two contrasting L. chinensis genotypes, herbicide tolerant (LC18002) and herbicide susceptible (LC17041) were polished by abrasive paper for 30 s, Pour the air-dried sand through a 20mesh sieve into a plastic cup with a diameter of 9.5cm and add water to completely moisten the soil. After soaking the seeds in water for 8~24 h, rinse them and air dry them naturally. The seeds are then spread evenly over the moist soil and a thin layer of soil is placed on top. After sowing, the seeds were cultured in an arti cial climate chamber. Two culture conditions were set. reads were generated. Readscontaining more than 10% poly-N and more than 50% low-quality reads (Q≤20) were removed from the raw data using Trimmomatic v0.33 [33]. Concurrently, the Q20 and Q30 values, GC content, and sequence duplication level of the clean data were calculated. All downstream analyses were based on clean, high quality data.

Identi cation of Differential Expressed Genes (DEGs) and Transcription Factors
In order to quantify the gene expression, count-based method, FeatureCounts was used. Finally, the transcript counts were used for pairwise differential gene expression analysis using edgeR package [34].
A cut-offvalue of log2 ratio ±1 and q-value 0.05 were used to lter out the signi cant transcripts in each case. Further, transcription factors were identi ed in L. chinensis, Transcription factor prediction was analyzed by iTAK software and HMMSCAN was used to identify TF, Transcription factor was compared to plant speci c transcription factor database PlnTFDB.

Functional annotation and classi cation
Gene Ontology (GO, http://www.geneontology.org/) enrichment analysis of the DEGs was implemented using GOseq (v. 1.22) [30] using Wallenius'noncentral hypergeometric distribution, which can adjust for gene length bias in DEGs. Affected pathways were determined using Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.kegg.jp) [31]. We used KOBAS [32] software to annotate and identify the enriched KEGG pathways of the DEGs

Selection of candidate non-target genes
Genes that were commonly expressed between treated S and T, and non-treated S and T, and between treated T and treated S plants were selected for further evaluation based on their gene ontologies (GO).Genes that were assigned with GO molecular function and biological process related to metabolism and signaling pathways (oxidoreductase activity, nuclear acid binding transcription factor activity, hydrolase activity, transferase activity, transmembrane transporter activity, transferase activity, protein transporter activity, biosynthetic process, small molecule metabolic process, signal transduction, homeostatic process, immune system process, cell wall organization, secondary metabolic process, nitrogen cycle metabolic process) were evaluated based on UniProt and their foldchange. Contigs with predicted annotations related to stress response, signaling, transcription factors, and herbicide metabolism were selected as potential candidate NTS genes involved in cyhalofop-butyl tolerance.

Phylogenetic analysis
To build a neighbor-joining tree, we selected amino acid sequences of candidate genes and obtained amino acid sequences of homologous genes in A. thaliana, rice, sorghum, maize, and soybean from NCBI (https://www.ncbi.nlm.nih.gov/) and CYP databases (https://drnelson.uthsc.edu/CytochromeP450.html), as well as amino acid sequences of identi ed genes involved in herbicide metabolic pathways. We constructed a phylogenetic tree using MEGAX with 1000 bootstrap replicates. Total number of TF predicted 1,353 Table 3 Summary of the number of differential expression genes in each group