The plant material consisted of O. glaberrima acc. IRGC102206, PR121, BC₁F₁, BC₂F₁, and BC₂F₂ populations derived from the cross of PR121 and O. glaberrima acc. IRGC102206 O. glaberrima acc. IRGC102206 were identified as highly resistant and susceptible to M. graminicola, respectively, in a previous study (Kaur, 2020). The resistance was confirmed over 2 years of screening both in the nematode-infested pots and sick plot under field conditions, as the number of galls per plant was significantly lower as compared to susceptible checks. O. glaberrima has characteristics of early maturity, moderate to tall height, seed shattering, lower yield, and more resistance to various diseases and pests. However, PR121 has short stature, bacterial blight resistance, and better lodging tolerance.

The BC₁F₁, BC₂F₁, and BC₂F₂ populations were grown in the nursery during cropping seasons 2018–2020. The 25-day-old seedlings were transplanted in irrigated field conditions with plant-to-plant and row-to-row distances of 20 and 30 cm, respectively. Each plant of BC₁F₁, BC₂F₁, and BC₂F₂ after 15 days of transplanting was split into four plantlets, and three replicas of each BC₁F₁, BC₂F₁, and BC₂F₂ plant were transferred to a nematode-infested sick plot with an initial population density of 1 J2/g of soil while one replica was raised under non-infested conditions. Standard agronomic practices were implemented during the raising of the crop, except that the soil was not flooded during screening in the nematode-infested sick plot. Each plant of all populations was uprooted from the nematode-infested sick plot after 60 days of transplanting. The roots were washed immediately under running tap water to count the galls per root system. Root gall index was calculated on a scale of 1–5 as given by Gaur et al. (2001). Rating was done as follows: 1 for 0–1 gall (highly resistant), 2 for 1–10 galls (resistant), 3 for 11–30 galls (moderately resistant), 4 for 30–100 galls (susceptible), and 5 for >100 galls (highly susceptible). The segregation pattern for nematode resistance was checked in each generation using standard chi-square analysis for the goodness of fit. Data for each BC₁F₁ plant on different morphological characters like plant height (cm), root length (cm), fresh shoot weight (g), and fresh root weight (g) were measured immediately after the uprooting of the plants. The shoots and roots of each plant were packed separately in brown paper bags for drying, to achieve constant weight for measuring dry shoot and root weight (g). Data were analyzed using a generalized linear model (GLM) of SAS software version 9.4 (SAS Institute, Cary, NC, United States).

Genomic DNA was isolated from the young leaves of parents and each BC₁F₁ plant using the CTAB method (Doyle and Doyle, 1987). The degradation and contamination of DNA were checked on 0.8% agarose gel while DNA was quantified using Thermo scientific NanoDropTM 1000 spectrophotometer. A total of 512 simple sequence repeat (SSR) markers from the universal core genetic map (Orjuela et al., 2010) spanning all 12 rice chromosomes were used for the parental polymorphic survey on PR121 and O. glaberrima acc. IRGC102206. Primer sequences were retrieved from the Gramene database (http://www.gramene.org/; IRGSP, 2005). In vitro amplification using polymerase chain reaction (PCR) was performed in a 96-well PCR plate in Eppendorf and Applied Biosystems master cyclers. The total PCR reaction of 20 μl was prepared using the following components: 100 ng template DNA, 0.50 μM each of forward and reverse primers and 2× Emerald Amp^® GT PCR Master Mix containing an optimized buffer, PCR enzyme, dNTP mixture, gel-loading dye (green), and a density reagent. A negative control (without template DNA) was included in each amplification reaction. PCR profile of 95°C for 5 min, followed by 30 cycles of 1 min at 94°C, 1 min at 55–60°C, and 1 min at 72°C with a final extension of 10 min at 72°C was used for amplification. The amplified products were resolved in 3.0% agarose gel, and amplicons were scored by comparing them to parental alleles. The linkage map was constructed using the Kosambi mapping function of QTL IciMapping version 4.1 (Meng et al., 2015) through MAP functionality in a graphic form representing the position of markers within linkage groups by using a threshold LOD score of 3.0. Composite interval mapping (CIM) at a 95% threshold level was used for the identification of QTLs based on 1,000 permutation tests using Windows QTL cartographer version 2.5 (Wang et al., 2012). The position of putative QTLs corresponded to the location (in centiMorgans) of peak LOD scores in the scan of individual chromosomes and was designated according to the chromosome position. The proportion of observed phenotypic variance attributable to the QTL was estimated by the coefficient of determination (R ²) using the maximum likelihood of CIM.

The DNA concentration of each sample was normalized to 500 ng/μl, and 5 μg of the total DNA of the individual plant was used for making extreme bulks. An equal amount of DNA from 10 individual plants with few galls (1–2.33 galls) was mixed to generate the resistant DNA bulk. Similarly, an equal amount of DNA from 10 individual plants with higher gall numbers per plant (40–75 galls) was mixed to generate the susceptible DNA bulk (Figure 1). Paired-end sequencing libraries from the two extreme bulk and parents using 2 μg of DNA were prepared according to the Illumina manufacturer’s instructions. The sequencing of libraries using the Illumina HiSeqTM 2500 platform was outsourced from NGB Diagnostics Pvt Ltd., India. Illumina sequencing of genomic libraries for each of the parents and the two bulks (2 × 150 bp) produced 33–35 million reads per sample for a total of 136 million raw reads. The quality of raw reads was assessed using FASTQC (version 0.11.8; Andrews, 2010) with default parameters. Poor-quality sequences were filtered and removed while contaminated adapter sequences and any unwanted bias from their ends were trimmed using Trimmomatic (version 0.39; Bolger et al., 2014). A Phred score of 30 was kept as the overall quality threshold for raw reads. The filtered reads were further re-checked for quality using FASTQC.

High-quality sequences were aligned and mapped to the Oryza sativa Indica Group ASM465v1 reference sequence of cultivar 93-11, available at Ensembl plants (https://plants.ensembl.org/Oryza_indica), using Bowtie 2 algorithm with default parameters (version v2.0.0; Langmead and Salzberg, 2012). The Bowtie 2 default mode is faster than all Burrows-Wheeler Aligner (BWA) modes and more than 2.5 times faster than the BWA default mode. All Bowtie 2 modes aligned a greater number of reads than either BWA or short oligonucleotide alignment program 2 (SOAP2). To keep only uniquely mapping reads, the output SAM files were converted into BAM files, then read groups were added, sorted, and indexed using SAMtools (version 0.1.19; Li et al., 2009). The output BAM files containing uniquely mapped reads were used for SNP Calling through the GATK (Genome Analysis Toolkit) pipeline. Subsequently, GATK’s HaplotypCaller component (version 4.0; Mckenna et al., 2010) was used to perform a joint variant calling of all samples. The indels and missing data were filtered out using the variant filtration parameter of vcftools (https://github.com/vcftools/vcftools). The filtered VCF file in table format was used as an input file for QTLseqr (https://github.com/bmansfeld/QTLseqr) package developed by Mansfeld and Grumet (2018). SNPs with a reference allele frequency of 0.2 from both the bulks were filtered out as these might be due to sequencing or alignment error. In G′ approach, run GprimeAnalysis first counted the number of SNPs within the sliding window, and then a tricube-smoothed ∆SNP index was calculated within a window size of 1.0-Mb genomic region. The ∆SNP index (>0.1) was used to calculate p-values (<0.05) and an FDR(q) of 0.01 to identify potential QTLs associated with root-knot nematode resistance. The G′ determines the statistical significance of QTLs as background noise is less and also addresses the linkage disequilibrium (LD) between SNPs. One important advantage of this method is that p-values can be estimated for each SNP using non-parametric estimation of the null distribution of G′ (Magwene et al., 2011). Comparison of the QTL-seq method (Delta-SNP index) and G′ method by Mansfeld and Grumet (2018) showed that a confidence interval of 99% with the QTL-seq method was not as stringent as using an FDR of 0.01 in the G′ method. All commands and codes for BSA-seq analysis are available at https://github.com/bmansfeld/QTLseqr.

The QTL regions harboring the candidate genes based on the annotation of Oryza sativa indica Group ASM465v1 (https://plants.ensembl.org/Oryza_indica) were identified. To identify non-synonymous SNPs among the candidate genes of the two parents (PR121 and O. glaberrima), nucleotide changes were investigated using ExPASy translate tool (http://web.expasy.org/translate/). Potential candidate genes and their corresponding Ensembl IDs were further subjected to the ShinyGO v0.74 database (Ge et al., 2019) to obtain gene ontology (GO) annotation against O. sativa subsp. indica. GO enrichment was calculated by a p-value cut-off (FDR) at 0.05 for the genes.

The F₁s (with medium-sized ligule) generated from the cross of PR121 (S) and O. glaberrima acc. IRGC102206 (R) were partial to completely sterile; therefore these F₁s were backcrossed with PR121 to develop the BC₁F₁ population. A total of 10,800 spikelets of F₁ plants were cross-pollinated and a 0.95% seed setting of BC₁F₁ (103 seeds) was obtained. Out of the 103 seeds, 69 seeds (67%) were germinated in the nursery and transplanted in controlled (normal) conditions. After 25 days of transplanting, three replicas of each BC₁F₁ plant along with their parental genotypes were screened in the nematode-infested sick plot. Based on gall number, 39 and 30 BC₁F₁ plants were categorized as resistant and susceptible, respectively, for root-knot nematode resistance that corresponded to a single locus segregation ratio statistically (Table 1). The resistant parent IRGC102206 showed gall number from 0 to 1 while the susceptible parent PR121 exhibited gall number in the range of 64–82 (Figure 2). The average gall number and gall index among BC₁F₁ individuals ranged from 1.00 to 70.33 and 1.33 to 4.00, respectively (Table 2). The resistant parent O. glaberrima accession IRGC102206 had an overall mean gall index of 1.0 while the susceptible parent PR121 had a mean gall index of 4.0. The segregation pattern in successive generations, that is, BC₂F₁ and BC₂F₂, was authenticated, thereby confirming that nematode resistance is governed by a single dominant gene (Table 1). Overall, BC₁F₁ plants showed stunted growth under nematode infestation. Plant height and root length among BC₁F₁ plants ranged from 61 to 118 cm and 7.57 to 21.03 cm, respectively. Significant differences were observed among BC₁F₁ plants for all the traits (Table 2). There was a nominal decrease in growth parameters in the resistant individuals compared to the susceptible plants indicating that the plant growth parameters were affected by nematode infestation.

We also carried out rough QTL mapping using 69 individuals of the BC₁F₁ population. A total of 512 microsatellite markers spanning 12 chromosomes of rice were used for the parental polymorphic survey and the markers per chromosome varied from 53 (chromosomes 1 and 2) to 31 (chromosome 9). Parental polymorphism among all chromosomes ranged from 7.6% to 50.9% with an average of 31.25%. A comprehensive list of polymorphic SSR markers is given in Supplementary Table S1. A total of 100 polymorphic SSR markers were genotyped and a genetic linkage map was generated with a total map length of 1,901.21 cM, with an average distance of 22.63 cM. One putative QTL associated with gall numbers was detected on chromosome 6 designated as qGN6.1, between the marker interval of RM3183 and RM27001 explaining the phenotypic variance of 41% at an LOD score of 3.95 (Table 3). Two QTLs for dry root weight and dry shoot weight (qDRW3.1 and qDRW3.2) were co-localized on chromosome 3 whereas QTL for fresh root weight, qFRW8.1, was located on chromosome 8 with 17% of total phenotypic variance (Table 3).

Based on the gall number of BC₁F₁ population, 10 plants from each of the extreme values of frequency distribution were selected and pooled as resistant bulk (RB) and susceptible bulk (SB), respectively. The whole-genome re-sequencing data from PR121, IRGC102206, resistant, and susceptible bulks were aligned with the O. sativa cultivar 93-11 reference genome using Bowtie 2 algorithm with default parameters. A total of 35.0, 33.6, 34.3, and 34.0 million paired-end reads were generated from IRGC102206, PR121, RB, and SB, respectively (Table 4), and 32–34 million reads of paired-end sequences were retained which were equivalent to 12.8× to 14.9× coverage of the rice genome indicating the high quality of the sequencing data (Supplementary Figure S1). Approximately, 31.23 (96.11%), 29.78 (87.66%), 31.07 (93.33%), and 31.05 (93.86%) million read pairs of PR121, IRGC102206, RB, and SB were uniquely mapped to the reference genome, respectively (Table 4). The GC content of raw reads ranged from 42% to 43% for all the samples.

A total of 3,692,066 variants (SNPs/indels) were identified by a joint variant calling with the reference genome through the GATK pipeline. RB had the highest variant (1,641,906) followed by IRGC102206 (1,049,225), PR121 (695,856), and SB (305,079). Maximum SNPs (284,012) were specific to chromosome 1 while chromosome 11 had the least SNPs (123,492). Likewise, the variant rate (SNPs/Mb) also varied among chromosomes, and chromosome 6 possessed the highest variant rate (1 SNP/195 bases) while the lowest variant rate (1 SNP/136 bases) was detected for chromosome 7 (Supplementary Table S2). A total of 1,416,115 SNPs were detected between RB and SB. Based on reference allele frequency (0.20) and maximum total depth (300), 174,651 SNPs were used for QTL identification through BSA-seq analysis.

The BSA-seq analysis detected nine putative QTL regions on chromosomes 2, 3, 4, 6, 11, and 12, based on the calculation of Gʹ values of SNP within a window size of 1.0-Mb genomic region across the entire length of chromosomes. QTLs on chromosomes 3, 6, and 12 with large G′ peaks above FDR(q) of 0.01 were considered as major QTL regions responsible for root-knot nematode resistance in rice. Despite the major G′ peaks, five minor G′ peaks at the margin of significance have been identified on chromosomes 2, 3 (left to qNR3.1), 4, and 12 (on both sides of qNR12.1). Two adjacent QTLs on chromosome 3, namely, qNR3.1 and qNR3.2, had peaked in genomic intervals of 22.00–23.3 Mb and 24.0–25.6 Mb, respectively (Figure 3B). Similarly, QTLs on chromosome 12, designated as qNR12.1 and qNR12.2, were located at 5 Mb distance between intervals of 10.5–12.7 Mb and 17.5–22.9 Mb, respectively. qNR12.2 was present at the proximal region on the long arm of chromosome 12 displaying two adjacent component peaks (Figure 3F). The region covered by significant QTLs varied from 0.2 Mb (qNR4.1) to 5.4 Mb (qNR2.2). The direction of the ΔSNP index value and the allele frequency difference (AFD) indicated that qNR2.1, qNR4.1, qNR6.1, and qNR6.2 originated from the donor parent IRGC102206 (Table 5). However, QTLs on chromosomes 3, 11, and 12 originated from a susceptible parent, PR121. The number of SNPs present in qNR3.1, qNR3.2, qNR11.1, qNR12.1, and qNR12.2 region was 192, 330, 328, 691, and 322, respectively.

The genomic interval between qNR6.1 and qNR6.2 was 2.4 Mb (11.2–13.6 Mb) and 0.6 Mb (15.4–16.0 Mb), respectively. The main Gʹ peak exhibited a subpeak in the QTL, qNR6.1, region whereas qNR6.2 exhibited a clear sharp peak region and fall near the centromere region (ranging approximately 17.00–22.59 Mb) of the short arm of chromosome 6. These two QTLs must be in continuity suggesting nematode resistance locus spanned within 3.0-Mb regions indicating a significant association for nematode resistance and hence considered as the major locus. The peaks of qNR2.1 and qNR4.1 were at the margin of significance, that is, qNR2.1 peak slightly exceeding and qNR4.1 nearly reaching the threshold; and, therefore, referred to as minor QTLs. The number of SNPs present in qNR2.1, qNR4.1, qNR6.1, and qNR6.2 regions was 345, 103, 264, and 27, respectively. Among all QTLs, qNR6.1 and qNR12.2 showed the highest Gʹ peaks that greatly exceeded the threshold level indicating the candidate regions for root-knot nematode resistance and susceptibility, respectively. Also, the genomic interval for root-knot nematode resistance on chromosome 6 overlapped with the identified gall number QTL (qGN6.1) through SSR markers. Therefore, qNR6.1 and qNR6.2 were selected as promising QTL regions originating from the donor parent for the identification of candidate genes and their sequence variations allied with root-knot nematode resistance.

The candidate genes within the QTL regions were identified based on BGI annotation. The genomic interval of candidate genes and their annotation have been summarized in Supplementary Table S3. The total number of candidate genes annotated in QTLs, qNR3.1, qNR3.2, qNR11.1, qNR12.1, and qNR12.2, was 7, 10, 16, 10, and 66, respectively. A total of nine annotated genes were identified within a 2.4-Mb interval of qNR6.1 and three genes within 0.6-Mb interval of qNR6.2. The less number of annotated genes within the QTL (qNR6.1 and qNR6.2) interval is because of low recombination rate in the centromere region. In the qNR2.1 region, seven genes were annotated whereas none of the candidate genes was annotated in the 240-Kb region of qNR4.1. We emphasized on those QTLs which originated from the donor parent to identify non-synonymous SNPs in the candidate genes. Among the nine annotated genes of qNR6.1, eight genes contained a total of 55 non-synonymous SNPs. In the candidate region of qNR6.2, a total of 10 non-synonymous SNPs were detected (Figure 4). The candidate gene encoding glycosyltransferase (BGIOSGA022720) had a maximum number of non-synonymous substitutions by SNPs (15), while peroxidase (BGIOSGA022766) and thioredoxin reductase (BGIOSGA021328) had the least number of non-synonymous SNPs (1). The examination of sequence reads of peroxidase and thioredoxin reductase on Integrative Genome Viewer (IGV) showed that non-synonymous SNPs were specific to the resistant parent; however, the susceptible parent also shared identical alleles to the reference genome (Figure 5). The majority of SNPs have been identified in the coding regions of candidate genes; however, mitogen-activated protein kinase (BGIOSGA022897), fructose-6-P-1-phosphotransferase subunit alpha (BGIOSGA022810), peroxidase (BGIOSGA022766), and thioredoxin reductase (BGIOSGA021328) had intronic SNPs (Table 6). Comparison of amino acid sequences of glycosyltransferase (BGIOSGA021421) and auxin-responsive protein (BGIOSGA022887) revealed that none of the SNPs caused non-synonymous substitutions.

In the qNR2.1 region, a total of 85 SNPs were identified within seven candidate genes and only 16 of them had non-synonymous substitution. Among these genes, the miR164 gene (ENSRNA049493468) had no SNP, and none of the SNPs in growth-regulating factor 1 (BGIOSGA005497) and hexosyl transferase (BGIOSGA005471) caused non-synonymous substitutions. Pectin esterase (BGIOSGA005479) genes had 36 synonymous SNPs and 8 non-synonymous SNPs. The maximum non-synonymous SNPs were present in pectinesterase (BGIOSGA005479) followed by a UMP-CMP kinase gene (BGIOSGA009139) among the candidate genes of qNR2.1 (Table 6). Therefore, candidate genes with non-synonymous SNPs could be the preferred genes for further studies. GO annotation showed the involvement of these genes in various biological processes and molecular functions. The candidate genes were predicted to be involved in different types of biosynthetic, metabolic, and catabolic processes (Figure 6). Transferase activity, protein metabolic process, and catabolic process appeared as the top enriched GO terms.

M. graminicola is emerging as a devastating pest and affects both upland and irrigated rice causing up to 80% yield losses (Mantelin et al., 2017). It has been observed that intensely galled roots were shorter than the non-infested plants because terminal galls inhibited the roots from further elongation and prevented the plant’s root system to absorb and translocate water (Kaur, 2020). Therefore, the severity of galling has a direct consequence on plant yield. To alleviate the plant damage and yield losses, breeding for rice root-knot nematode resistance is one of the efficient and most economical strategies. With the advent of molecular markers, several nematode resistance genes and QTLs have been identified in different crop species such as tomato, potato, cotton, sugarbeet, and cowpea (Cai et al., 1997; Vander-Vossen et al., 2000; Wang et al., 2006; Wu et al., 2009; Ndeve et al., 2018).

In the past, most studies reported that rice root-knot nematode resistance was a quantitative trait (Amoussou et al., 2004; Prasad et al., 2006; Shrestha et al., 2007; Jena et al., 2013; Galeng-Lawilao et al., 2018). Here, we report that nematode resistance in O. glaberrima acc. IRGC102206 is under dominant monogenic control. Ideally, in BC₁F₁, 1:1 segregation is expected for heterozygous alleles and homozygous alleles of the recurrent parent. Due to the presence of heterozygous alleles, the trait under consideration should be moderately resistant. However, we got plants of varying resistances, from resistant to moderate resistant, suggesting that there might be non-allelic minor modifying genes that might act in an additive manner such that the phenotype is a resultant of the copy number of the alleles present in the genotype toward resistance. Similarly, the BC₂F₂ progenies showed wide variation for gall numbers indicating that there are a few modifying minor genes which affect the expression of a major gene to alter the resultant disease reaction. Previously, O. glaberrima variety CG14 has been identified as resistant to M. graminicola (Plowright et al., 1999; Soriano et al., 1999). But O. glaberrima has not been explored further to study the genetics of root-knot nematode resistance due to hybrid sterility in the interspecific crosses. However, introgression of O. glaberrima genes is possible by repeated backcrossing and doubled haploid breeding, even though there is a risk of losing the desirable traits from the parents (Jones et al., 1997). However, Neelam et al. (2020) have successfully transferred the bacterial blight resistance gene xa45(t) into Pusa 44 background and mapped the gene on chromosome 8 from O. glaberrima acc. IRGC102600 by backcross breeding.

So far, no information is available for rice root-knot nematode resistance mapping in African rice cultivars. To the best of our knowledge, this is the first report on the identification of a major locus and a few minor QTLs for root-knot nematode resistance using interspecific BC₁F₁ population derived from a cross of O. sativa cv. PR121 and O. glaberrima acc. IRGC102206 by traditional QTL mapping and BSA-seq approach. Traditional QTL mapping identified one major large effect QTL for gall number, qGN6.1, explaining 41% of the phenotypic variance localized in the 12.29–16.36 Mb region with increasing alleles from the donor parent. This finding was consistent with the results of BSA-seq as the major locus for nematode resistance overlapped the genomic region between 11 and 16 Mb of chromosome 6, thereby validating the existence of a major QTL for root-knot nematode resistance in this region. We also detected QTLs on chromosomes 2 and 4 with minor effects using the BSA-seq approach. Based on the BC₁F₁ population, we anticipated that the SNP index of the RB had a mixture of PR121 and IRGC102206 alleles, and SB had PR121 alleles for the genes/QTLs associated with nematode resistance. The SNP index of the QTL regions on chromosomes 2, 4, and 6 agreed with this expectation, but the SNP index of the QTLs on chromosomes 1, 3, 11, and 12 displayed a contrasting blueprint. Thus, the present results inferred that QTLs on chromosomes 2, 4, and 6 probably conferred the nematode resistance, whereas QTLs on chromosomes 1, 3, 11, and 12 likely contributed to the susceptibility. Identification of the higher number of genomic regions contributed by PR121 is likely because we used a backcross population having a higher proportion of recurrent parent (PR121).

Previous studies identified QTLs for root galling, eggs per total root system and eggs per gram of roots, nematode reproduction, and nematode tolerance on chromosomes 1, 3, 4, 5, 7, 8, 9, 11, and 12 (Jena et al., 2013; Dimkpa et al., 2016; Galeng-Lawilao et al., 2018). To date, there is only one report inferring that a single dominant resistance gene (Mg1(t)) located on chromosome 10 in Asian rice cultivar Abhishek confers resistance against M. graminicola (Mhatre et al., 2017). So far, the Hsa-1Og nematode resistance gene has been identified from O. glaberrima TOG5681 on the long arm of chromosome 11 between markers RM206 and RM254 that confers resistance to the cyst nematode (Lorieux et al., 2003). Lahari et al. (2019) identified genomic region from 23 Mb to the bottom of rice chromosome 11 for root-knot nematode resistance through the QTL-seq approach, which might be co-localized with the results reported by Lorieux et al. (2003). Similarly, Dimkpa et al. (2016) also identified an SNP (id11008353) associated with gall number on the long arm of chromosome 11. Also, QTLs, qNR12.1 and qNR12.2, on chromosome 12 have different genomic locations for gall number–associated SNPs (Dimkpa et al., 2016). However, the BSA-seq analysis in our study identified QTL for root-knot nematode resistance on the short arm of chromosome 11 (1.01–2.38 Mb) contributed by the susceptible parent. A root-knot nematode resistance locus from 23 Mbp to the bottom of rice chromosome 11 was identified by Dimkpa et al. (2016). This might be because different sources of root-knot nematode resistance were used to identify QTL regions, as the nature of gene actions varies with the background.

Shrestha et al. (2007) identified a QTL for the gall number on the long arm of chromosome 6 flanked by R2654 and RG778 markers, explaining 9.6% of the total phenotypic variation. The results indicated that QTL, qGN6.1/qNR6.1-6.2, and QTL reported by Shrestha et al. (2007) were present on either sides of the centromere of chromosome 6, suggesting that gall number QTL identified in our population was controlled by a different genetic locus. QTLs/linked markers for root-knot nematode resistance identified in previous studies have been summarized in Supplementary Table S4. In the present study, one QTL for fresh root weight and two QTLs for dry shoot weight and dry root weight were identified on chromosomes 8 and 3 (two QTLs), respectively. The QTLs for dry shoot weight and dry root weight were co-localized. Earlier, the QTL for fresh root weight was mapped on chromosomes 2 and 12 (Galeng-Lawilao et al., 2018). It has to be noteworthy that in the present study, genomic regions identified for different traits through either traditional mapping or BSA-seq are unique. The major QTL could be fine mapped to identify the putative candidate gene for nematode resistance.

Among the 19 candidate genes, 13 genes harbored non-synonymous mutations. Three glycosyltransferase genes were present in the qNR6.1 region, two of which (BGIOSGA022720 and BGIOSGA022727) had a total of 18 non-synonymous substitutions. Glycosyltransferases are involved in carbohydrate biosynthesis and cell-wall synthesis (Egelund et al., 2004). The upregulation of glycosyltransferase genes provided resistance against Pseudomonas syringae and M. incognita and M. hapla in Arabidopsis and tomato, respectively (Langlois-Meurinne et al., 2005; Schaff et al., 2007). However, its role in M. graminicola resistance is not understood currently in rice. In the present study, the peroxidase (BGIOSGA022766) gene located in qNR6.1 plays an important role in ROS detoxification by regulating the H₂O₂ level and the oxidation of toxic reductants during pathogen attack. A higher level of peroxidase activity was reported in resistant genotypes of tomato during M. incognita infection (Chawla et al., 2013) and sweet potato (Yeon-Woo et al., 2019). Transcriptome profiling of M. graminicola infested rice roots showed the transcript abundance of peroxidases compared to uninfected roots (Kyndt et al., 2012). Rice root-knot nematode interactions revealed higher transcripts for nucleotide binding, catalytic, phosphatase, hydrolase, and ATPase activity after 3 and 7 days post inoculation (dpi) of M. graminicola (Haegeman et al., 2013).

Thioredoxin reductase has a role in ROS signaling and the protection of antioxidant enzymes to establish plant immunity during pathogen attacks. The ROS accumulation during pathogen attack results in oxidative modification of reactive free thiols in S-sulfenic acids (–SOH) of signaling proteins and antioxidant enzymes such as catalase (CAT). Thioredoxin reductase reduces these thiol modifications to enable ROS signaling and protection of antioxidant enzyme activities (Mata-Perez and Spoel, 2019). Similar to thioredoxin reductase, two genes of the qNR6.2 region, acyl-coenzyme A oxidase (BGIOSGA021280) and mitogen-activated protein kinase (BGIOSGA022897), might have a potential role in signaling pathways. Acyl-CoA oxidase metabolizes the pheromone (ascaroside) secreted by plant-parasitic nematode into chemical signals through the peroxisomal β-oxidation pathway that might act as a repellent for nematodes and thus reduces infection (Manohar et al., 2020). Mitogen-activated protein kinases (MAPKs) are intracellular signaling molecules that induce a defense response against root-knot nematode. The differential expression of the MAPK genes occurred after the activation of the effector-triggered immunity (ETI) pathway in rice upon M. graminicola infection as reported by Hatzade et al. (2020). It has been postulated that the genes involved in defense responses, phenylpropanoid, and hormone pathways were induced in response to M. graminicola infection in O. glaberrima line TOG5681, compared to Nipponbare (Petitot et al., 2017).

In the QTL region qNR2.1, the kinesin-like protein gene identified has diverse roles in biotic and abiotic stresses ranging from environmental to developmental processes like cell division, cell expansion, tropisms, and hormonal signaling. Kinesin modulates the cell wall structure and function by affecting the orientation and structure of cellulose microfibrils within the cell wall (Abdelkhalek et al., 2019). Therefore, kinesin-like protein and glycosyltransferase might be involved in the synthesis and maintenance of cell wall integrity during plant–nematode interactions. Growth-regulating factors are transcription factors regulated post-transcriptionally by miRNA396 (miR396) in different plant species (Omidbakhshfard et al., 2015). The growth-regulating factors, GRF1 and GRF3, in Arabidopsis target a plethora of genes having roles in programmed cell death, hormone signaling, and basal defense responses during biotic stresses (Piya et al., 2020). GRF1and GRF3 control the development and differentiation of syncytia during nematode (Heterodera schachtii) infection in Arabidopsis. In the present study, we also identified GRF1 (BGIOSGA009139) in the qNR6.2 region. The presence of non-synonymous SNPs in these candidate genes suggests their potential role for root-knot nematode resistance in rice.