Circular RNAs Repertoire and Expression Profile during Brassica rapa Pollen Development

Circular RNAs (circRNAs) are covalently closed RNA molecules generated by the back-splicing of exons from linear precursor mRNAs. Though various linear RNAs have been shown to play important regulatory roles in many biological and developmental processes, little is known about the role of their circular counterparts. In this study, we performed high-throughput RNA sequencing to delineate the expression profile and potential function of circRNAs during the five stages of pollen development in Brassica rapa. A total of 1180 circRNAs were detected in pollen development, of which 367 showed stage-specific expression patterns. Functional enrichment and metabolic pathway analysis showed that the parent genes of circRNAs were mainly involved in pollen-related molecular and biological processes such as mitotic and meiotic cell division, DNA processes, protein synthesis, protein modification, and polysaccharide biosynthesis. Moreover, by predicting the circRNA–miRNA network from our differentially expressed circRNAs, we found 88 circRNAs with potential miRNA binding sites, suggesting their role in post-transcriptional regulation of the genes. Finally, we confirmed the back-splicing sites of nine selected circRNAs using divergent primers and Sanger sequencing. Our study presents the systematic analysis of circular RNAs during pollen development and forms the basis of future studies for unlocking complex gene regulatory networks underpinning reproduction in flowering plants.


Introduction
Transcriptome-wide sequencing studies have shown that the genomic information is extensively transcribed into protein-coding RNAs and RNAs with no protein coding potential but play widely important regulatory roles. Although alternative splicing and post-transcriptional events can generate numerous mRNA isoforms to increase the proteincoding and regulatory capacity of the eukaryotic genomes, these RNAs comprise only a small percentage (usually about 1-2%) of all RNAs in the cell [1]. The majority of expressed RNAs are non-coding RNAs that, based on their size, can be divided into categories, such as small non-coding and long non-coding RNAs (lncRNAs) [2,3]. These non-coding RNAs have been known to be key regulatory factors controlling gene expression at the posttranscriptional, transcriptional, and epigenetic levels [4]. CircRNAs are another layer of transcriptome complexity that came into focus about a decade ago and are considered noncoding RNAs; however, the protein-coding ability for a subset of them has been reported in recent studies [5,6]. CircRNAs are covalently closed single-stranded RNA molecules (with no 5 caps or 3 poly-A tails) following transcription from a wide range of genomic positions, including exons, introns, and intergenic regions, that undergo a non-canonical splicing event termed as back-splicing [7] (Figure 1). In humans and animals, it was revealed that base pairing from reverse complementary sequences in the flanking introns of circularized exons, as well as some RNA-binding proteins such as Quaking (QKI), promote circRNA production [8]. However, the molecular mechanism behind back-splicing in plants hind back-splicing in plants is still poorly understood [9,10]. Despite being known f decades, circRNAs have been considered as aberrant byproducts of mis-splicing even [10]. With the advances in RNA-seq technology and bioinformatics algorithms, thousan of circRNAs have been identified in the eukaryotic tree of life and revealed to have ce and tissue-specific expression patterns [9,11]. Earlier studies have demonstrated th lncRNAs play important roles in reproductive-related processes such as floral transitio and flower development [12]. Cold-Assisted Intronic Noncoding RNA (COLDAIR) an Cold-Induced Long Antisense Intragenic RNA (COOLAIR) are lncRNAs that control flow ering in Arabidopsis by repressing Flowering Locus C (FLC) through epigenetic mech nisms [13,14]. In rice, it has been revealed that more than 700 long intergenic non-codin RNAs are essential factors for the biogenesis of phased small interfering RNAs that ass ciate with germline-specific Argonaute protein MEL1, which itself is involved in the d velopment of pre-meiotic germ cells, suggesting that rice lncRNAs are important elemen in reproduction [15]. Cell-specific expression of thousands of protein-coding transcripts underpins poll developmental processes [16][17][18][19]. Cis-regulatory elements, tissue-specific promoters, h tone modifications, and DNA methylation have also been reported to be involved in t regulatory network of pollen development [20][21][22][23][24]. LncRNAs have also been implicat in controlling pollen reproductive developmental processes [12]. For example, Long-Da Specific Male-Fertility-Associated RNA (LDMAR) and Zm401 are lncRNAs required f pollen development in rice and maize, respectively [25,26]. A recent study in Brassica ra found that more than 12,000 lncRNAs were expressed during pollen development an fertilization [27]. In addition, recent studies have implicated circRNAs in the regulation pollen development. For instance, in rice, 186 differentially expressed circRNAs have be identified in different pollen developmental stages during the fertility transition of t photo-thermosensitive genic male sterile line [28]. In soybean, 2867 circRNAs were ide Cell-specific expression of thousands of protein-coding transcripts underpins pollen developmental processes [16][17][18][19]. Cis-regulatory elements, tissue-specific promoters, histone modifications, and DNA methylation have also been reported to be involved in the regulatory network of pollen development [20][21][22][23][24]. LncRNAs have also been implicated in controlling pollen reproductive developmental processes [12]. For example, Long-Day-Specific Male-Fertility-Associated RNA (LDMAR) and Zm401 are lncRNAs required for pollen development in rice and maize, respectively [25,26]. A recent study in Brassica rapa found that more than 12,000 lncRNAs were expressed during pollen development and fertilization [27]. In addition, recent studies have implicated circRNAs in the regulation of pollen development. For instance, in rice, 186 differentially expressed circRNAs have been identified in different pollen developmental stages during the fertility transition of the photo-thermosensitive genic male sterile line [28]. In soybean, 2867 circRNAs were identified in flower buds of the cytoplasmic male sterile line (NJCMS1A) and its maintainer (NJCMS1B), of which 1009 circRNAs were differentially expressed between two lines, pointing towards the role of circRNAs in flower and pollen development [29]. Similar work in B. campestris revealed 31 differentially expressed circRNAs between cytoplasm male sterile and fertile lines involved in anther development [30]. To further improve our understanding of the role of circRNAs in pollen development, we performed a time series of RNA-seq experiments during pollen development in B. rapa (AA, 2n = 20). Using three different bioinformatics prediction tools, we identified 1180 circRNAs expressed in five pollen developmental stages. We performed differential gene expression analysis, functional annotation, and pathway enrichment analysis to explore the potential function of circRNAs during pollen development. We also predicted miRNA-circRNA interactions to investigate the potential role of circRNAs as competing endogenous RNAs (ceRNAs) in post-transcriptional gene regulation.

Plant Materials, Growth Conditions, and Sample Collection
In this study, Brassica rapa, accession no. ATC 92270 Y.S (AND)-168, were grown at 21/18 • C day/night under 16/8 h' light/dark (200 µmol m −2 s −1 light intensity) with 60% humidity. Based on the size of floral buds, five pollen developmental stages were identified: pollen mother cells: ≤1 mm buds, tetrad: 1.2-2 mm buds, uninucleate pollen: 2-2.6 mm buds, binucleate pollen: 3-4 mm buds, and mature pollen: 5-6 mm buds ( Figure S1) [27]. The samples were pooled from 15 different plants for each stage of pollen development. Collected buds were immediately transferred to a Petri dish containing modified half strength B5 medium (13% sucrose) on ice. Pollen from the buds in their binucleate and mature pollen stages was released by dissecting out the anthers and squashing them in the modified B5 medium. Because the buds in the first three groups (pollen mother cells, tetrad, and uninucleate pollen) were too small, the whole buds were squashed in the modified B5 medium to release the pollen. The resulting suspensions were then filtered using a 40 µm mesh (pluriStrainer Mini 40µm, pluriSelect, Leipzig, Germany) into 2 mL tubes. Tubes were then centrifuged at 4 • C for 3 min at 150 g, and the supernatant was discarded. The pellet was washed with modified half strength B5 medium, and finally, pure pollen grains were collected from the medium using centrifugation at 4 • C for 3 min at 150 g. After discarding the medium, the pellet was immediately frozen in liquid nitrogen and stored at −80 • C for later use [31].

RNA Extraction and Sequencing
The total RNA was isolated from collected pollen samples using the mirVana™ miRNA Isolation Kit (Thermo-Fisher; Part Numbers AM1560, AM1561, Carlsbad, CA, USA) according to the manufacturer's instructions. To remove contaminated DNA, all the isolated RNA samples were treated with TURBO™ DNase (Ambion, Carlsbad, CA, USA) and then submitted for sequencing. The RNA-seq libraries were then constructed using Illumina TruSeq stranded kit following the rRNA depletion step, and the single-end reads (100 bp) sequencing was performed at the Australian Genome Research Facility (AGRF), Melbourne. The sequencing data were deposited to NCBI's Sequence Read Archive (SRA) under accession number PRJNA763698.

Prediction of Potential circRNA-miRNA Interactions and Visualization
The potential interaction between known miRNAs and our identified differential expressed circRNAs was carried out using web tool psRNATarget (v2, http://plantgrn. noble.org/psRNATarget/home, accessed on 21 August 2020) [50] with default parameters, except we set the Expectation value to 3.0. We selected "Brassica rapa, 157 published miRNA" as input miRNAs and the sequence of our identified differentially expressed circRNAs as the target sequences. The circRNA-miRNA network was visualized using Cytoscape software (v3.8.2) [51]. We also predicted the mRNA targets of miRNAs using Targetfinder [52]. Then, GO enrichment analysis and KEGG pathway enrichment analysis were performed on predicted targets as described in Section 2.5.

Circular RNA Validation
To validate our identified circRNAs that were predicted by bioinformatic approaches, PCR with divergent primers followed by Sanger sequencing was used. Briefly, total RNA from all five pollen developmental stages was reverse transcribed into complementary DNA (cDNA) using SuperScript™ III Reverse Transcriptase (Invitrogen, Carlsbad, CA, USA) and random hexamers according to the manufacturer's protocol. The sequences of 15 in silico predicted circRNAs that were predicted to be expressed in developing pollen were used to design 15 pairs of divergent primers using Primer3web (v4.1.0) [53]. The PCR procedure was as follows: initial step at 94 • C for 3 min; followed by 40 cycles at 94 • C for 30 s, appropriate annealing temperature (Table S1) for 45 s, and 72 • C for 30 s; and then 1 cycle at 72 • C for 7 min. The PCR products were then visualized by agarose gel (1.5%) electrophoresis, and the desired bands were recovered from the gel using Wizard ® SV Gel and PCR Clean-Up System (Promega, Madison, WI, USA). Finally, the PCR products were cloned into a pJET1.2/blunt vector using the CloneJET PCR Cloning Kit (Thermo Scientific, Vilnius, Lithuania) and subjected to Sanger sequencing to confirm the back-spliced junction sites.

Identification and Characterization of circRNAs
To explore the involvement of circRNAs in B. rapa pollen development, rRNA-depleted RNA libraries for five pollen developmental stages (pollen mother cells, tetrad, uninucleate pollen, binucleate pollen, and mature pollen) were constructed. A total of about 23 gigabytes of high-quality clean data with more than 233 million reads was obtained (Table S2). After analyzing the sequencing data with three circRNA identification tools, more than 2000 circRNAs were obtained (Table S3). Out of the total number of circR-NAs, we identified 1180 unique circRNAs in all five pollen developmental stages with 824 circRNAs for find_circ, 380 circRNAs for CIRI2, and 137 circRNAs for CIRCExplorer2 ( Figure S2 and Table S4). When comparing the number of circRNAs found by all three tools, only six circRNAs were common between all, and more circRNAs were common between CIRI2 and find_circ ( Figure S2). Among 1180 identified circRNAs, 1034 (87.63%) circRNAs were predicted to have at least one exon from protein-coding genes (Figure 2A), 85 (7.20%) circRNAs were generated from intergenic regions, and of the remaining cir-cRNAs, 61 (5.17%) were intronic. While 133 genes could produce two or more types of circRNAs through alternative back-splicing, 737 genes produced only one circRNA isoform ( Figure 2B). Moreover, the parent genes of circRNAs were unevenly distributed on different chromosomes; chromosome nine with 184 circRNAs and chromosome four with 63 cir-cRNAs produced the most and the least circRNAs, respectively, and 21 circRNAs were produced from four scaffolds ( Figure 2C). The length distribution of circRNAs was mainly distributed between 100 to 600 bp ( Figure 2D). Finally, we visualized circRNA expression vs. pollen developmental stages using a clustering heatmap (

The Expression Patterns of circRNAs in B. rapa Pollen Development
To investigate whether the circRNAs are expressed in a specific manner during B. rapa pollen development, we compared the expression patterns of circRNAs in each pollen developmental stage with its previous stage (Table S5). The results show that 966 circRNAs had a significant differential expression, of which 367 circRNAs displayed stage-specific expression patterns. Comparing the expression of circRNAs between stages, a more distinct expression was observed when pollen developed from binucleate pollen to mature pollen or developed from pollen mother cell to tetrad. ( Figure 4A). Additionally, more circRNAs tended to up-regulate during pollen development, except when pollen developed from uninucleate pollen to binucleate pollen, where more circRNAs were downregulated ( Figure 4B). Finally, we visualized the expression profile of the 60 most significant (α ≥ 0.001) differential expressed circRNAs (Tables S5 and S6) among different stages of pollen development ( Figure S3). The expression analysis results indicate that circRNAs have a distinctive expression pattern during pollen development, pointing towards their defined stage-specific roles in pollen developmental progression.

The Expression Patterns of circRNAs in B. rapa Pollen Development
To investigate whether the circRNAs are expressed in a specific manner d rapa pollen development, we compared the expression patterns of circRNAs in eac developmental stage with its previous stage (Table S5). The results show that 96 stages of pollen development ( Figure S3). The expression analysis results indicate that circRNAs have a distinctive expression pattern during pollen development, pointing towards their defined stage-specific roles in pollen developmental progression.

Conservation Analysis of circRNAs between B. rapa and other Plant Species
The comparison of identified circRNAs in this study with the collection of circRNAs obtained from PlantcircBase showed that 38.71% of B. rapa circRNAs were homologous to the circRNAs in the database of which more than 35% were homologous to the Arabidopsis circRNAs (Table 1).

Functional Annotation Analysis of CircRNAs Parent Genes
To explore the potential functions of circRNAs during pollen development, we performed GO enrichment analysis on the parental genes of all the differential expressed circRNAs. The results show that circRNA parent genes belonged to three GO categories: biological process, molecular function, and cellular component (Table S7). For the biological process, the circRNA-host genes were enriched in a variety of metabolic, catabolic, and cell processes such as protein folding, glucan catabolic process, regulation of phosphorus metabolism, meiotic DNA double-strand break formation process, meiotic cell cycle process, meiosis I cell cycle process, and ncRNA metabolic process. In the molecular function, the enriched GO terms included GTPase activity, calmodulin binding, catalytic activity acting on RNA, kinase regulator activity, and metal ion binding. In the cellular component category, only four terms, cytoplasm, endoplasmic reticulum, microtubuleassociated complex, and katanin complex, were enriched ( Figure 5A,B).
To further investigate the function of circRNA host genes, Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was conducted. The results reveal the enrichment of 23 significant pathways (Table S8), including cyanoamino acid metabolism, protein processing in the endoplasmic reticulum, peroxisome, and aminoacyl-tRNA biosynthesis, among others ( Figure 5C).

Conservation Analysis of circRNAs between B. rapa and other Plant Species
The comparison of identified circRNAs in this study with the collection of circRNAs obtained from PlantcircBase showed that 38.71% of B. rapa circRNAs were homologous to the circRNAs in the database of which more than 35% were homologous to the Arabidopsis circRNAs (Table 1).

Functional Annotation Analysis of CircRNAs Parent Genes
To explore the potential functions of circRNAs during pollen development, we performed GO enrichment analysis on the parental genes of all the differential expressed circRNAs. The results show that circRNA parent genes belonged to three GO categories: biological process, molecular function, and cellular component (Table S7). For the biological process, the circRNA-host genes were enriched in a variety of metabolic, catabolic, and cell processes such as protein folding, glucan catabolic process, regulation of phosphorus metabolism, meiotic DNA double-strand break formation process, meiotic cell cycle process, meiosis I cell cycle process, and ncRNA metabolic process. In the molecular function, the enriched GO terms included GTPase activity, calmodulin binding, catalytic activity acting on RNA, kinase regulator activity, and metal ion binding. In the cellular component category, only four terms, cytoplasm, endoplasmic reticulum, microtubule-associated complex, and katanin complex, were enriched ( Figure 5A,B).
To further investigate the function of circRNA host genes, Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was conducted. The results reveal the enrichment of 23 significant pathways (Table S8), including cyanoamino acid metabolism, protein processing in the endoplasmic reticulum, peroxisome, and aminoacyl-tRNA biosynthesis, among others ( Figure 5C).

Prediction of miRNA Target Sites in circRNAs
To further evaluate the function of circRNAs in the post-transcriptional regulation of genes in pollen development, we predicted the binding sites of miRNAs in the sequences of differential expressed circRNAs (Table S9). The results show that 88 circRNAs contained putative miRNA-binding sites for 73 miRNAs. Of these 88 circRNAs, 49 (55.68%) had only one miRNA-binding site, followed by 22 circRNAs (25%) with two miRNA-binding sites, and the remaining circRNAs (19.32%) contained binding sites for three to nine miRNAs ( Figure 6). For example, circRNA A06:606103-667861 and A09:10484746-10513688 had binding sites for eight and nine miRNAs, respectively. MiRNAs could be targeted by one or more circRNAs as well. Among 74 identified miRNAs, 27 (36.49%) were targeted by a single circRNA, 24 miRNAs (32.43%) were targeted by two circRNAs, and the remaining miRNAs (31.08%) were targeted by three to twelve circRNAs ( Figure  6). For instance, bra-miR9569-5p, bra-miR5716, and bra-miR9563a-3p were targeted by seven, ten, and twelve circRNAs, respectively.

Prediction of miRNA Target Sites in circRNAs
To further evaluate the function of circRNAs in the post-transcriptional regulation of genes in pollen development, we predicted the binding sites of miRNAs in the sequences of differential expressed circRNAs (Table S9). The results show that 88 circRNAs contained putative miRNA-binding sites for 73 miRNAs. Of these 88 circRNAs, 49 (55.68%) had only one miRNA-binding site, followed by 22 circRNAs (25%) with two miRNA-binding sites, and the remaining circRNAs (19.32%) contained binding sites for three to nine miRNAs ( Figure 6). For example, circRNA A06:606103-667861 and A09:10484746-10513688 had binding sites for eight and nine miRNAs, respectively. MiRNAs could be targeted by one or more circRNAs as well. Among 74 identified miRNAs, 27 (36.49%) were targeted by a single circRNA, 24 miRNAs (32.43%) were targeted by two circRNAs, and the remaining miRNAs (31.08%) were targeted by three to twelve circRNAs ( Figure 6). For instance, bra-miR9569-5p, bra-miR5716, and bra-miR9563a-3p were targeted by seven, ten, and twelve circRNAs, respectively. We also predicted miRNA targets and performed functional enrichment analysis to investigate what processes in the cell could be affected by the interaction of circRNAs and miRNAs. The results show that the 73 predicted miRNAs could interact with 741 genes (Table S10) belong to various GO terms, mainly DNA processes and replication, amino acid metabolic processes, and gene expression processes (Table S10, Figure 7A,B). KEGG pathway enrichment analysis showed that miRNA target genes are involved in multiple pathways such as metabolic pathways, biosynthesis of secondary metabolites, and RNA transport, among others (Table S11, Figure 7C). We also predicted miRNA targets and performed functional enrichment analysis to investigate what processes in the cell could be affected by the interaction of circRNAs and miRNAs. The results show that the 73 predicted miRNAs could interact with 741 genes (Table S10) belong to various GO terms, mainly DNA processes and replication, amino acid metabolic processes, and gene expression processes (Table S10, Figure 7A,B). KEGG pathway enrichment analysis showed that miRNA target genes are involved in multiple pathways such as metabolic pathways, biosynthesis of secondary metabolites, and RNA transport, among others (Table S11, Figure 7C).

Validation of circRNAs Using Divergent Primers
To confirm the accuracy of the identified circRNAs, we selected 15 circRNAs for experimental validation using PCR and Sanger sequencing. Using a set of divergent primers designed for each circRNA, we successfully amplified nine circRNAs using nine pairs of divergent primers. Then, all the amplified PCR products were further validated by sequencing to confirm the presence of the back-spliced junctions (Figure 8). Our validated circRNAs are composed of one to three exons, except for circRNA A03:5712452-5713124, which has an intron in its conjunction site.

Validation of circRNAs Using Divergent Primers
To confirm the accuracy of the identified circRNAs, we selected 15 circRNAs for experimental validation using PCR and Sanger sequencing. Using a set of divergent primers designed for each circRNA, we successfully amplified nine circRNAs using nine pairs of divergent primers. Then, all the amplified PCR products were further validated by sequencing to confirm the presence of the back-spliced junctions (Figure 8). Our validated circRNAs are composed of one to three exons, except for circRNA A03:5712452-5713124, which has an intron in its conjunction site.

Discussion
There is increasing evidence that circular RNAs are not simply rare splicing events but are highly regulated molecules with unique selected sequences from different genomic regions [54]. In animals, circRNAs can carry out important biological functions in

Discussion
There is increasing evidence that circular RNAs are not simply rare splicing events but are highly regulated molecules with unique selected sequences from different genomic regions [54]. In animals, circRNAs can carry out important biological functions in a cell-type, tissue-, or developmental-stage-specific manner [55]. Although most of the circRNAs are expressed at a low level compared with their cognate mRNAs, they can play important roles because they are naturally resistant to diverse exonucleases, and their half-life is much longer than mRNAs [54,56]. To date, thousands of distinct circRNAs have been found in several taxa, including different species of plants [57,58]. Various studies have suggested that plant circRNAs are closely associated with growth and development [59][60][61].
In the present study, we investigated the expression and potential function of circRNAs during pollen development in B. rapa. According to the comparison studies, circRNA detection tools have been designed with different detection algorithms, and no single method provides all the metrics needed for circRNA identification [36,[62][63][64]. Here, we used CIRI2, find_circ, and CIRCExplorer2, the top three identification tools for circRNA analysis in the literature [36], to achieve a better overview of the circRNA expression profile in our datasets. After combining the results, only six identified circRNAs were common between the three tools, which is normal as each tool uses a different aligner and metrics for identifying circRNAs in a given dataset [64]. Different tools, especially when compared pairwise, may share a subset of identified circRNAs, but in most cases, this number decreases after combining the results from more than two tools [64]. In total, we identified 1180 circRNAs in all five stages of pollen development. Compared with previous reports, this number is similar to other plant species such as Arabidopsis (970) [65], Gossypium arboretum (1041), and Gossypium raimondii (1478) [66]. We found that during pollen development, circRNAs originated from different genomic regions, but most of our identified circRNAs (87.63%) were generated from exonic regions, which is consistent with previous studies on many plant species such as Arabidopsis (85.1%) [67], Chinese cabbage (70.49%) [68], and Brassica campestris (63.13%) [30]. Most of the genes can produce only one isoform of circRNAs; however, multiple isoforms of circRNAs produced through alternative back-splicing from a single gene have also been observed [69][70][71]. Here, we found that most of the circRNA parent genes (about 83%) generate only one circRNA isoform, with the remaining 17% of the genes able to produce two or more circRNA isoforms. Conservation between different plant species is another feature of circRNAs that has been shown in multiple studies [67,[72][73][74]. Although more than 37% of our identified circRNAs were found to be conserved between B. rapa and other plant species, the majority of them were conserved between B. rapa and Arabidopsis, the two species from the same family with closely related genomes [75]. When Arabidopsis was not taken into account, only about 3% of B. rapa circRNAs were conserved with circRNAs from other plant species. The low conservation rate could be due to the differences between datasets, as they were prepared from various tissues and/or under stress conditions, or the various bioinformatics tools that have been used for analyzing those datasets [73].
The expression of circRNAs changes during different stages of growth and development. For instance, by analyzing the expression profile of circRNAs in various plant species, it was revealed that circRNAs are closely involved in leaf growth in Arabidopsis [59], fruit ripening in pepper [60], and pollen development in B. campestris [30], rice [28], and soybean [29]. Herein, our differential expression analysis revealed that the stage-specific expression pattern of circRNAs changes during pollen development. We found 130 differentially expressed circRNAs while pollen mother cells developed to tetrad, and 191 differentially expressed circRNAs while binucleate pollen developed to mature pollen. We noted that the most common circRNAs were shared between two closely related stages. For example, 272 shared differentially expressed circRNAs when tetrads developed to microspores and then to binucleate pollens. These findings indicate that circRNA expression is related to different pollen developmental stages, and they have the potential to play important roles during pollen development in B. rapa.
Previous studies showed that circRNAs could regulate various processes such as chromatin structure, gene expression, translation, and even cell division [76][77][78]. For example, centromeric-retrotransposon-derived circRNAs in maize can bind to centromeres via R-loops to promote chromatin looping in centromere regions [76]. Exon-intron-containing circRNAs in human cells can increase the transcription of their parental genes by interacting with U1 small nuclear ribonucleoproteins at the promoters [77], or in Arabidopsis, circRNAs can form a strong RNA-DNA hybrid, resulting in transcriptional pausing [78]. CircRNAs also can form circRNA-protein structures and affect translation and cell cycle regulation. One such example is CircPABPN1, which can suppress HuR from binding to its cognate mRNA, resulting in reduced translation of PABPN1 mRNA [79]. Another example is circFOXO3, which has multiple binding sites for cell-cycle-regulating proteins such as p53, CDK-2, and p21. It was reported that circFOXO3, p21, and CDK-2 could form a ternary structure and inhibit CDK-2/Cyclin-E complex activation, which is necessary for G1 to S phase transition [80]. During pollen development, several mitotic and meiotic cell divisions occur, which involve the expression of a few thousand genes [16,81]. Studies showed that during pollen development, the tapetum endoplasmic reticulum was highly involved in biosynthesis, folding, and secreting proteins [82,83]. For pollen wall formation, lipidic components and polysaccharides are required as tapetum secretes lipid components onto the pollen surface [84,85], or β-glucosidase, which is involved in the regulation of polysaccharide metabolism, downregulated in the sterile floral buds of B. rapa [86]. As circRNAs' functions may be related to their parent gene's function, we annotated the biological roles of circRNAs' parent genes using GO and KEGG analysis to better understand the potential function of circRNAs in pollen development. We noted that circRNA parent genes could be involved in various important processes and functions related to pollen development such as protein synthesis and fate, cell cycle and DNA processing, kinase and phosphorus activities, polysaccharide metabolism, gene silencing, RNA processing, and antiporter activities. KEGG analysis showed that the circRNAs' parent genes are involved in 23 pathways, mainly amino acid metabolism and protein processing, lipid biosynthesis and metabolism, metabolic pathways, carbon metabolism, and RNA degradation. Similar results were observed when the function of miRNA target genes was investigated. Based on our results and previous findings, we can assume that circRNAs might play vital functions during pollen development in B. rapa.
Small RNAs such as miRNAs can regulate gene expression through mRNA cleavage, translational repression, or gene silencing via miRNA-directed DNA methylation [87]. Studies on animal and human cells have reported that circRNAs can bind to miRNAs and sequester them from their target mRNAs to regulate gene expression at the posttranscriptional level [88]. For example, it has been shown that circSry, associated with testis development in mice, contains 16 binding sites for miR-138 [89], and CDR1as, which is a highly expressed circRNA in mammalian brains, contains more than 70 binding sites for miR-7. Studies showed that CDR1as regulates gene expression by acting as miR-7 storage (sponge) and ensuring the release of an appropriate amount of miR-7 to the target mRNAs [8,90,91]. Prediction of circRNA-miRNA networks in plant species revealed that plant circRNAs, compared with animals and humans, have fewer interactions with miRNAs, suggesting that the main function of plant circRNAs might not be an miRNA decoy [92]. Nevertheless, several studies proposed that plant circRNAs could target miRNAs by acting as competing endogenous RNAs in post-transcriptional regulation of the genes [28][29][30]93,94]. For example, circRNAs were proposed to function as miRNA sponges in flower development in Arabidopsis [93], anther development in B. campestris [30], and pollen development in rice and soybean [28,29]. A database named "GreenCircRNA" collected all the identified circRNAs in different plant species and predicted circRNA-miRNAs interactions for each species [95]. In the present study, we found 88 circRNAs with putative miRNA binding sites. A few of them had multiple binding sites for the same or different miRNAs, suggesting that circRNAs could interact with miRNAs or act as miRNA sponges to regulate pollen development. We also noted some well-known and important miRNAs in our predicted network, such as miR156, miR157, miR158, miR161, miR162, miR164, miR172, miR395, and miR396. These miRNAs proved to be involved in many biological and developmental processes such as vegetative growth, flowering, fertility, and fruit ripening [87,[96][97][98]. For example, we found nine circRNAs with binding sites for miR156, seven circRNA with binding sites for miR172, four circRNAs with binding sites for miR396, and three circRNAs with binding sites for miR164, which all are known to regulate flower development in Arabidopsis [99,100], barley [101], tomato [102], Brassica napus [103], and strawberry [104]. In addition, we found nine circRNAs with binding sites for miR158 that were previously reported as key regulatory miRNA during pollen development in Brassica campestris [105]. Accordingly, our results suggest that circRNAs could act as ceRNAs during pollen development in B. rapa.

Conclusions
In conclusion, by studying the profile of circRNA expression during pollen development in B. rapa, we identified 1180 circRNAs, of which 367 circRNAs showed stage-specific expression. Functional characterization of circRNA host genes revealed that circRNAs were mainly related to biological and molecular processes in pollen development such as mitosis and meiosis cell cycles, protein biosynthesis, protein modification, and polysaccharide processes. Moreover, 88 circRNAs were found to contain miRNA-binding sites, suggesting the role of circRNAs in gene expression as post-transcriptional regulatory elements. Our study revealed the potential functions of circRNAs during pollen development and paves the way for further experiments on studying the molecular mechanism of these new regulatory RNA molecules.