The 203 Oryza species used for the subgroup classification of LM8 included 180 cultivated varieties and 23 weedy rice varieties. The O. sativa groups japonica and indica were designated as GJ and XI, respectively. Using the 3K database (Wang et al., 2018; https://registry.opendata.aws/3kricegenome/), 20 cultivated varieties were randomly selected from the GJ-adm (admixture components within japonica and indica), GJ-trp (Southeast Asian tropical), GJ-sbtrp (Southeast Asian subtropical), GJ-tmp (primarily East Asian temperate), XI-1A (East Asia), XI-1B (modern varieties of diverse origins), XI-2 (South Asia), XI-3 (Southeast Asia), and XI-adm groups. The NIP and R498 genomes were selected as the representative japonica and indica reference genomes for the whole-genome comparison with LM8.

The SNP databases for the 180 cultivated varieties are available online (SNP-Seek; http://snp-seek.irri.org). Information for 10 straw hull (SH) and 10 black hull awned (BHA) US weedy rice lines and 3 Chinese weedy rice lines was downloaded from the National Center for Biotechnology Information (NCBI) database (SRR4334499). The NIP, R498, and Arabidopsis genome data at the chromosome level were obtained from online databases (www.rice.uga.edu/, www.mbkbase.org/R498/, and www.ncbi.nlm.nih.gov/data-hub/genome/GCF_000001735.4/, respectively). Details regarding the LM8 genome are available in the NCBI BioProject database (accession number PRJNA754271).

Variant SNPs and insertions/deletions (InDels) were detected using the GATK (4.2.2.0) software, which was followed by (1) identification of active regions; (2) assembly of plausible haplotypes; (3) estimation of the per read likelihoods using the PairHMM algorithm; and (4) determination of the sample genotype. First, the GATK-HaplotypeCaller tool (versions 4.2.2.0) was used to analyze the bam files to obtain the sample GVCFs, after GATK-GenotypeGVCFs were used to genotype the variants with CombineGVCFs. Second, raw variants generated from the calling steps were filtered to identify the high-quality variants using the following parameters of GATK VariantFiltration: QD <2.0, MQ <40.0, FS >60.0, MQRankSum <−12.5, ReadPosRankSum <−8.0, and SOR >4.0. After combining the 3K-SNP data, the 5.12 Mb raw variants were filtered again (mis <10% and maf >5%) to generate the 4.03 Mb variant files (Table S1).

On the basis of the whole-genome SNPs among the 201 Oryza species included in this study, a neighbor-joining phylogenetic tree was constructed using the Treebest (versions 1.9.2) software, with 1,000 bootstrap replicates. As previously described (Lee et al., 2014), a typical method to construct trees has been: 1) calculating p-distance from all SNP data between two samples, 2) making the p-distance matrix for all samples, 3) constructing and drawing the phylogenetic tree image. The principal component analysis (PCA) of the SNPs was completed using the GCTA software to cluster the principal components into different subsets according to the differences in the individual genome SNPs. The population structure was analyzed using the ADMIXTURE software. Low cross-entropy values reflected high-quality runs. Independent runs were performed for each simulated K value (from 2 to 10). The K value for which the cross-entropy curve reflected a sensible model was chosen (i.e., based on the likelihood value).

Two comparisons were completed (LM8 versus NIP and LM8 versus R498) using the MUMmer (4.0.0rc1) software with the following parameters: nucmer -l 50 -c 100 -mum; delta-filter -i 90 -l 100 -1. According to the results of the collinearity analysis, presence-absence variations (PAVs) and structural variations (SVs) were identified using the SVUM and SYRI programs, respectively. In the final sequence variation files, variant sequences longer than 50 kb were retained, whereas variant sequences in the gap region were eliminated during the original variation test. Each SV was calculated for 1 Mb windows to generate the final Circos results.

To identify the WRKY and NB-LRR genes, the WRKY and TIR, CC, NBS, or LRR amino acid motifs were searched. The WRKY and NLR sequences in the LM8, NIP, and Arabidopsis genomes were obtained from database (accession number PRJNA754271 in the NCBI BioProject, www.rice.uga.edu, and www.ncbi.nlm.nih.gov/data-hub/genome/GCF_000001735.4/, respectively). The Molecular Evolutionary Genetics Analysis (MEGA; v11) program was used to construct a maximum-likelihood gene family phylogenetic tree, with 1,000 bootstrap replicates. Next, the WRKY and NB-LRR genes were mapped to 12 chromosomes using MapChart (2.3).

In 2021, Li et al. published the LM8 high-quality genome sequence and revealed LM8 is a japonica-type weedy rice line. To further classify LM8 within the japonica group, its SNPs as well as those in other weedy rice lines and accessions from the 3K database selected on the basis of the 3K classifications were downloaded to construct a neighbor-joining phylogenetic tree. The indica group (XI) served as the control. The phylogenetic tree constructed according to pairwise Nei’s genetic distances showed that the 10 subpopulations could be clustered into two groups (japonica and indica) (Figure 1A). Additionally, LM8 was located in the ‘adm’ branch (Figure 1A), implying that it was admixed (i.e., between japonica and indica). Moreover, LM8 was also classified in the japonica group. To verify the classification of LM8 in GJ-adm, a PCA was performed using the genome-wide SNPs in domesticated and weedy rice varieties (Figures 1B, C, and S1). The first three principal components (PCs), which explained 36.68% of the variance (Table S2), corresponded to genetic evolutionary factors as follows: PC1 separated the japonica (GJ-adm, GJ-trp, GJ-sbtrp, and GJ-tmp) and indica (XI-1A, XI-1B, XI-2, XI-3, XI-adm, SH, and BHA) subgroups (Figures 1B, C), whereas PC2 and PC3 distinguished between the weedy rice lines. Consistent with this finding, the weedy rice lines SH and BHA formed two clusters. The centralized subpopulation classification suggested the population clustering according to 3K was accurate. In terms of LM8, PC1 revealed the obvious differences between LM8 and the cultivated (japonica and indica groups) or weedy rice varieties (Figure 1B).

Next, the STRUCTURE program and an increasing subpopulation (K) value (2 to 10) were used to assign individuals to population structures and to explain the PC1-based clustering of LM8. The final population subgroups were determined according to (1) the likelihood value of these models and (2) the previously reported classification of accessions on the basis of the 3K database. When K was set to 2, the likelihood value was highest (0.570) and most of the rice accessions were clearly divided into the indica and japonica groups (Figure S3). Additionally, 50.43% and 49.57% suggested the “misplaced” LM8 was derived from a japonica and indica mixture (Table S3). Of the nine runs for K = 9, the run with the lowest likelihood value (0.462) was optimal for grouping samples and was selected for assigning the posterior membership coefficients to each accession (Figure S3). These rice accessions were grouped into the following 11 subpopulations: SH, BHA, GJ-adm, GJ-trp, GJ-sbtrp, GJ-tmp, XI-1A, XI-1B, XI-2, XI-3, and XI-adm (Figure 1D). In addition, the indica/japonica admixture of LM8 was easily identified and was in accordance with the phylogenetic relationships, with 22.98% from XI-2, 22.86% from GJ-trp, 17.76% from GJ- sbtrp, 9.86% from GJ-tmp, 8.63% from XI-3, 8.44% from BHA, 6.59% from XI-1A, 2.36% from XI-1B, and 0.52% from SH (Table S3). Nearly 10% of the genome contained genetic remnants of weedy rice, but in terms of geographical locations, LM8 tended to be similar to the tropical or subtropical accessions. These results suggested that the LM8 admixed genetic background was derived mostly from indica/japonica hybridizations rather than from weedy rice lines. Furthermore, LM8 is a japonica-type weedy rice line.

Considering a considerable proportion of its admixture components was from XI and GJ, LM8 may be a novel rice resource with a complementary gene pool. To further explore genome-wide diversity, the PAVs and SVs in the genomes of NIP (a typical japonica variety), R498 (a typical indica variety), and LM8 were detected and examined. Next, the deletion (DELs), insertion (INSs), inversion (INVs), translocations (TRANS) and duplication (DUPs) were identified between LM8 and NIP (as reference) and between LM8 and R498 (as reference; Figure 2). No matter comparing to NIP or R498, DELs and INSs were distributed relatively evenly across 12 chromosomes in the LM8 genome (Figures 2A, B). In contrast, INVs, TRANS and DUPs were unevenly distributed across 12 chromosomes in LM8, with many of them detected on chromosomes 2, 5, 6, 7, 8, 10, 11, and 12 (Figures 2A, B). Comparisons to NIP, the higher variation density of INVs were detected in the position of Chromosome 3 (20000001-21000000), Chromosome 5 (16000001-17000000), Chromosome 6 (12000001-19000000), Chromosome 8 (6000001-8000000 and 12000001-13000000), Chromosome 11 (18000001-19000000) and Chromosome 12 (14000001-15000000; Figure 2A, Table S4). Besides, there was a sharp peak in the 5000001-6000000 region of chromosome 8 and 78 variation density of TRANSs (Figure 2A, Table S4). Also, the obvious variation density peak of DUPs were detected in the position of Chromosome 3 (20000001-21000000), Chromosome 5 (15000001-16000000), Chromosome 6 (12000001-20000000), Chromosome 7 (11000001-14000000), Chromosome 8 (5000001-15000000), Chromosome 11 (9000001-10000000), Chromosome 12 (14000001-18000000; Figure 2A, Table S4). Differently, the maximum value of TRANSs located the 12000001-13000000 position on Chromosome 12 in the comparison of LM8 and R498 group (Figure 2B, Table S5). And comparisons to R498, the distinct peak value of INVs and DUPs were detected in the 14000001-16000000, 9000001-12000000, 12000001-14000000, 14000001-15000000, 11000001-15000000 region of chromosome 6, 7, 8, 10, 11, respectively (Figure 2B, Table S5). According to the result of the above description, there were significant correlation between INVs and DUPs of distribution in both LM8 versus NIP and LM8 versus R498 comparisons.

As expected, there were fewer SVs between LM8 and NIP than between LM8 and R498. The LM8 vs NIP comparison revealed 27,125 SVs, of which 9,806 were DELs, 11,209 were INSs, 1,040 were INVs, 3,783 were TRANSs, and 1,287 were DUPs in LM8 (Figures 2A, C). Among the 36,760 SVs detected by the LM8 vs R498 comparison, 12,121 were DELs, 15,105 were INSs, 984 were INVs, 6,879 were TRANSs, and 1,671 were DUPs in LM8 (Figures 2B, C). We defined the SVs exclusive to NIP or R498 as specific SVs, whereas the SVs in both NIP and R498 were designated as common SVs. Of the 63,885 SVs revealed by the two comparisons, 9,739 were common SVs (15.24%), including 4,773 DELs, 3,608 INSs, 240 INVs, 991 TRANSs, and 127 DUPs in LM8. In contrast, 27.21% (17,386/63,885) and 42.29% (27,021/63,885) of the SVs were specific to the LM8 vs NIP and LM8 vs R498 comparisons, respectively.

The WRKY TFs have core functions affecting plant hormone crosstalk and processes influencing plant growth and development (e.g., stress responses). It is unclear whether the abundant SVs in LM8 affect WRKY gene responses to various biotic and abiotic stresses under natural conditions. To identify WRKY genes and determine their chromosomal distribution, 96 WRKY genes in the LM8 and NIP genomes and 73 WRKY genes in the Arabidopsis genome were analyzed (Figure S4). Details regarding these genes, such as their gene IDs and WRKY domains and positions, are listed in Table S6. In LM8, WRKY genes were detected on 12 chromosomes, some of which were included in four clusters comprising at least four genes (Figure 3A). These clusters were distinct from those in Arabidopsis (Figure S5), with cluster 1 on chromosome 1 (35,398,161-35,468,751), cluster 2 on chromosome 5 (28,617,903-28,933,603), cluster 3 on chromosome 11 (756,297-803,420), and cluster 4 on chromosome 12 (854,961-904,182; Table S6 and Figure 3A). Interestingly, 15 of the 18 WRKY genes in the four clusters belonged to group III. More specifically, the comparison with Arabidopsis indicated 13 of these 15 genes were clustered in subgroup IIIb (Figure 3B). These 13 WRKY genes were identified as follows: ORUFILM01g003806, ORUFILM01g003807, ORUFILM01g003809, ORUFILM01g003811, ORUFILM11g002299, ORUFILM11g002300, ORUFILM11g002301, ORUFILM11g002302, ORUFILM11g002303, ORUFILM12g001262, ORUFILM12g001263, ORUFILM12g001266, and ORUFILM12g001268.

The maximum-likelihood phylogenetic tree constructed using the WRKY gene sequences in LM8, NIP, and Arabidopsis and MEGA (v11) (Figure 3B) divided the genes into subgroups Ia, Ib, IIa, IIb, IIc, IId, IIIa, IIIb, and IV, which was in accordance with previously reported classifications (Ross et al., 2007). The group IV WRKYs lack a complete zinc-finger motif (Eulgem et al., 2000). An earlier study indicated japonica contains relatively few group IV WRKY TFs (Ross et al., 2007). In another study, Xie et al. (2005) proposed that group IV genes may be incorrectly annotated because of genome sequencing errors or they might be pseudogenes that lack biological functions. Thus, the three main groups of WRKY family genes (i.e., I, II, and III) are indicated by three red circles in Figure 3B. The phylogenetic analysis suggested group I diverged into two clades (subgroups Ia and Ib), whereas group II diverged into four clades, with subgroups IIa and IIb clustered in one clade and subgroups IIc and IId in another (Figure 3B). Group III was the largest, with subgroup IIIa and IIIb in two clades (Figure 3B). One of the subgroup IIIb members in LM8 was related to a gene in NIP (Figure 3B). In the phylogenetic tree, the following three genes were clustered with Arabidopsis genes and not NIP genes: ORUFILM02g002693 (chromosome 2: 26,334,166-26,346,714), ORUFILM05g002725 (chromosome 5: 1,733,309-1,734,147), and ORUFILM05g001757 (chromosome 5: 15,083,202-15,087,922) (Figure 3B). Moreover, they were specifically detected in the LM8 genome and not in the NIP genome (Figure 3A). According to the constructed tree, ORUFILM02g002693 is most closely related to AT2G44745, whereas ORUFILM05g002725 and ORUFILM05g001757 are most closely related to AT2G23320 and AT2G40740, respectively (Figure 3B).

To clarify the pathogen resistance associated with sequence variants that may be applicable for breeding, the NBS-encoding genes in LM8, NIP, and Arabidopsis were downloaded and analyzed. The NIP genome had the most NBS genes (496), followed by the genomes of LM8 (464) and Arabidopsis (165) (Table S7 and Figure S6). The NBS genes in Arabidopsis were detected on 5 chromosomes, while the NBS genes in LM8 and NIP were on 12 chromosomes. In Arabidopsis, the NBS genes were mainly designated as CNL and TNL types (Mondragón-Palomino et al., 2002), according to Xiang et al. (2017), which followed by CN, NBS, NL, RN, RNL, and TN types (Figure 4). With the exception of one RNL-encoding gene (ORUFILM12g000772) in LM8, the NBS genes in the LM8 and NIP genomes were mainly CN, CNL, NBS, NL, and TN types (Figure 4). These results are consistent with the findings of an earlier study, in which CN- and CNL-encoding genes were identified in Oryza species (i.e., non-TNLs) (Zhou et al., 2004). In the LM8 and NIP genomes, there were three TN-encoding genes that lacked the LRR domain (Figure 4). Compared with the NIP genome, the LM8 genome had fewer CN-encoding genes, but more CNL-encoding genes (Figure 4). To elucidate the evolutionary relationships among the predicted NBS genes, a phylogenetic tree were constructed based on the fact the diversity of the protein domain was generally consistent with that of the NBS region (Figure 5). The tree revealed a lack of specific clustering, with the NBS genes divided into almost 200 groups. One of the clades included most of the CN-encoding genes in LM8 and NIP and the TNL-encoding genes in Arabidopsis, implying they may share a common origin and may be functionally similar. Furthermore, the phylogenetic tree indicated the NBS genes in Arabidopsis and rice (LM8 and NIP) formed distinct clusters, unlike the WRKY genes (Figures 3, 5).