Genome-wide association analysis identifies resistance loci for bacterial blight in a diverse collection of indica rice germplasm

Bacterial blight, which is caused by Xanthomonas oryzae pv. oryzae (Xoo), is one of the most devastating rice diseases worldwide. The development and use of disease-resistant cultivars have been the most effective strategy to control bacterial blight. Identifying the genes mediating bacterial blight resistance is a prerequisite for breeding cultivars with broad-spectrum and durable resistance. We herein describe a genome-wide association study involving 172 diverse Oryza sativa ssp. indica accessions to identify loci influencing the resistance to representative strains of six Xoo races. Twelve resistance loci containing 121 significantly associated signals were identified using 317,894 single nucleotide polymorphisms, which explained 13.3–59.9% of the variability in lesion length caused by Xoo races P1, P6, and P9a. Two hotspot regions (L11 and L12) were located within or nearby two cloned R genes (xa25 and Xa26) and one fine-mapped R gene (Xa4). Our results confirmed the relatively high resolution of genome-wide association studies. Moreover, we detected novel significant associations on chromosomes 2, 3, and 6–10. Haplotype analyses of xa25, the Xa26 paralog (MRKc; LOC_Os11g47290), and a Xa4 candidate gene (LOC_11g46870) revealed differences in bacterial blight resistance among indica subgroups. These differences were responsible for the observed variations in lesion lengths resulting from infections by Xoo races P1 and P9a. Our findings may be relevant for future studies involving bacterial blight resistance gene cloning, and provide insights into the genetic basis for bacterial blight resistance in indica rice, which may be useful for knowledge-based crop improvement.


Introduction
Bacterial blight, which is caused by Xanthomonas oryzae pv.oryzae (Xoo), is one of the most devastating diseases of cultivated rice (Oryza sativa L.) in tropical and temperate regions worldwide [1].This disease is frequently prevalent in southern China and southeast Asia, resulting in heavy rice yield losses [2,3].Developing and deploying resistant cultivars carrying major resistance (R) genes has been the most effective approach for managing bacterial blight [4].
Novel R genes and QTLs associated with bacterial blight resistance have been identified and used in breeding programs.However, the rapid loss of bacterial blight resistance in rice varieties carrying a single R gene remains a problem for breeders.This has been in part because the mechanisms mediating the pathogenesis of Xoo and the genetic basis for bacterial blight resistance in rice have not been fully characterized.To date, most of the studies on bacterial blight resistance have been based on a single resistant parent or bi-parental genetic mapping populations.Additionally, the genetic variability of bacterial blight resistance in rice accessions has yet to be addressed because only one or two parental lines have been used in previous studies.A genome-wide association study (GWAS) involving high-density single nucleotide polymorphisms (SNPs) based on next-generation sequencing may be useful for detecting genetic variants that can be directly applied to improve rice cultivars [26].In this study, 172 global indica accessions were inoculated with representative strains of six Xoo races from China and the Philippines to evaluate their resistance reactions.We also conducted a GWAS of bacterial blight resistance based on the genotyping of 700,000 SNPs with a high-density rice array [27].Our study objectives were to elucidate the genetic basis of bacterial blight resistance and identify loci related to bacterial blight resistance in indica rice lines, which may provide useful information for improved rice production.

Plant materials and bacterial inoculations
We examined 172 indica rice accessions from 26 countries (S1 Table ).To evaluate bacterial blight resistance, the seeds of all plant materials were sown in a seedling nursery, and 30-dayold seedlings were transplanted to the experimental farm at the Institute of Crop Sciences, Chinese Academy of Agricultural Sciences, Beijing, China.There were nine plants in each row (20 × 17 cm).The following six representative Xoo strains were used to artificially inoculate plants: GD1358 (race C5) and V (race GV) from China, and PXO61 (race P1), PXO340 (race P3c), PXO339 (race P9a), and PXO99 (race P6) from the Philippines.The bacterial races were grown on peptone sucrose agar medium at 30˚C for 2 days, and each inoculum was prepared by suspending the bacterial mass in sterile water at a concentration of 10 8 cells ml −1 .Five central plants for each line were inoculated with each Xoo race (with two replicates) at the tillering stage (plant age: 65 days), and four or five of the uppermost leaves of each plant were inoculated with each Xoo race using the leaf-clipping method [28].The lesion lengths (LLs) were measured on all inoculated leaves 3 weeks post-inoculation when lesions were obvious and stable.The average LL of one accession was calculated based on 15 of the longest lesions from five individual plants (i.e., three lesions per plant) for each replicate.The average LL of two replicates for each accession was used as the phenotype for the GWAS.

Genotyping and population structure analysis
All 172 accessions were genotyped for 700,000 SNPs using a high-density rice array [27].To avoid the influence of linked SNPs during the population structure analysis, we used the linkage disequilibrium (LD) pruning tool of PLINK [29] to obtain a subset of 45,311 independent SNPs with a minor allele frequency (MAF) > 5% and a missing data ratio (MDR) < 0.1 according to 'indep-pairwise 50 10 0.5'.We used PHYLIP (version 3.6) [30] to construct an unrooted neighbor-joining tree with 100 bootstrap replicates.The genetic structure of the whole population was predicted with the ADMIXTURE program [31].EIGENSOFT [32] was used to conduct principal components analysis (PCA) for estimating the number of subpopulations in the GWAS panel.

Association mapping
A total of 279,266 SNPs with a MAF > 5% and MDR < 0.2 were filtered for association analyses of the whole panel.The GWAS was completed using a linear mixed effects model to determine the associations between each SNP and the phenotype (i.e., mean lesion length caused by Xoo).We applied the efficient mixed model analysis (EMMA) feature of the EMMA eXpedited (EMMAX) software [33].We used the Balding-Nichols matrix based on genome-wide SNP data to develop the kinship matrix, which measured the genetic similarity between individuals.The first seven principal components explaining 80% of the genetic variation were used as covariates.
Based on 1000 permutation tests, the genome-wide significance thresholds at a significance level of 0.05 were P = 3.49 × 10 −7 , 1.82 × 10 −7 , 2.58 × 10 −6 , 3.15 × 10 −7 , 3.50 × 10 −7 , and 1.34 × 10 −7 for Xoo races P1, P3c, P6, P9a, C5, and GV, respectively.The Manhattan and quantilequantile plots for GWAS results were created using the R package qqman [34].We obtained independent association signals using a published method [35].Briefly, multiple SNPs exceeding the threshold in a 5-Mb region were clustered using an LD r 2 value !0.25, and SNPs with the minimum p value in a cluster were considered as the lead SNPs.The pairwise LD r 2 values were calculated with PLINK [29], and the R package LDheatmap [36] was used to draw the heatmap of pairwise LD.The SNPs spanning a region from 2 kb beyond the 3 0 end of a gene to the 5 0 end were concatenated as the haplotype.We excluded SNPs with no functional effect according to the haplotype analysis, and only used the haplotypes shared by at least five accessions.Phenotypic variations explained by multiple lead SNPs of association loci were estimated by stepwise regression using SAS [37].

Analysis and annotation of significant signals
Synonymous and nonsynonymous SNPs and SNPs with large-effect changes were annotated based on the gene models of the annotated version of the Nipponbare reference genome (IRGSP 1.0) [38] using snpEff software (version 4.0) [39].Enriched gene ontology terms were identified using the agriGO online tool [40].We listed all significant SNPs located within genes and annotation information for the Nipponbare reference genome (IRGSP 1.0) [38].

Population structure of rice accessions
We genotyped a diverse global collection of 172 O. sativa ssp.indica accessions (S1 Table ) using a high-density rice array [27], resulting in 700,000 high-quality SNPs.A total of 45,311 independent SNPs with MAF > 5% and MDR < 0.1 were used for genetic structure analyses.The neighbor-joining tree developed using PHYLIP version 3.6 [30] with 100 bootstrap replicates and the ADMIXTURE results [31] revealed that the 172 accessions were grouped into three clusters (S1A and S1C Fig) .A similar result was observed using PCA, in which most of the genetic variation in the accessions was explained by the first two principal components.When we plotted the first two components against each other, most accessions clustered in three groups with a few admixed accessions (S1B Fig) .We then divided the 172 indica accessions into three subgroups, namely ind-I (78 accessions), ind-II (36 accessions), and ind-III (58 accessions) (S1 Table ).

Reactions of accessions to different Xoo races
All accessions were inoculated with six Xoo races at the primary tillering stage (S1 Table ).Based on the LL values calculated using the unweighted pair-group method with arithmetic mean, the six bacterial races were divided into three groups.One group consisted of race P9a, another group comprised races P1 and P6, and the third group included the remaining three races, which were more pathogenic than the races in the other groups (Fig 1A).Additionally, the accessions were clustered into three major groups (i.e., indicated in green, red, and blue in Fig 1A) based on the LL values for the six Xoo races.The accessions in the red group mainly consisted of indica-I accessions, and were equally resistant to race P1.The resistance level of each accession was based on the average LL for ten plants (two replicates) as follows: LL < 5 cm, 5 cm LL < 10 cm, 10 cm LL < 15 cm, and LL ! 15 cm represented resistant, moderately resistant, moderately susceptible, and susceptible levels, respectively.A considerable proportion (i.e., from 94.8% for P6 to 98.2% for P3c) of accessions were moderately susceptible or susceptible to Xoo races P6, P3c, C5, and GV, respectively (Fig 1B).However, there were similarly high proportions of accessions that were resistant and susceptible to both P1 and P9a, suggesting the R genes with large effect on P1 and P9a resistance were carried by only some accessions.The wide range of LLs observed in different rice germplasm accessions indicated there was substantial genotypic variability associated with resistance to Xoo races, especially for P9a and P1 (Fig 1C).For all Xoo races (except P9a), the average LL of subgroup indica-I was significantly smaller than that of the other two subgroups.For P9a, the average LL of subgroup indica-II was significantly higher than that of subgroups indica-I and indica-III (p < 0.001 for both; Fig 1C).

Genome-wide association analysis for Xoo resistance in indica rice
We used association mapping to identify genome-wide associated loci underlying the resistance to six Xoo races.A total of 279,266 SNPs with a MAF > 0.05 and MDR < 0.2 in the whole population and subpopulation (i.e., whole population without the admixed accessions highlighted in S1B Fig) were used for the association analyses with a linear mixed effects model in the EMMAX program [33].We constructed the empirical kinship (K) matrix using the BN matrix with 211,943 SNPs (MAF > 0.01 and MDR < 0.05), and covariates (Q) were generated based on the PCA of the population structure.The similarities in the Manhattan plots and quantile-quantile plots for the whole population (Fig 2 ) and the subpopulation (S2  ), explaining 55.3%, 59.9%, and 13.3% of the variation in LL, respectively.Significantly associated SNPs were not identified for bacterial races P3c, C5, and GV.The 121 detected SNPs were distributed in 96 annotated genes as follows: 60 in intergenic regions, 14 in introns, one in a promoter, 13 in missense variants, 27 in synonymous variants, two in 5 0 untranslated regions, two in 3 0 untranslated regions, and two in stop codons (S2 Table ).Of the 121 significant signals, only SNP rs11_ 28436056 was detected for both P1 and P6.This SNP was located in a Xa26 family gene (LOC_Os11g47290), which encodes a receptor kinase.Gene ontology analysis indicated that these 96 genes were significantly enriched for defense responses, programmed cell death biological processes, and purine ribonucleotide-bindingrelated functions (S3 Table ).

Regions strongly associated with bacterial blight resistance
The physical distance between neighboring SNPs was 1-131,763 bp (average: 1349 bp) in our GWAS panel.Thus, we combined adjacent significantly associated SNPs within an LD block (i.e., multiple SNPs exceeding the significance threshold in a 5-Mb region were clustered using an LD r 2 value !0.25) as a resistance region.Additionally, the SNPs with the lowest p values in the LD block were considered the lead SNPs.We detected 12 strongly associated regions containing 120 significant SNPs (including one common SNP for P1 and P6) for three of six Xoo races (Table 1).Among these, one hotspot region was associated with resistance to P9a, and spanned an approximately 1.0-Mb interval (i.e., 16,502,066-17,531,046 bp) on chromosome 12, which included 58 (85.3%) of 68 significantly associated SNPs (Fig 3A).Nine significant SNPs (lead SNP rs12_17304839, P = 1.8 × 10 −12 ) in this region overlapped with xa25 (LOC_Os12g29220) (S2 Table ), which is a recessive R gene that is responsible for race-specific resistance to P9a in rice [17].Another hotspot region (i.e., L11), which included 48 significant SNPs associated with resistance to P1, was located at the end of chromosome 11.We determined that 43 (89.6%) of the 48 SNPs in L11 were located in an approximately 1.45-Mb interval (i.e., 27,252,984-28,704,769 bp) (Fig 4A).The corresponding region contained 33 genes that encoded proteins mainly involved in defense responses (e.g., NB-ARC domain-containing protein and protein kinase) according to the Nipponbare reference genome.This region carried the well-known R genes Xa4 and Xa26 [16,41].The lead SNP (i.e., rs11_28142810, P = 2.2 × 10 −15 ) of region L11 was located in the coding sequence (CDS) of LOC_Os11g46870, which encodes a protein kinase.This location corresponded to the Xa4 gene region based on the physical position of the BAC clone 3H8 that carries Xa4 [41].The second leading SNP (i.e., rs11_28437434, P = 1.3 × 10 −14 ) of region L11 was located in the CDS of LOC_Os11g47290 (i.e., encoding a receptor kinase), which is the Xa26 paralog MRKc [42] (S2 Table ).
We did not detect any significant association signals or a resistant haplotype associated with LOC_Os11g47210 (data not shown).However, we identified one strong association signal for P1 located at around 28,437,434 bp (i.e., rs11_28437434, P = 1.3 × 10 −14 ) on chromosome 11.There were also five significant SNPs for P1 and one significant SNP for P6 located in the CDS region of LOC_Os11g47290 (encoding a receptor kinase), which is a Xa26 paralog (i.e., MRKc) [42] (Fig 4A and S2 Table).We mainly analyzed the LOC_Os11g47290 haplotype, and  ).These haplotypes explained 45.3% and 32.2% of the phenotypic variations in the lesions  ).Our results suggest that chromosomes 11 and 12 were important for the evolution of rice disease resistance.At the genome level, loci associated with bacterial blight resistance exhibited race specificity, similar to some reported R genes and QTLs [1,24,25].In this study, only two SNPs with weak but significant signals associated to resistance for P6 were detected (Table 1 and Fig 2), no any SNPs significantly associated resistance for Xoo races P3c, C5, and GV was identified in the tested rice accessions due to limited phenotypic variation (Fig 1), suggesting lack of R gene/QTL with large effect on resistance for these races in the association panel.It is generally accepted that GWAS has low power to find associations for rare alleles [44].To clarify whether low allele frequency resulted in false positive associations in GWAS, further studies will be needed to ascertain the functions of resistance candidate genes by expression profiling or genetic transformation.
Our GWAS and haplotype results provided some interesting information on the differentiation and accumulation of resistance loci for Xoo in rice.Bacterial blight is prevalent in both tropical and temperate areas, but outbreaks and epidemics of this disease frequently occur in the tropics because of a favorable climate for the survival, propagation, and infection of Xoo [1,45].During the co-evolution of plant hosts and pathogens, the strong selection pressure resulting from diverse pathogen populations might have induced indica rice lines to generate novel disease resistance specificities or evolve novel effectors to overcome the effects of genes associated with disease susceptibility.).Nipponbare reference genome indicated that this japonica rice variety carries xa25 [43] as well as the associated disease resistance haplotype.These results suggested that the effects of the xa25 alleles in indica accessions were overcome by P9a during the co-evolution between the host and pathogen.
Breeding for bacterial blight resistance has been one of the most important objectives among rice breeders in Asia since the 1970s.Rice germplasm carrying Xa4 and Xa3/Xa26 exhibit a broad resistance to Chinese Xoo races, and have been obtained from the International Rice Research Institute and used as the bacterial blight resistance donor in rice breeding programs [3].The results in study suggest that breeding for bacterial blight resistance involved the end of the long arm of chromosome 11 (Fig 2 and Table 1), which carries several bacterial blight resistance genes (e.g., Xa4, and Xa26).The haplotype analysis indicated that accessions with P1 resistance haplotypes related to the Xa4 candidate gene (i.e., LOC_Os11g46870), Xa26 paralog gene (i.e., LOC_Os11g47290), and a P9a-specific susceptible haplotype of the xa25 gene (i.e., LOC_Os12g29220) were associated with differential disease susceptibilities among the indica subgroups (Fig 3 and Fig 4).The P1 resistance haplotypes associated with the Xa4 candidate gene (i.e., LOC_Os11g46870) and Xa26 paralog gene (i.e., LOC_Os11g47290) were more frequently detected in the indica-I subgroup than in the indica-II and indica-III subgroups (Fig 4B and 4C).This implies that conventional breeding has applied a strong selection pressure for certain R genes, and the breeding process for disease resistance of certain varieties has contributed to the accumulation of BB resistance loci.

Putative roles of transposable elements affecting bacterial blight resistance in rice
We identified several retrotransposons and transposons among the candidate genes anchored by the SNPs associated with bacterial blight resistance, including significant SNPs (S2 Table ).Transposable elements have been suggested to contribute to the evolution of genes by providing cis-regulatory elements leading to changes in expression patterns [46][47][48][49].Stressinduced changes to retrotransposons may play a role in generating host genetic plasticity in response to environmental stresses [50][51][52].In rice, Xa21 encodes a receptor-like kinase and confers resistance to Xoo.This gene is a member of a multi-gene family in which 17 transposon-like elements have been identified in the 5 0 and 3 0 flanking regions and introns [49,53,54].Researchers have attempted to diversify the rice bacterial blight resistance genes.Efforts have involved point mutations [55], deletion and duplication of intragenic DNA repeats encoding blocks of leucine-rich elements subject to adaptive selection [56], intragenic and extragenic recombinations [57], and insertion of transposable elements [49].Our results suggest that transposable elements have highly enriched the diversity of bacterial blight resistance genes during evolution.Additionally, a large-scale attempt to identify and annotate loci provides useful insights into the genetic control of indica rice traits.

Candidates for fine-mapped genes and novel loci
We identified 12 bacterial blight resistance loci containing 121 significantly associated signals using 317,894 SNPs.Two hotspot regions (i.e., L11 and L12) were located within or nearby cloned R genes (i.e., xa25 and Xa26) and a fine-mapped R gene (i.e., Xa4).These observations revealed the relatively high resolution of a GWAS involving a relatively large population, which increased our ability to investigate genetic diversity and a high-density SNP map.Additionally, our haplotype analysis results suggest that the Xa4 candidate gene (i.e., LOC_Os11g46870) and the Xa26 paralog (i.e., LOC_Os11g47290) likely confer resistance to P1.Moreover, two novel loci on chromosomes 8 [i.e., L6 (rs8_3317231, P = 3.3 × 10 −12 ) and L7 (rs8_9671192, P = 4.2 × 10 −8 )] were significantly associated with resistance to P9a and P1, respectively.These results provide valuable information for future studies involving bacterial blight resistance gene cloning.Moreover, the identified diverse rice accessions carrying fine-mapped and novel loci from nearly 20 countries will provide more available donors in breeding programs aiming at developing bacterial blight resistance in different rice growing regions.

Conclusions
Our study provides new insights into the genetic basis of the evolution of bacterial blight resistance in rice.The findings reported herein may be useful for knowledge-based crop improvement.Future research will focus on validating the effects of these candidate genes and their functional variants.We will use genetic transformations and DNA insertion mutant screens to verify that these genes confer bacterial blight resistance to rice.

Fig 1 .
Fig 1. Susceptible and resistant reactions of a whole population and six subpopulations of 172 indica accessions inoculated with representative strains of six Xanthomonas oryzae pv.oryzae (Xoo) races from China and the Philippines.(A) Hierarchical cluster of accessions and races based on lesion length (LL).(B) Number of accessions in the following reactions to six Xoo races: resistant (LL < 5 cm), moderately resistant (5 cm LL < 10 cm), moderately susceptible (10 cm LL < 15 cm), and susceptible (LL ! 15 cm).(C) Boxplots for LLs following infections by six Xoo races in three indica subgroups divided by population structure analysis results.Box edges represent the 0.25 and 0.75 quantiles with median values indicated by bold lines.'***' refers to a significant difference in the average LLs among the indica subgroups (p < 0.001).We used the following Xoo races: C5 (strain GD1358) and GV (strain V) from China, and P1 (strain PXO61), P3c (strain PXO340), P9a (strain PXO339), and P6 (strain PXO99) from the Philippines.https://doi.org/10.1371/journal.pone.0174598.g001

Fig 2 .
Fig 2. Manhattan and quantile-quantile plots for a genome-wide association study of bacterial blight resistance in indica rice.(A, G) P1 (strain PXO61).(B, H) P6 (strain PXO99).(C, I) C5 (strain GD1358).(D, J) P3c (strain PXO340).(E, K) P9a (strain PXO339).(F, L) GV (strain V).The strength of the associations for the lesion lengths caused by six Xanthomonas oryzae pv.oryzae (Xoo) strains is indicated as the negative logarithm of the p value for the linear mixed effects model.Based on 1000 permutation tests, the genome-wide significance thresholds (i.e., horizontal red lines in all Manhattan plots) at

Fig 3 .
Fig 3. Hotspot region for the resistance to Xanthomonas oryzae pv.oryzae race P9a and haplotype analysis of the peak associated with the gene on chromosome 12.(A) Local Manhattan plot (top) (16.5-17.5 Mb) and linkage disequilibrium heatmap (bottom) (17.2-17.5 Mb) surrounding the hotspot region on chromosome 12.The arrow indicates the position of the peak single nucleotide polymorphism (SNP) located in xa25 (LOC_Os12g29220).Dashed lines indicate the xa25 region.(B) Gene structure and haplotype analysis of xa25 in 140 accessions based on five significant SNPs in xa25.Haplotypes with fewer than five accessions are not shown.(C) Lesion lengths caused by P9a infections of accessions in three haplotypes of xa25 in different indica subgroups.Box edges represent the 0.25 and 0.75 quantiles with median values indicated by bold lines.Whiskers extend to data no more than 1.5-times the interquartile range, and the remaining data are represented by dots.'***' refers to a significant difference based on Duncan's multiple comparison tests (p < 0.001).https://doi.org/10.1371/journal.pone.0174598.g003

Fig 4 .
Fig 4. Hotspot region for the resistance to Xanthomonas oryzae pv.oryzae race P1 and haplotype analysis of the peak associated with the gene on chromosome 11.(A) Local Manhattan plot (top) (23-29 Mb) and linkage disequilibrium heatmap (bottom) (28.0-28.8Mb) surrounding the hotspot region on chromosome 11.Red arrows and points indicate the positions of the peak single nucleotide polymorphisms located in the Xa4 candidate gene (i.e., LOC_Os11g46870) and Xa26 paralog (i.e., LOC_Os11g47290), respectively.Dashed lines indicate the xa25 region.(B) Gene structure and haplotype analysis of the Xa4 candidate gene (i.e., LOC_Os11g46870).(C) Gene structure and haplotype analysis of the Xa26 paralog (i.e., LOC_Os11g47290).Lesion lengths caused by P1 infections of accessions in three haplotypes of LOC_Os11g46870 (D) and LOC_Os11g47290 (E) in different indica subgroups.'***' refers to a significant difference based on Duncan's multiple comparison tests (p < 0.001).https://doi.org/10.1371/journal.pone.0174598.g004 Considerable selection pressures from various Xoo races have contributed to the diversification of disease resistance loci in indica.This has resulted in the differentiation of bacterial blight resistance among the indica subgroups.In this study, 45.1% and 47.9% of the xa25 haplotype accessions resistant to P9a belonged to indica-I and indica-III subgroups, respectively (Fig 3B and S1 Table).However, 92.6% of the susceptible xa25 haplotype accessions (i.e., Hap2) belonged to the indica-I subgroup, in which 13 accessions resistant to P1 carried LOC_Os11g46870 or LOC_Os11g47290 resistance haplotypes (Fig 4 and S1 Table

Table and Table 1 )
. Another distinct peak comprising four SNPs (lead SNP rs8_3317231, P = 3.3 × 10 −12 ) in region L6 located on chromosome 8 was significantly associated with resistance to P9a.This region contained the EARLY FLOWERING 7 gene (LOC_Os08g06070) (Table 1, Fig 2 and S2

Table 1 . Twelve regions with signals significantly associated with bacterial blight resistance based on a genome-wide association study involving indica rice. Locus ID Chromosome Lead SNP position (bp) LD block No. of significant SNP in LD block Xoo races Most significant p value MSU ID of genes harboring lead SNP Known R genes
6% of Hap2 accessions belonged to the indica-I subgroup, while 64.3% of Hap1 accessions belonged to the indica-II subgroup.Meanwhile, 45.1% and 47.9% of the Hap3 accessions belonged to the indica-I and indica-III subgroups, respectively.These results explained the phenotypic variability among the three indica subgroups in response to Xoo race P9a (Fig 1C). https://doi.org/10.1371/journal.pone.0174598.t00192.