Evaluation of blast response in the germplasm collection
A total of 290 and 278 rices, belonging to a panel of 311 japonica accessions (including 11 aus, indica and aromatic cvs as outgroups), were respectively evaluated under field and growth chamber conditions for resistance to rice blast.
The distributions of blast disease scores for the field evaluations (means of the 2013 and 2014 evaluations) are presented in Fig. 1A. The Pearson correlation value among the average SES scores of the two years (means of the two replicates of 2013 and one replicate in 2014) was 0.88 (p-value < 2.2e-16), indicating a robust relationship among the two year data. Blast evaluation was based on a 1–9 scoring system on leaves (1 = absence of symptoms and 9 = diseased area > 85%); we considered all accessions with a score below or equal to 3 as resistant; these lines showed higher numbers of small lesions with some becoming elliptical (indicating capabilities in limiting the infection spread), and the remainder as susceptible. Using this threshold, only 5.5% (16/290) of accessions were classified as resistant, while 65.5% (190/290) were highly susceptible (SES scores 8–9; Fig. 1a; Additional file 2: Table S1). The broad sense heritability of blast resistance scores in field conditions was very high (h2 = 0.85 calculated on two repetitions; Additional file 3: Table S2) as measured in the 2013 experimental design. Thus, our scoring of blast disease resistance is considered accurate even in absence of replications during 2014.
The 278 rice accessions evaluated under growth chamber conditions were artificially inoculated with a mixture of three Italian blast isolates (IT02, IT03, IT10), showing a broad pattern of virulence towards several blast R genes (Roumen et al. 1997; Urso et al. 2016). Based on a 1–5 scoring system on leaves (where accessions with a blast disease score equal to- or below 2 were considered as resistant and the remaining as susceptible), 59 accessions (21.2%) were classified as resistant (Fig. 1b; Additional file 2: Table S1). Similar to the field evaluations, a large percentage of accessions (41%) were classified as highly susceptible (score > 4). The broad sense heritability of blast resistance scores in growth chamber conditions was also high (h2 = 0.74; Additional file 3: Table S2).
Pair-wise comparisons of disease resistance between the field and growth chamber data indicate a robust relationship between the recorded phenotypes (Fig. 1c; r = 0.83), suggesting a similar response of the panel in the two conditions. With the exceptions of Brazos and Giza_177, all the accessions showing field score equal or lower than 3 also showed growth chamber score lower than 2 (in the 11 accessions for which both field and growth chamber data were available). However, of the 54 accessions resistant in growth chamber (score equal or lower than 2), only 11 showed resistance in the field (score equal or lower than 3), indicating that in field conditions additional blast strains contributed to a higher disease pressure, or that resistance expresses differently in the field. Notably, three accessions (Lamone, Dixiebelle and Giano) with growth chamber score equal or lower than 2, showed field infection scores of 7 (Additional file 2: Table S1).
Based on the evaluation of the disease in the field and growth chamber, the differences in rice blast resistance among genetic sub-populations of the panel were also analyzed (Fig. 1D). The Structure analysis, conducted at increasing K values (as proposed by Courtois et al. 2012) on the whole panel, showed K = 2 as the most probable value (Additional file 1: Fig. S1). At this K value the group of temperate japonica accessions (208) was defined, and the percentage of varieties classified as admixed was 12.5% (Additional file 2: Table S1). A principal component analysis (PCA) on genotypic data confirmed this clustering, as the two groups were separated with admixed accessions located in between, while higher K values in the Structure analysis, implemented with a Jukes and Cantor model-based neighbor-joining tree, allowed the expected distinction of tropical japonica, and a few indica, aromatic and aus accessions (Additional file 1: Fig. S1; Additional file 2: Table S1). The median of blast scores in the temperate japonica group (K1) were 8.75 and 4.00 for the field and growth chamber evaluations, respectively; while the K2 cluster, containing mainly tropical japonica, but also indica, aus and aromatic accessions,showed median scores of 5.63 and 2.17 for field and growth chamber evaluations respectively; admixed varieties scored in between the two main groups (6.00 and 2.17, respectively for field and growth chamber). The differences among subpopulations were statistically significant, as resulted from Kruskal-Wallis and Dunn’s tests (Additional File 3: Table S2). These results suggest that even though genetic distance among temperate and tropical japonica is low (Additional file 1: Fig. S1d), the latter sub-population is more resistant to Italian strains of the blast fungus.
Analysis Of Linkage Disequilibrium And Gwas For Blast Resistance
The SNP panel used in the GWAS consisted of 246,084 markers. After filtering for callrate (> 95%) and minor allele frequency (> 5%) the final set included 37,423 SNPs, with an average value of 3,119 markers/chromosome. Assuming an estimated genome size of 373 Mbp (Kawahara et al. 2013), the average density was 9.97 kbp/marker (ranging from 14.9 for chromosome 3 to 5.9 for chromosome 11). In the GWAS panel, the mean LD decay over the physical distance, computed as r2, was 934.2 kbp, ranging from 445 kbp for chromosome 11 to 1,355 kbp for chromosome 10 (Additional file 4: Fig. S2).
Although such values are higher than to those generally present in literature (Mather et al. 2007), similar or even higher LD values in rice populations were previously highlighted (Kumar et al. 2015; Xu et al. 2011). Therefore it was assumed that the panel has sufficient detection power for an association mapping approach.
GWAS was conducted considering the disease resistance scores in both field and growth chamber conditions. A total of 14 significant MTAs were identified, with -log10(p) of the peak SNP markers ranging between 3.87and 14.59 and passing a False Discovery Rate threshold of 0.05 (Table 1; Fig. 2). Of these, 8 MTAs were identified for the field data, while 6 MTAs were detected under growth chamber conditions. In almost all the MTAs, the resistant phenotype was associated with the allele in lower frequency (minor allele). However, an opposite situation was discovered for MTAs BRF06-2 and BRF10, for which the accessions carrying the most frequent allele (major allele) displayed lower average blast scores with respect to the genotypes with the least frequent allele.
Table 1
Significant associations between SNP genotypes and blast resistance in the analyzed rice accessions.
Category | MTA ID | Peak Marker/Region (bp) | Chr. | -log10(p) | #SNP | Associated region (bp) | Peak marker R2* | LD block (bp) | SNP alleles at peak marker** | Average blast score | Allelic frequence (%) |
| | | | | | Start | End | Size | | Start-End | | | Total | K1 | K2 | Adm. |
Field | BRF01 | S1_32855693 | 1 | 4.369 | 1 | - | - | - | 6.23% | 32,855,693 − 33,101,889 | A | 7.39 | 93.5 | 92.0 | 100.0 | 90.9 |
G | 5.26 | 6.5 | 8.0 | 0.0 | 9.1 |
Field | BRF04 | S4_31566281 | 4 | 14.594 | 31 | 31,257,415 | 32,476,401 | 1,218,986 | 19.54% | - | T | 7.59 | 91.4 | 96.4 | 75.9 | 91.7 |
G | 3.23 | 8.6 | 3.6 | 24.1 | 8.3 |
Field | BRF06-1 | S6_10184497 | 6 | 5.889 | 17 | 10,108,408 | 11,927,672 | 1,819,264 | 7.77% | - | C | 7.60 | 92.4 | 94.4 | 92.2 | 74.2 |
A | 4.52 | 7.6 | 5.6 | 7.8 | 25.8 |
Field | BRF06-2 | S6_29136506 | 6 | 3.963 | 1 | - | - | - | 5.08% | - | A | 7.06 | 58.8 | 57.4 | 59.6 | 64.9 |
G | 7.44 | 41.2 | 42.6 | 40.4 | 35.1 |
Field | BRF07 | S7_23126442 | 7 | 4.723 | 1 | - | - | - | 6.18% | - | A | 7.35 | 94.5 | 92.3 | 100.0 | 100.0 |
G | 4.95 | 5.5 | 7.7 | 0.0 | 0.0 |
Field | BRF10 | S10_14442878 | 10 | 4.450 | 2 | 14,421,776 | 14,442,878 | 21,102 | 5.79% | - | G | 6.75 | 58.6 | 47.9 | 79.3 | 83.8 |
A | 7.86 | 41.4 | 52.1 | 20.7 | 16.2 |
Field | BRF11-1 | S11_27464620 | 11 | 5.476 | 28 | 26,482,681 | 27,672,972 | 1,190,291 | 7.26% | 26,090,712 − 28,791,205 | C | 7.57 | 91.7 | 99.5 | 72.4 | 83.8 |
T | 4.12 | 8.3 | 0.5 | 27.6 | 16.2 |
Field | BRF11-2 | S11_28957771 | 11 | 6.304 | 2 | 28,900,124 | 28,957,771 | 57,647 | 8.43% | 28,807,359 − 29,018,344 | T | 7.64 | 83.8 | 91.3 | 56.9 | 86.5 |
G | 5.05 | 16.2 | 8.7 | 43.1 | 13.5 |
Growth chamber | BRGC01 | S1_1719121 | 1 | 3.872 | 1 | - | - | - | 5.65% | | A | 3.48 | 80.4 | 97.9 | 26.4 | 68.6 |
G | 2.71 | 19.6 | 2.1 | 73.6 | 31.4 |
Growth chamber | BRGC06 | S6_10390066 | 6 | 6.068 | 20 | 9,906,359 | 10,685,711 | 779,352 | 9.55% | | C | 3.51 | 89.1 | 94.7 | 86.5 | 63.9 |
T | 2.00 | 10.9 | 5.3 | 13.5 | 36.1 |
Growth chamber | BRGC07 | S7_22505151 | 7 | 4.167 | 1 | - | - | - | 6.14% | | G | 3.41 | 92 | 90.4 | 92.5 | 100.0 |
A | 2.47 | 8 | 9.6 | 7.5 | 0.0 |
Growth chamber | BRGC11-1 | S11_25555354 | 11 | 4.413 | 5 | 24,811,878 | 25,555,354 | 743,476 | 6.03% | 24,068,338 − 26,090,282 | T | 3.54 | 87 | 97.3 | 67.3 | 61.1 |
G | 1.95 | 13 | 2.7 | 32.7 | 38.9 |
Growth chamber | BRGC11-2 | S11_27672972 | 11 | 5.975 | 25 | 26,482,681 | 28,167,315 | 1,684,634 | 8.34% | 26,090,712 − 28,791,205 | C | 3.47 | 92.4 | 97.9 | 75.5 | 88.9 |
T | 1.65 | 7.58 | 2.1 | 24.5 | 11.1 |
Growth chamber | BRGC11-3 | S11_28957771 | 11 | 7.143 | 34 | 28,809,430 | 29,007,377 | 197,947 | 11.54% | 28,807,359 − 29,018,344 | T | 3.54 | 84.1 | 91.0 | 58.5 | 86.1 |
G | 2.25 | 15.9 | 9.0 | 41.5 | 13.9 |
For each trait. the peak marker (SNP with the highest p-value) is reported. * variance explained by the marker/region; ** Major allele is written in bold letters |
Among the MTAs identified in field, the largest allele effects (i.e. the highest difference in average blast scores among the major and minor allele of the peak marker) were observed for BRF04, BRF11-1 and BRF06-1 (4.36, 3.45 and 3.08, respectively), suggesting that gene(s) underlying these MTAs contributed most effectively to blast resistance under these conditions. The same was observed for BRGC11-2, BRGC11-1 and BRGC06 in growth chamber conditions (1.82, 1.59, and 1.51, respectively). For these loci, the highest incidence of favourable alleles at the peak markers was observed in either the K2 (identifying the tropical japonica cluster) or the admixed group (Table 1). This supports the idea that the blast resistance loci identified, in both field and growth chamber conditions, were mainly derived from tropical and admixed accessions, as already suggested from the analysis of phenotypic scores.
The SNPs showing the most significant p-values for the two chromosome 6 MTAs BRF06-1 (from 10,108,408 to 11,927,672 bp) and BRGC06 (from 9,906,359 to 10,685,711 bp) largely overlapped, even considering that the two peak markers were separated by about 205 Kbp from each other (Table 1; Fig. 3; Additional file 5: Fig. S3). On chromosome 11, two regions are defined by largely overlapping MTAs in the two experiments (Table 1; Fig. 3): BRF11-1/BRGC11-2 and BRF11-2/BRGC11-3, showing the same peak marker S11_28957771. Also, in these cases, the data suggest a common origin of resistance for the field and growth chamber experiments.
Analysis of candidate genes at the blast resistance loci
To search for candidate genes underlying the MTAs identified in this study, we considered all annotated genes available from the Oryza sativa reference sequence (Os-Nipponbare-Reference-IRGSP-1.0; RAP database: http://rapdb.dna.affrc.go.jp/download/irgsp1.html) and located within the regions defined by -/+ 100 kbp sequences from the MTA peak marker, according to the extension of the mean r2 decay). Moreover, the positions of the significant MTAs were compared with those of known blast R genes or QTLs. Results are summarized in Fig. 3 and Additional file 6: Table S3.
The genomic intervals identified by BRF06-1/BRGC06, BRF11-1/BRGC11-2 and BRF11-2/BRGC11-3 overlapped in both field and growth chamber assays. For BRF06-1/BRGC06, the interval between the two peak markers (from 10,184,497 to 10,390,066 bp) was treated as a single peak. The search on RAP-DB detected a total of 26 genes in the chromosome 6 interval from 10,084,497 to 10,490,006 bp. Of these, twenty were functionally annotated and included Os06g0286700 (from 10,387,793 to 10,390,465 bp), corresponding to the cloned blast R gene Piz and its allelic variants Piz-t, Pi2, Pi9 and Pi40 (Zhou et al. 2006; Wang et al. 2017), which localizes in close proximity of the BRGC06 peak marker (at 10,390,066 bp; Fig. 3). Two additional cloned blast R genes, Pi59(t) and Pigm(t), were also located in this region (Fig. 3).
Thirty-six genes were detected in the chromosome 11 region analyzed for the MTAs BRF11-1/BRGC11-2 (27,364,620 − 27,772,972 bp; Additional file 6: Table S3). Three of them (Os11g0678400, Os11g0684700, Os11g0684100) corresponded or were similar to NB-ARC containing proteins. These types of sequences have been described as part of the catalytic domain of R proteins, and show a high grade of polymorphism in different alleles of the blast resistance gene Pib (van Ooijen et al. 2008; Vasudevan et al. 2015). In addition, literature searching highlighted that a large cluster of Pi genes is localized in a region overlapping the BRF11-1/BRGC11-2 MTAs (Fig. 3). This group includes, among the others, Pik (Os11g0689100; 27,984,697 − 27,989,128 bp), as well as its allelic forms (Pik-m, Pik-p, Pik-s, Pik-g and the two recently described alleles Pik-e and Pik-x; Chen et al. 2015; Li et al. 2019a). The region investigated for the MTAs BRF11-2/BRGC11-3 (28,857,771 − 29,057,771 bp) included 22 annotated genes, none of which showing functions that could be clearly related to blast resistance.
Overall, the genomic interval defined by the two chromosome 11 MTAs, BRF11-1/BRGC11-2, is overlapping or in close proximity to a dense cluster of Pi genes (Fig. 3), although none of the considered R genes was coincident with the peak markers of the MTAs.
The remaining MTAs were specific to either growth chamber or open field conditions. In the BRF01 region (32,755,693 − 32,955,693 bp), no genes having a function related to pathogen resistance were detected using RAP-DB, with the exception ofa MATH domain containing protein (Os01g0775300; Additional file 6: Table S3). This type of protein is a pathogen-responsive membrane bound receptor kinase, involved in the early phases of fungal interaction in Arabidopsis (Peškan-Berghöfer et al. 2004; Shahollari et al. 2005). However, BRF01 is a single-marker MTA located approximately 242–290 Kbp from a cluster of Pi genes: Pi64 (Os01g0781200), Pish (Os01g0782100) and Pi37 (Os01g0781700). Local LD decays slowly in this region, as evidenced by a 246 Kbp block, spanning from 32,855,693 to 33,101,889 bp, and covering the entire region between BRF01 (32,855,693 bp) and Pi64 (33,098,082 − 33,103,904 bp) (Fig. 3). The remaining two genes lie in a region spanning approximately 101 Kbp, covered only by monomorphic markers in our dataset, therefore no further information can be provided; nevertheless, we cannot exclude the possibility that the LD block spans the whole Pi gene cluster in this region of chromosome 1.
Twenty-six genes, 18 of which were functionally annotated, were detected in the BRF04 region (31,157,415 − 32,576,401 bp; Additional file 6: Table S3). This list included an ABC-transporter (Os04g0620000) and a disease resistance domain containing protein (OsO4g0621500). The latter was located at ca. 10 Kbp from the peak marker, raising its possible involvement in blast resistance governed by BRF04. Moreover, co-positional relationships were observed for the MTA region defined by BRF04 and the Pi39(t) and Pi63 blast resistance loci (Fig. 3).
The BRF06-2 region contained 27 genes, 21 with a known function. Among them, Os04g0695600 encodes a RING protein with ligase E3 function. This class of enzymes is potentially involved in Magnaporthe grisea resistance in rice by influencing cell wall properties and the accumulation of hydrogen peroxide and phenolic compounds (Li et al. 2011). Moreover,the BRF06-2 region was positionally related to Pitq1 as defined by Sharma et al. (2012; Fig. 3).
Thirty-five genes were detected in the region investigated for BRF07 candidates, including Os07g0573100 and Os07g05731200, both coding for adenosine-sulphate kinases. These enzymes play a role in the synthesis of sulfur compounds, such as glutathione, which are involved in the stress response and resistance mechanisms (Wang et al. 2016b). On the same region BRGC07 (22,405,151 − 22,605,151 bp) is also localized and contains Os07g0563300, encoding for a B3 domain transcriptional repressor which regulates mechanisms of resistance to phytofagous insects in rice, by activating the salicylic acid signalling pathway and suppressing the jasmonic acid/ethylene-dependent pathway (Wang et al. 2015). These represent mechanisms that might also have a possible role in resistance to fungal pathogens. Both BRF07 and BRGC07 partially overlap with the region where the blast R gene Pi17 was mapped (22,250,443 − 24,995,083 bp; Sharma et al. 2012; Fig. 3).
Within the 200 Kbp interval surrounding the BRF10 peak marker (14,442,878 bp), 21 genes were annotated, 17 of which have a known function (Additional file 6: Table S3). Among them, Os10g0413400 encodes a Glycerol–3–phosphate acyltransferase, an enzyme homologous to AtGPAT6, which plays a role in Arabidopsis resistance to fungal pathogens (Chanda et al. 2008). No cloned or mapped blast R genes were previously identified in this region, raising the possibility that it represents a new source of blast resistance.
The chromosome 1 BRGC01 included Os01g0132100, encoding for a leucine-rich repeat containing protein. No cloned or mapped blast R genes map within this region; however, an LD block spanning the whole interval (approximately 500 Kbp) contained both BRGC01 and the gene Pit (Fig. 3).
Finally, the investigated region surrounding BRGC11-1 contained 24 genes. Fifteen of them have a known function. Four genes have an annotation correlated with disease resistance: Os11g0644100 (Leucine-rich repeat N-terminal domain containing protein), Os11g0647600 (Plant disease resistance response-related protein) and the two NB-ARC domains containing genes Os11g0648600 and Os11g0649100. Moreover, the associated regions is positionally related with Pikh (and its allelic form Pi54; Sharma et al. 2005).
Validation of BRGC11-1 on chromosome 11 in a biparental segregating population
BRGC11-1 was an MTA providing a high allelic effect (i.e. blast score difference between major and minor allele: 1.59) under growth chamber conditions. To further dissect and validate this region, a biparental F2/F3 segregating population derived from a cross between the resistant accession Salvo (showing the resistant allele G at BRGC11-1) and the susceptible accession Maratelli (showing the susceptible allele T; Additional File 2: Table S1), was used. The disease reactions of Salvo and Maratelli and the 130 F2/F3 lines derived from the cross were subjected to artificial inoculation with a mixture of three blast isolates (IT02, IT03, IT10) under growth chamber conditions. Infection and disease scores were evaluated following the procedures adopted for the GWAS of the 278 accessions. The analysis of variance, calculated using the phenotypic data revealed a highly significant effect of the genotype (P < 0.0001) on the disease incidence, with no significant differences observed between replications (P = 0.38). No trangressive lines (with extreme resistance/susceptibility phenotypes) were observed, suggesting that resistance factors were contributed only by Salvo. The estimated broad sense heritability for resistance within the Salvo x Maratelli population was 92.5%, indicating that the majority of the phenotypic variance was due to genetic effects.
A total of 3,724 SNP markers together with 46 SSRs were used to assemble the genetic linkage map. After elimination of the unlinked loci, the genotypic data relating to 1,669 informative marker loci were assembled into 15 linkage groups corresponding to the 12 rice chromosomes (Additional file 8: Table S5A). Two linkage groups were identified for chromosomes 2, 4 and 5. The overall length of the map was 1,798.33 cM with individual chromosome genetic length ranging from 120.2 cM (chromosome 6) to 195.9 (chromosome 11) and average chromosome length of 149.9 cM. The total number of mapped loci per chromosome ranged from 67 (chromosome 8) to 226 (chromosome 1) with an average of 139.1 loci per chromosome. The genome-wide marker density was 1.08 cM/marker, varying from 0.68 cM/marker (chromosome 1) to 2.25 cM/marker (chromosome 8) (Additional file 8: Table S5B).
The analysis of blast resistance allowed the identification of one highly significant QTL (SxMBRQTL-1) on chromosome 11, with a LOD score of 29.81 and explaining 74% of the phenotypic variation. The QTL peak was identified by the SSR marker RM224 and, based on Darvasi and Soller (1997), the confidence interval covered a region of 5.85 cM, from S11_25478006 to S11_26583984. The negative value of the additive effect suggests that the resistance allele at this locus was derived from Salvo (Additional file 8: Table S5C). The region defined by this QTL overlapped with the peak marker of BRGC11-1 (S11_25555354, 25,555,354 bp), and partially overlapped with regions defined by BRGC11-2 (26,482,681 − 28,167,315 bp) and BRF11-1 (26,482,681 − 27,672,972 bp).