Identification of spontaneous mutation for broad-spectrum brown planthopper resistance in a large, long-term fast neutron mutagenized rice population

Background The development of rice varieties with broad-spectrum resistance to insect pests is the most promising approach for controlling a fast evolving insect pest such as the brown planthopper (BPH). To cope with rapid evolution, discovering new sources of broad-spectrum resistance genes is the ultimate goal. Results We used a forward genetics approach to identify BPH resistance genes in rice (Oryza sativa L.) using double digest restriction site-associated DNA sequencing (ddRADseq) for quantitative trait loci (QTL)-seq of the backcross inbred lines (BILs) derived from a cross between the BPH-susceptible cultivar KDML105 and BPH-resistant cultivar Rathu Heenati (RH). Two major genomic regions, located between 5.78–7.78 Mb (QBPH4.1) and 15.22–17.22 Mb (QBPH4.2) on rice chromosome 4, showed association with BPH resistance in both pooled BILs and individual highly resistant and susceptible BILs. The two most significant candidate resistance genes located within the QBPH4.1 and QBPH4.2 windows were lectin receptor kinase 3 (OsLecRK3) and sesquiterpene synthase 2 (OsSTPS2), respectively. Functional markers identified in these two genes were used for reverse screening 9323 lines of the fast neutron (FN)-mutagenized population developed from the BPH-susceptible, purple-pigmented, indica cultivar Jao Hom Nin (JHN). Nineteen FN-mutagenized lines (0.24%) carried mutations in the OsLecRK3 and/or OsSTPS2 gene. Among these mutant lines, only one highly resistant line (JHN4) and three moderately resistant lines (JHN09962, JHN12005, and JHN19525) were identified using three active, local BPH populations. The 19 mutant lines together with three randomly selected mutant lines, which did not harbor mutations in the two target genes, were screened further for mutations in six known BPH resistance genes including BPH9, BPH14, BPH18, BPH26, BPH29, and BPH32. Multiple single nucleotide polymorphisms (SNPs) and insertion-deletion (Indel) mutations were identified, which formed gene-specific haplotype patterns (HPs) essential for broad-spectrum resistance to BPH in both BILs and JHN mutant populations. Conclusion On the one hand, HPs of OsLekRK2–3, OsSTPS2, and BPH32 determined broad-spectrum resistance to BPH among RH-derived BILs. On the other hand, in the JHN mutant population, BPH9 together with seven significant genes on chromosome 4 played a crucial role in BPH resistance. Electronic supplementary material The online version of this article (10.1186/s12284-019-0274-1) contains supplementary material, which is available to authorized users.


Background
Rice (Oryza sativa L.) is an important cereal crop that feeds almost half of the world's population (Mohanty 2013), and is mostly grown in Asia (Muthayya et al. 2014). Insect pests, such as brown planthopper (BPH), white-backed hopper, green leafhopper, stem borer, and gall midge, cause severe damage to the rice crop across the world (Ane and Hussain 2016). Among these, BPH (Nilaparvata lugens Stål) is one of the most economically damaging insect pests of rice, and causes substantial losses in yield in the major rice growing areas in Asia each year.
During plant development from the seedling to the reproductive stage, BPH sucks the phloem sap, causing whole plant senescence called hopper burn (Lou et al. 2005;Dale 1994). A significant threat to rice production was reported in the early 1970s, when new high yielding varieties and extensive use of fertilizers and pesticides caused the rapid evolution of pesticide-tolerant biotypes of insect pests. These newly introduced rice cultivars, lacking broad-spectrum resistance, exhibit short life spans. Therefore, the best approach to prolong BPH resistance is to develop rice cultivars with durable resistance against multiple BPH biotypes. More than 30 BPH resistance genes or quantitative trait loci (QTL) have been genetically mapped (Ling and Weilin 2016). However, only eight have been cloned, including BPH14 (Du et al. 2009), BPH26 (Tamura et al. 2014), BPH3 (Liu et al. 2014), BPH29 , BPH 9 (Zhao et al. 2016), BPH18 (Ji et al. 2016), BPH32 (Ren et al. 2016), and BPH31 (Prahalada et al. 2017).
Rathu heenati (RH), a Sri Lankan landrace, is one of the most strongest sources of durable resistance against BPH in Southeast Asia. This landrace has been widely used as a BPH resistance donor in breeding programs in the Great Mekong Sub-region, including Thailand, Cambodia, Myanmar, and Laos (Wang et al. 2012). Several BPH resistance loci have been identified in RH, including BPH3 (Sun et al. 2005) and sesquiterpene synthase II (OsSTPS2) (Kamolsukyunyong et al. 2013) on chromosome 4, and BPH32 (Jairin et al. 2007) on chromosome 6. The BPH3 locus has been cloned from RH, and comprises the lectin receptor kinase (OsLecRK1-3) gene cluster (Liu et al. 2014). The CDSs of OsLecRK1, 2, and 3genes are 2442, 2436, and 2436 bp, in length, respectively (GenBank accession nos. KF748957,KF748965,and KF748973,respectively). Sequence comparison between rice varieties showing differences in BPH resistance identified 17, 7, and 31 polymorphic sites in OsLecRK1, OsLecRK2, and OsLecRK3, respectively (Liu et al. 2014). The CDS of the OsSTPS2 gene cloned from RH is 1518 bp in size (GenBank accession no. KC511027), and contains a 21-bp in-frame insertion in exon 5, adding 7-aa to the encoded protein, compared with the protein sequence in the BPH susceptible variety KDML105 (Kamolsukyunyong et al. 2013). The BPH32 gene contained 20 single nucleotide polymorphisms (SNPs) and several insertion-deletion (Indel) mutations in the 585-bp CDS (Ren et al. 2016).
Rice varieties with durable resistance to BPH are key to sustainable rice production in Asia because of rapidly evolving BPH biotypes and high input farming practices. Natural variation in BPH resistance is limited, as rice is an obligate self-pollinating crop. Chemical and physical mutagenesis has been used to induce mutations and create novel sources of resistance in various genotypes, which are then used as germplasm for the development of new resistant varieties. In rice, induced mutations have been reported using gamma irradiation (Till et al. 2007;Wu et al. 2005), ethyl methanesulfonate (EMS) (Wu et al. 2005), fast neutron (FN) (Wu et al. 2005;Ruengphayak et al. 2015), and ion beam (Yamaguchi et al. 2009). The effects of FN mutagenesis on structural changes in the genome have been reported in soybean (Bolon et al. 2011(Bolon et al. , 2014Campbell et al. 2016), Arabidopsis (Li et al. 2002), and rice (Ruengphayak et al. 2015;Li et al. 2017). Deletions are the main structural rearrangements induced by FN mutagenesis (Bolon et al. 2011(Bolon et al. , 2014Li et al. 2002). In soybean, FN treatment has been shown to induce chromosomal rearrangement near the target gene (Campbell et al. 2016). Induction of single nucleotide variation (SNV) by FN has been reported in the indica rice cultivar Jao Hom Nin (JHN) (Ruengphayak et al. 2015) and japonica rice cultivar Kitaake (Li et al. 2017). In Kitaake, SNV is the most abundant mutation, accounting for 48% of the total number of mutations, and 58% of the SNVs are located within rice genes (Li et al. 2017). Structural variation in CDSs may create useful mutations; for example, tandem duplication in the waxy gene of rice (Wanchana et al. 2003), structural rearrangement in NAP1 gene in the gnarled trichrome mutant of soybean (Campbell et al. 2016), and haplotype change in OsFRO1 in a rice mutant tolerant to iron (Fe) toxicity (Ruengphayak et al. 2015).
In the current study, we combined QTL-seq (Takagi et al. 2013) and double digest restriction site-associated DNA sequencing (ddRADseq) (Peterson et al. 2012) to scan the entire rice genome for the identification of genes responsible for BPH resistance. The derived SNPs were used for reverse screening of a large FN-derived mutant population. We demonstrated that haplotype patterns (HPs) of OsLecRK2, OsLecRK3, OsSTPS2, and BPH32 contributed significantly to broad-spectrum BPH resistance in RH-derived backcross inbred lines (BILs). Additionally, HPs of BPH9 together with HPs of seven significant genes on chromosome 4 played a crucial role in broad-spectrum BPH resistance in the FN-derived mutant population.

Distribution of BPH-damage AUC in BILs
Based on the average area under the curve (AUC) of BPH damage scores from six BPH populations (Additional file 1: Table S1), the distribution of BILs, KDML105, and RH are shownin Fig.1. According to the AUC scores, KDML105 (AUC = 34.6) and RH (AUC = 11.3) parental lines were marked as the most Fig. 1 Frequency distribution of average AUC of damage scores on BILs from 6 BPH populations susceptible and resistant genotypes, respectively, whereas the AUC scores of their progeny ranged from 31.8 to 13.9. On the extreme tails, 16 and 19 BILs selected for susceptible and resistance pools have the average AUCs ranged from 23.4-31.8 and 13.9-18.4, respectively (Additional file 2: Figure S1). Also, within each pool, 4 susceptible and 5 resistance BILs selected as extreme individuals have the average AUCs ranged from 28.0-34.3 and 13.9-16.0, respectively (Additional file 2: Figure S2). These selected BILs were used in the next part for QTL-seq/ddRAD analysis.

Genome-wide SNP index analysis and SNP mining in BILs
A total of 33.75 million paired-end reads were generated from the sequencing of 11 ddRADseq libraries (NCBI accession nos. SRR7851264 to 69, SRR7851271, SRR7851273 to 76), and SNPs in each library were called for QTL-seq analysis (Additional file 1: Tables S2, S3). Using a 2 Mb-sliding window for genome-wide scanning, 37,212 windows were allocated to the reference genome. Based on the maximum average ΔSNP-indexes (max ΔSNP) of individuals and pools, two QTLs for BPH resistance were mapped on genomic windows between 5.78-7.78 Mb (QBPH4.1) and 15.22-17.22 Mb (QBPH4.2), on the short and long arms of chromosome 4, respectively (Fig. 2a, b). Linkage disequilibrium analysis of polymorphic SNPs in QBPH4.1 and QBPH4.2 in the BIL population revealed three haplotype blocks (Fig. 2d). The first block covered a 396-kb genomic region containing F-box118 (LOC_Os04g 11450), F-box119 (LOC_Os04g11660), resistance protein LR10 (LOC_Os04g11780), and an unknown protein (LOC_Os04g12110). The second haplotype block contained OsLecRK3 (LOC_Os04g12580), whereas the third block contained OsSTPS2 (LOC_Os04g27430).
A total of 107 and 94 SNPs were identified in QBPH4.1 and QBPH4.2 genomic windows, respectively. In the QBPH4.1 window, 31 SNPs were identified in 11 genes, resulting in 22 non-synonymous mutations and one premature stop codon (Table 1). In the F-box119 protein, two amino acids change from aspartic acid (D) to glutamic acid (E) and lysine (K) to asparagines (N) were found at position 253 and 279 respectively. In the resistance protein LR10, six amino acid changes were identified at positions 371 (isoleucine [I]  located at amino acid position 387 in the protein encoded by the terpene synthase gene (LOC_Os04g 27720) (Table 2). In total, five SNPs and one indel marker were developed from the QBPH4.2 region. Additionally, the region on chromosome 6 in RH spanning 0.03-2.03 Mb, which harbors the previously cloned BPH resistanceconferring BPH32 gene, (Ren et al. 2016), showed the max ΔSNP. This window contained nine non-synonymous mutations in eight annotated genes. Unfortunately, no SNP was identified in the BPH32 gene using ddRADseq data (Additional file 1: Table S4). Therefore, the SNP marker identified previously in BPH32 by Ren et al. (2016) was used as a marker in this study. All of the SNPs identified in QBPH4.1, QBPH4.2, and BPH32 were used for QTL detection in the BIL population using single marker analysis.

Roles of OsLecRK3 and OsSTPS2 in broad-spectrum BPH resistance
To determine candidate genes that were exclusively involved in BPH resistance, single marker QTL analysis was performed using 28 candidate SNPs and BPH damage AUC of the selected BILs infested by five BPH populations. The results revealed a significant association between the BPH resistant QTL and OsLecRK3SNP in all five BPH populations with phenotypic variance ranging from 16% to 44% (Table 3). The second most important marker was the 21-bp indel in OsSTPS2, which accounted for 8%-30% of the total phenotypic variance in BPH resistance Specific associations were detected between some SNPs and specific BPH populations. The SNP markers in F-Box118 protein, oxidoreductase (LOC_Os04g26910), and BPH32 showed major associations, while SNP markers specific to the Verticillium wilt resistance (VWR) protein (LOC_Os04g28210) showed a minor association, specifically with the Huai Thalaeng (HTL) BPH population. With the Kamphaeng Phet (KPP) BPH population, SNP markers from the LR10 and unknown protein (LOC_Os04g12110) accounted for 4.5-7% of the total phenotypic variation, whereas terpene synthase gene explained approximately 7.6% of the total phenotypic variation in BPH resistance with the Phitsanulok (PSL) BPH population. Therefore, QTL-seq analysis revealed the two most important regions of chromosome 4, which contributed the RH-derived broad-spectrum BPH resistance in BILs.
Reverse screening of JHN mutant population for BPH resistance using SNP markers from OsLekRK3 and OsSTPS2 To confirm the functional impact of SNP and 21-bp indel markers in OsLecRK3 and OsSTPS2 genes, respectively, we screened 9323 JHN mutant lines for polymorphisms in OsLecRK3 and OsSTPS2 genes. This analysis identified 11   Table 4 and Additional file 2: Figure S3). The UBN BPH population was used as the representative of BPH from the rainfed rice farming area in the Northeastern region of Thailand while the TPY and CNT BPH populations were used as the representative of BPH from the irrigated rice farming area in the Central region of Thailand. Four mutant lines (JHN4, JHN12005, JHN19525, and JHN09962) showed resistance to BPH in the BPH infestation test (Table 4). The mutant line JHN4, like RH, showed resistance to all three BPH populations. The mutant line JHN12005 showed strong resistance to CNT and TPY populations but only moderate resistance to UBN population. The third mutant line, JHN19525, showed strong resistance to TPY and moderate resistance to CNT and UBN. The last mutant line, JHN09962, was resistant to TPY and UBN but susceptible to CNT. Interestingly, JHN19874, JHN07766, and seven other mutant lines carried the same mutations in OsLecRK3, and OsSTPS2 genes as RH but were susceptible to all three BPH populations. This observation indicates that OsLecRK3 and OsSTPS2 genes may not fully explain the broad-spectrum resistance to BPH in the FN-induced mutants based on the polymorphic markers, or either RH or FN-mutants may harbor additional R gene besides these two genes.
Sequence analysis of the OsLecRK3 gene, the target of reverse screening, showed that four resistant mutant lines (JHN19525, JHN09962, JHN12005, JHN4), three susceptible mutant lines (JHN19578, JHN19874, JHN11183), and RH shared nine synonymous amino acid substitutions ( Fig. 3 and Additional file 3: Data S1). Furthermore, extending SNP mining into the adjacent genes (OsLecRK1 and OsLecRK2) revealed similar trends as OsLecRK3; the amino acid sequence of resistant and susceptible lines exhibited greater homology with RH than with the WT (Fig. 3). In OsLecRK1, 16 amino acid substitutions, and three amino acid insertions were shared between the resistant lines, susceptible lines, and RH (Additional file 3: Data S2), whereas in OsLecRK2, six amino acid substitutions were shared among these genotypes. Interestingly, novel amino acids were identified at positions 386,407,429,433,495,502, and 534 of the OsLecRK2 protein, with mutants differing from both the WT and RH (Additional file 3: Data S3). These seven novel amino acids, found only in OsLecRK3-positive mutants, ruled out the possibility of contamination or outcrossing with RH during seed propagation. In particular, amino acid change at positions 534 was located within the protein kinase domain. Greater homology was detected in the nucleotide sequence of OsLecRK1-3 and the encoded amino acid sequence among mutants and RH compared with the WT. However, only four mutant lines were BPH resistant. Thus, we speculated that other genetic factors contributed to the BPH resistance of the mutants. Furthermore, the mutant line JHN4 was identified as a promising new BPH resistance donor for pyramiding genes to obtain broad-spectrum BPH resistance in rice.

Mutant lines carrying BPH32 are susceptible to local BPH populations
To determine whether BPH32 imparts BPH resistance in the JHN mutant population, the SNP in BPH32 was used for reverse screening 9323 randomly selected mutant lines in the M 6 generation. A total of 16 mutant lines carrying the resistant allele of BPH32 and susceptible alleles of OsLecRK3 and OsSTPS2 were identified. To determine if these BPH32-positive mutant lines showed broad-spectrum resistance, seeds of all 16 lines were retrieved from the mutant seed stock for BPH-infestation test. Unfortunately, only five lines could be germinated and exposed to three BPH populations. As a result, all the lines tested were found to be susceptible to all BPH-populations tested (Additional file 1: Table S5). To find out if the selected mutants carried new SNPs in the BPH32, the genomic sequence of the BPH32 gene of JHN wild-type, selected mutants, and RH were compared. However, no new polymorphic SNPs were found (data not shown), suggesting that the gene had no contribution to BPH resistance in the mutant population. These data suggested that the genetic basis of broad-spectrum BPH resistance in the JHN mutant is complex and different from that in RH-derived BIL populations. Thus, genes other than OsLecRK, OsSTPS2, and BPH32 likely confer resistance to BPH in the JHN mutant population.
(See figure on previous page.) Fig. 2 The genome-wide scanning of average ΔSNP-indexes graphs plotted along the rice genome. a The individual extreme BILs, b the pooled extreme BILs. c Candidate QTL locations were identified on rice chromosome 4 with statistical confidence intervals under the null hypothesis of no QTL (P < 0.05). d Linkage disequilibrium (LD) plot for SNPs in the QBPH4.1 and QBPH4.2 regions. The D' values (normalized linkage disequilibrium scores) for SNPs covered the LD blocks were shown. The SNPs were genotyped in the BIL population by KASP genotyping      Table S6). These SNP markers were used to genotype BPH resistant and susceptible mutants to identify HPs associated with broadspectrum BPH resistance in mutant lines.
Among the 15candidate genes identified in QBPH4.1 and QBPH4.2 windows and six known BPH resistance genes, 62 SNPs were identified, with 1-8 polymorphic SNPs per gene (Additional file 4: Table S7). No SNPs were identified in BPH14 or BPH32 gene, which implied that the two genes were fixed in the selected mutants. The WT HP was designated as HP1, and alternative HPs were then identified. The number of HPs varied from 2 to 10; two HPs were identified in the gene encoding inorganic phosphate transporter protein (LOC_Os04g 10800), and 10 HPs were identified in the OsLecRK3 gene. These mutants were evaluated for BPH damage AUC against three BPH populations. Therefore, it was possible to associate the HPs of candidate genes with BPH resistance (Table 5).
Single marker analysis of 21 candidate gene HPs from QBPH4.1, QBPH4.2, and BPH resistance genes revealed eight genes in which HPs were significantly associated with BPH resistance in the selected mutants, with phenotypic variance ranging from 13.5% to 83.3% (Table 5). Interestingly, the most significant association with BPH resistance was identified in HP2 of BPH9 with phenotypic variance in BPH resistance of 83.3%, 71.0%, and 40.6% against UBN, TPY, and CNT populations, respectively. SNPs at nucleotide positions 274 (A:G, WT:mutant) (SNP id. BPH9_Chr12_22,886,067) and 258 (G:T, WT:mutant) (SNP id. BPH9_Chr12_22,886,105) in exon 1 of BPH9 generated four HPs (Additional file 1: Table S6). The A-G and G-T HPs were associated with susceptibility and resistance, respectively, to all BPH populations. By substituting the susceptible A-G with the resistant G-T HPs, the damage reducing ratios (DRRs) were 0.53, 0.46, and 0.25 against CNT, UBN, and TPY populations, respectively (Table 5). The replacement of A-G HP with G-T HP resulted in one amino acid substitution from arginine (R) in WT and susceptible mutants to glycine (G) in resistant mutants, as the mutation at nucleotide position 258 resulted in synonymous substitution. The second highest impact HP was linked to the two F-Box genes, F-Box118 and F-Box119. For F-Box119, two SNPs, resulting in non-synonymous substitutions, were identified at nucleotide positions 837 (SNP id. F-box_R04006388017) and 759 (SNP id. F-box_R04006388095) of the single exon gene, generating four HPs (Additional file 1: Table S6). These SNPs caused two amino acid substitutions from lysine (K) and aspartate (D) to asparagine (N) and glutamic acid (E), respectively. The C-G HP was associated with susceptibility, while the A-C HP was associated with resistance against TPY and H (heterozygous)-G HP was associated with resistance against CNT (Table 5). In the F-Box118 gene, three SNPs generated three HPs (Additional file 1: Table S6). These SNPs caused three non-synonymous amino acid substitutions from glutamine (Q), asparagine (N), and leucine (L) to glutamic acid (E), tyrosine (Y), and valine (V), respectively. The C-A-C and G-T-G HPs were associated with susceptible and resistance, respectively, to the TPY population. Consequently, the DRR of haplotype substitutions in F-Box118 and F-Box119 genes were 0.58 and 0.50, respectively, against the TPY population (Table 5).
The third most significant HPs were specifically associated with the TPY population. Three promoter-specific SNPs and the functional 21-bp indel marker created four HPs in the OsSTPS2 gene (Additional file 1: Table S6). The G-T-C-21 bp_insertion were associated with BPH resistance against the TPY population. The 21-bp insertion in the OsSTPS2 gene together with the WT HP in its promoter increased resistance to BPH, with a DRR of 0.60 against the TPY population (Table 5). Other significant genes which reduced BPH damage of TPY population were the gamma thionin protein (LOC_Os04g11165), LR10, and VWR protein with DRRs of 0.56, 0.60, and 0.69, respectively. However, a mutation in the gene encoding inorganic phosphate transporter protein increased the BPH damage AUC by DRR of 1.37 against TPY infestation.

Genes that play a dominant role in broad-spectrum BPH resistance
The BPH9 gene was considered the most significant gene in the selected mutant lines. Four HPs were identified in mutant lines. The HP2 of BPH9 had the most substantial impact on the DRR. Three resistant mutant lines (JHN4, JHN19525, and JHN09962) were identified as carrying BPH9 HP2 (Table 6). The JHN4 mutant was resistant to all BPH populations tested (CNT, TPY, and UBN). On the other hand, JHN19525 showed strong resistance to TPY and moderate resistance to CNT and UBN populations, whereas JHN09962 was resistant to TPY and UBN but susceptible to CNT. Four mutant Table 4 Genotyping data at two loci used for reverse screening of selected mutant lines and BPH damage AUC tests using 3 BPH populations. The "Ins" indicated insertion allele while "Del" indicated deletion allele of the OsSTPS2-21 bp-indel marker. Letter a, b, c, d, e, f, g indicated similarity of the traits in multiple comparison tests Five mutant lines (JHN07766, JHN05678, JHN19572, JHN19577, and JHN19578) harbored HPs at significant genes on chromosome 4 similar to those in four resistant lines (JHN4, JHN12005, JHN19525, and JHN09962); however, these five mutants showed susceptibility to all BPH populations tested because their HP on BPH9 gene was WT (HP1). This observation confirmed BPH9 as the most important BPH resistance gene in the mutant population (Additional file 1: Table S8).
Therefore, mutations in BPH9 as well as in inorganic phosphate transporter, gamma thionin, F-Box118, F-Box119, LR10, OsSTPS2, and VWR genes (all on chromosome 4) were crucial for obtaining broad-spectrum resistance in the JHN mutant population. RH, the donor of broad-spectrum BPH resistance, carried HP2 at BPH32 and HP2 or HP3 at all 7 significant genes on chromosome 4. Furthermore, HPs of significant genes for five extreme resistant BILs and susceptible BILs were reported (Table 7) and showed that OsLecRK2-3 and OsSTPS2 played significant roles in BPH resistance in these BILs, as every BPH resistant BIL carried HPs 3, 2, and 2 in OsLecRK2, OsLecRK3, and OsSTPS2 respectively. On the other hand, the BPH32 gene possibly played interactive roles depending on the HPs of OsLecRK2-3 and the OsSTPS2 genes conferring broadspectrum BPH resistance in RH and RH-derived BILs. Moreover, genes encoding inorganic phosphate transporter protein, the gamma thionin protein, F-box118, F-box119, and LR10 may not play a role in the BPH resistance of BILs since the HP at these genes in the resistance BIL 423(4) was different from that in RH. Therefore, mutations in OsLecRK2-3 and OsSTPS2 were crucial for broad-spectrum resistance to BPH in RH and BIL populations; this suggests that breeding for broad-spectrum BPH resistance must aim to accumulate multiple BPH resistance genes from chromosomes 4, 6, and 12.

The origin of the BPH resistant JHN mutant line
To investigate whether BPH resistance found in mutant lines were induced by FN or contamination during mutant population advancement, two pieces of evidence were generated. Based on 8928 genome-wide SNPs, phylogenetic analysis clearly showed four clusters related to JHN, BILs, RH, and the japonica out-group (Fig. 4). The sequence information for OsLecRK1-3 of RH are from the GenBank database (accession no. AIE56222, AIE56230, and AIE56238 respectively). Amino acid changes were highlighted in yellow, and amino acid that JHN WT differed from RH were highlighted in green. The amino acid positions that are corresponding to KASP-SNP markers used for haplotype pattern analysis of mutant lines were highlighted by lighted-gray      (Fig. 5a). The total span of genomic region introgressed from RH ranged from4.9 Mb in BIL49 to11.7 Mb in BIL24. By contrast, very small homologous sequences were identified in the resistant (JHN4) and moderately resistant (JHN09962) mutants (Fig. 5b). Additionally, mutants carried either the same haplotype as that in JHN WT or unique HP within the OsLecRK2 CDS (Fig. 5c). Four SNPs (positions 6,956,441; 6,956,628; 6,956,639; and 6,956,769) created novel amino acids (positions 495, 433, 429, and 386 respectively) in the OsLecRK2 gene of JHN4, which were not found in either RH or JHN WT (Fig. 3). Together with its distinguish seed color of dark purple, BPH resistance identified in JHN4 is not contaminated from a pollen or seed source but induced by FN.

QTL-seq/ddRADseq identified candidate genes for BPH resistance
In this study, we successfully combined QTL-seq with ddRADseq to identify genomic region harboring BPH resistance genes from rice cultivar RH. The two genomic regions, QBPH4.1 and QBPH4.2, were identified on rice chromosome 4. For the QBPH4.1 region, two reportedly candidate genes and one cloned gene for BPH resistance were localized (Kamolsukyunyong et al. 2013;Liu et al. 2014). These genes were LOC_Os04g11660 and LOC _Os04g11780, encoding F-box 119 and LR10 proteins (Kamolsukyunyong et al. 2013) and LOC_Os04g12580, encoding the OsLecRK3, one of the members of the gene cluster (OsLecRK1-3) responsible for BPH resistance in RH (Liu et al. 2014). For the QBPH4.2 region, the most interesting candidate gene is LOC_Os04g27430, encoding the OsSTPS2 gene, reportedly to play essential roles in antixenosis resistance mechanisms against BPH (Kamolsukyunyong et al. 2013). The BPH32 gene (LOC_Os06g03240) has been mapped to an area between the simple sequence repeat markers RM589 and RM588 and served as one of the primary target loci in RH as BPH resistance donor for marker-assisted backcrossing in this BIL population (Jairin et al. 2009). Our QTL-seq analysis could not identify significant SNP with statistical confidence in this gene (Fig. 2a, b). This observation may reflect the minor effect of BPH32 gene on BPH resistance in the RH-derived BIL population. By using The JHN mutant lines were analyzed with JHN-WT and RH, c Haplotype patterns at seven novels amino acids of mutant lines. SNPs at position 6,956,325 and 6,956,441 were developed as the KASP markers for HP analysis. Filled squares represented gene-specific SNPs, filled arrowheads represented intergenic SNPs, and asterisks represented the recombination breakpoint known BPH32-specific SNP (Ren et al. 2016), only 5 out of 16 BPH32-positive lines could be tested for BPH resistance in this study and found that all tested mutants were BPH susceptible (Additional file 1: Table S5). Sequence comparison confirmed that no other new mutation was found in the BPH32 gene. As such the gene is not essential for BPH resistance in the JHN-mutagenized population.
Our results confirmed that for the durable BPH resistance rice variety RH, not only single genes/ gene cluster of OsLecRKs, OsSTPS2, or BPH32 is essential for BPH resistance, but all of them working together in response to BPH infestation. As we shown in Table 3, SNP marker from OsLecRK3 and indel marker from OsSTPS2 associated with resistance to all the BPH populations tested but SNP from BPH32 associated with resistance to only one BPH population. This result may suggest that the gene cluster OsLecRKs, and the gene OsSTPS2 play a major role while the gene BPH32 play a minor role for BPH resistance in RH. As the single marker analysis may not fully explain the effect of these genes, our HP analysis of the extreme resistance and susceptible BILs may better support this speculation. It showed that HP2 of BPH32 is useful when HPs at OsLecRK2, OsLecRK3, and OsSTPS2 were 3, 2, and 2 respectively (Table 7). On the other hand, HP2 of BPH32 will not be effective when the HPs of those three genes were not in the said pattern.
Breeding for broad-spectrum resistance requires multiple gene haplotypes Improving broad-spectrum resistance creates a solid foundation for cultivars with durable resistance against highly variable insects such as BPH. To date, at least 25 major BPH genes located as gene clusters on chromo-somes3, 4, 6, and 12 have been identified (Jing et al. 2017). Durable BPH resistance comprises the recognition of herbivore/damage-associated molecular patterns (HAMPs/DAMPs), and the activation of resistance genes encoding the coiled-coil nucleotide-binding leucine-rich repeat (CC-NB-LRR) domain proteins (Jing et al. 2017).To develop new rice varieties with durable BPH resistance, molecular markers must be developed for effective pyramiding of genes with strong resistance to BPH. Pyramiding BPH14 (a CC-NB-LRR protein) and BPH15 (a LecRK protein) enabled the successful development of a hybrid rice cultivar with broad-spectrum BPH resistance (Hu et al. 2012). In this study, two major QTLs on rice chromosome 4 (QBPH4.1 and QBPH4.2) harboring OsLecRK genes (Liu et al. 2014) and OsSTPS2 genes (Kamolsukyunyong et al. 2013), respectively, and a minor QTL harboring BPH32 on chromosome 6 (Jairin et al. 2009;Ren et al. 2016) were shown to be involved in broad-spectrum BPH resistance derived from RH. Pyramiding QBPH4.1, QBPH4.2, and BPH32 into KDML105 by marker-assisted backcrossing has successfully improved broad-spectrum BPH resistance against diverged BPH populations (Vanavichit et al. 2018). Reverse genetics using the JHN mutant population revealed comprehensive mutations within the two gene clusters on chromosome 4 together with BPH9, enabling the induction of a mutant line exhibiting broad-spectrum BPH resistance. We have previously demonstrated the role of insect-inducible volatile factors, including sesquiterpenes and monoterpenes (Pitija et al. 2014), in antixenosis mechanisms supporting durable resistance against BPH populations. In this study, we identified the role of at least eight genes harboring 22 SNPs in broad-spectrum BPH resistance in JHN mutant lines (Table 5). These mutant lines can be used as the donor in the breeding program. However, pyramiding these eight genes or 22 SNPs into one variety using marker-assisted analysis may not be efficient. Our results showed the organization of these target SNPs into eight gene-specific HPs. These results would empower breeders to design a more efficient gene pyramiding program for developing rice cultivars with strong BPH resistance.

Genome background relationship of JHN and RH
Based on our phylogenetic analysis using SNPs derived from ddRAD, it was clearly illustrated that JHN, BILs, and RH were separated into different clusters. While RH was grouped with Basmati370, N22, and FR13A from South Asia, JHN and 22 JHN M 6 mutant lines were grouped together. Also, the RH-derived BILs were grouped with KDML105, the recipient parent. The phylogenetic analysis reflected no relationship on the genome background between RH and JHN mutants and thus harboring a different set of BPH resistant genes in selected JHN mutant lines is possible. As reported by Zhao et al. (2016), RH contained HP1 susceptible allelotype at BPH9 gene, and we also showed that the BPH9 (HP1) was not essential for BPH resistance in RH-derived BILs (Table 7). On the other hand, as we showed in Table 5, the HP2 of BPH9 was associated with BPH resistance of mutants in all three BPH populations tested by reducing BPH damage as DRR of 0.535, 0.245, and 0.462 for CNT, TPY, and UBN respectively.
To confirm if the HP2 of BPH9 is genuinely crucial for BPH resistance in JHN mutant population, the 9313-mutant-lines were screened for HP2 mutation of BPH9 gene. As a result, there were three more mutant lines found to harbor HP2 of BPH9 gene. These mutant lines were JHN19727, JHN22247, and JHN22995. HPs of 21 candidate gene from QBPH4.1, QBPH4.2, and BPH resistance genes of these mutant lines were also characterized (Additional file 4: Table S7). Multiple BPH resistance of these lines will be evaluated in the future.

Rapid induction of spontaneous mutations using FN mutagenesis
FN mutagenesis typically generates deletions and chromosome rearrangements at a genome-wide scale (Gilchrist and Haughn 2010). The induction of multiple nucleotide substitutions is another essential characteristic of FN-induced mutagenesis (Belfield et al. 2012). In the model rice cultivar Kitaake, FN at 20 Gy resulted in single base substitution (SBS), deletion, insertion, translocation, and tandem duplication in the M 2 population (Li et al. 2017). Because of the high frequency of deletions, loss-of-function mutations are the primary cause of mutation detection in most FN-mutagenized populations (Li et al. 2001(Li et al. , 2002. In the M 2-3 population developed from FN mutagenesis of Kitaake, SBSs account for 47.5% of all mutations (Li et al. 2017). Similarly, in JHN mutant population generated using FN mutagenesis, SBSs, deletions, and insertions accounted for 56%, 23%, and 21% of the total mutations, respectively (data not shown). These data suggest that the genetic basis of mutagenesis was similar in the JHN and Kitaake mutant populations.
Considering the high homology of the induced mutations in the QBPH4.1, QBPH4.2, and BPH9 to RH and BPH donors, we believe that the genetic variation detected in M 6 JHN mutant population was caused by FN-induced spontaneous mutations. Our FN-mutagenized population survived two significant forces during development. First, during the early phase, a major loss of lethal mutants may have eliminated all large chromosomal rearrangements from the mutagenized population. Second, the JHN mutant population was supposedly exposed to environmental stresses during generation advancement by self-pollination in an organic paddy field in Thailand. This exposure of an FN-mutagenized population to biotic and abiotic stresses in an open paddy field may impose natural selection and homologous recombination during the prolonged generation advancement time frame. In Arabidopsis, abiotic stresses lead to heritable changes in the frequency of recombination, point mutation, and microsatellite mutation (Yao and Kovalchuk 2011). Additionally, the rate of homologous recombination is strongly affected by high temperature, short day length, and moderate change in environmental conditions (Boyko et al. 2005, Kloosterman et al. 2013). On an average, 9.6 spontaneous mutations per line in the M 2-3 Kitaake rice population contained 7.6 SBSs and two small indels ). On the other hand, in the M 6 JHN mutant population, only 4% of SBSs were novel, based on the comparison with OrzaSNP database. Therefore, higher mutation rate in the JHN mutant population could be due to 1) synergistic effects of selection against lethal mutations at the beginning of the generation advancement, 2) several generations of genetic recombination during long-term self-pollination, and 3) natural selection during generation x advancement. Spontaneous mutations may arise from local genomic lesions and repair, resulting in SBSs and single base frame shift mutations, which occur in a non-random manner at mutation hotspots throughout the genome (Maki 2002). The induced broad-spectrum BPH resistance identified in JHN4 may represent several mutation hotspots induced by FN mutagenesis and natural exposure to local BPH populations during long-term generation advancement via self-pollination.
Transposable elements (TE) are also sensitive to irradiation and environmental stresses. Stress-responsive Arabidopsis mutant lines acquired exapted TE-derived genes coding for new proteins, thus gaining new roles for host adaptation to stressful environments, such as high phosphate and arsenic levels, salinity, and freezing (Hoen and Bureau 2015;Joly-Lopez et al. 2017). In mutagenized Kitaake rice, FN-induce TE mutations at 58.6% compared to 25.7% in the same flanking sequence tag population (Hong and Jung 2018). It would be interesting to further explore the possibility whether FN-induced exapted TE transposition in the JHN mutant population. Thus, FN-mutagenized population contains newly induced suites of SBSs and indels as a basis for creating new genetic variability.

Conclusion
In this study, we used ddRADseq with QTL-seq analysis to identify SNPs and candidate resistance genes associated with BPH resistance in RH-derived BILs. Two major genomic regions associated with broad-spectrum BPH resistance were localized between 5.78-7.78 Mb (QBPH4.1) and 15.22-17.22 Mb (QBPH4.2), forming three linkage disequilibrium (LD) blocks on the rice chromosome 4. Twenty-one significant SNPs were organized into three haplotype blocks: two blocks in QBPH4.1 and one block in QBPH4.2. Functional markers associated with OsLecRK3 and OsSTPS2 genes were used for reverse screening of 9323 FN-mutagenized lines generated from the BPH susceptible rice cultivar JHN in the M 6 generation. Twenty-two mutants, including 19 identified mutants, were evaluated for BPH resistance using three active BPH local populations representing important rice growing areas in Thailand. As a result, one resistant mutant (JHN4) and three moderately resistant mutants (JHN12005, JHN09962, and JHN19525) were identified. Further screening with six known BPH resistance genes revealed that BPH9 was dominant to OsLekRK2-3 and OsSTPS2 genes in reducing BPH damage in the mutant population. On the other hand, OsLekRK2-3, OsSTPS2, and BPH32 determined the broad-spectrum BPH resistance in RH-derived BILs. Significant gene-specific HPs involved in broad-spectrum BPH resistance in both BILs and JHN mutant population were identified. We demonstrated that long-term FN mutagenesis is a useful tool for generating not only novel but also a natural genetic variation for functional genetics and molecular breeding. Together, our data suggest that RH and the BPH-resistance mutant lines can be used as the donors in a breeding program to improve broad-spectrum resistance of rice crop against diverged BPH populations. Four gene-specific haplotypes including OsLecRK2, OsLecRK3, OsSTPS2, and BPH32 are needed in the breeding program using RH as a donor. On the other hand, eight gene-specific haplotypes including BPH9, inorganic phosphate transporter, gamma thionin, F-Box118, F-Box119, LR10, OsSTPS2, and VWR genes are needed in the breeding program using the BPH-resistance mutants as a donor.

Plant materials
The BILs were developed from a cross between KDML105 and the BPH resistant donor, RH, by backcrossing using marker-assisted selection (Jairin et al. 2009). The F 1 plants were backcrossed with the recurrent parent. BPH resistant BC 1 and BC 2 plants were selected, and the BPH3-linked marker RM589 on chromosome 6 where BPH32 was localized (Ren et al. 2016), was used to generate BC 2 F 1 and BC 3 F 1 generations, respectively. Two individual BC 3 F 3 plants were developed from the resistant BC 3 F 2 progenies (n = 2343) that were heterozygous at the linked marker on chromosome 6. A total of 105 BC 3 F 5 BILs were generated by self-pollination of selected heterozygous BC 3 F 3 plants. These BILs were infested with six local BPH populations collected from intensively irrigated and rain-fed rice production areas. BILs those were highly resistant (19 and 5 for pool and individuals respectively) and susceptible (16 and 4 for pool and individuals respectively) to BPH were identified and used for QTL-seq/ ddRADseq analysis both as pools and as individuals.
The FN-mutagenized population was developed from the rice cultivar JHN, a photoperiod insensitive, semidwarf purple rice variety (Ruengphayak et al. 2015). Contamination of the JHN population can be readily recognized based on the purple grain color. Initially, 100,000 purified seeds derived from a panicle-to-row of a single JHN plant were treated with 33 Gy of FN by the Office of Atoms for Peace, Thailand, in 2012. Starting with the first generation (M 1 ), the whole mutant population was field grown; however, with the impact of FN on germination, more than 70% of the population was lost in subsequent generations. In the M 4 generation, 21,024 lines were identified as stable. Generation advancement of the mutant population followed the same procedure as that used for the development of a random inbred line population, but with slight modifications. For each line, eight plants were randomly selected, and panicles were individually bagged with a pollination glassine envelope, before pollination, and until seed set. Seed-containing panicles were processed individually, and stored individually or as pools. In total, 9313 M 5-6 mutant lines were used in this study. The reverse screening was conducted using three selected polymorphic SNP/indel markers in OsLecRK3, OsSTPS2, and BPH32 genes. The M 5 mutants that showed the same alleles as RH for these markers were selected, and M 6 families derived from these M 5 mutants were used for the BPH infestation test.

BPH populations and screening
Six local BPH populations collected from six critical rice-growing provinces in Thailand (HTL, KPP, Nan [NAN], PSL, TPY, and UBN) were used for BPH damage screening of BILs while three BPH population (CNT, TPY, and UBN) were used to evaluate the selected mutants using a modified standard seedbox screening method (SSBS) (Heinrichs et al. 1985).
The BILs and candidate mutant lines were assessed for BPH resistance at the seedling stage. Three replications of each line were under greenhouse conditions and arranged in a randomized complete block design. The second or third instar nymphs of the BPH were released on 7-10-day-old rice seedlings at the rate of 8-10 insects per plant. Seedling reaction to BPH was recorded daily as a damage score for five consecutive days, which started 7-10 days after infestation or when the susceptible check Taichung Native1 (TN1) was completely dead. The 5-day-long damage scores were presented as AUC for rating the damage caused by each BPH population. The AUC of damage rating of BPH populations were calculated using the trapezoid rule (Litsinger 1991

Genotyping-by-sequencing
Genotyping-by-sequencing of BILs and mutants was performed as an individual plant or as a pool of plants using ddRADseq (Peterson et al. 2012, Pootakham et al. 2015. Genomic DNA was isolated from leaves tissue using DNeasy Plant Mini Kit (Qiagen, Carlsbad, CA, USA), and quantified using NanoDrop ND-8000 Spectrophotometer (Thermo Scientific, Wilmington, DE, USA). Resistance and susceptible pools contained 200 ng of DNA from selected BPH resistant lines (n = 19) and susceptible lines (n = 16). The two pools were used for ddRADseq library construction. In addition, highly resistant lines (n = 5) and highly susceptible lines (n = 4) were sequenced individually as internal controls. Paired-end sequencing libraries (mean read length of 112 bp) with an insert size of approximately 200 bp were prepared. The DNA libraries were labeled using PstI-adapters containing specific 9-bp barcodes. Libraries were sequenced using the Ion Proton PITM Chip (Life Technologies, Grand Island, NY, USA), according to the manufacturer's instructions.

QTL-seq analysis
Short sequence reads from BPH resistant and susceptible pools and extreme BILs were locally aligned to the Nipponbare reference genome (Os-Nipponbare-Reference-IRGSP-1.0 pseudo-molecules), with a fragment size of 80-300 bp, using Bowtie2 Software (Langmead and Salzberg 2012). Description of sequencing data was added by using the command line AddOrReplaceReadGroups of the Picard command line tools (http://broadinstitute.github.io/picard/). Sequence data were improved in the region spanning the indel polymorphisms using the Genome Analysis Toolkit (GATK) (https://software.broadinstitute.org/gatk/), including indel positioning by GATK Realigner Target Creator and re-alignment by GATK IndelRealigner. An improved overlapping sequence of restriction associated DNA (RAD tag), was used for SNV calling using GATK UnifiedGenotyper, with a minimum confidence threshold and emit confidence score of 0.5 and 0.3, respectively. Low-quality RAD tag with a Phred quality score lower than 15 and base quality score less than 20 (Q > 20, 1:100; 90%) were excluded.
The ddRADseq libraries were sequenced at more than 6X coverage, and the identified SNVs were used for SNP-index calculation. The bi-allelic depths for the reference and alternate allele score at each position were calculated (Takagi et al. 2013). The ΔSNP-index was calculated by subtracting the SNP-index of the resistant pool from that of the susceptible pool. Alternatively, subtracting the SNP-index of the resistant lines from that of the susceptible lines was used for ΔSNP-index analysis. To identify candidate genomic regions responsible for BPH resistance, the ΔSNP-index between susceptible pool vs. resistant pool and susceptible individuals vs. resistant individuals were calculated as absolute values. The average ΔSNP-index of all 2-Mb sliding genomic windows were calculated at 10 kb increments and plotted with a statistical confidence interval.

LD analysis
LD analysis was performed in the BC 3 F 5 population (KDML105 × RH) by pairwise comparisons among kompetitive allele-specific PCR(KASP)-SNP markers distributed across the two QTLs for BPH resistance identified by QTL-seq analysis using the HAPLOVIEW software version 4.2 (Barrett et al. 2005) using the following parameters: minor allele frequency (MAF) > 0.05; Hardy-Weinberg P-value cut-off, 0; and percentage of genotyped lines > 0.50. LD was estimated using squared allele frequency correlations (r 2 ) between pairs of loci. P-value < 0.001 was used as a criteria for significant LD, the remaining r 2 values were considered as uninformative. The pattern and distribution of LD were visualized and studied from LD plots generated for each chromosome using HAPLOVIEW software version 4.2. LD blocks were defined using confidence intervals (Gabriel et al. 2002).

KASP-SNP design and single marker analysis of BILs
SNP markers from sliding windows QBPH4.1 and QBPH4.2 and previously reported BPH resistance genes were designed using 100 bp flanking sequence on either side of the polymorphisms. Two allele-specific forward primers were designed with differences at the 3′ ends where the target SNP was located, and one common reverse primer was designed following KASP genotyping technology (Additional file 4: Table S7). Single marker analysis was performed by one-way analysis of variance (ANOVA) and simple linear regression using GenStat software (18th edition) (https://www.vsni.co.uk/software/ genstat/).

Reverse screening of JHN mutant population
Screening of 9323 JHN mutant lines was conducted using an SNP in OsLecRK3, 21-bp indel in OsSTPS2, and SNP in BPH32. DNA was extracted using the cetyltrimethylammonium bromide (CTAB) method. Approximately 2-5 ng of DNA was genotyped using the KASP genotyping platform (KBioscience/LGC, Middlesex, UK). The M 6 seeds of selected mutant lines were retrieved from the mutant seed bank, germinated, and tested for BPH resistance.
Gene sequences were assembled using CAP3 assembly (Huang and Madan 1999). Amino acid gene sequences were translated by using the ExPAsy Translate Tool (https://web.expasy.org/translate/). Nucleotide sequences of the genes from mutant and the predicted amino acid sequences of rice lines were compared using Clustal Omega (https://www.ebi.ac.uk/Tools/msa/clustalo/).

HP mining
Nucleotide sequences of known BPH resistance genes (BPH14, BPH29, BPH18, BPH26, and BPH9) from different BPH resistant rice cultivars were compared with the genomic sequence of rice cultivar Nipponbare to identify SNPs. These SNP markers were used together with SNPs identified in QBPH4.1 and QBPH4.2 to generate the HP of candidate genes. HP mining was performed on each gene individually (Toivonen et al. 2000). The impact of each HP on BPH resistance was determined using a single haplotype analysis by one-way ANOVA and simple linear regression using GenStat software (18th edition) (https://www.vsni.co.uk/software/genstat/). HPs that affected BPH damage AUC significantly (P < 0.05) were associated with BPH resistance. The DRRs of significant HPs were calculated by dividing the mean AUC of a significant HP by the mean AUC of the JHN WT HP.

Phylogenetic analysis
To determine whether selected mutant lines were contaminated by outcrossing, phylogenetic analysis was performed. The BPH resistant mutant lines, JHN4 and JHN09962, together with five BPH susceptible mutant lines and 15 WT-like mutant lines were genotyped by sequencing. Nine extreme BILs, RH, KDML105, Pokkali, and other germplasm were included in this analysis. Phylogenetic analysis was conducted in MEGA7 (Kumar et al. 2016). A total of 8928 RAD-derived SNPs were used as inputs in the UPGMA method at 1000 replicates bootstrap tests (Sneath and Sokal 1973). Genetic distances were computed using the p-distance method (Nei and Kumar 2000).

Graphical genotyping
To dissect the genomic region flanking QBPH4.1 region, 84 polymorphic RAD-derived SNPs located in intergenic regions from 1.45-16.45 Mb on chromosome 4 were analyzed. 32 and 34 Thirty-two intergenic SNPs were used for graphical genotyping analysis of selected BILs while thirty-four intergenic SNPs were used for selected mutants. The KASP-SNP gene-specific markers were used within the OsLecRK1-3region.

Additional files
Additional file 1: Table S1. Evaluations for brown planthopper (BPH) resistance of the progeny using 6 BPH populations; KPP = Kamphaeng Phet, NAN = Nan, PSL = Phitsanulok, UBN = Ubon Ratchathani, TPY = Ta Phraya, HTL = Huai Thalaeng. Table S2. Summary of ddRADseq data. Table S3. distribution on the rice chromosomes of identified SNPs in each ddRADseq library. Table S4. List of SNP identified on the BPH32gene-containing sliding window (0.03-2.03 Mb) on chromosome 6. SI = SNP index; QS = QBPHS; QR = QBPHR; IndS = BPHS; IndR = BPHR; AAC = amino acid change; (−) = synonymous variant; (+) = missense variant or stop codon. DNA markers used for molecular breeding programs to improve BPH resistance were highlighted in gray. Table S5. Genotyping data at BPH32-SNP of selected mutant lines by reverse screening and BPH damage AUC tests using 3 BPH populations. Table S6. KASP-SNP marker from QBPH4.1, QBPH4.2, and BPH-R genes. Table S8. HP analysis of significant genes in 4 resistant mutant lines and 5 susceptible mutant lines. Table S9. List of overlapping primers for sequence analysis of OsLecRK1-3 genes in mutant lines. (DOCX 61 kb) Additional file 2: Figure S1. The aggressiveness of three BPH populations used for BPH resistance validation of mutant lines. Figure S2. Average AUC of 105 BILs and their parental varieties, KDML105 and RH. The selected extreme lines for susceptible and resistance individuals were lebeled in yellow and red respectively. Figure S3. Aggressiveness of three BPH population used for BPH resistance validation of mutant lines. (DOCX 23 kb)