Multi-omics of Circular RNAs and Their Responses to Hormones in Moso Bamboo (Phyllostachys edulis)

Abstract

Circular RNAs (circRNAs) are endogenous non-coding RNAs with covalently closed structures, which have important functions in plants. However, their biogenesis, degradation, and function upon treatment with gibberellins (GAs) and auxins (1-naphthaleneacetic acid, NAA) remain unknown. Here, we systematically identified and characterized the expression patterns, evolutionary conservation, genomic features, and internal structures of circRNAs using RNase R-treated libraries from moso bamboo (Phyllostachys edulis) seedlings. Moreover, we investigated the biogenesis of circRNAs dependent on both cis- and trans-regulation. We explored the function of circRNAs, including their roles in regulating microRNA (miRNA)-related genes and modulating the alternative splicing of their linear counterparts. Importantly, we developed a customized degradome sequencing approach to detect miRNA-mediated cleavage of circRNAs. Finally, we presented a comprehensive view of the participation of circRNAs in the regulation of hormone metabolism upon treatment of bamboo seedlings with GA and NAA. Collectively, our study provides insights into the biogenesis, function, and miRNA-mediated degradation of circRNAs in moso bamboo.

Keywords: Alternative splicing; Circular RNA; Degradome; Phyllostachys edulis; Phytohormone.

Figure 1

Characterization of circRNAs in moso bamboo A. Flow chart for multi-omics sequencing including mRNA sequencing, circRNA sequencing, degradome sequencing, and small RNA sequencing upon hormone treatment. B. Venn diagram showing the number of circRNAs detected in different samples (upper panel). The density plot of back-splicing reads supporting circRNAs (lower panel). C. Validation of circRNAs using RT-PCR after RNase R treatment with linear RNA (NTB) as control (left panel). Sanger sequencing validated back-splicing junctions of circ-RAD16, which was not resistant to RNase R (right panel). D. Number (left) and overlap (right) of circRNAs in Arabidopsis thaliana and Oryza sativa L. that show homology to circRNAs in moso bamboo. E. Multiple sequence alignment of conserved circ-GSL1. Asterisk symbols indicate highly conservative nucleotides. F. Diagrams of different AS types (left panel) and percentage of these circRNAs (right panel). Gray bars and black lines represent exons and introns, respectively. Dotted curves and colored bars indicate AS events. Colored arced lines represent back-splicing (circularization). G. GO enrichment analysis of circRNAs with AS events. Darker orange node color indicates more significant P values. The circle size is proportional to the number of genes enriched in the terms. The arrows represent hierarchical relations between GO terms. circRNA, circular RNA; AS, alternative splicing; GO, Gene Ontology; A3SS, alternative 3′ splice site; A5SS, alternative 5′ splice site; ExonS, exon skipping; IntronR, intron retention; A3BS, alternative 3′ back-splice site; A5BS, alternative 5′ back-splice site; MIC, mutually inclusive circRNA; mRNA, messenger RNA; RT-PCR, reverse transcription-polymerase chain reaction; GA, gibberellin; NAA, 1-naphthaleneacetic acid.

Figure 2

Biogenesis of circRNAs in moso bamboo A. Distribution of ICSs (brown arrows) in flanking introns of ecircRNAs. Circularized exon and flanking intron are indicated by green bar and gray line, respectively. B. Schematic drawings show expression vectors including flanking ICSs (brown arrows), circularized exon (green bar), and linear exon (yellow bar). PCR primers for circRNAs are indicated by black arrows. Semi-quantitative RT-PCR in the lower panel shows circularization efficiency for circularized exon and linear exon, respectively. C. Experimental detection of the function of introns in multiple exon circularization of circ-PKL1. Schematic drawings of expression vectors including or excluding interior introns. Semi-quantitative RT-PCR in the lower panel shows circularization efficiency of circ-PKL1 with or without interior intron. D. Percentage of circRNAs showing high PCC (r ≥ 0.5 or r ≤  −0.5) with proteins SF3A, SF3B, FUS, HNRNPL, QKI, NF90, and DHX9. E. RNAi knockdown of PedQKI. Schematic drawing in the upper panel shows RNAi expression vector to knock down expression of PedQKI. The lower panel shows PCR-based validation of the expression of six randomly selected circRNAs in the QKI RNAi sample. ATCB is a housekeeping gene used as a control.  F. PCCs between circRNAs and 17 proteins including 1, 2, and 14 core spliceosomal factors, splicing factors, and other RBPs, respectively. Number in the right panel indicates percentage of circRNAs with high PCC to those proteins. G. Distribution of PCCs between ecircRNAs/ciRNAs and linear RNAs generated from the same host genes. ecircRNA, exonic circRNA; ICS, inverted complementary sequence; EV, empty vector; OV, overexpression vector; NOS, nopaline synthase; PCC, Pearson correlation coefficient; RNAi, RNA interference; RBP, RNA-binding protein; ciRNA, circular intronic RNA.

Figure 3

circRNAs regulate AS in moso bamboo A. Overview of overlapping regions between AS events and ecircRNAs. White bars indicate exons, black bars indicate introns, blue arrows indicate AS events, and colored arced lines indicate back-splicing junctions. B. Overview of overlapping regions between AS events and ciRNAs. C. The histograms show the number of AS events in the transcribed regions of ecircRNAs (brown bars) and random RNA segments (gray bars) in three species. D. The histograms show the number of AS events in the transcribed regions of ciRNAs (light yellow bars) and random RNA segments (gray bars) in three species. E. Visualization of IntronR and ExonS events in transcribed regions of circ-BRE1 and circ-BRE2, respectively. F. RT-PCR validation of circ-BRE1-1 and circ-BRE1-2 and their corresponding AS events, linear BRE1-1 and linear BRE1-2. Divergent arrows represent divergent primers, and convergent arrows represent convergent primers. Linear RNA of ACTB was used as a control. G. The histogram plot shows the overlap between R-loop and AS events in transcribed regions of ecircRNAs. Black arrow in the left panel indicates the observed number of ecircRNAs located in R-loop regions. H. Root phenotype upon DMSO and CPT treatment (left panel). Semi-quantitative RT-PCR shows the expression of circ-BRE1-1 and linear BRE1-1 upon CPT treatment (right panel). R-loop, RNA:DNA hybrid; AltA, alternative acceptor; AltD, alternative donor; DMSO, dimethyl sulfoxide; CPT, camptothecin.

Figure 4

circRNAs regulate miRNA-associated genes in moso bamboo A. The upper panel shows the miRNA-associated genes involved in several major steps in miRNA biogenesis and modes of action in plants. The table in the lower panel shows the number of miRNA-associated genes and their corresponding circRNAs. B. RT-PCR and Sanger sequencing for validation of circ-DCL4 with RNase R treatment. C. The left panel shows the vector construction for circ-DCL4 and linear DCL4. RT-PCR validation in the right panel shows the overexpression of circ-DCL4 and expression of linear DCL4. miRNA, microRNA; HEN1, small RNA 2′-O-methyltransferase; CDKF, cyclin-dependent kinase F-1; HST, shikimate O-hydroxycinnamoyltransferase.

Figure 5

miRNA-mediated cleavage of circRNAs in moso bamboo A. Flow chart for construction of customized degradome libraries for enriching the decaying circRNA without poly(A) tails. B. The box line plot shows the distribution of reads in 5′ UTR, CDS, and 3′ UTR of the gene from poly(A)+ and poly(A)− degradome libraries. C. The UpSet plot shows the intersection of cleavage sites among different types of libraries and different hormone treatments. The histogram plot in the upper panel shows the number of cleavage sites shared by different libraries. The histogram plot in the lower panel and the line plot in the right panel show total cleavage sites from each library and the combination of different libraries. D. Plot in the upper panel shows the distribution of cleavage sites from two types of libraries in genes divided into three groups: the upstream region of circRNA body, downstream region of circRNA body, and circRNA body. Bar in lower panel shows the percentage of cleavage sites from two types of libraries in transcript regions of circRNAs. E. Schematic overview of decaying reads of circRNAs spanning back-splicing sites. Gray bars represent exons, black lines represent introns, black arrows represent cleavage sites, and colored arced lines represent back-splicing. Blue bars and dash arced lines indicate the mapped back-splicing junction reads. F. Length distribution between cleavage sites and the 3′ end of transcribed regions of circRNAs. G. Venn diagram shows the number of three subgroups. Subgroups 1, 2, and 3 represent the number of circRNAs including cleavage sites based on degradome sequencing, the number of circRNAs with a distance of < 47 nt between cleavage sites and the 3′ end of transcriptional regions of circRNAs, and the number of circRNAs determined by degradome reads spanning back-splicing sites. H. Number of circRNAs and linear RNAs including miRNA-mediated cleavage. I. Schematic overview of miR166-mediated cleavage sites in circ-NHLRC2. J. The upper panel shows the vector construction for overexpressing circ-NHLRC2 only (OV1) or both circ-NHLRC2 and miR166 (OV2). The lower panel indicates the expression of circ-NHLRC2 and miR166 detected by RT-PCR. Divergent and convergent arrows represent divergent and convergent primers, respectively. rRNA, ribosomal RNA; CDS, coding sequence; UTR, untranslated region; TSS, transcription start site; TTS, transcript termination site.

Figure 6

Hormone-induced circRNAs in moso bamboo A. The scatterplot in the left panel shows RPMs for H₂O (X-axis) and GA treatments (Y-axis). Semi-quantitative PCR using divergent primers in the right panel validated the differential expression data of circRNAs in response to GA treatments with linear ACTB as internal reference gene using convergent primers. B. The scatterplot shows the differential levels of circRNAs in response to NAA, and semi-quantitative PCR shows the validation of the levels of circRNAs based on sequencing. C. Venn diagram shows overlapping and unique differential circRNAs in response to GA and NAA. D. The histogram plot shows the top 10 GO terms enriched for hormone-induced circRNAs. E. Heatmap showing expression levels of circular transcripts related to cell wall, cellulose, and lignin. Red or blue represents high and low abundance of circRNAs, respectively. Image and bar chart in the right panel show the phenotype and height of seedlings upon GA and NAA treatments after 2 weeks. F. The PCCs of circRNAs and their host linear RNAs related to fast growth, NAA, and GA. G. The vector construction (left panel) and RT-PCR validation (right panel) of circ-CSLA1 and linear CSLA1. Divergent and convergent arrows represent divergent and convergent primers, respectively. H. RT-PCR using divergent primers revealed the expression levels of circRNAs in six transformed lines (T1 generation). I. The histogram shows the plant height of 24-day-old rice seedlings transformed by circ-SPY, circ-MYBS3, circ-WRKY4, circ-CSLA1, circ-AGO1A, and circ-GID1, respectively. J. The phenotypes of plant height of rice seedlings transformed by circ-AGO1A and circ-GID1. RPM, reads per million; FPKM, fragments per kilobase of exon model per million mapped fragments; WT, wild type; OE, Over expression.

Figure 7

The potential interplay between hormone and circRNA metabolism Module I: circRNAs originating from hormone-related genes exhibit dynamic expression upon exposure to GA and NAA. Module II: hormones could affect circRNA biogenesis by regulating several splicing factors. Module III: hormones could regulate the degradation of circRNAs though modulating RSC. Module IV: functional circRNAs might regulate hormone metabolism by regulating splicing and transcription of their linear cognates or generating potential proteins from translatable circRNAs. RISC, RNA-induced silencing complex; RSC, RNA silencing complex.