Transcriptome Analysis Identied Core Genes Involved in Maize Resistance to Rhizoctonia Solani

Abstract


Background
Maize (Zea mays L.) is a major cereal crop grown throughout the world for food, feed and fuel. However, maize production is frequently affected by biotic and abiotic stresses with devastating effects on yield and product quality. Among these stresses, banded leaf and sheath blight (BLSB) caused by Rhizoctonia solani is considered to be one of the emerging and serious diseases limiting maize production under climate change scenarios [1]. The dominant anastomosis group (AG) on maize is AG1-IA [2]. This pathogen usually causes necrotic spot on the leaf, sheath, and stem. Under high temperature and humidity conditions, it colonizes almost all the aboveground parts of the plants, including cobs and ears.
Control of R. solani is a di cult problem in maize, and the traditional control methods are to use agricultural prevention measures and chemical fungicides to reduce the damage caused by R. solani. With the improvement of the awareness of food safety and environmental protection and requirements for high quality maize products, it is becoming increasingly important to cultivate disease-resistant varieties. However, most of the existing maize varieties are susceptible to R. solani, and there are few maize germplasms with high resistance, which greatly limits the study of resistance inheritance [3].
BLSB resistance is a quantitative trait regulated by multiple genes. To our knowledge, no major resistant gene has been mapped or cloned in maize. Up to now, many quantitative trait loci (QTLs) associated with BLSB resistance have been mapped on the maize genome. For example, 11 QTLs for resistance to BLSB were identi ed using the F 2 individuals derived from the cross of inbred lines R15 and 478 [3], while four BLSB resistance QTLs were mapped using the F 2 individuals derived from the cross of inbred lines CML429 and DM9 [4]. However, very few candidate genes were cloned in these QTLs.
Global gene expression analysis can not only be used to screen differentially expressed genes (DEGs), but also an important means to explore the molecular mechanisms of plant-pathogen interaction. RNAsequencing (RNA-seq) technology provides a more accurate measurement of transcription levels to identify major changes in gene expression, and numerous DEGs have been identi ed in different plantpathogen interactions using this approach [19][20][21]. RNA-seq also has been successfully used in investigating transcriptional changes of plants during interaction with R. solani [22][23][24]. However, all these studies focused on the expression changes of host genes at the early stage of R. solani penetration. Analysis of maize transcriptional changes during R. solani expansion has not been reported yet. Because there is no major resistant gene conferring BLSB resistance, genes involved in restricting R. solani expansion also have high application value in resistance improvement.
In this study, gene expression pro le of maize infected by different virulent strains of R. solani for 3 and 5 d was investigated by RNA-sEq. A total of 388 core genes involved in maize resistance to R. solani were identi ed. Among these genes, differentially expressed transcription factors (TFs), hormone-related genes and pathogenesis-related (PR) genes were collected and analyzed. Furthermore, we identi ed and validated two novel genes involved in R. solani resistance. Taken together, these results provide a comprehensive insight into the resistant mechanisms of maize against R. solani expansion and potential gene resources for the improvement of BLSB resistance.

Results
Comparative analysis of B73 infected by different virulent R. solani strains To compare the different responses of maize to R. solani, the maize inbred line B73 was inoculated with high virulent strain (HVS) and low virulent strain (LVS) of R. solani, respectively. At 5 d post inoculation (dpi), B73 exhibited longer lesion (5.57 ± 0.25 cm) infected by HVS than by LVS (3.97 ± 0.34 cm) (Fig.   1a). Because R. solani was the same pathogen that caused BLSB on maize and sheath blight (ShB) on rice, and it would be much easier and faster to do the gene function research using rice transformation than maize transformation. The pathogenicity of HVS and LVS on the rice cultivar Zhonghua 11 (ZH11) was also investigated. After inoculation for 5 d, HVS caused longer lesion (5.03 ± 0.31 cm) than LVS (3.85 ± 0.23 cm) (Fig. 1b). These results indicated that the two strains showed the same difference in pathogenicity on maize and rice. Furthermore, the biomass of R. solani was measured by histological staining and qPCR. No signi cant difference in hyphae amount was observed at 1 dpi. At 3 and 5 dpi, sheaths infected by HVS showed more hyphae than by LVS (Fig. 1c, d). Because there were no major resistance genes for BLSB, we were more interested in genes associated with restriction of hyphae expanding instead of penetration. Time points of 3 and 5 dpi were selected for the transcriptome analysis.

Identi cation of DEGs in B73 inoculated with LVS and HVS
To explore the transcriptional changes of maize under R. solani attack, sheaths of B73 plants were inoculated with LVS and HVS, respectively. Un-inoculated sheaths were used as control. For each R. solani strain, the sheath samples from 3 and 5 dpi were pooled together as one biological replicate. Finally, the sheath samples with three biological replicates were collected for RNA-seq. Nine libraries were sequenced, and more than 54 million raw reads were obtained with an average of six million in each library. The percentage of clean reads ranged from 97.96% to 98.76%. More than 78% of clean reads were mapped to the reference genome (Table S1). The raw data was deposited in the National Center for Biotechnology Information Sequence Read Archive (NCBI SRA) under the accession number PRJNA739273.
Next, DEGs were identi ed in LVS-and HVS-infected B73 compared to un-inoculated B73 with the threshold of |fold change| ≥ 2 and FDR ≤ 0.001. In LVS-and HVS-infected B73, 1545 and 1166 DEGs were up-regulated, and 1470 and 462 were down-regulated, respectively (Fig. 2a). More DEGs were regulated by LVS infection than by HVS. Venn diagram showed that 49.8% of the up-regulated DEGs were speci cally expressed in LVS_DEGs, while 72.8% of the down-regulated DEGs were speci cally expressed in LVS_DEGs. In HVS_DEGs, only 33.5% and 13.2% of up-and down-regulated DEGs were speci cally expressed, respectively (Fig. 2b).
Interestingly, GO term "response to stimulus" has the most DEGs both in LVS_DEGs and HVS_DEGs, and only defense-related GO terms, such as "response to stimulus", "defense response", "cell wall organization", "response to stress", "response to biotic stimulus" and "response to hormone", were commonly enriched in LVS_DEGs and HVS_DEGs. So we speculated that a set of common DEGs might constitute a major part of the immune transcriptional reprogramming engine.
Identi cation of core defense-related genes during maize and R. solani interaction To prove the above hypothesis, we focused on the 625 DEGs commonly expressed both in LVS_DEGs and HVS_DEGs (Fig. 2b, Table S3). Among these DEGs, 388 genes were up-regulated and 237 were downregulated (Fig. 2b). Meanwhile, hierarchical clustering of the fold change of the 625 DEGs showed that they were divided into two groups: genes in the larger group were up-regulated and genes in the small group were down-regulated both by LVS and HVS (Fig. 4a). Then, the expression patterns of several DEGs were veri ed by quantitative Real-Time PCR (qRT-PCR). As shown in Fig. S1, the expression patterns of 10 DEGs (including 5 up-regulated and 5 down-regulated) were consistent with the RNA-seq analysis, indicating that the RNA-seq results were reliable. GO (biological process) enrichment of the 388 upregulated genes and 237 down-regulated genes was performed, respectively. Defense-related GO terms were signi cantly enriched in the up-regulated genes, including "defense response", "salicylic acid biosynthesis", "phenol biosynthesis" and "response to stress" (Fig. 4b, Table S4), suggesting that the 388 up-regulated genes were important components in maize resistance to R. solani. However, no defenserelated GO terms were signi cantly enriched in the down-regulated genes (Fig. 4c, Table S4). Taken together, these results indicated that most of the 388 up-regulated genes might be associated with maize immunity and are potential candidates for further functional studies. Therefore, the 388 genes were considered to be the core of maize immune genes.
Identi cation of differentially expressed TFs in the core gene list TFs usually play positive or negative roles in plant defense response. In the core gene list, 26 genes encoded TFs, including six ZF, ve bHLH, two WRKY, two bZIP, two MYB, two AP2/EREBP, two HSF, one NAC, one TCP, one JMJ, one VQ and one GRAS (Fig. 5a, Table S5). These TFs were further classi ed into two clusters in a hierarchical clustering based on the fold change. Most of the TFs (20 TFs) were moderately induced, while the minority TFs (6 TFs) were highly induced by R. solani (Fig. 5b). Four of them have been reported to be involved in maize defense response. For example, ZmbZIP111 was upregulated during the induced systemic resistance activation by Trichoderma atroviride [25], while ZmNAC41 was transcriptionally induced by the hemibiotrophic ascomycete fungus Colletotrichum graminicola [26]. ZmJMJ22 was identi ed to confer resistance to maize rough dwarf disease [27]. ZmbHLH124 played an important role in the maize resistance-speci c response to Puccinia sorghi [28].
These functionally characterized TFs implied that the core gene list contains important regulators associated with maize immunity.
TFs bind to the cis-elements in the promoter region to regulate target gene expression. Therefore, the TFbinding motifs in the promoters of the 388 core genes were investigated. Promoter analysis showed that binding sites of WRKY, bHLH, ZF, HSF, MYB, AP2/EREBP, bZIP, NAC and TCP TFs were found in the core genes (Table S6). Noticeably, the binding motifs for WRKY, bHLH, ZF, HSF, MYB and AP2/EREBP TFs were mostly identi ed among these core genes, and WRKY12, bHLH3, AZF1, HSFB2A, MYB105 and RAV1 were represented in these TFs (Fig. 5c, Table S6). In plants, these TFs play various roles in plant development and response to biotic and abiotic stresses. For example, overexpressing BrWRKY12 enhanced Arabidopsis and Chinese cabbage resistance to Pectobacterium carotovorum [29]. AtbHLH3 acts as a transcription repressor to regulate jasmonate (JA)-mediated defense and development [30]. AtAZF1 negatively regulates abscisic acid (ABA)-responsive and auxin-inducible genes under abiotic stress [31]. . These results suggested that these identi ed TFs might regulate maize resistance to R. solani by targeting the core genes.
Identi cation of hormone signaling-related and PR genes in response to R. solani infection Previous studies indicated that salicylic acid (SA), JA and ABA signaling pathways are involved in R. solani resistance [7,35,36], so hormone signaling-related DEGs were examined in the core gene list. A total of 28 genes were identi ed, including six SA-related, six JA-related and 16 ABA-related genes (Table S7).
PR genes are markers in plant defense response and overexpressing some PR genes could enhance plants resistance to R. solani [37,38]. A total of 12 up-regulated PR genes were identi ed in the core gene list, including one ZmPR1 gene, two ZmPR2 genes, three ZmPR4 genes, one ZmPR8 gene, three ZmPR9 genes and two ZmPR10 genes (Fig. 6d, Table S8). In addition to these hormone-related and PR genes, 87 other genes have been characterized to function in plant defense response (Table S9).
ZmNAC41 and ZmBAK1 positively regulate rice resistance to R. solani Among the 388 core genes, the NAC transcription factor 41 (ZmNAC41, GRMZM2G439903) gene and BRASSINOSTEROID INSENSITIVE 1-associated receptor kinase 1 (ZmBAK1, GRMZM2G015933) gene exhibited very high transcription levels when infected by R. solani. To determine the role of ZmNAC41 and ZmBAK1 in R. solani resistance, we overexpressed these two genes in the rice cultivar ZH11 and obtained 23 ZmNAC41 overexpression (OE) and 25 ZmBAK1 OE lines, respectively. In T 1 generation, three OE lines of the two genes were selected and sheaths were inoculated with R. solani. At 5 dpi, the ZmNAC41 and ZmBAK1 OE lines were more resistant to R. solani than ZH11, with average 3.66 cm and 3.72 cm reduction in lesion length, respectively (Fig. 7a, b). Meanwhile, the expression of OsPR1a, OsPR1b and OsPR10 was also detected in OE lines post inoculation with R. solani. Compared to ZH11, these PR genes were more strongly induced in the OE lines (Fig. 7c, d). Taken together, these results indicated that ZmNAC41 and ZmBAK1 positively regulate rice resistance to R. solani.

Discussion
BLSB is an important maize disease that severely threats yield and quality. Although some DR genes have been identi ed in maize and other plants [7][8][9][10][11][12][13][14][15][16][17][18], these genes still cannot meet the needs of BLSB resistance breeding. Meanwhile, the defense mechanisms of maize against R. solani are not very clear. In this study, we performed RNA-seq to explore the transcriptional changes in maize at the R. solani expansion stage. Thousands of DEGs were identi ed under LVS and HVS infection. Noticeably, less DEGs were regulated by HVS (1166 up-regulated and 462 down-regulated) than by LVS (1545 up-regulated and 1470 down-regulated) (Fig. 2), suggesting that HVS might suppress partial resistance that antagonizes LVS. A set of DEGs commonly regulated by LVS and HVS was identi ed and most of them were involved in plant defense response. Therefore, these genes were considered core genes in maize-R. solani interaction.
Among the identi ed core genes, a total of 26 TFs encoding genes were induced by R. solani (Fig. 5a, Table S5). In addition to ZmbZIP111, ZmNAC41, ZmJMJ22 and ZmbHLH124, which have been functionally characterized in maize defense response [25][26][27][28], other types of TFs, such as ZF, WRKY, MYB, AP2/EREBP and VQ have been also reported to function in plant immunity to R. solani. For example, GhZFP1 enhances R. solani resistance by interacting with CZIRD21A and GZIPR5 [39]. OsWRKY80, OsWRKY4, BdWRKY38 and ZmWRKY79 positively regulates plants resistance to R. solani [8, 13,16]. Overexpressing OsMYB-R1 in rice showed endurance against R. solani [40]. A cell-wall-degrading enzyme-induced AP2/EREBP TF, OsAP2/ERF152, leads to resistance to R. solani by activating OsMPK3/6 [41], while another AP2/EREBP TF, NbOPBP1, confers rice resistance to R. solani [42]. VQ TFs usually interact with WRKY to regulate plant defense response [43], however, their roles in R. solani resistance have not been illustrated. Previous and our studies showed that VQ could be highly induced by R. solani [44], suggesting that VQ TFs may also have the R. solani resistant function. Taken together, these results demonstrated that the identi ed TFs in this study may have potential disease resistance to R. solani, and their functions should be further investigated.
SA and JA are typical phytohormones involved in plant defense response to pathogens, and previous studies showed that SA functioned in resistance to biotrophic pathogens, while JA was associated with resistance to necrotrophic pathogens [45]. Current evidence indicated that SA-and JA-mediated defense response is mutually antagonistic [46]. R. solani is a necrotrophic pathogen, while some ndings suggested a hemi-biotrophic nature of R. solani [36]. Therefore, the mechanisms of resistance against R. solani induced by SA and JA are more complex. In this study, 12 up-regulated genes of the core genes were identi ed to be related to SA and JA signaling pathways (Fig. 6a, b, Table S7). In SA signaling pathway, SA synthesis-related genes were mostly found. CBP60g and SARD1 are two master TFs in plant immunity [47]. During pathogen attack, SA biosynthetic genes such as ISOCHORISMATE SYNTHASE 1 (ICS1), ENHANCED DISEASE SUSCEPTIBILITY 5 (EDS5) and AVRPPHB SUSCEPTIBLE 3 (PBS3) were coordinately activated by CBP60g and SARD1 [48,49]. PAL is an important enzyme in SA synthesis [50] and it mediates plant defense response by regulating SA accumulation [51,52].
AAE activates JA biosynthetic precursors as one of the key JA biosynthetic enzymes [53]. In this study, ZmAAE11 was identi ed to be up-regulated by R. solani. In Medicago truncatula, knocking-down MtAAE3 increased the plant susceptibility to Sclerotinia sclerotiorum [54]. NPF proteins were originally identi ed as nitrate or di/tri-peptide transporters [55]. Recent studies revealed that this transporter family also can transport JA [56]. ZmNPF5.2 was dramatically induced by R. solani in the RNA-seq data. However, there is no direct evidence suggesting resistant functions of this transporter family. Additionally, JA-induced genes ZmKCS and ZmTPS2 were up-regulated post R. solani inoculation. KCS was reported to positively regulate abiotic stress tolerance and defense against herbivore wounding [57,58], while in M. truncatula, silencing MtTPS10 resulted in a susceptibility phenotype under Aphanomyces euteiches infection [59].
ABA is well known for its involvement in response to abiotic stress, and its functions in defense response have been also illustrated [60,61]. In this study, 16 ABA-related genes were identi ed among the core genes (Fig. 6c, Table S7). NCED is the key ABA biosynthetic enzyme that catalyzes the oxidative cleavage of cis-epoxycarotenoid [62]. Here, ZmNCED6 was identi ed to be up-regulated by R. solani. However, the roles of NCED in plant defense response have not been reported yet. Our previous studies showed that DR gene Os2H16 regulated R. solani resistance via ABA signaling pathway [7], and the TF OsASR2 bound to the GT-1 cis-element to regulate Os2H16 expression [9]. In this study, two TF genes, ZmASR6 and ZmASR10, were induced by R. solani, implying their potential involvement in R. solani resistance. Additionally, some ABA regulatory genes which have been functionally characterized in plant immunity, such as SnRK, SYR, PP2C and expansin, were also identi ed in the core genes.
PR proteins are synthesized in response to pathogen attack and limit the growth of pathogens. According to the molecular mass, isoelectric point, localization and biological activity, PR genes can be classi ed into 17 families [63]. In this study, 12 PR genes among the core genes were up-regulated post inoculation with R. solani, including one ZmPR1 gene, two ZmPR2 genes, three ZmPR4 genes, one ZmPR8 gene, three ZmPR9 genes and two ZmPR10 genes (Fig. 6d, Table S8).
PR1 genes are generally considered as marker genes for systemic acquired resistance (SAR) and SA plays an important role in the signal transduction of SAR [64]. Several PR1 promoters contain SAinducible cis-elements [65]. Previous studies showed that overexpressing PR1a increased tolerance to oomycete pathogens in transgenic tobacco [66], and some PR1 proteins had an antimicrobial activity against Phytophythora infestans [67], indicating the promising resistant function of PR1 genes to oomycetes. Here, one ZmPR1 gene was signi cantly up-regulated by R. solani, implying that it may result in resistance in maize against R. solani.
PR2 genes are predicted to encode beta-1,3-glucanase, which hydrolyzes beta-1,3-glycosidic bonds in beta-1,3-glucans, and have been proved to degrade the fungal cell walls [68]. The PR2 proteins may also degrade endogenous substrates to generate elicitors and induce defense response [63]. In this study, the transcription levels of two ZmPR2 genes were obviously increased challenged by R. solani, suggesting that these two genes may regulate R. solani resistance through the above pathways.
PR4 protein family consists of two classes of proteins. Class PR4 proteins contain a conserved cysteine-rich chitin-binding domain at the N-terminal, while Class PR4 proteins lack the chitin-binding domain [69]. Of the three identi ed ZmPR4 genes, two genes (GRMZM2G145461 and GRMZM2G358153) belong to Class and the other one (GRMZM2G465226) belongs to Class . PR4 genes are also considered as JA pathway markers [70]. Studies showed that knocking-out VvPR4b increased susceptibility to Plasmopara viticola [71] and OsPR4b protein exhibited in vitro antifungal activity against R. solani [72].
PR8 genes encode chitinases that have well-understood antifungal activity due to their hydrolytic action [73,74]. Plant chitinases are widely distributed and structurally heterogeneous enzymes, and they mainly function in plant-fungus interaction. In the process of preventing fungal pathogen infection, plants secret chitinases to digest chitin in the fungal cell wall, produce pathogen-associated molecular patterns and nally lead to immune recognition [75]. In this study, one ZmPR8 gene was signi cantly up-regulated by R. solani, suggesting its possible roles in blocking the invasion or expansion process by degrading the cell wall of R. solani.
PR9 proteins are reported to have peroxidase activity and are involved in lignin biosynthesis, which act as the cell wall barrier against pathogens [76]. Peroxidase is also involved in the oxidation of isophenols to forms that are more toxic to pathogens [77]. The two identi ed ZmPR9 genes in this study were enriched to the "cell wall organization" GO term (Table S4). Overexpressing rice peroxidase 114 gene enhance disease resistance to Alternaria radicina in transgenic carrot [78]. Hence, we predicted that up-regulation of the two ZmPR9 genes may limit the in ltration of R. solani toxins into cells and thereby the spread of infection.
PR10 proteins are classi ed as intracellular proteins and are closely related to a group of major food allergens, which belong to the Bet v 1 family [79]. Although several studies have reported that PR10 proteins have ribonuclease or antimicrobial activities [80, 81], their biological functions remain unclear. In this study, two ZmPR10 genes were induced after R. solani inoculation, and one ZmPR10 protein encoded by GRMZM2G112488 has been reported to possess broad spectrum antifungal activity against Botrytis cinerea, S. sclerotiorum, Fusarium oxysporum, Verticillium dahlia and Alternaria solani [82]. So their roles in resistance to R. solani need to be further con rmed.

Conclusions
We used RNA-seq to explore the transcriptional changes of maize at the expansion stage of different virulence strains of R. solani. Selection and GO enrichment analysis of DEGs identi ed a core gene list that played important roles in limiting the expansion of R. solani in maize tissues. The core gene list mainly included TF genes, PR genes, hormone-related genes and genes which have been characterized to function in plant defense response. Furthermore, the resistant function of two genes in the core gene list, ZmNAC41 and ZmBAK1, was investigated in transgenic rice plants.

Plant materials and inoculation
Maize inbred line B73 and rice cultivar ZH11 were grown in a growth chamber at 28 ± 2°C with a 16/8 h light/dark cycle. R. solani strains were cultured on potato-dextrose-agar medium at 25°C. For the sequencing library construction, sheaths of B73 were inoculated with LVS and HVS, and harvested at 3 and 5 dpi for RNA extraction. For expression analysis of PR genes, sheaths of ZH11 and transgenic rice plants were inoculated with R. solani and harvested at 24 h post inoculation for RNA extraction.

RNA extraction and mRNA library construction
Total RNA was isolated from sheaths inoculated with LVS and HVS at 3 dpi and 5 dpi using TRI reagent (Sigma-Aldrich, USA). The RNA of 3 dpi and 5 dpi was mixed as one biological replicate (three biological replicate, respectively). The mRNA was puri ed using oligo (dT)-attached magnetic beads and then fragmented. The rst-and second-strand cDNAs were synthesized and end repaired. The cDNA fragments with adaptors were ampli ed by PCR, and the products were puri ed by Ampure XP Beads. Quality and quantity of the library was assessed using the Agilent 2100 bioanalyzer. The quali ed library was ampli ed on cBot to generate the cluster on the owcell. The ampli ed owcell was sequenced single end on the HiSeq4000 platform (BGI-Shenzhen, China).
Validation of RNA-seq by qRT-PCR Ten DEGs were randomly selected and expression levels were validated by qRT-PCR. Primers were designed in qPrimerDB (https://biodb.swu.edu.cn/qprimerdb/) (Table S10). cDNA was synthesized with a ReverTra Ace qPCR RT Master Mix with gDNA Remover Kit (Toyobo Life Science, Japan). qRT-PCR was performed with a KOD SYBR qPCR Mix Kit (Toyobo Life Science, Japan) on a qTOWER3G Touch thermal cycler (Analytikjena, Germany). ZmActin was used as internal control to normalize the results. The relative expression was calculated using the 2 -ΔΔCt method.

Vector construction and rice transformation
For generating ZmNAC41 and ZmBAK1 overexpression vectors, the cDNAs of the two gene were cloned into the XcmI-digested pCXUN vector. Primers were listed in Table S10. The constructs were transformed into Agrobacterium tumefaciens strain EHA105. For rice transformation, embryonic callus derived from mature embryos was infected by EHA105 strain.

Disease assays
Sheaths of ZH11 and transgenic rice plants were inoculated with R. solani for 5 d. The disease symptoms were taken photographs and the lesion length was measured on R. solani-infected sheaths. Three independent T 1 transgenic lines and 10 plants for each line were assessed.

Statistical analysis
Means and SD values were calculated using Microsoft Excel 2019. The two-tailed student's t test was used to calculate the signi cance between samples. All statistical analysis was performed using SPSS software version 19.

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Availability of data and materials
The sequencing raw data for all libraries presented in this study are openly available in NCBI with the accession PRJNA739273 at https://www.ncbi.nlm.nih.gov/bioproject/PRJNA739273.

Competing interests
The authors declare that they have no competing interests.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.