Selection and determination of sampling time points
To determine the time points for transcriptome sampling, we used fugus cake cultured with R. solani AG1 IA for 2 days, inoculated the middle of the leaves of the different varieties at the tillering stage and wrapped the leaves with plastic wrap. The results showed that 24 h after inoculation, mycelia were attached to the leaf surface of Yanfeng47. In addition, 48 h after inoculation, disease spots began to form on the leaf surface of Yanfeng47, and Gangyuan8 also started to show disease spots. Seventy-two hours after inoculation, the lesion area of Yanfeng47 had expanded, and Gangyuan8 also showed obvious lesion. Ninety-six hours after inoculation, the leaves of Yanfeng47 had wilted, and the lesion area of Gangyuan8 was larger. Consistent with the results of previous studies, the resistance of Gangyuan8 and Yanfeng47 to R. solani AG1 IA showed significant differences (Supplementary Figure S1).
RNA-seq results
To investigate the changes in the gene expression profiles of Gangyuan8 and Yanfeng47 leaves during the initial stage of AG1 infection, we performed a transcription analysis of leaves after AG1 IA inoculation by high-throughput sequencing. Transcriptome data were generated from the leaves at different time points. A transcriptome analysis of each variety was performed using 4 time points and 3 biological replicates of each time point. Thirty samples yielded 236.04 Gb of clean data, and 6.18 Gb of clean data were obtained from each sample. The Q30 base percentage was at least 93.91%. The clean reads from each sample were sequenced using the designated reference genome, and the alignment efficiency ranged from 76.72–95.72%. To examine the quality of the biological replicates, we calculated Pearson correlation coefficient (PCC) values for each pair of samples and performed a cluster analysis, which showed that the 96-h samples were clustered far from their replicated samples (Supplementary Figure S2). We also found that the 48- and 72-h samples were clustered together within the same species, which might indicate that the changes in gene expression within the samples from 48 to 72 h were not obvious. This result suggested that the clustering order of all the samples was roughly consistent with the sampling time. These samples were also subjected to a principal component analysis (PCA) (Supplementary Figure S3). The first and second principal components (PC1 and PC2) showed that the 96-h samples were separated from the other samples. The 0- and 24-h Yanfeng47 samples were also separated by a long distance, which indicated that the susceptible varieties were more responsive to stress, and this this finding was consistent with the results shown in Figure S1.
Differentially expressed genes
Differentially expressed genes (DEGs) were identified by comparing the gene expression data obtained from Gangyuan8 and Yanfeng47 leaves at five time points. At all tested time points after inoculation, the number of upregulated genes exceeded the number of downregulated genes in both cultivars, and the number of upregulated genes in Yanfeng47 leaves was higher than that in Gangyuan8 leaves (Table 1). In both cultivars, the number of DEGs continued to increase from 0 to 96 h after inoculation. In Gangyuan8, the number of DEGs identified 96 h after inoculation was significantly higher than that found at other time points, and the number of DEGs showed an increasing trend from 0 to 96 h after inoculation. The results showed that the influence of the pathogen on the Gangyuan8 leaves decreased from 0 to 96 h after inoculation through a self-regulation mechanism. This study compared the DEGs of the two varieties at the same time points (GG24-YY24, GG48-YY48, GG72-YY72 and GG96-YY96; Figure 1) and compared the DEGs of each variety at the different time points (GG24 h, GG48 h, GG72 h, GG96 h, YY24 h, YY48 h, YY72 h and YY96 h; Figure 2). As shown in Figure 3, 1656 and 3688 DEGs were identified in the two varieties at 72 and 96 h after inoculation, respectively, and the number of these DEGs was higher than that at other time points. In Gangyuan8 leaves, 555 DEGs were shared among the DEGs at 24, 48, 72 and 96 h after inoculation, and these included 409 upregulated genes and 122 downregulated genes. A total of 2460 DEGs were found in the Yanfeng47 leaves, and these included 1135 upregulated genes and 1159 downregulated genes. As indicated in Figure 4, from 24 to 96 h after inoculation, the number of DEGs in Yanfeng47 leaves was higher than that in Gangyuan8 leaves.
Table 1
Statistics on the number of differentially expressed genes (DEGs)
DEG Set | Number of DEGs | Upregulated | Downregulated |
GG24 h | 1560 | 1186 | 374 |
GG48 h | 2780 | 2180 | 600 |
GG72 h | 3171 | 2406 | 765 |
GG96 h | 5953 | 3269 | 2684 |
YY24 h | 5783 | 3419 | 2364 |
YY48 h | 5804 | 3483 | 2321 |
YY72 h | 6216 | 3994 | 2222 |
YY96 h | 7112 | 4020 | 3092 |
Functional enrichment analysis
The enriched DEGs between the control and inoculation treatments were classified by GO annotation, and we found that the enriched DEGs were involved in multiple biological activities (Figure S4). All unigenes were annotated in the GO database and classified into three main categories. Among biological processes, “chloroplast envelope”, “chloroplast stroma” and “plasma membrane” were the three processes showing the highest degree of gene enrichment. The analysis of molecular functions revealed that “iron ion binding” was the only process with the highest degree of enrichment. In the cellular component category, “pentose-phosphate shunt”, “response to cadmium ion”, “response to cold”, and “response to salt stress” were the four processes with the highest enrichment degree.
In this study, some GO pathways were enriched in only a few DEGs (Figure S5). GO: 0016168, GO: 0018298, GO: 0016023 and GO: 0009523 were enriched only in the GG_0_VS_24 and GG_0_VS_48 h comparisons. This result indicated that the functional pathways involving these genes may be related to the early response of the varieties to pathogens. Our results suggest that at early time points in the rice-AG1 IA interaction, the genes related to chlorophyll binding, protein-chromophore linkage, cytoplasmic membrane-bounded vesicle and photosystem II are more active in Gangyuan8 leaves than in Yanfeng47 leaves.
Analysis of the gene expression profiles of rice cultivars
We also wanted to further analyze whether a certain gene is only expressed in the early stage of GG but not in YY; thus, these genes are likely related to the difference in the responses of the varieties to AG1 IA. Therefore, we extracted all DEGs and used them to generate expression profiles (Figure 3). In the diagram, the red colors represent high expression, and the blue colors represent low expression. The GG sample is shown on the left, and the YY sample is presented on the right. A total of 11 clusters were obtained, and the expression information of each cluster is shown in the line chart on the right of Figure 3. In this line chart, the red color represents the GG sample, and the black color represents the YY sample, as in cluster 9. The expression of these genes was high at 24 h and then decreased, and this profile is probably related to the early response of rice varieties to the pathogen. At the same time, the expression level of cluster 9 in the GG sample was higher than that in the YY sample, which may be due to the high expression of genes involved in resistance to infection in GG; thus, GG was identified as a disease-resistant variety. A similar trend was found for cluster 10. The analysis of clusters 9 and 10 obtained for GG revealed that the expression level of the genes related to cluster 9 increased continuously before 24 h and then decreased continuously after 24 h. In cluster 10, the increase and decreased occurred at 24 and 72 h, respectively. The GO enrichment results showed that the related genes in cluster 10 were mainly concentrated in the integrity of the membrane, metabolic process and cytoplasmic membrane-bound vesicle (Figure S6). The expression of cluster 11 in GG and YY remained highly after inoculation. In addition, the genes in cluster 11 were enriched in the response to salt stress, heme binding and extracellular matrix.
Metabolic pathway analysis after AG1 IA inoculation
We selected the DEGs of the two varieties at different time points for KEGG enrichment analysis (Figure S7). We found that some of the metabolic enrichments associated with resistance, namely, phenylpropanoid biosynthesis, alpha-linolenic acid metabolism, phagosomes, and the TCA cycle, were only enriched in Gangyuan8, and plant-pathogen interactions were only enriched in Yanfeng47. Based on the results of the KEGG enrichment analysis of the DEGs in Gangyuan8 at different stages, we found that immune-related metabolism was significantly enriched at the initial stage of infection, and phenylpropanoid biosynthesis was significantly enriched after 24 h and was only enriched in the resistant variety. Alpha-linolenic acid metabolism, which is highly related to jasmonic acid metabolism, was significantly enriched at 48 h. At the same time point, other resistance-related metabolisms, such as diterpenoid biosynthesis, were also significantly enriched. Based on the previous expression profile analysis, we preliminarily confirmed that the genes in clusters 9, 10 and 11 were involved in the disease resistance pathway of rice. To identify potential regulatory genes closely related to phenylalanine metabolic pathways, we identified DEGs by comparing the expression changes between the two varieties. At 96 h, the expression level of genes related to phenylalanine metabolism in the GG samples was significantly higher than that in the YY samples, and the difference between the two varieties can be seen in the heatmap (Figure 4). Jasmonic acid is an important signaling molecule in plants, and the biosynthesis of jasmonic acid in plants originates from linolenic acid. To explore the role of the jasmonic acid pathway in rice sheath blight resistance, we drew a heatmap of the related genes involved in alpha-linolenic acid metabolism and compared the differences between the two varieties. The results showed that the expression of genes related to alpha-linolenic acid metabolism was higher in the GG samples than in the YY samples. As shown in Figure 5, the expression level of related genes in the GG samples was always significantly higher than that in the YY samples. The expression levels of genes closely related to plant hormone signal transduction in the leaves of the two cultivars are shown in Figure 6. From 24 to 72 h, the expression levels of related genes in the GG samples were significantly higher than those in the YY samples.
Quantitative RT-PCR (qRT-PCR) validation of DEGs
To verify the reliability of the sequencing results, we selected 12 genes that were expressed in both the GG and YY samples for qPCR verification. The results showed that the expression trends of the 12 genes were similar to the sequencing results, which indicated that the sequencing results used in this study were relatively reliable (Figure S8).