Transcriptome analysis reveals the molecular mechanisms of the defense response to gray leaf spot disease in maize

Background Gray leaf spot (GLS), which is caused by the necrotrophic fungi Cercospora zeae-maydis and Cercospora zeina, is one of the most impactful diseases in maize worldwide. The aim of the present study is to identify the resistance genes and understand the molecular mechanisms for GLS resistance. Results Two cultivars, ‘Yayu889’ and ‘Zhenghong532,’ which are distinguished as resistant and susceptible cultivars, respectively, were challenged with the GLS disease and a RNA-seq experiment was conducted on infected plants at 81, 89, 91, and 93 days post planting (dap). Compared with the beginning stage at 81 dap, 4666, 1733, and 1166 differentially expressed genes (DEGs) were identified at 89, 91, and 93 dap, respectively, in ‘Yayu889,’ while relatively fewer, i.e., 4713, 881, and 722 DEGs, were identified in ‘Zhenghong532.’ Multiple pathways involved in the response of maize to GLS, including ‘response to salicylic acid,’ ‘protein phosphorylation,’ ‘oxidation-reduction process,’ and ‘carotenoid biosynthetic process,’ were enriched by combining differential expression analysis and Weighted Gene Co-expression Network Analysis (WGCNA). The expression of 12 candidate resistance proteins in these pathways were quantified by the multiple reaction monitoring (MRM) method. This approach identified two candidate resistance proteins, a calmodulin-like protein and a leucine-rich repeat receptor-like protein kinase with SNPs that were located in QTL regions for GLS resistance. Metabolic analysis showed that, compared with ‘Zhenghong532,’ the amount of salicylic acid (SA) and total carotenoids in ‘Yayu889’ increased, while peroxidase activity decreased during the early infection stages, suggesting that increased levels of SA, carotenoids, and reactive oxygen species (ROS) may enhance the defense response of ‘Yayu889’ to GLS. Conclusion By combining transcriptome and proteome analyses with comparisons of resistance QTL regions, calmodulin-like protein and leucine-rich repeat receptor-like protein kinase were identified as candidate GLS resistance proteins. Moreover, we found that the metabolic pathways for ROS, SA, and carotenoids are especially active in the resistant cultivar. These findings could lead to a better understanding of the GLS resistance mechanisms and facilitate the breeding of GLS-resistant maize cultivars. Electronic supplementary material The online version of this article (10.1186/s12864-018-5072-4) contains supplementary material, which is available to authorized users.


Background
Gray leaf spot (GLS), which is caused by the necrotrophic fungi Cercospora zeae-maydis and Cercospora zeina, is one of the most serious foliar diseases in maize among almost all maize-growing regions worldwide [1,2]. Cercospora spores overwinter on maize plant debris in diseased fields, and conidia begin to germinate under suitable climatic conditions (i.e., warm temperatures and high humidity) in the next planting season [3,4]. The disease symptoms are characterized by distinct rectangular shapes parallel to the veins [5], and the severity of GLS is highly dependent upon climate conditions and cultivars, which could cause yield losses of over 70% owing to severe blight, stalk deterioration, and lodging [2,6]. This disease has also occurred commonly in Sichuan Province, China since 2007 [7] owing to a lack of resistant cultivars and the temperate mountain climate, which favors the spread of this disease.
Compared with regular methods for disease control by chemicals, cultivation of GLS-resistant hybrids is urgently needed, as it is a cost-effective and environmentally friendly approach for reducing yield loss caused by GLS. Thus, a better understanding of the genetic architecture and the underlying genes and mechanisms could contribute to crop improvement and disease management. In general, plant immunity to necrotrophic infections is associated with complex biological and physiological processes [8], including host cell death [9,10] and the production of various secondary metabolites [11,12] and hormones [13][14][15] as well as the accumulation of reactive oxygen species (ROS) [16,17] and a variety of cell wall modifications [14,18]. Previous studies have been focused on GLS resistance QTL mapping and a variety of genetic loci were identified [19][20][21]. QTL mapping of a maize recombinant inbred line (RIL) population derived from subtropical CML444 × SC Malawi maize inbred lines identified seven GLS-resistance QTL in total from five sites, and each of two QTLs (located in bin 4.08 and bin 6.06/6.07) explained more than 11% of the phenotypic variation [22]. The maize nested association mapping (NAM) population was used to identify three QTLs that reduced GLS disease incidence by more than 10%, and experiments with near isogenic lines were conducted to develop hypotheses about resistance mechanisms [19]. Bi-parental QTL mapping was combined with GWAS to refine the map positions of GLS-resistance QTLs on chromosomes 1, 6, 7, and 8. One QTL was fine-mapped to a region with 15 candidate genes [21]. An integrated systems genetics approach recently identified four hotspots that coincide with GLS severity QTLs located on chromosomes 4, 9, and 10 [23]. Experiments with a maize RIL population in South Africa supported a hypothesized network associated with susceptibility in the maize-GLS pathosystem. However, no GLS resistance genes have been cloned, and the knowledge of maize resistance to GLS is still limited.
RNA-Seq can provide an efficient way to assess the global expression variation in coding genes and identify gene-based markers at the whole genome level. Recently, RNA-Seq analysis of resistant and susceptible sub-tropical maize lines showed that zealexins and kauralexins are accumulated after infection by C. zeina of both germplasms adapted to the southern African climate [24]. Several candidate genes for the gibberella ear rot resistance were identified by screening global gene expression profiles for differentially expressed genes mapped to within resistance QTL regions [25]. RNA-Seq was also used to identify SNPs that provide a useful resource for genetic and genomic studies and are a valuable asset for development of functional markers in soybean [26].
'Yayu889, ' as a maize cultivar with high proven resistance to GLS in field conditions of China, renders itself an excellent material for understanding maize resistance to GLS. To facilitate our analysis, the susceptible cultivar 'Zhenghong532' was used as a reference strain. The transcriptional response of these two maize cultivars to GLS was analyzed during four stages of infection by using RNA-Seq, and both common and specific differentially expressed genes (DEGs) associated with the defense response to GLS in the two cultivars were investigated. To the best of our knowledge, this study is the first to report the induced defense response to GLS in maize by using a comparative transcriptome analysis. By screening global gene expression profiles for DEGs followed by quantitative analysis of protein and comparisons to resistance QTL regions previously reported, we have identified candidate genes for GLS resistance. Moreover, we found several pathways involved in maize responses to GLS by metabolic analysis. This may ultimately help to reveal the molecular mechanisms of maize-GLS interactions, which forms an important step towards developing molecular markers-assisted selection (MAS) for GLS-resistant genetic materials.

Germplasm and field trials
Two commercial maize lines, 'Yayu889' and 'Zhenghong 532, ' were used in our study, and both of them were provided by Sichuan Nongda Zhenghong Bio. Co., Ltd., China. Both of them were planted in Guangyuan (N 32°3 6′ 19.96″, E 106°05′ 59.34″; 1354 masl), Sichuan Province, China, where GLS outbreaks occur every year. For each line, seeds were sown in three replicate rows of 10 plants. Before this study, a natural infection experiment was conducted for three consecutive years in the field, confirming that 'Yayu889' and 'Zhenghong532' were resistant and susceptible, respectively.
Leaf samples of three biological replicates were collected at 81, 89, 91, and 93 days post planting (dap) and stored in liquid nitrogen. Disease symptoms were evaluated at 81, 89, 91, 93, and 123 dap, respectively, based on a 1-5 scale with 0.25 increments, according to disease progression on the ear leaf [27].

cDNA library construction and sequencing
Total RNA was isolated and purified from leaf tissues with TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. To ensure that the RNA meets the requirements for transcriptome sequencing, a Nanodrop spectrophotometer (Thermo Fisher Scientific Inc., Wilmington, DE, USA), Qubit RNA Kit (Life Technologies, Carlsbad, Ca, USA), and 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) were used to determine the purity, concentration, and integrity of samples, respectively.
Both cDNA library preparation and sequencing were conducted by the Biomarker Technology Company in Beijing, China. All 24 libraries were sequenced using an Illumina HiSeq™4000 (Illumina, San Diego, CA). The raw data were filtered with the FASTQ_Quality_Filter tool from the FASTX-toolkit. After preprocessing the RNA-Seq data, the reads were mapped to the maize reference genome version 4 (B73 RefGen_v4) using Tophat2 [28], which was download from ftp://ftp.ensemblgenomes.org/pub/release-32/plants/fasta/zea_mays/dna/. The sequence Alignment generated by Tophat2 was then processed by the software package Cufflinks to assemble the Sequence Alignment/Map file into transcript fragments. FPKM [29] was used as the unit of measurement to estimate transcript abundance.

Differential expression analysis
Differential expression analysis was performed using the DESeq R package (1.10.1) [30]. DESeq provides statistical routines for determining differential expression in RNA-Seq gene expression data using a model based on a negative binomial distribution. The resulting P-values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate (FDR). Genes with an adjusted P-value < 0.05 and log 2 (fold change) > 1 (according to DESeq) were assigned as differentially expressed.
Gene Ontology (GO) enrichment analysis of the DEGs was implemented using the GOseq R packages based on the Wallenius non-central hyper-geometric distribution [31], which can control for gene length bias in DEGs. KOBAS [32] was used to test the statistical enrichment of DEGs in KEGG pathways.

Weighted gene co-expression network analysis
Weighted Gene Co-expression Network Analysis (WGCNA) is a statistical tool for clustering transcripts that have a similar expression pattern across a group of samples [33,34]. WGCNA was originally developed to analyze microarray data but can also be used to examine RNA-Seq data [35]. The input data for the WGCNA were the normalized values for each transcript. First, all available samples from each time point in 'Yayu889' samples were collected to identify modules that had different expression patterns. Next, a soft threshold was chosen to create networks with a scale-free topology using the method described by Langfelder and Horvath [33]. After the networks were built, modules of transcripts with similar expression patterns were created and eigengenes for these modules were calculated. Finally, correlations between these eigengenes and the factor of interest (at 81 dap, 89 dap, 91 dap, and 93 dap) were calculated.

Gene validation and expression analysis
To validate RNA-seq results, qPCR was performed. RNA samples were reverse-transcribed into cDNA using the PrimeScript® RT reagent Kit (Takara Code RR047B; Takara, Shiga, Japan). Expression profiles of genes were examined in triplicate using SYBR® Premix Ex TaqTM II (Tli RNaseH Plus; Takara Code RR820A) in LightCycler 480 (Roche Applied Science, Switzerland) following the 20-μL Real Time system, including 10 of μL SYBR® Premix Ex TaqTM II (2X), 0.6 μL of forward and reverse primer (6 μM) with a final concentration of 0.3 μM, 1.0 μL of cDNA (50 ng/μL), and 7.8 of μL sterile distilled water. Two-step PCR was performed according to the manufacturer's procedure, and the initial denaturation step was conducted at 95°C for 30 s, followed by 40 cycles of 95°C for 5 s and 60°C for 30 s, and then 72°C for 60 s. All primers are listed in Additional file 1.

Protein precipitation and trypsin digestion
Proteins dissolved in solution were precipitated by adding samples to three sample volumes of pre-cooled acetone (320,110, Sigma). White protein precipitates were observed upon mixing. To maximize the efficiency of protein precipitation, each sample was incubated at − 20°C overnight (approximately 12 h). Precipitated proteins were then collected by centrifugation with a benchtop centrifuge at max speed (approximately 14 k) for 15 mins. The protein pellet was washed twice with equal volumes of pre-cooled acetone. For the washing step, a shorter centrifugation time of 5 mins was used. Residual acetone was removed by drying in a biosafety cabinet or light blowing with air.
Proteins were re-suspended in 8 M urea and reduced with 20 mM DTT at 60°C for 1 h. Proteins were then alkylated with 40 mM IAA at room temperature for 30 min in darkness. Alkylation reactions were quenched with 10 mM DTT. Samples were diluted to 2 M urea concentrations with HPLC grade water. Protein concentrations of the samples were determined with a modified Lowry's assay (DC Assay Kit, Cat. 500-0111, Bio-Rad, Hercules, CA, USA) with the use of bovine serum albumin (BSA) to construct a calibration curve. An appropriate amount of trypsin was then added to the samples to reach an enzyme-to-substrate ratio of 1:100. Digestion was performed with 100 mM triethylammonium bicarbonate (pH 8, T7408, Sigma-Aldrich, St. Louis, MO, USA) at 37°C for 18 h. Digested proteins were desalted with ZipTip (Cat. ZTC18S960, Millipore, Billerica, MA, USA) for LC-MS/MS analysis.
Protein quantification using the multiple reaction monitoring (MRM) method Separation of protein digests was performed on an Eksigent NanoLC-Ultra® 2D System and cHiPLC® system (Eksigent, Dublin, CA, USA) in serial column mode. For each injection, the sample was desalted on a 75 μm × 15 mm analytical column, then eluted onto a second analytical column to create a 30-cm column length for separation. Both column chips were filled with ChromXP™ C18-CL 3 μm 120 Å phase (Eksigent, Dublin, CA, USA). Peptides were separated using a linear gradient formed by A (2% ACN, 0.1% FA) and B (98% ACN, 0.1% FA), from 12 to 32% of B over 60 min at a flow rate of 250 nL/min. Each injection was performed using a full loop injection with a 1-μL sample loop. The MS analysis was performed on a TripleTOF® 5600 system (AB Sciex, Framingham, MA, USA) using both the MRM HR workflow and Scheduled MRM HR workflow. High resolution mode (> 30,000) was used to capture MS spectra with the accumulation time being 250 ms per spectra. We performed full scan MS/MS under high-sensitivity mode, and optimized the accumulation time per-cycle. Rolling collision energy was used to set collision energy that spread of 5 V. The retention time windows used for datasets were 2.5 min.
ProteinPilot (Version 4.5.0.0 by AB Sciex) was used for sample profiles with the Paragon method to source possible peptide matches from the UniProt protein database (release 2013_08). The following search parameters in the Paragon algorithm in ProteinPilot were used: Sample type, identification; Cys alkylation, Iodoacetamide; Digestion, trypsin; Instrument, TripleTOF 5600; Special Factors, none; and Search Effort, Thorough ID. FDR analysis in the ProteinPilot software was performed, and FDR < 1% was used for the protein identification threshold. The candidate peptides were verified by using PeptideAtlas (www.peptideatlas.org), and the XICs for specific m/z values of different peptides were viewed by PeakView (version 1.1.1.2 by AB Sciex). The MS/MS spectrum detected by the mass spectrometer and the theoretical fragmentation spectrum for the peptides were manually checked to confirm the presence of the peptide in the XIC profile. Multiquant (version 2.0.2 by AB Sciex) was used to calculate peak area integration for the top three transitions selected from the MS/MS spectra.

Validation of SNP markers in candidate proteins
The SAMtools software package was used to simultaneously call SNPs across all samples [36]. The results were filtered to discard SNPs with quality scores below 70. To ensure reliability of these SNPs, those identified in at least eight replicates (out of 12 libraries) for a cultivar that had a total read depth of at least six were assigned as SNPs for that cultivar. The identified SNPs from RNA-Seq were selected for PCR amplification and Sanger sequencing. Flanking sequences of selected SNPs were extracted from the reference genome and PCR primers were designed with Premier Primer 5 (Premier Biosoft, Palo Alto, CA, USA). All primers are listed in Additional file 1.

SA content determination assay
Leaf tissues with fresh weights of 100 mg were homogenized and extracted in 900 μL of 0.05 M phosphate saline buffer (pH, 7.4). The homogenate was centrifuged at 1 × 10 4 rpm at 4°C for 20 min, and the supernatant was collected. Then, the levels of salicylic acid (SA) were assessed using enzyme linked immunosorbent assay (ELISA) test kits (Hermes Criterion Biotechnology, Vancouver, Canada) according to the manufacturer's instructions.

Peroxidase enzyme activity assay
The crude enzyme extraction method was the same as those outlined above. Peroxidase (POD) activity was determined using the method of Upadhyaya et al. [37] with minor modification. The reaction mixture contained 3.9 mL of 0.3% guaiacol, 50 μL of enzyme extract, and 50 of μL 0.3% H 2 O 2 . The absorbance was monitored at 470 nm for at least 5 min at 1-min intervals; an absorbance change of 0.01 represented one unit of POD activity.

Carotenoid content determination assay
The total carotenoids were determined by the method of Yang et al. [38] with slight modification. Leaf tissues with fresh weights of 30 mg were ground and extracted in 1 mL of an acetone-water mixture (4:1). Total carotenoids were also calculated using the Yang et al. method [38].

Statistical analysis
Data obtained in the experiments were used to evaluate the disease index and MRM were subjected to an analysis of variance (ANOVA) using the Statistical Package for Social Science (SPSS; SPSS Inc., Chicago, IL, USA) version 17.0. The statistical significance was judged at a threshold of P < 0.05.

Evaluation of maize cultivars against GLS
'Yayu889' and 'Zhenghong532' were cultivated in 2012 and 2011, respectively, in southwestern China. They showed similar growth rates and developmental progressions throughout the sampling period. The growing cycle from sowing to tassel lasted 98 days for 'Yayu889' and 100 days for 'Zhenghong532.' The growing cycle from sowing to maturity lasted 170 days for 'Yayu889' (plant height, 280 cm; ear length, 19.8 cm) and 163 days for 'Zhenghong532' (plant height, 240 cm; ear length, 19.5 cm), while different resistance levels were exhibited in response to GLS. The development of disease symptoms of the two maize cultivars was investigated at 81, 89, 91, and 93 dap, respectively, and GLS symptoms were scored at 123 dap.
Of the two cultivars evaluated against GLS, 'Yayu889' was resistant (disease index, 14.9%), while 'Zhenghong532' was susceptible (disease index, 78.4%; Fig. 1). We first observed small chlorotic spots on the bottom leaves at 89 dap. Then, spots developed into small tan spots with chlorotic borders at 91 dap. At 93 dap, mature lesions were formed, which were gray to tan in color and distinctly rectangular in shape. According to the observed phenotypic change, we sampled tissues at three time points (89, 91, and 93 dap) with typical lesions and at the beginning time point (81 dap) with no disease symptoms to investigate the differential transcript changes after challenging 'Yayu889' and 'Zhenghong532' with GLS, respectively. Three biological replicates were sequenced for each genotype and time point. The beginning time point (81 dap) was used as control group for profiling the induced transcription changes under GLS in both cultivars.

Differentially expressed genes
A total of 628.95 million reads were generated via 150-bp paired-end sequencing from 24 cDNA libraries. The raw data were filtered with the FASTQ_Quality_Filter tool from the FASTX-toolkit (Additional file 2). The clean data were used for further analysis. After preprocessing the RNA-Seq data, the reads were mapped to the maize reference genome version 4 (B73 RefGen_v4) using Tophat2 [28] (Additional file 3). Differential expression analysis was performed using the DESeq R package (1.10.1).
By comparing the RNA-Seq expression profiles of each pair of three biological replicates (Additional file 4), the correlation coefficients were evaluated. The results showed that three biological replicates in each group had high consistency.
Functional annotations of the DEGs available from seven public databases were used to investigate the functions and roles in response to GLS. Several categories of defense-related genes, which were already reported in plant immunity to necrotrophic fungi infection [39][40][41][42][43][44], are described in the following sections.

Quantitative RT-PCR validation
To validate the RNA-Seq data, 12 DEGs were randomly selected for qRT-PCR assays based on their expression patterns at three time points. The results of the selected DEGs showed that the qRT-PCR results were consistent with the RNA-Seq results, as both RNA-Seq and qRT-PCR analyses showed similar expression patterns of up-and down-regulation, indicating that the RNA-seq analysis was well suited for analysis of GLS-induced maize leaf transcriptomes (Fig. 4, Additional file 1).

GO and KEGG enrichment analyses
GO category enrichment analysis was applied to elucidate the functional enrichment of the DEGs and Fisher's exact test when a corrected P-value < 0.05 was used. A total of 22 GO terms were enriched in the resistant 'Yayu889' cultivar, of which 11 common GO terms were significantly enriched at 91 and 93 dap (Fig. 5a). Several of the 22 common GO terms in 'Yayu889' were associated with general plant responses to fungal attacks, including 'response to salicylic acid' (GO:0009751), 'protein phosphorylation' (GO:0006468), 'transmembrane receptor protein serine/threonine kinase signaling pathway' (GO:0007178), and 'chitin binding' (GO:0008061). In the susceptible 'Zhengong532' cultivar, there were 31 significantly enriched GO terms, while only two common GO terms were enriched at 91 and 93 dap (Fig. 5b). Although 'Zhenghong532' displayed a higher number of GO terms, many of them were not associated with a specific plant response to biotic stress. Furthermore, a different battery of response mechanisms to GLS was activated by 'Zhenghong532' compared to 'Yayu889, ' e.g., 'oxidation-reduction process' (GO:0055114) and 'carotenoid biosynthetic process' (GO:0016117).
In order to compare and summarize the response of the two cultivars to GLS, the DEGs were mapped to the reference canonical pathways in KEGG. Based on the plant-pathogen interaction map, most DEGs were involved in the hypersensitive response and Ca 2+ signaling pathways. 'Yayu889' expressed more Ca 2+ -dependent protein kinase (CDPK) genes and CML genes from these two pathways compared to 'Zhenghong532' (Additional file 11). As summarized in the plant-pathogen interaction pathway, four CML genes (Zm00001d005766, Zm00001d042056, Zm00001d029028, and Zm00001d0

Weighted gene co-expression network analysis (WGCNA)
To determine whether coordinated transcriptional responses to GLS were correlated with resistance in 'Yayu889, ' the datasets of 'Yayu889' expression across individual days (81, 89, 91, and 93 dap) were subjected to WGCNA. A total of 5046 transcripts were assigned to 13 co-expression modules, named after randomly assigned colors (Fig. 6a). Six of them were significantly relevant to GLS disease responses (Fig. 6b) as the eigengene was associated with the development of necrotic lesions.
Here, three modules (DG, DM, and PU) were listed for further analysis. We applied GO enrichment analysis to investigate the functional enrichment of transcripts in these three modules, respectively. In module DG, several biological processes related to GLS disease were enriched, including 'protein phosphorylation, ' 'response to chitin, ' 'hormone-mediated signaling pathway, ' and 'ionotropic glutamate receptor signaling' (Additional file 12). Several genes belonging to the DG module, including two MAPK (Zm00001d048027 and Zm00001d047758) and LYK5 (Zm00001d053695), were up-regulated in 'Yayu889' and expressed at the baseline level in 'Zhenghong532' (Fig. 7, Additional file 13). In the DM module, both 'protein phosphorylation' and 'cellular response to hormone stimulus' GO terms were overrepresented (Additional file 12). Two leucine-rich repeat receptor-like protein kinase genes (LRR-RLK; Zm00001d041476 and Zm00001d038522) of this module were significantly up-regulated at 91 and 93 dap compared to 81 dap, while baseline expression of these genes was observed at these two stages in 'Zhenghong532' (Fig. 7, Additional file 13). The PU module was significantly enriched for 'vitamin E biosynthetic process, ' 'hormone biosynthetic process, ' 'protein phosphorylation, ' and 'L-phenylalanine catabolic process' (Additional file 12).
Next, multiple pathways involved in plant immunity, including 'response to salicylic acid, ' 'protein phosphorylation, ' 'response to chitin, ' 'oxidation-reduction process, ' and 'carotenoid biosynthetic process, ' which were common or specifically enriched by differentially expressed analysis and WGCNA, were more relevant to the response to GLS in maize and thus summarized for further analysis.

MRM quantitation of candidate proteins
GLS resistance pathways were identified by transcriptome profiling and data from the published literature, revealing a total of 56 DEGs that were associated with disease resistance (Fig. 7, Additional file 13). To develop our analysis further, we used MRM to quantify the expression of candidate proteins. As shown in Table 1, a total of 12 proteins involved in 'response to salicylic acid, ' 'protein phosphorylation, ' and 'response to chitin' were quantified. In 'Yayu889, ' two of these proteins (CML25 and CML42) were significantly up-regulated at 93 dap compared to 81 dap, while baseline expression was viewed at all three stages in 'Zhenghong532, ' confirming the role of CML in GLS responses. The remaining 10 proteins showed a trend of up-regulation in both genotypes, but the expression level of 'Yayu889' was significantly higher than that in 'Zhenghong532' at all four time points, including two CML proteins, two LRR-RLK proteins, RING-H2 finger protein ATL2A, nematode resistance protein-like HSPRO2, putative serine/threonine-protein kinase-like protein CCR3, protein LYK5, MAP kinase kinase (MKK), and glutathione-S-transferase (GST).

Candidate gene identification
Co-localization of QTLs with DEGs facilitates the screening of candidate genes for GLS resistance, and the impact of  Table 2, among the 12 candidate proteins, four (Zm00001d038522, Zm00001d038061, Zm00001d036951, and Zm00001d035467) were mapped within the GLS QTLs that were previously reported [19,20,22,23]. Interestingly, all four proteins were located on chromosome 6. Then, we screened SNPs in these proteins by using the RNA-Seq data and validated them by PCR assay. As shown in Table 2, four SNPs were identified within two proteins, a non-synonymous SNP (SNP4-A) located in the exons of CML25 (Zm00001d038061) and unique to 'Yayu889, ' while the LRR-RLK (Zm00001d038522) protein carries two unique non-synonymous SNPs (SNP1-G, SNP2-A) and a unique synonymous (SNP3-C) in 'Zhenghong532.' Considering that they were co-localized with resistance QTLs, these two proteins are highly likely resistance proteins for GLS disease. Moreover, another two proteins that mapped to a resistance QTL region, i.e., GST (Zm00001d036951) and CCR3 (Zm00001d035467), are also noteworthy.

Salicylic acid, peroxidases, and carotenoids are induced by GLS
In our study, the GO term 'response to salicylic acid' (GO:0009751) was specifically enriched in 'Yayu889' at 91 and 93 dap, and several DEGs, including two MYB genes, alpha-dioxygenase 1, and a HSPRO2 gene, were identified; 'Yayu889' expressed a putative MYB gene  Expressions with the same lower-case letter are not significantly different (P > 0.05) (Zm00001d051149) at all three phases and another MYB gene (Zm00001d017268) at 91 and 93 dap. However, in 'Zhenghong532, ' both of the genes only showed a weaker up-regulation at 91 dap (Fig. 4, Fig. 7, Additional file 13). Additionally, two HSPRO2 genes (Zm00001d012321 and Zm00001d042811) that belong to this GO term were up-regulated in 'Yayu889, ' but either down-regulated or expressed at baseline levels in 'Zhenghong532' (Fig. 4, Fig.  7, Additional file 13). Next, the SA content of the two maize genotypes was measured using ELISA. As shown in Fig. 8a, the SA content of 'Yayu889' was significantly increased at 91 and 93 dap compared to 81 dap, while no significant changes were observed in 'Zhenghong532.' The GO term 'oxidation-reduction process' (GO:0055114) was specifically found in 'Zhenghong532.' Protein SRG1 and four important redox enzymes of cytochrome P450, acyl-desaturase 2, geraniol 8-hydroxylase, and indolin-2-one monooxygenase were enriched. In order to verify this GO term from the metabolic level, the POD activity associated with this pathway was measured. As shown in Fig. 8b, the POD activity showed a trend of decreasing first and then increasing in both of the cultivars, but 'Yayu889' exhibited a greater decrease at the early infection stage compared to that of 'Zhenghong532.' The GO term 'carotenoid biosynthetic process' (GO:0016117) was overrepresented solely in 'Zhen-ghong532.' The total carotenoid content of the two cultivars was thus measured. As shown in Fig. 8c, the carotenoid level was significantly increased during early infection stages and had stabilized at 93 dap in 'Yayu889, '  These results indicate that SA, carotenoids, and POD may participate in the regulation of maize's defense response to GLS.

Discussion
GLS is a significant threat to maize cultivation in China owing to the current lack of resistant cultivars and the suitability of the climate of many regions to GLS. In this study, sets of genes were identified that respond to fungal infection in susceptible and resistant maize cultivars and define how each cultivar perceived the intruder and programmed its defense system. A limitation of this study is the genetic diversity between both maize lines. To reduce the impact of diverse genetic backgrounds on the results of the analysis, the datasets were filtered step by step. First, six major gene families, reported to be related to plant immunity [39][40][41][42][43][44], were found to be greatly expressed in the resistant genotype at all infection stages compared to the susceptible genotype (Fig. 3, Additional file 5). This indicates that the screened DEGs were mainly related to the disease resistance of maize. Then, bioinformatics analysis and public databases were utilized to enrich the metabolic pathways related to the plant-pathogen interactions and identify DEGs that belong to these pathways and that were potentially associated with GLS. At the same time, we further corroborated our results through metabolic and protein quantification experiments. Finally, to further narrow the screening range of resistance genes, candidate DEGs were co-localized with the resistance QTL loci for previously reported GLS [25].
Three expression modules associated with GLS responses were observed by WGCNA. The RING-H2 finger protein ATL2A, which belongs to module DG (Table 1, Additional file 12), was an early pathogen-associated molecular pattern (PAMP)-responsive gene and was reported to participate in the E3 ubiquitination degradation pathway in Arabidopsis [45,46]. Another protein (LYK5) in module DG, (Table 1, Additional file 12), was reported to be a chitin receptor that could induce plant innate immunity in Arabidopsis [47,48] by regulating the production of ROS and MAPK phosphorylation by interacting with the E3 ubiquitin ligase PUB13 [49]. We hypothesized that in 'Yayu889, ' the PAMP-triggered immunity (PTI) system was activated after being challenged by GLS infection; once the pathogen invaded, LYK5 and ATL2A were activated followed by an induction of the E3 ubiquitination degradation pathway. In addition, the L-phenylalanine catabolic process in module PU, which contains four phenylalanine ammonia-lyase (PAL) proteins, was reported to participate in the biosynthesis of defense-related phytohormone SA and phytoalexin in rice [15]. Additionally, a previous study found that SA could induce the resistance of maize to GLS by enhancing PAL activity [50].
SA plays an important role in resistance and defense induction in response to pathogen attack [51]. In this study, two MYB genes (Zm00001d051149 and Zm00001d017268) were identified (Fig. 7, Additional file 13), which were reported to regulate the SA-dependent signaling pathway [52] and modulate antagonistic interactions by activating SA-mediated defenses and repressing JA-mediated defenses [53]. Additionally, HSPRO2 (Zm00001d012321) was specifically activated in 'Yayu889' (Fig. 7), which was consistent with a report that HSPRO2 positively induced the basal resistance in Arabidopsis against Pseudomonas syringae through the downstream function of SA [34]. A resistance induction test for maize exposed to GLS was previously reported, and a 30.89% induced resistance was produced by spraying SA on the leaves of maize plants [50]. Similar results were produced in our study, indicating that only the resistant cultivar 'Yayu889' exhibited induced SA expression (Fig. 8a). In summary, these results indicate that through modulation of the production of SA, 'Yayu889' regulated the level of tolerance to GLS disease, possibly through interactions of the SA signaling molecules with ROS. This was also reported under other stress conditions in other plants [54] and it was confirmed in our work by induction of POD activity in 'Yayu889.' A change in cellular reduction-oxidation status is one of the earliest responses detected in challenged cells [55]; Several redox enzymes were involved in this pathway. For example, maize cytochrome P450 was reported to produce phenylacetaldoxime and indole-3-acetaldoxime in heterologous systems and might contribute to plant defense [56]. Similarly, acyl-[acyl-carrier-protein] desaturase may enhance the resistance of Arabidopsis and soybean to multiple pathogens [57]. Additionally, GST belongs to the GSH redox system that suppresses cercosporin [58]. Lastly, production of reactive oxygen species (ROS) occurs rapidly in response to an attempted pathogen invasion of potential host plants [59]. Based on the analysis of POD activity in the present study (Fig. 8b), we hypothesized that 'Yayu889' induced ROS levels by inhibiting POD activity at early stages of GLS infection, such as H 2 O 2 production and then, the inhibition of H 2 O 2 levels by enhancing POD activity after the infection was established. This hypothesis corroborated the concept that ROS plays a role in resistance during early infection while promoting disease during the later stages of infection [8,60].
Protein phosphorylation has been regarded as essential to plant immunity [61]. Genes within this group include various protein kinases, such as receptor kinases and MAPKs. In this study, two serine/threonine protein kinases were identified, which were reported to participate in a phytohormone-activated signaling pathway, and the regulation of plant-type hypersensitive response in Arabidopsis [62]. 'Yayu889' is known to encode two LRR-RLK genes (Zm00001d041476 and Zm00001d038 522; Additional file 13) and two surface-localized pattern recognition receptors (PRRs). These genes were reported to constitute an important layer of innate immunity in plants. One of the genes (Zm00001d038522) contained a SNP within a major QTL region (Table 2); thus, it may be a candidate resistance gene with high potential to be applied in MAS to produce GLS-resistant genetic materials [63]. Additionally, a MKK7 protein (Zm00001d048027) in this pathway was consistently expressed at higher levels ( Table 1, Additional file 13) in 'Yayu889, ' and previous studies have found that ectopic expression of MKK7 could elevate SA levels and enhance basal resistance to Arabidopsis pathogens [64].
Ca 2+ ions are recognized as a crucial second messenger in signaling pathways, combining the perception of environmental stimuli with plant adaptive responses [65]. This conversion into biological responses requires the presence of calcium sensor relays such as calmodulin (CaM) and CML proteins [66,67]. CML9 and CML8 were found to act as positive regulators of defense mechanisms against microbial pathogens in Arabidopsis [68,69]. In the present study, protein expression of CML25 and CML42 was only induced in 'Yayu889' (Table 1). Notably, CML25 was co-localized with a major QTL region and the SNP4-A may be the reason for the difference in protein expression among the two maize cultivars ( Table 2). This indicates that the Ca 2 + -dependent signaling pathway may participate in 'Yayu889' resistance to GLS by activation of SA and ROS, and this was verified in plant immunity to other pathogens [68,70,71].
Cercosporin, the non-host-selective toxin produced by Cercospora zeina, is a photo-activated perylenequinone that converts molecular oxygen to active oxygen species, including singlet oxygen [72]. Plant carotenoids play a key role in quenching singlet oxygen molecules and are hypothesized to play a role in defense against toxins that produce reactive oxygen species [73,74]. The defense role of carotenoids against GLS in maize was confirmed by experiments of total carotenoid quantitation for the two maize genotypes in this study (Fig. 8c).
The use of molecular markers and MAS can facilitate gene stacking for more durable resistance, since multiple genes that are effective against the pathogens can be combined into a single breeding line or a variety in a manner that would not be possible with phenotypic selection only [75]. We have narrowed the potential protein candidates by co-localizing with previously mapped resistance QTL regions and screened the SNPs in each protein ( Table 2). This method can promote the screening of potential resistance genes or molecular markers for GLS. In this study, the four proteins located in the resistance QTL regions (i.e., CML25, LRR-RLK, GST, and CCR3) were tightly associated with GLS resistance. QTL mapping of a maize recombinant inbred line (RIL) population derived from subtropical CML444 × SC Malawi maize inbred lines found that the QTL region in chromosome 6 (located in bin 6.06/6.07) explained more than 11% of phenotypic variation [22], CML25 and LRR-RLK, which also contains SNPs, were located in this region, indicating that they are most likely candidate resistance proteins. Additionally, SNP1-G, SNP2-A, and SNP4-A are potential molecular markers for MAS. However, the functions of these proteins or markers and their contribution to the enhancement of resistance to GLS need to be further verified. Moreover, GST (located in bin 6.04) and CCR3 (located in bin 6.01) mapped to other QTLs on chromosome 6 and are candidate proteins, because this region (qGLS6.01-6.04) was associated with 6.85% of variation in GLS resistance [20]. The remaining eight proteins were not located in any GLS-resistance QTL regions, which may have little effect on phenotypic variation. However, we failed to develop markers associated with these proteins by using RNA-Seq owing to the limited number of sampled populations; thus, further analysis could use genomewide association prediction models to establish molecular markers associated with these proteins. Since we have developed markers for all relevant proteins, genomic selection is a better approach for the capture of little-effect genes than MAS. The reason for this is that the practical implementation of MAS for stacking multiple disease resistance genes is difficult; therefore, the population size needed for the maintenance and fixing of multiple genes quickly exceeds the reasonably available resources for MAS.

Conclusions
To develop resistant cultivars and control GLS, understanding host-pathogen interactions is essential. In this study, novel GLS resistance genes and metabolites were identified by combining transcriptional, protein, and metabolic expression levels, and an interaction network of maize resistance to GLS was also constructed, which was regulated by ROS, SA, and Ca 2+ signaling pathways as well as carotenoids. In particular, four genes on chromosome 6 that were consistently expressed at higher levels in the resistant cultivars and occurring within GLS resistance QTL regions seemed most promising, namely a CML25, a LRR-RLK, a GST, and a CCR3 protein. Further validation of the association of these proteins is necessary to improve our understanding of maize resistance to GLS. This information will provide new insights for molecular marker-assisted breeding of GLS-resistant maize cultivars.