Rapid identification of rice blast resistance gene by specific length amplified fragment sequencing

ABSTRACT

Excavation of resistance genes is one of the most effective and environment-friendly measures to control the devastating rice disease caused by Magnaporthe oryzae. Many resistance genes have been mapped and characterized in the last century. Nevertheless, only a few of the total resistance genes could be really applied in the rice breeding program. Huazhan (HZ) is a new native rice restorer line developed in China and widely used in hybrid rice in recent years. HZ and its crossed combinations usually show a broad spectrum of resistance against rice blast in different rice ecosystems in China. Dissection of the genetic background of HZ is very useful for its further application. In this study, a combined method based on bulked segregation analysis (BSA) and specific length amplified fragment sequencing (SLAF-seq) was used to identify blast resistance gene(s) in HZ. A total of 56,187 SLAFs labels were captured and 9051 polymorphic SLAFs markers were analysed and procured in this study. One trait associated with candidate resistance genes region on chromosome 12 overlapping 10.2–17.6 Mb has been identified, in which 10 NBS-LRR (nucleotide-binding site-leucine-rich repeat) coding genes were used as resistance gene candidates. Our result indicated that SLAF-seq with BSA is a rapid and effective method for initial identification of blast resistance genes. The identification of resistance gene in HZ will improve its molecular breeding and resistance variety application.

KEYWORDS:

• Rice blast
• resistance gene
• gene identification
• SLAF-seq

Introduction

Many genetic markers including morphological markers, biochemical markers and DNA markers such as RFLP, RAPD, AFLP, STS, SSR, CAPS, etc. have been used to identify different qualitative genes and quantitative trait locus (QTL).[1–3] DNA markers based on different polymorphism-detecting techniques have promoted the gene mapping and marker-assisted selection (MAS) breeding.[4] Among the genetic markers, Single nucleotide polymorphisms (SNPs) are more serviceable than other conventional markers because of its abundance and stability.[5] Utilization efficiency of DNA markers needs to improve since most of the classical DNA markers are ordinarily used as single locus to map qualitative loci. High-throughput sequencing technologies offer more efficient methods for sequence-based SNP genotyping, which is more helpful for genetic mapping and genome analysis.[6] Combination of specific length amplified fragment sequencing (SLAF-seq) and reduced representation library (RRL) construction has provided an efficient method for large-scale de novo SNP discovery and genotyping.[7] For example, Chen et al. (2013) developed 89 SLAF markers in Thinopyrum elongatum.[8] Xia et al. (2015) identified a new maize inflorescence meristem mutant based on association analysis using the SLAF-seq method.[9] Bulked segregant analysis (BSA) method was considered to be an elegant method of mapping resistance genes.[10,11] Many blast resistance genes such as Pi37, Pi39, Pi41, Pi50, Pi58 and Pi61, etc. were mapped by this traditional method.[12–17] However, the availability of molecular markers limits gene fine mapping in some cases, because of the poor polymorphism among similar genetic background. Furthermore, genotyping of each marker in the whole population is still costly and time-consuming. Whole genome sequencing of individuals pools (Pool-seq) provides a cost-effective alternative method for gene identification.[18] Pool-seq method has been increasingly used for population genomic research on both model and non-model organisms.[19–21] The Pool-seq technique has also been applied for rice blast resistance QTL mapping.[22]

Rice blast caused by the ascomycete fungus Magnaporthe oryzae could infect the leaves, leaf collars, panicles, culms and culm nodes, which is known to occur in more than 85 countries worldwide.[23,24] It is estimated that destroyed rice by blast could be enough feed for more than 60 million people each year.[25] Blast disease has never been eradicated from a region around the world. In order to control blast disease effectively, an integrated management program was ordinarily implemented to avoid overuse of a single control method. Application of resistant varieties is always considered to be one of the most effective and environmental methods to control blast disease, since chemical control such as carpropamid-type fungicide could pollute the ecosystem.[26,27] Between the qualitative and quantitative blast resistance, QTLs are usually deemed to be more durable than qualitative resistance genes, but the qualitative resistance genes were utilized more easily than that of quantitative resistance in rice breeding,[28] thus breeders are inclined to exploit resistance genes in the rice breeding program.

Although more than 100 blast resistance genes have been identified or mapped based on genetic or molecular marker method [28–31] and 27 of which have already been characterized,[32–34] only a few resistance genes such as Pi2, Pi9, Pi39, Pi40, Pi50 and Pikh, etc. have been reported to possess broad-spectrum resistance,[35–38] while most of them confer race specific resistance including Pia, Pib, Pii, Pikm, Pit, Pi12 and Pi19, etc.[38,39] It is still indispensable to excavate more durable resistance genes from more germplasm, since blast resistance would often break down for the reason that rare pathotypes increase or novel pathotypes emerge.[40–42] In this study, a rice blast resistance gene derived from the restorer line Huazhan (HZ) was identified and mapped based on SLAF-seq method.

Materials and methods

Plant material

HZ is a new restorer line developed in China, which has been used widely in recent years. The hybrid rice combinations derived from HZ show high blast resistance in China. An F₂ population with 219 individuals derived from Lijiangxintuanheigu (LTH, susceptible parent) and HZ was used for genetic analysis and resistance gene mapping. All the seeds including two parents and their corresponding 8 F₁ and 219 F₂ plants were evenly sown in porcelain trays (30 × 20 × 5 cm) filled with yellow mud in greenhouse. Soil was kept moist with water by auto-spraying equipment. Plants were fertilized with 3–5 g of (NH4)₂SO₄/m² before inoculation at 3–4 leaf stage.

Fungal material

The M. oryzae isolate GD10-3001 which was collected from South China in 2010 showed stable pathogenicity. GD10-3001 was separated from the original strain through single spore separation method and then conserved on sterilized rice stems in a refrigerator. The isolate was transferred to yeast culture and incubated in a biochemical incubator at 25 °C for 7–10 days referred to Zhu et al.[15] Then the mycelium was shifted to sterilized maize grains and cultured under purple light at 25 °C for 10–13 days. The isolate sticked on maize was transferred to sterilized porcelain plates (25 cm × 19 cm × 2 cm) covered with a piece of wet gauze and incubated at 25 °C for 3–4 days. The conidium was washed from the maize by sterilized water referred to Chen et al.[29] The concentration of spores was adjusted to 1 × 10⁵/mL mixed with 0.2% of Tween-20 before inoculation.

Rice pathogen inoculation and resistance evaluation

Rice pathogen inoculation was performed as described previously.[15] Seven days later, blast symptoms were evaluated using a 0–9 scale standard based on the lesion types described by Standard Evaluation System for rice (SES, International Rice Research Institute (IRRI) 1996). Scores identified as 0–3 were considered to be resistant, scores with 4–9 were considered as susceptible.

Rice DNA extraction

Total genomic DNA was extracted from young leaves of F₂ population and their parents according to the cetyltrimethyl ammonium bromide (CTAB) method.[43] The total DNA was detected with NanoDrop ND-1000 Spectrophotometer (NanoDrop Technology, Wilmington, DE, USA) and diluted to 50 ng/ul. RNase A was used to remove RNA contamination from genomic DNA.

SLAF library construction and high-throughput sequencing

The reference genome of japonica rice Nipponbare, with the updated version build 5.0, was downloaded (http://rgp.dna.affrc.go.jp/IRGSP/Build3/build3.html). A pre-design SLAF experiment was performed to build up conditions to maximize SLAF-seq efficiency. First, we estimated the restriction lengths of fragments using the drilling data derived from the rice reference genome. The criteria were consulted as followed: suitable SLAFs numbers for gene mapping, average distribution of SLAFs and avoidance of repeated SLAFs. To ensure the sequencing depth consistency, a tight length range of about 30–50 bp was selected. A pilot polymerase chain reaction (PCR) amplification was performed to detect the RRL features in target length range. The pre-design step of the scheme was adjusted when non-specific amplified bands detected on the gel. The pilot experiment was to assure that the SLAFs were evenly distributed across the rice genome. The SLAF library was constructed according to the optimized pre-designed scheme. The procedure was carried out based on Sun et al. (2013) with minor adjustments.[7] All of the purified DNA samples were incubated at 37 °C with MseI (New England Biolabs, NEB), T4 DNA ligase (NEB), ATP (NEB) and MseI adapter. Restriction-ligation reactions were heat-inactivated at 65 °C first, and then the samples were digested with another restriction enzyme EcoRI at 37 °C. PCR reaction was performed with mixtures of diluted restriction-ligation samples, dNTP, MseI primer containing barcode1 and Taq DNA polymerase (NEB). PCR products were purified by an E.Z.N.A.® Cycle Pure Kit (Omega) before pooled. All the purified products were run on a 2% agarose gel. Fragments of about 330–380 bp in size were isolated using a Gel Extraction Kit (Qiagen, Hilden, Germany). The fragment products were subjected to PCR reaction with Phusion Master Mix (NEB) and Solexa amplification primer mix before adding barcode2. The target products were gel-purified and diluted for pair-end sequencing through an Illumina high-throughput sequencing platform (Illumina, Inc. San Diego, CA, USA).

SLAF-seq data grouping and genotype definition

In sequencing data assay, P, M, aa and ab value represented resistant parent HZ, susceptible parent LTH, resistant pool (Rp) and susceptible pool (Sp), respectively. All SLAF single-end reads were grouped based on sequence similarity as detected by BLAT (-tile size = 11, -step size = 11).[44] Sequences with 90% or more similarity were clustered into one SLAF locus. Alleles were specified to each SLAF based on the minor allele frequency evaluation. Four kinds of SLAF tags are in one locus, since rice is a diploid species. Repetitive SLAFs were filtered out. SLAFs with sequence depth less than 213 were considered as low-depth SLAFs, which were unsuited for further analysis. SLAFs with two, three or four tags were considered to polymorphic SLAFs and regarded as potential markers for further polymorphic analysis. The potential markers were classified into three segregation patterns. We confirmed the origin of alleles according to sequencing depth of the parents and R/S pools. The potential markers with one genotype deriving from P and the other from M were verified as polymorphic markers.

Resistance gene association analysis

SLAF markers with Ratio_ab ≥ 2 were considered to be Diff_markers. Regions with three or more consecutive Diff_markers were regarded as candidate resistance gene(s) regions. Difference ratio was counted as follows: Paa represented depth of Rp (aa) population derived from P, and Maa indicated depth of Rp (aa) population derived from M; Pab represented depth of Sp (ab) population derived from P, and Mab referred to depth of Sp (ab) population derived from M. The associated value was calculated as follows: Ratio_aa = Paa/Maa; Ratio_ab = Mab/Pab. The Diff_markers’ statistics procedure was as follows: two kinds of genotypes derived from parents were selected according to suitable depth (60–300) after standardizing the parents’ markers data. Diff_markers were screened according to the difference values of Rp and Sp. The gene mapping candidate regions were ascertained based on the distribution and difference ratio of Diff_markers.

Resistance spectrum analysis of HZ and eight near isogenic lines (NILs)

A set of monogenic lines carrying different major single blast R genes, kindly provided by IRRI, including four broad spectrum major genes Pi2,Pii,Pi5,Pi7 and four blast resistance genes Pita,Pita-2,Pi12,Pi19 located on chromosome 12, were used as the monogenic lines. Fifty M. oryzae representative isolates collected from South China were used to evaluate the blast resistance spectrum (RS) of restorer line HZ carrying Pi-hz(t) gene and the monogenic lines mentioned above (Table S2). RS was calculated as followed fomula: RS = (number of incompatible isolates/number of total tested isolates) × 100%.

Results and discussion

Genetic analysis of resistance gene and construction of gene pools

The resistant donor HZ was crossed with susceptible parent LTH. Their F₁ individuals all showed resistant phenotype against M. oryzae isolate GD10-3001. An F₂ population derived from HZ/LTH was used to analyse resistance inheritance of HZ against GD10-3001. Segregation of resistant (R) and susceptible (S) progenies of 219 F₂ individuals fitted a 3:1 ratio (171 R: 48 S, x² = 0.42, P = 0.52). The segregated result suggested that HZ harbours a dominant resistance gene designated as Pi-hz(t) tentatively. The plants with disease scales 0–2 were considered to be distinctly resistant, while the plants with disease scales between 6 and 9 were recognized to be distinctly susceptible. In this study, two parents HZ and LTH, 30 distinctly resistant and 30 clearly susceptible F₂ individuals DNA were extracted and mixed as Rp and Sp for further SLAF-seq and super-BSA analysis.

SLAF-sequencing and SLAF tags data analysis

Total DNAs of four samples HZ, LTH, Rp and Sp were digested into fragments for SLAF library construction and SLAF-sequencing. After SLAF library construction and high-throughput sequencing, a total of 3,670,254,150 Mb data including 18,169,575 single-end reads were procured with each read length of ∼202 bp originally. Among them, only 1.69 Gb data including 6,190,100 valid single-end reads, which were constructed from the samples of HZ (1,121,190), LTH (1,704,652), Rp (2,050,068) and Sp (1,314,190), respectively (Table 1). The number of SLAF tags was 56,187, and the average coverage for each tag was 108.79-fold. All the SLAF tags were drew evenly in the distribution diagram on each rice chromosome (Figure 1), which illustrated that the rice genome has been successfully simplified, since the SLAF tags were positioned on the rice genome equally.

Table 1. Sequencing data of each sample.

Sample	Sample ID	Useful reads
HZ	P	1,121,190
LTH	M	1,704,652
Rp	Aa	2,050,068
Sp	Ab	1,314,190
Total		6,190,100

Figure 1. SLAF tag distribution on each rice chromosome.

SLAF tags polymorphic analysis

The 56,187 high-quality SLAF tags were divided into six kinds of markers. Among the SLAF markers, 9051 polymorphic markers were identified based on polymorphic analysis of allele numbers and variation among gene sequences. Among the 56,187 SLAF tags, the polymorphic tags markers were divided into three kinds: INDEL (insertion or deletion), SNP and EDP (endonuclease digestion polymorphism), of which their polymorphism rates were 0.24%, 15.62% and 0.23%, respectively. The other useless types were: no polymorphism (81.29%), unknown type (2.11%) and repeat sequence (0.48%), etc. The SLAF tags were distributed explicitly and evenly on each rice chromosome (Figure 2), which proved that the SLAF-seq was prosperous.

Figure 2. Markers distribution on each rice chromosome.

Association analysis and gene annotation of the target region

The genetic analysis of the F₂ population indicated that HZ contained at least a dominant resistance gene. The F₂ population of HZ/LTH was divided into distinctly resistant and susceptible pools with 30 individuals each. We simplified the genetic model for association analysis. A total of 5648 markers with one genotype deriving from P and the other from M were selected to compute the difference ratio of Rp and Sp. 91 Diff_markers were obtained and mapped on chromosome 12 according to association analysis. The Diff_markers distributed on chromosome 12 were visualized in Figure 3 statistically. The associated assay indicated that target gene was mapped on 10M–20M on chromosome 12, since this region had clear concentration of Diff_markers based on the Ratio_ab above green lines and Ratio_aa under purple lines (Figure 3). Further signal amplification of this region narrowed the target gene in 10M–18M interval (Figure 4), since the Diff_markers were mapped on this region. The intervening genomic region (about 7.4 Mb) of chromosome 12, where the target gene was mapped, contains at least 143 genes annotated according to the International Rice Genome Sequencing Project (IRGSP) build 5.0 (Table 2; supplementary Table S1; Figure 4).

Figure 3. Ratio distribution of Diff_markers on chromosome 12. The Ratio_ab above the green line and Ratio_aa under the purple line indicated that corresponding markers are linked with target gene. There was a obvious concentrated region on 10M–20M based on the Diff_marker above the green and under the purple line.

Figure 4. Ratio distribution of Diff_markers on 10M–18M of chromosome 12. The 91 Diff_markers on chromosome 12 were selected to carry out association analysis. The regions with three or more consecutive Diff_markers were considered to be related candidate gene regions. Nine related candidate gene regions overlapped 10.2M–17.6M on chromosome 12 were obtained and there were 42 Diff_markers and 143 annotated candidate genes in which based on the conference genome sequence version IRGSP build 5.0 (Table 2). Taken together, the candidate gene regions of the Pi-hz(t) gene covered a size of 7.4 Mb, approximately.

Table 2. Information of the association region.

ID	Chr_ID	Start	End	Size (M)	Marker numbers	Gene
1	Chromosome 12	10,214,000	10,231,000	0.017M	3	2
2	Chromosome 12	10,709,000	10,982,000	0.273M	3	27
3	Chromosome 12	12,680,000	12,797,000	0.117M	4	4
4	Chromosome 12	13,183,000	13,632,000	0.449M	3	23
5	Chromosome 12	13,785,000	14,227,000	0.442M	6	24
6	Chromosome 12	14,341,000	14,539,000	0.198M	3	3
7	Chromosome 12	14,648,000	14,912,000	0.264M	5	6
8	Chromosome 12	15,312,000	15,494,000	0.182M	3	3
9	Chromosome 12	16,071,000	17,629,000	1.558M	12	51

Note: Start and end represent position on Nipponbare reference genome; Diff_marker represents the polymorphic SLAF tag with a difference ratio (Ratio_ab ≥ 2).

Resistance characterization of HZ and eight NILs

Resistance evaluation of HZ and eight NILs against 50 M. oryzae isolates collected from South China indicated that HZ conferred a high RS with 94% of tested isolates, the NIL IRBLz5-CA had higher RS with 90%, two NILs IRBL7-M and IRBLta2-Pi had middle RS with 58% and 46%, and five NILs IRBL5-M, IRBLi-F5, IRBLta-CP1, IRBL12-M and IRBL19-A had lower RS with 28%, 26%, 12%, 6% and 6%, respectively (Table 3; Table S2). The resistance characterization of HZ and eight NILs against M. oryzae isolates derived from South China suggested that HZ can be used as an excellent resistant donor for MAS breeding in South China.

Table 3. Resistance spectrum of HZ and eight near isogenic lines against 50 M. oryzae isolates from South China.

Lines	Resistance spectrum (%)
HZ	94
IRBLz5-CA (Pi2)	90
IRBLi-F5 (Pii)	26
IRBL5-M (Pi5)	28
IRBL7-M (Pi7)	58
IRBLta-CP1 (Pita)	12
IRBL12-M (Pi12)	6
IRBL19-A (Pi19)	6
IRBLta2-Pi (Pita-2)	46

Conclusions

In this study, 56,187 of SLAF tags had been developed with SLAF method. Nine thousand and fifty-one polymorphic markers were identified with high quality and quantity that met the requirements. The integral and accurate SLAF tags were well-distributed on each rice chromosome. The 91 Diff_markers provided useful data for precisely narrowing down the candidate trait-associated regions into a 7.4 Mb interval which had 143 genes annotated by IRGSP build 5.0. Annotation of the 143 candidate genes showed that there are 10 NBS-LRR (nucleotide-binding site-leucine-rich repeat) coding genes as resistant candidates in the target region, which could contribute to isolate the Pi-hz(t) gene in further research. BSA combined with SLAF-seq method could be an efficient and cost-effective way for resistance genes or QTLs mapping. Resistance characterization analysis of HZ and eight NILs showed that HZ and IRBLz5-CA had higher RS with 94% and 90%, which suggested that HZ and IRBLz5-CA would be excellent resistance source for improving blast resistance breeding program in South China. The lower level RS of NILs IRBL5-M, IRBLi-F5, IRBLta-CP1, IRBL12-M and IRBL19-A indicated that resistance genes Pi5, Pii, Pita, Pi12 and Pi19 are not recommended to apply in South China.

Acknowledgments

The authors greatly appreciate the kind help of Biomarker Technologies Company Limited in providing the technology support.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This research is supported by grants from the Guangdong Science and Technology Projects [project number 2013B020301011]; National Natural Science Foundation [grant number 31101403], [grant number 31301304]; Guangzhou Sciences and Technology Project [project number 2012J2200066]; Earmarked Fund for Modern Agro-Industry Technology Research System [CARS-01-24]; National Transgenic Research Projects [project number 2014ZX0800904B]; Science and Technology Project of Guangdong, China (the Construction of Green Management R&D Center for Crop Pests in Guangdong Province).