A genome-wide association study in Indian wild rice accessions for resistance to the root-knot nematode Meloidogyne graminicola

Rice root-knot nematode (RRKN), Meloidogyne graminicola is one of the major biotic constraints in rice-growing countries of Southeast Asia. Host plant resistance is an environmentally-friendly and cost-effective mean to mitigate RRKN damage to rice. Considering the limited availability of genetic resources in the Asian rice (Oryza sativa) cultivars, exploration of novel sources and genetic basis of RRKN resistance is necessary. We screened 272 diverse wild rice accessions (O. nivara, O. rufipogon, O. sativa f. spontanea) to identify genotypes resistant to RRKN. We dissected the genetic basis of RRKN resistance using a genome-wide association study with SNPs (single nucleotide polymorphism) genotyped by 50K “OsSNPnks” genic Affymetrix chip. Population structure analysis revealed that these accessions were stratified into three major sub-populations. Overall, 40 resistant accessions (nematode gall number and multiplication factor/MF < 2) were identified, with 17 novel SNPs being significantly associated with phenotypic traits such as number of galls, egg masses, eggs/egg mass and MF per plant. SNPs were localized to the quantitative trait loci (QTL) on chromosome 1, 2, 3, 4, 6, 10 and 11 harboring the candidate genes including NBS-LRR, Cf2/Cf5 resistance protein, MYB, bZIP, ARF, SCARECROW and WRKY transcription factors. Expression of these identified genes was significantly (P < 0.01) upregulated in RRKN-infected plants compared to mock-inoculated plants at 7 days after inoculation. The identified SNPs enrich the repository of candidate genes for future marker-assisted breeding program to alleviate the damage of RRKN in rice.

feeding site (3-5 giant cells) that serves as the permanent nutrient sink for the nematode. Typical hook-shaped galls are formed around the giant cells that affect water and nutrient translocation of plant and in turn plant growth and several yield-contributing traits (such as root and shoot weight/length, tiller number, panicle/spikelet length, grain weight and percentage filled grains) are impaired [3][4][5][6]8].
In order to mitigate the loss by RRKN in rice, use of chemical nematicides, crop rotation (with mung bean, mustard and sesame as poor hosts) and continuous flooding are the recommended practices [3,6]. Although prolonged flooding may reduce the RRKN populations by inhibiting the second-stage juveniles (J2) to invade rice roots, increasing scarcity of water for agricultural use may render this practice impractical to adopt. Rotation with other crops during the growing season may incur additional costs to bear for the small-scale rice farmers who constitute a major part of the farming community in Asia. Considering the continual withdrawal of nematicides from the market due to its residual effect on the environment, growing resistant/tolerant rice cultivars offer an economically and environmentally viable option to manage RRKN problems in rice.
Resistance to RRKN has been unraveled in Oryza longistaminata, O. glaberrima (African rice) and even in O. sativa (Asian rice). Nevertheless, very few of these are truly resistant and majority of Asian rice genotypes are susceptible to M. graminicola [8,9]. Efforts were made to introgress RRKN resistance from O. glaberrima to O. sativa, however, interspecific progenies did not exhibit the same degree of resistance documented with O. glaberrima [10,11]. Resistance to RRKN in rice is reportedly quantitative in nature and governed by several genes with additive effect. Five QTLs (quantitative trait loci, in chromosomes 1, 2, 6, 7 and 9) were found to be associated with root galling by RRKN in RILs (recombinant inbred lines) derived from the cross of O. sativa accessions 'Bala' and 'Azusena' [12]. Jena et al. [13] reported QTLs associated with number of galls and eggs per root system on  [7,8].
Genome-wide association study (GWAS) offers a viable strategy to complement traditional bi-parental linkage mapping in order to discern the genetic basis of trait variation [14]. Compared to classical biparental QTL mapping, GWAS takes less time to dissect quantitative traits at higher resolution because of the greater recombination rate and distribution of SNP (single nucleotide polymorphism) markers across the genome of natural plant populations. Additionally, prior knowledge of genotype pedigree/crosses is not required in GWAS as bi-allelic SNP markers can be used to estimate the population structure [15][16][17]. To date, only in a few studies, GWAS has been used to identify novel QTLs for phytonematode resistance or susceptibility in different plants including Arabidopsis [18,19], wheat [20], soybean [21] and rice [22]. QTLs associated with root galling were mapped to the chromosomes 1, 3, 4, 5, 11 and 12 in a global panel of Asian rice [22]. Nevertheless, additional studies on geographically different rice panels (such as wild rice accessions that harbor greater genetic diversity) will enrich the resistant gene pool which can be deployed in future rice breeding programs. Due to ever increasing population pressure and rapid urbanization wild rice resources are depleting at an alarming rate. The Indo-Burma region is considered the biodiversity hotspot of wild rice populations [23,24].
In the present study, we investigated the genetic basis of RRKN Figures 1a and 1b show the frequency distribution of the gall numbers and nematode multiplication factor (MF), respectively, while all values (including numbers of total endoparasites, egg masses and eggs per egg mass) are detailed in Supplementary Table S1. According to Fig. 1a, 49 of the accessions showed highly resistant response (0-2 mean number of galls), 8 exhibited highly susceptible response (7-12 galls), while 87 accessions showed moderately resistant response (4-5 galls). The mean number of galls was recorded to be 7-8 in reference genotypes Taipei 309 and Pusa 1121. A similar trend was observed in terms of MF; a moderate MF value of 12-20 was calculated in maximum of 89 accessions (Fig. 1b).
The results of a repeat run with 50 randomly selected accessions revealed a strong correlation with the initial gall (r = 0.780, P < 0.001) and MF (r = 0.720, P < 0.001) data, indicating the high reproducibility of PF-127-based assay.
When all the genotypes were categorized into different taxonomic groups of wild rice, accessions from O. nivara (2.9 mean numbers of galls) showed the least number of galls as compared to the accessions that are classified as O. rufipogon (4.1) and O. spontanea (5.2). A significant difference (P < 0.01, n > 30) in gall numbers between different species (Fig. 1c) was found, which explains the genetic variation between wild rice species in our study. When MF values were taken into account a significant difference (P < 0.01, n > 30) was observed between O. spontanea and O. nivara/O. rufipogon but not between O. nivara and O. rufipogon (Fig. 1d). When accessions were stratified into different agro-climatic zones (independent of species), no significant difference was observed. This may be because each geographic location contained all the three species of wild rice (Supplementary  Island Regency (IR). Cluster V nested O. nivara accessions from MGP (Fig. 4c). Our results suggested that population structure analysis may not always conform with diversity analysis as in former case inherent genetic structure are inferred based on their ancestral lineage while in latter genotypes can be grouped into additional clusters according to their geographical origin.

Genome-wide association analysis (GWAS)
A total of 17 significant SNPs associated with M. graminicola resistance (based on four different parameters including gall number, egg mass, eggs per egg mass and MF) to rice were identified by GWAS using mixed linear model (MLM), which controlled for population structure and familial relatedness (Fig. 5a, 5b; Table 1). The 17 association SNPs explained 5.33-12.61% of the total phenotypic variation (Table 1). We identified two SNPs to be associated with gall number, appearing to represent two QTLs on chromosome 2 and 6 in rice (threshold for significance -log 10 (P) > 4; Fig. 5; Table 1). Two SNPs each were associated with egg mass (on chromosome 2 and 4) and eggs per egg mass (in chromosome 1 and 11). No significant SNPs for total endoparasites were observed (only one SNP was detected at -log 10 (P) > 2.75). The greatest number of SNPs (9) was detected (in chromosome 1, 2, 3, 4, 6, 10 and 11) to be associated with the most important trait, i.e. nematode multiplication factor (-log 10 (P) > 4; Fig. 5; Table 1). Overall, two QTLs on chromosome 1 and 4 were found to be associated with all the traits (Table 1). Percent heritability explained by corresponding SNPs for phenotypic traits such as gall number, egg mass, eggs per egg mass and MF were 78, 71, 65 and 89%, respectively ( Table 2). were located on different chromosomes (Table 1). Several resistance proteins such as leucine-rich repeat (LRR), NBS-LRR, Cf2/Cf5, RGA3, Lr34 analog ABC transporter, zinc finger motifs, WD repeat, Fbox (analogous to WD and LRR) were located in different QTLs (Table 1). Additionally, a number of protein kinases, Exo70 exocyst complex and kelch motifs were identified which are not listed in Table 1. Gibberellin, pathogenesis-related (PR) proteins, HEAT repeats and 14-3-3 proteins which are required for plant immunity to pathogen stress were also identified. Aquaporin, nodulin and auxin regulatory proteins were identified which are known to regulate the giant cell initiation and maintenance during nematode infection (Table 1).

Discussion
Broad and novel genetic resources for M. graminicola resistance in rice The wild relatives of rice harbor rich and novel genetic resources which can be used to improve pest and disease resistance in cultivated rice [30,31]. For example, resistance to grassy stunt virus, bacterial leaf blight, neck blast and brown plant hopper was successfully introgressed into the cultivated rice from their wild relatives [32][33][34]. A number of beneficial traits including tolerance to biotic and abiotic stresses have been lost in cultivated rice which possesses a narrow genetic base because of domestication and breeding bottlenecks [23,24]. According to an estimate, modern rice varieties have retained only 20% of the genetic diversity present in their wild relatives [35]. India has unprecedented diversity of wild rice germplasm and landraces that are spread over fifteen diverse agro-climatic zones [23,24]. Therefore, these untapped genetic resources were taken as the prime candidates for the association panel in our study in order to unravel the conserved loci that govern M. graminicola parasitism. The higher susceptibility of indica and especially japonica rice is not surprising maybe because no selection is in function to introduce resistance in japonica as RRKN is not a pest of temperate rice cultivation [22]. Intriguingly, a majority of O. nivara accessions from middle Gangetic plains (a suspected genetic diversity hotspot; [23]) has shown the least susceptibility to RRKN infection. This may be because of their eco-geographical isolation in accordance with the rationale that gene flow is inversely related to geographic distance of natural populations [45].
Genetic loci determining the parasitic success of M. graminicola in rice By analyzing the marker-trait associations in 272 rice accessions via GWAS, we identified 17 significant SNPs that governed the RRKN resistance-related traits such as numbers of galls, egg mass, eggs per egg mass and MF ratio. By applying the -log 10 (P) score of 4 as a threshold for significance, we identified two QTLs at chromosome 1 and 4 of rice were associated with all the traits. A similar stringency level or a less stringent threshold (-log 10 (P) > 3) was adopted while reporting significant SNPs associated with nematode resistance/susceptibility in rice [22], wheat [20], soybean [21] and Arabidopsis [46]. In particular, four SNPs in chromosome 4, three SNPs each in chromosome 1, 2, 11, two SNPs in chromosome 6, and one SNP each in chromosome 3 and 10 were found to be associated with different traits. The trait such as total endoparasite counts in rice accessions was not taken into consideration as only one SNP was detected for this trait at -log 10 (P) > 2.75. Lowering the threshold for significance in GWAS study may increase the false discovery rate which in turn reveals more common alleles with smaller effect size in test populations [19].
A number of transcription factors (bZIP, SCARECROW, MYB, MADS-box, WRKY, ARF, GRAS etc.) that regulate plant defense responses to nematode and various biotic stress [47][48][49][50][51][52][53][54], were located on different chromosomes within 200-kb LD of the suspected QTLs in this study. The 200 kb window was adopted because of the fact that LD decay occurs in rice at 50-500 kb [25][26][27][28]. Although in wild rice the decay is much higher at 1-200 kb, O. nivara maintains LD over larger distance due to their high level of self-pollination [29].  [65]. The nucleotide-binding state of NBS-LRR proteins regulates the activity of plant resistance (R) proteins that are involved in activation of plant innate immunity upon pathogen recognition [66]. The R gene in tomato contains Cf2/Cf5 locus (encode LRR) that confer resistance in tomato to M. incognita [55]. Upon pathogen invasion in host plants, the activated membrane-bound NBS-LRR immune receptors translocate to the cell nucleus and interact with specific transcription factors (although not exclusive, comprise members of ERF, bHLH, bZIP, MYB, NAC and WRKY families) which modulate the plant immunity [50]. In our qRT-PCR analysis, compared to uninfected plants the expression of NBS-LRR and Cf2/Cf5 protein was positively regulated in M. graminicola-infected Pusa 1121 at 7 dpi. In coherence, other defense response regulatory genes such as MYB TF, SCARECROW, zinc fingers, bZIP TF and 14-3-3 were upregulated in nematode-infected plants than the uninfected ones.
Overexpression of a rice ADP-ribosylation factor induced pathogen resistance in tobacco by regulating the transcript accumulation of pathogenesis-related (PR) genes and salicylic acid (SA) [61]. OsRac1 (GTP binding Rac protein) was shown to be a component of disease resistance pathway acting downstream of R gene when OsRac1 transformed japonica rice was infected with the blast fungus, Magnaporthe oryzae [63]. Increased resistance to tobacco mosaic virus was documented when tobacco was transformed with OsRac1 [64,67]. In our qRT-PCR study, both ADP-ribosylating and Rac protein was substantially upregulated in M. graminicola infected plants than the uninfected ones at 7 dpi. Together, our data suggest changes in expression level of these candidate genes can be a contributory factor to confer resistance in rice to M. graminicola. However, expression of ARF18 TF (auxin response factor) was downregulated in nematode-infected plants compared to uninfected plants at 7 dpi in our qRT-PCR study. Auxin-regulated proteins are known to coordinate the balance between plant root growth and disease resistance by promoting the auxin biosynthesis and suppressing the benzoxazinoid-based defense compound formation [68]. The direct role of auxin influx (AUX1, LAX3) and efflux (PIN3) proteins in giant cell formation were unraveled during Arabidopsis-M. incognita compatible interaction [69].

Conclusions
In conclusion, forty accessions displaying a high degree of resistance to M. graminicola were identified. This data uncovers the potential of O. nivara, O. rufipogon and O. sativa f. spontanea as novel resources for RRKN management. GWAS was successfully applied to dissect the genetic architecture of resistance to RRKN in Indian wild rice populations. The present study provides an example of exploring the untapped genetic resources and novel genes which enrich the repository of candidate genes utilizable for future marker-assisted rice breeding program for RRKN resistance.

Rice accessions screening assay
Seeds of different rice accessions were surface-sterilized via 70% ethanol prior to soaking overnight in distilled water. Seeds were germinated in Petri dish containing wet filter paper in a growth chamber at 28 ± 2 °C and 4-5 days old seedlings were used for infection bioassays in Pluronic gel medium (PF-127, Sigma-Aldrich) as described previously [36,37]. Briefly, 60 ml of 23% PF-127 was poured into 150 × 20 mm Petri dish containing 7-9 uniformly distributed seedlings of identical genotype at < 15 °C and approximately 30 RRKN J2s were inoculated at the root tip of each seedling by a pipette tip. Dishes were incubated, plantlets were harvested at 16 dpi and nematode infection potential was assessed as described above.

Genome-wide association analysis
The set of 272 rice accessions were genotyped using in-house developed "OsSNPnks" 50K genic

Ethics approval and consent to participate
Not applicable on the manuscript.

Consent for publication
Not applicable on the manuscript.

Availability of data and materials
The data sets supporting this article are included in the article and in the additional files.     qRT-PCR-based differential expression patterns of 12 candidate genes in infected and uninfected root tissues of Pusa 1121 plants at 7 days post inoculation (dpi) with M. graminicola. Pusa 1121 was used as the reference cultivar. Data reflect gene expression levels in whole roots collected at the time of nematode inoculation (0 dpi control), in whole roots at 7 days after mock inoculation (7 dpi control) and at 7 days after nematode inoculation (7 dpi infected). Gene expression levels were normalized using the internal reference gene of O. sativa, i.e. 18S rRNA. Data are expressed as the fold change in expression.