Systematic identification and characterization of exon–intron circRNAs

The FEICP pipeline can be used to effectively identify EIciRNAs from HTS data

To identify EIciRNAs from HTS data, we developed finding EIciRNAs from paired-end RNA sequencing (FEICP) (Methods) (Fig. 1A). The first step of FEICP was to detect circRNAs from paired-end RNA sequencing (RNA-seq) data using the CIRI2 pipeline, a sensitive and accurate tool commonly used to identify circRNAs generated by backsplicing (Hansen et al. 2016; Zeng et al. 2017; Vromman et al. 2023). We used paired-end, but not single-end, RNA-seq reads from the FEICP pipeline, as paired-end RNA-seq covers longer regions and would be more reliable in detecting EIciRNAs. Next, based on the BSJs, the annotated host transcripts of the circRNAs were identified, and the HTS reads containing both BSJs and intronic sequences between circularized exons were used to identify candidate EIciRNAs. Selection of HTS reads containing both intronic sequences and BSJs ruled out the possibility of including intronic sequences from linear transcripts. In the bioinformatic identification of LIRs, introns are considered to be retained only when the following criteria are met in the HTS reads: the detection of exon–intron and exon–exon junctions and >90% read coverage of intronic sequences (Braunschweig et al. 2014; Middleton et al. 2017). The read coverage was calculated using all the sequencing reads, as the fragment length of BSJ reads from HTS is limited and generally cannot cover the whole intron. We then applied the same criteria for intron retention (IR) to detect EIciRNAs from the circRNAs (Fig. 1A).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 1.

FEICP pipeline and genome-wide analyses of EIciRNAs. (A) Workflow of the pipeline finding EIciRNAs from paired-end RNA sequencing (FEICP). Briefly, circRNAs were identified from paired-end high-throughput sequencing (HTS) data, followed by detecting exon–intron junctions or intron sequences from paired reads of backsplicing junctions (BSJs), and then retained introns and EIciRNAs were validated. (B) Correlation of Ct value examined by RT-qPCR and BSJ reads predicted by FEICP for 20 randomly selected EIciRNAs in HEK293 cells (left). Divergent primers to amplify EIciRNAs were indicated (right). (C) Pearson’s correlation of the BSJ counts of EIciRNAs commonly detected by HTS and Oxford Nanopore Technology (ONT) sequencing. (D) Correlation of EIciRNA BSJ reads (backspliced reads per million [BRPM]) from published HTS or ONT data sets of eight human tissues. (E) Distribution of full-length of EcircRNAs and EIciRNAs in HEK293 cells. The full-length of EcircRNAs and EIciRNAs was calculated by psirc and FEICP, respectively. (F) Boxplots showing the length of flanking introns of EcircRNAs and EIciRNAs (left) and Alu counts in their flanking introns (right). All introns annotated in the human genome and 2000 randomly selected introns were used as controls. (G) Proportion of EcircRNAs and EIciRNAs in cytosol and nucleus across seven human cell lines. For E and F, P-values were calculated using the Wilcoxon rank-sum test.

We assessed the EIciRNAs from sequencing reads with different fragment lengths and read lengths (Supplemental Fig. S1A,B). The number and length of detected EIciRNAs increased as the fragment size increased; however, the difference in the number of EIciRNAs among different fragment sizes was relatively small, and the range of EIciRNA lengths was similar across the various fragment sizes of HTS sequencing (Supplemental Fig. S1A). The number of EIciRNAs increased with increasing read length, whereas the increase was small when the read length reached 150 bp. The distribution of EIciRNA lengths remained similar across different read lengths (Supplemental Fig. S1B). We recommend constructing libraries with fragment sizes >250 bp and sequencing with read lengths no less than 150 bp to detect EIciRNAs. Using FEICP, we analyzed four biological replicates of paired-end HTS data from HEK293 cells (Supplemental Table S1). We randomly selected 20 FEICP-identified EIciRNAs and quantified them using RT-qPCR (Fig. 1B; Supplemental Table S2). The results showed that the Ct values and counts of BSJs of EIciRNAs were strongly correlated (Pearson’s correlation = −0.87) (Fig. 1B; Supplemental Table S2). In comparison, CIRI-full identified 287 EIciRNAs, ∼71% (203) of which overlapped with FEICP (Supplemental Fig. S1C). One of the key differences between CIRI-full and FEICP was that CIRI-full was based on BSJ and RO featured to reconstruct full-length circRNAs; thus, the overall length of the EIciRNAs captured by CIRI-full (median length 343 nt) was obviously shorter than that captured by FEICP (median length 821 nt) (Supplemental Fig. S1D).

We next examined the performance of the FEICP pipeline in distinguishing between IR events and cryptic exons. It was found that the proportion of EIciRNAs overlapping with circRNAs containing cryptic exons was similar between FEICP and previous ONT data also from HEK293 cells (isoCirc) (Supplemental Fig. S1E). Out of the total 1824 EIciRNAs, there were 1001 EIciRNAs with exons greater than two; although larger exons might be missed, even they might be included in the EIciRNAs. FEICP was unable to predict EIciRNAs with multiple introns, presumably owing to the relatively larger length of introns. FEICP was unable to detect exons or introns far away from the BSJ; even they might be included in the EIciRNAs. Another pipeline, CYCLeR, identified 2104 EIciRNAs in total, among which 130 and 64 overlapped with FEICP and isoCirc, respectively (Supplemental Fig. S1F). CYCLeR identified circRNA-specific features by comparing circRNA-enriched samples with control total RNA-seq samples, followed by inferring full-length circRNAs using a graph-based algorithm (Stefanov and Meyer 2023). This approach allowed CYCLeR to detect IR in circRNAs independently of BSJ reads, and identified EIciRNAs with a larger length spanning more exons and introns compared with the other three pipelines (Supplemental Fig. S1D,G). The distinct algorithm of CYCLeR could account for the relatively small overlap in EIciRNA identification between CYCLeR and the other software (Supplemental Fig. S1F). Comparable numbers of EIciRNAs, circRNAs, and proportions of BSJ reads were detected between our HTS data and isoCirc data, and 348 EIciRNAs were shared from FEICP identification and isoCirc (Supplemental Figs. S1F, S2A; Supplemental Table S1; Xin et al. 2021). Similar to the approach used by circtools to infer exon inclusion in circRNAs, we extracted the full-length of circRNAs detected by isoCirc, followed by quantifying their abundance in our RNase R+ HTS data and public RNase R– HTS data (Jakobi et al. 2019). It was found that ∼70% of all circRNAs and ∼75% of EIciRNAs were enriched significantly by RNase R treatment, and the proportion of enrichment increased as the BSJ count increased, indicating a good agreement of internal exon–intron structures of circRNAs between the two data sets (Supplemental Fig. S2B). Detection of circRNAs by isoCirc was dependent on rolling circle, which might be prone to miss relatively longer EIciRNAs, and it was noticed that isoCirc-detected EIciRNAs had a shorter length compared with those identified by FEICP (Supplemental Fig. S1D). Among 317 out of 343 EIciRNAs with larger lengths (1500–4000 nt) whose expression was detected by FEICP but not by isoCirc, we successfully validated 10 EIciRNAs with the highest BSJ reads via RT-PCR experiment (Supplemental Fig. S2C; Supplemental Table S2). We compared the correlation of the expression levels (the BSJ counts) of the 348 overlapping EIciRNAs and found that the correlation among four HTS replicates was 0.76–0.96, and the correlation between the HTS data and isoCirc was 0.31–0.60 (Fig. 1C), indicating the reliability of FEICP in quantifying EIciRNAs.

Next, we detected EIciRNAs in 17 human tissues using FEICP in a published study (Ji et al. 2019). EIciRNAs of eight out of these 17 tissues were also characterized in isoCirc data (Xin et al. 2021). We found that the correlation between the number of EIciRNAs in the eight tissues from the two data sets was 0.88 (Supplemental Fig. S2D), and the correlation between the normalized number of EIciRNA BSJs between these data sets was 0.73 (Fig. 1D). Taken together, these findings suggest that FEICP is reliable for identifying EIciRNAs from paired-end HTS data.

Features of EIciRNAs compared with those of EcircRNAs

A total of 1824 EIciRNAs and 14,445 EcircRNAs were identified in HEK293 cells (Supplemental Fig. S3A). Most of the EIciRNAs (1771, 97.09%) and EcircRNAs (14,165, 98.06%) were generated from protein-coding genes (PCGs) (Supplemental Fig. S3A). Among the 13,510 expressed PCGs in HEK293 cells, 9.81% (1326) could generate EIciRNAs, and 91.63% of EIciRNA-generating genes produced only one or two isoforms of EIciRNAs (Supplemental Fig. S3B). In comparison, 37.69% (5092) of the genes could generate EcircRNAs, and 63.75% of the EcircRNA-generating genes produced one or two isoforms of EcircRNAs (Supplemental Fig. S3C).

We next compared the sequence features of EIciRNAs and EcircRNAs. The length of the EIciRNAs (median length of 821 nt) was significantly longer than that of the EcircRNAs (median length of 512 nt) (Fig. 1E), which was reasonable considering that the EIciRNAs harbored retained introns. The flanking introns of both the EcircRNAs and EIciRNAs were significantly longer and contained more Alu elements than did the controls; these introns were all introns in the human genome or 2000 randomly selected introns (Fig. 1F). Previous studies have revealed that longer and more Alu sequences in flanking introns are features of backsplicing-generated circRNAs (Jeck et al. 2013; Westholm et al. 2014; Li et al. 2015). The flanking introns of EIciRNAs were significantly shorter and had fewer Alu elements than did those of EcircRNAs (Fig. 1F). EIciRNAs and EcircRNAs had stronger strength than the controls of random introns for both of the 3′ SSs of the upstream flanking introns and the 5′ SS of the downstream flanking introns (Supplemental Fig. S3D), whereas the strength of either the 5′ SS of the upstream flanking introns or the 3′ SS of the downstream flanking introns of the EIciRNAs and EcircRNAs was not significantly different (Supplemental Fig. S3D). Furthermore, the 3′ SS strength of the upstream flanking introns of EIciRNAs was significantly greater than that of EcircRNAs (Supplemental Fig. S3D).

To characterize the transcriptome-wide cellular localization of the EIciRNAs and EcircRNAs, we applied FEICP and CIRI2 to the nuclear and cytoplasmic poly(A)-minus RNA-seq data of seven cell lines from The ENCODE Project (The ENCODE Project Consortium 2012) to detect EIciRNAs and EcircRNAs, respectively (Fig. 1G). In these seven cell lines, essentially all the EIciRNAs were localized in the nucleus (Fig. 1G). In addition to showing the innate property of CIR, EIciRNAs displayed several features distinct from those of EcircRNAs.

The expression of EIciRNAs is highly tissue specific and finely regulated by RBPs

To systematically analyze EIciRNAs in humans, we collected 244 public paired-end RNA-seq data from 102 human cell lines and tissue types, and a total of 18,968 EIciRNAs from these data were identified with the FEICP pipeline, which overlapped with 2120 of the 8477 EIciRNAs identified by the FL-circAS database (Fig. 2A,B; Chiang et al. 2024). These samples represented ∼95% of the human PCGs, and the number of PCGs reached a saturation point (Supplemental Fig. S4A). By random sampling, we found that, overall, ∼16% of human introns contributed to the CIRs in EIciRNAs (Fig. 2A), and ∼37% of human PCGs could give rise to EIciRNAs (Fig. 2A). In addition, brain tissues and neural cells tended to have more EIciRNAs than did nonneural cells and tissues (Supplemental Fig. S4B). These data sets were generated by different laboratories and varied in many aspects, including sequencing library preparation, read length, sequencing depth, and data quality. Therefore, these data were not suitable for further in-depth analyses that require parallel comparisons. We then focused on the 38 ENCODE data sets, which were obtained from 19 human tissues with two biological replicates and had comparable and sufficient sequencing depths (373 million reads in depth on average). We found that both the number of EIciRNAs and the number of EIciRNA BSJs (normalized using backspliced reads per million [BRPM]) were generally greater in brain tissues than in nonbrain tissues (Fig. 2C; Supplemental Fig. S4C). Analysis of both the HTS and ONT data from eight human organs also showed that the brain and testis had more EIciRNAs than did the other tissues (Supplemental Fig. S4D). Seven brain tissues showed a distinct cluster of EIciRNAs, and the correlations of EIciRNA expression among brain tissues were greater than those among nonbrain tissues (Fig. 2D).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 2.

Identification and analyses of EIciRNAs in the human transcriptome. (A) The percentage of intron retained in circRNA (CIR) and genes generating EIciRNAs with the accumulating number of cell or tissue samples. The circles represent the mean values from 1000 iterations, and exponential curve fitting was applied. The corresponding equation and limit value are shown. (B) Overlap of EIciRNAs detected by FEICP and FL-circAS. (C) The number of EIciRNAs detected from ENCODE HTS data of 19 human fetal tissues. Twelve nonbrain tissues and seven brain tissues are labeled as gray and red, respectively. (D) Hierarchical clustering of EIciRNA BRPMs from 19 human fetal tissues. (E) Hierarchical clustering of Spearman's correlations of expression levels of RNA-binding proteins (RBPs) and EIciRNAs in 19 human fetal tissues. Cluster 1 and Cluster 2 are two row clusters for EIciRNAs, and Cluster A and Cluster B are two column clusters for RBPs. (F) Boxplots showing the tissue specificity index (tau) of PCGs, lncRNAs, circRNAs, and EIciRNAs, as well as parental genes of EIciRNAs and all circRNAs, from available HTS data sets (PRJCA000751) of 17 human tissues. P-values were calculated using the Wilcoxon rank-sum test. (PCGs) Protein-coding genes. (G) GO analysis of parental genes of EIciRNAs. The most significant five top GO terms are shown for six human tissues.

To analyze the potential roles of RBPs in the differential expression of EIciRNAs among tissues, we computed the Spearman's correlation coefficients of expression levels between all pairs of 1031 RBPs and 1628 EIciRNAs identified from the 19 tissues, followed by hierarchical clustering (Fig. 2E; Supplemental Methods). Two distinct RBP and EIciRNA clusters were revealed (Fig. 2E). A total of 1231 (75.6%) EIciRNAs formed cluster 1, which was negatively correlated with the expression of RBP cluster A (Fig. 2E). Cluster 2 had 397 (24.4%) EIciRNAs, and these genes were positively correlated with the expression of the RBPs in cluster A (Fig. 2E). Among the 1031 RBPs, binding profiles of 98 RBPs were characterized in the ENCODE eCLIP-seq data. Among them, 80 RBPs belong to cluster A and 18 belong to cluster B. We then calculated the binding density in both flanking introns, circularizing exons, and retained introns of EIciRNAs for each of the 98 RBPs (Supplemental Fig. S4E). We found that, compared with the RBPs in cluster B (positively regulated), those in cluster A (negatively regulated) had stronger binding density in both the circular exons and retained introns of these EIciRNAs, with no apparent difference in flanking introns (Supplemental Fig. S4E). These results suggested that RBPs complexly regulate the expression of EIciRNAs and that the majority of EIciRNAs are negatively correlated with RBPs. Protein–protein network analysis of the top 100 RBPs with the most negative correlations with EIciRNAs (Supplemental Methods) showed that they were involved in pathways such as RNA splicing, mRNA metabolic process, and post-transcriptional regulation of gene expression (Supplemental Fig. S4F). Several RBPs, including HNRNPK, SNRNP70, SRSF1, and SRSF2, were core nodes of this network (Supplemental Fig. S4F).

In one study, HTS data from 17 human tissues were obtained (Ji et al. 2019), and we analyzed the expression of PCGs, lncRNAs, and circRNAs (Fig. 2F). FEICP identified 4837 EIciRNAs from these tissues. Most of the EIciRNAs (3625, ∼75%) were expressed in only one tissue type, indicating that the EIciRNAs showed high tissue specificity (Supplemental Fig. S5A). We used the widely used tau method to quantify the tissue specificity of PCGs, lncRNAs, circRNAs, EIciRNAs, and the parental genes of EIciRNAs and circRNAs (Supplemental Methods) (Fig. 2F; Yanai et al. 2005). Compared with those of all PCGs, the other three types of RNAs showed greater tissue specificity, among which EIciRNAs showed the strongest tissue specificity (Fig. 2F). In contrast, we also found that the expression of parental genes of both EIciRNAs and circRNAs showed lower tissue specificity than did those of the all PCGs control (Fig. 2F). We then performed a gene enrichment analysis of those genes, whose mRNAs were expressed without strict tissue specificity but generated EIciRNAs in only one tissue type (Fig. 2G). We found that these parental genes were enriched in GO terms highly correlated with tissue identity (Fig. 2G). We further analyzed 70 ENCODE RNA-seq data sets from 18 tissues of four adults (Supplemental Fig. S5B). At least 50 EIciRNAs detected eight tissues from all individuals, after which hierarchical clustering of correlations between all pairs of these tissues was performed based on the expression of the EIciRNAs (Supplemental Fig. S5B). Biologically similar tissues were clustered together, which showed that EIciRNAs could be applied as markers of tissue identity (Supplemental Fig. S5B).

Genes generating EIciRNAs are more actively transcribed

Next, we examined the gene expression features of the EIciRNA-generating genes. We detected the CIRs and quantified gene expression from the 38 RNA-seq data sets used in Figure 2C. For comparison, we analyzed genes with LIRs using 71 poly(A)-plus RNA-seq data from 30 human tissues from ENCODE. All the expressed PCGs (TPM ≥ 1) were divided into 10 bins according to their TPM, and the fraction of LIRs or CIRs was calculated for each bin (Fig. 3A). The overall proportion of CIRs increased with increasing parental gene expression (Fig. 3A). However, the overall proportion of LIRs tended to negatively correlate with parental gene expression (Fig. 3A), and this feature of LIRs was reported in a previous study (Braunschweig et al. 2014). We compared the proportions of genes producing EIciRNAs or EcircRNAs in different expression-level bins and found that the proportions of genes generating either EIciRNAs or EcircRNAs generally increased with increasing gene expression levels (Fig. 3B). However, there was a notable decline for the proportion of EcircRNA-producing genes in the highest expression-level bin (Fig. 3B). To quantify the increasing rate for proportions of genes generating EcircRNAs or EIciRNAs, we calculated mean values of the proportions for each bin and then obtained their increment compared with the preceding bin (Fig. 3B). It was found that with the increase of gene expression levels, the accelerating trend increased gradually for genes generating EIciRNAs, whereas the trend decreased from the second bin on for genes generating EcircRNAs (Fig. 3B). To investigate the expression levels of genes that give rise to EIciRNAs versus EcircRNAs, we analyzed RNA-seq data from 12 cell lines and 17 tissues, whose sample preparation allowed analyses of EIciRNAs and EcircRNAs together (Supplemental Methods) (Fig. 3C; Supplemental Fig. S6A). In all these cell lines and tissues, the expression levels of EIciRNA-generating genes were significantly higher than those of genes generating only EcircRNAs but not EIciRNAs (Fig. 3C; Supplemental Fig. S6A). The expression levels of EIciRNA-generating genes were also greater than those of either all expressed PCGs (overall-control) or 2000 randomly selected PCGs (random-control) (Fig. 3C; Supplemental Fig. S6A).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 3.

Parental genes of EIciRNAs are actively transcribed. (A) Boxplots displaying the distributions of percentages of CIRs and introns retained in linear RNA (LIRs) in 10 groups (deciles) with ranked expression levels. (B) Boxplots displaying the distributions of percentage of genes generating EcircRNAs or EIciRNAs in 10 groups (deciles) with ranked expression levels. (Left) The mean proportion in each group is indicated as a point. (Right) The increment of mean proportion of genes generating EcircRNAs or EIciRNAs compared with the preceding bin is represented as a point for each bin. (C) Cumulative distribution of expression levels (transcripts per million [TPM]) of all genes, random genes, parental genes of EcircRNAs, and EIciRNAs in eight cell lines and eight tissues. (All genes) All expressed PCGs with TPM ≥ 1, (random genes) 2000 genes randomly selected from all genes, (parental genes of EIciRNAs) PCGs generating EIciRNAs, and (parental genes of EcircRNAs) PCGs generating EcircRNAs but no EIciRNAs. (D) Distributions of RNAPII and RNAPII Ser2P ChIP-seq signals around the transcription start site (TSS) regions of the indicated groups of genes in the K562, GM12878, and HeLa-S3 cell lines. (E) DNase-seq and ChIP-seq signals of H3K9ac, H3K27ac, H3K4me1, H3K4me2, H3K4me3, and H3K79me2 around the TSS regions of the indicated group of genes in the K562, GM12878, and HeLa-S3 cell lines. In C–E, P-values were calculated with the Wilcoxon rank-sum test to compare the parental genes of EIciRNAs and parental genes of EcircRNAs.

To investigate whether the relatively higher expression levels of EIciRNA-generating genes were associated with more active transcription, we assessed the RNAPII and RNAPII Ser2P signals on the corresponding genes in the K562, GM12878, and HeLa-S3 cell lines (Fig. 3D). These three cell lines were the only ones with chromatin-immunoprecipitation sequencing (ChIP-seq) data for both RNAPII and RNAPII Ser2P that were available in The ENCODE Project. We found that in these three cell lines, the ChIP-seq signals of both RNAPII and RNAPII Ser2P around the transcription start site (TSS) of EIciRNA-generating genes were significantly greater than those of genes generating only EcircRNAs but not EIciRNAs and were significantly greater than those of the overall-control and random-control genes (Fig. 3D). We also analyzed six cell lines for which only RNAPII ChIP-seq data were available (Supplemental Fig. S6B). We found that in these cell lines, the ChIP-seq signals of RNAPII at the TSSs of EIciRNA-generating genes were significantly greater than those of genes generating only EcircRNAs but no EIciRNAs and were greater than those of the overall-control and random-control groups (Supplemental Fig. S6B). We further assessed nascent RNAs in four cell lines using public transient transcriptome sequencing of nascent transcripts (TT-seq) data (Supplemental Fig. S6C). In these four cell lines, EIciRNA-generating genes were more actively transcribed than were the other genes (Supplemental Fig. S6C).

Six histone modifications (H3K9ac, H3K27ac, H3K4me1, H3K4me2, H3K4me3, and H3K79me2) and chromatin accessibility associated with active transcription were subsequently evaluated in K562, GM12878, and HeLa-S3 cells (Fig. 3E). We observed that, in all three cell lines, the DNase-seq and ChIP-seq signals of all six histone modifications were significantly greater around the TSSs of EIciRNA-generating genes than around those generating only EcircRNAs but not EIciRNAs and were also greater than those of the overall-control and random-control groups (Fig. 3E). However, the difference was weaker, or not significant, or not consistent among cell types for four histone modifications (H3K9me3, H3K27me3, H3K36me3, and H4K20me1) associated with repressed transcription (Supplemental Fig. S6D). To further investigate the correlation between the generation of EIciRNAs and the levels of these epigenetic marks, we compared the EIciRNA-generating genes with genes not producing EIciRNAs but having expression levels comparable to those of EIciRNA-generating genes (EIciRNA-match genes) (Supplemental Fig. S7A). We found that the signals around the TSSs of EIciRNA-generating genes were greater than those of EIciRNA-match genes for most activation marks (Supplemental Fig. S7B). When the top 10% of genes with high expression levels were examined as two groups with or without EIciRNA-generating, the signals around the TSSs of EIciRNA-generating genes were again generally higher than those of without EIciRNA-generating genes for most activation marks (Supplemental Fig. S7B). Overall, the genes generating EIciRNAs tended to be associated with active transcription and the corresponding epigenetic marks.

CIRs are distinct from LIRs and other introns

We wondered whether CIRs had some specific sequence features, and we analyzed three groups of introns, namely, the CIRs, LIRs, and non-circRNA-retained introns (NCIs), in HEK293 cells (Supplemental Methods) (Fig. 4A; Supplemental Tables S1, S3). NCIs refer to introns located between circularized exons in the genome but are not retained during backsplicing to form EcircRNAs. The NCIs and LIRs had little overlap, and a small portion of the CIRs overlapped with the LIRs (164 of 1333 CIRs, ∼12.30%) (Fig. 4B). Similar results were found by analyzing published RNA-seq data from four cell lines, namely, HEK293, HeLa, K562, and SH-SY5Y cells (Supplemental Fig. S8A). Therefore, CIRs were generally distinct from LIRs, and the sequence features of CIRs were subjected to further investigations.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 4.

Features of CIR compared with LIR and the other introns. (A) A schematic demonstration of the identification of LIR, NCI, and CIR. LIR was detected from poly(A)-plus RNA-seq data through IRFinder with the cutoff (IRratio ≥ 0.1). CIR was detected from RNase R–treated RNA-seq data using FEICP with the cutoff (percent intron retention [PIR] ≥ 0.1). NCI represents spliced introns of EcircRNAs with the cutoff (PIR ≤ 0.02). (B) Venn diagram showing the overlap of LIR, NCI, and CIR in HEK293 cells. (C) Illustrative diagram for training the neural network (NNetwork) to distinguish LIR, NCI, and CIR. (D) Performance of the NNetwork in distinguishing three groups of introns from each other. ROC curves and precision-recall curves were plotted for pairwise comparisons, with the corresponding area under the ROC curve (AUROC) and area under the precision-recall curve (AUPRC) values shown in parentheses. (E–G) Scatter plots showing the results of t-SNE analysis of the top 50 genomic features for distinguishing NCI and LIR (E), NCI and CIR (F), and LIR and CIR (G). The colors indicate the group of introns with which a feature was positively correlated in the pairwise comparison. The point size indicates the feature importance. (H) Density curves showing the genomic distribution of LIR, CIR, and overlapped IR across their host transcripts.

We analyzed the introns belonging to only one group (15,060 NCIs, 5729 LIRs, and 1169 CIRs) and extracted 1309 sequence features of these introns and the corresponding upstream exons (E1) and downstream exons (E2) (Fig. 4C; Supplemental Table S4). These features included five aspects, namely, transcript features, sequence conservation, RNA secondary structure, sequence motifs, and nucleosome positioning (Supplemental Methods) (Fig. 4C; Supplemental Table S4). All these features were subsequently used to train a three-layer neural network (NNetwork) (LeCun et al. 2015; Yeom et al. 2021), which aimed to predict whether an intron belonged to the NCI, CIR, or LIR (Fig. 4C; Supplemental Methods). We then evaluated the performance of this NNetwork model using receiver operating characteristic (ROC) curves and precision-recall curves (Fig. 4D). By computing the area under the ROC curve (AUROC) and the area under the precision-recall curve (AUPRC), we found that the NNetwork model showed the highest distinction between the NCIs and LIRs (AUROC = 0.93, AUPRC = 0.85), and the CIRs could also be distinguished from the NCIs and LIRs (AUROC = 0.77 and 0.79, AUPRC = 0.30 and 0.45, respectively) (Fig. 4D). These results showed that, for these sequence features, the NCIs and LIRs were the most different, whereas the CIRs and NCIs or the CIRs and LIRs were also significantly different, although to a lesser extent. We also extracted the 1309 genomic sequence features of introns for linear RNAs, which were used to predict whether an intron was retained by training a three-layer NNetwork (NNetwork-linear), which achieved an AUROC of 0.89 and an AUPRC of 0.37 (Supplemental Fig. S8B). When applying this two-classifier based on LIRs to distinguish CIRs from NCIs, we observed a noticeable decrease in both the AUROC (0.72) and AUPRC (0.27) compared with those of the three-class classifier (Supplemental Fig. S8C). We thus confirmed that CIRs were inherently different from LIRs.

We next assessed the importance of each of the 1309 features in distinguishing groups of introns from each other. The top 50 features distinguishing the two groups were extracted and further analyzed (Fig. 4E–G; Supplemental Table S4). Consistent with the findings of previous reports (Braunschweig et al. 2014; Yeom et al. 2021), a lower SS strength, higher intron GC content, and shorter introns were highly predictive of LIRs over NCIs and of CIRs over NCIs (Fig. 4E,F; Supplemental Fig. S9A). We compared the SS strength of the top 500 CIRs or LIRs with that of the corresponding bottom 500 IRs according to their percent intron retention (PIR) (Supplemental Fig. S9B). We found that, for LIRs, the top 500 had lower 5′ and 3′ SS strengths than did the bottom 500 (Supplemental Fig. S9B), whereas for CIRs, the top 500 had lower 5′ SS strength than did the bottom 500, and there was no significant difference in the 3′ SS strength (Supplemental Fig. S9B).

CIRs had higher sequence conservation scores at both the 5′- and 3′-ends of the introns compared with NCIs, which might imply the conserved regulatory functions of CIRs (Fig. 4F). The motif scores of several RBPs, including SRSF1, FXR2, and ELAVL1, had either positive or negative correlations with the CIRs. (Fig. 4F). We further integrated the eCLIP-seq data of 134 RBPs from ENCODE and compared their binding in CIRs over NCIs (Supplemental Methods) (Supplemental Fig. S9C; Van Nostrand et al. 2020). Eighty-nine of the 134 RBPs had significantly more binding sites in the CIR group than in the NCI group (Supplemental Fig. S9C). Additionally, among the 50 features predictive of CIRs over NCIs, two were for exons that flanked the retained introns in EIciRNAs; these exons were shorter in length and had lower translatability (Fig. 4F).

Several features, including a lower GC content, a longer length, a greater propensity to form a local secondary structure, a greater sequence conservation of the 5′-end of the intron, and a smaller position (toward the 5′-end of the gene) of the intron, were predictive of CIRs over LIRs (Fig. 4G; Supplemental Fig. S9A). By plotting the retained introns along their host transcripts, we observed an evident trend that CIRs peaked at the 5′ portion of the transcripts, whereas LIRs peaked at the 3′ portion of the transcripts in five data sets from four cell lines (Fig. 4H; Supplemental Fig. S9D). Consistently, the 164 introns shared by CIRs and LIRs in our HEK293 data set peaked in the middle of the transcripts (Fig. 4H). We also plotted the relative positions of the start exons and end exons of the EIciRNAs along the transcripts and found that the start exons of the EIciRNAs peaked at the 5′ region of the transcripts, whereas the end exons showed a relatively uniform distribution across the middle of the transcripts (Supplemental Fig. S9E). This finding was consistent with the reported localization pattern of start and end exons for circRNAs (Guo et al. 2014; Westholm et al. 2014; Rybak-Wolf et al. 2015). LIRs tend to locate to the 3′ region of genes, and a position closer to the 3′-end might be associated with the functional role of LIRs in regulating mRNA translation and decay (Braunschweig et al. 2014; Mauger et al. 2016; Schmitz et al. 2017; Wong and Schmitz 2022).

A genome-wide CRISPR-Cas9 screen for genes regulating CIRs identifies SRSF1

The introns retained in CIR were unique, and we subsequently sought to identify regulators of IR in EIciRNAs by performing a genome-wide CRISPR-Cas9 screening. To this end, we constructed an intron-split-GFP circRNA reporter construct (EIciGFP) in which the intron retained in EIciEIF3J or EIciPAIP2 was inserted between the split GFP sequence and the internal ribosome entry site (IRES) (Fig. 5A; Yang et al. 2017). EIciGFP could be produced when backsplicing occurred and the intron was retained, and it was found to localize in the nucleus (Fig. 5A; Supplemental Fig. S10A). When the intron was not retained, EcircGFP was produced; EcircGFP was localized mainly in the cytoplasm to serve as a translational template for the GFP (Fig. 5A; Supplemental Fig. S10A). As a control, we constructed an mCherry reporter, which was composed of the same promoter and used IRES as the EIciGFP reporter (Fig. 5A). The EIciGFP reporter with the intron retained in EIciEIF3J or EIciPAIP2 and the mCherry reporter were used to construct a stable HEK293 cell line, which was termed E-In cells or P-In cells (Fig. 5A; Supplemental Fig. S10A–E).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 5.

Genome-wide CRISPR screen identifies regulators of EIciRNA biogenesis. (A) Constructs for expressing GFP and mCherry protein reporters. The reverse complementary sequence (red arrows) in the flanking introns was used for the circularization of split GFP. The retained introns of two reported EIciRNAs (EIciPAIP2 and EIciEIF3J) were inserted between the split GFP (FP) and IRES sequences. Two GFP reporters (P-In for EIciPAIP2 and E-In for EIciEIF3J) were stably integrated into the genome of HEK293 cells along with the mCherry construct (Supplemental Methods). (B) Procedure of genome-wide CRISPR screen. Two reporter cells (P-In or E-In) were infected with the lentivirus sgRNA library and selected with 1 μg/mL puromycin for 7 d. The infected cells expressing mCherry proteins were gated with FACS, and 10% cell populations of high-GFP and low-GFP were used for genomic extraction, PCR amplification, and HTS. (C) Venn diagram showing the overlapped genes with significantly enriched sgRNAs in the screens of P-In or E-In reporters. (D) GO analysis of 80 genes overlapped in high-GFP of both P-In and E-In groups. (E) Ranks of eight RBPs enriched in the biological process of RNA splicing. Ranks were calculated using MAGeCK (Li et al. 2014). (F) Heatmap showing the relative fold change of seven EIciRNAs examined by RT-qPCR when eight RBPs were individually knocked down by shRNA. (G) Flow cytometric analysis showing the GFP expression in P-In (left) and E-In (right) cells upon SRSF1 knockdown. (shCtrl) shRNA with scrambled sequences, (shSRSF1-1, shSRSF1-2) two independent shRNAs against SRSF1. (H) Semiquantitative RT-PCR gels and quantification of the EIciGFP levels in P-In and E-In cells after SRSF1 knockdown. 18S rRNA was used as the loading control. Error bars represent standard deviation (SD) in triplicate experiments. P-values from two-tailed Student's t-test. (**) P < 0.01, (***) P < 0.001.

A whole-genome CRISPR-Cas9 knockout lentivirus library (Brunello), which was composed of 76,441 single-stranded guide RNAs (sgRNAs) targeting 19,114 genes in the human genome (Doench et al. 2016), was transfected into either E-In cells or P-In cells (Fig. 5B). Using fluorescence-activated cell sorting (FACS), we isolated cells that had the same mCherry fluorescence intensity as nonlentivirus-infected cells and collected the top 10% and the bottom 10% of the cells according to GFP expression (high-GFP or low-GFP group) (Fig. 5B). sgRNAs obtained from the input (bulk cells without sorting), high-GFP, and low-GFP groups were subjected to HTS (Fig. 5B; Supplemental Table S5).

Genes targeted by sgRNAs that were significantly enriched in high-GFP or low-GFP group compared with the input group were identified. From the high-GFP group, 854 and 709 candidate genes were identified from P-In and E-In cells, respectively (Fig. 5C), which had an overlap of 80 genes (Fig. 5C). For the low-GFP group, we obtained 304 and 36 candidate genes in P-In and E-In, respectively, which had an overlap of nine genes (Fig. 5C). In this study, we focused on 80 genes, whose expression was knocked out by CRISPR-Cas9, resulting in GFP expression. GO analysis of these 80 genes revealed enrichment in biological processes related to RNA metabolism and gene expression, such as RNA splicing and ribonucleoprotein complex biogenesis (Fig. 5D). Among these genes, 26 were known RBP genes (Gerstberger et al. 2014), and eight of them encoded factors associated with RNA splicing (Fig. 5E; Supplemental Fig. S10F). All these eight RBPs were also detected to be negatively correlated with EIciRNA expression (Fig. 2E; Supplemental Fig. S10G). We then examined these eight RBPs experimentally by knocking down with the corresponding shRNAs and investigated the expression levels of seven EIciRNAs, including EIciPAIP2 and EIciEIF3J (Fig. 5F). SRSF1 was the only gene whose knockdown led to significantly increased levels of all seven EIciRNAs (Fig. 5F). Knockdown of SRSF1 with two independent shRNAs in the reporter cell lines P-In and E-In resulted in increased EIciGFP and GFP levels (Fig. 5G,H; Supplemental Fig. S10H,I). Therefore, SRSF1 was identified as a candidate of negative regulators for CIRs through genome-wide CRISPR-Cas9 screening with two EIciRNA reporters.

SRSF1 regulates EIciRNA biogenesis

To further investigate the roles of SRSF1 in CIRs and EIciRNAs, we generated a stable SRSF1-knockdown (shSRSF1) HEK293 cell line (Supplemental Fig. S11A). Both poly(A)-plus and RNase R–treated RNA-seq were performed, after which the expression levels of mRNAs and EIciRNAs were quantified in the control and shSRSF1 cells (Fig. 6A,B; Supplemental Tables S6, S7).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 6.

SRSF1 binds to CIR to inhibit EIciRNA biogenesis. (A) Cumulative distribution of normalized EIciRNA expressions in HEK293 cells treated with shRNA targeting SRSF1. (shCtrl) shRNA with scrambled sequences, (shSRSF1) shRNA against SRSF1. (B) Changes of EIciRNA parental genes upon SRSF1 knockdown. (EIciRNA-up) Up-regulated EIciRNAs upon SRSF1 knockdown, (EIciRNA-down) down-regulated EIciRNAs upon SRSF1 knockdown, and (EIciRNA-stable) unaltered EIciRNAs upon SRSF1 knockdown. (C) Cumulative distribution of the IRratio of LIR in HEK293 cells treated with shRNA targeting scramble sequence or SRSF1. (D) UpSet plot showing the overlap of genes harboring LIR or CIR regulated by SRSF1. (E) Boxplots displaying the counts of SRSF1 binding sites in CIR and the length of CIR in the indicated groups. (F) Predicted binding motif from SRSF1 iCLIP-seq reads mapped to introns retained in up-regulated EIciRNAs upon SRSF1 knockdown. (G) Integrative Genomics Viewer (IGV) (Robinson et al. 2011) snapshot showing the SRSF1 iCLIP-seq signals in the EIciPAIP2 locus. Three SRSF1 binding GA-rich regions in the retained intron of EIciPAIP2 are framed as red dotted lines and labeled as #1, #2, and #3, respectively. (H) Semiquantitative RT-PCR gels and the quantification of EIciGFP levels in HEK293 cells after transfection of indicated mutation plasmids. The eukaryotic resistance gene (NeoR) of the plasmids was used as the loading control. In A and C, P-values were calculated with the Wilcoxon rank-sum test. In B, P-values calculated with the chi-squared test were shown for the comparison of proportion between the indicated group of genes against all regulated genes by SRSF1 knockdown. In E, P-values were calculated with the two-tailed Student's t-test. In H, error bars represent SD in triplicate experiments, and P-values were calculated with the two-tailed Student's t-test. (***) P < 0.001, (****) P < 0.0001.

The global expression levels of EIciRNAs significantly increased upon SRSF1 knockdown (Fig. 6A; Supplemental Table S6). In detail, 56% (658) of the EIciRNAs were up-regulated (EIciRNA-up), 33% (382) of the EIciRNAs were down-regulated (EIciRNA-down), and 11% (127) of the EIciRNAs were not significantly changed (EIciRNA-stable) (Supplemental Fig. S11B). When the expression of EIciRNA parental genes was examined, we found that, compared with those of EIciRNA-stable parental genes, significantly more EIciRNA-up parental genes had increased mRNA levels, whereas significantly more EIciRNA-down parental genes had decreased mRNA levels (Fig. 6B; Supplemental Table S7). By performing global run-on sequencing (GRO-seq) upon SRSF1 knockdown, we found that the global levels of nascent RNAs from all PCGs expressed in HEK293 cells were decreased (Supplemental Fig. S11C). Those EIciRNA-stable parental genes showed no significant difference to all PCGs in the nascent RNA levels (Supplemental Fig. S11C). EIciRNA-down parental genes showed significantly larger decreases, whereas EIciRNA-up parental genes had significantly fewer decreases in nascent RNA levels, compared with all PCGs (Supplemental Fig. S11C). It seemed that the overall effect of SRSF1 knockdown would lead to down-regulation in the transcription, whereas EIciRNA-up parental genes had a tendency to resist the down-regulation. Overexpression of FLAG-tagged SRSF1 resulted in a significant decrease in the nascent production of EIciRNAs and the corresponding mRNAs of their parental genes in the three EIciRNA/mRNA pairs examined (Supplemental Fig. S11D,E).

We then investigated the alterations in LIRs and CIRs upon the SRSF1 knockdown. We found that the IR ratio of LIRs significantly increased upon SRSF1 knockdown (Fig. 6C; Supplemental Table S8), consistent with the findings of previous reports in HepG2 cells (Middleton et al. 2017). The majority (1741/1879, ∼93%) of genes with significant changes in IR showed only regulated LIRs or regulated CIRs (Fig. 6D). Although most of the genes generally possessed multiple introns, the IR of most of the genes was unidirectionally regulated by SRSF1, with LIRs in 982/1041 (∼94%) genes being either up-regulated or down-regulated and CIRs in 642/700 (∼92%) genes being either up-regulated or down-regulated (Fig. 6D).

Reanalysis of published SRSF1 iCLIP-seq data from HKE293 cells revealed that there were more SRSF1 binding sites in the retained introns of EIciRNAs regulated than in those unregulated by SRSF1 (Fig. 6E; Howard et al. 2018). Moreover, there was no significant difference in SRSF1 binding to exons or flanking introns of EIciRNAs (Supplemental Fig. S11F). The length of introns was significantly longer in EIciRNAs regulated by SRSF1 than in those unregulated by SRSF1 (Fig. 6E). Compared with introns retained in up-regulated EIciRNAs, those in down-regulated EIciRNAs were longer and had more SRSF1 binding sites (Fig. 6E). When the number of SRSF1 binding sites was normalized to the intron length, there was no significant difference in the density of SRSF1 binding sites among the introns retained in the stable, up-regulated, or down-regulated EIciRNAs (Supplemental Fig. S11G). Consistent with these findings, the metaprofiles of SRSF1 binding to the three groups of EIciRNAs were not significantly different (Supplemental Fig. S11H). Therefore, we speculated that the complicated changes in CIRs resulting from SRSF1 knockdown were more related to the length and might also be regulated in combination with some other RBPs required further investigation.

Analysis of SRSF1 iCLIP-seq reads mapped to retained introns in the up-regulated EIciRNAs identified a purine-rich motif (Fig. 6F). We found via iCLIP-seq that introns of EIciPAIP2 and EIciEIF3J harbored several SRSF1 binding sites, and among them, three and two were found to be purine-rich in introns of EIciPAIP2 and EIciEIF3J, respectively (Fig. 6G; Supplemental Fig. S11I). We then experimentally examined the SRSF1 binding sites in introns of EIciEIF3J and EIciPAIP2 (Fig. 6G; Supplemental Fig. S11I). Random mutation of purines to pyrimidines in these sequences with EIciGFP reporters led to significantly increased levels of the corresponding EIciGFP (Fig. 6G,H; Supplemental Fig. S11I,J). Overall, we concluded that SRSF1 could suppress CIR and thus the biogenesis of a portion of EIciRNAs by binding to specific motifs in introns of parental genes, although the overall regulatory mechanisms of EIciRNA biogenesis are complicated.

Regulation of CIRs and EIciRNAs has functional roles in neuronal differentiation

To investigate the biological relevance of CIR and EIciRNAs, we applied an established neuronal differentiation model by inducing SH-SY5Y cells with retinoic acid (RA) (Supplemental Fig. S12A–C; Ross et al. 1983; Rybak-Wolf et al. 2015). HTS of RNAs from noninduced and RA-induced cells was performed to evaluate circRNAs and mRNAs, each with two biological replicates (Supplemental Table S9). We detected 384 high-confidence EIciRNAs in total (Supplemental Fig. S12D; Supplemental Table S9), among which 138, 164, and 299 EIciRNAs were identified on days 0 (noninduced), 3, and 6 of SH-SY5Y differentiation, respectively (Fig. 7A). Temporal clustering of these EIciRNAs revealed six clusters (Fig. 7B). To investigate the EIciRNAs that could promote SH-SY5Y differentiation, we focused on cluster 1, in which most EIciRNAs increased continuously during differentiation (Fig. 7B). Gene enrichment analysis revealed that the parental genes of these EIciRNAs were enriched in pathways associated with neurological activities and cell polarity (Fig. 7C). Moreover, the overall expression levels of the parental genes of these EIciRNAs increased significantly from day 0 to day 3, whereas the increase was not significant from day 3 to day 6 (Fig. 7D), implying that these EIciRNAs might play regulatory roles in the early stages of SH-SY5Y differentiation. When the expression of the 141 well-defined splicing factors was analyzed (Papasaikas et al. 2015), we found that their overall expression levels decreased during differentiation (Supplemental Fig. S12E,F).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 7.

Functional roles of EIciRNAs in neuronal differentiation. (A) Circos plots (Krzywinski et al. 2009) showing genomic distribution (hg38) and backspliced reads per million (BRPM) of EIciRNAs at day 0, day 3, or day 6 after retinoic acid (RA)–induced SH-SY5Y differentiation. Two biological replicates were used for analysis. (B) Fuzzy clustering showing the temporal expression patterns of EIciRNAs during SH-SY5Y cell differentiation. (C) Enriched GO terms of parental genes of EIciRNAs (cluster 1 in B). (D) Violin plots displaying the expression levels (TPM) of EIciRNA (cluster 1) parental genes during SH-SY5Y cell differentiation. (E) Heatmap showing the expression levels of 32 serine-rich (SR) and HNRNP proteins during SH-SY5Y cell differentiation (Papasaikas et al. 2015). (F) RT-qPCR analysis of EIciLIMK1 and LIMK1 mRNA expression in SH-SY5Y cells after SRSF1 knockdown. (shCtrl) shRNA with scrambled sequences, (shSRSF1-1 and shSRSF1-2) two independent shRNAs against SRSF1. (G) Nuclear run-on assay showing the nascent level of LIMK1 mRNA in SH-SY5Y cells after EIciLIMK1 knockdown with shRNA or ASO against the EIciLIMK1 BSJ. (Scramble ASO) ASO with scrambled sequences. (H) RT-qPCR analysis of EIciLIMK1 and LIMK1 mRNA expression in SH-SY5Y cells after overexpression of LIMK1 (left) and EIciLIMK1 (right). (I) Western blot showing the expression levels of LIMK1 and phosphorylated cofilin (p-cofilin) protein levels in SH-SY5Y cells in uninduced and RA-induced SH-SY5Y cells after EIciLIMK1 knockdown. GAPDH was used as the loading control. (J) Representative immunofluorescence (IF) images of F-actin in RA-induced SH-SY5Y cells after EIciLIMK1 knockdown (left). The F-actin fluorescence intensity was quantified with Fiji (right). N = 40. In D and J, P-values were calculated with the Wilcoxon rank-sum test. In F–H, error bars represent SD in triplicate experiments, and P-values were calculated with the two-tailed Student's t-test. (*) P < 0.05, (**) P < 0.01, (***) P < 0.001, (****) P < 0.0001.

Specifically, among the 32 serine-rich and HNRNP splicing factors that are known to regulate splicing by direct intronic binding, 27 showed a decreasing trend during differentiation (Fig. 7E). The opposite trend in expression between the levels of these splicing factors and EIciRNAs was consistent with the finding that the expression levels of the majority of RBPs were negatively correlated with those of EIciRNAs (Fig. 2E). We noticed that the serine/threonine kinase LIMK1, which was enriched in the GO terms “positive regulation of GTPase activity” and “positive regulation of axonogenesis,” could generate an EIciRNA (EIciLIMK1) in cluster 1 (Supplemental Fig. S13A). Considering that SRSF1 has previously been shown to play a role in regulating EIciRNA biogenesis (Figs. 2E, 4F, 5, 6; Supplemental Figs. S4F, S9C), we experimentally examined whether SRSF1 regulated EIciLIMK1 in SH-SY5Y cells (Fig. 7F). SRSF1 knockdown significantly increased in the expression of the EIciLIMK1 and LIMK1 mRNA (Fig. 7F; Supplemental Fig. S13B). RNA immunoprecipitation (RIP) experiments showed that SRSF1 bound to the flanking introns, the circular exons, and the retained intron of the EIciLIMK1 locus (Supplemental Fig. S13C,D).

To determine whether EIciLIMK1 plays a role in neuronal differentiation, we knocked down EIciLIMK1 with shRNAs or RNase H–based antisense oligonucleotides (ASOs) (Supplemental Fig. S13E) and found that the expression level of LIMK1 mRNA significantly decreased upon EIciLIMK1 knockdown (Supplemental Fig. S13E). Further, nuclear run-on assays showed that knockdown of EIciLIMK1 resulted in significantly decreased LIMK1 transcription levels (Fig. 7G). However, the expression level of EIciLIMK1 did not change upon LIMK1 mRNA knockdown or overexpression (Supplemental Fig. S13F; Fig. 7H). The mRNA levels of LIMK1 were unaffected by overexpression of EIciLIMK1 with plasmid (Fig. 7H). These results indicated that EIciLIMK1 could promote the transcription of the LIMK1 gene in cis.

Western blot analysis revealed that the levels of both LIMK1 protein and Ser-3-phosphorylation of cofilin (p-cofilin) significantly decreased after knockdown of either LIMK1 mRNA or EIciLIMK1 in SH-SY5Y cells (Supplemental Fig. S13G,H; Fig. 7I). Cofilin, a substrate of LIMK1, can be phosphorylated to p-cofilin, which is involved in the formation of F-actin during neuronal differentiation (Yang et al. 1998; Ben Zablah et al. 2021; Villalonga et al. 2023). Consistently, the mean F-actin intensity decreased significantly after knockdown of either the EIciLIMK1 or LIMK1 mRNA (Fig. 7J; Supplemental Fig. S13I,J). Taken together, these data suggested that SRSF1 could suppress the biogenesis of EIciLIMK1 and that the expression level of SRSF1 decreased to release this suppression during neuronal differentiation. EIciLIMK1 regulated the expression of its parental gene in cis, and up-regulated LIMK1 led to increased levels of p-cofilin, which facilitated the formation of F-actin during neuronal differentiation.