RNA‐seq bulked segregant analysis combined with KASP genotyping rapidly identified PmCH7087 as responsible for powdery mildew resistance in wheat

Powdery mildew causes considerable yield losses in common wheat (Triticum aestivum L.) production. Mapping and cloning powdery mildew‐resistant quantitative trait loci can benefit stable yield production by facilitating the breeding of resistant varieties. In this study, we used the powdery mildew resistance introgression line ‘CH7087’ (harboring the resistance gene PmCH7087) and developed a large F2 population and a corresponding F2:3 segregation population containing 2,000 family lines for molecular mapping of PmCH7087. Genetic analysis demonstrated that the resistance phenotype was controlled by a single dominant gene. According to the performance exhibited by the F2:3 lines, 50 resistant lines and 50 susceptible lines without phenotype segregation were chosen for pooling and bulked segregant RNA sequencing (BSR‐Seq) analysis. A region spanning 42.77 Mb was identified, and genotyping of an additional 183 F2:3 lines with extreme phenotypes using 20 kompetitive allele‐specific polymerase chain reaction (KASP) markers in the BSR‐Seq mapping regions confirmed this region and narrowed it to 9.68 Mb, in which 45 genes were identified and annotated. Five of these transcripts harbored nonsynonymous single nucleotide polymorphisms between the two parents, with the transcripts of TraesCS2B01G302800 being involved in signal transduction. Furthermore, TraesCS2B01G302800.2 was annotated as the closest homologue of serine/threonine‐protein kinase PBS1, a typical participant in the plant disease immune response, indicating that TraesCS2B01G302800 was the candidate gene of PmCH7087. Our results may facilitate future research attempting to improve powdery mildew resistance in wheat and to identify candidate genes for further verification and gene cloning.


INTRODUCTION
Common wheat (Triticum aestivum L) is one of the major crops providing food energy to humans, and its stable reproduction is threatened by various diseases (Wellings, 2011). Powdery mildew, one of the most devastating diseases affecting wheat, is caused by Blumeria graminis f. sp. tritici and can cause not only dramatic yield losses but also decreased wheat grain quality in common wheat worldwide (Singh et al., 2016). Identifying and utilizing new genes to breed cultivars with improved disease resistance have long been a crucial and effective strategy formed for overcoming these biotic stresses (Nelson et al., 2018). Marker-assisted selection can accelerate these breeding programs (Anderson, 2009;Yaniv et al., 2015). Increasing numbers of powdery mildew resistance genes and quantitative trait loci (QTLs) have been reported in the gene pool of elite wheat varieties and its wild relatives (Li, Dong, Ma, et al., 2020;Lu et al., 2020;Qie et al., 2019;Xie et al., 2020). However, major resistance genes generally are overcome rapidly when they are applied over larger areas and for longer periods of time because of changes in virulence and the appearance of novel races in the pathogen population. Therefore, it is necessary to investigate novel resistance gene pools and use long-lasting disease resistance genes to breed resistant wheat cultivars.

Core Ideas
• A span of 42.77 Mb and 385 genes on chromosome 2B were annotated. • Based on the KASP genotypes, seven SNP markers were associated with PmCH7087. • TraesCS2B01G302800 was the most likely candidate gene for PmCH7087.
reported two cloned genes Pm2 and Pm4, the latter of which encoded a putative chimeric protein of a serine/threonine kinase and multiple C2 domains and transmembrane regions (Sánchez-Martín et al., 2016;Sánchez-Martín et al., 2021). Pm21, encoding a typical coiled-coil nucleotide binding siteleucine-rich repeat protein, confers broad-spectrum resistance to powdery mildew disease in wheat (He et al., 2018). Through homology-based cloning, the Pm8 and Pm17 were isolated on the 1RS introgression in wheat cultivar 'Amigo' (Singh et al., 2018). Through a map-based strategy, Pm24 was cloned; this gene encodes the WTK3 protein with putative kinase-pseudokinase domains and is a broad-spectrum powdery mildew resistance gene . Pm41, located on chromosome 3BL, encodes the CNL protein and confers powdery mildew resistance Wang et al., 2014). Bulked segregant RNA sequencing (BSR-Seq) is a modified bulked segregant analysis for mapping genes and was first developed and applied in maize (Zea mays L.), a species with a relatively larger genome size (Liu et al., 2012), and subsequently introduced to wheat as described above. Segregant RNA sequencing used RNA-seq data to call single nucleotide polymorphisms (SNPs), which balances the sequencing cost and molecular marker density for species with a large genome. In addition, expression information, especially regarding differentially expressed genes (DEGs), can provide the necessary information for gene screening, especially for adverse stress-related traits, such as disease resistance. The mapping resolution of BSR-Seq has also been determined to be satisfactory (Liu et al., 2012). Recently, many genes from different crops that are responsible for disease resistance, abiotic stress, and development have been reported. Pm4b was reported to be located on chromosome 2AL by BSR-Seq and further mapped in a 3.0-cM genetic interval spanning 6.7 Mb (Wu et al., 2018). A two-step BSR-Seq approach identified Pm5e, which encodes a leucine-rich-repeat-containing protein, and a rare nonsynonymous single nucleotide variant in which the C-terminal leucine-rich repeat domain was determined to confer the resistance to powdery mildew . Yr15, a wheat yellow rust disease resistance gene, was mapped to a 0.77-cM interval by implementing BSR-Seq and its corresponding SNP markers were developed (Ramirez-Gonzalez et al., 2015). Through a combination of BSR-Seq and linkage analysis, Pm61 was located in a 0.71-cM genetic interval in common wheat . In the Chinese cabbage (Brassica rapa L. ssp. pekinesnis) DH line 'FT', the male sterile gene ftms was mapped onto chromosome A05 via BSR-Seq analysis, and the region was narrowed to approximately 1.7 Mb via linkage analysis (Tan, Liu, et al., 2019). To isolate waterlogging-tolerance genes in maize, however, a set of 10 waterlogging sensitive and eight waterlogging-tolerant inbred lines were adopted for BSR-Seq, and six alleles were identified by integrating the DEG data and genomic mapping positions (Du et al., 2017).
Tall  (Li & Wang, 2009). The resistance genes can be introgressed into wheat cultivars by wheat-T. ponticum distant hybridization. The introgression CH7087 has demonstrated high resistance to powdery mildew in wheat. In the present study, a large F 2 population was developed by crossing CH7087 with 'Taichung 29', a powdery mildew-susceptible cultivar. A subsequent F 2:3 family population with 2,000 lines was also developed. After phenotyping, BSR-Seq of bulked pools was used to identify SNPs associated with powdery mildew resistance. Subsequently, we developed kompetitive allele-specific polymerase chain reaction (KASP) markers to genotype the extreme F 2:3 lines to further confirm and narrow down our mapping region. Nonsynonymous SNP annotations, gene enrichment analyses, and DEGs were combined to screen thecandidate gene. The markers linked to PmCH7087 may be useful for marker-assisted selection in wheat breeding programs and map-based cloning of PmCH7087.

Materials
The powdery mildew-resistant accession 'R431', derived from T. ponticum and the common wheat cultivar Taichung 29, was provided by Shanxi Agricultural University, Taiyuan. CH7087 is an introgression line with high powdery mildew resistance, whereas Taichung 29 is highly susceptible to powdery mildew disease. To develop high-resolution molecular mapping of the resistant QTLs underlying CH7087, an F 2 segregation population with 2,372 individuals was developed by crossing CH7087 with Taichung 29 (equivalent to a backcross). The corresponding F 2:3 families were also developed for powdery mildew resistance phenotyping, and the performance of every F 2:3 family individual was used to evaluate the corresponding F 2 resistance.

Powdery mildew resistance testing
Powdery mildew resistance was evaluated at the seedling stage. The CH7087, Taichung 29, 'SY95-71', F 1 hybrid and F 2:3 families were grown in a greenhouse with a light time of 18 h and a dark time of 6 h at 26°C. The resistance of CH7087, Taichung 29, SY95-71, and part of the F 1 hybrid populations were evaluated in 2016. The parents and F 2 populations were evaluated in 2017. Four seeds were planted discretely in a 4by 4-by 4-cm plug tray. Twenty seeds from all three parents (CH7087, Taichung 29 and SY95-71) and the F 1 hybrid were sown into mixed peat with vermiculite (2:1). Approximately 15 seeds were sown for every F 2:3 family line. At the oneleaf stage, the spores of B. graminis f. sp. tritici isolate E09 were inoculated artificially and water was sprayed to maintain wetness. To ensure successful infection, we performed inoculations two times in the morning on two continuous days. The wheat cultivar Mianyang 11 was used as a susceptible control. After the susceptible control was heavily infected (approximately 7 d after inoculation), the infection types (ITs) were evaluated on a 0 to 4 infection type scale according to a previously described method (He et al., 2009). The appearance of F 2:3 lines was recorded three times to ensure the reliability of the disease resistance measurements. For seed reproduction of the F 2 and F 2:3 families, sowing, watering, and fertilizer management were the same as in field production. In order to understand the response to different strains, the performance of CH7087 were evaluated as well in the light incubator with different B. graminis f. sp. tritici strain (E20, E21, and E26). The evaluation conditions and methods were the same as above.

RNA library construction and Illumina sequencing
Individuals of one family line with a consistent phenotype in the F 2:3 population were chosen to guarantee accurate IT pooling. Fifty homozygous resistant lines (IT = 0) and 50 homozygous susceptible lines (IT = 4) of the F 2:3 families were screened. Leaves of 10 randomly selected F 2:3 lines after inoculation for 7 d were collected and mixed equally to extract the total RNA, which was treated as the corresponding genotypes of F 2 . The leaves of the parents CH7087 and Taichung 29 after inoculation for 7 d were also sampled. Total RNA was extracted by the total RNA extractor (Sangon Biotech). Equivalent RNA of every single F 2:3 line was pooled to construct resistant and susceptible bulks. Finally, RNA samples were divided into four groups: the resistant parent CH7087 (R01), the susceptible parent Taichung 29 (S01), homozygous resistant bulks (R02), and homozygous susceptible bulks (S02) for subsequent sequencing library construction. After measurement of the concentration with a Nanodrop Drop 2000 spectrophotometer (Thermo Scientific) and quality control with an Agilent Bioanalzer 2100 (Agilent Technologies), cDNA libraries were generated. All the libraries were loaded into one lane of a flow cell and sequenced by an Illumina X-TEN platform with the PE150 model according to the manufacturer's instructions.

Differentially expressed gene analysis
Clean data of the four groups were obtained by filtering sequencing reads with more than 50% of bases having a Q value of ≤10 or an ambiguous sequence content exceeding 10%. Next, the remaining reads were aligned to the wheat genome sequence (https://urgi.versailles.inra.fr/download/ iwgsc/IWGSC_RefSeq_Assemblies/v1.0/, accessed 28 May 2021) by TopHat2 software with the default settings (Kim et al., 2013). The components of Cuffquant and Cuffnorm were used to determine the quantitative expression as fragments per kilobase of transcript per million fragments mapped (Trapnell et al., 2010). DESeq (Anders & Huber, 2010) was used for calling DEGs at the level of fold change ≥2 and a false discovery rate of <0.01. Further detailed information for typical RNA-seq analysis was obtained according to previously described methods with slightly modified points (Xue et al., 2017).

Single nucleotide polymorphism calling and molecular mapping
After TopHat2 mapping was performed, the GATK Hap-lotypeCaller module was used for SNP calling in parallel (McKenna et al., 2020) and SnpEff was used for SNP annotation (Cingolani et al., 2012). Single nucleotide polymorphisms with a sequencing depth of <4 in all four groups were filtered. The SNP-index algorithm (Takagi et al., 2013) was used to estimate the association between high-quality SNPs and powdery mildew. SNPNUM was used to fit the SNPindex values . Powdery mildew resistanceassociated loci were identified when the fitted values of the SNP-index were above the threshold at the 99% confidence interval . A circos graph containing chromosomes, genes, and SNP density was generated by CIRCOS 0.66 software (http://circos.ca/, accessed 28 May 2021).

2.6
Bulked segregant RNA sequencing mapping confirmation and candidate gene screening Bulked segregant RNA sequencing analysis resulted in an region associated with powdery mildew resistance. To con-firm and narrow the region, we used KASP markers in an additional panel of individuals with extreme phenotypes from the F 2:3 population. We selected 30 SNPs evenly distributed throughout the region and converted them into KASP markers for population genotyping. The primer information is shown in Supplemental Table S1. In total, 183 BSR individual-free F 2:3 families, including 136 resistant families and 47 susceptible families, were genotyped by KASP markers, and a further t-test was used to generate new regions harboring significantly different SNPs. Moreover, nonsynonymous SNPs and DEGs between the two parents (R01 vs. S01) and between the resistant and susceptible bulks (R02 vs. S02) were identified in the new region. Gene Ontology (Ashburner et al., 2000), Kyoto Encyclopedia of Genes and Genomes (Kanehisa & Goto, 2000), eggNOG (Huerta-Cepas et al., 2016), and NR (Deng et al., 2006) database enrichment analyses were conducted to help identify potential candidate genes.

Genetic analyses of the resistant phenotype
The performance of CH7087, Taichung 29, and the F 1 hybrid and F 2:3 families inoculated with B. graminis f. sp. tritici isolate E09 was evaluated. To simplify phenotyping, plants with scores between 0 and 2 points were regarded as resistant; the rest were marked as susceptible. Typical performance scores of 0 to 4 are shown in Supplemental Figure S1. As shown in Table 1, CH7087 and the F 1 population were highly resistant (IT = 0), whereas Taichung 29 and SY95-71 were highly susceptible (IT = 4). The data from F 2:3 families segregated as 533 homozygous resistant: 1,102 segregating: 520 homozygous susceptible (χ 2 1:2:1 = 1.271, P df2 = .530) ( Table 1). The segregation ratio was in keeping with previous results that powdery mildew resistance was indicative of dominant inheritance in CH7087. The resistant gene was temporarily designated as PmCH7087.

Bulked segregant RNA sequencing analyses
The BSR-Seq strategy was used to map the resistance gene PmCH7087. RNA sequencing produced 12. 74, 13.24, 40.54, and 40.02 Gb of clean data in R01, S01, R02, and S02, respectively. The percentage of the Q30 value (the Q30 valued s is equivalent to the probability of an incorrect base call 1 in 1,000 times) ranged from 93.04 to 93.61%, and the GC content of each sample ranged between 56.55 and 58.08%. In addition, the mapped ratio was between 84.27% and 86.67% for the four samples (Table 2). These results indicated that a relatively reliable dataset was generated for T A B L E 1 The segregation patterns of F 1 plants, F 2 : 3 families, and their parents inoculated with powdery mildew The percentage of bases for which the quality value of the clean data was ≥30; Mapped reads: the number of clean reads mapped on the reference genome; Mapped ratio: the percentage of mapped reads in the clean reads.

No. of plants or lines
follow-up analysis. The clean sequencing data were uploaded to the Sequence Read Archive database (accessions: SRR10019739, SRR10019740, SRR10019741, and SRR10019742). We further identified 64,999 expressed genes from the four sequencing libraries and 64,709 SNPs were called between the parent CH7087 and Taichung 29 in total ( Figure 1, Table 3). After filtering, we obtained 10,187 high-quality SNPs. The SNPs were roughly enriched on the long arms and short arms of the chromosome (Figure 1). Table 3 shows the information of these SNPs. The average SNP density was 0.68 SNPs per Mb, with the highest density being in chromosome 3B (1.10 SNPs per Mb) and the lowest density being in chromosome 4A (0.47 SNPs per Mb). The SNP sequencing depth of the parents (R01 and S01) fluctuated between 7.75× to 9.49×, and the two pools ranged from 20.44× to 22.94×; 2,608 of the 10,187 SNPs were nonsynonymous (Table 3).
To map PmCH7087, individuals with extreme phenotypes were pooled for BSR-Seq analysis. Two regions on chromosome 2B (425,795,642-453,022,545; 666,652,840-682,194,230) were estimated to be associated with powdery mildew disease (Figure 2). These regions spanned 42.77 Mb, and 385 genes were annotated in the two regions by the reference genome (Chinese Spring reference genome V1 version), including 25 genes harboring nonsynonymous SNPs between the two parents.

Confirmation and further mapping of PmCH7087
Since many genes were annotated in the BSR-Seq mapping region, 30 KASP markers spanning from 425,795,642 to 682,194,230 bp on chromosome 2B were used to genotype additional F 2 individuals for confirmation and further mapping of PmCH7087. Out of 30 SNPs, 20 were genotyped successfully in the two parents, which was in keeing with the SNPs called by RNA sequencing, and were further used for genotyping.
On the basis of the the KASP genotypes between the contrasting phenotyping individuals, a t-test demonstrated that seven continuous SNP markers (from BSNP96 to BSNP461) showed significant differences (Figure 3), indicating that these SNP markers were associated with PmCH7087. This region was between 425,795,642 and 435,479,154, a smaller region of 9.68 Mb, in which 45 transcripts were annotated; five of these transcripts harbored nonsynonymous SNPs between the two parents (Supplemental Table  S2). Furthermore, the expression status of this 9.68 Mb F I G U R E 1 Circos plot of genome-wide genes and single nucleotide polymorphisms (SNPs) distributions. The outer circle indicates chromosomes; the middle circle indicates gene distribution, and the inner circle indicates the SNP density distribution region in the Chinese Spring reference genome was investigated. The DEG analysis demonstrated that no gene was upregulated or downregulated between R01 and S01 and between R02 and S02. Gene enrichment analyses showed that TraesCS2B01G302800.2 and TraesCS2B01G302800.1 were annotated in signal transduction mechanisms by the eggNOG database. Moreover, TraesCS2B01G302800.2 was annotated as a serine/threonine-protein kinase PBS1 by the NR database (Supplemental Table S2). These results indicated that TraesCS2B01G302800 was the candidate gene regarding PmCH7087, but this possibility warrants further confirmation and verification of the deep molecular functional.

PmCH7087 is a single dominant gene for wheat powdery mildew resistance
In the present study, we investigated the genetic characteristics of PmCH7087 and found that all F 1 plants were diseaseresistant, and a typical Mendelian 3:1 segregation ratio of the F 2 population demonstrated that the resistance of powdery mildew was controlled by a single dominant gene (Table 1). More importantly, PmCH7087 was mapped to chromosome 2B by BSR-Seq (Figure 2), which was different from the locus reported earlier (Bin 2BL: 0.89-1.00) (Zhan et al., 2014), indicating different underlying powdery mildew resistance mechanisms in T. ponticum.
Several previously identified QTLs have been mapped onto the chromosome 2B, such as Pm6, Pm33, Pm51 (Zhan et al., 2014), Pm52, Pm63, Pm64, MlZec1, PmJM22 (Yin et al., 2009), and MlAB10 (Maxwell et al., 2014. To investigate the relationship between PmCH7087 and these genes, we mapped the corresponding linked markers of these QTLs onto the Chinese Spring reference genome V1 (Supplemental Table S3). No shared region was observed between PmCH7087 and these QTLs except for Pm33, Pm52, and MlZec1, which had been mapped previously. In addition, different disease resistance patterns exist between PmCH7087 and other gene series. The source, resistance phenotype, and position of the resistance genes indicate that PmCH7087 may represent a novel resistance gene resource.

4.2
TraesCS2B01G302800 is a putative candidate gene for PmCH7087 Since 385 genes were annotated in the primary mapping region, we further developed 20 SNP-based KASP markers that were evenly distributed throughout the two regions, and we genotyped 183 individuals (136 resistant individuals and 47 susceptible individuals) of the F 2 population with extreme performance. The t-test resulted in seven KASP markers associated with disease resistance, which covered 9.68 Mb (425,795,479,154 bp). No associated SNP was detected in the region from 666,652,840 to 682,194,230 bp identified by BSR-Seq, suggesting that these two regions were false positives. The confirmed region (425,795,642-435,479,154 bp)  The protein kinase AvrPphB susceptible1 (PBS1) has long been considered to play an important role in the plant immune response (Kim et al., 2016;Pottinger et al., 2020;Qi et al., 2014). We screened a putative candidate gene, TraesCS2B01G302800, with a nonsynonymous SNP. This gene was annotated as encoding the serine/threonine-protein kinase PBS1. However, the expression of this gene exhibited no significant differences between the two parents or between the two extreme pools. We also compared the coding sequences of TraesCS2B01G302800 among different plant species. High sequence similarity was observed in two durum wheat (Triticum durum Desf.) varieties ('Svevo' and 'Zavitan') and nine common wheat varieties ('Arinalrfor', 'Jagger', 'Julius', 'Lancer', 'Landmark', 'Mace', 'Norin61', 'Stanley', and 'SyMattis') (Supplemental Figure S1), indicating that the homologs were conservative and might endure strong selection pressure. To further verify the function of TraesCS2B01G302800, the mutant libraries developed by UC Davis (Krasileva et al., 2017) and Yantai University (10th National Congress on Wheat Genomics and Molecular Breeding, unpublished data, 2019) can be requested to screen candidate gene-mutated mutants and the performance of these mutants under powdery mildew stress. Protocols in molecular biology, such as transient transfection (virus-induced gene silencing) or hereditable transformation (transgenosis and gene editing), can provide solid evidence to uncover the mystery of TraesCS2B01G302800.
Map-based cloning is considered to be an effective technique to identify and isolate target genes. Pm3b (Yahiaoui et al., 2004), Pm4 (Sánchez-Martín et al., 2016), Pm21 (Ying et al., 2018, Pm24 , and Pm41  were reported to have been cloned by this strategy. In light of our identification of the resistance gene PmCH7087 in this study, further backcrossing to develop a backcross population or a near-isogenic line might be necessary to map PmCH7087 into a smaller region, which could strengthen the mapping results obtained and isolate the gene. In that case, PmCH7087 can be widely used to enhance powdery mildew resistance in wheat breeding.

Segregant RNA sequencing combined with KASP in a large F 2:3 family can accelerate gene mapping for wheat
With development of high-throughput sequencing technologies, BSR-Seq has become an effective method for rapid discovery of novel genes and genetic markers linked to target genes (Pearce et al., 2014). Although BSR-seq has many advantages in the study of target genes, there are still some limitations. For example, the SNPs from transcriptome sequencing data are only from exon sequences, whereas the majority of polymorphisms reside in intergenic regions. A combination analysis of BSR-Seq and resequencing of the two parents might deal with this deficiency. Furthermore, there were strict requirements for the time and method of sampling in BSR-Seq. Despite all these, BSR-seq still provides a large number of effective markers for gene mapping and isolation (Nishijima et al., 2017).
Recently, a protocol combining bulked segregant analysis (pooling DNA) and local QTL mapping via KASP genotyping was developed to rapidly map the Cf-10 gene to a 790-kb region in an F 2 population . In the present study, we used BSR-Seq for primary mapping and KASP markers to further map the powdery mildew gene PmCH7087 in a large F 2:3 population. As disease resistance always gives rise to transcription factor responses, BSR can not only the map the function but can also help us screen DEGs and transcription factors in mapping regions, especially for regions with limited combinations. Similarly, TmBr1, a gene responsible for stem brittleness, was isolated by the MapRseq approach and was present in recombinant cold spots (Deng et al., 2019). Therefore, we believe that combining BSR-Seq analysis and KASP in the primary mapping region in a large F 2:3 population can considerably accelerate the isolation of the genes responsible for stress-related traits, such as disease resistance and resistance to abiotic stresses (drought stress, temperature stress, or mineral deficiency stress).

A C K N O W L E D G M E N T S
This work was supported by the National Natural Science Foundation of China (No. 31601303) and the Shanxi Provincial Program of key R&D projects (China) (No. 201903D221053).

C O N F L I C T O F I N T E R E S T
The authors declare that there is no conflict of interest.