Identication and Characterization of Rice Circular RNAs Responding to Xanthomonas Oryzae pv. Oryzae Invasion

The emerging role of circular RNAs (circRNAs) in various biological processes have advanced our knowledge of transcriptional and post-transcriptional gene regulation. The number and expression of plant circRNAs vary with species and treatments. However, the expression prole and the potential role of circRNAs during plant response to pathogen invasion are still elusive. In this study, we identied 3517 circRNAs from PXO99 A -infected rice leaves using the ribosomal RNA (rRNA) depleted RNA-Sequencing technique coupled with the CIRI2 and CIRCexplorer2 pipeline. Among them, 2994 (85.13%) circRNAs arised from the exons of their parent genes, 1214 circRNAs were previously unknown and 276 circRNAs exhibited differential expression proles upon PXO99 A infection over time. In addition, 31 differentially expressed circRNAs (DEcircRNAs) were predicted as the corresponding 121 miRNAs sponges. Functional analysis of both host genes and target mRNAs suggested that these identied circRNAs might play an important role in reprogramming rice responses to PXO99 A invasion, mainly by mediating photorespiration, chloroplast, peroxisome and diterpenoid biosynthesis associated pathways. These results inferred a potential provide novel clues PRRs: ncRNAs: Non-coding RNAs; DEcircRNAs: Differentially expressed circRNAs; RNA-seq: RNA sequencing; Hpi: Hour post inoculation; KEGG: Kyoto encyclopedia of genes and genomes; GO: Gene Ontology; RT-qPCR: Real-time quantitative PCR.

bioinformatics tools signi cantly facilitate the discovery of novel noncoding RNAs (ncRNAs) (Memczak et al. 2013;Ponting et al. 2009;Ravasi et al. 2006). Diverse types of ncRNAs with varied functions have been found in both unicellular and multicellular organisms, which have shifted the attention in the elds of molecular biology from protein-coding RNAs to ncRNAs (Morris and Mattick 2014).
As a novel class of ncRNAs, circular RNAs (circRNAs) has received increasing attention (Memczak et al. 2013). Structurally different from the linear RNA molecules, their 3' and 5' ends are covalently joined to form closed loops which arise from back-splicing of pre-mRNAs (Chen and Yang 2015) and can protect them from RNA exonuclease-mediated degradation (Cocquerelle et al. 1993;Jeck et al. 2013).
Characterization analyses show that circRNAs are mainly divided into exonic, exon-intronic, intronic, and intergenic region (Salzman et al. 2012;Chu et al. 2018) and have various isoforms elicited from a circRNA locus Chu et al. 2018). Over the past few years, thousands of circRNAs have been identi ed in animals. Their signi cance is being addressed frequently in literature for their role in transcriptional and post-transcriptional regulation, such as functioning as miRNA sponges, modi ers of their host genes expression, and regulators of parental genes transcription and splicing (Ebbesen et al. 2016;Qu et al. 2015). However, their function remains elusive in plants, although they have been widely identi ed in both monocot and dicot plants (Lu et al. 2015;Sun et al. 2016;Wang et al. 2017;Zhou et al. 2018;Zuo et al. 2016). Recently, there has been an increasing focus on characterizing their relevance and function in plant stress responses. For example, in tomato fruit under chilling stress, 163 circRNAs were identi ed with chilling-responsive expression, of which 102 contain miRNA-binding sites and have the potential to act as miRNA sponges (Zuo et al. 2016). Likewise, 584 circRNAs were differentially expressed in kiwifruit plants during canker pathogen infection, and the expression of particular circRNAs is altered upon the stage of infection . A total of 2932 circRNAs were identi ed in rice leaf transcriptomes of blast-resistant and -susceptible accession, of which 636 were detected speci cally subjected to Magnaporthe oryzae (M. oryzae) infection, responsible for circRNAs diversity in rice (Fan et al. 2020). These reports suggested that the number and the expression of plant circRNAs vary with treatments and plant species. However, the potential role of circRNAs in rice responding to Xoo invasionat least to our knowledge -is largely missing. Our study primarily focused on identi cation and characterization of circRNAs involved in rice-Xoo interaction, and uncovered their potential regulation function in plant responses to bacterial pathogen invasion.

Discovery and characteristic of circRNAs in rice leaves
OsPR1a and OsPR1b (pathogenesis-related proteins coding genes) were regarded as reliable marker genes of plant defense responses (Agrawal et al. 2001). To observe whether the transcriptional regulatory events occurred in rice-Xoo interactions, we checked the expression pro les of OsPR1a and OsPR1b by RT-qPCR using speci c primers (Table S6) at different time points post inoculation (hpi) with PXO99 A . The levels of the OsPR1a and OsPR1b transcripts showed a similar induction pro le ( Figure S1), which is consistent with the defense-responsive expression in previous report (Jia et al. 2020). Their transcript levels slightly induced at 2 hpi and 6 hpi or markedly increased at 12 hpi and decreased drastically at 24 hpi ( Figure S1). Combining the infection strategy and time points chosen in the study of pathogenresponsive lncRNAs in rice (Yu et al. 2020), 0 h (control), 2, 6, 12 and 24 h were chosen as the ideal sampling time points for the identi cation of circRNAs responding to Xoo in rice. 5 rRNA-depleted RNA libraries were constructed for 3-week-old seedlings of Nipponbare infected by PXO99 A . The deep-sequencing using Illumina Hiseq X-ten platform yielded 45057324, 47264368, 42590676, 41416705 and 51015354 raw reads from 0, 2, 6, 12 and 24 hpi, respectively, followed by a mapping to the rice reference genome (Table S1). After removing the redundant circRNAs by CD-HIT-EST, a total of 3517 circRNAs were identi ed using CIRCexplorer2 (589 circRNAs) and CIRI2 (2648 circRNAs) from all samples (Fig. 1A, Figure S2), only 280 circRNAs were shared by both algorithms ( Figure S2). Among them, 512, 486, 622, 621 and 733 were respectively identi ed at 0, 2, 6, 12 and 24 hpi (Fig. 1A).
According to their genomic loci, the identi ed 3517 circRNAs were grouped as exon, intergenic, antisense and intron regions. Most circRNAs (85.13%) were derived from exons, and a few from intergenic regions (13.48%) (Fig. 1B). Moreover, we compared the circRNAs detected in our study with the publicly available rice circRNAs from PlantcircBase. It revealed that 2303 (65.48%) of circRNAs in Nipponbare have already been identi ed, yet a lot of novel circRNAs (1214, 34.52%) were only detected in this study (Fig. 1C).

CircRNAs are differentially expressed in response to PXO99 A infection
To reveal circRNAs involved in rice-Xoo interaction, we investigated the expression pro le of circRNAs and identi ed differentially expressed circRNAs (DEcircRNAs) from rice ag leaves infected by PXO99 A at different hpi. The expression pro le revealed a total of 276 circRNAs were detected as signi cantly DEcircRNAs in response to PXO99 A infection ( Fig. 2A). Among them, 91, 129, 129 and 138 DEcircRNAs were respectively detected at 2, 6, 12 and 24 hpi. The 112 circRNAs were shared by at least two postinfection time points and 30 circRNAs were overlapped among different time points ( Fig. 2A). As to the up-and down-regulated DEcircRNAs, the down-regulated circRNAs were obviously more than the upregulated genes at 2 hpi, and less than the up-regulated genes at other time points (Fig. 2B). Like all circRNAs in rice, exonic DEcircRNAs were predominant circRNAs (65.22%), followed by intergenic (31.88%) and antisense circRNAs (2.90%) (Fig. 2C), moreover, most DEcircRNAs were also shorter than Functional categorization of circRNA parental genes Since circRNA biogenesis has close relationship with their functions in plants, GO and KEGG pathway analysis were performed for the host genes of the circRNAs collected from all the time points. The host genes are involved in three main functional categories, including cellular component, molecular function, and biological process (Fig. 3A). The major enrichments for the cellular component terms were 'cell', 'organelle' and 'membrane' while 'catalytic activity' and 'binding' were the most enrichments for molecular function terms with abundant circRNAs. The GO terms 'cellular process', 'metabolic process' and 'response to stimulus' were signi cantly enriched and accounted for the most of biological process, indicating the putative involvement of circRNAs in regulating rice response to pathogen infection. KEGG pathway analysis showed that the host genes were signi cantly enriched in carbon xation in photosynthetic organisms and glyoxylate and dicarboxylate metabolism (p.adjust < 0.05), followed by photosynthesis proteins and photosynthesis ( Figure S4, Table S2).
To uncover whether the DEcircRNAs-producing genes possess any bacterial blight resistance function, the parental genes of all the DEcircRNAs were also predicted and annotated. The photorespiration and chloroplast-related terms (like chlorophyll biosynthetic process, chloroplast membrane, glycolate oxidase activity, thylakoid) were signi cantly enriched by DEcircRNA-producing genes, indicating their vital role in rice-pathogen interaction (Fig. 3B, Table S2). In addition, we found involvement of DEcircRNA-producing genes in regulation of 'glyoxylate and dicarboxylate metabolism', 'arginine biosynthesis' 'autophagy', and 'peroxisome' (Fig. 3C, Table S2), inferring the relevance between circRNAs and rice response to pathogen infection.
circRNA-miRNA-mRNA regulating network Given the functional regulatory evidence of some circRNAs with binding capabilities to miRNAs , we constructed the circRNA-miRNA-mRNA network in this study to reveal the potential role of these DEcircRNAs in post-transcriptional regulation. As a result, 31 circRNAs out of the 276 DEcircRNAs harbored one or more miRNA binding sites for 121 miRNAs (Table S3). Furthermore, target prediction showed 4192 mRNAs to be the potential targets of 36 miRNAs, including several conserved miRNAs like miR398, miR444 and miR812 (Table S4). Functional enrichment analysis was subjected to uncover any resistance functions possessed by the targeted mRNAs. The GO enrichment analysis suggested that the targeted mRNAs were involved in various stress-responsive, such as pattern speci cation process and gibberellin metabolic process (Table S5). The KEGG pathways revealed signi cant enrichment of diterpenoid biosynthesis (p.adjust < 0.05), which was followed by starch and sucrose metabolism and transcription factors (TFs) (Fig. 4A, Table S5). To predict the functions of circRNAs more accurately, circRNA-miRNA-mRNA regulating network corresponding to the top 3 enriched pathways was further constructed using Cytoscape (v3.8.0) (Fig. 4B). In this study, the predicted miRNA target mRNAs included various TFs, like MYB, ABF, WRKY and TGA, which were predicted to participate in plant-pathogen interaction and plant hormone signal transduction (Table S5).
Experimental validation of the identi ed circRNAs To con rm these circRNAs detected in this study, PCR and Sanger sequencing were performed to validate the back-splicing sites of 7 randomly selected circRNAs. 7 pairs of convergent primers were used to detect linear transcripts from gDNA and cDNA as positive controls (Fig. 5). In contrast, all 7 pairs of circRNAs can only be ampli ed in cDNA samples, but not in gDNA samples (Fig. 5), suggesting the presence of the head-tail splicing junctions. PCR products with expected length of divergent primers were further veri ed to span the back-splicing junction by Sanger sequencing and sequence mapping (Fig. 5). Moreover, the expression of four DEcircRNAs were detected by RT-qPCR using cDNA of total RNAs as templates.
Obviously, in comparison with 0 h, the expression of all the four circRNAs were signi cantly up-regulated at different hpi (Fig. 6), which is basically consistent with the deep sequencing results, indicating that the circRNAs identi ed in our study are reliable (Fig. 6). The above results revealed that those identi ed circRNAs may participate in rice innate immunity or the regulation of rice bacterial leaf streak disease development.

Discussion
As an emerging potential class of regulatory RNAs, circRNAs have attracted the attention of many researchers. Widespread and various circRNAs have been detected in both animals and plants with the progress in the development of sequencing and bioinformatics technology (Kristensen et al. 2019;Zhao et al. 2019). Recent studies also highlighted the tight correlation between circRNAs and plant development and stress-induced responses Litholdo and da Fonseca 2018). In case of rice, current evidence has pointed out the potential biological roles of circRNAs in response to fertility transition , phosphate starvation (Ye et al. 2015) and fungus infection (Fan et al. 2020). However, a little is known regarding the regulation roles of circRNAs in rice defense responses against pathogens.
To obtain comprehensive and reliable circRNA predictions, two state-of-the-art tools, CIRCexplorer2 ) and CIRI2 , were adopted for circRNAs characterization in this study. Obviously, there are quite different in terms of the number of detected circRNAs and little overlap between the two prediction algorithms when they were adopted for rice circRNA libraries. Such differences are largely due to different algorithms possessing different advantages and shortcomings regarding the prediction accuracy and sensitivity of detecting circRNAs , no relation to animal or plant species. For example, CIRCexplorer identi ed the most grape circRNAs, followed by nd_circ and CIRI (Gao et al. 2019). While in Arabidopsis thaliana, nd_circ identi ed the most circRNAs among the three methods (Zhang et al. 2020b), however, CIRI identi ed the most circRNAs when they were applied to human circRNA libraries (Zeng et al. 2017). Thus, it is necessary to combine several algorithms to achieve reliable and comprehensive predictions in both animal and plant species Gao et al. 2019;Hansen et al. 2016). Several rice circRNA back-splicing sites and expression pattern were experimentally validated by divergent PCR combined with Sanger sequencing and RT-qPCR, respectively. All these data further con rmed the reliability of prediction results in this study.
With respect to the characteristics of circRNAs identi ed in this study, 2303 (about 65.5 %) were homologous to the rice circRNA sequences in PlantcircBase database and 1214 were novel in rice. Consistent with previous studies, all these circRNAs derived from diverse genomic regions, and the exon circRNAs accounted for the largest proportion (Fan et al. 2020;Wang et al. 2019;Ye et al. 2015;Gao et al. 2019;Tong et al. 2018), which can be explained as the exon skipping events contributing to the most circRNA formation in these plants (Conn et al. 2017). However, most circRNAs were intergenic in kiwifruit ) and potato . In soybean, the highest proportion of intronic circRNAs occurred in root, while exonic in the stem and leaf (W. . Integrating these ndings showed that plant circRNAs biogenesis are tissue-and species-speci c. A recent study of cucumber also found that the circRNA numbers are signi cantly positive in correlation with chromosome length (Zhu et al. 2019), which is also con rmed by our results and previous study in rice , the longest chromosome in rice genome (chromosome 1) producing more circRNAs than others.
Moreover, diverse non-GT/AG splicing signals, except the canonical GT/AG, were also used by the most circRNAs in our study, which is consistent with the results in rice (Ye et al. 2015) and cucumber (Zhu et al. 2019), and inconsistent with those in grape (Gao et al. 2019) and cotton (T. , inferring a circular preference of circRNAs processing in most plants and highlighting the complexity of backsplicing in rice. Plant stress responses are very complex processes and involve many speci c-inducible genes and signaling transduction pathways. We observed that PXO99 A infection could speci cally induce the formation of some circRNAs in rice, of which the host genes were signi cantly related to various primary metabolic pathways, mainly involving photorespiration, chloroplast, and peroxisome-associated pathways (Table S2). Similar results were also indicated in the GO enrichment analysis for protein-coding genes in kiwifruit induced by bacterial canker pathogen invasion ). More importantly, evidences from the literature have con rmed the involvement of these metabolic pathways in plant defense responses (Kangasjärvi et al. 2012;Rojas et al. 2014;Ho us et al. 2017;Mammarella et al. 2015). For instance, glycolate oxidase (GOX), the rst enzymes in the photorespiratory pathway, could be induced in response to fungal pathogen in Brassica napus, Arabidopsis (Bohman et al. 2002) and Hordeum vulgare (Schäfer et al. 2004). While gox mutants of Arabidopsis and Nicotiana benthamiana showed susceptibility to non-host pathogens (Rojas et al. 2012). Chloroplast Cu/Zn superoxide dismutase (CSD2) participates in controlling ROS levels in infected tissues (Mateo et al. 2004).
Peroxidases PRX33 and PRX34 play important roles in Arabidopsis basal resistance via regulating apoplastic ROS production (Daudi et al. 2012). As outlined above, the positive correlations between circRNAs expression and that of their cognate linear genes have been observed in many plant species, such as Arabidopsis (Ye et al. 2015), rice (Ye et al. 2015), kiwifruit , grape (Gao et al. 2019), etc. These ndings con rmed our hypothesis that the DEcircRNAs have valuable functions in the synthesis of metabolic processes through regulation of their parent genes expression during rice-PXO99 A interaction.
Besides cis-regulating parent genes, circRNAs have been proposed to sponge to speci c miRNAs and to further prohibit them from regulating their target mRNAs . Here, we found 31 DEcircRNAs containing one or more miRNA binding sites in circRNA-miRNA-mRNA network. Unlike the various functions of circRNAs producing genes, mRNAs predicted in this network were only signi cantly involved in diterpenoid biosynthesis, which has been con rmed to be relevant to rice plant defense against bacterial pathogen Xoo recently (Lu et al. 2018). One possible explanation was that miRNA sponges might not be the main mechanism for plant circRNAs (Chu et al. 2018). In addition, several conserved miRNAs were predicted to be the targets of certain rice DEcircRNAs, including miR156, miR398, miR414, miR444, and miR812, which are well known to regulate a wide range of developmental and physiological processes (Baldrich and San Segundo 2016;Kansal et al. 2015;Kruszka et al. 2012;Li et al. 2014;Wang et al. 2016). Particularly, rice plants overexpression miR398b showed enhanced resistance to M. oryzae by in uencing H 2 O 2 accumulation at the infection site and defense gene expression (Li et al. 2014). These results further implied that these circRNAs might make a regulatory contribution to rice immunity, whether control rice immunity via acting as miRNA sponges need experimental evidence in the future.

Conclusions
In this study, all these ndings collectively con rmed rice circRNAs may also regulate rice defense against pathogens at multiple levels. Functional characterization of rice circRNAs and discovery of their interaction machinery between host genes or miRNAs will aid to develop appropriate strategies for rice protection.

Plant materials and pathogen inoculation
The rice cultivar Nipponbare (Oryza sativa L. ssp. japonica) seeds were soaked in water and placed in an incubator (28 °C) for 2 days. The germinated seeds were sown in pots containing nursery soil, placed in an illuminated growth chamber with a 14 h light/10 h dark cycle at 70% RH and 28-32 °C. One-month-old rice seedlings were inoculated with wild Xoo strain PXO99 A (laboratory collection) at 10 9 cfu ml −1 by the leaf-clipping method (Kauffman 1973). Brie y, the tips of rice leaves are cut with scissors previously dipped in bacterial suspension. Approximately 2 cm-length leaf fragments next to the inoculated sites were harvested from two stages: before inoculation (0 h, used as a control) and after inoculation (2, 6, 12 and 24 h), respectively. At each time point, samples were collected from at least eighteen different plants and immediately frozen in liquid nitrogen.
CircRNA-seq library preparation and sequencing A total of 5 samples from 5 time points were harvested for library construction. The RNA-Seq library construction and sequencing were performed by the Novogene Corporation (Beijing). RNA purity and concentration were assayed by using the NanoPhotometer® spectrophotometer (IMPLEN, CA, USA) and Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA). The integrity of total RNAs was further evaluated by gel electrophoresis and an Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA).
A total amount of 5 μg RNA per sample was used for circRNA-seq library preparation. Ribosomal RNAs were depleted using an Epicentre Ribo-ZeroTM Kit (Epicentre Technologies Corp., Chicago, IL), whereas linear RNAs were digested with RNase R (Epicentre, USA). The remaining RNAs were broken into small fragments under elevated temperature and converted into double stranded DNAs. The cDNA fragments (250 ~ 300 bp in length) were puri ed and used for PCR ampli cation. Finally, PCR products were puri ed and libraries were quality-controlled using the BioAnalyzer 2100 system. The resulting libraries were sequenced using an Illumina HiSeq 2000 system and 150 bp paired-end reads were generated.

Identi cation of circRNAs
Low-quality reads (Q ≤ 20) and adapters were removed from the sequencing data by Trim galore (https://github.com/FelixKrueger/TrimGalore) before circRNA identi cation. The remaining clean reads were mapped to the rice reference genome (IRGSP-1.0), generating a sequence alignment map le by bwa (v0.7.17, mem-T 19) with default parameters (Li and Durbin 2009). The output of BWA was used to identify the circRNAs by CIRIquant (V1.1) (Zhang et al. 2020a) . CIRI2 and CIRCexplorer2 are the two tools that were used in the pipeline Zhang et al. 2016). One or more base differences might be present between the results of different prediction software programs. Then, repeated circRNAs in the prediction results were removed by CD-HIT-EST following a method proposed in a previous report (Fu et al. 2012;Zhang et al. 2020b): (1) a length difference between the two sequences less than 10 bp and (2) an alignment sequence exceeding 99.7% of the shorter sequence. The other CD-HIT-EST parameters were the default parameters.

Differential expression and annotation analysis
All the differentially expressed circRNAs (DEcircRNAs) were examined in the comparison after PXO99 A infection (2, 6, 12 and 24 h) versus 0 h using the CIRIquant (V1.1) with default parameters. In addition, the circRNAs that DE_score value equal to zero were removed by host python script. The functions of some DEcircRNAs may depend on their host genes; thus, the function of DEcircRNAs host genes were annotated by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases. Gene enrichment analysis of these host genes were performed using clusterPro ler R package (Yu et al. 2012). All detected DEcircRNAs were compared against the PlantcircBase database with BLASTN. We did BLAST (with a threshold of E-value < 1E-5) for nding the conserved circRNAs using our identi ed DEcircRNAs as the query.
Prediction for circRNA-miRNA-mRNA relationships The functions of some circRNAs may be independent of their host genes and may act as miRNA decoys to achieve their effects (Hansen et al. 2013). To construct the circRNA-miRNA-mRNA network, miRNA sequences were obtained from miRBase (http://mirbase.org/) and the plant microRNA database (PMRD: http://bioinformatics.cau.edu.cn/PMRD/) (Kozomara et al. 2019;Zhang et al. 2010), while mRNA sequences were downloaded from the rice annotation project database (https://rapdb.dna.affrc.go.jp/index.html). The DEcircRNAs sequences were extracted with an in-house python script. Then, the miRNAs identical alignment regions in circRNAs and mRNAs were predicted by GSTAr.pl (https://github.com/MikeAxtell/GSTAr), and the minimum free energy (MFE) of miRNA-circRNA or miRNA-mRNA duplexes was calculated with the RNAhybrid program. Then, the miRNA-targeted mRNA and miRNA-decoyed circRNA were predicted following a method proposed in a previous report (Tang et al. 2018;Ivashuta et al. 2011). The general criteria used to de ne a miRNA decoy were as follows: no more than six mismatched or inserted bases present between the 9 th to 20 th nucleotides of the miRNA 5' end, perfect matching of the 2 nd to 8 th bases of the miRNA 5' end sequence, and no more than four mismatches or indels in other regions. The circRNA-miRNA-mRNA regulating network was generated by Cytoscape (v3.8.0) (Shannon et al. 2003).

Validation of circRNAs and quantitative real-time PCR (RT-qPCR)
Genomic DNA was extracted from fresh rice leaves using Dzup (plant) Genomic DNA Isolation Reagent (Sangon Biotech, Shanghai). In accordance with the manufacturer's protocol, total RNA of all samples was isolated using EZ-10 Total RNA Mini-Preps Kit (Sangon Biotech, Shanghai) and reverse-transcribed into cDNA with a combination of random and oligo (dT) primers using the Hifair® III 1st Strand cDNA Synthesis Kit (gDNA digester plus) (YEASEN) in accordance with manufacturer's protocol.
To validate the identi ed circRNAs in this study, a pair of divergent primers and convergent primers (Table  S6) were designed for each circRNA using Primer5 software, and both cDNA and gDNA (negative control) were used as template for polymerase chain reactions (PCRs). For each PCR ampli cation, 100 ng of cDNA or genomic DNA was used with TaKaRa Ex Taq ® DNA Polymerase. The PCR results were con rmed by agarose gel electrophoresis and Sanger sequencing.
The expression level of circRNA candidates in rice leaves were experimentally detected by RT-qPCR with Magic SYBR Green qPCR Mix (Magic Biotech, Hangzhou, China) and the ABI 7500 quantitative PCR system (Applied Biosystems, Foster City, CA). Data were calculated using the 2 −∆∆CT method and were normalized to the expression of the reference gene OsActin. Each assay was performed three times in triplicate and the results are displayed as the means ± SD of three biological replicates.