Systematic dissection of disease resistance to southern corn rust by bulked-segregant and transcriptome analysis

Xiaohuan Mu, Zhuangzhuang Dai, Zhanyong Guo, Hui Zhang, Jianping Yang, Xinke Gan, Jiankun Li,Zonghua Liu, Jihua Tang, Mingyue Gou

State Key Laboratory of Wheat and Maize Crop Science, Collaborative Innovation Center of Henan Grain Crops, College of Agronomy, Henan Agricultural University,Zhengzhou 450002, Henan, China

Keywords:Maize Southern corn rust BSA-seq RNA-seq Lignin

ABSTRACT Southern corn rust (SCR) is a destructive maize disease caused by Puccinia polysora Underw. To investigate the mechanism of SCR resistance in maize,a highly resistant inbred line,L119A,and a highly susceptible line, Lx9801, were subjected to gene mapping and transcriptome analysis. Bulked-segregant analysis coupled with whole-genome sequencing revealed several quantitative trait loci (QTL) on chromosomes 1, 6, 8, and 10. A set of 25 genes, including two coiled-coil nucleotide-binding site leucine-rich repeat (CC-NBS-LRR) genes, were identified as candidate genes for a major-effect QTL on chromosome 10. To investigate the mechanism of SCR resistance in L119A, RNA-seq of P. polysorainoculated and non-inoculated plants of L119A and Lx9801 was performed. Unexpectedly, the number of differentially expressed genes in inoculated versus non-inoculated L119A plants was about 10 times that of Lx9801,with only 29 common genes identified in both lines,suggesting extensive gene expression changes in the highly resistant but not in the susceptible line.Based on the transcriptome analysis,one of the CC-NBS-LRR candidate genes was confirmed to be upregulated in L119A relative to Lx9801 independently of P.polysora inoculation.Gene Ontology and Kyoto Encyclopedia of Genes and Genomes analyses indicated that transcription factors, as well as genes involved in defense responses and metabolic processes, were dominantly enriched, with the phenylpropanoid biosynthesis pathway most specifically activated. Consistently, accumulation of phenylpropanoid-derived lignin, especially S lignin, was drastically increased in L119A after P. polysora inoculation, but remained unchanged in Lx9801, suggesting a critical role of lignin in SCR resistance. A regulatory network of defense activation and metabolic change in SCR-resistant maize upon P. polysora infection is described.

1. Introduction

Maize (Zea mays L.) provides food, feed, energy, and forage worldwide [1]. With rapid growth in global population and livestock production, the demand for maize is continually increasing.However, maize production is threatened by various diseases.Southern corn rust (SCR) is a globally destructive maize disease caused by Puccinia polysora Underw. [2,3], which attaches to leaf surfaces and inhibits photosynthesis. Yield losses caused by SCR may reach 45%–100% under optimal environmental conditions[4,5].To combat SCR,it is essential to identify SCR resistance genes and characterize the molecular basis of SCR resistance.

Several studies [2,3,6–9] have been conducted to identify genetic and molecular mechanisms of maize resistance to SCR.Candidate loci and/or genes for SCR resistance were identified on the short arm of chromosome 10[2,6,9–13].Recently,two quantitative trait loci (QTL), qSCR4.01 and qSCR6.01 on chromosomes 4 and 6, were identified [7,14]. Using a proteomic approach in SCRresistant and -susceptible maize inbred lines, ZmREM1.3, a differentially expressed remorin protein, was identified. Resistance to SCR was increased in ZmREM1.3-overepressing plants but reduced in a mutant of ZmREM1.3, confirming the gene’s role in SCR resistance [3]. However, no SCR resistance gene has been functionally confirmed by forward-genetic studies. Systematic dissection of SCR resistance mechanisms using resistant and susceptible germplasm awaits investigation.

Bulked-segregant analysis coupled with whole-genome sequencing(BSA-seq)is an effective method for quickly identifying candidate genes or genomic regions controlling a given phenotype[15,16].Candidate genes for Fusarium wilt and sterility mosaic disease in pigeonpea [Cajanus cajan (L.) Millsp.] were identified by BSA-seq[17].RNA sequencing(RNA-seq)is an efficientway to assess global gene expression profiling and has been used to identify differentially expressed genes (DEGs), thereby dissecting the maize defense response to infections of Colletotrichum graminicola,Fusarium verticillioides, F. graminearum, Cercospora zeae-maydis, C. zeina,and Setosphaeria turcica [18–23]. Although proteomics has been employed to identify differentially expressed proteins in maize response to P. polysora [3], there is a lack of systematic transcriptomic analysis of maize genes following infection by the fungus.

In this study, integrated BSA-Seq, transcriptomic and physiological analyses were performed to uncover the mechanism underlying the difference of SCR resistance in the highly resistant line L119A versus the susceptible line Lx9801.The candidate SCR resistance genes and critical defense and metabolic pathways conferring SCR resistance were identified.

2. Materials and methods

2.1. Plant materials, genetic analysis and BSA-seq

A susceptible maize line Lx9801, a resistant line L119A, and BC4F1and BC4F2populations were tested for SCR resistance in field experiments at the experimental station in three locations: Xinxiang and Shangqiu in Henan province and Ledong in Hainan province,China.Plants were evaluated as resistant or susceptible on a five-point rating scale[6]:1,immune to SCR without visible infection on leaves; 3, highly resistant to SCR with chlorotic flecks but without uredinia on leaves; 5, resistant to SCR with few uredinia on leaves;7,susceptible to SCR with moderate numbers of uredinia on leaves; 9, highly susceptible to SCR with large numbers of uredinia on the whole plant.

Leaves of 25 highly SCR-resistant (with score 1) and 25 highly susceptible(with score 9)lines from BC4F2and parental lines were bulked for DNA extraction by the cetyltrimethylammonium bromide method [24]. Equal amounts of bulked DNA of resistant plants(R-bulk)and susceptible plants(S-bulk),as well as their parents,were subjected to whole-genome resequencing using an Illumina HiSeq (Illumina, Inc., San Diego, CA, USA) platform (Illumina NovaSeq6000), producing 75 Gb data for the R- and S-bulks,respectively, and 25 Gb data for each parent. After removal of low-quality regions and adapter sequences, clean data were obtained and then aligned to the maize B73 RefGene_v4 reference genome with Burrows-Wheeler Aligner [25]. Grouping and duplicate-read identification were performed with Picard(https://broadinstitute.github.io/picard/). Finally, GATK [26] was used to realign suspicious intervals, and to call and filter variants.QTLseqr[27]was performed to detect QTL.Annotation of SNPs and InDels was performed with SnpEff [28]. The statistics of BSA-seq were calculated as previously described [29]. The threshold of the G statistic(G’)was defined as the genome-wide false discovery rate (FDR) of 0.001. Peak regions above the threshold value were defined as association regions. Preliminary QTL regions on each chromosome were first identified by sliding-window analysis of G’ with a 5 Mb window, and the QTL region on chromosome 10 was then refined by scanning the G’ and delta SNP-index with a 200 kb window [30].

2.2. RNA-seq and data analysis

For RNA-seq samples, L119A and Lx9801 were grown in pots(25 cm wide × 17 cm high) under natural conditions at Henan Agricultural University, Zhengzhou, Henan province. Twelve pots of each inbred line with two maize plants per pot were cultivated.Pots were arranged randomly. A P. polysora strain isolated from field-grown leaves in Xinxiang was used for plant inoculation. At V6 to V7 stage (six to seven expanded leaves), half of the plants were inoculated by spraying with 10 mL solution consisting of 0.02% Tween 20 (v/v) and approximately 105–106spores mL-1.The remaining plants as control were mock-inoculated with 10 mL water containing 0.02%Tween 20.At 14 days after inoculation, the middle parts of inoculated leaves were sampled. Four replicates with 2 plants per replicate were collected for each treatment. The samples were frozen in liquid nitrogen for RNA extraction and cell wall preparation. B73 was cultivated at the Xinxiang experimental station. At silking stage, plants were separated into leaf, tassel, silk, husk, and cob, immediately frozen in liquid nitrogen, and stored at –80 °C. Three replicates with three plants per replicate were collected.

RNA was extracted from samples with Trizol reagent and examined with a Nanodrop spectrophotometer(Thermo Fisher Scientific Inc.,Wilmington,DE,USA)and a 2100 Bioanalyzer(Agilent Technologies,Santa Clara,CA,USA).RNA samples of three biological replicates with sufficient quality were used for RNA-seq library construction on the Illumina HiSeq platform (Illumina Hiseq Xten)(Illumina,Inc.,San Diego,CA,USA).Clean reads were mapped to the maize reference genome using HISAT2[31].Alignments were processed with StringTie [32] software to assemble transcript isoforms and quantify expression value as fragments per kilobase of exon model per million mapped reads(FPKM)of known and novel genes.Gene expression profile data for principal components analysis (PCA) plot was preprocessed by Pareto scaling. Differentially expressed genes with fold change(FC)＞2 and false discovery rate(FDR)＜0.05 were identified with the R package DESeq2[33].Venny 2.1 (http://bioinfogp.cnb.csic.es/tools/venny/index.html) was used to construct a Venn figure. AgriGO_v2.0 (http://systemsbiology.cau.edu.cn/agriGOv2/)[34]was used for Gene Ontology(GO)analysis.The R packages‘‘clusterProfiler”and‘‘pathview”[35]were used for Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis.Shared DEGs in Lx9801 and L119A before and after P.polysora inoculation were defined as common genes(CGs),whereas DEGs unique to either Lx9801 or L119A were defined as specific genes(SGs).Subcellular location analysis was performed with SUBA4 (http://suba.live/) [36]. Transcription factors (TFs) were identified with PlantTFDB(http://planttfdb.cbi.pku.edu.cn/index.php)[37].Bubble plots were drawn with the‘‘ggplot2”package[38].Gene classification was performed with MapMan Data (https://mapman.gabipd.org/mapmanstore). For any input gene list, Plant Regulomics(http://bioinfo.sibs.ac.cn/plant-regulomics/) retrieves the factors,treatments,and experimental/environmental conditions regulating the input from integrated omics data.

2.3. Cell wall preparation and lignin measurement

Fresh leaves sampled for the RNA-seq experiment were used to measure lignin. Four replicates with two plants for each replicate were collected for each treatment. Leaf samples were ground into powder in liquid nitrogen,extracted with 70%ethanol,and centrifuged at 10,000 r min-1for 10 min. The pellet was resuspended and extracted with chloroform:methanol (1:1) three times. The residues were resuspended and extracted with acetone.The resulting pellet was dried at 35 °C. Total lignin was measured by the acetyl bromide method and G and S lignin monomers were measured by the thioacidolysis method, as previously described [39].

2.4. Validation of RNA-seq and candidate gene expression pattern by quantitative real time-PCR (qRT-PCR)

To validate the RNA-seq results, three biological replicates of RNA used in the RNA-seq analysis were subjected to qRT-PCR. To assess the expression patterns of candidate genes, total RNA extracted from leaf, tassel, silk, husk, and cob of B73 were used in qRT-PCR. Total RNA was reverse-transcribed into cDNA using RNase-free DNase Kit (Cat. #RR047A, TaKaRa, Otsu Shiga, Japan).Expression profiles of genes were quantified with a SYBR Green system in a LightCycler 480(Roche Applied Science,Basel,Switzerland).The primers used in qRT-PCR analyses are listed in Table S1.The maize ZmUBQ gene was used as an internal control.

2.5. Statistical analysis

Chi-squared (χ2) tests were performed to determine the goodness of fit of the segregations in the BC4F1(1:1) and BC4F2(3:1)populations. Data for the lignin contents and qRT-PCR were subjected to variance analysis using one-way ANOVA procedure implemented in SPSS 19.0(SPSS,Inc.,Chicago,IL,USA).Differences were compared by least significant difference test at P ＜0.05.

3. Results

3.1. Characterization of the SCR-resistant phenotype in L119A and its backcross populations

Line L119A was highly resistant and Lx9801 highly susceptible to SCR (Fig. 1A). A BC4F1population was generated by backcross processes using L119A as donor and Lx9801 as recurrent parent.The BC4F1plants were then self-pollinated to generate the BC4F2population. Both BC4F1and BC4F2plants were divided into resistant (score 1–5) and susceptible (score 7–9) groups. Among the 449 BC4F1plants scored, 224 were resistant and 225 susceptible(Fig. 1B; Table S2), fitting a 1:1 ratio well (χ2= 0.0022,P = 0.9624). For the 438 BC4F2plants scored, 319 were resistant and 119 susceptible (Fig. 1C; Table S2), fitting a 3:1 ratio(χ2= 1.0989, P = 0.2945). These results suggested the action of a major-effect gene conferring SCR resistance in L119A.

Fig. 1. The phenotypic distributions in BC4F1 and BC4F2 populations from the cross between L119A and Lx9801. (A) Phenotypes of P. polysora-infected plants of L119A,Lx9801, and their BC4F2 populations. Disease level is based on the five-point rating scale. (B, C) Plant frequency distribution of disease level in BC4F1 (B) and BC4F2 (C)populations.Numbers of level 1,3 and 5 are assigned to the resistant group(R),and numbers of level 7 and 9 to the susceptible group(S).The BC4F1 population consisted of 224 resistant and 225 susceptible plants,fitting a 1:1 ratio(χ2=0.0022,P=0.9624).The BC4F2 population consisted of 319 resistant and 119 susceptible plants,showing 3:1 Mendelian segregation (χ2 = 1.0989, P = 0.2945).

3.2.Identification of genetic loci conferring SCR resistance in L119A by BSA-seq

Eight QTL on chromosomes 1,6,8,and 10 were identified in the G’ plot (Fig. 2A; Table S3). The highest peak was located on chromosome 10, and the locus (QTL8) was accordingly assigned as a major-effect QTL (Fig. 2A). QTL8 was further confined to a 400 kb region (1,397,359–1,797,359) on chromosome 10 with a peak at 1,597,359 by scanning the G’ and delta SNP-index with a 200 kb window. Twenty-five candidate genes were identified in this region based on the B73 reference genome, including two coiledcoil nucleotide-binding site leucine-rich repeat (CC-NBS-LRR)genes, Zm00001d023265 and Zm00001d023267 (Table S4). qRTPCR analysis indicated that the expression of both Zm00001d023265 and Zm00001d023267 was higher in L119A than in Lx9801 independently of P.polysora inoculation(Fig.2B,C).Both genes were expressed in all tested organs,with the highest expression in leaves (Fig. 2D, E). We speculated that Zm00001d023265 and Zm00001d023267 are the most reasonable candidate genes conferring SCR resistance in L119A.

Fig. 2. Identification of genetic loci conferring SCR resistance in L119A by BSA-seq. (A) Detection of SCR-resistance QTL using BSA-seq. (B) Relative gene expression level of Zm00001d023265 in non-inoculated and inoculated plants of L119A and Lx9801. (C) Relative gene expression level of Zm00001d023267 in non-inoculated and inoculated plants of L119A and Lx9801.(D)Relative gene expression level of Zm00001d023265 in five maize organs.(E)Relative gene expression level of Zm00001d023267 in five maize organs. Different letters above error bars in (B–E) indicate significant difference at P ＜0.05.

Fig. 3. Number of differentially expressed genes (A) and Venn diagram (B) for inoculated versus non-inoculated plants of L119A and Lx9801.

3.3. Transcriptome analysis of inoculated and non-inoculated L119A and Lx9801 by RNA-seq

To understand the molecular mechanism of SCR resistance, we performed the transcriptome analysis of L119A and Lx9801 inoculated and non-inoculated with P.polysora.Generally,similar numbers of transcripts were detected by RNA-seq in non-inoculated L119A (32360) and Lx9801 (31876) as in inoculated L119A(32411)and Lx9801(32891).PCA plot indicated that the replicates of different treatments were clustered into distinct groups, suggesting the reliability of the RNA-seq data (Fig. S1).

In L119A, 727 and 89 genes were significantly up- or downregulated in inoculated versus non-inoculated plants (Fig. 3A;Table S5). In contrast, only 68 genes were up-regulated and 18 down-regulated in inoculated versus non-inoculated plants of Lx9801(Fig. 3A; Table S6). Thus,＞10-fold more DEGs were identified in L119A than in Lx9801 (Fig. 3A; Tables S5 and S6). qRT-PCR analysis of six randomly chosen DEGs (Zm00001d031815,Zm00001d049217, Zm00001d038049, Zm00001d008862,Zm00001d046234, and Zm00001d051934) showed good agreement with the RNA-seq data(Fig.S2).None of the 25 candidate genes on chromosome 10 (Table S7) showed differential expression in inoculated versus non-inoculated plants in either L119A or Lx9801.However,the expression of Zm00001d023265 was upregulated in L119A versus Lx9801 independently of P.polysora inoculation (Table S7), in agreement with the qRT-PCR data (Fig. 2B).Different from the qRT-PCR data (Fig. 2C), the expression of Zm00001d023267 was not significantly changed in Lx9801 versus L119A. Zm00001d023248 and Zm00001d023259 were not detected in non-inoculated plants of Lx9801 and L119A, but were upregulated and downregulated, respectively, in inoculated Lx9801 versus inoculated L119A (Table S7).

Only 29 DEGs were shared by L119A and Lx9801, and were assigned as CGs (Fig. 3B; Table S8). There were 57 Lx9801-specific DEGs and 787 L119A-specific DEGs, which were classified as SGs (Fig. 3B; Tables S5 and S6). GO-term enrichment analysis detected 7 biological processes (BP) and 1 molecular function(MF) that were enriched in CGs (Table S9), with mainly the terms‘‘carboxylic acid metabolic process”,‘‘defense response to fungus”,‘‘oxoacid metabolic process”,and‘‘organic acid metabolic process”represented (Fig. 4A; Table S9). Seven BP were enriched in SGs in Lx9801, with mainly the terms ‘‘response to fructose/ hexose/mo nosaccharide/sucrose/disaccharide” represented (Fig. 4A; Table S10). In contrast, 229 BP, 50 MF, and 4 cellular component (CC)GO terms were enriched in the SGs of L119A (Fig. 4A; Table S11).The SGs of L119A were preferentially associated with defenserelated GO terms such as ‘‘response to stimulus/other organism/oxygen-containing compound/hormone/chitin/chemical”, ‘‘defense response”, ‘‘cellular response to acid chemical” (Fig. 4A; Table S11). Six genes encoding pathogenesis-related (PR) proteins(Zm00001d018734, Zm00001d018738, Zm00001d009296,Zm00001d023811, Zm00001d028816, and Zm00001d029558) and four genes encoding NBS-LRR resistance (NLR) proteins(Zm00001d011737, Zm00001d014650, Zm00001d047952, and Zm00001d052927) were specifically up-regulated in L119A,whereas only one NLR gene (Zm0001d014650) was up-regulated in both L119A and Lx9801 (Fig. 4B; Tables S5 and S6). A set of 787 SGs of L119A were tested with Plant Regulomics online software [40]. Of 743 listed genes, 712 were shared with published genes induced by pathogen infection (Table S12), confirming the extensive activation of defense responses in L119A. In general,genes involved in defense responses were specifically induced in L119A but not in Lx9801, a finding consistent with the difference in disease symptoms (Fig. 1A). We thus investigated these SGs in L119A in view of their potential involvement in SCR resistance.Most proteins encoded by the SGs of L119A were predicted by SUBA4 to be located in cytosol and nucleus (Fig. 5A), implying the involvement of a large number of TFs in the SGs of L119A.Consistently, 74 TFs were identified in the SGs of L119A by PlantTFDB(Fig. 5B). These TFs belong to 18 gene families, among which WRKY, MYB, NAC, and ERF were the four most highly enriched TF families. Seventeen WRKY and 11 MYB TFs were upregulated in inoculated versus non-inoculated L119A plants (Fig. 6; Table S13). However, expression of those genes was not significantly changed in Lx9801 (Fig. 6; Table S13). Many genes involved in redox state regulation were upregulated in L119A, while almost all were downregulated or unchanged in Lx9801(Fig.6;Table S14).

Fig.4. GO and KEGG analysis of differentially expressed genes in inoculated versus non-inoculated plants of Lx9801 and L119A.(A)GO term analysis of common genes(CGs)and specific genes(SGs)in L119A and Lx9801.(B)Expression changes in representative SGs involved in disease resistance in inoculated versus non-inoculated plants of L119A and Lx9801. Colors correspond to log2 values of fold change (FC). (C) KEGG pathway analysis of SGs in inoculated versus non-inoculated plants of L119A.

Fig. 5. Prediction of protein localization and transcription-factor enrichment. (A) The proportion of predicted subcellular location of proteins encoded by specific genes in L119A. (B) Gene number of each transcription factor family enriched in specific genes in maize L119A.

Fig.6. Differential changes of defense-associated gene expression in inoculated versus non-inoculated plants of L119A and Lx9801. Colors correspond to log2 values of fold change(FC)of inoculated versus non-inoculated plants.The upper and lower(separated by blanks)cells represent the expression data from L119A and Lx9801,respectively.Each cell represents one gene in each gene category.

3.4. SGs of primary and secondary metabolisms in L119A

The KEGG analysis indicated that the SGs of L119A were enriched mainly in pathways involved in biosynthesis of primary metabolites, especially amino acids (tryptophan, cysteine, and methionine) and secondary metabolites (phenylpropanoid, flavonoid, zeatin, and diterpenoids) (Fig. 4C). Based on the gene classification by MapMan Data, most of the phenylpropanoid-lignin biosynthetic genes including those encoding two phenylalanine ammonia lyases (Zm00001d003016 and Zm00001d051161), four hydroxycinnamoyl coenzyme A: shikimate hydroxycinnamoyl transferase (Zm00001d022592, Zm00001d050455,Zm00001d030542, and Zm00001d037073), one 4-coumarate-CoA ligase (Zm00001d032103), one caffeate O-methyltransferase(Zm00001d048087), and one cinnamoyl-CoA reductase(Zm00001d019669) were upregulated in L119A after P. polysora inoculation (Fig. 6; Table S15). Flavonoid biosynthetic genes including two chalcone synthase (Zm00001d052673 and Zm00001d007403), two dihydroflavonol-4-reductases(Zm00001d031802 and Zm00001d020961), and one flavonone-3-hydroxylase (Zm00001d001960) were upregulated in L119A after P.polysora inoculation(Fig.6;Table S15).However,the expression of most of these genes was not significantly changed in inoculated versus non-inoculated Lx9801 (Fig. 6; Table S15).

Fig.7. Quantification of lignin content in inoculated versus non-inoculated plants of L119A and Lx 9801.(A)Total lignin content measured by acetyl bromide method.(B,C)Content of G lignin(B)and S lignin(C)measured by the thioacidolysis method.(D)Ratio of S lignin to G lignin content presented in(B)and(C).Different letters above error bars in (A–D) indicate significant difference at P ＜0.05.

3.5. Lignin accumulation in L119A and Lx9801

Total lignin content was generally higher in L119A than in Lx9801 independently of P.polysora inoculation(Fig.7A).Total lignin content was 19% higher in inoculated versus non-inoculated L119A but remained unchanged in inoculated versus noninoculated Lx9801 (Fig. 7A). Consistently, both G and S lignin monomer levels increased in inoculated versus non-inoculated L119A,but the increase in S lignin(71%)was higher(P ＜0.05)than that in G lignin (17%) (Fig. 7B, C), resulting in a higher S/G ratio in inoculated L119A(Fig.7D).However,there appeared to be a subtle decrease in G lignin but no significant change in S lignin in inoculated versus non-inoculated Lx9801 (Fig. 7B, C). The S/G ratio in Lx9801 was generally higher(P ＜0.05)than that in L119A independently of P.polysora inoculation(Fig.7D).Thus,lignin biosynthesis was induced specifically in the SCR-resistant but not in the SCRsusceptible line,and the induction rates of G and S lignin monomer differed.

4. Discussion

BSA-seq was performed in this study to identify the causal locus of the SCR resistance in the resistant line L119A,8 QTL were identified, and a major-effect QTL was located on chromosome 10.Many QTL or genes conferring SCR resistance in maize have been reported previously [2,6–9,12,14]. On chromosome 1, a stable SCR-resistance QTL was identified in a 100.9 Mb region (chromosome 1: 93,193,678–194,118,967) between simple sequence repeat (SSR) markers Umc2025 and Umc1919 based on the B73 V4 reference genome[8].In our study,QTL1 and QTL2 were identified in respectively the 10.8 Mb (chromosome 1: 43,399,099–54,212,595) and 6.4 Mb (chromosome 1: 289,124,388–295,526,4 80) regions on chromosome 1 and thus represented novel QTL for SCR resistance. On chromosome 6, a SCR-resistance QTL,qSCR6.01, was mapped to a 2.95 Mb region based on the B73 RefGen_v3 reference genome corresponding to a 3.08 Mb region(chromosome 6:77,998,607–81,075,445)on B73 RefGen_v4 reference genome [7]. QTL3 and QTL4 were identified in the 27.5 Mb(chromosome 6: 100324867–127501574) and 11.1 Mb (chromosome 6:133026374–144099079)regions,respectively,with a physical distance of about 20–50 Mb from qSCR6.01. On chromosome 8, a minor SCR-resistance QTL was previously identified between SSR markers Umc1360 and Umc1034 with genetic positions of 51.7 and 70.3 cM, corresponding to a 51.84 Mb region (chromosome 8:20,552,200–72,390,224)on B73 RefGen_v4 reference genome[12].QTL5,QTL6,and QTL7 detected in this study were located in 31.4 Mb(chromosome 8:6,455,717–37,891,932),19.7 Mb(chromosome 8:73,899,522–93,569,844),and 21.8 Mb(chromosome 8:134,888,636–156,706,016) regions, respectively, with an approximately 10 Mb overlap with the published QTL on chromosome 8[12]. Several studies have identified major QTL for SCR resistance on chromosome 10. RppP25, a SCR-resistance gene, was anchored to a 40 kb region on chromosome 10 based on the B73 RefGen_v2 reference genome corresponding to a 96 kb region (chromosome 10:2,651,981–2,748,521)on the B73 RefGen_v4 reference genome[6].RppCML496,a major QTL for resistance to SCR,was anchored to a 128 kb region on chromosome 10 (chromosome 10: 2,640,817–2,768,842) [13]. Another SCR-resistance locus, RppM, was recently anchored to chromosome 10 (chromosome 10:1,586,659–1,697,392)on the B73 RefGen_v4 reference genome based on BSAseq combined with fine-mapping in an F2population derived from a cross between Jing2416K and Jing2416[9].The major-effect QTL(QTL8) was mapped to a 6.7 Mb region of chromosome 10 (chromosome 10: 7615–6,766,893), overlapping those of RppP25,RppCML496, and RppM. By further refining the preliminary QTL8 on chromosome 10 to a 400 kb region(1,397,359–1,797,359)with 200 kb window according to G’and delta SNP-index,25 candidate genes were identified in this region based on B73 reference genome, among which Zm00001d023265 and Zm00001d023267, the two genes highly expressed in leaves, encode CC-NBS-LRR type R proteins. Those genes have also been proposed [9] as candidate SCR-resistance genes for RppM.Although both genes were downregulated in Lx9801 versus L119A independently of P.polysora inoculation, only expression of Zm00001d023265 was downregulated in Lx9801 versus L119A in the RNA-seq data (Table S7). While those information highlights the two CC-NBS-LRR genes, especially Zm00001d023265, as major candidate genes, we could not exclude other candidate genes. It remains to be determined whether additional genes are present in this region of the L119A and Lx9801 genome, given that the B73 genome was used as the reference for BSA-seq analysis.

Unexpectedly, based on the transcriptome data, the number of DEGs in inoculated versus non-inoculated L119A was about 10 times that in Lx9801, with only 29 CGs identified, suggesting extensive gene expression change in the highly resistant but not the susceptible line. We suspect that the mild change of gene expression in Lx9801 upon P.polysora inoculation is due to the biotrophic characteristics of this pathogen, which does not cause obvious leaf morphological and physiological change after infection.Besides,plant immune response did not appear to be induced in inoculated Lx9801 according to the SGs identified in Lx9801,probably owing to a lack of functional immune receptors like the CC-NBS-LRR protein described above. Consistently, SGs in Lx9801 were mostly genes involved in carbon metabolism, likely a consequence of photosynthesis inhibition by the disease itself. In contrast, the SGs in L119A were enriched with defense-responsive genes (Fig. 4A). ＞90% of the SGs in L119A were shared with published genes induced by pathogen infection,confirming the intense activation of defense response in L119A.

Plant defense response is regulated by an elaborate regulatory network consisting of immune receptors, TFs, and genes controlling reactive oxygen species accumulation and secondary metabolite biosynthesis etc. [41–44]. Upon perception of pathogen infection by the immune receptors like CC-NBS-LRR proteins, TFs,especially those of WRKY and MYB families, function in transcriptional regulation of defense response genes[45–49].After inoculation of P. polysora, expression levels of 74 TFs including 17 WRKY TFs and 11 MYB TFs were changed in L119A,indicating the activation of a transcription regulatory network in L119A in response to P.polysora.Along with the activation of upstream TFs,downstream defense-associated genes were upregulated in L119A, indicating the broad activation of downstream defense genes. Twenty-two redox-associated genes, including 11 genes encoding glutathione S-transferases and 11 genes encoding peroxidases, were upregulated in L119A after inoculation(Fig.6;Table S14),suggesting that the redox state in L119A is beneficial to SCR resistance.

As well documented in previous studies [50–52], MYB transcription factors have also been implicated in the transcriptional regulation of secondary metabolite biosynthesis,including phenylpropanoid, lignin, and flavonoid biosynthesis. Among those MYB TFs, AtMYB15 specifically regulates pathogen-triggered lignin biosynthesis to combat disease [48]. Along with the expression changes of the 11 MYB TF genes, a series of genes involved in secondary metabolism were upregulated in L119A upon P. polysora infection, but remained unchanged in Lx9801 (Fig. 6).Phenylpropanoid-lignin biosynthetic genes encoding phenylalanine ammonia lyase, hydroxycinnamoyl coenzyme A: shikimate hydroxycinnamoyl transferase, 4-coumarate-CoA ligase, caffeate O-methyltransferase, and cinnamoyl-CoA reductase were all upregulated in L119A after P. polysora inoculation. In our previous studies [39,53], genes encoding the membrane steroid-binding protein (Zm00001d017380) and the cytochrome b5 protein(Zm00001d017425) functioned in lignin biosynthesis by modulating P450 protein complex formation and electron transfer. Both genes were specifically upregulated in inoculated versus noninoculated L119A but not in Lx9801 (Tables S5 and S6). In agreement with the specific induction of lignin biosynthetic genes in L119A,total lignin content in L119A was generally higher than that in Lx9801, and both G and S lignin levels were increased in inoculated versus non-inoculated L119A but slightly decreased or remain unchanged in inoculated versus non-inoculated Lx9801.The lignin biosynthetic gene ZmCCoAOMT2 was involved in multiple disease resistance to southern leaf blight and gray leaf spot,and ZmCAD was essential for leaf and sheath blight resistance in maize[54,55].Those studies and our findings collectively suggest that lignin plays a vital role in resistance to SCR. Besides, the finding that the increase in S lignin was dramatically higher than that in G lignin in L119A suggests that S-lignin is more responsive to SCR.This is consistent with previous studies on the effect of lignin composition on Pseudomonas syringae resistance in Arabidopsis [56].

In conclusion, we have identified several QTL contributing to SCR resistance in L119A, and observed extensive gene expression change in L119A but not in Lx9801. Based on integration of BSAseq and RNA-seq, we propose that P. polysora is recognized likely by immune receptors such as the CC-NBS-LRR protein in L119A but not in Lx9801. The immune receptors may then activate some key TFs, which further regulate the expression of a series of defense-associated genes, including PR genes, redox-associated genes,and secondary-metabolite biosynthetic genes.Upregulation of lignin biosynthetic genes, as well as the specifically induced accumulation of lignin in L119A, suggest that lignin plays vital roles in SCR resistance. Identification of the causal gene and effect of lignin in SCR resistance awaits further studies to fully understand the mechanism of SCR resistance.

Data availability

Data of RNA-seq in this study have been deposited in the SRA database (https://www.ncbi.nlm.nih.gov/sra; accession ID PRJNA689981).

CRediT authorship contribution statement

Xiaohuan Mu:data curation, formal analysis, visualization,writing of original draft;Zhuangzhuang Dai:data curation,investigation;Zhanyong Guo:data curation for BSA-seq;Hui Zhang:resources,creation of BC4F1 population;Jianping Yang:data curation for lignin content;Xinke Gan:data curation for qRT-PCR;Jiankun Li:resources, cultivation of maize;Zonghua Liu:resources,creation of L119A line;Jihua Tang:resources, creation of BC4F1 population, conceptualization;Mingyue Gou:conceptualization,funding acquisition, project Administration, writing, review, and editing.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work was supported by the Zhongyuan Thousand Talents Program (ZYQR201912168, to MG), the National Natural Science Foundation of China (U2004207, to MG), Fund for Distinguished Young Scholars in Henan (212300410007) and the Startup Grant of Henan Agricultural University (30601732, to MG and 30500926, to XM).

Appendix A. Supplementary data

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2021.07.001.