circRNAs Are Involved in the Rice-Magnaporthe oryzae Interaction^{1  [OPEN]}

Abstract

Circular RNAs (circRNAs) play roles in various biological processes, but their functions in the rice (Oryza sativa) response to Magnaporthe oryzae remain elusive. Here, we demonstrate that circRNAs are involved in the rice-M. oryzae interaction using comparative circRNA-sequencing and transgenic approaches. We identified 2932 high-confidence circRNAs from young leaves of the blast-resistant accession International Rice Blast Line Pyricularia-Kanto51-m-Tsuyuake (IR25) and the blast-susceptible accession Lijiangxin Tuan Heigu (LTH) under M. oryzae-infected or uninfected conditions; 636 were detected specifically upon M. oryzae infection. The circRNAs in IR25 were significantly more diverse than those in LTH, especially under M. oryzae infection. Particularly, the number of circRNAs generated per parent gene was much higher in IR25 than in LTH and increased in IR25 but decreased in LTH upon M. oryzae infection. The higher diversity of circRNAs in IR25 was further associated with more frequent 3′ and 5′ alternative back-splicing and usage of complex splice sites. Moreover, a subset of circRNAs was differentially responsive to M. oryzae in IR25 and LTH. We further confirmed that circR5g05160 promotes rice immunity against M. oryzae. Therefore, our data indicate that circRNA diversity is associated with different responses to M. oryzae infection in rice and provide a starting point to investigate a new layer of regulation in the rice-M. oryzae interaction.

Rice (Oryza sativa) is a staple crop that feeds more than half of the global population (Liu et al., 2014). Rice blast, caused by the fungal pathogen Magnaporthe oryzae, is one of the most devastating diseases threatening rice production worldwide (Liu and Wang, 2016). Rice has multiple layers of defense against the invasion of blast fungus (Liu et al., 2013). First, several rice pattern-recognition receptors, including Chitin Elicitor Binding Protein (CEBiP), Chitin Elicitor Receptor Kinase1 (CERK1), Lysin Motif-Containing Protein4 (LYP4), and LYP6, recognize chitin fragments of the blast fungus and trigger defense responses (Shimizu et al., 2010; Liu et al., 2012) called pathogen-associated molecular pattern-triggered immunity (PTI). PTI can be suppressed by effector proteins secreted by M. oryzae (Khang et al., 2010; Mentlak et al., 2012). To counteract the pathogen, rice deploys resistance (R) proteins to recognize avirulence effectors from M. oryzae, leading to a second layer of defense called effector-triggered immunity (ETI). More than 100 R genes have been identified, and about 30 of them have been functionally characterized to act as on-off switches in regulating rice blast resistance (Wang et al., 2017a). Some of the R genes, such as Pigm and Pi2, have been widely exploited in blast disease-resistance programs (Shi et al., 2015; Deng et al., 2017). Both PTI and ETI can be modulated by microRNAs (miRNAs; Padmanabhan et al., 2009).

miRNAs are a class of noncoding RNAs that act in various biological processes and stress-induced responses. Increasing evidence supports the role of miRNAs in fine-tuning rice immunity against M. oryzae. Genome-wide small RNA analyses have revealed a number of candidate miRNAs responsive to M. oryzae infection or elicitors (Campo et al., 2013; Li et al., 2014; Li et al., 2016; Wang et al., 2018). Transgenic approaches further confirm miRNAs, such as miR7695, miR398b, miR160a, and miR166k-166h (Campo et al., 2013; Li et al., 2014, 2019; Salvador-Guirao et al., 2018), positively regulate rice blast resistance; on the contrary, miR169a, miR164a, miR319b, miR396, and miR167d (Li et al., 2017; Wang et al., 2018; Zhang et al., 2018; Chandran et al., 2019; Zhao et al., 2019) act as negative regulators of rice immunity against the blast pathogen. miRNAs exert their biological functions via targeting mRNAs with sequence complementarity (Seitz, 2009); meanwhile, binding of miRNAs to their targets could be interfered with by endogenous noncoding RNAs, such as target mimics (Wu et al., 2013) and circular RNAs (circRNAs; Hansen et al., 2013). Target mimics have been reported to regulate rice immunity against the blast pathogen (Li et al., 2017; Chandran et al., 2019; Li et al., 2019), but it is unknown whether circRNAs are involved in the rice-M. oryzae interaction.

circRNAs are produced from precursor messenger RNAs through back-splicing in which an upstream 3′ splicing acceptor site is joined to a downstream 5′ splicing donor site (Ashwal-Fluss et al., 2014). According to their genomic origins, circRNAs are mainly classified as exonic, intronic, exon-intronic, and intergenic (Bolha et al., 2017). circRNAs regulate gene expression in animals through different mechanisms, such as miRNA sponges, specific RNA-RNA interactions, and affecting alternative splicing (Li et al., 2018). Since the identification of their roles in regulating gene expression in animals (Memczak et al., 2013), circRNAs have also been identified in different plant species, including Arabidopsis (Arabidopsis thaliana), rice, wheat (Triticum aestivum), barley (Hordeum vulgare), maize (Zea mays), tomato (Solanum lycopersicum), potato (Solanum tuberosum), soybean (Glycine max), cotton (Gossypium hirsutum), and kiwifruit (Actinidia deliciosa; Wang et al., 2014, 2017b, 2017c; Lu et al., 2015; Ye et al., 2015, 2017; Darbani et al., 2016; Zuo et al., 2016; Chen et al., 2017, 2018; Zhao et al., 2017).

While reports identifying plant circRNAs are rapidly increasing, investigation of their roles is far more challenging. Plant circRNAs have been found to function during fiber development, flowering, and fruit coloration based on spatial- and/or tissue-specific expression patterns (Tan et al., 2017; Wang et al., 2017c; Tong et al., 2018). The expression of plant circRNAs also responds to abiotic and biotic stimuli, including pathogen invasion (Zuo et al., 2016; Wang et al., 2017c; Zhou et al., 2017; Ghorbani et al., 2018; Xiang et al., 2018; Gao et al., 2019). For example, the accumulation of a set of circRNAs is altered in kiwifruit upon canker pathogen infection (Wang et al., 2017c). However, it is unknown whether circRNA production and expression changes are functionally linked to plant disease resistance.

To explore whether circRNAs are involved in rice blast resistance, we performed genome-wide identification of circRNAs in leaf transcriptomes of the resistant accession International Rice Blast Line Pyricularia Kanto51-m-Tsuyuake (designated as IR25) and the blast-susceptible accession Lijiangxin Tuan Heigu (LTH) before and after M. oryzae infection via the ribosomal RNA-depleted RNA-sequencing technique (Lu et al., 2015). circRNAs were identified with the CIRI2 pipeline, and their numbers and abundance were compared between IR25 and LTH, and/or between the samples during M. oryzae infection and without infection. circRNAs in IR25 were more diverse than those in LTH regardless of M. oryzae infection, which was associated with the increased complexity of alternative splicing sites. We also obtained transgenic lines overexpressing circR5g05160 and confirmed that up-regulation of this circRNA led to enhanced blast disease resistance. Overall, our results indicate that increased alternative splicing complexity contributes to the diversity of circRNAs, and circRNAs are involved in the rice-M. oryzae interaction.

RESULTS

Identification of High-Confidence circRNAs from the Leaf Transcriptomes of Blast-Resistant and Blast-Susceptible Rice Accessions

To explore the different responses to M. oryzae between blast-resistant and blast-susceptible rice accessions, we performed RNA-sequencing (RNA-seq) on the leaves from 3-week-old seedlings of IR25 and LTH before (IR25-0h, LTH-0h) and after M. oryzae infection (IR25-12h, LTH-12h, IR25-24h, LTH-24h). In total, we obtained 18 transcriptome data sets with each containing an average of about 200 million paired-end reads for analysis (Supplemental Table S1). Then, we performed principle component analysis based on the expression level of each gene and found a clear separation of both the LTH and IR25 rice samples between before (0 h) and after M. oryzae infection (12 h and 24 h) in two experiment benches (Supplemental Fig. S1A). Samples from bench 1 and bench 2 were also obviously separated (Supplemental Fig. S1A). However, there was no clear separation between LTH and IR25 samples at each time point from either experimental bench (Supplemental Fig. S1A). These results indicated the overall gene expression profiles were similar between IR25 and LTH, regardless of M. oryzae infection.

Next, we performed edgeR analysis to identify differentially expressed genes (DEGs; ≥2-fold increase or decrease, false discovery rate [FDR] < 0.05) between the two accessions at each time point. From one bench, a total of 914, 1456, and 1758 DEGs were found in IR25 compared to the LTH at 0, 12, and 24 h, respectively (Supplemental Fig. S1B, left). From the other bench, the number of DEGs were 4352, 2687, and 2855 at 0, 12, and 24 h, respectively (Supplemental Fig. S1B, right). For expression levels of all unique DEGs, unsupervised hierarchical clustering analysis revealed a clear difference between LTH and IR25 rice samples in response to M. oryzae infection (Supplemental Fig. S1C). As to the up- and down-regulated DEGs between IR25 and LTH, gene ontology (GO) analysis uncovered a consistently significant enrichment in “response to stress” and “response to biotic stimulus” (Supplemental Fig. S2). These data indicate that the transcriptomes of IR25 and LTH are differently modulated in response to infection with M. oryzae and thus are suitable to screen circRNAs differentially responsive to M. oryzae.

To identify circRNAs from leaf transcriptome data of both IR25 and LTH, we removed low-quality reads, and the rest of the reads were run with a custom pipeline based on CIRI2, followed by a validation mapping step (Supplemental Fig. S1D, left). In the validation step, we aligned all reads onto an artificial circRNA sequence database built for all circRNAs retrieved from the CIRI2 pipeline, and only circRNAs with at least two mapped reads covering the back-splicing junction were kept as candidate circRNAs (Supplemental Fig. S1D, right). A total of 374 to 1527 circRNAs were identified in each of the IR25 and LTH transcriptome samples (Fig. 1A). Overall, 2932 circRNAs were identified in all samples, of which 72.74%, 22.21%, and 5.05% were exonic, intergenic, and intronic circRNAs, respectively (Supplemental Table S2). These circRNAs were distributed across the whole genome and were more commonly found at both ends of chromosomes (Fig. 1B). Such distribution is similar to a recent report in maize (Chen et al., 2018). A total of 1646 circRNAs were previously reported and deposited in PlantcircBase for rice (release 3; Chu et al., 2017), while 1286 circRNAs were only detected in this study (Supplemental Table S2). Moreover, 636 circRNAs were specifically triggered by M. oryzae infection (Supplemental Table S3). Out of them, 411 and 144 circRNAs were detected only in M. oryzae-infected IR25 and LTH, respectively; 81 circRNAs were common in both M. oryzae-infected IR25 and LTH (Supplemental Table S3). Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis indicated that the parent genes producing M. oryzae-triggered circRNAs were involved in multiple stress-responsive pathways, such as biosynthesis of secondary metabolites, protein export, and terpenoid backbone biosynthesis (Supplemental Table S4). These results indicate that circRNAs widely exist in leaf transcriptomes of both IR25 and LTH and M. oryzae infection induces more production of circRNAs in IR25 than in LTH.

To validate these circRNAs, we randomly selected 12 circRNAs and validated them by PCR using pairs of divergent and convergent primers (Supplemental Table S5). All 12 pairs of convergent primers successfully amplified the expected length of fragments from both cDNAs of total RNAs and genomic DNA (Fig. 1C, Supplemental Fig. S3A). By contrast, the divergent primers could only amplify fragments from cDNA of total RNAs, but not from genomic DNA (Fig. 1C; Supplemental Fig. S3A). The amplification products from the divergent primers were confirmed to span the back-splicing junction of circRNAs via Sanger sequencing and sequence mapping (Fig. 1D; Supplemental Fig. S3B). Moreover, all 12 pairs of divergent primers yielded amplification products from cDNAs of total RNAs with or without treatment of RNase R that digests linear RNAs (Fig. 1E). These data indicate the validation of all 12 randomly selected circRNAs. Thus, the circRNAs identified by our analysis pipeline are highly reliable.

The Diversity of circRNAs Is Constantly Higher in IR25 Than in LTH Regardless of M. oryzae Infection

We noticed more circRNAs were repeatedly identified in IR25 than in LTH (Fig. 1A). This prompted us to determine the difference between blast-resistant and susceptible rice accessions in producing circRNAs. A total of 2633 and 1787 circRNAs were identified from all samples of IR25 and LTH, respectively, of which 1488 circRNAs overlapped (Fig. 2A, top), indicating most of the LTH-produced circRNAs are included in those of IR25. Focusing on the circRNAs detected in at least two samples, 1204 and 688 circRNAs remained in IR25 and LTH, respectively, of which 588 overlapped (Fig. 2A, bottom), indicating more diverse circRNAs in IR25 than in LTH. To remove the potential effect of sequencing bias on the number of identified circRNAs, we obtained the normalized number of circRNAs in each million of high-quality sequencing tags (Fig. 2B), which confirmed IR25 generated more circRNAs than did LTH.

To explore whether the increased circRNA diversity in IR25 depended on the M. oryzae infection, we plotted the normalized number of circRNAs to each postinfection time point. The results demonstrated that M. oryzae infection led to a marginal increase in the circRNA diversity; however, the circRNA diversity in IR25 was constantly higher than that in LTH (Fig. 2C). We then analyzed circRNAs located in the gene regions and found that the diversity increase became more pronounced in IR25 but not in LTH (Fig. 2D), indicating the biogenesis of diverse circRNAs from genic loci might be specifically triggered in the blast-resistant accession by the M. oryzae infection.

To further dissect the potential origins of the increased circRNA diversity in IR25, we analyzed the number of circRNAs generated by each gene. The average number of circRNAs per gene was constantly higher in IR25 than in LTH (Fig. 2E), implying genes in IR25 are more capable of generating circRNAs than those in LTH. Moreover, the numbers of circRNA-generating genes were 275 in common, 360 specifically in IR25, and 57 specifically in LTH (Fig. 2F; Supplemental Table S6), indicating there were more circRNA-generating genes in IR25 than in LTH, and most (∼82.9%) of the circRNA-generating genes in LTH were included in IR25. Then, we analyzed the average number of circRNAs generated by each gene of the 275 circRNA-producing genes shared by IR25 and LTH. The average number of circRNAs generated per gene was significantly higher in IR25 than in LTH (Fig. 2G). Upon infection with M. oryzae, the number of circRNAs generated per gene was obviously reduced in LTH; by contrast, it was increased in IR25, and the increase became significant at 24 h post inoculation (hpi) with M. oryzae (Fig. 2G). Taken together, our results indicate that IR25 is capable of producing more circRNAs than LTH, particularly upon M. oryzae infection.

It has been reported that circRNAs can regulate the expression of their parent genes, therefore exerting a biological function associated with their parent genes (Li et al., 2018). To uncover whether the circRNA-generating genes possess any blast resistance functions, we subjected the circRNA-generating genes to functional enrichment analysis according to the circRNAs specifically detected in IR25 and LTH and those shared by IR25 and LTH. We found regulation of autophagy and peroxisome were among the top functional terms enriched by circRNA-generating genes and were particularly more enriched by IR25-specific genes (Fig. 2H; Supplemental Fig. S4). As autophagy and peroxisomes have been documented to function in plant immunity (Mammarella et al., 2015; Hofius et al., 2017), our findings imply a link between circRNAs and blast resistance.

Alternative Splicing Contributes to the Increased circRNA Diversity in IR25

Since alternative back splicing is required for the biogenesis and diversification of circRNAs in rice (Lu et al., 2015; Ye et al., 2017), we hypothesized more diverse circRNAs in IR25 may be attributed to more complicated alternative splicing in IR25 than in LTH. Therefore, we annotated the downstream 5′ back-splice sites and upstream 3′ back-splice sites from all the identified circRNAs. As expected, the numbers of both alternative 5′ back-splice sites and 3′ back-splice sites were significantly higher in IR25 than in LTH regardless of M. oryzae infection (Fig. 3, A and B). We further viewed the alternative back-splicing events by calculating their percent circularized-site usage (PCU); each isoform of an alternative back-splicing event had its own PCU (Supplemental Table S7). The heatmap plot of PCUs showed alternative back-splicing events were more prevalent in IR25 samples than in LTH samples (Fig. 3C). These results indicated alternative back-splicing patterns were more diverse in IR25 than in LTH, which could contribute to the increased circRNA diversity in IR25.

Internal alternative splicing inside circRNAs is also essential for the biogenesis and diversity of circRNAs in animals (Gao et al., 2016; Zhang et al., 2016). A number of rice circRNAs consist of multiple exons that could be derived from internal splicing events of circRNAs (Ye et al., 2017). Therefore, we hypothesized that more diverse circRNAs in IR25 may also be attributed to more internal alternative splicing in circRNAs of IR25 than those of LTH. To this end, we adopted CircRNA Identifier-Alternative Splicing (CIRI-AS; Gao et al., 2016) to identify and analyze the internal alternative splicing in circRNA. Indeed, internal alternative splicing events inside circRNAs were significantly more prevalent in IR25 than in LTH, although all alternative splicing events were not obviously different between IR25 and LTH (Fig. 3D). In addition, all four basic types of alternative splicing were higher in IR25 than in LTH (Fig. 3E).

Taken together, our results demonstrated both alternative back-splicing and internal alternative splicing inside circRNAs may contribute to the higher diversity of circRNAs in IR25 than in LTH.

circRNAs Are More Likely Derived from Complex Splice Sites

Recent studies in mammals have shown the presence of very complex splice sites that generate splice junctions with more than two donors or acceptors (Gao et al., 2016; Zhang et al., 2016; Feng et al., 2019; Zheng et al., 2019). We speculated complex splice sites should have a higher probability of successful back-splicing. To test this hypothesis, we analyzed the number of splice junctions formed by each detected splice site in all samples. The number of splice junctions displayed similar profiles in both IR25 and LTH (Supplemental Fig. S5). Most splice sites underwent only one splicing event, which yielded one splice junction. Nevertheless, a significant fraction of splice sites underwent two or three splicing events, and thousands of splice sites underwent four or more splicing events (Supplemental Fig. S5).

We then analyzed whether the linear and circRNA splice sites were associated with their complexity. The total number of circRNA splice sites was less than that of the linear splice sites (Fig. 4A, left). However, the percentage of complex splice sites (≥2 junctions) in circRNAs was more than that in linear transcripts; by contrast, the percentage of simple splice sites (junction = 2) in circRNAs was significantly lower than that in linear transcripts (Fig. 4A, right). We further explored whether the association between circRNAs, and their splice site complexity was different in IR25 and LTH. Both the number (Fig. 4B, left) and percentage (Fig. 4B, right) of complex splice-site-associated circRNAs were higher in IR25 than in LTH, although the differences were not significant. Using LOC_Os11g36030 and LOC_Os2g148000 as examples, the reads distribution demonstrated diverse circRNAs were generated from more complicated splice sites in IR25 than in LTH (Fig. 4C). These results indicate IR25 has a high capacity for supporting complex splicing, which may facilitate back-splicing and circRNA production.

The Abundance of Some circRNAs Is Significantly Altered upon M. oryzae Infection

To explore whether circRNA abundance responded to M. oryzae infection, we violin-plotted circRNA expression level (total circRNA junction reads per mapped million reads). Total circRNA abundance was marginally different among samples (Supplemental Fig. S6). We then screened circRNAs differentially expressed in IR25 and LTH in response to M. oryzae infection and identified 20 and 31 up-regulated circRNAs and eight and 16 down-regulated circRNAs in IR25 at 12 and 24 hpi, respectively. By contrast, we identified six and five up-regulated circRNAs and 21 and 13 down-regulated circRNAs in LTH at 12 and 24 hpi, respectively (Fig. 5A; Supplemental Table S8). Apparently, the number of up-regulated circRNAs was more than that of down-regulated circRNAs in IR25 upon M. oryzae infection; on the contrary, more circRNAs tended to be down-regulated in LTH (Fig. 5A). There were 18 circRNAs constantly up-regulated and six circRNAs constantly down-regulated at both 12 and 24 hpi in IR25, and two circRNAs constantly up-regulated and eight circRNAs constantly down-regulated in LTH (Fig. 5B). These circRNAs could be the ones involved in rice response to M. oryzae. As an example, the expression of a circRNA from LOC_Os12g19381 was constantly down-regulated upon M. oryzae infection in both IR25 and LTH (Fig. 5C). LOC_Os12g19381 was annotated with the putative function of ribulose bisphosphate carboxylase, implying M. oryzae infection might affect rice photosynthesis.

Next, we screened circRNAs differentially expressed between IR25 and LTH. We identified 38, 40, and 51 circRNAs with higher expression and 15, 12, and 14 circRNAs with lower expression in IR25 than in LTH at 0, 12, and 24 hpi, respectively (Fig. 5D). These circRNAs may be associated with regulation of blast disease resistance. Moreover, a significant number of circRNAs overlapped among different time points (Fig. 5E; P < 0.05). As an example, a circRNA from LOC_Os11g11890 was shown to have constantly higher expression in IR25 than in LTH (Fig. 5F). LOC_Os11g11890 was annotated as an ortholog of the disease-resistance gene RPG1 in barley (Brueggeman et al., 2002), suggesting IR25 might deploy circRNAs for blast resistance.

To validate these differentially expressed circRNAs, we randomly selected 12 circRNAs for reverse transcription quantitative PCR (RT-qPCR). All 12 circRNAs were differentially expressed after M. oryzae infection in both IR25 and LTH, which was consistent with the circRNA-seq results (Fig. 5G; Supplemental Fig. S7). These results demonstrated the abundance of most rice circRNAs was higher in IR25 than in LTH, and infection of M. oryzae led to more up-regulated circRNAs in IR25 but more down-regulated circRNAs in LTH. Therefore, the differentially expressed circRNAs should be a priority for functional characterization in the future.

Overexpression of circR5g05160 Leads to Enhanced Blast Resistance

To confirm the involvement of circRNAs in regulation of the rice-M. oryzae interaction, we selected a circRNA Chr5:2512798|2514806 derived from LOC_Os05g05160 (circR5g05160) for functional analysis. We selected circR5g05160 because it was preferentially expressed in IR25 (Supplemental Table S6), and its parent gene encodes a putative plant immunity regulator MPK14 (Tena et al., 2011); moreover, the predicted circR5g05160 sequence contained potential targeting sites of miRNAs (Supplemental Fig. S8C). We confirmed by qPCR that circR5g05160 was differentially responsive to M. oryzae with much higher induction in IR25 than in LTH, especially at early infection stages (Fig. 6A). We then validated this circRNA via PCR with divergent and convergent primers (Supplemental Fig. S8A; Supplemental Table S5) and confirmed the junction site by Sanger sequencing (Supplemental Fig. S8B). We also amplified the full length of circR5g05160 using the divergent primers and found circR5g05160 is a 607 nucleotide exonic circRNA (Supplemental Fig. S8C). There was an alternative splicing event within circR5g05160 resulting in a 10-nucleotide shift upstream of the second exon (Supplemental Fig. S8C). Sequence analysis confirmed one potential targeting site for osa-miR168-5p and one for osa-miR2103 (Supplemental Fig. S8, C and D).

To overexpress circR5g05160 in rice, we first cloned the DNA fragment containing predicted circR5g05160 and its endogenous flanking introns into the pCAMBIA1300 vector under the control of the Cauliflower mosaic virus 35S promoter (Supplemental Fig. S8E). We then tested whether this construct (p35S-circR5g05160) could efficiently produce circR5g05160 in the Nicotiana benthamiana transient expression system by qPCR with divergent primers and Sanger sequencing. circR5g05160 was highly expressed in N. benthamiana, and the junction site and full-length sequence were exactly the same as in rice (Supplemental Fig. S8, C and F), indicating the construct p35S-circR5g05160 can express the full length of circR5g05160. Next, we introduced p35S-circR5g05160 into the rice accession Nipponbare (NPB) via Agrobacterium tumefaciens-mediated transformation and obtained multiple independent transgenic rice lines that overexpressed circR5g05160 (Fig. 6B). In these lines, linear transcript accumulation of LOC_Os05g05160 was not greatly changed (Fig. 6B). Subsequently, two transgenic lines were used for blast disease assay. Both lines formed much smaller disease lesions and supported significantly less fungal growth of M. oryzae than did NPB (Fig. 6, C–E), indicating enhanced resistance to the blast disease.

To understand how circR5g05160 is involved in blast resistance, we examined the expression of some defense-related genes in two transgenic lines overexpressing circR5g05160. PTI-related genes, such as NAC4, PBZ1, and PR10b in rice (Li et al., 2014), were expressed at higher levels in overexpressing lines OX20 and OX24 than in wild-type NPB; upon M. oryzae infection, their expression levels were induced much higher in OX20 and OX24 than in NPB, especially at 12 hpi (Fig. 6, F–H). ETI-related genes, such as HSP90, SGT1, and RAR1 (Thao et al., 2007; Shirasu, 2009), also displayed differential expression in OX20 and OX24 compared to NPB. The basal expression levels of HSP90, SGT1, and RAR1 were 10- to 80-fold higher in OX20 and OX24 than in NPB, although their expression levels were mostly similar between overexpressing lines and NPB after M. oryzae infection (Fig. 6, I–K). These results indicate that circR5g05160 may be involved in both PTI and ETI to regulate rice immunity against M. oryzae.

DISCUSSION

circRNAs are widespread and diverse in both animals and plants and have potential regulatory functions (Li et al., 2018). In plants, circRNAs are closely associated with development and stress-induced responses (Lu et al., 2015; Ye et al., 2015, 2017). Here, we demonstrated that circRNAs are involved in the rice-M. oryzae interaction. First, out of 2932 high-confidence circRNAs identified in this study, 636 circRNAs were specifically generated upon M. oryzae infection (Supplemental Table S3). Second, more circRNAs were produced in the blast-resistant accession IR25 than in the blast-susceptible accession LTH, and circRNA diversity was significantly increased by M. oryzae infection in IR25 but not in LTH (Fig. 2; Supplemental Table S3). The higher capability of circRNA biogenesis in IR25 was further associated with more alternative back-splicing, internal alternative splicing inside circRNAs, and usage of complex splice sites (Figs. 3 and 4). Third, a subset of circRNAs was differentially responsive to M. oryzae in IR25 and LTH (Fig. 5), and functional analysis showed that circR5g05160 contributed to rice immunity against M. oryzae (Fig. 6).

As a posttranscriptional process in eukaryotic organisms, alternative splicing leads to production of multiple, distinct functional transcripts and diverse proteins from a single gene (Black, 2000; Nilsen and Graveley, 2010). Both alternative back-splicing and alternative splicing inside circRNA are required for the biogenesis and diversity of circRNAs in animals (Gao et al., 2016; Zhang et al., 2016; Feng et al., 2019; Zheng et al., 2019). In fact, alternative circularization of circRNA is reported as a common feature in rice, which results in a set of circRNA isoforms from the same locus (Lu et al., 2015; Ye et al., 2017). In this study, we found more circRNA production was associated with more 5′ and 3′ alternative back-splice sites and more internal alternative splice sites inside circRNAs in IR25 than in LTH (Fig. 3), implying that internal alternative splicing, as well as alternative back-splicing, may contribute to different responses to M. oryzae in rice.

Complex splice sites can generate splice junctions with more than two donors or acceptors (Gao et al., 2016; Zhang et al., 2016; Feng et al., 2019; Zheng et al., 2019). As alternative back-splicing events compete with the canonical splice site and alternative splice site, complex splice sites may have more chances to produce circRNAs. In this study, the percentage of complex splice sites was significantly higher in circRNAs than in linear transcripts (Fig. 4). Moreover, both the number and percentage of complex splice sites were obviously higher in IR25 than in LTH, which is positively associated with the capability of producing circRNAs in rice accessions (Fig. 4). Therefore, complex splice sites may be involved in circRNA biogenesis in response to M. oryzae in rice.

We observed that M. oryzae infection could specifically trigger the production of some circRNAs in rice, of which the parent genes were enriched in defense-related pathways, such as “Biosynthesis of secondary metabolites” and “Terpenoid backbone biosynthesis” (Piasecka et al., 2015; Supplemental Table S4). Particularly, a number of circRNAs (411) were specifically induced in the resistant accession IR25 upon M. oryzae infection. Their parent genes were significantly enriched in the KEGG pathway “Splicesome” (Supplemental Table S4), implying splicing-related genes were strongly modulated, likely contributing to increased circRNA diversity in IR25. Interestingly, for circRNA parent genes shared by IR25 and LTH, those in IR25 showed a higher ability to generate circRNAs than those in LTH, especially under M. oryzae infection (Fig. 2G). Moreover, IR25-specific parent genes were significantly enriched in peroxisome- and autophagy-related pathways (Fig. 2H), which are involved in plant immunity and disease resistance (Hofius et al., 2009; Daudi et al., 2012). For instance, peroxidases PRX33 and PRX34 are required for apoplastic reactive oxygen species production and basal resistance to pathogens in Arabidopsis (Daudi et al., 2012). Autophagic components such as ATG6 function in hypersensitive cell death and plant immunity (Hofius et al., 2009; Yue et al., 2015). Taken together, our data support that production of rice circRNAs is responsive to M. oryzae infection, and some circRNAs may be involved in rice immunity against M. oryzae.

Numerous circRNAs have been identified by genome-wide analysis in plants, but relatively few have been functionally characterized. For instance, an exonic circRNA from Arabidopsis SEPALLATA3 causes floral homeotic phenotypes via increasing the expression of the cognate exon-skipped alternative splicing isoform (Conn et al., 2017). A grape (Vitis vinifera) circRNA Vv-circATS1 enhances cold tolerance when ectopically expressed in Arabidopsis (Gao et al., 2019). The rice circRNA Os08circ16564 was successfully overexpressed in rice accession NPB, but the resulting transgenic plants displayed no obvious phenotypes, and the biological function of Os08circ16564 is still unknown (Lu et al., 2015). In this study, we demonstrated the function of circR5g05160 in rice resistance to blast disease (Fig. 6), providing the first line of evidence that circRNAs function in the rice-M. oryzae interaction. In addition, we found potential miRNA-targeting sites in circR5g05160 (Supplemental Fig. S8, C and D); whether this circRNA modulates rice immunity through sponging miRNAs could be a future research focus. Alternatively, circR5g05160 may interfere with the function of its parent gene OsMPK14, which belongs to a gene family critical for plant immunity (Tena et al., 2011).

CONCLUSION

In summary, we have identified a much larger number/diversity of expressed circRNAs in leaves of blast-resistant IR25 than in that of blast-susceptible LTH. The difference may be associated with different responses to M. oryzae and attributed to differences in efficiency of alternative back-splicing and usage of complex splice sites. We also demonstrated a strong link between circRNA production and rice blast resistance by functional analysis. These results provide new insights into circRNA biogenesis in rice and uncover regulatory factors in rice immunity, which will help us to understand the complicated regulatory network in the rice-M. oryzae interaction.

MATERIALS AND METHODS

Rice Accessions and Fungi

Rice (Oryza sativa) plants used in this study were the blast-susceptible japonica accession LTH and a corresponding blast-resistant accession IR25. IR25 is a monogenic line containing the blast resistance gene Pikm that was introduced into LTH by backcrossing with Tsuyuake (Tsunematsu et al., 2000). All rice plants were grown at 26°C ± 2°C and 70% relative humidity under 12 h light:12 h darkness. Magnaporthe oryzae strains (NC-24, NC-34, B9, Zhong8-10-14, 97-95-1, F1, NC-14, B13, and E37) were cultured in complete media at 28°C under 12 h light:12 h darkness for 2 weeks. Spores were collected and adjusted to a concentration of 2 × 10⁵ spore mL⁻¹. Equal volume of spore suspension for each strain were mixed and sprayed onto 3-week-old seedlings of LTH and IR25. Leaves were collected at 0, 12, and 24 hpi, immediately frozen in liquid nitrogen, and stored at −80°C until use.

circRNA-seq Library Construction and Sequencing

Leaf tissues of both rice accessions from the following three stages were collected for library construction: before inoculation (0 h) and after inoculation (12 and 24 h). Total RNA of each sample was isolated using TRIzol (Ambion) according to the manufacturer’s instructions. To remove the noise from differences in individual plants and experimental variation, we generated a total of three sets of RNA-seq data from two independent benches. In one experimental bench, we mixed total RNAs from different plants of two biological repeats for library preparation. In the other bench, we used total RNAs from each biological repeat for library construction (Supplemental Table S1).

Total RNA was treated with RQ1 DNase (Promega) to remove DNA. The quality and quantity of the purified RNA were determined by measuring the absorbance at 260 and 280 nm using SmartSpec Plus (Bio-Rad). RNA integrity was further verified by 0.015 g mL⁻¹ agarose gel electrophoresis.

For each sample, 25 µg of total RNA was used for circRNA-seq library preparation. Total RNA was depleted of ribosomal RNAs using the RiboMinus kit (Illumina). The remaining RNAs were fragmented at 95°C followed by end repair and 5′ adaptor ligation. Then reverse transcription (RT) was performed with the RT primer harboring a 3′ adaptor sequence and randomized hexamer. The cDNAs were purified and subjected to PCR amplification. PCR products corresponding to 300 to 500 bp were purified and quantified and then stored at −80°C until they were used for sequencing.

For high-throughput sequencing, the libraries were prepared following the manufacturer’s instructions (Illumina), and the Illumina NextSeq 500 system and the HiSeq X Ten system were used for 150-nucleotide paired-end sequencing by ABlife.

Identification and Quantification of circRNAs

For genome-wide identification of circRNAs, reads containing adapter or poly-N and low-quality reads were filtered from the raw sequencing reads by in-house Perl scripts. The resultant clean reads were mapped to the rice reference genome (IRGSP v5.0), generating a sequence alignment map file. Sequence alignment map files were then used to identify circRNAs with CIRI2 as described previously (Gao et al., 2015). Based on their genomic origins, circRNAs were classified into three major types: exonic, intronic, and intergenic circRNAs. For quantification of circRNAs (Song et al., 2016), the numbers of back-spliced reads of each circRNA was normalized to the total sequencing reads in a corresponding sample data set, defined as reads per million mapped reads.

DEGs and circRNA (DEC) Analysis

The R Bioconductor package edgeR (Robinson et al., 2010) was utilized to screen out DEGs and DECs. FDR < 0.05 and fold change ≥ 2 were set as the cut-off criteria for identifying DEGs, and P < 0.05, fold change ≥ 1.5 for identifying DECs.

To sort out functional categories of DEGs and genes hosting DECs, GO terms and KEGG pathways were identified using the KOBAS 2.0 server (Xie et al., 2011). A hypergeometric test and the Benjamini-Hochberg FDR control procedure were used to define the enrichment of each term.

Analysis of Alternative Splicing in circRNAs

Analysis of alternative splicing was performed using CIRI-AS as previously described (Gao et al., 2016). CIRI-AS is a novel algorithm for detecting internal components of circRNAs based on split alignment of back-splice junction read pairs and distribution of sequencing depth (Gao et al., 2016).

Genomic DNA, Total RNA Extraction, and RNase R Treatment

Genomic DNA was extracted from fresh tissue using the cetyltrimethylammonium bromide method (Murray and Thompson, 1980). Genomic DNA was used as a negative control for PCR with divergent primers. Total RNA was isolated from all samples using TRIzol reagent (Ambion) according to the manufacturer’s protocol. Total RNA samples were treated with DNase I (NEB) and purified by RNA Clean & Concentrator-25 spin columns (ZYMO Research) to remove DNA contamination and salts. RNA was evaluated using 0.02 g mL⁻¹ Tris-acetate-EDTA-agarose gel electrophoresis. For RNase-R treatment, the purified DNase I-treated total RNA was incubated for 15 min at 37°C with 3 units/µg of RNase R (Epicentre) and subsequently purified by phenol-chloroform extraction and reprecipitated in three volumes of ethanol.

PCR Amplification and Sanger Sequencing

The divergent and convergent primers were designed for circRNA validation (Supplemental Table S5). Convergent primers were used as positive controls for linear transcripts, and divergent primers were used to detect the candidate circular template. For each PCR amplification, 20 ng of cDNA or genomic DNA was used with 2× Phanta Master Mix (Vazyme Biotech). PCR products of the expected length were dissected from a gel and purified using the GeneJET Gel Extraction kit (Invitrogen). Purified PCR products were Sanger sequenced at Sangon Biotech or TsingKe Biological Technology Company.

RT-qPCR

For RT-qPCR, first-strand cDNA was retro-transcribed from total RNA with or without RNase R treatment with random six mers and SuperScript III reverse transcriptase (Invitrogen) and subjected to PCR reactions with the SYBR Green PCR MasterMix (Takara) on a Bio-Rad CFX Connect Real-Time PCR Detection System. The expression level of circRNA was normalized to endogenous linear rice ubiquitin (OsUbi, AK059011) transcripts. Each set of experiments was repeated three times. The primers used for RT-qPCR are listed in Supplemental Table S5.

Plasmid Construction and Overexpression of circR5g05160

The genomic region of circR5g05160 containing the endogenous flanking introns (Supplemental Fig. S8E) was amplified from LTH genomic DNA with primers Chr5:2512798|2514806-KpnI-F/SalI-R (Supplemental Table S5) and cloned into a pCAMBIA1300 vector under the control of Cauliflower mosaic virus 35S promoter leading to a recombinant plasmid p35S-circR5g05160. For transient expression, Agrobacterium tumefaciens GV3101 cells containing p35S-circR5g05160 were adjusted to optical density at 600 nm = 0.1 to 0.3 with 10 mm MgCl₂ and infiltrated into the fully expanded leaves of Nicotiana benthamiana. At 2 d postinfiltration, leaves were collected for total RNA extraction and qPCR. Leaves infiltrated with10 mm MgCl₂ were used as the control (CK). The production and expression levels of circR5g05160 were examined with divergent primers of Chr5:2512798|2514806 (Supplemental Table S5). Full length of circR5g05160 was validated by Sanger sequencing.

A. tumefaciens strain EHA105 containing plasmid p35S-circR5g05160 was subjected to genetic transformation in rice NPB. Positive transgenic lines were confirmed by a hygromycin sensitivity test as described previously (Li et al., 2014). Transgenic plants overexpressing circR5g05160 (circR5g05160OX) were further confirmed by qPCR with the divergent primers of Chr5:2512798|2514806. RT-qPCR was also performed to examine the accumulation of linear transcripts from the parent gene LOC_Os05g05160 with the convergent primers located out of the circR5g05160 region. Primer sequences are listed in Supplemental Table S5.

Blast Disease Assay

Fully expanded leaves of 3-week-old seedlings were challenged with spores (2 × 10⁵ spore mL⁻¹) of M. oryzae eGFP-tagged Zhong8-10-14 (GZ8) through punch inoculation or spraying inoculation methods (Park et al., 2012; Li et al., 2019). The disease phenotype was photographed at 7 d postinoculation. Disease lesions were measured with ImageJ software. Relative fungal biomass was quantified as previously described (Park et al., 2012; Li et al., 2019).

Accession Numbers

The raw sequencing data were submitted to NCBI’s Gene Expression Omnibus and are accessible through GEO series accession number GSE131641.

Supplemental Data

The following supplemental materials are available.

• Supplemental Figure S1. Analysis of RNA-seq data from leaves of IR25 and LTH.

• Supplemental Figure S2. GO enrichment analysis of DEGs between IR25 and LTH.

• Supplemental Figure S3. Successfully amplified and sequenced circRNAs.

• Supplemental Figure S4. KEGG analysis on parent genes of circRNAs identified in IR25 and LTH.

• Supplemental Figure S5. The distribution of the number of splice sites in different junctions of circRNAs identified in IR25 and LTH.

• Supplemental Figure S6. Expression levels of circRNAs in IR25 and LTH.

• Supplemental Figure S7. RT-qPCR validation of differentially expressed circRNAs in IR25 and LTH.

• Supplemental Figure S8. Full-length validation and overexpression strategy for circR5g05160.

• Supplemental Table S1. Summary of circRNA-seq data from leaves of blast-resistant and -susceptible rice accessions with or without M. oryzae infection.

• Supplemental Table S2. High-confidence circRNAs identified in all samples.

• Supplemental Table S3. The 636 circRNAs specifically triggered by M. oryzae infection in IR25 and LTH.

• Supplemental Table S4. KEGG enrichment analysis of parent genes producing circRNAs in IR25 and LTH only upon M. oryzae infection.

• Supplemental Table S5. Primers used in this study.

• Supplemental Table S6. circRNAs detected in at least two samples and the corresponding parent genes indicated in Figure 2F.

• Supplemental Table S7. PCU of alternative 5′ and 3′ back-splicing sites.

• Supplemental Table S8. Differentially expressed circRNAs in IR25 and LTH upon M. oryzae infection.

ACKNOWLEDGMENTS

We thank Dr. Cai-Lin Lei (Institute of Crop Science, Chinese Academy of Agricultural Sciences) for providing the seeds of Pikm-monogenic line IR25 and Dr. Li-Huang Zhu (Institute of Genetics and Developmental Biology, Chinese Academy of Sciences) for providing the M. oryzae strain Zhong8-10-14.

LITERATURE CITED

Author notes