Circular RNA expression profiles during the differentiation of mouse neural stem cells

Background Circular RNAs (circRNAs) have recently been found to be expressed in human brain tissue, and many lines ofevidence indicate that circRNAs play regulatory roles in neurodevelopment. Proliferation and differentiation of neural stem cells (NSCs) are critical parts during development of central nervous system (CNS).To date, there have been no reports ofcircRNA expression profiles during the differentiation of mouse NSCs. We hypothesizethat circRNAs mayregulate gene expression in the proliferation anddifferentiation of NSCs. Results In this study, we obtained NSCs from the wild-type C57BL/6 J mouse fetal cerebral cortex. We extracted total RNA from NSCs in different differentiation stagesand then performed RNA-seq. By analyzing the RNA-Seq data, we found 37circRNAs and 4182 mRNAs differentially expressedduringthe NSC differentiation. Gene Ontology (GO) enrichment analysis of thecognate linear genes of these circRNAsrevealed that some enriched GO terms were related to neural activity. Furthermore, we performed a co-expression network analysis of these differentially expressed circRNAs and mRNAs. The result suggested a stronger GO enrichmentin neural features for both the cognate linear genes of circRNAs and differentially expressed mRNAs. Conclusion We performed the first circRNA investigation during the differentiation of mouse NSCs. Wefound that12 circRNAs might have regulatory roles duringthe NSC differentiation, indicating that circRNAs might be modulated during NSC differentiation.Our network analysis suggested the possible complex circRNA-mRNA mechanisms during differentiation, and future experimental workis need to validate these possible mechanisms. Electronic supplementary material The online version of this article (10.1186/s12918-018-0651-1) contains supplementary material, which is available to authorized users.


Background
Canonical RNA had been known asa linear molecule that terminates in 5′ cap and 3′ poly(A) tail. However, in 1976 a new class of endogenous circular RNA (circRNA) was discovered in a viroid [1]. Since the first discovery, many cir-cRNAshave been identified and they present in all eukaryotes ranging from yeast, fruit flies (Drosophila), to humans [2,3]. As a class of non-coding RNAs, circRNAs are characterized by covalently closed continuous loop. Due to the 5′ and 3′ end deficiency, circRNAs cannot be degraded by RNase R and present permanent stability in cells [4]. For many years after 1976, there had been only a few circRNAs whose expression was experimentally validated in cells [5]. However, the recent application of innovative, high-throughput sequencing technologies allowed investigators to discover tens of thousands of circRNAs from poly (A)-minus or RNase R treated RNA-Seq data [6,7]. In the past a few years, circRNAs increasingly attractedresearchers' attention and re-emerged as a burgeoning class of non-coding RNAs [8].
With the accelerating generation of RNA-Seq data, more and more cellular and functional circular RNAshave been recently reported in different cell types and tissues [9]. In human neural tissue, circRNAs were foundto play a functional role in maintenance, plasticity and abnormality of neural circuitry [10,11]. For instance, Yang and his colleagues reported that circZNF292 silencing could suppress the proliferative and angiogenic potential of glioma cells, suggesting that circRNAsmight regulate many biological processes [12]. Recently, many lines of evidence have unveiled that circRNAs are enriched in brain tissuewhen compared with their expressions in other tissues [13]. Investigators also found that while most circRNAs are inlow expression level, they are conserved and expressed in a time-, cell type-andgene-specific manner [14]. Cognate genesrefer to those linear genes have the same sequences with circRNAs. Some circRNAs are present at higher copynumbers when compared to their cognate genes (> 10-fold) [15].So far, the function of most cir-cRNAshas still been elusive. One common function is that some circRNAs can act their roles as microRNA (miRNA) 'sponges'sponge RNAs containing complementary binding sites to a miRNA of interest so that theyprevent miRNAs from the connections with the target mRNA [16,17]. Some circRNAs can also regulate miRNAs through RNA-binding proteins (RBPs) and affect cellular function, indicating another important function in the interaction with miRNAs [18]. Although functional studies have been reported in literature, investigation ofcircRNA'sfunction in neural tissue has not been systematically done yet.
Central nervous system (CNS) is a significant part of neurodevelopment, which contains brain and spinal cord. The CNS development starts in early embryo and plays an extremely important role in human development from fetal to adult stages. In recent years, studies showed that circRNAs have close connection with neuronal development. CircRNAs are not expressed equally in the CNS, rather, they are differentially distributed in different brain regions and different developmental stages [19]. For example, Vanet al. found a striking regulation of circRNAs during neuronal development [20,21]. In the meanwhile, several circRNAswere found to be expressed at high level in the cerebellum when compared withthe brainstem [17]. A large amount ofcircRNAswere detected in synapses and their regulations were also confirmed [22,23]. Many cir-cRNAswereexperimentally verified in different cell typesin mice [24]. Zhao et al. reported that circRNAsmight have physiological functions in neuron, andsome of these cir-cRNAswere related tobrain disease [25]. Many brain diseases are attributed to the abnormal development of CNS [26]. For example, dysplasia in CNS leads to many cerebral diseases, such as CNS tumor in childhood, and Alzheimer's and Parkinson's in elder people [27,28]. Some studies indicated that over-expressed circRNAsmight initiate the occurrence of neurological diseases [29]. However, the expression pattern and function of circRNAs in the brain diseases or CNS development are still largely unclear.
Recently, several circRNAs were found to involve in neural function.For example, twocircRNAs, circRNAsexdetermining region Y(SRY) and cerebellar degenerationrelated protein 1 antisense(CDR1as) have several binding sites with miR-138 and miR-7, respectively [30,31]. By binding miRNA, these circRNAscould regulate the expression of miRNA and furthermore suppress their function, which is known as 'sponging RNA' [32]. And circ-HIPK3 was shown to bind with miR-24 and furthermore modulate human cell growth [33]. While these cir-cRNAs act as competitive RNAs for miRNAs, other studies suggested their connection with RBPs as well. CircRNAs from the muscle blind (mbl) and FOXO3 genes were reported to bind,sequester andtransport RBPs [34]. The researchers thought these circRNAs might regulate the interaction of RBPs with their RNAtargets [35]. Some other studies revealed that alternativesplicing of circRNAs might also lead to the newbinding sequences for some RBPs and,thus, influence functions [36]. So far, it is not clear whether miRNA sponging andRBP binding are shared functions of circRNAs.
In this work, we examined circRNA expression pattern during the NSC differentiation by analyzing the RNA-Seq data. Previous experimental results reported that at the time course of differentiation day 2, the NSC activation reached their peaks and the NSC differentiation was also strongly activated [37,38]. In this study, we identified differentially expressed circRNAs and their cognate linear mRNAs in different differentiation stages.There wasa total of 37 circRNAs differentially expressed in this process. It is likely some of thesecircR-NAsinvolved in differentiation and resulted in the corresponding expression profiles. Further analysis of the expression profiles during NSCdifferentiation could help us uncover the possible regulatory circRNAs. We further found out the differentially expressed mRNAs and then constructed a co-expression network between these potentially regulatory circRNAs and differentially expressed mRNAs by the same binding miRNA. The Gene Ontology (GO) enrichment analysis on them suggested stronger enrichment in neural features and pointed outa possible regulation of circRNAsduring the differentiation.Furthermore, the opposite expression patterns between circRNAs and mRNAs suggested complex circRNA-mRNA mechanisms in the NSC differentiation.

CircRNAsare abundant and highly expressed in NSC differentiation
Total RNA was collected from mouse NSCs cultured in differentiation-suppressed medium or induced to differentiation with two replicates (Fig. 1). One group of NSCs wasculturedand kept undifferentiated with the different iation-suppress ingredient bFGFin 6 days as 0d.nsc group. In the 2d.nsc group, NSCs were first kept undifferentiated in 4 days and then induced to differentiation in 2 days without adding bFGF.And in the 6d.nsc group, NSCs were induced at the beginning of culturing and were kept in differentiation state in 6 days. Paired-end ribominus RNA sequencing (RNA-Seq) was performed, and the UROBORUS computational pipeline was applied to detect potential circRNAs in differentiation [39]. Firstly, RNA-Seq data were mapped to reference genomeusing toolTopHat (version 2.1.1). The transcriptome was then reconstituted with Cufflinks(version 2.2.1) (Additional file 1: Table S1). Welisted the top 200 circular RNA genes expressed in eachgroup, and presented their pair-wise expressioncorrelationmatrix resultsofthese 200circular RNAgenes (shown in Additional file 1: Figure S1). The correlation values ofthese three groupswere allabove 0.94, indicating the strong consistencyand reproducibility.Moreover, the moderate correlations (0.68-0.84)in2d.nsc versus 0d.nsc and 6d.nsc versus 2d.nscindicated their close relationship but withdifferent expressionpatterns. In the comparison of the 6d.nsc and 0d.nsc groups, the correlations were among 0.51-0.57. The weaker correlations suggested strongerdifferent #iation-related circular RNA expression variation between groups6d.nsc and 0d.nsc than that between groups 2d.nsc and 0d.nsc, or between groups6d.nsc and 2d.nsc. Differential linear RNAexpression analysis was further performed to check the quality of the RNA-Seq data in NSCdifferentiation experiments. The comparison between groups2d.nsc and0d.nscidentified2917linear RNAs weresignificantly upregulated or downregulated in the NSC differentiation process, and2360differentially expressedlinear RNAsbetween groups6d.nsc and 2d.nsc (Additional file 1: Table S2 and Figure S2). We performed GO enrichment analysis onthese linear RNAs. The resultsindicated the consistent and significant enrichment of these genes in function that is related to neurodevelopment, neural activity, and differentiation (Additional file 1: Figure S3). Moreover, the expression of well-known specific markers of neurodevelopment, neural activity or differentiation was also found in the differentially expressedgene list, furtherhighlighting the regulatory function of these differentially expressed mRNA during the differentiation.
To further find the specific mRNAs with regulatory function during the differentiation, we found1173 mRNAs that weredifferentially expressed in both differential expression analysis between2d.nsc versus 0d.nsc and 6d.nsc versus 2d.nsc.That is, these mRNAs were differentially expressed during the whole NSC differentiation. Then, we separatedthese differentially expressedmRNAs into four expressionpatterns for better illustration in further analysis. Pattern 1 consisted ofthose mRNAs with their expression upregulated;we named this expression profile as0d.nsc < 2d.nsc and 2d.nsc < 6d.nscduringthe differentiation. Pattern 2 contained those mRNAs with the opposite expression against pattern 1.Pattern 3 consisted of mRNAs with expression profileas0d.nsc < 2d.nsc and 2d.nsc > 6d.nsc. And similarly, pattern 4 showed the opposite trend ofPattern 3. We separately applied GO and Kyoto Encyclopedia of Genes and Genomes (KEGG)pathway enrichment analyses on these differentially expressed mRNAsduring the differentiation ( Fig. 2 and Additional file 1: Figure S4). As summarized in Additional file 1: Table S3, these analyses consistently revealed the enriched GO terms much related to neural activitywith one exception (the results between Pattern 4 with other Patterns). Patterns 1, 2 and 3 had several irrelevant GO enrichment terms in their top lists while the top GO enrichment terms in Pattern 4 were most connected to the neural activity or development of CNS. This observation suggested thatmore mRNAs with lower expression in 2d.nsc had regulatory function in differentiation. When we traced back to the number of expressed mRNAs in the different groups, we found that the number of mRNAs is the lowest when the threshold was set to FPKM> 10 or FPKM > 100 (Additional file 1: Table S1). In the previous GO enrichment analysis, we also discovered that theenriched GO terms from down-regulated mRNAs in 2d.nsc versus 0d.nsc and up-regulated mRNA in 6d.nsc versus 2d.nscweremore associated with neural activity than those up-regulated mRNAs in 2d.nsc versus 0d.nsc and down-regulated mRNA in 6d.nsc versus 2d.nsc, which coincided with the variation in Pattern 4 (Additional file 1: Figure S3). This phenomenon suggested that there were fewer mRNAs expressed in the medium differentiation stage, and those mRNAs in Pattern 4(i.e., expression profile in 0d.nsc > 2d.nscand 2d.nsc < 6d.nsc) might more involve in neuronal related functions during the differentiation of NSCs.
In the next step, RNA-Seqshort reads were aligned to the reference genome (version GRCm38/mm10) and the UROBORUS computational pipeline was applied to detect circRNAs in differentiation. More than 4000 potential circRNAs were found in every group (Additional file 1: Table S4). We listed the top 100 expressed cir-cRNAs and found that all these circRNAswere expressed in tissues from mouse brain based on our search of the database circPedia. This additional information preliminarily confirmed the expression feature of circRNAs detected in the NSC differentiationin our study. We further checked the function of their cognate genesusing UCSC (The University of California Santa Cruz) Genome Browser. We found many of them have regulations in development and neural activity. In the next, we further explored whether these circRNAshave regulatory function during NSC differentiation.

CircRNAexpression is modulated during NSC differentiation
To perform quantitative analysis ofcircRNA expression profiles, we firstly filtered out those circRNAs with expression level < 0.1 RPM. Then, we excluded those cir-cRNAs with only expressed in onereplicate. The process resulted in80 circRNAsfor further studies (Fig. 3).All the replicates had an overall better correlation for both cir-cRNAsandtheir cognate mRNAs.CircRNAs had overall correlation coefficients above 0.93 in replicates. And the expressions between these circRNAs and their cognate genes in replicates were compared and the correlation coefficientswere all above 0.92 (Additional file 1: Figure S5). We tested the potential functions of these circRNAs by gene set enrichment analysis. We first listed the cognate genes of these circRNAs and then applied GO enrichment analysis of them (Fig. 4).
For GO analysis, genes were organized into hierarchical categoriesto discover their regulatory networks based onthree domains: Biological Process, Cellular Component and MolecularFunction. The results showed thatmost of the significant enriched GO terms wererelated to CNS development.We also found that these cognate genes were associated with the following GO Biological Process terms: regulation of synapse assembly, axo-dendritic transport, receptor localization to synapse, cerebellar granule cell differentiation and neuron remodeling; GO Cellular Component terms: dendritic shaft, dendritic spine, synapse, cell projection, cell junction, neuronal cell body and axon; and GO Molecular Function terms: RNA binding and protein binding. These significantly enriched GO terms are closely associated with CNS development, suggesting the potential connections between these circRNAs and CNS development. Taken together, these results suggested that the circRNAs we identified might have close connection with the neural activity.
The comparison of circRNA expression profiles between 2d.nsc and 0d.nsc and between 6d.nsc and 2d.nsc revealed a global change in circRNA expression, that is, these molecules were likely altered during NSC differentiation (Tables 1 and 2). We identified 28 and 22 cir-cRNAs that were differentially expressed between 2d.nsc and 0d.nsc and between 6d.nsc and 2d.nsc, respectively. We found that 16 circRNAs were upregulated and 12 circRNAs were downregulated in comparison between 2d.nsc versus 0d.nsc. And in the comparison between 6d.nsc and 2d.nsc, 21 circRNAs were downregulated and only 1 circRNA was upregulated. GO enrichment analysis was applied to the cognate genes of these circRNAs: 21 downregulated circRNAsin6d.nsc versus 2d.nsc, 16 upregulated and 12 downregulated circRNAsin2d.nsc versus 0d.nsc (Fig. 5a, b and c). As the results, both downregulated circRNAs in 6d.nsc versus 2d.nsc and upregulated circRNAsin 2d.nscversus 0d.nschad the similar significant enrichment biological process terms, which were also strongly related to NSC differentiation. The GO analysis on downregulated circRNAsin 2d.nscversus 0d.nscindicated the significant enrichment for gene differentiation and translationwhich were alsoenriched in NSC differentiation. Collectively, the resultsabove suggested that circRNA expression is broadly modulated during NSC differentiation and these altered circRNAs are related to neural function.

Possible regulation of circRNAsin NSC differentiation
To find the possible regulatory circRNAsduring NSC differentiation, we examinedthe circRNAdifferential expression changesbetween conditions 2d.nsc and 0d.nsc, as well as,between conditions 6d.nsc and 2d.nsc. The results showed that11 circRNAs were upregulated in group 2d.nsc versus 0d.nsc and downregulated in group 6d.nsc versus 2d.nsc while only 1 circRNA had the opposite expression profile, i.e., expression profile in 0d.nsc > 2d.nsc and 2d.nsc < 6d.nsc (Table 3). Interestingly, our GO enrichment analysis on the cognate genes from these 11 circRNAs were consist with the previous GO results (Fig. 5d). These 11circRNAsmay perform function related with NSC differentiation. In our search of thecircPedia database, we found all these circRNAs expressed in tissues of mouse brain. We noted thatonly three circRNAs,Circ-APP,Circ-Arhgap5and Circ-Bnc2, shared the same expression patterns with their cognate linear genes.The remainingcircRNAsmatched only Pattern 4 (expression profile in 0d.nsc < 2d.nscand 2d.nsc > 6d.nsc) and their cognate linear genes were only withPattern 1 (expression profile in 0d.nsc < 2d.nscand 2d.nsc < 6d.nsc).
We next investigated these 11 circRNAsfor their cognate linear RNAs' functionfrom previous reports.It has been reported that gene App regulates neuronal stem cell differentiation andinduces glial differentiation [40]. The APPswe/PS1DE9 (APP/PS1) mice have been widely used for the study of Alzheimer's disease due to its lack of App protein that contributes to Alzheimer's disease [41]. Arhgap5 has no direct function on brain or other neural activity, but it is critical for embryonic mammary bud development [42]. Previous studies have shown that there is significant expression ofPrrc2b protein in developing rat brains [43]. Pak et al. found that Nrxn1expression is high during human neocortical development and aging process. Heterozygous mutations in Nrxn1can lead to synaptic transmission defects [44]. Kim et al. recently identified that Glis3could regulate neurogenin-3 expression in pancreatic β-cells [45]. Acosta et al. applied experiment and confirmed thatAdgrl3 variants are associated with a refined phenotype of attention-deficit/hyperactivity disorder(ADHD) in the MTA study [46]. An ultraconserved brain-specific enhancer in ADGRL3were reported with-ADHD susceptibility [47].STK39has been reported as associated with Parkinson's Disease and has regulation in cell invasion and proliferation [48]. Besides, several groups reported that Dcbld2 is one identification of novel neuroendocrine-specific tumor genes and regulates gastric cancer cell proliferation and invasion [49].Collectively, these previous studies supported our analysis results that circRNAs can performs similar function because the circular and linear RNA were homologous. There has been no direct evidence ofGnptg,Ralgps1 and Bnc2in CNS development or involvement in other neural activity, but our resultsindicated that the circular isoforms from these three genes have potential function during NSCs differentiation. This novel finding warrants further investigation in future and we will keep studying on them.

Construction of the circRNA-mRNA co-expression network
To further detect the regulation of these circRNAs in NSC differentiation, we constructed a co-expression network for thedifferentially expressed circRNA and mRNA molecules. We first downloaded the database on mirTarBase and found the miRNA target sites verified by experiment for the differentially expressedmRNA. For the lack of the validatedfunction of circRNAs, we could not find their biological function such as their targeted miRNA or co-expressed mRNA in the existing database. To find their potential connection, miRanda was applied to detected the differentially expressedcircRNAs for their possible miRNA target sites for further study. According to the knowledge, a circRNAhas a junctionsite between the head and tail of its sequence. We need to consider the possibility of binding site on the junction site. Foridentifying possible miRNA target sites for circRNAs, we traced back to the sequences of their cognate genes. Then, we copied 100-bp reads from the start site and added them to the end of the sequence to assemble a pseudo circRNA sequence because of the cyclic further of circRNAs.Next, we appliedmiRanda for predicting cir-cRNAbinding sites with miRNA. The results were filtered by setting stricter threshold with miRanda score > 150 and minimum free energy < − 20. By this analysis pipeline, we obtained 492miRNA target sites for these circRNAs. Based on the above results, we constructed a circRNA-mRNA co-expression network (Fig. 6). We listed the mRNA which shared the same miRNA target sites with circRNA. Then, we built the co-expression network by using Cytoscape.
To further find the possible regulatory roles of the RNA (both circRNA and mRNA), we first applied GO enrichment analysis on all the mRNA genes in the reconstructed co-expression network (Additional file 1: Figure S6). We found that these mRNA genes were enriched in cell adhesion, multicellular organism development, long-term synaptic potentiation, regulation of cell proliferation, regulation of osteoblast proliferation, long-term memory, neuron migration, neuromuscular junction development, cell differentiation and dozens of GO Biological Process terms closely associated with NSC differentiation.And these mRNA genes were also enriched in axon, cell surface, neuronal cell body, dendrite, cell junction, neuromuscular junction, neuron projection and cell-cell junction as GO Cellular Component terms. These significant enriched GO terms are associated with NSC differentiation, suggesting mRNAs in the co-expression network involved in the regulation during NSC differentiation. As for the KEGG Pathway analysis, Axon guidance(mmu04360) was most significantly enriched. These results verified the function of mRNA in the co-expression network, which further indicates the regulatory function of circRNAs in neural activity. However, to better study on these circRNAs, we need to classify the mRNA in the network for their different expressed patterns during the differentiation.
These GO termsare much related to neural activity. It is worth noting that theseenriched GO termsin Pattern 4ranked the topin the analysis rather than some irrelevant terms appeared in the top list for other patterns. The result implies that the mRNAs with expression profile in0d.nsc > 2d.nscand 2d.nsc < 6d.nscmight regulate the differentiation, which coincided with the phenomenon that group 2d.nsc had the least mRNAs at 0.1, 10 and 100 thresholds (Additional file 1: Table S1).
In the GO enrichment analysis and KEGG pathway analysis, mRNAs expressed in group 2d.nsc had better results than the other two groups;again, this isconsistentwith the above phenomenon.These results further validated that mRNAs in the Pattern 4 might have regulatory function in differentiation and their co-expression network with circRNAsmight be critical in the process as well. In the NSCs, circRNAs had an overall high expression in the medium differentiation stage while mRNA had the opposite trend.Their specific function duringdifferentiation will need be experimentally investigated in future. Interestingly,for the Pattern 1 (expression in 0d.nsc < 2d.nsc < 6d.nsc), the top enriched GO termsin Biological Process included negative regulation of ossification, negative regulation of cell proliferation, negative regulation of mitotic cell cycle, negative regulation of osteoblast proliferation,central nervous system developmentpositive regulation of axonogenesis and cell differentiation. These enriched GO terms suggestedmRNA genes in Pattern 1 had strong connection with differentiation and CNS development. They have regulatoryroles in suppressing proliferation and promoting differentiation of NSCs. However, this resultdid not stand out from the differentially expressed mRNAs of the previous Pattern 1 (Additional file 1: Figure S4). These GO terms did appear in the top listbut the bottom from enriched Biological Process terms.In Pattern 4, we observed that most enriched GO terms in Biological Process had a better performance with stronger connection with differentiation. Taken together,circRNAs and mRNA in our co-expression network had closer relationships with NSCs differentiation and showed potential regulatory functions. In future work, we will generate more experimental datasetto gain higher confidence in the expression patterns of thesecircRNAs and the linear mRNA co-expression network.
Due to the limitation of technology and funding problems, we did not apply further biological experiments to validate the results. The judgement of some GO enriched analyses was set as P-value rather than adjusted P-value because the number of RNAs is small. Besides, we also need perform more replicates for high statistical confidence in the further study.

Conclusions
In this study, we first found that 12 circRNA expression was modulated during mouse NSC differentiation. Then, weconstructed one specific co-expression network for the differentially expressed circRNAs and mRNAs. The GO enrichment analysisrevealed that circRNAs mightplay regulatory roles during NSCs differentiation. In addition, we found that the expression profiles of circRNAswere different from mRNAs. There were more circRNA in 2d.nsc than in 0d.nsc and 6d.nsc while there were fewer mRNAs in 2d.nsc than in 0d.nsc and 6d.nsc. In addition, circRNA in Pattern 3 had closer connections with neural activity while mRNA in the opposite Pattern 4 had closer relationships. Based on this observation, we hypothesized that there are complexcircRNA-mRNA regulation mechanisms that act on NSC differentiation.This study is the first to investigate the circRNA expression profiles during the NSCs differentiation in mice and we will keep studying on this.

NSC differentiation and RNA extraction
To identify the biological function of circRNAsduring the differentiation of NSCs, we first separated and cultured NSCs from the wild-type C57BL/6 J mouse fetal cerebral cortex. NSCs from mouse fetal cerebral cortex at E15.5 days were then dissociated by trituration. After centrifugation at 1000 r/min for 4 min, cells were cultured and seeded in culture flask treated with regent Poly-L-ornithine hydrochloride (GIBCO) and Fibronectin (GIBCO). This   procedure could help NSCs adherent to the surface of the flask. Cells were maintained at 37°C in the incubator for 2 weeks, during which most alive cells in the flask were NSCs with the help of the specific medium selection. The recipe of proliferation medium is DMEM/ F-12(1:1) (GIBCO) with 2% B27 (GIBCO), 1% Penicillin-Streptomycin (GIBCO) and 25 ng/ml bFGF (Mouse basic fibroblast growth factor, GIBCO). Next, approximately 2 × 10 6 cells/mL were seeded in six separate flasks. We divided these flasks into 3 groups. Two flasks of NSCs were kept their undifferentiated state during a six-dayculturing. They did not differentiate and we named this group as group 0d.nsc. Another replicate of flasks was induced to differentiate on day 4. We induced them to differentiate into neuron and neurogliocyte cells by wiping off the differentiation-suppress regent bFGF. These NSCs differentiated in 2 days and we marked these as group 2d.nsc. The rest of the flasks differentiated at the first time of the experiment and kept differentiating in 6 days. We named this group as 6d.nsc. We extracted the total RNA of these six flasks on day 6. Finally, we extracted total RNA from the differentiate cells at three different time points. Total RNA was extracted by using Total RNA Kit (TIANGEN, co. LTD., China), which used RNase-free spin column to isolate RNA from the cells. These RNAs were stored in liquid nitrogen and sent to Beijing Genomics Institute (BGI) for transcriptomic assay (RNA-Seq, paired end, 100-bp). A total of six groups of data were obtained andthere were approximately 43 million reads for each sample.

Library preparation and sequencing
End-repair and adaptor ligation were processed by using Illumina's TrueSeq Total RNA Library Prep Kit. Size of 250-300 bp cDNA fragments were separated and then PCR-amplified for about 20 cycles to build the library. After further purification and detection, the RNA-Seq was generated by using Illumina Hi-Seq 4000 sequencer.

Identify differential expression genes
RNA-Seqshort reads were aligned to the mouse reference genome (GRCm38 mm10) using TopHat (version 2.1.1). On average, approximately 85.7% (minimum 82.3%; maximum 88.9%) of RNA sequencing reads could successfully mapped to the mouse reference genome. By using Cufflinks, we assembledthe transcripts, which had been annotated in database GENCODE. Then we employed the tool Cuffdiff (version 2.2.1) to identify the differential expression genes and transcripts between the different time courses. Stricter standard Q-values were used and the cutoffs were set as below 0.05.

CircRNA detection pipeline
The RNA-Seq data were initially filtered through TopHat by mapping to the reference genome. The canonical splicing sites were detected and most reads were mapped. Then, the unmapped reads, were processed by URO-BORUS (version 1.0.0, http://uroborus.openbioinformatics.org/en/latest/).

Functional pathway analysis
GO terms have been commonly used forenrichment analysis of genes of interest (e.g.,differentially expressed genes in this study). We performed functional annotation clustering by using the Database for Annotation, Visualization and IntegratedDiscovery (DAVID version 6.8) [50]. The background gene set consists of all expressed in mouse muscles andthe P-value was set< 0.05. We also performed enrichment analysis using KEGGpathways and the P-value was set < 0.05 [51].
Prediction of miRNA target sites on circRNAs By using miRanda, the database for predicted microRNA targets and target downregulation scores, we obtained the miRNA target sites for circRNAs. We assembled the sequences of circRNAs by following two steps: (1) assembled the annotated exons from the differentially expressedcircRNAs, and (2) copied 100-bp reads from the startof the sequence and added them to the end of the sequence. Then, we applied miRanda on these sequences to predict miRNA target sites.We set the complementarity score threshold to 150.0 and set the complementarity energy threshold to − 20.0 kcal/mol to improve the prediction accuracy.
Visualization of the circRNA-mRNA co-expression network By combing the verified miRNA targets from mirTar-Base and the predicted miRNA targets from miRanda, we constructed the circRNA-mRNA co-expression networkby using Cytoscape (version 3.5.1), an open source software platform for visualizing complex networks. The network contains three types of nodes (circRNA, miRNA, mRNA), and the nodes mRNA and circRNAshared the same miRNA targets.The construction of the co-expression network included three steps: (1)identification of the those up-or down-regulatedcircRNAs, (2)appliedmiRanda to predict miRNA target sites on differentially expressedcircRNAs, and (3)obtained the list ofdifferentially expressedmRNAs and selected their miRNA targetsfrom mirTarBase database with experiment support.