The wheat breeding line LS5082 was bred and provided by Prof. Sishen Li, Shandong Agricultural University, China, and maintained in our lab as a resistant germplasm to powdery mildew. The wheat cultivars Shannong 29, Shimai 22 and Huixianhong, which are susceptible to all the Bgt isolates tested, were used as susceptible parents and crossed with LS5082 to obtain F1 hybrids, F2 populations and F2:3 families for genetic analysis of the powdery mildew resistance in LS5082. Wheat cultivar Huixianhong was used as the susceptible control for phenotypic assessment. Five resistant donors, Coker 747 (with Pm6) (Bennett et al. 1984), Liangxing 99 (with Pm52) (Wu et al. 2019), CH7086 (with Pm51) (Zhan et al. 2014), WE35 (with Pm64) (Zhang et al. 2019), Qingxinmai (with PmQ) (Li et al. 2020b) were used to compare their phenotypic responses to different Bgt isolates with those of LS5082 through multi-race response experiments. Twenty-six susceptible cultivars collected from different ecological regions of China were used to evaluate the applicability of closely linked markers for MAS.
Reactions to different Bgt isolates and genetic analysis
To evaluate the resistant spectrum of the Pm gene in LS5082 and compare its reaction pattern with the catalogued Pm genes on the same chromosome arm, LS5082 was tested against 11 Bgt isolates collected from different wheat production regions of China and with different virulence patterns using genotypes Coker 747 (with Pm6), Liangxing 99 (with Pm52), CH7086 (with Pm51), WE35 (with Pm64), Qingxinmai (with PmQ) as controls. For each Bgt isolate, five seeds of each of the genotypes above were sown in 72-cell rectangular trays with each tray 3 × 3cm. The susceptible control Huixianhong was sown randomly in each tray. Each tray was then put in an independent space with high humidity to avoid cross-infection. The growing conditions were set at a daily cycle of 14 h light at 24oC and 10 h of darkness at 18oC. When the seedlings had grown to the one-leaf stage, they were inoculated with fresh conidiospores, previously multiplied on Huixianhong seedlings. For the first 24h after inoculation, the seedlings were put in a dark place with high humidity at 18oC. From the second day, the growing conditions were set at a daily cycle of 14 h light at 22oC and 10h of darkness at 18oC. When the spores were fully developed on the first leaves of the susceptible control, at 10-14 days post inoculation, infection types (ITs) were surveyed using the 0-4 scale described by An et al. (2013), in which ITs 0, 0;, 1 and 2 are regarded as resistant, and ITs 3 and 4 as susceptible. For all the experiments above, three independent repeats were carried out.
To investigate the inheritance of the powdery mildew resistance in LS5082, the Bgt isolate E09 was firstly selected to inoculate one-leaf seedlings of LS5082, Shannong 29, Shimai 22 and Huixianhong, and their crossed F1 plants, F2 population and F2:3 families for genetic analysis. For the resistant and susceptible parents and their F1 hybrids, 10 seeds were selected for sowing; for each F2 population, no fewer than 200 seeds were sown; and for each F2:3 family, 20 seeds were sown. As the segregation ratio against isolate E09 was consistent with that of monogenetic inheritance, other isolates avirulent on LS5082 were all used to inoculate randomly selected 10 homozygous resistant, 10 segregating and 10 homozygous susceptible F2:3 families for genetic analysis using the same procedures as above to clarify if the same gene conferred resistance to all the Bgt isolates, i.e., to confirm that no other genes were included in LS5082 conferring resistance to powdery mildew. When the phenotypic data was obtained, goodness-of-fit analysis was carried out using a chi-squared (χ2) test to assess deviations of the observed phenotypic data from theoretically expected segregation ratios (SPSS 16.0 software, SPSS Inc, Chicago, USA) at P ≤ 0.05.
Preparation of samples for BSR-Seq, library construction and RNA sequencing
After scoring the ITs of F2:3 families of LS5082 × Shannong 29 using Bgt isolate E09, 50 resistant plants from 50 homozygous resistant F2:3 families, and 50 susceptible plants from 50 homozygous susceptible F2:3 families were selected for isolation of their total RNA using the Spectrum Plant Total RNA kit (Sigma-Aldrich) following the manufacturer’s protocol. Equal amounts of RNA from each of the 50 resistant plants were mixed, and similarly for the 50 susceptible plants to construct resistant and susceptible RNA bulks, respectively. The RNA integrity of the resistant and susceptible bulks was assessed using the Agilent 2100 Bio analyzer (Agilent Technologies, Santa Clara, CA, USA). An RNA Integrity Number (RIN) ≥ 7 was considered to meet the sequencing standard. The cDNA libraries were constructed using TruSeq Stranded mRNA LTSample Prep Kit (Illumina, San Diego, CA, USA) following the manufacturer’s protocol, and the quality of the cDNA libraries was assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The eligible cDNA libraries were sequenced on the Illumina HiSeq sequencing platform (Illumina HiSeq4000), sequenced according to the standard protocol by Beijing Biomarker Technologies Company Limited (Beijing, China). The sequencing indicator was set as 20 Gb clean data for each bulk.
SNP calling and BSR association mapping
After the cDNA libraries were sequenced, the raw data were filtered. Firstly, the reads with adapter, ribosome RNA sequences, reads with N proportion ≥ 10 and poor-quality reads with more than 50% bases with a Q value ≤ 10 were eliminated, and high-quality clean data was obtained; secondly, the STAR software (https://github.com/alexdobin/STAR) was used to assemble the clean data using the reference genome of Chinese Spring (v1.1); thirdly, SNP calling was carried out following the reference flowchart aimed at RNA-Seq by the software GATK (v3.1-1), and the SNP index values in the two bulks was calculated using the MutMap method (Abe et al., 2012) with SNPs in the susceptible bulk as a reference; fourthly, the ΔSNP index between resistant and susceptible bulks for each SNP was calculated (Takagi et al. 2013) using the following formula: ΔSNP_index = (SNP_index of resistance bulk) − (SNP_index of susceptible bulk); fifthly. Euclidean distance (ED) algorithm (Trapnell et al. 2010) was carried out to lock the associated interval based on the ΔSNP_index using the associated threshold value of 0.21. The interval with an associated threshold value > 0.21 was considered to be the candidate interval.
Development of new markers and molecular mapping of the Pm gene
After confirming the candidate interval, the markers linked to catalogued Pm genes in the candidate region were firstly used to screen polymorphic markers to preliminarily map the Pm gene in LS5082. To saturate the marker density, the SSR, Indel and other different sequences between resistant and susceptible bulks obtained from the BSR-Seq were used to develop new markers available for gel electrophoresis detection in BMK Cloud (supported by Biomarker Technologies Corporation). PCR amplification and separation and visualization of the PCR products followed Ma et al. (2015).
In addition to the markers available for gel electrophoresis detection, Kompetitive Allele-Specific PCR (KASP) markers were also developed for mapping the Pm gene(s) in LS5082. Firstly, distinctive SNPs in the candidate interval was screened; secondly, sequences of 100 bp upstream and downstream of the distinctive SNP were acquired for KASP marker development using both the Polymarker website (http://www.polymarker.info/) and Premier 5 software; thirdly, the amplification sequences of the designed primers were aligned once again in the Triticeae Multi-omics Center (http://126.96.36.199/) to ensure specificity of the sequences. The KASP primers were used in a Bio-Rad CFX real-time PCR system (Bio-Rad Laboratories, Inc, California, United States) with a final volume of 10 μL containing 3.0 μL of genomic DNA (125 ng), 5.6 μL of 2 × KASP Master Mix (provided by LGC), and 0.17 μL of primer mix (balanced mix of three pairs of primers for each marker). The amplification procedure was set as follows: 94°C for 15 min, followed by 10 touchdown cycles of 94°C for 20 s, 64 to 58°C (decreasing 0.6°C per cycle), and 38 cycles of regular amplification (94°C for 20s, 58°C for 60s), and the final ﬂuorescence was detected at 20°C using software Bio-Rad CFX Manager 3.1.
All the developed markers above were tested for polymorphisms between the resistant and susceptible parents and their derived resistant and susceptible bulks. The resulting markers were then genotyped on the F2:3 population of LS5082 × Shannong 29. After obtaining the genotyping data, a χ2 test was carried out to assess deviations of the observed phenotypic data of the F2:3 families from theoretically expected segregation ratios for goodness-of-fit analysis. The linkage map of the powdery mildew resistance gene(s) in LS5082 was constructed based on Lincoln et al. (1993) and Kosambi (1943) using the MAPMAKER 3.0 and the Kosambi function.
Differential expression analysis
After mapping clean reads (obtained from BSR-Seq) to the reference genome, FPKM (Fragments per kilo bases of exon per million fragments mapped) was carried out to calculate the expression level of the reads (Garber et al., 2011). Software EBSeq was used to detect differentially expressed genes (DEGs) using Fold Change ≥ 2 and FDR (False discovery rate) <0.01 as standard. Statistical significance of DEGs was determined using a combination of multiple tests and false discovery rate (FDR) (Reiner et al., 2003). Statistics and clustering analysis of DEGs between resistant and susceptible bulks were carried out to present the expression patterns in a genome-wide scale.
Functional annotation of the DEGs was performed using IWGSC (the International Wheat Genome Sequencing Consortium) database (v1.1). COG and KEGG pathway enrichment analyses were performed using an R package for DEGs (Yu et al., 2012). For COG analysis, the Unigene sequences were aligned to the COG database to predict possible functions, and to determine the gene function distribution characteristics (http://www.ncbi.nlm.nih.gov/COG). For KEGG pathway analysis, the KEGG database was used to blast against the metabolic pathway (https://www.kegg.jp/kegg/pathway.html).
After analysis of the DEGs using COG and KEGG, the DEGs relating to disease resistance and/or stress tolerance were selected to profile their expression patterns using RT-qPCR. To detect transcript levels, primers for specific DEGs were designed based on the coding sequences of the selected genes (Table S1). Firstly, the seedlings of LS5082 and Shimai 22 were inoculated with the Bgt isolate E09 at the one-leaf stage, and the first leaves of LS5082 and Shimai 22 seedlings were sampled 3, 6, 12, 24, 36, 48 and 72 h after inoculation, with three parallel experiments; secondly, total RNA was extracted using the Spectrum Plant Total RNA kit (Sigma-Aldrich) following the manufacturer’s recommendations, and then quantified by measuring absorbance at the wavelengths of 260 and 280 nm using a Nano Drop 1000 spectrophotometer (Thermo Scientific); thirdly, Promega DNase I-treated RNA was used for cDNA synthesis using Invitrogen SuperScript-II reverse transcriptase following the manufacturer’s guidelines. Finally, the cDNA was analyzed using RT-qPCR, following the procedure described by He et al. (2016; 2018) using SYBR green master mix (Applied Biosystems) with a Rotor-Gene-Q (Qiagen). Amplification was followed by melt curve analysis. Relative quantification was carried out using the 2−ΔΔCt method (He et al. 2018). Oligonucleotides amplifying ACTIN were used for normalization.
Evaluation of the closely linked markers for MAS
To evaluate markers for MAS, the closely linked markers were genotyped for LS5082 and 27 susceptible wheat cultivars from different regions of China. Polymorphisms between LS5082 and 27 susceptible wheat cultivars were analyzed, and the markers that stably amplified polymorphic band(s) between LS5082 and the susceptible cultivars were regarded as effective for MAS in these genetic backgrounds. To transfer the Pm gene to applicable backgrounds, these cultivars were crossed with LS5082 to construct F2 and F3 segregation populations for MAS.