Xoo-responsive transcriptome reveals the role of the circular RNA133 in disease resistance by regulating expression of OsARAB in rice

Abstract

Upon Xanthomonas oryzae pv. oryzae (Xoo) infection of rice leaves, the invasion induces systematic expression changes for both the coding genes and the non-coding genes, allowing the plant to make corresponding responses. However, the roles of circular RNAs (circRNAs) in rice defending against Xoo remain largely unknown. To address this question, we conducted a whole-transcriptomic analysis to systematically screen the differentially expressed (DE) mRNAs and non-coding RNAs (ncRNAs) in rice responding to Xoo infection. Our results revealed a total of 4076 DE mRNAs, 89 DE long non-coding RNAs (lncRNAs), 82 DE microRNAs (miRNAs), and 14 DE circRNAs identified from Xoo-infected rice plants at 48 h post inoculation. Three circRNAs (ciR52, ciR298, and ciR133) were found to be able to form circular RNAs, and their expression was induced by Xoo infection. ciR133 was found to repress the expression of its parental gene OsARAB (putative arabinofuranosidase gene) during Xoo infection. Overexpression of ciR133 and mutation of OsARAB enhanced rice resistance against Xoo, without compromising main agronomic traits. Our data suggest that circRNAs are associated with rice response to Xoo infection, providing a potential strategy for breeding Xoo-resistant rice plants by manipulating ciR133 and OsARAB.

Background

Oryza sativa L. (Rice) is one of the most important crops, feeding more than half of the world’s population. Xanthomonas oryzae pv. oryzae (Xoo) is a devastating bacterial disease on rice, causing yield loss up to 50% in Asia and Africa (Adhikari et al. 1995; Ahmed et al. 2020). Xoo first enters rice leaves through wounds or margins, and then moves to the leaf xylem, leading to the plant wilting and yield loss (Wang et al. 2017; Oliva et al. 2019). Transferring resistance (R) genes into different rice cultivars has been demonstrated to be the most effective and economical strategy to control Xoo infection (Sun et al. 2004; Buddhachat et al. 2022). To date, more than 40 R genes have been identified in rice, and 14 of which have been cloned, including Xa1 (Yoshimura et al. 1998), Xa21 (Pruitt et al. 2015), and Xa23 (Wang et al. 2015).

Endogenous non-coding RNAs (ncRNAs) have been implicated in the interaction of rice resistance to Xoo infection (Lu et al. 2021; Song et al. 2021). Overexpression of miR169o in rice reduced the disease resistance to Xoo by decreasing the transcription of defense genes (Yu et al. 2018). miR1861k may play a role in broad spectrum resistance to Xoo (Zhang et al. 2015). Down-regulation of miR156 and overexpression of its two targets IPA1 and OsSPL7 enhanced rice resistance to Xoo infection (Liu et al. 2019). The Xoo-induced long ncRNA (lncRNA) ALEX1 up-regulates the expression of jasmonate (JA)-responsive genes, such as JAZ8 and MYC2, and mediates the increased endogenous levels of JA to enhance resistance against Xoo (Yu et al. 2020).

The covalently closed circular RNAs (circRNAs) were first reported in virus in 1975 (Takahashi and Diener 1975), and later they have been identified in almost all eukaryotes and prokaryotes (Li et al. 2018). Advances in next-generation sequencing have enabled the identification of circRNAs in Arabidopsis, rice, soybean, wheat, maize, and many other plant species (Lai et al. 2018). circRNAs are formed by canonical back-splicing of gene introns or exons from precursor mRNA, and their biogenesis is catalyzed and modulated by RNA binding proteins and flanking intronic complementary sequences (Li et al. 2018). Because of their special topological structure, circRNAs are more stable than other RNAs (Chen et al. 2019a). Generally, circRNAs are expressed at lower levels than linear RNAs or other ncRNAs, but they show more variants due to alternative splicing (Feng et al. 2019). Although the functions of most circRNAs remain largely unexplored, they are known to sequestrate miRNAs or proteins, and modulate gene transcription. The first reported functional circRNA in eukaryotes is ciRS-7, which was identified to act as a miRNA sponge to inhibit miR-7 function in human brains (Hansen et al. 2011; Memczak et al. 2013). Circ-MALAT1 in mammalian cancer stem cells was observed to function as a PAX5 mRNA translation brake by competing with the ribosomes entry site (Chen et al. 2019a). In Arabidopsis, a circRNA derived from exon 6 of the SEPALLATA3 (SEP3) gene was demonstrated to pause transcription of its cognate DNA locus by forming an RNA:DNA or R-loop (Conn et al. 2017). In rice, CircR5g05160 was found to be induced by Magnaporthe oryzae infection, which consequently promoted the immunity in rice (Fan et al. 2020). However, the response of circRNAs to Xoo infection in rice is not known yet.

Plants must protect their cell walls from penetration by bacterial and fungal pathogens that seek access to the nutrient-rich environment within the cell. The plant cell wall is dynamic, and its structure and composition are changed during pathogen infection, partially due to modification by microbial cell wall-degrading enzymes (Zhao and Dixon 2014). Plants have developed a dedicated mechanism to maintain cell wall integrity which could be utilized to enhance crop disease resistance (Bacete et al. 2018). Primary cell walls are composed of three polysaccharides: cellulose, hemicelluloses, and pectin, which determine the cell shape and mechanical strength. Lignin, a phenylpropanoid polymer, is deposited in the cell wall during secondary wall thickening (Brown et al. 2005). Arabinofuranosidases (ARABs) are enzymes that catalyze the hydrolysis of terminal, non-reducing α-L-arabinofuranoside residues, and contribute to reduce the number of links between lignin and cell wall polysaccharides, and to degrade lignin-carbohydrate complexes (Hashimoto et al. 2011). Transgenic rice plants expressing ARAB (ARAF) from Coprinopsis cinerea displayed the reduced ferulate-lignin cross-links by detaching arabinose side chains from arabinoxylan and the increased abundance of alkaline-resistant benzyl ether cross-links (Maruyama et al. 2021). Overexpression of rice OsARAF1 and OsARAF3 led to an increase in cellulose accumulation and saccharification efficiency, including an increase in arabinose content and a decrease in glucose content; however, it does not affect modifications of cell wall polysaccharides (Sumiyoshi et al. 2013). Polygalacturonase depolymerizes pectin by hydrolysis, thereby altering pectin composition and structure and activating cell wall defense. Loss-of-function mutation of ltn-212 destroyed cell wall integrity, resulting in spontaneous cell death and an auto-activated defense responses, including reactive oxygen species burst and pathogenesis-related gene expression, as well as enhanced resistance to Xoo infection (Cao et al. 2021). Knocking out of transcription factor WRKY53 thickens sclerenchyma cell walls of the vascular bundle, conferring the resistance to Xoo (Xie et al. 2021).

In the present study, our transcriptomic analysis revealed that a total of 4076 mRNAs, 82 miRNAs, 89 lncRNAs, and 14 circRNAs displayed differential expression in rice seedlings after 48 h of Xoo infection. An analysis for transgenic rice showed that ciR133 could repress expression of its parental gene OsARAB, and plays a positive role in defending against Xoo infection. Therefore, ciR133 and OsARAB may be the candidate genes that can be used for Xoo-resistance breeding in rice plants.

Results

Transcriptomic analysis of rice plants infected by Xoo

One-month-old rice leaves were treated with Xoo strain AXO1947 and H₂O, and were harvested two days later for transcriptomic sequencing (Fig. 1a). Pair-end RNA-Sequencing was used to identify mRNAs, lncRNAs, and circRNAs, while small RNA-Sequencing was used to identify miRNAs (Fig. 1b). Quality control of the 6 RNA samples yielded 108.80 Gb of clean data, with Q30 percentages higher than 93.47% and GC contents ranging from 43.86% to 45.25%. Subsequently, 740 million clean reads were generated and mapped to the rice genome (IRGSP-1.0), resulting in unique alignment rates ranging from 27.79% to 37.11% for total RNA-seq (Additional file 1: Table S1). To validate the RNA-seq data, real-time quantitative PCR (RT-qPCR) was used to examine the expression of five up-regulated and five down-regulated genes on differentially expressed (DE) mRNAs, and non-coding RNAs (lncRNAs and miRNAs) (Additional file 2: Figure S1). The results showed that the DE mRNAs, lncRNAs, and miRNAs identified from the RNA-sequencing could be verified by RT-qPCR.

Fig. 1

Flow chart of whole transcriptomic RNA-sequencing. a The workflow of wild type rice ZH11 infected by Xanthomonas oryzae (Xoo African strain AXO1947). Thirty-day-old plants were treated with Xoo and H₂O for 48 h for RNA extraction. The experiment was repeated with three biological replicates. b The schematic diagram for the generation of two cDNA libraries, RNA sequencing, and analysis of differentially expressed (DE) genes

A total of 31,698 mRNAs, 1,060 lncRNAs, 691 miRNAs, and 1,341 circRNAs candidates were expressed in the rice seedlings treated wtih Xoo (Additional file 1: Table S2). All of them are listed in Additional file 1: Tables S3–S6. Of them, 1358 new genes are sequenced (Additional file 1: Table S3), and 253 new miRNAs are found (novel miRs in Additional file 1: Table S5). After filtering the protein-coding RNAs (mRNAs) based on prediction by multiple software, we predicted that the majority of lncRNAs (624 lncRNAs) were intergene-lncRNAs, while the rest were antisense-lncRNAs, intronic-lncRNAs, and sense-lncRNAs.

The transcriptome of Xoo-infected rice plants

Analysis of DE mRNAs and ncRNAs after Xoo infection revealed a total of 4076 mRNAs, 89 lncRNAs, 82 miRNAs, and 14 circRNAs with altered expression, including 2041, 52, 39, and 10 up-regulated mRNAs, lncRNAs, miRNAs, and circRNAs respectively, and 2035, 37, 43, and 4 down-regulated mRNAs, respectively (Additional file 1: Table S2). All of the DE RNAs are listed in Additional file 1: Tables S7–S10.

Combined analysis of DE lncRNAs with DE miRNAs and mRNAs revealed that 78 DE lncRNAs co-expressed with the target DE Cis-mRNAs, 4 DE miRNAs co-expressed with their target DE lncRNAs, and three DE lncRNAs (MSTRG.26116.1, MSTRG.38793.1, and MSTRG.9627.3) co-expressed with the target DE Cis-mRNAs and DE miRNAs simultaneously (Additional file 1: Table S11 and Additional file 2: Figure S2a). All 89 DE lncRNAs co-expressed with target DE trans-mRNAs, and 4 DE lncRNAs (MSTRG.11997.2, MSTRG.26116.1, MSTRG.38793.1, and MSTRG.9627.3) co-expressed with target DE trans-mRNAs and DE miRNAs simultaneously (Additional file 2: Figure S2b and Additional file 1: Table S11). Based on the targeted genes of DE lncRNAs and DE miRNAs, the KEGG annotation suggests five categories, including cellular process, environmental information processing, genetic information processing, metabolism, and organismal systems (Additional file 2: Figures S2c, S3).

The DE mRNAs include some known pathogenesis-related genes, such as PR10a and OsSWEET14 (Additional file 1: Table S12), while the DE miRNAs include rice immunity related miRNAs, such as miR169a and miR528 (Additional file 1: Table S13) which have been reported in pathogen-infected rice (Li et al. 2014). These data indicate that Xoo infection can cause significant changes in gene expression. The combined analysis of DE miRNAs and DE mRNAs, among the 82 DE miRNAs and 4,076 DE mRNAs after Xoo infection, indicated that there are 33 down-regulated miRNAs and their putative target genes are up-regulated; however, there are 7 up-regulated miRNAs and their putative target genes are down-regulated (Table 1). This result suggests that these miRNAs play important roles in response to Xoo infection by regulating the putative target genes.

Table 1 All differentially expressed (DE) miRNAs in rice ZH11 plants by combined analysis with the eligible predicted target DE mRNAs after 48 h of Xoo infection

Expression change of circRNAs after Xoo infection

A total of 1341 circRNA candidates were identified from rice seedlings using the find_circ and CIRI2 software packages (Fig. 2a and Additional file 1: Tables S2, S6). The majority (85%) of the parental genes of circRNAs produced only one circRNA, while some generated two to eight (Fig. 2b). The length of circRNAs ranged from 200 to over 3000 nucleotides (nt). Exonic circRNAs were mostly abundant in 400 nt, while a few intergenic circRNAs were longer than 3000 nt (Fig. 2c). The average expression level of circRNAs before Xoo infection was higher than that after 48 h of Xoo infection (Fig. 2d).

Fig. 2

Summary of circRNAs in ZH11 seedlings after Xoo infection. a Identification of differentially expressed (DE) circRNAs after Xoo infection. A total of 1,341 circRNAs were detected, and 14 of which showed differential expression (FC > 2 and P < 0.05) after Xoo infection. Seven of the DE circRNAs were predicted to derive from the known functional parental genes. b Number of genes capable of producing 1–8 circRNAs. c Number of circRNAs with different lengths and circRNAs are produced from different locations of genes. d Average expression level (transcripts per million, TPM) of all circRNAs in ZH11 and Xoo-infected ZH11

The find_circ and CIRI software were used to analyze the DE circRNAs in response to Xoo-infection. A total of 14 DE circRNAs were found from rice seedlings after 48 h of Xoo infection (Fig. 2a, b and Table 2). Only one up-regulated and one down-regulated circRNA were identified by two software, which are chr03:11550445|11551137 (ciR133), and chr06:16279995|1680386 (Fig. 3a). Of the 14 DE circRNAs, and seven genes were generated from parental genes, namely Os03g0320100 (a-L-arabinofuranosidase), Os04g0684500 (OsWSL5), Os04g0686650, Os06g0480200, Os07g0672500, Os09g0268300 (Hexose carrier), and Os09g0343200. circRNAs and other noncoding RNAs cannot function independently in plants, so we further studied and analyzed the relationship between circRNAs, DE miRNAs, and DE mRNAs through bioinformatics (Fig. 3c). We found that 225 circRNAs potentially target DE mRNAs, and 55 circRNAs have putative miRNA binding sites. Further analysis predicted that three DE circRNAs may interact with DE miRNAs (chr4:2014063|24014278, chr4:31276588|31301714, chr8:27930611|27963872) to form a DE circRNA-DE miRNA network, and seven DE circRNAs may interact directly with DE mRNAs (3:11550445|11551137, 4:35093089|35095669, 9:10618503|10619042, 9:5140146|5140580, 7:28426090|28427266, 6:16279995|16280386, 4:34965141|34969807) to form a DE circRNA-DE mRNA network.

Table 2 Differentially expressed circRNAs in rice ZH11 plants after 48 h of Xoo infection

Fig. 3

Validation of the differentially expressed (DE) circRNAs in rice plants after Xoo-infection. a Number of DE circRNAs after Xoo-infection were identified by find_circ and CIRI2 software. b Number of the up- and down-expressed circRNAs after Xoo-infection. c Number of the potential network prediction of circRNA-miRNA-mRNA. The DE circRNAs were used for predicting the target mRNAs and miRNAs. d–f Validation of circular RNA formation of 3 circRNAs by (RT)-PCR assays. d PCR diagram using genomic DNA (left) and circular cDNA (right) of circRNAs as template. e The electrophoretogram of circular RNA formation of 3 circRNAs by (RT)-PCR assays. The circular cDNA of circRNAs (ciR52, ciR298, and ciR133) could be amplified by divergent primers (first column), but not their genomic DNA (gDNA) (third column). Circular cDNA from circRNAs and gDNA of the parental genes could be amplified by the convergent primers (second and fourth columns). f Expression change of the DE circRNAs in ZH11 after Xoo infection using divergent primers in RT-qPCR assays. OseEF1-a was used as a reference gene. Means ± SD represent the value with standard deviation (n = 3 biological repeats)

Because the divergent primer pair for RT-PCR can amplify bands from the circular RNAs but not from linear RNAs and the convergent primer pair can detect the linear RNAs but not circular RNAs (Fig. 3d), these primer pairs were used to distinguish circRNAs from linear RNAs or linear DNAs by PCR and RT-qPCR (Fig. 3e, f). Three circular RNAs (ciR52, ciR298, and ciR133) could be amplified using the divergent primer pair when cDNAs were used as templates, but not gDNAs as templates. However, their parental genes can be amplified with the convergent primer pair when cDNAs and gDNAs were used as templates (Fig. 3e). The result showed that ciR52, ciR298, and ciR133 can form circular RNAs in rice plants. RT-qPCR was also performed with the convergent primer pair to examine the circRNAs responding to Xoo treatment (Fig. 3f). The result showed that the expression of the three circRNAs (ciR52, ciR298, and ciR133) could be induced by Xoo infection.

ciR133 is a circular RNA in rice plants

To further verify whether ciR133 is a circular RNA in rice, the back splicing junction site was first validated by Sanger sequencing (Fig. 4a). ciR133 is a 400 nt exonic circRNA, formed by the back-splicing of the 4th to 6th exon of its parental gene OsARAB (Os03g0320100/LOC_Os03g20420), which encodes a putative a-L-arabinofuranosidase (Fig. 4a).

Fig. 4

ciR133 is a non-coding circular RNA derived from OsARAB. a Schematic diagram of ciR133 RNAs formed by back-splicing from the 4th, 5th and 6th exons of its parental gene OsARAB. b Validation of circular RNA formation of ciR133 by RT-PCR assays. Total RNAs of rice seedlings were treated with RNase R, and then used for the RT-PCR assays. c RT-qPCR analysis for ciR133 expression level in the Xoo infected rice leaves. RT for total RNAs were performed with random or oligo(dT) primers. Os-eEF1-1a was used as a reference gene. Means ± SD represents the value with standard deviation (n = 3 biological repeats, **P ≤ 0.01)

Then, the circularity of ciR133 was further confirmed by RT-PCR using the divergent primer or the convergent primer (Fig. 4b). RNase R is known to degrade linear RNAs but not circRNAs, total RNAs were digested with RNase R for reverse transcription (RT) to further confirm circular ciR133. Therefore, only the circular ciR133 can be amplified using the RNase R-digested RNAs as template and the divergent primers for RT-PCR (Fig. 4b). However, even using the convergent primers for RT-PCR, the linear mRNA of OsARAB cannot be amplified using the RNase R-digested RNAs as a template (Fig. 4b).

Thirdly, because the circular RNAs lack poly (A) tail, the efficiency of RT using oligo (dT) as RT primer is lower than that using random primer. Therefore, RT-qPCR results showed that the expression level of the cDNA reverse-transcribed ciR133 by random primer was significantly higher than that using the cDNA reverse-transcribed ciR133 by oligo (dT) primers (Fig. 4c). This result indicates that ciR133 lacks a poly (A) tail which is usually present in protein coding genes. The above three results demonstrate that ciR133 can form a circular RNA in rice plants.

ciR133 represses the expression of its parental gene OsARAB

Does ciR133 regulate the expression of its parental gene OsARAB? We first analyzed the expression of ciR133 and OsARAB during Xoo infection (Fig. 5a). The results showed that the expression of ciR133 was induced, but the expression of OsARAB was repressed by Xoo (Fig. 5a), which indicates that ciR133 may repress expression of its parental gene OsARAB during Xoo infection. We then used a dual-luciferase system to analyze the interaction of ciR133 and OsARAB mRNA in tobacco plants (Fig. 5b). The fluorescence value of ciR133 and OsARAB co-transformed tobacco leaves was significantly lower than that of empty vector and OsARAB co-transformed tobacco leaves, indicating that ciR133 could negatively regulate expression of OsARAB. To further confirm this observation, we constructed the transgenic rice plants with over-expression (ciR133-ox) and down-expression (Si-ciR133) of circular ciR133, and with overexpression of its corresponding linear 133 (Linear133-ox) (Fig. 5b and Additional file 2: Figure S4). Overexpression of circular ciR133 in the ciR133-ox rice (Fig. 5c) resulted in repression of OsARAB expression (Fig. 5d). Conversely, repression of circular ciR133 in the Si-ciR133 rice (Fig. 5c) led to up-regulation of OsARAB expression (Fig. 5d). However, overexpression of linear ciR133 in Linear133-ox (Additional file 2: Figure S5b) did not alter expression of ciR133 (Fig. 5c) or OsARAB (Fig. 5d). These results indicated that only circular ciR133 can repress the expression of its parental gene OsARAB.

Fig. 5

ciR133 represses expression of OsARAB and responds to Xoo infection. a Expression change of ciR133 and its parental gene OsARAB in four weeks old ZH11 plants after Xoo (AXO1947)-infection. b Dual-luciferase assay of interaction between ciR133 and OsARAB in tobacco leaves. Schematic diagram of interaction construction was displayed in upper panel. Interaction of ciR133 and OsARAB mRNA was displayed in left panel, and the luciferase activity was shown in right panel. N. benthamiana were grown at 30°C. The plants were infected for two days at 25°C. c Expression level of circular ciR133 in the transgenic rice lines with overexpression of circular ciR133 (ciR133-ox), down-expression of circular ciR133 (Si-ciR133), and overexpression of linear ciR133 (Linear133-ox). To confirm whether the circular ciR133 RNAs were altered in these transgenic rice seedlings (ciR133-ox and Si-ciR133), their total RNAs were digested with RNase R for RT, and the divergent primers were used for qPCR assays. Os-eEF1-1a was used as a reference gene. d Expression change of OsARAB in ciR133-ox, Si-ciR133, and Linear133-ox seedlings. Os-eEF1-1a was used as a reference gene. Means ± SD represents the value with standard deviation (n = 3 biological repeats)

Overexpression of ciR133 or knocking out of OsARAB enhanced rice resistance against Xoo infection

The expression of ciR133 was highly induced by Xoo among the 3 DE circRNAs (Figs. 3f and 5a), prompting us to investigate the role of ciR133 that is derived from OsARAB in rice-Xoo interaction. To this end, the transgenic rice lines with overexpression (OsARAB-ox) and mutation (Cri-ARAB) of OsARAB were generated (Additional file 2: Figure S5a, c). These transgenic rice lines were then selected for Xoo infection assays (Fig. 6). The results showed that the bacterial lesion lengths of ciR133-ox and Cri-ARAB transgenic rice plants were shorter than that of ZH11 (Fig. 6a, b). However, Linear133-ox, Si-ciR133, and OsARAB-ox transgenic rice plants did not show any differences in lesion length compared to ZH11. This result indicated that overexpression of ciR133 or mutation of OsARAB could enhance rice resistance to Xoo infection. Nevertheless, these transgenic rice plants did not show any differences in major agronomic traits with ZH11 (Additional file 2: Figure S6). To further elucidate how ciR133 is involved in Xoo resistance, we examined the expression of some defense-related genes. Three genes, the OsPR10a, OsPR10b, and NAC4 that were reported by Li et al. (2014), were up-regulated in ciR133-ox and Cri-ARAB plants (Fig. 6c–e).

Fig. 6

ciR133 and OsARAB mediate disease resistance against Xoo infection in rice. a, b Phenotype (a) and lesion lengths (b) of the ciR133 and OsARAB transgenic rice leaves after Xoo-infection. The 1.5 month-old rice leaves were inoculated with Xoo strain AXO1947 for 14 days. Scale bar = 1 cm in a. Means ± SD in b represents value with standard deviation (n = 3 biological repeats), **P < 0.01 in Student’s t-test (n > 15). c–e Relative expression of the pathogenesis related genes (OsPR10a, OsPR10b, and NAC4) in the transgenic rice plants. OseEF1-a was used as a reference gene. Means ± SD represent the value with standard deviation (n = 3 biological repeats)

Discussion

Transcriptomic assay of rice leaves in response to Xoo infection

Bacterial leaf blight can cause a significant reduction to rice yield. Both mRNAs and ncRNAs function in regulating plant growth, development, and response to biotic and abiotic stresses. In this study, the rice variety ZH11 was monitored for responsive changes in transcriptome, including mRNAs and ncRNAs after infection of the Xoo strain AXO1947, which resulted in the identification of 31,698 protein-coding genes, 1341 circRNAs, 1060 lncRNAs, and 691 miRNAs in the rice plants (Additional file 1: Tables S2–S6). Of these, 4076 DE genes, 89 DE lncRNAs, 82 DE sRNAs, and 14 DE circRNAs were identified in Xoo-infected rice leaves (Additional file 1: Tables S2, S7–S10). The accuracy of the RNA-seq data was verified by the previously reported DE mRNAs (OsSWEET14 and OsPR10a) (Additional file 1: Table S12) and DE miRNAs (osa-miR169a and osa-miR172) (Additional file 1: Table S13) (Yu et al. 2018; Oliva et al. 2019), and confirmed by the results of some DE gene expression (Additional file 2: Figure S1). Combined analysis with all the DE mRNAs, circRNAs, lncRNAs, and miRNAs indicate that there was no co-targeted DE mRNAs of DE lncRNAs, miRNAs, and circRNAs, but there are 7 DE lncRNAs with the co-targeted DE mRNA and DE miRNAs, such as Os06g0712200 which has been identified in a Xoo infection analysis (http://systbio.cau.edu.cn).

Transcriptional analysis of circRNAs in response to Xoo infection

circRNAs have been found to be involved in various stages of plant development, stress responses, and immune responses (Chen et al. 2019b). According to PLANTCIRCBASE v.7.0 (Xu et al. 2022), rice has been bioinformatically predicted to yield > 40 k circRNAs, and RNA-seq was used to assembly the full-length sequences of 6519 circRNAs in rice seedlings (Chu et al. 2022). Additionally, the find_circ toolkit was used to predict 12,037 circRNAs in rice roots (Ye et al. 2015). Two years later, the same group used TopHat-Fusion software to predict 1046, 1328, and 1225 full-length sequences of circRNAs from three different public datasets (Ye et al. 2017). CIRI2 was used to observe 5723 and 3480 circRNAs in rice varieties IR25 and LTH, respectively (Fan et al. 2020). In IR25 young leaves, Fan et al. (2020) identified 2932 circRNAs and found that 636 circRNAs were specifically detected upon M. oryzae infection. Here, our results predicted 1341 circRNAs in rice seedlings (Additional file 1: Tables S2, S6), and 14 circRNAs showed differential expression after Xoo infection (Additional file 1: Tables S2, S10). Three DE circRNAs (ciR52, ciR298, and ciR133) were confirmed to form circular RNAs and respond to Xoo infection (Fig. 3d, e). These results imply that some circular RNAs play important roles in rice response to Xoo infection.

Circular ciR133 may affect cell wall to attenuate Xoo infection

ciR133 is an exonic circRNA derived from its parental gene OsARAB (Fig. 4a), which is predicted to encode an α-L-arabinofuranosidase, an enzyme that is involved in the degradation of lignocelluloses (Mondher and Narayan 2006). Previous studies have reported a certain extent of correlation between circRNAs and their parental genes (Pan et al. 2018). Although numerous circRNAs have been identified in plants by genome-wide analysis, few have been functionally characterized. Xoo infection induced expression of ciR133 and repressed expression of its parental gene OsARAB (Fig. 5a). Overexpression of ciR133 enhanced rice resistance against Xoo infection, and knocking out of OsARAB also increased rice resistance to Xoo infection (Fig. 6). This result suggests that ciR133 mediates disease resistance to Xoo infection through repression of OsARAB in rice, thereby maintaining cell wall integrity during Xoo infection (Sumiyoshi et al. 2013; Bacete et al. 2018; Drula et al. 2022). Therefore, we proposed that Xoo-induced ciR133 expression may repress expression of OsARAB. As a result, ciR133-overexpression and the knocking out of OsARAB could reinforce cell wall integrity in response to Xoo infection. The arabinofuranosidase is also considered as a virulence factor of Sclerotinia sclerotiorum to cause canola necrosis in leaves and stems, and disruption of the encoding gene of its arabinofuranosidase reduced S. sclerotiorum virulence on canola tissue (Yajima et al. 2009). Likewise, Xoo-induced ciR133 may also repress the Xoo-arabionofuranosidase expression and reduced its virulence on rice. In addition, three pathogenesis-related genes (OsPR10a, OsPR10b, and NAC4) (Li et al. 2014) were induced in ciR133-ox and Cri-ARAB plants (Fig. 6c–e), which may also contribute to the resistance to Xoo.

Methods

Plant materials and growth conditions

Rice cultivar Zhonghua 11 (Oryza sativa L. ssp. Japonica, ZH11) was used in this study. Seeds were germinated in the dark for 3 days, then the germinated seeds were planted in pots or grown under optimum conditions in the green house in South China Botanical Garden, Chinese Academy of Sciences, under natural light (12 h light/12 h darkness) and temperature (25–35°C). Four weeks old plants were used for Xoo inoculation.

Bacterial blight strains and inoculation assays

Virulent Xoo strain AXO1947 that is obtained from Africa (Gu et al. 2004) were cultivated on PSA medium (10 g/L peptone, 10 g/L sucrose, 1 g/L glutamic acid, 16 g/L agar, pH 7.0) at 28°C under darkness for 2 days. Bacteria were collected and re-suspended in sterile water at an OD₆₀₀ = 0.5 for inoculation of rice leaves. The bacteria were injected into the leaves with a syringe. The inoculated rice leaves were harvested 48 h later and were immediately frozen in liquid nitrogen, then stored at –80°C for further analysis. The leaves with the mock treatment were harvested and used as a control. Each treatment has three replicates. Bacterial blight inoculation was performed with the leaf-clipping method as previously described (Zeng et al. 2020). Bacterial lesion lengths were measured at two weeks post inoculation.

Total RNA isolation, library preparation, and RNA sequencing

A total of 0.5 g leaves was taken from three plants and mixed together as one sample. The total RNAs were isolated and extracted using RNA TRIzol DP411 (Tiangen, Beijing, China), followed by the manufacturer’s instructions. Quantity of total RNAs were conducted by NanoDrop 2000 (Thermo Fisher Scientific, New York, USA) and its integrity was detected by LabChip GX (Agient2100, CA, USA) with RIN ≥ 6.5. High quality of total RNAs from leaves were mixed and sent to BioMarker Technologies (Beijing, China) for library construction and sequencing.

Library preparation for mRNA, lncRNA, and circRNA sequencing

For mRNA, lncRNA, and circRNA sequencing, about 1.5 μg RNA per sample was used as input and the rRNA was removed using the Ribo-Zero rRNA Removal Kit (Epicentre, Madison, WI, USA). Sequencing libraries were generated using NEBNextR UltraTM Directional RNA Library Prep Kit for IlluminaR (NEB, USA). Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5×). First strand cDNA was synthesized using random hexamer primer and reverse transcriptase. Second-strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3’ ends of DNA fragments, NEBNext adaptor with hairpin loop structure were ligated to prepare for hybridization. In order to select insert fragments of preferentially 150 ~ 200 bp in length, the library fragments were purified with AMPure XP Beads (Beckman Coulter, Beverly, USA). Then 3 μL USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37°C for 15 min before PCR. Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. Finally, PCR products were purified with AMPure XP system (Beckman Coulter, Suzhou, China) and library quality was assessed on the Agilent Bioanalyzer 2100 (Agilent, Berlin, German).

Library preparation for sRNA sequencing

A total of 2.5 ng RNA per sample was used as input material for the RNA sample preparation. Sequencing libraries were generated using SRNA Library Construction SOP by BMK-2020 version (BioMarker Technologies, Beijing, China). First, the 3’SR Adaptor was ligated with the RNA. Mixed 3’SR adaptor, RNA and nuclease-free water were incubated for 2 min at 70°C in a preheated thermal cycler. The mixture was then transferred on ice. Then, 3’Ligation Reaction Buffer (2×) and 3’Ligation Enzyme Mix were added to ligate the 3’SR adaptor. In order to prevent the formation of adaptor-dimer, the SR RT Primer was utilized to hybridize with the excess of 3’SR adaptor that remains unbound after the 3’ligation reaction. This hybridization process transforms the single-stranded DNA adaptor into a double-stranded DNA molecule. The double-stranded DNAs are not suitable substrates for ligation-mediated reactions. Then the ligation of the 5’SR Adaptor, followed by the synthesis of the first chain through reverse transcription. Finally, PCR amplification and size selection were performed. For fragment screening, PAGE gel electrophoresis was employed, and rubber cutting was used to recycle the pieces that contained small RNA libraries. Finally, PCR products were purified using AMPure XP system and library quality was assessed on the Agilent Bioanalyzer 2100 system.

Clustering and sequencing

Clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kitv3-cBot-HS (Illumina) according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina platform and reads were generated.

Reads mapping and transcriptomic assembly

Quality control was performed using FastQC (de Sena Brandine and Smith 1874). RNA-seq paired-end reads were aligned to the rice IRGSP-1.0 genome (RAP-DB|HOME (affrc.go.jp) via TopHat (Trapnell et al. 2009) and assembled via Cufflinks (Trapnell et al. 2010). Prediction and identification of mRNAs, lncRNAs, and circRNAs were performed with visual pipeline described in Fig. 1b. The transcriptome was assembled using the StringTie (https://ccb.jhu.edu/software/stringtie/index.shtml), based on the reads that are mapped to the reference genome. The assembled transcripts were annotated using the gffcompare program (http://ccb.jhu.edu/software/stringtie/gffcompare.shtml).

The unknown transcripts were used to screen putative lncRNAs. Three computational approaches including CPC/Pfam/CPAT were combined to sort non-protein coding RNA candidates from putative protein-coding RNAs in the unknown transcripts (Bateman et al. 2004; Kong et al. 2007; Wang et al. 2013). Putative protein-coding RNAs were filtered using minimum length and exon number threshold. Genes with lengths more than 200 nt and contain more than two exons were selected as lncRNA candidates and were further screened using CPC/Pfam/CPAT that can distinguish the protein-coding genes from the non-coding genes. The different types of lncRNAs, including lincRNA, intronic lncRNA, anti-sense lncRNA, sense lncRNA, were selected using cuffcompare program. The unmapped reads were used to predict circRNAs.

We used CIRI2 (CircRNA Identifier) tools (Gao et al. 2015) and find_circ (Memczak et al. 2013) to identify circRNA. It scans SAM files twice and collects sufficient information to identify and characterize circRNAs. Briefly, during the first scan of SAM alignment, CIRI detects junction reads with PCC signals that reflect a circRNA candidate. Preliminary filtering is implemented using paired-end mapping (PEM) and GT-AG splicing signals for junctions. After clustering junction reads and recording each circRNA candidate, CIRI scans the SAM alignment again to detect additional junction reads and meanwhile performs further filtering to eliminate false positive candidates resulting from incorrectly mapped reads of homologous genes or repetitive sequences. Finally, the identified circRNAs are output with annotation information. We used TargetFinder tools (https://github.com/carringtonlab/TargetFinder) to predict target miRNA.

Small RNA-Seq data were analyzed using Bowtie (https://bowtie-bio.sourceforge.net/bowtie2/manual.shtml#command-line). The generated reads were subjected to sequence alignment against the Silva (https://bio.tools/silva#!), GtRNAdb (http://gtrnadb.ucsc.edu/), Rfam (https://rfam.org/), and Repbase databases (https://www.girinst.org/repbase/), followed by filtering of ribosomal RNA (rRNA), transfer RNA (tRNA), small nuclear RNA (snRNA), small nucleolar RNA (snoRNA), other non-coding RNA (ncRNA), and repeats. The remaining reads were used to detect known miRNAs and novel miRNAs predicted by comparing with known miRNAs from miRBase (https://www.mirbase.org/). Randfold tools soft (https://hpc.ilri.cgiar.org/randfold-software) was used for novel miRNA secondary structure prediction.

DE mRNAs and ncRNAs analysis

DESeq R package (1.10.1) (Robinson et al. 2010) was utilized to establish the expression levels of all transcripts. The P-value < 0.05 and log2 (fold change, FC) ≥ 1.0 were set as the cut-off criteria for identifying the DE-mRNAs, DE-miRNAs, DE-lncRNAs, and DE-circRNAs.

RT-qPCR

Total RNAs from leaves before and after inoculation with Xoo was extracted using a TRIzol kit (Invitrogen, USA), and were reverse transcribed into cDNA using an M-MLV RT Kit (Accurate Biotechnology, AG11711, China). For RT of mRNAs and lncRNAs, the oligo(dT) primers were used. For RT of miRNAs, the specific stem-loop primers were used. For RT of circRNAs, total RNAs were treated at 37°C with 3 U/µg of RNase R (Epicentre Biotechnologies, Madison, WI, USA) for 30 min, and the random 6-mer primers were used.

For qPCR, the divergent primers for detecting circRNAs and the convergent primers for detecting linear RNAs were designed using Primer 3 (version 0.4.0) (http://bioinfo.ut.ee/primer3-0.4.0/). The above cDNAs were used as templates for qPCR with SYBR qPCR Real Master Mix (ChamQ Universal, China) on LightCycler 480 (Roche, USA). OseEF1-a was used as the reference gene for mRNAs, lncRNAs, and circRNAs, and U6 was used as the reference gene for miRNAs.

Validation of the circular RNAs by RT-PCR assays

To verify circular RNAs, the total RNAs were treated with or without RNase R for RT assay, and the divergent primers and convergent primers were designed for PCR assays. The divergent primers of PCR could only amplify bands using the circular cDNAs as the template, but not using genomic DNAs (gDNAs) as template; while the convergent primers could amplify bands using both linear cDNAs and gDNAs as templates (primer pairs as in Fig. 3d). The PCR reaction was performed using 2 × Taq Plus Mix (Vazyme, China) according to the manufacturer’s protocol. Sanger sequencing of these PCR bands was performed after purifying PCR products using DNA Gel Extraction Kit (Tsingke Biotechnology, China) according to the manufacturer’s protocol.

Construction of transgenic rice with overexpression of circular and linear ciR133

Due to the circular nature of circRNAs, the construction methods of overexpression and interference vectors are slightly different from those of the coding genes. For overexpression of circular ciR133 and liner133, the vector construction was performed as following process. The full length of ciR133 and the target fragments from the start site of the upstream intron to the end site of the downstream intron of ciR133 were cloned using the specific primers. The genomic region of ciR133 was amplified from ZH11 gDNA, and the cognate linear fragment of ciR133 was amplified from ZH11 cDNA. The specific primer pairs for cloning ciR133 and linear fragment are both with restriction endonuclease of Bgl II and BstE II, respectively. The fragments were then cloned into pCAMBIA1301 vectors under the control of Cauliflower mosaic virus 35S promoter to generate the recombinant plasmids of pCAMBIA1301-ciR133 and pCAMBIA1301-35S:linear133, respectively (Additional file 2: Figure S4). The full length of ciR133 and its cognate linear fragment were validated by Sanger sequencing.

For RNAi of circular ciR133, the vector construction was performed as following process. The circRNA back-to-back splice site is the only difference from the parental gene transcript. RNAi technology can not only effectively inhibit the expression of circRNA, but also ensure that the parental gene transcription is not affected. Taking 20 nt before and after the back-to-back splicing site after ciR133 looping as the target sequence, the primers were designed and cloned into the intermediate vector pNW55 using the artificial microRNA technology (http://wmd3.weigelworld.org/cgi-bin/webapp.cgi), and then sub-cloned into vector pCAMBIA1301 at BglII and AflII digestion sites.

Agrobacterium tumefaciens strain EHA105 containing plasmids were transformed into ZH11. Rice transformation of ZH11was performed as previously described (Hiei et al. 1994). Positive transgenic lines were confirmed by a hygromycin sensitivity test as well as by amplifying the Hpt. Expression levels of circRNA133 in the transgenic rice lines that overexpress ciR133 (ciR133-ox) or linear ciR133 (35S:linear 133-ox) were further confirmed by RT-qPCR assays.

Mutation and overexpression of OsARAB in rice

Mutation of OsARAB (LOC_Os03g20420) were implemented using the CRISPR/Cas9 genome editing system (Ma et al. 2015). The anchor sequence was targeted to the 7th exon of OsARAB. gDNA templates of the transformed rice plants were extracted, and the editing site was amplified by PCR using a pair of flanking primers. The PCR products were sequenced and analyzed using DSDecode (Ma et al. 2015) to confirm OsARAB deletion. To further confirm the edited sequence, RT-PCR products from the mutated OsARAB seedlings were sequenced again. To construct OsARAB -overexpressing plants (OsARAB-ox), a 822 bp OsARAB cDNA fragment containing the open reading frame (ORF) was inserted downstream of the 35S promoter in pCambia1301 at AflII and BglII digestion sites to produce the plasmid p35S-OsARAB, and then the vectors were transferred into A. tumefaciens strain EHA105 to transform ZH11.

All the primers used in this study are listed in Additional file 1: Table S14. IBM SPSS 20.0 software was used to determine the statistical significance of the data.