Whole-transcriptome analysis and construction of an anther development-related ceRNA network in Chinese cabbage (Brassica campestris L. ssp. pekinensis)

Anther development is precisely regulated by a complex gene network, which is of great significance to plant breeding. However, the molecular mechanism of anther development in Chinese cabbage is unclear. Here, we identified microRNAs (miRNAs), mRNAs, long non-coding RNAs (lncRNAs), and circular RNAs (circRNAs) related to anther development in Chinese cabbage (Brassica campestris L. ssp. pekinensis) to construct competitive endogenous RNA (ceRNA) regulatory networks and provide valuable knowledge on anther development. Using whole-transcriptome sequencing, 9055, 585, 1344, and 165 differentially expressed mRNAs (DEmRNAs), miRNAs (DEmiRNAs), lncRNAs (DElncRNAs), and circRNAs (DEcircRNAs) were identified, respectively, in the anthers of Chinese cabbage compared with those in samples of the vegetative mass of four true leaves. An anther-related ceRNA regulatory network was constructed using miRNA targeting relationships, and 450 pairs of ceRNA relationships, including 97 DEmiRNA–DEmRNA, 281 DEmiRNA–DElncRNA, and 23 DEmiRNA–DEcircRNA interactions, were obtained. We identified important genes and their interactions with lncRNAs, circRNAs, and miRNAs involved in microsporogenesis, tapetum and callose layer development, pollen wall formation, and anther dehiscence. We analyzed the promoter activity of six predominant anther expression genes, which were expressed specifically in the anthers of Arabidopsis thaliana, indicating that they may play an important role in anther development of Chinese cabbage. This study lays the foundation for further research on the molecular mechanisms of anther growth and development in Chinese cabbage.


Results
Identification of miRNAs and their functions. To fully understand the miRNA repertoire related to anther development in Chinese cabbage, ' Ant' and 'Mix' sRNA libraries were constructed and sequenced. A total of 24,706,757 and 22,291,255 raw reads were obtained from the ' Ant' and 'Mix' sRNA libraries, and after filtering, 13,351,418 and 12,257,469 valid reads were obtained, respectively (Table S1). Based on original sequencing data analysis and statistics, we calculated the length distribution of the filtered valid data to be distributed in the range of 20-24 nt, which conforms to the typical characteristics of Dicer enzyme cutting (Table S2). A total of 1360 miRNAs (709 known and 651 novel) were identified (Table S3). Family analysis of detected miRNAs to explore the existence of miRNA families in other species aids in understanding the conservation of miRNAs in evolutionary relationships. A total of 65 miRNA families conserved in plants were identified in this study (Table S4).
To further understand the potential function of miRNAs, the target genes of miRNAs were analyzed using GO and KEGG (Tables S7 and S8). GO annotation revealed that most DEmRNAs were annotated to the 'biological process' (BP) ontology. Additionally, common annotations were 'regulation of transcription, DNA-templated, ' 'transcription' , and 'protein phosphorylation' in the BP ontology; 'nucleus, ' 'integral component of membrane' , 'cytoplasm' , and 'plasma membrane' in the 'cellular component' (CC) ontology; and 'molecular function' , ' ATP binding' , 'sequence-specific DNA binding transcription factor activity' , and 'DNA binding' in the 'molecular function' (MF) ontology (Fig. 1c). KEGG enrichment analysis revealed that the targets of DEmiRNAs were significantly enriched in several pathways that may be related to anther development, including 'plant hormone signal transduction' , which plays an important regulatory role in plant anther development 31 , 'endocytosis' , which is essential for the development of male reproductive organs in Arabidopsis 32 , 'plant-pathogen interaction' , which plays a dominate role in recessive genic male sterility of cabbage 33 , 'protein processing in the endoplasmic reticulum' , and 'oxidative phosphorylation' , which were also significantly enriched in the comparative transcriptome analysis of tobacco sua-cytoplasmic male sterility, indicating that they may be involved in anther development 34 (Fig. 1d).  (Table S9). The mRNA and lncRNA expression levels were estimated using the FPKM value of Cuffdiff. A total of 2384 lncRNAs and 33,821 mRNAs were identified (Table S3). Compared with those of 'Mix, ' 1344 lncRNAs and 9055 mRNAs were differentially expressed in the anthers of Chinese cabbage. Among the DElncRNAs and DEmRNAs, 1230 lncRNAs and 5177 mRNAs were upregulated, while 114 lncRNAs and 3878 mRNAs were downregulated in the anther (Table S5, Fig. 2a, d). Next, we compared the identified lncRNAs with the genomic characteristics of Chinese cabbage protein-coding genes. The structural characteristics and expression levels of lncRNAs and mRNAs were significantly different (Fig. S1). Most transcript lengths of mRNA were longer than 1000 bp, while those of lncRNAs were less than 500 bp. The average ORF length of mRNAs was longer than that of lncRNAs. The ORF length of most lncRNAs was 0-50 aa, while the ORF length of most mRNAs was 100-500 aa. Compared with those of mRNAs, lncRNAs had fewer exons, on average, and most lncRNAs had only one or two exons. Additionally, mRNA expression levels were higher than those of lncRNAs.
GO enrichment and KEGG analyses were performed to explore the functions of the DEmRNAs (Tables S7 and  S8). DEmRNAs were commonly annotated to the GO terms 'regulation of transcription, DNA-templated' , 'oxidation-reduction process' , 'protein phosphorylation' , and 'defense response' in the BP ontology; 'nucleus' , 'integral component of membrane' , 'plasma membrane' , and 'chloroplast' in the CC ontology; and 'protein binding' , ' ATP binding' , 'DNA-binding transcription factor activity' , and 'DNA binding' in the MF ontology (Fig. 2e). Additionally, many GO terms related to pollen development, including 'pollen tube growth' , 'carbohydrate-binding' , 'plant-type cell wall' , 'lipid transport' , and 'cell wall modification' , were significantly enriched. The KEGG analysis  www.nature.com/scientificreports/ revealed that the following pathways were highly enriched in the DEmRNAs: 'photosynthesis' , 'photosynthesisantenna proteins' , 'phenylpropanoid biosynthesis' , 'glucosinolate biosynthesis' , 'pentose and glucuronate interconversions' , 'glycerophospholipid metabolism' , 'amino sugar and nucleotide sugar metabolism' , 'cutin, suberin, and wax biosynthesis' , 'glycerolipid metabolism' , and 'alpha-Linolenic acid metabolism' (Fig. 2f). These GO terms and KEGG pathways of DEmRNAs may be involved in anther development in Chinese cabbage. lncRNAs can cis-regulate neighboring genes to regulate the expression of transcriptional or post-transcriptional genes 35,36 . lncRNA cis-regulated target genes were predicted based on positional relationships. We searched for lncRNAs and mRNAs with differential expressions 100 Kb upstream and downstream of the chromosomes; these lncRNAs were considered cis-regulated. By analyzing the cis-regulatory function of lncRNAs, a co-expression network of DElncRNAs and DEmRNAs was constructed. Most DEmRNAs and DElncRNAs were one-tomany matching, while several were one-to-one matching (Table S10). Considering that the graph cannot display much network regulation information between lncRNAs and mRNAs, more mRNAs for co-expression with lncRNAs were selected to establish their interaction network graph. As displayed in Fig. 3 GO and KEGG enrichment analyses of genes targeted by DElncRNAs were performed (Tables S7 and S8). Similar to DEmRNAs, common GO terms of DElncRNAs were 'nucleus, ' 'integral component of membrane, ' 'plasma membrane' , 'cytoplasm' , 'chloroplast' , 'protein binding' , 'regulation of transcription' , 'DNA-templated' , and ' ATP binding' (Fig. 2b). We analyzed the KEGG enrichment of the target genes of DElncRNAs and found 127 pathways, of which 16 were significantly enriched KEGG metabolic pathways, including 'photosynthesis-antenna proteins' , 'phenylpropanoid biosynthesis' , 'flavonoid biosynthesis' , 'plant-pathogen interaction' , 'photosynthesis' , 'protein export' , 'pentose and glucuronate interconversions' , and 'amino sugar and nucleotide sugar metabolism' , which may be involved in the regulation of anther development (Fig. 2c).
Identification of circRNAs and their functions. A total of 106 and 65 circRNAs were identified from the ' Ant' and 'Mix' RNA libraries, respectively (Table S3). We identified 165 DEcircRNAs, including 62 upregulated and 103 downregulated circRNAs (Table S5, Fig. 4a). circRNAs composed of exons or introns were counted    www.nature.com/scientificreports/ using the CIRCexplorer2 software. The results revealed that more than 97% of the circRNAs were full-exon circRNAs (Fig. 4b).
Constructing a ceRNA-miRNA-target gene regulatory network related to anther development in Chinese cabbage. It has been reported that lncRNAs and circRNAs can interact with miRNAs through miRNA reaction elements in the ceRNA network 27 . To reveal the global regulatory network of proteincoding RNAs and ncRNAs that are related to anther development based on gene annotation, the DEmRNAs that play an important role in anther development were selected to establish the ceRNA-miRNA-target gene regulatory network based on the ceRNA theory consisting of 49 DEmRNAs, 55 DEmiRNAs, 196 DElncRNAs, and 17 DEcircRNAs. We established candidate ceRNA relationships through targeting miRNA relationships and obtained 450 pairs of ceRNA relationships, including 97 DEmiRNA-DEmRNA, 281 DEmiRNA-DElncRNA, and 23 DEmiRNA-DEcircRNA interactions (Table S11). Perl script was used to integrate lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA networks, and the Cytoscape software (https:// cytos cape. org) was used to visualize the regulatory relationship (Fig. 5). In the ceRNA network, the target genes were mainly involved in plant hormone signal transduction, starch and sucrose metabolism, pentose and glucuronate interconversions, RNA degradation, plant-pathogen interactions, carbon metabolism, pyruvate metabolism, and carbon fixation in photosynthetic organisms, based on the KEGG analysis (Table S11). This advanced the understanding of the mechanism of anther development in Chinese cabbage. www.nature.com/scientificreports/ Some important transcription factors, enzymes, and regulatory proteins involved in anther and pollen development were identified in this study. These included AP2-like ethylene-responsive transcription factors (BraA01g001240.3C, regulated by 27 DElncRNAs, four DEcircRNAs, and two DEmiRNAs; BraA07g016770.3C, regulated by nine DElncRNAs and ten DEmiRNAs). It has been reported that the expression of many AP2 genes in anther deficient mutants also decreased, which may be involved in the regulation of anther development 38 . Other transcription factors, enzymes, and regulatiry proteins included ethylene-responsive transcription factors (BraA09g056770.3C, regulated by three DElncRNAs and one DEmiRNA; BraA08g016860.3C, regulated by two DElncRNAs and one DEmiRNA; BraA02g036660.3C, regulated by four DElncRNAs, two DEcircRNAs, and one DEmiRNA), which play an important role in the programmed cell death of papilla cells, compatible pollination, and reproductive growth in Brassica rapa 39 ; phosphoenolpyruvate carboxylase (PEPC) (BraA04g029950.3C, regulated by ten DElncRNAs and two DEmiRNA), which is involved in accelerating the accumulation of storage substances during pollen maturation 40 ; squamosa promoter-binding-like protein (SPL) (BraA04g003010.3C and BraA06g043820.3C, regulated by 15 DElncRNAs, one DEcircRNA, and 13 DEmiRNAs), which is a plant specific transcription factor involved in anther development 41 ; and nitrate reductase (BraA02g024390.3C, regulated by five DElncRNAs, one DEcircRNA, and one DEmiRNA), whose deficiency has an important effect on flower transformation and nitric oxide release during flower development in Arabidopsis thaliana 42 . The above were identified in the ceRNA network, indicating their specific synergistic regulation roles in anther development. More importantly, several well-studied miRNA-target gene pairs related to anther development were identified in this network, including miR156-SPL, miR159-AP2, and miR164-NAC domain-containing protein pairs.
The following important genes and their interactions with lncRNAs, circRNAs, and miRNAs in the main anther development process of Chinese cabbage were also identified in this study: BraA09g053910.3C, which affect anther dehiscence. Our results revealed that the regulatory network of ceRNA-miRNA-target genes may be involved in anther development in Chinese cabbage. Further functional studies of ceRNAs, miRNAs, and target genes will help us better understand the mechanisms of anther development and sexual reproduction.

Verification of RNA-Seq results.
To confirm the quality of the RNA-Seq and the expression patterns of miRNAs, lncRNAs, circRNAs, and mRNAs in the ' Ant' and 'Mix' groups, 13 DEmRNAs (eight known mRNAs and five novel mRNAs), seven DEmiRNAs (five known miRNAs and two novel miRNAs), seven DElncR-NAs, and seven DEcircRNAs were randomly selected for qRT-PCR (Fig. S2). Except for MSTRG.1758.1 and MSTRG.6184.1, the expression patterns of other genes were concordant with the RNA-Seq data, suggesting that the RNA-Seq results were reliable. Then, the expression patterns of six predominant anther expression genes (BraA09g005520.3C, BraA01g044980.3C, BraA06g036460.3C, BraA05g041890.3C, BraA04g030880.3C, and BraA02g023130.3C) in different parts of flower organs, including sepals, petals, pistils, anthers, and filaments, were analyzed. The expression levels of these six genes in the anthers were significantly higher than those in the other parts, indicating that these six genes may play an important role in anther development (Fig. S3).
Promoter activity analysis of predominant anther expression genes. To further explore the function of the identified genes that may be involved in Anther development and analyze their promoter activity, we selected six predominant anther expression genes, including four unknown functional genes (BraA09g005520.3C, BraA01g044980.3C, BraA06g036460.3C, and BraA04g030880.3C), one encoding polcalcin Bra n 2-like protein (BraA05g041890.3C), and one encoding GDSL esterase/lipase (BraA02g023130.3C). RNA-Seq data and qRT-PCR revealed that the expression of these genes in ' Ant' was much higher than that in 'Mix' (Fig. S2), and their expression levels in the anther were significantly higher than those in other floral organs (Fig. S3). Their detailed spatial expression patterns were examined in Arabidopsis plants expressing the promoter-driven β-glucuronidase gene. Histochemical staining demonstrated that GUS signals of the six genes only appeared in the anther (Fig. 6), suggesting that BraA09g005520.3C, BraA01g044980.3C, BraA06g036460.3C, BraA05g041890.3C, BraA04g030880.3C, and BraA02g023130.3C are anther specific genes and may perform important functions in regulating anther development.

Discussion
lncRNAs and circRNAs have been reported to act as ceRNAs and regulate each other by competitive binding to common miRNA response elements. The ceRNA regulatory network formed by lncRNAs, circRNAs, miRNAs, and mRNAs through miRNA response elements is of great significance to post-transcriptional gene regulation in various physiological and pathophysiological processes 43 . In the ceRNA network, miRNAs play an important role in the connection and regulation of different ceRNAs (lncRNAs and circRNAs). The mRNAs can be translated into proteins with direct functions, while lncRNAs and circRNAs can indirectly affect the expression of mRNAs by competitively binding to common miRNAs 27 . Anther development is an important biological process in plant reproduction, and its gene regulatory network has been gradually revealed. However, a comprehensive identification and analysis of lncRNAs and circRNAs as ceRNAs in the anthers of Chinese cabbage is unavailable. To investigate the function of lncRNAs and circRNAs in anther development in Chinese cabbage, we performed www.nature.com/scientificreports/ sequencing and whole-transcriptome analysis of lncRNAs, circRNAs, miRNAs, and mRNAs from the vegetative mass of four true leaves ('Mix') and anther (' Ant') samples. A total of 9055 mRNAs, 585 miRNAs, 1344 lncRNAs, and 165 circRNAs were identified as differentially expressed in the anthers. Compared with mRNAs, lncRNAs have fewer exons, lower expression levels, and shorter transcript lengths. To date, most of the functions of lncRNAs are not fully understood. Co-expression networks of lncRNAs and mRNAs aid in predicting the lncRNA function 44 . By analyzing the cis-regulatory function of lncRNAs, a coexpression network of lncRNAs and mRNAs was constructed. Most mRNAs and lncRNAs were one-to-many matched, but several were one-to-one matched, similar to those in cucumber 45 . It has been demonstrated that circRNAs play an important role in regulating miRNA-mediated post-transcriptional gene expression 46 . In this study, all 165 circRNAs identified were novel circRNAs, including 62 upregulated and 103 downregulated cir-cRNAs in the anther. Here, we first constructed a ceRNA-miRNA-target gene regulatory network related to anther development in Chinese cabbage based on ceRNA theory to obtain 450 pairs of ceRNA relationships. Our findings provide a systematic recognition of mRNAs and ncRNAs and lay the foundation for further work on the regulation mechanism of anther development and the utilization of hybrid breeding in Chinese cabbage.
Anthers are important male reproductive organs of plants that can produce pollen grains. Anther development is a complex process, involving a series of biological events 46 . During the normal development of microspores, many nutrients, including sugars, starch, amino acids, and proteins, play an important role in the development of anthers 47  The key regulatory genes in the process of anther development have been widely reported in the model plant Arabidopsis, including those involved in microsporogenesis, tapetum and callose layer development, pollen wall formation, and anther dehiscence 48 . In this study, we identified a series of genes involved in anther development and their interactions with miRNAs, lncRNAs, and circRNAs (Fig. 7). Studies have demonstrated that the expression of peptidyl-prolyl cis-trans isomerase (PPIase) can affect the pollen fertility of rice by affecting microspore  The tapetum is composed of a secretory cell, which is an important part of pollen formation. In the early stage of anther development, the tapetum surrounds the anther, providing various nutrients for microspore development. After meiosis of microspore cells, the tapetum secretes callose degrading enzymes, which decompose the callose wall of the cyst tetraspore, thus releasing microspores 50 . The F-box protein has been reported to play an important role in regulating tapetum degeneration and pollen maturation during Arabidopsis anther development 51 . Autophagy is an important cellular process that can break down cellular components during aging, starvation, and stress. In Arabidopsis, autophagy plays an important role in tapetum programmed cell death and pollen development 52 . Zinc finger proteins play an important role in the late development of the tapetum. Silencing of the tapetum-specific zinc finger gene (TAZ1) leads to premature degeneration and pollen abortion in Petunia 53 . Ca 2+ -binding proteins play an important role in regulating cell degradation and pollen formation in the rice tapetum. Ca 2+ -binding protein mutations cause delayed degeneration of tapetum cells, reduced degeneration of the callose wall around microspores, pollen abortion, and complete male sterility 54 9089.1, MSTRG.4677.1, MSTRG.4306.1, MSTRG.25939.1, MSTRG.21171,  MSTRG.19614.1, MSTRG.19265.1, MSTRG.18069.1, MSTRG.17974.1, MSTRG.14213.1, MSTRG.1403.1,  MSTRG.13821.1, MSTRG.12550.1, MSTRG.11661, MSTRG.11253.1, MSTRG.10845.1, and MSTRG.10406.1) and  circRNAs (circRNA161, circRNA146, and circRNA106). BraA05g031610.3C, encoding NAC domain-containing protein, is regulated by 21 DElncRNAs and ath-miR829-p5. Five other DEmRNAs encoding NAC domain-containing proteins, namely BraA10g029650.3C, BraA03g044390.3C, BraA02g043100.3C, BraA08g000170.3C, and BraA02g002460.3C, were the target genes of ath-miR164a, bra-miR164a, and mtr-miR164a, and further interacted with MSTRG.18058.1. These regulatory networks may participate in pollen wall formation in Chinese cabbage. Anther dehiscence determines the success of sexual reproduction in flowering plants by releasing pollen grains during pollination and fertilization. Many AP2 genes were downregulated in the anther deficient mutants, which may be related to the regulation of anther development. The putative AP2 interacting protein is involved in many functions of development and stress response, including flower development, pathogenesis, photomorphogenesis, and plant hormone signal transduction 38  www.nature.com/scientificreports/ and development. In Arabidopsis, the normal development of anthers depends on the interaction between cells so as to regulate cell proliferation and differentiation, and LRR-RLK and MAPK are involved in anther dehiscence and differentiation 57 . In rice, cellulose synthase-like (CSL) genes play a pivotal role in pollen formation, anther dehiscence, cell proliferation, stomatal formation, and cell arrangement in different tissues 58 . E3 ubiquitinprotein ligase affects the fertility of Arabidopsis by affecting anther dehiscence and pollen development 59  , and their related ceRNAs are possibly involved in anther dehiscence and differentiation. These results extend the regulatory network of anther development by adding new regulatory ceRNA-miRNA-target gene relationships. Here, all identified lncRNAs, circRNAs, miRNAs, target genes, and their interactions provide a meaningful database for future research. However, to understand the biological process of anther development, further experimental research and computational analyses are required.

Materials and methods
Plant materials. In this study, the Chinese cabbage double haploid (DH) line 'FT' , obtained from an isolated microspore culture, was used as the material 60 . The samples consisted of the following two parts. (1) On August 10, the DH 'FT' seeds germinated at room temperature were sown in a greenhouse of Shenyang Agricultural University after vernalization at 4 °C for 20 days. During flowering, the flower buds of different sizes were collected, the anthers were stripped out, and the anthers at different development stages were mixed and named "Ant. " (2) On August 10, the DH 'FT' seeds germinated at room temperature were sown in the greenhouse of Shenyang Agricultural University. When the seedlings grew to the four true-leaf stage, the entire plant was taken out and cleaned, and then the whole vegetative mass (including the root, stem, and leaf) was sampled and named "Mix. " The samples were quickly frozen in liquid nitrogen and completely ground for RNA extraction. For mRNA, lncRNA, and circRNA, the total RNA representing a specific fat type was used to consume ribosomal RNA according to the Epicenter Ribosomal RNA Zero Gold kit (Illumina). Purified poly (A)-or poly (A) + RNA fragments were cut into small pieces with divalent cations. Then, according to the mRNA-Seq sample preparation kit (Illumina), the RNA fragments were reverse transcribed to create the final cDNA library, with an average insertion size of 300 bp (± 50 bp) for the paired-end library. The paired ends were sequenced using an Illumina Hiseq 4000 (LC-BIO). miRNA identification and target gene prediction. After removing bases with adapter dimer and low complexity, common RNA families (tRNA, rRNA, sonRNA, and snRNA) and duplication from the raw reads, which were unique sequences of 18-25 nucleotides were mapped to specific species precursors in miRBase 21.0 to identify known and novel miRNAs. We analyzed the expression of miRNAs in ' Ant' and 'Mix' libraries and compared the two libraries in pairs to identify differentially expressed miRNAs (DEmiRNAs). When |log2(Fold change)|> 1 and P value ≤ 0.05, miRNA was considered to be differentially expressed between ' Ant' and 'Mix' libraries. To predict the miRNA-targeted genes and determine the miRNA binding sites, we used a target finder algorithm (Target Finder 50), and the alignment score value was used as the screening standard 61 . The miRNA and mRNA were compared, and complementary pairing in the comparitive position was scored 62 . The default maximum alignment score in this study was 4 points. GO terms and KEGG pathways of miRNA targets were also annotated. GO terms were enriched using Blast2GO by referring to the GO database 63 . KEGG pathway analysis was performed with reference to the KEGG pathway database (www. kegg. jp/ kegg/ kegg1. html) 64.
Identification and differential expression analysis of mRNA, lncRNA, and circRNA. The reads containing adapter contamination, undetermined bases, and low-quality bases were removed using Cutadapt 65 . FastQC was used to verify sequence quality (http:// www. bioin forma tics. babra ham. ac. uk/ proje cts/ fastqc/). The clean reads were mapped on to the Brassica rapa v3.0 reference genome according to Bowtie2 66 and Tophat2 67 , and the mapped reads were assembled using StringTie 68 . After the final transcriptome generation, the expression levels of all transcripts were estimated using StringTie and Ballgown 69 . Differentially expressed mRNAs (DEm-RNAs) were identified when |log2 (fold change)|> 2 and P value ≤ 0.01.
Transcripts that overlapped with known mRNAs and those less than 200 bp were discarded. CNCI 70 , CPC 71 , and Pfam 72 were used to predict transcripts with encoding potential. Transcripts with a CNCI score < 0 and a CPC score < − 1 were deleted. The rest of the class codes (I, j, o, u, x) were considered lncRNAs, which may play cis roles in adjacent target genes. To investigate the function of lncRNAs, the cis-target genes of lncRNAs were predicted. In this study, Perl script was used to screen 100 Kb upstream and downstream of coding genes. Then, we determined the functional analysis of lncRNA target genes using internal scripts. A standard of | log2 (fold change) |> 1 and P value ≤ 0.05 were used to screen differentially expressed lncRNAs (DElncRNAs).
A circexplorer 69,71 was used to assemble the mapped reads to circRNAs; following this, tophat-fusion and CIRCExplorer recognized the reverse splicing reads in the unmapped reads. Based on the structural characteristics and splicing sequence characteristics of circRNAs, we used the following criteria to identify circRNAs: (1) mismatch ≤ 2; (2) back-spliced junction reads ≥ 1; and (3)  www.nature.com/scientificreports/ expression of circRNAs was calculated using in-house scripts. The R package edgeR 70 was used to identify differentially expressed circRNAs (DEcircRNAs) between ' Ant' and 'Mix' libraries with the criteria of P-value ≤ 0.05 and | log2 (fold change) |> 1.
The GO terms and KEGG pathways of mRNAs, lncRNAs, and circRNAs were also annotated 63,64 .
Construction of ceRNA network. In this study, ceRNA analysis was divided into two parts: miRNA-mRNA and miRNA-lncRNA/circRNA. Target finder software was used for miRNA-mRNA analysis, and the alignment score value was used as the screening standard. The prediction results were scored by comparing the strands of miRNA and mRNA and complementary pairing in comparative positions was scored-one point for a mismatch or missing nt, 0.5 points for a G:U pairing, and double penalty point for a G:U pairing in the core area, starting from the 2 nd to the 13 th nt of miRNA in the core area. The alignment score indicated the degree of matching between the target gene and miRNA. The lower the score, the more complete and credible the match. The default maximum penalty score in this study was four points. Search 36 (36.3.6) was used for miRNA-lncRNA/ circRNA analysis with the following rules 73 : (1) bulges must be on the ncRNA and in the middle of the miRNA; (2) at most, four mismatches were allowed in positions other than the middle of the miRNA, and continuous mismatches could not exceed two; (3) bulges were not allowed in positions other than the middle of the miRNA. The software predicted the target binding relationship between miRNA and mRNA (3' UTR), lncRNA, and circRNA, and then constructed the ceRNA relationship pairs with circRNA/lncRNA-miRNA-mRNA based on common miRNA binding. The network diagram was constructed using Cytoscape (https:// cytos cape. org) to visualize the regulatory relationship.
Quantitative real-time PCR (qRT-PCR). The expression levels of selected DEmRNAs, DEmiRNAs, DEcircRNAs, and DElncRNAs were validated by qRT-PCR. Total RNAs and small RNAs of ' Ant' and 'Mix' were extracted using an RNApure Total RNA Kit (Aidlab, RN03, China). RNA was reverse transcribed using an Evo M-MLV RT Kit II (Accurate Biotechnology, AG11711, China). Oligo dT primers and random 6-mer primers were used to construct the 1st strand cDNA for quantitative analysis of mRNAs and lncRNAs. The miRNAs were reverse-transcribed using the downstream primers of U6 endogenous reference gene and the specific stem-loop primers. The circRNAs were reverse-transcribed using the qRT-PCR downstream primers and random 6-mer primers. qRT-PCR was performed using a SYBR® Green Premix Pro Taq HS qPCR Kit (Accurate Biotechnology) with an ABI 7300 instrument (Thermo Fisher Scientific, Waltham, MA, USA). The qRT-PCR procedure was as follows: 50 °C for 2 min; 95 °C for 30 s; 40 cycles at 95 °C for 5 s and 60 °C for 30 s; 95 °C for 15 s; 60 °C for 1 min; and 95 °C for 15 s. U6 was used as an endogenous control for the miRNAs. Actin was used as an endogenous control for circRNAs, lncRNAs, and mRNAs. All qRT-PCR reactions were performed in three technical replicates and three biological replicates. The primer sequences are listed in Table S12. RT is the specific stem-loop primers used to amplify miRNAs.
Promoter activity assay. The promoter region, which was 2000 bp upstream of the anther predominant expression gene, was amplified and a recombinant plasmid was obtained by connecting the target fragment with a pC1301IgT vector and transferred into Escherichia coli DH5 α competent cells. The expression vector was transformed into Agrobacterium tumefaciens GV3103 competent cells and was introduced into wild-type Arabidopsis according to the Agrobacterium-mediated floral-dip method 74 . The positive transgenic Arabidopsis was separately screened on MS solid medium containing 30 μg/mL hygromycin and then validated using PCR. Finally, GUS staining was performed on transgenic plants 75 .
Ethics approval and consent to participate. The current study complies with relevant institutional, national, and international guidelines and legislation for experimental research and field studies on plants (either cultivated or wild), including the collection of plant material.

Data availability
The datasets supporting the conclusions of this article are included within the article and its additional files. The transcriptome sequencing data were deposited in the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) Database under accession number GSE171705 and GSE171706. Genomic sequences and gene annotation information of Brassica rapa were downloaded online at http:// brass icadb. cn/#/.