Seeds of 19 populations of L. multiflorum were collected. Among them 18 populations were sampled from different locations of Denmark, whereas one population (ID-27) was obtained from France. The majority of populations were collected from fields where unsatisfactory control of ALS and ACCase inhibitor herbicides had been observed. The L. multiflorum is not an endangered species and no specific permit is required for collection of seeds from agricultural fields. Two known susceptible populations ID-290, a population maintained by the Department of Agroecology, Aarhus University and the commercial forage grass variety Sikem were used as standard reference populations.

Two herbicides, iodosulfuron-methyl-sodium and clodinafop-propargyl, belonging to the groups of ALS and ACCase inhibitors, respectively, were applied through foliar application to L. multiflorum plants at the 3–4 leaves stage. The common trade names of the herbicides are Hussar OD (100 g L⁻¹ iodosulfuron-methyl-sodium + 300 g L⁻¹ mefenpyr, Bayer Crop Science, Denmark) as an ALS inhibitor and Topik (100 g L⁻¹ clodinafop-propargyl + 25 g L⁻¹ cloquintocet-mexyl, Syngenta, Denmark) as an ACCase inhibitor. All spray solutions were prepared in deionized water and applied in mixture with 0.5 L ha⁻¹ Renol (vegetable oil) using a laboratory pot sprayer fitted with two Hardi ISO F-110-02 nozzles. The nozzles were operated at a pressure of 310 kPa and a speed of 5.5 km ha⁻¹ delivering a spray volume of 150 L ha⁻¹. Herbicides solutions were prepared at 2 g a.i. ha⁻¹ iodosulfuron-methyl sodium (≈0.0013%) and 25 g a.i. ha⁻¹ clodinafop-propargyl (≈0.0167%). Plants were harvested 4 weeks after herbicide application and fresh weights were recorded.

In order to get idea about resistance level of collected populations, initial screening was carried out. The populations were chosen keeping in view the availability of seeds for further experiments. Fifteen populations were tested with ALS and ACCase inhibitors at a rate known to provide full control of susceptible L. multiflorum plants grown in pots (2 g a.i. ha⁻¹ iodosulfuron-methyl sodium (≈0.0013%), 25 g a.i. ha⁻¹ clodinafop-propargyl (≈0.0167%). For these purpose three trays of each population were prepared, one serving as control while the other two trays were treated with iodosulfuron-methyl sodium and clodinafop-propargyl. Plants were grown in a potting mixture consisting of soil, peat, and sand (2:1:1 w/w) containing all necessary macro- and micronutrients. After seedling emergence in a glass-house, trays were placed outdoor and thinned to one plant per hole in the tray. In total 104 individual plants per tray were maintained except for population ID-290 with only 45 individual plants due to a low germination. The herbicide efficacies were expressed as percentage survival and percentage decrease in fresh weight compared to control 4 weeks after herbicide application (Table 1).

Keeping in view initial screening and previous bioassays conducted in our laboratory, 17 populations were included in the dose-response studies. Due to low germination of the reference population ID-290, the commercial variety Sikem was also included as reference control. Seeds were sown in 2 L pots filled with the potting mixture used for the herbicide screening trials. Pots were placed in a glasshouse and sub-irrigated automatically one to three times per day depending on plant size. After germination the number of seedlings per pot was maintained at six and experiment was laid out in a completely randomized block design with three replications. Statistical analysis were done using R statistical software (R Development Team, http://www.r-project.org) with add on package drc (Ritz and Streibig, 2005; Knezevic et al., 2007) and a three-parameter log-logistic model were fitted to the data applying the drc modelFit function: Y=D1+exp[b(log(x)-log(ED50))] where Y is the fresh weight, x represents herbicide dosage (g a.i. ha⁻¹), D is mean response when herbicide dose is close to zero, ED₅₀ is the dosage (g a.i. ha⁻¹) that reduces fresh weight by 50% and b is the slope of the curve around ED₅₀. The ED₅₀ values for each population and herbicide are shown in Table 2 along with a resistance index (RI), which is calculated as the ratio between the estimated ED₅₀ value of the tested and susceptible populations (ED_50R/ED_50S).

Leaf materials from all population were harvested for DNA extraction. Briefly, genomic DNA was extracted from plants using a DNeasy Plant Kit (Qiagen, USA). PCR was performed using the universal primers of ACCase genes as described by Délye et al. (2011) and specific primers of ALS gene of L. multiflorum from the GenBank accession number AF310684.1 (Table S1). PCR products of respective genes were amplified using initial denaturing for 5 min at 94°C, followed by 35 cycles consisting of 94°C for 40 s, 62°C for 35 s, 72°C for 30 s, and 72°C for 5 min for final extension. The PCR products were purified, sequenced and analyzed for single nucleotide polymorphism (SNP) using CLC main workbench according to manufacturer instructions (CLC bio, Aarhus, Denmark).

Total RNA was extracted from 50 mg of leaf material of individual plants of six selected populations using RNeasy Plant Mini Kit (Qiagen, Stanford, California, USA). RNA was extracted from three independent biological samples of selected populations. To eliminate genomic DNA from RNA samples they were treated twice with RNase-Free DNase (Qiagen, USA). Concentration and purity of total RNA were measured through NanoDrop 1000 spectrophotometer (Thermo Scientific, USA). Samples with concentration more than 100 ng μL⁻¹ and optical density absorption ratios of A₂₆₀/A₂₈₀ and A₂₆₀/A₂₃₀ between 1.8–2.2 and 2.0–2.2, respectively, were taken for cDNA synthesis. For cDNA synthesis and differential gene expression we followed the guidelines proposed by Bustin et al. (2010). The qPCR reactions were performed using the StepOnePlus™ Real-Time PCR System Thermal Cycling Block (Applied Biosystems Foster City, USA). Two internal control genes were chosen based on their stable expression to herbicide treatment as found by Gaines et al. (2014). Same gene specific primers as described by Gaines et al. (2014) were used and presented in Table S2. Moreover, amplicon sizes were checked on 2.5% (w/v) agrose gels. Reactions were done in duplicate and a negative control consisting of template without primers were also included for each primer. Reactions were also conducted as described by Gaines et al. (2014) with some modification. Briefly, 20-μL volume of reaction included 10 μL of SyberGreen Master Mix, 2 μL of 1:10 diluted cDNA, 5 μL of 0.5 p mol μL⁻¹ primers (1:1 mix of forward and reverse primers), and 3 μL of nuclease-free distilled water. Reaction conditions included 15 min incubation at 95°C, then 40 cycles of 95°C for 30 s and 60°C for 1 min, followed by a melt-curve analysis to confirm single PCR product amplification. Threshold-cycle (ΔC_T) values were calculated for each reaction. Gene-specific PCR efficiency was used to calculate the expression of target genes relative to the expression of internal reference genes. Equivalent slopes for target and internal control genes were observed in amplification plots. The ΔC_T value was calculated as follows: ΔC_T (target genes) = C_T (target gene) − C_T (reference gene), where C_T is the cycle number at which PCR product exceeded a set threshold. Relative transcript level (RTL) was calculated through = 1 × 2^−ΔCT. Transcript levels were calculated from three independent biological samples using three technical replicates of each biological sample.

The orthologs of HMR genes were identified in A. thaliana and O. sativa. The names of these orthologous are glycosyl-transferase (GT, At2g15490), nitronate monooxygenase (NMO, At5g64250), cytochrome P450 (CYP72A-1, At3g14680), and cytochrome P450 (CYP72A-2, At3g14690) with accession numbers NM_127109, NM_180932, NM_112329 and NM_112330 respectively, in A. thalina. Orthologs of HMR genes in the O. sativa were identified as Os01g0638600 (GT), Os08g0485400 (NMO), Os01g0627500 (CYP72A-1), and Os01g0728300 (CYP72A-2) with accession numbers NM_001050207, NM_001068622, NM_001050167, and NM_001050664, respectively. Rab GTPase (RGTP) and isocitrate dehydrogenase (IDE) were used as internal reference genes and identified as Os01g0276100 (NM_001049260) and Os01g0558600 (NM_001049871) in O. sativa, and At5G03290 (NM_120407) and At5G47200 (NM_124091) in A. thaliana. All the identified accession numbers of HMR orthologs were analyzed manually and validated through BLAST search using NCBI database (http://blast.ncbi.nlm.nih.gov/Blast.cgi).

The behaviors of HMR genes under various experimental conditions and at different developmental stages were explored using the Genevestigator (http://www.genevestigator.com). The Genevestigator is a manually curated and well-annotated database of expression profiling from 11 different plant species with more than 26,000 exclusive plant samples (www.genevestigator.com, Feb. 2016). In silico analysis of the expression of identified transcription factors (TFs) was also explored using the condition tool of this database.

The 1 kb upstream sequences from translational initiation codon of HMR orthologous in A. thaliana and O. sativa were obtained. In order to check occurrence of transcription factor (TF) in the promoter region of HMR genes, the promoters of HMR orthologous were analyzed using the web-based tool “The Plant Promoter Analysis Navigator (PlantPAN 2.0; http://PlantPAN2.itps.ncku.edu.tw).” This tool provides resources for detecting corresponding transcription factors and other regulatory elements in a promoter of a gene (Chow et al., 2015). The analysis was carried out separately on the promoters of each HMR orthologous of A. thaliana and O. sativa using Promoter analysis function of PlantPAN 2.0. This function recognizes the combinatorial cis-regulatory elements and associates them with corresponding transcription factors (TFs). Venn diagrams were plotted using online tools Venny 2.0 http://bioinfogp.cnb.csic.es/tools/venny/. All identified TFs were analyzed for biological functions using the AgriGO through singular enrichment analysis tool (http://bioinfo.cau.edu.cn/agriGO/). This analysis was performed with FDR correction and Fisher's exact test <0.05. Furthermore, identified TFs in A. thaliana and O. sativa were analyzed to find common and conserve TFs in both species. The InterPRO analysis was carried out to find conserve domain in TFs exclusively linked to the promoters of HMR genes.

The overrepresented cis-motif consensus patterns were identified using the Multiple Expectation maximization for Motif Elicitation (MEME) analysis tool (Bailey et al., 2006). The MEME was used to search best 5 cis-motif based on E-value, consensus patterns of 6–50 bases width, only on the forward strand of the input sequences and the distribution model used was the default Zero Or One Per Sequence (ZOOPS). A motif is a sequence pattern that occurs repeatedly in a group of DNA sequences. The motifs identified using MEME were analyzed through GOMO (Gene Ontology for Motifs; Buske et al., 2010). The purpose of GOMO is to identify possible roles (Gene Ontology terms) for DNA binding motifs. GOMO takes a motif and determines which GO terms are associated with the (putative) target genes of the binding motif. GOMO was run using the “multiple species category” which gave access to the plant database with significant threshold q < 0.01 and the number of score shuffling rounds over 1000.

In our primary screening assay, all tested populations showed higher tolerance to either one or both of the herbicides compared to our reference susceptible population (ID-290). However, we found different levels of resistance among the tested populations. For example, population ID-1023 was the most resistant population to iodosulfuron-methyl-sodium with a survival percentage of 96.2% and decrease in fresh weight of only 9.9% (Table 1). Population ID-903 also had a high survival rate (96.2%) and exhibited a low decrease in fresh weight (20.5%) following iodosulfuron-methyl-sodium exposure. The most resistant populations to clodinofop-propargyl were ID-903 and ID-06 with 88.5 and 94.3% survival, and 24.4 and 32.7% decrease in their fresh weights, respectively, while ID-1023 was found to be moderately resistant with 75% survival and 54.4% decrease in fresh weight (Table 1). Other populations had a very interesting resistance pattern to both herbicides such as ID-04, which had 97.1% survival and 29.5% decrease in fresh weight after exposure to iodosulfuron-methyl-sodium, whereas clodinofop-propargyl showed 60.6% survival and 67% decrease in fresh weight. In order to obtain a deeper insight into the resistance pattern of the populations, a dose response experiment was carried out.

The susceptible reference populations (ID-290 and Sikem) were found to be susceptible to both herbicides. Populations ID-1023 and ID-07 were found to be highly resistant to iodosulfuron-methyl-sodium with a RI of more than 47.4 (Table 2). The exact ED₅₀ of ID-1023 and ID-07 could not be estimated, as fresh weight was higher than 50% over the whole dose ranges. The populations ID-06, ID-904, ID-20, and ID-04 exhibited 2.3, 3.3, 3.7, and 6.4 fold higher tolerance, respectively compared with the susceptible reference (Table 2).

In contrast to iodosulfuron-methyl-sodium, ID-1023 was susceptible to clodinafop-propargyl and the resistance index of ID-07 was only 1.4-fold higher than the standard reference used in this study. The populations ID-974, ID-20, and ID-06 showed 1.6, 2.3, and 3.9 fold higher tolerance, respectively (Table 2). Resistance index of ID-23, ID-04, ID-691, and ID-27 to clodinafop-propargyl were not determined due to low seed germination. The population ID-27 was obtained from France and found resistant to ACCase herbicides in our earlier bioassays (data not shown). ID-06 and ID-20 were moderately resistant with RIs of 3.9 and 2.3 to clodinafop-propargyl (Table 2).

Comparison of the PCR fragments of the ALS gene from various L. multiflorum populations revealed a single nucleotide change from TGG to TTG at position 574 in 78 and 89% individuals of populations ID-1023 and ID-903, respectively (Table 3). This change leads to substitution of a tryptophan to leucine at 574 (Trp-574-Leu) and previous studies suggest that this mutation results in resistance to all classes of ALS inhibitors (Bernasconi et al., 1995; Bi et al., 2016).

PCR fragments of the ACCase gene also revealed a single nucleotide change from ATT to AAT leading to substitution of isoleucine to asparagine at 2041 (Ile-2041-Asn) in one population ID-903 (Table 4). This Ile-2041-Asn substitution was reported to cause resistance to ACCase inhibitors (Délye et al., 2003, 2005; Martins et al., 2014).

The expression patterns of four herbicide metabolism genes were tested in six populations of L. multiflorum without any herbicide exposure. These populations were selected based on their resistance patterns to ALS and ACCase inhibitors. Among the selected populations two (ID-23 and Sikem) were susceptible, two (ID-1023, ID-07) were resistant to ALS inhibitors, one (ID-27) was resistant to both ALS and ACCase inhibitors and one (ID-20) was found to be moderately resistant to both inhibitors. Before expression analysis, individual plants of these populations were tested for the presence of ALS and ACCase target site mutations. This revealed the presence of a mutation Trp-574-Leu in population ID-1023. In qPCR analysis, all four genes were exhibiting higher expression in some plants of the populations ID-20, ID-27, and ID-07. For example ID-20-3, ID-27-3, and ID-07-1 showed significant over-expression of all four genes (Table 5). ID-27-3 exhibited the highest expression of NMO, GT, and P450 among all populations. NMO was also highly expressed in all individual plants of population ID-27 and ID-07 (Table 5). However, there is considerable variation in the expression of GT, P450s, and NMO in populations and even in individuals of the same populations (Table 5). For example NMO was highly expressed in ID-20-3, but less expressed in the other plants of the same population, whereas expression of NMO was significant increased in all plants of ID-27, but P450s and GT were not significantly altered in ID-27-1. Nonetheless, some individuals of the three populations (ID-20, ID-27, and ID-07) exhibited significant constituent higher expression of HMR genes compared to the susceptible population (Table 5). Based on these results, we can conclude that populations ID-20, ID-27, and ID-07 showed signs of evolution of metabolic based herbicide resistance.

In order to check the behavior of HMR genes under various stresses, we utilized the publicly available gene expression database (Genevestigator). We checked the differential expression of HMR genes during different stages of plant development. NMO and CYP72A-1 had similar and high level of expression throughout developmental stages with maximum expression at senescence (Figure 1A). GT and CYP72A-2 had medium to low expression during different developmental stages except senescence, where CYP72A-2 tended toward higher expression (Figure 1A). Reference genes used in the study consistently showed a very high and stable expression pattern during different growth stages pointing them as good internal reference genes. In order to check behavior of these genes under various stress conditions, in silico expression was performed using the condition tool of the Genevestigator. This analysis revealed that 34 different perturbations (out of 3283 available in database) cause significant change in the expression (p < 0.001 and fold change >3.00) and the expression was most pronounced under various chemicals stress (Figure S1, Figure 1B). The strongest response was seen for the herbicide safener fenclorim, where expression increased 16-fold for NMO, 114-fold for GT, 9-fold for CYP72A-1, and 2.5-fold for CYP72A-2. Similarly herbicides such as dicamba and plant secondary metabolites phytoprostane strongly induce expression of these genes, whereas our control genes remain unaffected (Figure 1B). These results suggest that these herbicide metabolism genes have a role to detoxify herbicides and validate our observations in L. multiflorum.

Discovery of regulatory cis-elements in the promoter regions is essential to understand the spatial and temporal expression pattern of involved genes in a specific function. These genes may be regulated by a common set of transcription factors, and can be detected by the occurrence of specific cis-regulatory motifs in the promoter region. Hence we analyzed the promoters regions of HMR genes and reported top five motifs common in both A. thaliana and O. sativa. Among the top five common motifs present in the investigated promoters of both plants, only the first motif was significant based on E-value (Table 6). In GOMO analysis, this motif was found to be involved in making CUL4-RING ubiquitin ligase complex with sigma factor activity and localized in chloroplast and mitochondrial inner membrane. Cullin-RING (CUL) complexes represent a predominant group of ubiquitin E3 ligases with vital role in coping with environmental stresses and maintaining genome stability (Guo et al., 2013). Remaining top four common motifs were not considered significant based on E-value. However, these were all involved in transcription factor activity, helicase activity, regulation of transcription, kinase activity and circadian rhythm (Table 6).

In order to check the occurrence of TFs associated to the promoter regions of HMR genes, the PlantPAN promoter analysis tool was utilized. The analysis elucidate TFs associated to the promoters of individual gene and can determine common TFs among the four HMR genes. Interestingly, 136 TFs were common in all four HMR genes in A. thaliana and 85 TFs were found to be commonly associated to the promoters of HMR genes of O. sativa (Figure 2, Table S3). Among them, 14 TFs were conserved and exact homologs between A. thaliana and O. sativa (Table 7). In order to get deeper insights about these conserved TFs, their expression patterns were analyzed at different developmental stages and under various stress conditions. This revealed that two C2H2 type zinc finger transcription factors (ZAT6 and ZAT10) were up regulated in response to herbicides, whereas squamosa promoter-binding (SPB) type did not respond to chemicals (Figure S3). The expression of ZAT6 and ZAT10 is very similar to HMR genes during different growth developmental stages and under chemical stress (Figure S3). These zinc finger transcription factors had the highest expression level throughout developmental stages with top most expression at senescence and highly up regulated in response to herbicides and herbicide safeners (Figure S3).

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.