Genome-wide identification and characterization of circRNAs in wheat tiller

This study found that the intergenic circRNAs of wheat were more abundant than those of other plants. More importantly, a circRNA-mediated network associated with tillering was constructed for the first time. Circular RNAs (circRNAs) are a class of endogenous non-coding RNAs with covalently closed circular structures, which play an important role in transcriptional and post-transcriptional regulation. Tiller is an important agronomic trait that determines plant morphological architecture and affects spike number in wheat. However, no studies on the characteristics and functions of circRNAs involved in the regulation of wheat tiller. Here, we performed a genome-wide identification of circRNAs using ribosomal-depleted RNA-seq from wheat tiller of two pairs near-isogenic lines. A total of 686 circRNAs were identified and distributed on 21 chromosomes of wheat, of which 537 were novel circRNAs. Unlike other plants, the majority of these circRNAs (61.8%) were derived from intergenic regions. One circRNA-mediated network associated with tillering was constructed through weighted gene co-expression network analysis, including 323 circRNAs, 117 miRNAs, and 968 mRNAs. GO and pathway enrichment analysis of mRNAs suggested that these circRNAs are involved in cell cycle, ncRNA export from nucleus, developmental process, plant hormone signal transduction, MAPK signaling pathway, RNA degradation. Of these circRNAs, ten circRNAs are associated with known tillering/branching genes in rice or Arabidopsis thaliana, including OsCesA7, EBR1, DTE1, CRD1, LPA1, PAY1, LRK1, OsNR2, OsCCA1, OsBZR1. In summary, we present the first study of the identification and characterization of circRNAs in wheat tiller, and the results suggest these circRNAs associated with tillering could play an important role in wheat tiller formation and development.


Introduction
CircRNAs are short of the 5′ cap and 3′ poly(A) tail; they can arise from exons, introns, and intergenic regions. Cir-cRNAs are a class of RNA molecules with a closed-loop structure formed by back splicing, are stable in expression, and can't be easily degraded. Recently, studies have shown that circRNAs are closely related to regulation of gene expression and disease occurrence Patop and Kadener 2018). CircRNAs were proved to be the role of microRNA decoys by preventing miRNA from binding corresponding mRNAs, which could provide ideas for later research (Hansen et al. 2013). Meanwhile, circRNAs can participate in biological processes such as ubiquitination and cell senescence and affect the function and expression of protein while interacting with it. But conversely, some proteins also play a role in the production and degradation of circRNAs (Huang et al. 2020). Although circRNAs have been widely reported in human diseases, there is little study about the biological functions of plant circRNAs and it's kind of difficult to identify the plant circRNAs due to the low abundance. CircRNAs and the homology of their host genes identified in Arabidopsis thaliana, rice, and maize demonstrate the widespread and importance of circRNA in plants (Ye et al. 2015;Zhang et al. 2019).
Wheat (Triticum aestivum. L) is one of the most important crops in the world, providing plenty of calories and protein for human diets (Consortium et al. 2018). Wheat tillers are special branches that arise from axillary buds at the basal nodes, determining spike number and plant architecture and affecting yield . Previous studies have shown that tillering is a complex quantitative trait regulated by transcriptional and post-transcriptional, and involved in multiple regulatory networks related to growth and development (Wang et al. 2016bWorland et al. 1998). For example, overexpression TaD27-B can significantly reduce the tiller number and enhance strigolactones biosynthesis in wheat . Furthermore, some non-coding RNAs have been confirmed to play an important role in post-transcriptional regulation of tillers in rice, containing miR156, miR393, and miR160 (Liu et al. 2017;Wang et al. 2021;Zhao et al. 2019). Circle is a widespread miRNA sponge in animals and plants, or a precursor of some genes, and still has great potential for tillering regulation. However, the function of circRNAs in tiller is still unclear.
In this study, we used high-throughput RNA-sequencing technology to identify circRNAs in wheat tiller, further examined the distribution of circRNAs on chromosomes and understand the features and possible functions in the regulation of wheat gene expression. The present study is the first report of genome-wide identification of circRNAs in wheat tiller, confirming the existence of circRNAs in wheat tiller and providing the basis for studying the potential function of circRNAs in wheat tillering.

Genetic materials
In this study, H461, CN16, and CM107 were used to develop near-isogenic lines. Among them, H461 is a low tillering wheat variety carry Qltn.sicau-2D; CN16 and CM107 are commercial cultivar with free tillering (Wang et al. 2016a, b). Two populations of recombinant inbred lines (RILs) are generated using single-seed descent, including (1) H461/CN16 249 F 5 lines, (2) H461/CM107 200 F 4 lines. Heterozygous lines were selected with the closely linked markers to Qltn.sicau-2D and then selfpollinated, and the process of heterozygote screening was repeated three times (Wang et al. 2016a(Wang et al. , b, 2019. In the last round of selection, two homozygous isolines were isolated, one with and the other without the low tillering genotype (Fig. 1). Finally, NIL7A and NIL7B were successfully developed from the H461/CN16 population (Fig. S1A); NIL11A and NIL11B were developed from the H461/CM107 papulation (Fig. S1B).
The two NIL pairs were planted in 2016 in a field at the Sichuan Agricultural University, Wenjiang District (30° 97′ N, 103° 81′ E), Sichuan Province, China, under natural conditions with conventional management. Each line consisted of 20 rows with a length of 1.5 m and a spacing of 30 cm between rows, while the sowing density was 15 seeds per row and the spacing was 10 cm between seeds.

RNA extraction, library construction, and RNA sequencing
Tillering nodes (Fig. S2) were separated and collected from 10 to 11 a.m. in the four stages (Z20-Z23) as described by Zhou (Zhou et al. 2020;Zadoks et al. 1974), immediately frozen in liquid nitrogen, and stored at − 80 °C until use. Each sample contained three independent plants at the same stage and was pooled together for RNA isolation with biological replicates.
Total RNA was extracted by using Total RNA Purification Kit according to the manufacturer's protocol (LC-Bio, CA), and DNA was removed using DNase I. The quantity of RNA was measured using the NanoDrop 2000 spectrophotometer (Thermo Scientific, USA). The Agilent 2100 Bioanalyzer (Agilent Technologies, USA) was used to measure RNA integrity number (RIN); only samples with RIN > 9.0 were used in the next step. Then, the total RNA was depleted rRNA using the Riobo-Zero Magnetic Kit (Epicentre, USA) and eliminated linear RNA with RNase R (Epicentre, USA). The remaining RNAs were used to generate circRNA-seq libraries as described by Zuo (Zuo et al. 2016). A total of 32 libraries were sequenced on Illumina Hiseq X Ten platform by BGI (Shenzhen, CA), and paired-end reads (PE150) were generated according to the standard Illumina protocol.
The identified circRNA should satisfy the following conditions: (1) circRNA identification was performed in each of the 32 samples using find_circ and circRNA_finder; the circRNAs with junction reads ≥ 2 were retained, (2) with canonical splice site signals (GT-AG), (3) for find_circ or circRNA_finder, circRNAs were identified in at least two samples, (4) finally, the intersection of circRNAs identified by the two software tools was formed.

CircRNA full-sequence assembly
The sequence of circRNAs was extracted with TBtools ) and repeated once as a hypothetical reference. All of the clean data aligned the hypothetical reference by hisat2 (Kim et al. 2015). Mapped reads were used for assembly of circRNA full sequence by StringTie (Pertea et al. 2015) and removed transcripts without back-splice junctions. Transcripts from different samples were integrated and removed the redundant data.

Polymerase chain reaction amplification
To validate circRNAs identified in this study, polymerase chain reaction (PCR) was performed using divergent primers and convergent primers. Genomic DNA was extracted from leaves of plants using the cetyl trimethyl ammonium bromide (CTAB) method as previously described (Wang et al. 2016b). cDNA was synthesized from the total RNA using the ABScript III RT Master Mix (ABclonal, CA). PCR amplification conditions were as follows: 94 °C for 5 min; 35 cycles of 94 °C for 30 s, 60 °C for 30 s, 72 °C for 30 s; 72 °C for 10 min. Then, Sanger sequencing was performed on PCR products. All the primers used in the present study are shown in Table S1.

Functional analysis circRNA
Classify circRNAs according to their position in genome by bedtools (Quinlan and Hall 2010) and the sequence of circRNA host genes were extracted by TBtools. To predict miRNA-circRNA interactions, 1063 wheat miRNAs downloaded from PmiREN (https:// www. pmiren. com/); psRobot (Wu et al. 2012) was used to predict the miRNA target sites of circRNA with default conditions. Based on ceRNA theory, circRNA-mRNA with the same miRNA target sites was identified.
The genes connected with circRNA acting as miRNA decoys and the circRNA host genes were consolidated into a set. A total of 1809 genes in the set were quantified by 102 Page 4 of 10 Kallisto using previously data (Zhou et al. 2020) and used to construct gene co-expression networks by weighted gene co-expression network analysis (WGCNA, Langfelder and Horvath 2008), with default parameters, except that network type was unsigned, R square was 0.8, and min module size was 30. The correlation coefficient (CC) between module and tillering was calculated, and the network of module associated tillering (CC > 0.7 or < − 0.7) was visualized by Gephi (Bastian et al. 2009). Gene ontology (GO) and KEGG pathways annotation using KOBAS v3.0 (http:// kobas. cbi. pku. edu. cn/).

Identification of circRNAs in wheat tiller
To identify the circRNAs in wheat tiller, 32 cDNA libraries were constructed between free tillering and low tillering wheat at the tillering stage. In total, 54,678,616-182,816,432 clean reads per library were obtained by removing lowquality reads, adapters, and poly-N (Table S2). 5877 circR-NAs and 990 circRNAs were identified by circRNA_finder and find_circ, respectively, with strict threshold of at least two unique back-spliced reads and identified in at least two libraries. A total of 686 circRNAs were identified in both software tools ( Fig. 2A). Among them, 149 circRNAs were already identified in the public database of plant circRNAs, and 537 circRNAs were novel identified (Fig. 2B). Genomic mapping revealed these circRNAs are distributed widely on 21 chromosomes of wheat, of which the A genome had the least number (177, 25.8%) and the B genome had the largest number (280, 40.8%) (Fig. 2C).

Conservation of the circRNAs
BLAST analysis was performed to evaluate circRNA conservation among different species in PlantcircBase and PlantCircNet; 142, 114, 53, and 7 identified circRNAs of wheat had similar sequence in rice, barley, maize, and Arabidopsis thaliana, respectively (Table S4); 75 of these circR-NAs were conserved across multiple plants species (Fig.  S3). For instance, Tae_circRNA00535 contained the gene TraesCS6B03G0595800 was found in three species (maize, rice and barley). Tae_circRNA00102, contained the gene TraesCS2A03G0556400 was found in three species (rice, barley, and Arabidopsis thaliana). These results suggest that circRNAs of wheat are conserved between other species, particularly in grass.

CircRNA-miRNA-mRNA network associated with tillering
CircRNAs are emerging as important regulators of physiological development by regulating its host gene expression, and the circRNA can act as miRNA sponges, weakening the inhibitory effect of miRNA on its targets. We found 305 of the 644 circRNAs derived from 275 host genes (Table S5), and 26 circRNAs could act as 117 miRNAs decoys and are associated with 1545 genes (Table S6). Expression levels of all the 1809 genes were calculated by Kallisto using the RNA-seq data from previous study (Zhou et al. 2020) and then were used for WGCNA analysis (Table S7).
All the genes were grouped into eight modules; the Mod-ule2 contained 968 genes and had significant correlation with trait 'tiller.' Module6 (109 genes) and Module1 (74 genes) had significant correlation with genetic background. No module was significant correlation with stage (Fig. 4A). GO enrichment and functional analysis of the module2 implies that module2 might be associated with multiple biological processes and molecular functions (Fig. 4B, C), such as ncRNA export from nucleus, miRNA transport, cell cycle, developmental process, flower development, protein binding, receptor activity, transporter activity, etc. Furthermore, KEGG enrichment analysis of module2 was performed; the genes in module2 were significantly enriched in pathways such as plant hormone signal transduction, ubiquitinmediated proteolysis, MAPK signaling pathway, starch and sucrose metabolism, RNA degradation.

Validation of the identified circRNAs
To confirm the identification of wheat circRNAs, ten circR-NAs (unique junction reads > 5) were randomly selected for experimental validation. Divergent primers and convergent primers were designed for each circRNA selected and then used to amplify both genomic DNA and cDNA. Of the 10 circRNAs, 8 were amplified successfully (Fig. 7A) and the back splice sites of Tae_circRNA00061 was also validated by Sanger sequencing (Fig. 7B).

Discussion
CircRNA is a novel endogenous and largely non-coding RNA, which functions as miRNA sponging, transcriptional modulation, and translating into protein. The consequences of these activities contain cell proliferation, altering the cell cycle process, which are essential roles in cellular processes and the development of plants (Cheng et al. 2018;Conn et al. 2017;Du et al. 2016;Zheng et al. 2016). Recent studies have shown that circRNAs are widely distributed in plants; however, there are few studies of circRNAs on plant branching/tillering. Therefore, revealing the characteristics of cir-cRNAs in wheat tiller is essential to further understanding non-coding regulators on branch/tiller in plants.

Characterization of circRNAs in wheat tiller
Various tools were developed to identify circRNAs, including circRNA_finder, find_circ, CIRCexplorer, MapSplice, CIRI (Memczak et al. 2013;Westholm et al. 2014). Combining two tools to identify circRNA would decrease the false positive rate and strengthen the output accuracy (Hansen et al. 2016). In this study, 686 wheat circRNAs were identified using both tools, circRNA_finder and find_circ, and 537 of them were novel circRNAs compared to the circRNA database. All the circRNAs could divide into three types, which the intergenic circRNAs account for the largest proportion. This result is different from previous studies in rice maize and Arabidopsis thaliana (Tang et al. 2018;Ye et al. 2015). This may be due to the tools that can identify circRNAs in the interregional independent of gene annotation and a relatively low gene density in wheat. Cereal crops circRNAs have a large number of isoforms, many circRNA isoforms are generated from the same locus by alternative circularization, and host genes with multiple exons are preferentially circularized (Ye et al. 2015;Zeng et al. 2018). Our experimental results also support the above views. Several studies reported that some circRNAs are conserved across species. For example, a study on rice and Arabidopsis thaliana, identified more than 300 circRNA, was similar in two species by comparing the homology of their host genes (Ye et al. 2015). Tang (Tang et al. 2018) reported 232 circRNAs were conserved between maize and rice; 122 circRNAs were conserved between maize and Arabidopsis thaliana. As expected, in this study, sequence conservation analysis indicated that the circRNAs of wheat tiller distributed in coding regions exist high similarity among rice, maize, and barley. However, the circRNAs are low similarity with Arabidopsis thaliana. These results suggested the sequence conservation of circRNAs in wheat tiller may exist widely across monocot, and the potential functions of these conserved circRNAs in monocot tiller development deserve further exploration.

Functions of circRNAs in wheat tillering
Many genes have been reported to be involved in branching/tillering of plants. In rice, overexpression of PAY1 or excessive expression LPA1 can alter tiller number by influencing the polar transport activity of auxin (Sun  Zhao et al. 2015). Nitrogen and boron are essentially required for plant growth; in Arabidopsis thaliana, T-DNA insertion lines of NIP5-1 showed decrease in the absorption of boric acid and increased sensitivity of shoot development to boron deficiency (Takano et al. 2006). OsNR2 and DTE1 can regulate tiller formation by affecting the uptake of nitrate and boron in rice, respectively (Gao et al. 2019;Liu et al. 2015). In addition, loss of function of EBR1 in ebr1 mutant leads to accumulation of OsBAG4 protein that inhibits rice growth and development, and tiller (You et al. 2016). OsCesA7 encodes cellulose synthase A catalytic subunit 7 has significantly impact on tiller number (Wang et al. 2016a). In the present study, to further reveal the functional importance of cir-cRNAs, 323 circRNAs, 117 miRNAs, and 968 genes were used to construct circRNA-miRNA-mRNA and circRNA-mRNA networks based on WGCNA analysis to study their functions (Fig. 5A). GO analysis was enriched in ncRNA export from nucleus, cell cycle, developmental process. And several KEGG pathways are involved in plant hormone signal transduction (auxin), MAPK signaling pathway, suggesting that circRNAs of wheat tiller may involve cell proliferation and tiller formation; 12 wheat genes of the network were homologous to known tillering/branching genes in rice or Arabidopsis thaliana, containing OsC-esA7, EBR1, DTE1, CRD1, LPA1, PAY1, LRK1, OsNR2, OsCCA1, OsBZR1. Taken together, we speculate these circRNAs could play an important role in wheat tiller formation and development.

Conclusion
In this study, we identified and characterized 686 circRNAs in wheat tiller containing 537 novel circRNAs through analyzing 32 transcriptome data. One Major module associated with wheat tiller was found by WGCNA analysis, containing 323 circRNAs, 117 miRNAs, and 968 genes. Furthermore, circRNA-miRNA-mRNA and circRNA-mRNA networks were constructed; the possible functions of these circRNAs were predicted using GO and KEGG analysis. In addition, the regulatory mechanisms and functions of these circR-NAs in wheat tiller require future experimental studies to elucidate.