Comparison of tolerant and susceptible cultivars revealed the roles of circular RNAs in rice responding to salt stress

As a newly characterized class of noncoding RNAs, circular RNAs (circRNAs) have been identified in many plant species, and play important roles in plant stress responses. However, little is known about how salt stress mediates the expression of circRNAs in rice. In this study, we identified circRNAs from root tissues of salt-susceptible recipient cultivar 93–11 and salt-tolerant introgression line 9L136. A total of 190 circRNAs were identified. Among them, 93 circRNAs were differentially expressed under salt stress in 93–11 (36 up- and 57 down-regulated) and 95 in 9L136 (46 up- and 49 down-regulated). Salt stress significantly decreased the average expression level of circRNAs in 93–11, but circRNA expression levels were slightly increased in 9L136, suggesting that circRNAs had different response patterns in these two cultivars. Function annotation and enrichment analysis indicated that, through cis-regulation and circRNA-miRNA-mRNA network regulation, those induced circRNAs were commonly involved in transcription, signal transduction, ion transportation, and secondary metabolism. Compared to 93–11, salt-induced circRNAs in line 9L136 targeted more stress response genes participating in transcription regulation, ion transportation, and signal transduction, which may contribute to the salt tolerance of 9L136. Summarily, this study revealed the common response of rice circRNAs to salt stress, and the possible circRNA-related salt tolerance mechanisms of 9L136.


Introduction
During the past decade, a novel type of noncoding RNAs, termed as circular RNAs (circRNAs), have been discovered and widely studied . Unlike formerly known noncoding RNAs, circRNAs are covalently closed RNA circles formed by head-to-tail splicing of transcripts that are devoid of the 5'-cap and 3'-tail (Zhu et al. 2019a). With the development of high-throughput sequencing technologies, numerous circRNAs have been discovered in both animals and plants in recent years . In general, circRNAs can be classified into three main categories based on their genomic location, exonic circRNAs, intronic circRNAs, and intergenic circRNAs . Some circRNAs are conserved, more abundant than their linear counterparts, and exhibit cell, and tissue-specific expression patterns in different plant species, suggesting the potential physiological functions of these circRNAs (You et al. 2015;Zhao et al. 2017).
Although circRNAs are generally expressed at low levels, they may be involved in many physiological and molecular processes. Several circular RNAs were recently demonstrated to possess microRNA (miRNA) binding sites, which may function as sponges to sequester miRNAs and prevent their interactions with target mRNAs (Yin et al. 2017;Zhou et al. 2018). However, the miRNA sponge functions of most circRNAs have not been proved experimentally. In addition, many studies have proposed that some circRNAs could be translated into proteins in a cap-independent manner. For example, Zhang et al. (2018) showed that a circRNA containing an open reading frame (ORF) driven by the internal ribosome entry site (IRES) can translate a functional amino acid protein termed as SHPRH-146aa. Moreover, cir-cRNAs can interact with functional proteins to exert their functions in multiple processes. Li et al. (2017) reported that in humans, double-stranded RNA-binding domains containing immune factors NF90 and NF110 were components of circRNPs in the cytoplasm and could enhance circRNA production in the nucleus through stabilizing intronic RNA pairs. However, the general functions of most circRNAs remain far from clear. In plants, circRNAs have been identified and studied in increasing number of species, including Arabidopsis (Liu et al. 2017a, b), rice (Ye et al. 2016), wheat , tomato (Tan et al. 2017;Yin et al. 2018), kiwifruit (Wang et al. 2017b), soybean , potato , tea (Tong et al. 2018) and cucumber Zhu et al. 2019a). Similar to animal circRNAs, plant circRNAs are abundant, conserved, and expressed at low concentrations and in a tissue-specific manner. However, plant circRNAs also have distinct features, such as a significantly lower proportion of both repetitive elements and reverse complementary sequences in introns flanking exonic circRNAs (Meng et al. 2016;Chu et al. 2020). Moreover, plant circRNAs might also be involved in a broad range of biological pathways . For example, some circRNAs exhibited different expression patterns during fruit ripening (Yin et al. 2017) and under biotic and abiotic stresses (e.g. dehydration, salt, and heat stresses) (Wang et al. , 2017bPan et al. 2018;Zhu et al. 2019a). In addition, circRNAs can affect their parental genes in cis-or trans-actions. In rice, circRNAs has been found to negatively regulate the expression of their parental genes (Lu et al. 2015). In Arabidopsis, Conn et al. (2017) demonstrated that certain circRNAs have the ability to skew splicing preference and increase abundance of their cognate exon-skipped alternative splicing (AS) mRNA variants, thereby revealing that specific circRNAs have the ability to regulate the splicing of their cognate AS mRNAs. Several studies in animals have shown that circRNAs could serve as "miRNA sponges" to regulate mRNA expression, but "miRNA sponges" function of circRNAs in plants have not been experimentally characterized. These findings suggest that biological roles of circRNAs in plants are diverse. However, to date, few experiments have analyzed the functional and formation mechanisms of circRNAs in plants.
Globally, rice (Oryza sativa L.) is one of the most important crops and is quite sensitive to salinity (Gong et al. 2006). In rice, Ye et al. (2015) conducted a genome-wide identification of circRNA for the first time, and revealed the potential regulatory functions of 27 differentially expressed exonic circRNAs in response to Pi-starvation stress. Furthermore, Lu et al. (2015) identified 2354 circRNAs through deep sequencing, and the computational analysis of transcriptome data and transgenic analysis suggested that circRNA and its linear form could reduce the expression of its parental gene. Ye et al. (2016) identified nearly 3000 circRNAs in rice and found that non-GT/AG splicing signals were common in rice circRNAs. However, to our knowledge, there is still limited data on the role of circRNAs in the salt stress resistance of rice . In this study, we aimed to provide insight into the possible roles of circRNAs in response to salt stress, which could provide a basis for further functional studies of rice circRNAs.

CircRNA identification
Rice RNA-seq dataset used in this study was accessible in NCBI (National Center for Biotechnology Information) under the SRA (Sequence Read Archive) database (Bio-Project PRJNA395311). The dataset includes the control (0 h, CK) and salt treatment (1 h after treatment, 125 mM NaCl) samples from roots of the recurrent parent 93-11 and introgression line 9L136 seedlings that were cultured in full-strength Yoshida's solution to two and a half leafstage. Different samples were sequenced in three biological replicates (Wang et al. 2017c). Low quality reads, including (1) unknown (N) bases that were greater than 5%, (2) contained adaptor sequences, and (3) contained more than 50% bases with Q ≤ 20, were removed. The remaining highquality reads were mapped to the previously reported Oryza sativa L. reference genome RGAPv7 (MSU Rice Genome Annotation Project Release 7) using BWA (v0.7.15, mem -T 19) (Zhu et al. 2019a). CIRI2 (v2.05) is considered to be the most reliable software, with remarkably balanced sensitivity and reliability (Gao et al. 2017), thus the output SAM files of BWA were inspected using CIRI2 to identify circRNAs. CircRNA expression levels were normalized to be SRPBM (SRPBM = number of circular reads/number of mapped reads (units in billion)/read length). Those known rice circR-NAs were download from the database PlantCircBase (Chu et al. 2017). Junction reads were manually extracted and blastn was performed against known rice circRNAs (-evalue 0.05).

CircRNA annotation and GO enrichment
Databases including Nr (NCBI non-redundant protein sequences), GO (Gene Ontology), KEGG (Kyoto Encyclopedia of Genes and Genomes), KOG/COG (Clusters of Orthologous Groups of roteins), Swiss-Prot (A manually annotated and reviewed protein sequence database, https:// www. unipr ot. org/), and Pfam (Protein family) were used to perform gene function annotation . The identification of differentially expressed circRNAs was performed using the R package EBSeq (Zhu et al. 2019a). The cut-off criteria were set as an SRPBM value larger than 0.5 and a fold change larger than 2. GO annotations of the parent genes of differentially expressed circRNAs were collected and used to perform GO enrichment using topGO with default parameters. Venn diagram was draw using online tools (http:// bioin forma tics. psb. ugent. be/ webto ols/ Venn/). R package "circlize" (Gu et al. 2014) and "vioplot" (Adler and Kelly 2021) were used to draw the circos and violin diagrams.

Prediction of circRNA-miRNA-mRNA relationships
To construct the circRNA-miRNA-mRNA network, we collected the mature sequences of rice miRNAs from the miRBase. Then the miRNAs sponged by circRNAs were predicted by TargetFinder using default parameters (Zhu et al. 2019a). Meanwhile, rice coding genes (mRNAs) targeted by collected miRNAs were predicted by searching the CDS sequences of the RGAPv7 reference genome with an online tool (http:// plant grn. noble. org/ psRNA Target/) (Fang et al. 2020). The circRNA-miRNA-mRNA regulating network was generated by Cytoscape (v3.6.0) . The eggNOG class annotations of these miR-NAs targeted mRNAs were collected and used to interpret the possible regulating roles of circRNAs.

Validation of circRNAs
To further confirm the circRNAs identified from rice, total RNA was extracted from root sample of 93-11 using the Plant RNA Kit (Omega, London, UK). The cDNA was then synthesized from 500 ng of total RNA using the RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific, USA) according to the manufacturer's protocol. The divergent PCR primers were designed using the "outfacing" strategy to exclude linear mRNA from amplification . Sanger sequencing was further used to confirm the presence of back-splicing junctions.

Identification of rice circRNAs
To explore candidate circRNAs related to salt tolerance, the RNA-seq dataset obtained from the root tissues of salt-susceptible recipient cultivar 93-11 and salt-tolerant introgression line 9L136 (seedlings after 0 (CK) and 1 h of NaCl treatment (Na)) was used to perform the identification, expression profiling, and function annotation of rice circRNAs. In total, 864 million reads were obtained from 12 libraries, and about 94.98% of the reads could be mapped to the reference genome by BWA (Supplementary  Table S1). Then, the mapping SAM files were scanned by the CIRI2 tool with the default parameters, resulting in the identification of 14,261 junction reads (Fig. 1A).  Table S2). Among these 190 circRNAs, 83 and 66 were detected from CK and Na in 93-11, respectively; meanwhile, 88 and 82 were found from CK and Na in 9L136, respectively. Moreover, 22 circRNAs were detected in all samples (Fig. 1B). Further validation revealed that the majority of junction reads (12,801, 90%) overlapped with the reported circRNAs, hence suggesting the reliability of detected circRNAs in this study.

Expression patterns of rice circRNAs under salt stress
To study the effect of salt stress on rice circRNAs expression, we compared the average expression levels of all detected circular RNAs. A Wilcoxon rank-sum test indicated that, generally, the average expression levels of circRNAs were significantly decreased by salt stress in recipient 93-11 (p-value < 0.05. Figure 1H), whereas the expression levels of circRNAs were slightly increased by salt stress in introgression line 9L136. In order to reveal rice circRNAs that may have biological functions in response to salt stress, we Data are represented as the log 2 (SRPBM + 1). The gray dot represents the median. *p < 0.05, Wilcoxon rank-sum test. I Number of differentially expressed circRNAs. 93-11 and 9L136, two rice cultivars; vs., versus compared the expression profiles of differentially expressed circRNAs between control and salt stress samples. Among the 190 circRNAs, 96 circRNAs were found to be significantly differentially expressed under salt stress, including 93 in recipient 93-11 (36 up-and 57 down-regulated, Supplementary Table S3) and 95 in introgression line 9L136 (46 up-and 49 down-regulated, Supplementary table 4) (Fig. 1I).

Prediction of the circRNA-miRNA-mRNA regulating network
In addition to regulating parent genes, circRNAs were found to competitively and specifically bind miRNAs to prohibit miRNAs from the epigenetic regulating of target mRNAs (Liu et al. 2017a, b). In the present study, 27 miRNAs were predicted to be sponged by 19 differentially expressed cir-cRNAs in 9L136, and 37 miRNAs were predicted to be sponged by 26 differentially expressed circRNAs in 93-11, suggesting that these circRNAs may function as miRNA sponges in rice (Supplementary Tables S5-8). Furthermore, 426 mRNAs were predicted to be the targets of 27 miRNAs in 9L136, and 327 mRNA were predicted to be the targets of 37 miRNAs in 93-11 (Supplementary Tables S5-8). Egg-NOG class annotation analysis suggested that 19 circRNAs of 9L136 play roles in the regulation of physiological and biochemical processes, including transcription, signal transduction, translation, secondary metabolism, and inorganic ion transport, etc. (Fig. 3A, Supplementary Figure S1, Supplementary Table S6). Similar results were found in 93-11. CircRNAs also play roles through involving EggNOG The pink, magenta, cyan, and coral rectangle nodes represent categories of transcription, transport, enzyme, and component, respectively. Detailed information, including binding sites, relationship among cir-cRNAs, miRNAs, and mRNAs, and annotation of mRNAs were listed in Table S4, S5 classes such as transcription, translation, signal transduction, and secondary metabolism (Fig. 3B, Supplemental Figure  S1, Supplementary Table S8).

Experiment confirmation of salt-induced circRNA in rice
To further confirm the existence of rice circRNAs detected by this study, divergent primers were designed to perform the divergent polymerase chain reaction (PCR). Ten selected circRNAs were confirmed to be processed from back-splicing, demonstrating the reliability of our cicRNA identification based on high-throughput sequencing (Fig. 4A, 4B. Supplementary Table S9). Besides, the expression levels and regulation patterns of differentially expressed circR-NAs and target mRNAs were further confirmed by qRT-PCR. The up-and down-regulation patterns of circRNAs and mRNAs determined by RNA-seq were consistent with qRT-PCR (Fig. 4C, 4D. R 2 = 0.90 **, Supplementary Figure  S2), suggesting that: (1) expression patterns of circRNAs determined by RNA-seq were reliable, and (2) the possibility of circRNA-miRNA-mRNA regulating networks were further confirmed.

CircRNA identification
Rice circRNAs had been identified in multiple tissues, including leaves, anthers, pistils, seeds, shoots, and roots, implying the roles of circRNAs in different rice growth and development stages (Lu et al. 2015;Ye et al. 2015;Chu et al. 2017Ye et al. 2017. However, it is still much unclear how circRNAs participate in the response of rice to salt stress. Previously, Wang et al. (2017c) developed salt-tolerant introgression line 9L136 using wild Oryza rufipogon accession as the donor, and Oryza sativa indica cultivar 93-11 as the recipient. RNA-seq profiling revealed more salt stress-induced differentially expressed genes (1391) in introgression line 9L136 than recipient 93-11. In this study, this RNA-seq profiling dataset was used to explore candidate salt tolerance-related circRNAs. In total, 190 circR-NAs supported by at least two unique back-spliced reads were detected, including 118 in 93-11 and 130 in 9L136 (Fig. 1B).

circRNAs related to salt stress
In plants, many circRNAs have been reported to exhibit stress-inducible expression patterns . For example, Zhu et al. (2019a, b, c) identified 2787 cir-cRNAs in cucumber, with 1934 in root and 44 in leaf being differentially regulated under salt stress. Annotation and enrichment analysis of both parental genes and target mRNAs of salt-induced circRNAs suggested that circRNAs may paly roles in salt stress response by regulating abiotic stress associated genes. In Arabidopsis, Pan et al. (2018) uncovered 1583 heat-specific circRNAs, suggesting these circRNAs may participate in heat stress response. Plant salt stress tolerance is a very complex regulating progress involving in signal transduction, substance and energy metabolism (Zhu et al. 2019b, c). In this study, Wilcoxon rank-sum test indicated that salt stress significantly decreased the average expression level of circRNAs in salt-susceptible cultivar 93-11 (p-value < 0.05. Figure 1H), but slightly increased it in salt tolerant introgression line 9L136.
In order to reveal rice circRNAs that may have biological functions in response to salt stress, we compared the expression profiles of differentially expressed circRNAs between control and salt stress samples. Among the 190 circRNAs, 96 circRNAs were significantly differentially induced by salt stress, including 93 in recipient parent Supplementary (Fig. 1I). Since circRNAs have been proposed to be positive or negative regulators on their parent coding genes (Lu et al. 2015;Cheng et al. 2018), we predicted and annotated the parent genes of these differentially expressed circRNAs. GO enrichment analysis implied that the parent genes of differentially expressed circRNAs in 9L136 were involved in more biological and biochemical processes than in salt-susceptible cultivar 93-11 ( Fig. 2A

CircRNA-miRNA-mRNA regulating network prediction
Besides regulating parent genes, previous studies suggested that circRNAs could bind specific miRNAs to repress the regulating ability of miRNAs (Liu et al. 2017a, b). In this study, 27 miRNAs and 37 miRNAs were predicted to be competitively bound by 19 and 26 differentially expressed circRNAs in 9L136 and 93-11, respectively (Supplementary  Tables S5-8). Among those miRNAs, some were reported to be closely related to plant stress responses in rice. For example, osa-mir2925 was reported to be up-regulated by salt stress in rice, and osa-mir2925 potentially targets stress-response-associated genes, such as LOC_Os05g24780 (calcium-binding protein CML21), LOC_Os09g28200 (Heat stress transcription factor B-4c) in 9L136, and LOC_ Os01g43320 (GABA transporter) in 93-11 (Fig. 3) (Wang et al. 2017c). Another osa-miR2925 was competitively Fig. 4 Validation of rice circRNA through PCR and Sanger sequencing. A M, marker; 1 to 10, ten circRNA validation using Sanger sequencing. (Detailed information, including the corresponding ID of 1 to 10 circRNAs and their divergent primers, can be found in Supplementary Table S8). B Using Chr10:20,345,844|20,346,873 as an example. Upper left panel (Model), a model used to show the divergent primers for circRNAs backsplicing sites amplification. Upper right panel, an example showing that amplification of divergent primers can be amplified from cDNA but not from genomic DNA. Two middle panels, respectively, black rectangle (DNA sequence) representing downstream and upstream sequences in the genome, and red and green rectangles (circRNA sequence) representing back-splicing circRNA sequence. Lower panel, a Sanger sequencing example of rice circRNA Chr10:20,345,844|20,346,873. Validation the expression levels and regulation patterns of circRNAs and target mRNAs in control and salt-treated 9L136 C and 93-11 D rice seedlings by comparing RNA-seq and qRT-PCR data. Right histograms with names initialed by LOC_Os represent mRNA, and left histograms with names initialed by Chr represent circRNA. Histogram height represents log2(fold change) of expression levels of circRNAs and target mRNAs between control and salt-treated samples bound by circular RNA Chr10:18,818,949|18,958,248 in 9L136 and by Chr6:2,046,190|2,059,Chr10:18,818,949|18,958,248 and Chr6:2,046,190|2,059,120 showed distinct regulation patterns in 9L136 and 93-11 under salt stress. In salt-tolerant line 9L136, Chr10:18,818,949|18,958,248 was down-regulated, which may decrease its binding effect to osa-mir2925. Meanwhile, in 93-11, Chr6:2,046,190|2,059,120 was up-regulated, which may promote the sponge effect of Chr6:2,046,190|2,059,120 to osa-mir2925. It seems that in this regulation network, the distinct regulation patterns of circular RNAs in salt-resistant line 9L136 and salt-susceptible cultivar 93-11 result in the opposite regulating effect on miRNA, which further leads to the differential expression of stress-response-associated mRNAs, and may contribute to the salt-resistant phenotype of 9L136 and salt-susceptible phenotype of 93-11.
Additionally, it should be noticed that several miRNAs competitively bound by circRNAs were markedly different between the two rice cultivars, thus indicating the different regulating networks of circRNAs in the two rice cultivars in response to salt stress (Fig. 3. Supplementary Tables S5-8).
For example, the salt-tolerant introgression line 9L136 showed several unique regulation networks. In 9L136, osa-miR2102-5p and osa-mir5809 were predicted to be competitively bound by circRNAs Chr10:18,818,949|18,958,248 and Chr7:14,201,482|14,303,571, respectively. The osa-miR2102-5p has been widely reported to be involved in salt and drought stresses in Gramineae plants, such as barley (Zare et al. 2019), maize (Wang et al. 2014a, b), and Spartina alterniflora (Qin et al. 2015). In rice, osa-mir5809 was reported to be involved in salt (Huang et al. 2019) and heat stresses (Mangrauthia et al. 2017), as well as leaf senescence (Xu et al. 2014). Both circRNA Chr10:18,818,949|18,958,248 and Chr7:14,201,482|14,303,571 were down-regulated by salt stress in 9L136, which could decrease the competitive binding of circRNAs to miRNA osa-miR2102-5p, and further enhance the suppression of osa-miR2102-5p to target mRNAs. These findings implied the involvement of circR-NAs in the response of rice to salt stress. CircRNAs may also have various biological functions during the growth and stress responses of rice. However, the function of circRNAs in rice needs further experimental validation. Further analysis demonstrated that 426 mRNAs were predicted to be the targets of 27 miRNAs in 9L136, and 327 mRNA were predicted to be the targets of 37 miRNAs in . EggNOG class annotation analysis revealed that these mRNAs participated in the regulation of physiological and biochemical processes, including transcription, signal transduction, and secondary metabolism (Fig. 3, Supplementary Figure S1, Supplementary Tables S5-8). Among those miRNA targets, more mRNAs were differentially expressed in 9L136 (108) than in 93-11 (83), suggesting a more complex regulation network in the salt stress tolerance introgression line 9L136.

Expression patterns of circRNAs and their corresponding parental genes in different cultivars
Previous studies in several species have suggested that most circRNAs regulate the expression level of their corresponding parental genes. In rice, Lu et al. (2015) proposed that circRNA and its linear form may act as a negative regulator of its parental gene. However, expression profiles of circR-NAs in tea and Arabidopsis showed a positive correlation between circRNAs and their parental genes (Tong et al. 2018;Cheng et al. 2018). In cucumber, both opposite and positive trends between the expression levels of circRNAs and parental mRNAs have been reported (Zhu et al. 2019a).
In salt-sensitive cultivar 93-11, circRNA Chr1:30,513,415|30,521,330 was down-regulated and its parental gene LOC_Os01g53090, a pathogen-related protein, was up-regulated by NaCl stress (Table 1). Previous studies reported that LOC_Os01g53090 was up-regulated by both biotic and abiotic stresses such as Magnaporthe oryzae and aluminum stresses (Vijayan et al. 2013;Arbelaez et al. 2017). Another circRNA, Chr2:18,305,592|18,305,907, and its parental gene LOC_Os02g30714, a putative 11-β-hydroxysteroid dehydrogenase that has been proposed to take part in perception and transduction for many environmental stimuli (Wang et al. 2014a, b), were up-regulated (Table 1 and Table 2). In rice, LOC_Os02g30714, which was up-regulated by OsWRKY13, was an activator of rice in resistance to both bacterial and fungal pathogens, and was up-regulated during subsequent recovery after cold stress (Qiu et al. 2008;Yang et al. 2015). Chr4:18,170,122|18,170,775 was down-regulated and its parental gene LOC_Os04g30420, a zinc-binding dehydrogenase, was down-regulated. A previous study showed that LOC_Os04g30420 participated in the response to heat and drought stresses (Wilkins et al. 2016). These analyses reveal that many stress related genes in rice are the potential targets of circRNAs and could be differentially regulated by associated circRNAs in salt-sensitive cultivars.
Similarly, in introgression line 9L136, LOC_ Os01g53090, LOC_Os02g30714, LOC_Os02g39550, and LOC_Os04g30420, the four parent genes of corresponding differentially regulated circRNAs, were also differentially expressed under NaCl treatment, implying the roles of these parent genes and circRNAs in salt stress response in both salt-susceptible and salttolerant rice cultivars (Table 2). Furthermore, many genes that function in ion transportation, such as LOC_ Os02g08010 (Ca 2+ -transporting ATPase, parental gene of Chr2:4,198,743|4,199,218) (Table 2).
Furthermore, LOC_Os05g31254 (calmodulin-related calcium sensor protein gene), a candidate gene for salt tolerance in rice (Wang et al. 2017c), was predicted to be the parent gene of circRNA Chr6:2,046,190|2,059,120 in 93-11, and the parent gene of circRNA Chr10:18,818,949|18,958,248 in 9L-163. In Arabidopsis, the calmodulin-related calcium  sensor protein gene has been reported to modulate stress responses. The functions of LOC_Os05g31254 and Chr10:18,818,949|18,958,248 need to be further studied, since they will be valuable for determining salt stress mechanisms and conducting salt-resistance breeding. In summary, our study reveals the possible roles of rice circRNA in response to salt stress, which will expand our understanding of the characteristics of plant circRNAs and facilitate the determination of salt stress regulatory mechanisms in rice. Furthermore, the complicated relationships between the abundances of circRNAs and their parent genes need to be further explored through molecular biology approaches.