The 421 cultivated rice varieties used in this study were selected from diverse germplasms collected domestically and internationally, and the names, geographical origin, and subpopulation information of which have been previously reported (Zhao et al., 2015). In 2022, these rice varieties were planted at the South Breeding Base of Sichuan Agricultural University in Lingshui County, Hainan Province, China. The seeds of each were individually harvested in April 2023, uniformly air-dried and placed in the seed storage facility in the Rice Research Institute of Sichuan Agricultural University in Chengdu, Sichuan Province.

Sixty grains of each accession were sown in a flowerpot (15 cm diameter and 16 cm height) in July 2023 with three replications for herbicide evaluation in the Chengdu Campus of Sichuan Agricultural University. Half-lethal dose of glyphosate, glufosinate and mesotrione solutions was sprayed, respectively at the 3-4 leaf stage as 2.25 L/hm² of 41% glyphosate (Fuhua Chemical), 3.00 L/hm² of 20% glufosinate (BASF) and 1.20 L/hm² of 10% mesotrione (Hubei Jiahui Xingcheng Biotechnology Co., Ltd.). The response of each pot to herbicide injury was scored with 0-4 scale five days after the application of glyphosate and glufosinate, but seven days after the application of mesotrione, respectively, where 0 was for the plants with no apparent injury, level 1 for slight effect on seedlings, level 2 for moderate effect, level 3 for severe effect and 4 with the most severe injury.

GWAS analyses for herbicide resistance traits were performed for the entire panel and for indica and japonica group separately, using mixed linear models to take population structure and relative kinship into consideration for statistical association. SNPs with a minor allele frequency (MAF) greater than 5% and a missing rate less than 15% were selected for the association analysis, and 6.3 million variants, including 5.5 million SNPs and 800,000 InDels, were used for GWAS, as well as estimating population structure and kinship coefficients. The thresholds of genome-wide significance were determined using a modified Bonferroni correction as described by Li et al. (2012), where the total number of SNPs (M) used for threshold calculation was replaced by the effective number of SNPs (Me), in mixed linear models provided by EMMAX program (Kang et al., 2010). The physical locations of the SNPs and InDels were identified based on the Rice Annotation version 7.0 of the variety Nipponbare (Kawahara et al., 2013) from the website http://rice.uga.edu/.

The Bayesian clustering program fast Structure (Raj et al., 2014) was used to calculate different levels of K (K = 2-5), and the command choose_K.py was employed to identify the model complexity for maximizing the marginal likelihood. LD heatmaps and gene haplotypes were constructed using SNPs and InDels (MAF >0.03) from the 1 Mb regions around each peak SNP using R package ‘LD heatmap’ (Shin et al., 2006).

We annotated the functional variations within a 1 Mb region surrounding the associated loci using snpEff software (Cingolani et al., 2012), and classified all polymorphisms in the candidate region into four groups, as for 1) those predicted to induce amino acid exchanges or to change splicing junctions (GT or AG at the beginning or end of an intron, respectively), 2) those located at the 5′ flanking sequences of genes (e.g., within 2kb of the first ATG, typically the promoter region), 3) those located within a gene but did not meet the criteria for Group I or II (for example, located in a coding region but not predicted to change an amino acid, an intron, or a 3′ noncoding sequence), and 4) those located outside of coding regions. By identifying the most significant variant in Group I, potential candidate genes were determined.

The haplotype of each candidate gene was analyzed by downloading the genotypes of variant sites located 2kb upstream of the transcription start site and within the coding sequence of the gene from the RiceVarMap 2.0 (ricevarmap.ncpgr.cn) (Zhao et al., 2021). To improve the efficiency of the haplotype analysis, those variant sites with a minor allele frequency less than 0.05 and a missing rate greater than 10% were filtered out. The genotypes of remaining high-quality variants were used for haplotype classification. To better display the distribution pattern of haplotypes, only those haplotypes containing more than 10 varieties were shown. The information on geographical distribution of the 421 accessions was also downloaded from the RiceVarMap website. The world map was plotted using the ‘maps’ package in R (Becker and Wilks, 1995), and the characteristics of different haplotypes on geographical distribution were displayed in pie chart form using the ‘ape’ package in R (Paradis and Schliep, 2019).

Violin plots, correlations, and multiple comparisons were constructed using phenotypic means from three replicates for each accession. The P values for Pearson’s correlation coefficients were calculated by two-sided t-tests using the cor.test() function in R (Ihaka and Gentleman, 1996). Multiple comparisons of phenotypic data among different subgroups or haplotypes were conducted using the duncan.test() command in the ‘agricolae’ package in R, with a significant probability at 0.01. The ‘vioplot’ package in R was used to create violin plots for the multiple comparisons (https://github.com/TomKellyGenetics/vioplot).

By grading the herbicide resistance, we observed significant differences among the 421 accessions responding to the three herbicides, indicating the feasibility of the grading standards ( Figures 1A ; Supplementary Figure S1 ). The cultivars generally exhibited low resistance to all three herbicides, with average resistance ratings around level 3. Less than 5% cultivars exhibited relatively high resistance (resistance ratings smaller than level 2 to the three herbicides, and a Temperate Japonica cultivar HB-6-2 was the relatively highest glufosinate and glyphosate resistant cultivar in the tested materials ( Supplementary Figure S1 ; Supplementary Tables S4 , S5 ). There is a significant positive correlation between the resistance to the three herbicides in cultivated rice (r² = 0.26,0.13,0.14), although the correlation coefficients are relatively low ( Figure 1B ). Among them, the positive correlation between glufosinate resistance and glyphosate resistance is the most significant, with a correlation coefficient of 0.26 (p= 7.5 × 10^-5).

A total of 6.3 million high-quality variants, including 5.5 million SNPs and 0.8 million InDels, queried from RiceVarMap (http://ricearmap.ncpgr.cn), were used to estimate population structure and Kinship coefficients, as well as GWAS. When K equaled 2, the 421 accessions in the whole panel were divided into 300 Indica, 106 Japonica and 15 Intermediate ones. And When K equaled 5, the 300 Indica accessions were further divided into 38 Aus, 143 Indica I (IndI), 119 Indica II (Ind II), while the 106 Japonica accessions were divided into 60 Temperate Japonica (TeJ) and 46 Tropical Japonica (TrJ). Based on this population structure, we noticed a varying resistance to different herbicides of different subgroups ( Figure 1C ). For glufosinate and glyphosate, the TeJ subgroup exhibited significantly higher resistance than other subgroups ( Figure 1C ); for halosulfuron-methyl, the Aus subgroup showed significantly lower resistance than other subgroups. Thus, the cultivars with different genetic backgrounds demonstrated diverse herbicide resistance.

The strict inbreeding mating habits of rice have led to significant subpopulation structure and considerable linkage disequilibrium (LD), which reduce the mapping resolution of association studies and increase the occurrence of type I errors (Zhou et al., 2017). To address these issues, we used a mixed model to correct the interferences from subpopulation structure by incorporating the first three principal components (PCs) as covariates in the model, and conducted genome-wide association studies (GWAS) in the Indica group (303 accessions), Japonica group (181 accessions), and the entire group (all varieties) (see Methods). The calculated genome-wide significance thresholds, based on a nominal level of 0.05, were P = 6.6 × 10^-8, 8.7 × 10^-8, and 2.0 × 10^-7 for the whole panel, Indica and Japonica, respectively ( Supplementary Table S1 ).

Five significant loci associated with resistance to glufosinate and glyphosate, but none to mesotrione were identified in the whole panel ( Figure 2 ; Supplementary Figures S2 , S3 ). Detailed association results for each herbicide in each subpopulation are presented in Table 1 . Three loci associated with glufosinate resistance in the whole panel totally explained 25.76% of the phenotypic variance, with effect of each varying from 15.95% to 19.20% ( Supplementary Table S2 ). For glyphosate resistance, two loci accounted for 26.69% of the phenotypic variance in total, and 5.15% and 20.06% alone, indicating a significant role of the additive effects on the variation of glyphosate resistance ( Supplementary Table S3 ).

Since RGlu6 explained the highest phenotypic variation (R² = 19.20%, Supplementary Table S2 ), we chose it as a major QTL for rice glufosinate resistance and conducted variant-function-based association analysis of the associated loci regions on RGlu6 (see Methods). By analyzing the function of variants within 500 kb upstream and downstream of the peak SNP and the linkage disequilibrium (LD) blocks in that segment ( Figures 3A, B ), we analyzed the candidate genes and functional variants of RGlu6. LD decay around RGlu6 was rapid with no obvious LD blocks, where the most significant functional variant (vg0623197432) was only 6 kb away from the peak SNP ( Table 1 ). A significant glufosinate resistance difference between type A and type G of vg0623197432 variant was observed (p= 1.8×10^-10, Figure 3C ). Since variant vg0623197432 was located on the exon of LOC_Os06g39070, which encodes a UDP-glucosyl transferase (UGT), we consider it as a candidate gene for RGlu6.

Using the genotypes of the 2 kb promoter and coding region variations of LOC_Os06g39070 with 61 SNPs and 13 InDels, we categorized RGlu6 of these 421 cultivars into 5 haplotypes, namely H1 - H5 ( Figure 3D ). Haplotype H1 was predominantly present in TeJ and TrJ, H2 and H3 mainly in IndI and IndII, H4 in Aus, and H5 only in TeJ, respectively. Multiple comparison demonstrated no significant differences were among haplotypes H1-H4, while haplotype H5 exhibited significantly higher resistance to glufosinate than others, indicating beneficial allelic variations of H5 haplotype ( Figure 3E ). Thus, we classified RGlu6 into two allelic genotypes with merging the indistinct H1-H4 into RGlu^A genotype, and considering H5 as RGlu^G genotype. By analyzing the geographic distribution of these allelic genotypes ( Supplementary Table S4 ), we observed that the accessions with the RGlu6^G genotype mainly originated from Europe ( Figure 3F ), which indicates that RGlu6^G may have arisen from a mutation and was retained locally through natural selection.

Analyzing candidate genes at the RGly8 locus for rice glyphosate resistance in a rapid LD decay around RGly8, we found that vg0808951557 was the most significant functional variant located only 7 kb from the peak SNP ( Figures 4A, B ). The T to C mutation of vg0808951557 which was in the fourth exon of LOC_Os08g14850, led to an amino acid changed from serine to proline ( Figure 4C ). Since LOC_Os08g14850 encodes a resistance protein and was reported to be involved biotic and abiotic stresses, we considered it as a candidate gene for RGly8.

Seventy four SNPs and five InDels were used to analyze the haplotypes of RGly8 and seven haplotypes, namely H1 - H7 were categorized based on the glyphosate resistance and genotypes of the 421 accessions ( Figure 4D ). Haplotypes H1-H5 were predominantly present in Indica, while H6 and H7 were mainly present in Japonica. The superior C genotype of vg0808951557 was exclusively found in the H6 haplotype. Multiple comparisons revealed significant differences for glyphosate resistance among these seven haplotypes, with H6 haplotype exhibiting a significantly higher glyphosate resistance than other haplotypes ( Figure 4E ). Therefore, we designated the H6 haplotype as the allelic genotype RGly8^C , and the other haplotypes as RGly8^T . Analysis of the geographic distribution of these alleles showed that the superior allele RGly8^C primarily originated from Europe, the same as that of RGlu6^G ( Figure 4F ; Supplementary Table S4 ).

These results suggested that resistance to both glufosinate and glyphosate may have been locally selected in Europe, leading to an increase in the frequency of superior resistance alleles in rice accessions originated from that region.