Abstract Circular RNAs (circRNAs) participate in complex diseases and have emerged as a novel biomarker and therapeutic target. It is challenging to identify circRNAs from total RNA-seq data and to validate and quantify their levels by qRT-PCR. We sought to characterize the circRNA transcriptome from RNA-seq with total RNA from septic peripheral blood mononuclear cells (PBMCs), assess/compare commonly used aligners, and optimize validation methods, considering sepsis is the leading cause of death in intensive care units (ICU) and there is no effective treatment. PBMCs were isolated from sepsis patients before and after intensive care, and RNA sequencing was conducted with total RNA. RNA-seq data for circRNAs were analyzed using four pipelines. A circRNA transcriptome has been generated, and four alignment and annotation pipelines have been evaluated. TopHat was found to be the most sensitive aligner while MapSplice was found to be the most accurate. Overall, circRNA expression in septic PBMCs were found to be more abundant and less diverse at ICU admission (ICU-AD) compared with ICU discharge (ICU-DC). There were three differentially expressed circRNAs between ICU-AD and ICU-DC. The addition of reverse primers to a reverse transcription reaction improved reproducibility and accuracy of qRT-PCR. Among these altered circRNAs, two different isoforms of circular RNA from the ASPH gene (circASPH(2,4), and circASPH(2,3,4) with the same backsplicing junction (BSJ) were identified and specifically quantified by qRT-PCR. All the two isoforms were highly expressed at ICU-AD and circASPH(2,3,4) levels positively correlated with the length of stay in ICU. The circRNAs predictably bind to miR-670, miR-7975, miR-6818, and miR-384 and interact with proteins TLR2 and ZC3H12D. In conclusion, it is feasible to identify circRNA transcriptome from RNA-seq with total RNA. Moreover, RT-PCR followed by gel electrophoresis is important to identify/ distinguish a new isoform of a circRNA with the same BSJ. CircRNA transcriptome in septic PBMCs evolves as infectious severity changes. CircASPH(2,3,4) has a clinical potential as an indicator of disease severity. Supplementary Information The online version contains supplementary material available at 10.1038/s41598-025-05652-3. Keywords: Sepsis, PBMCs, Circular RNA, Circular RNA ASPH, miRNA, RBPs Subject terms: Sequencing, Biotechnology, Molecular medicine, Predictive markers, Infectious-disease diagnostics Introduction Circular RNAs (circRNAs) are single-stranded RNA transcripts that are covalently linked at a back-splicing junction (BSJ), a non-conventional splice junction in which the downstream splice donor is linked to the upstream splice acceptor^[44]1. CircRNAs were found to be biologically significant, regulating gene expression through microRNA (miRNA) sequestration, protein interaction, and even translation^[45]1–[46]5. There is increasing interest in the putative role of circRNAs as a biomarker or therapeutic target in diseases for which no specific treatment currently exists. To accurately identify or discover circRNA expression profiles is the first key step for circRNA research. Currently, RNA sequencing (RNA-seq) with RNase R digestion is the commonly used method to enrich and identify circRNA, but this occurs at the expense of losing linear RNA. Additionally, RNase R digestion reduces the accuracy of circRNA quantification because some circRNAs are not resistant to RNase R, while some linear RNAs are resistant to RNase R due to secondary structures^[47]1. Total RNA-seq has been emerging as a platform for circRNA investigation. However, it poses a challenge to identify circRNAs from total RNA-seq data for most researchers without strong bioinformatics or computing knowledge. For circRNA identification algorithms, many studies have characterized their variation in sensitivity and specificity^[48]6; an increasing number of algorithms act as annotation programs that accept alignment data from multiple aligners, and the contribution of aligners to variation in circRNA detection is unknown. More research is required to evaluate the methodology of extracting circRNA reads from total RNA-seq, such as the optimal combination of aligners, annotation programs, parameters, and filtering conditions. Furthermore, due to low expression levels (1/10 to 1/1000 of their home genes) and lack of a poly-A tail, reproducibility and accuracy for quantifying circRNA levels by qRT-PCR becomes very challenging, despite being the most commonly used method to validate circRNA expression. Sepsis is a dysfunctional inflammatory response to infection, leading to multiorgan dysfunction or long-term immune suppression^[49]7. Sepsis affects up to 1–2% of all hospitalized patients, and is a major cause of death in the intensive care unit (ICU)^[50]8. Current treatments for sepsis rely on the detection of inflammatory biomarkers and timely administration of antibiotics, as there are no FDA-approved pharmacological interventions specific for sepsis^[51]8. Recently, a RNA-seq study reported that monocytes from peripheral blood mononuclear cells (PBMCs) of sepsis patients due to community-acquired pneumonia had higher average normalized circRNA counts compared with healthy controls, and that only two circRNAs from the versican gene (VCAN) and chromodomain helicase DNA binding protein 2 gene (CHD2) were significantly changed between the two groups^[52]9. However, there are no reports regarding the dynamic alteration of circRNA transcriptome in PBMCs from the same septic patients between the time of admission into ICU (ICU-AD) and discharge from ICU (ICU-DC), which could provide more precise and accurate data of circRNAs involved in sepsis. Moreover, whether circRNAs could be a predictor of clinical outcome, serving as potential biomarkers for diagnosis or prognosis of sepsis is not clear. The objectives of our study were to establish a method to identify circRNAs from total RNA-seq by comparing the sensitivity and accuracy of different aligners and to introduce a strategy for obtaining consistent and accurate qRT-PCR results. This study also aimed at determining dynamic circRNA transcriptome in PBMCs of sepsis patients between ICU-AD and ICU-DC to identify a biomarker for evaluation of disease severity. Results Dynamic circular RNA transcriptome identified by four alignment pipelines in septic PBMCs at ICU admission and discharge Principal component analysis and alignment with GRCh38.p13 transcriptome were first conducted using the total RNA-seq data to check RNA-seq quality. The results show clustering of samples from each condition (shown in Fig. [53]S1), and our RNA-seq data have an over 90% mapping rate to the transcriptome (Fig. [54]S2). To assemble a complete and comprehensive profile of the circRNA landscape in septic PBMCs, our RNA-seq data were analyzed by the four annotation pipelines (TopHat, STAR, MapSplice, and BBduk) to identify circRNAs. The BBduk pipeline used reference libraries of BSJs from circRNAs documented in circBase or circAtlas. A flowchart of bioinformatic pipelines is shown in Fig. [55]1. Fig. 1. [56]Fig. 1 [57]Open in a new tab A flowchart of bioinformatic pipelines. Overall, the total number of circRNA species in PBMCs aligned and annotated by the TopHat, STAR, and MapSplice pipelines were similar, while the BBduk pipeline reported around twice the number of circRNA species (Table [58]1). Across/Among all aligners, the greatest number of unique circRNA species identified was from the ICU-DC samples (Table [59]1, Fig. [60]2A). Most circRNAs identified by the fusion read aligners were documented in the circBase and circAltas, though there is a certain proportion of de novo circRNAs in septic PBMCs (Fig. [61]2A). In terms of overall abundance of circRNA reads from the 6 PBMC samples, we normalized the total BSJ reads to total mapped reads. We found there were no significant differences among the three aligners TopHat, STAR, and MapSplice, while BBduk reported around twice the circRNA abundance as the fusion aligners with circBase data, and around four times with circAtlas data (Fig. [62]2B). Two-way ANOVA revealed a significant difference between circRNA abundance at ICU-AD compared with ICU-DC across the different pipelines (p = 0.0147, Fig. [63]2C). We also found that BBduk reported greater circRNA abundance than the fusion aligners at ICU-AD. Table 1. Summary of circRNA identification from septic PBMCs. Aligner ICU admission^a ICU discharge^a Total^a Coverage (%) BBduk-circBase 714 928 1499 32.9 BBduk-circAtlas 919 1253 1875 TopHat 391 461 840 55.7 STAR 408 534 891 55.1 MapSplice 399 486 809 60.8 [64]Open in a new tab ^aCircRNA species identified with a read count ≥ 4 in at least two samples within the condition. Fig. 2. [65]Fig. 2 [66]Open in a new tab CircRNA annotation results by four annotation pipelines. (A) CircRNA species identified in each condition by the three fusion aligners, grouped into “known” circRNAs that are documented in circBase and circAtlas, or “de novo” circRNAs. Identified circRNAs were defined with a read count ≥ 4 in at least two samples within the group. (B,C) An abundance of circRNA reads, as a proportion of total mapped reads, graphed as individual samples. (C) Data in (B) graphed as mean and SD in each group. *, p < 0.05; **, p < 0.01; ***, p < 0.005. (D) Venn diagram of circRNA species identified by the four annotation pipelines in septic PBMCs in all 6 samples and within each group (n = 3). Identified circRNAs were defined with a read count ≥ 4 in at least two samples within the group. (E) Heatmap visualizing the expression of circRNA species identified by BBduk and at least one of TopHat, STAR, and MapSplice across all 6 PBMC samples. Data normalized as Z-scores across rows, and black represents a read of 0. (F) Correlation between read counts by BBduk and read counts by TopHat, STAR, and MapSplice. Dashed line represents simple linear regression analysis, with the R-squared value labelled. Venn diagrams were created to show the overlapping of identified circRNAs among the four pipelines (Fig. [67]2D). Notably, only 462 circRNA species were detected by all four pipelines, which is a little more than half of the total circRNA species identified by TopHat, STAR, and Mapsplice, but close to 1/3 of the total circRNA species identified by BBduk (Table [68]1). Taking a further comparison among three pipelines TopHat, STAR and Mapsplice, 513 circRNA species were identified by all three pipelines with the coverage of 55.7%, 55.1%, and 60.8%, respectively. Individual Venn diagrams for each time point of samples show similar trends to the total, regardless of group. These results indicate that different aligners, using different alignment algorithms, can lead to considerable variation in identified circRNAs. We also compared the expression of circRNAs identified through both the BBduk pipeline which used known circRNA BSJ sequences retrieved from circBase or circAtlas, and at least one of the fusion-read aligners. Overall, circRNA expression patterns are very similar across the four aligners as shown in the heatmap (Fig. [69]2E). Linear regression analysis comparing the expression patterns annotated by the individual fusion aligners to BBduk showed that circRNA expression characterized by MapSplice is the most similar to BBduk results (Fig. [70]2F, R^2 = 0.916). Since BBduk is assumed to select only reads that contain known circRNA BSJ sequences with no more than one mismatch, and MapSplice had the highest coverage among the pipelines compared, expression data from MapSplice were regarded as the most accurate results and thus used for the downstream analyses. CircRNA characteristics To characterize the circRNA landscape in PBMCs from septic patients, RNA-seq data aligned by MapSplice and annotated by CircExplorer2 were analyzed in the R statistical computing environment. Overall, the majority of circRNA species identified originated from chromosomes 1, 2 and 3, whereas chromosomes 22 and Y were found to contain the lowest amount of circRNA species (Fig. [71]3A). The vast majority of circRNA species were found to be uniquely expressed without other isoforms, and more genes in ICU-DC samples produced at least 2 circRNA isoforms compared with ICU-AD (Fig. [72]3B). About one third of circRNA species identified were back-spliced from exon 2 of their home genes (Fig. [73]3C, left panel), and the same pattern is observed when considering their abundance (Fig. [74]3C, right panel). Among circRNAs back-spliced from exon 2, about one quarter of them consisted of only exons 2 and 3 from their home genes (Fig. [75]3D). Identified circRNA species range from 84 to 5262 nucleotides (nt) in length, with 50% of circRNA species within 288–633 nt long (Fig. [76]3E). Fig. 3. [77]Fig. 3 [78]Open in a new tab CircRNA characteristics in PBMCs, as identified by MapSplice. (A) Chromosomal distribution of circRNA. circRNA species identified in septic PBMCs in all groups (Left), and in individual groups (Right). (B) Number of circRNA isoforms spliced from the same gene. (C) First exon of circRNAs containing the 3’ splice acceptor, grouped by circRNA species (left), and by circRNA abundance (right). (D) Last exon of circRNAs containing the 5’ splice donor, within circRNAs spliced from exon 2. (E) Frequency distribution of circRNA species by length, grouped into bins of 20 in width. The last bin combines all circRNAs from 2100 to 5500 nt. (F) Back-spliced-to-linear ratios of circRNAs, visualized as the FBP count of circRNA BSJ reads plotted against the FBP counts of the corresponding left or right linear spliced junctions, whichever is greater. Additionally, CIRCexplorer3-CLEAR was used to calculate the ratio of linear: circular RNA. Compared to their linear counterparts, most of the detected circRNAs were much lower expressed; however, circRNA species identified during ICU-AD generally had a higher normalized read count in relation to linear counterparts, than ICU-DC (Fig. [79]3F). Circular RNAs transcriptomes in septic PBMCs are different between the time of ICU admission and discharge Differential expression analysis was performed using the DESeq2 package in the R statistical computing environment and visualized as volcano plots (Fig. [80]4A). Three circRNA species were found to be differentially expressed between ICU-AD and ICU-DC: one circRNA species from the ASPH gene, which has a BSJ from exon 2 and exon 4, was overexpressed, while circRNAs from the VAPB and TUT4 genes were under-expressed at ICU-AD (Table [81]2, Fig. [82]4A). CircASPH is the circRNA with the highest mean reads. Fig. 4. [83]Fig. 4 [84]Open in a new tab Differential expression analysis of circRNAs in septic PBMCs. (A) Volcano plot depicting the differentially expressed circRNA between ICU admission and discharge. (B) Fold change of differentially expressed circRNAs and their linear counterparts. (C) Gene Ontology enrichment. Table 2. Differential expression of circRNA between ICU admission and discharge. circRNA baseMean log[2]FC lfcSE stat P value P adj circASPH(2,4) 226.48 3.111 0.672 4.632 3.62E−06 2.93E−03 circTUT4(2,3) 16.871 − 7.571 1.722 − 4.397 1.10E−05 4.45E−03 circVAPB(4,5) 13.796 − 7.277 1.764 − 4.126 3.70E−05 9.98E−03 circSUCLG2(8,9) 7.181 − 6.336 2.105 − 3.010 0.003 0.334 circABCA13(40,41) 11.775 5.001 1.672 2.990 0.003 0.334 circMETRNL(2) 9.830 − 6.801 2.286 − 2.975 0.003 0.334 circARHGAP18(6,7,8) 6.897 − 6.279 2.135 − 2.941 0.003 0.334 circSSH2(3,4,5) 8.727 6.654 2.274 2.926 0.003 0.334 circTMEM87A(3,4,5,6) 8.891 6.672 2.299 2.902 0.004 0.334 circPRDM2(7,8) 8.641 − 6.603 2.383 − 2.771 0.006 0.453 circASPH(4,5,6,7,8,9,10,11,12,13,14) 6.833 6.302 2.430 2.593 0.010 0.600 circPHKB(3,4,5,6,7) 5.335 − 5.907 2.312 − 2.555 0.011 0.600 circDOP1B(20,21,22,23,24) 5.251 5.916 2.351 2.516 0.012 0.600 circZBTB11(2,3) 5.189 − 5.863 2.358 − 2.487 0.013 0.600 circDNM3(7,8,9,10,11) 5.070 − 5.835 2.349 − 2.485 0.013 0.600 circELP3(11,12,13,14) 5.315 − 5.894 2.393 − 2.463 0.014 0.600 circPLXNC1(2,3,4) 78.066 1.798 0.734 2.451 0.014 0.600 circKMT2C(4,5,6) 13.044 − 3.739 1.528 − 2.447 0.014 0.600 circCSF2RA(4,5,6,7) 6.216 6.157 2.526 2.437 0.015 0.600 circCSF2RA(4,5,6,7) 6.216 6.157 2.526 2.437 0.015 0.600 circFBXO7(2) 4.798 − 5.759 2.392 − 2.408 0.016 0.619 circSLC39A10(2,3) 5.829 − 6.045 2.538 − 2.382 0.017 0.634 circBICD1(4) 4.871 − 5.775 2.449 − 2.358 0.018 0.642 circMETTL8(4,5) 5.848 − 6.035 2.584 − 2.335 0.020 0.642 circFUT8(3) 5.648 − 5.988 2.586 − 2.316 0.021 0.642 circRBM39(8,9,10) 4.501 − 5.662 2.446 − 2.315 0.021 0.642 circTASOR(2,3,4,5,6,7,8,9,10,11) 5.672 − 5.991 2.608 − 2.297 0.022 0.648 circRHOBTB3(6,7) 5.469 − 5.941 2.646 − 2.245 0.025 0.716 circSKA3(4) 5.012 − 5.828 2.660 − 2.191 0.028 0.783 circMAN2A1(2,3) 4.185 − 5.552 2.542 − 2.184 0.029 0.783 circAC068896.1(9) 4.073 − 5.528 2.569 − 2.152 0.031 0.797 circZFPM2(2,3) 4.016 − 5.495 2.555 − 2.151 0.032 0.797 circCCDC91(2,3,4,5,6) 28.099 − 2.636 1.238 − 2.130 0.033 0.812 circRASA3(2,3,4) 3.915 − 5.464 2.579 − 2.119 0.034 0.812 [85]Open in a new tab CircASPH(2,4) contains circASPH (2,4) and circASPH(2,3,4) , which are found in our validation experiments. Furthermore, we examined the relationship between the expression of circRNAs and their linear counterparts. It was found that the expression profile of most circRNAs at ICU-AD were similar to that at ICU-DC, and only about 0.4% circRNAs were significantly differentially expressed between these timepoints. A sizable proportion of circRNAs were found to be downregulated at ICU-AD, while very few linear RNAs were downregulated at the same timepoint (Fig. [86]4B). Gene Ontology analysis shows there are 4 pathways with > 100-fold enrichment for the differentially expressed circRNAs between ICU-AD and ICU-DC (Fig. [87]4C). Increasing reproducibility and accuracy of qRT-PCR for circRNA by an addition of reverse primers Since circRNA is a single-stranded RNA without a poly-A tail, cDNA was commonly synthesized with random primers, not oligo dTs. CircRNA is generally expressed at very low levels compared to its linear counterpart. To increase the cDNA synthesis of circRNA of interest, we added extra reverse primers of the circRNA along with hexamer random primers to the reaction. The resultant cDNA was used as a template for qPCR. As shown in Fig. [88]S3, the Ct value of the circRNA for the PCR using cDNA synthesized with an addition of reverse primers was dramatically decreased compared with the Ct value using cDNA synthesized with random primers only. More importantly, the variations among technical repeats became much smaller (< 0.5 Ct) in the PCR with cDNA synthesized with addition of reverse primers than those using cDNA made only by random primers, indicating improved reproducibility and accuracy of qPCR. Two isoforms of circASPH were identified by qRT-PCR, and circASPH(2,3,4) level positively correlated with length of stay in ICU To validate the RNA-seq results that circRNA back-spliced from the ASPH gene (circASPH) with the BSJ sequence of exon 2 and exon 4 was highly expressed in PBMCs at ICU-AD, qRT-PCR was conducted with divergent primers (F1 and R1) located at exon 2 (Fig. [89]5A, Table [90]S1), followed by gel electrophoresis with PCR products. It was found that there was a broad peak in the dissociation curve (Fig. [91]S4). Agarose gel electrophoresis of the PCR products showed there were two clear DNA fragments with different length around 200 bp (209 bp, and 155 bp) (Fig. [92]5B). Fig. 5. [93]Fig. 5 [94]Open in a new tab Validation of circASPH expression by qRT-PCR. (A) Diagram of two isoforms of circASPH and primers for qRT-PCR. Left: circASPH(2,4); and Right circASPH(2,3,4). (B) Gel electrophoresis of qPCR products. qPCR products were run on a 2% agarose gel and visualized by SYBG gold. The image was cropped, and its original image was presented in Supplementary Fig. [95]S7. (C) DNA sequencing with qPCR products specific to circASPH(2,4) and circASPH(2,3,4), respectively. Left: DNA sequence of the PCR products of circASPH(2,4) with reverse primer R2; Right: circASPH(2,3,4) with reverse primer R3. The BSJ and conjunction sequences between exons were underlined and pointed with an arrow. (D) A representative gel electrophoresis image of qPCR products for RNA with and without RNase R digestion. RT-PCR was conducted with 3 RNase R digested samples as well as their undigested controls using primers F1 and R1. PCR products were run on a 2% agarose gel and visualized by SYBG gold. The image was cropped and its original image was presented in Supplementary Fig. [96]S7. (E) Relative expression of circRNA isoforms circASPH(2,4) and circASPH(2,3,4), and linear ASPH mRNA by qRT-PCR. (n = 7). Data presented as a relative level of circRNA and ASPH normalized to ICU-AD. *, p < 0.05 (F) Correlation of circASPH(2,3,4) levels at ICU-AD and length of stay at ICU. n = 6; R^2 = 0.81, p < 0.05. To understand why there were two bands amplified, we searched and blasted with the transcriptome database. We found there might be two isoforms of circRNA having the same BSJ from ASPH: circASPH(2,4) and circASPH(2,3,4), named according to the nomenclature guide suggested by Chen et al.^[97]10. To specifically amplify each isoform, we re-designed qPCR primers F2/R2 specific for circASPH(2,3) where F2 spans across the 3’ end of exon 2 and 5’ end of exon 4, primers F3/R3 for circASPH(2,3,4) where F3 spans across the 3’ end of exon 2 and 5’ end of exon 3, shown in Fig. [98]5A and Table [99]S1. qPCR was conducted using each isoform specific primers. Agarose gel electrophoresis and DNA sequencing were applied to the resultant PCR products using each isoform specific primers. Each pair of isoform-specific primers had a clear and unique PCR band with the expected length (Fig. [100]5B) and a nice melting curve (Fig. [101]S4). DNA sequencing results showed that the two isoforms have the same BSJ sequence as well as their specific exon junction sequences, confirming the existence of two isoforms and PCR specificity (Fig. [102]5C). RNase R digestion was also performed with total RNA to remove linear RNA. The resultant RNase R digested RNA was subjected to cDNA synthesis, qPCR with primers F1 and R1, and gel electrophoresis. The results show that RNase R digested RNA samples had the same band pattern as their undigested controls (Fig. [103]5D), indicating the amplified two bands by primers F1 and R1 were true PCR products from the circRNAs. Using specific primers, two isoforms of circASPH and linear ASPH were validated to be significantly upregulated in ICU-AD compared with ICU-DC (Fig. [104]5E) with a larger sample size. We also performed correlation analysis between circRNA levels and the time stayed in ICU and hospital. We found that the level of circASPH(2,3,4) at ICU-AD positively correlates to the length of stay in ICU (Fig. [105]5F, R^2 = 0.81, p < 0.05), whereas there was no correlation between the circRNA and total hospitalization time. We also conducted qRT-PCR to determine the expression of circASPH(4,5,6,7,8,9,10,11, 12,13,14) which was not changed by RNA-seq with circExplore 2. As shown in (Fig. [106]S5), there was no significant change in the expression of circASPH(4,5, 6,7,8,9,10,11,12,13,14) in samples between ICU-AD and ICU-DC by qRT-PCR, which is consistent with the RNA-seq data. Moreover, we conducted qRT-PCR to determine the expression of unaltered circRNAs in RNA-seq including circFAM13B(8,9,10), circPCMTD1(3), and circVPS13C(10,11,12,13,14,15), which are spliced from the FAM13B, PCMTD1, and VPS13C genes. qRT-PCR results showed that the expression of circFAM13B(8,9,10), circPCMTD1(3), and circVPS13C(10,11,12,13,14,15) was not significantly changed between samples at ICU-AD and ICU-DC (Fig. [107]S6), aligned with the RNA-seq data. CircASPH interaction with miRNA and RNA binding proteins (RBPs) CircRNA could function as a miRNA sponge or bind to RBPs. To understand whether circASPH(2,4) and circASPH(2,3,4) potentially serve as a miRNA sponge, we conducted bioinformatic analysis to search miRNA binding sites on the two circRNA isoforms using different programs (miRDB, mirBase, and CircInteractome) (Table [108]S2). Bioinformatic analysis found that hsa-miR-6818-3p targets the BSJ and that each of hsa-miR-7975, hsa-miR-384, and hsa-miR-670-3p has two binding sites on circASPH(2,4) and circASPH(2,3,4), as predicted by miRDB. To construct a network of circASPH-miRNA-mRNA targets, potential target genes of the above four miRNAs were analyzed using TargetScan, miRDB, microT-CDS, and miRTarbase. The targeted genes predicted by all four programs were selected for network construction (Fig. [109]6A). The tight junction, antigen processing and presentation, Wnt signaling, Hedgehog signaling, and tuberculosis pathways were predicted to be regulated by at least two of the miRNAs. Fig. 6. [110]Fig. 6 [111]Open in a new tab Interactions of selected circRNA to miRNA and RBPs. (A) CircRNA-miRNA-mRNA network for circASPH, visualized in Cytoscape. CircASPH stands for circASPH(2,4) and circASPH(2,3,4) here. (B) Putative interactions of selected circRNAs to differentially expressed RBPs, predicted by RPISeq. To find potential RBPs, we extracted differentially expressed RBPs from our RNA-seq data. The mRNA expression for three RBPs, Toll-like receptor 2 (TLR2), cytoplasmic polyadenylation element binding protein 4 (CPEB4), and Zinc Finger CCCH-Type Containing 12D (ZC3H12D), were significantly changed in PBMCs between ICU admission and discharge (Table [112]3). Using RPISeq, we evaluated the possibility of interaction between the above three RBPs and circASPH(2,4) and circASPH(2,3,4). As shown in Table [113]3 and Fig. [114]6B, all two isoforms have a high likelihood of interaction with the RBPs. Table 3. Differentially expressed mRNA of RBPs between ICU admission and discharge and Probabilities of interaction between RBPs and circRNA exonic regions. RBP species Log2Fold change p-adj value Probabilities of interaction between RBPs and circRNA* circASPH(2.4) circASPH(2,3,4) CPEB4 1.5 0.0009416 0.8 0.9 TLR2 2.0 0.0036346 0.815 0.855 ZC3H12D − 2.7 0.0003226 0.78 0.795 [115]Open in a new tab *Number listed represents averaged RF & SVM Classifier scores calculated by RPISeq, the higher the score above 0.5, the higher the likelihood for interaction. On the other hand, some RBPs bind to the intronic regions flanking a circRNA, facilitating the joining of the splicing donor site and acceptor site, leading to back-splicing. We used CircInteractome and RPISeq to assess the likelihood of RBPs CPEB4 and TLR2 binding to the intronic flanking regions of the ASPH mRNA from which circASPH(2,4) and circASPH(2,3,4) are back-spliced since the two RBPs were up-regulated in ICU-AD (Table [116]3). The results show that up-regulated RBPs CPEB4 and TLR2 have very high probabilities of interaction with the 300-1000 bp flank regions of circASPH(2,4) and circASPH(2,3,4) (Table [117]4), indicating that CPEB4 and TLR2 may be responsible for the change of circASPH(2,4) and circASPH(2,3,4) levels. Table 4. Probabilities of interaction between TLR2 and pre-mRNA or flanking regions. RNA flanking region sequence Protein Probability RF classifier* SVM classifier* Upstream 1000nt of exon2 CPEB4 0.8 1 Downstream 1000nt of exon 4 CPEB4 0.8 1 Upstream 500nt of exon2 CPEB4 0.7 1 Downstream 500nt of exon 4 CPEB4 0.85 1 Upstream 300nt of exon2 CPEB4 0.55 1 Downstream 300nt of exon 4 CPEB4 0.85 1 Upstream 1000nt of exon 4 TLR2 0.65 0.99 Downstream 1000ntof exon2 TLR2 0.7 1 Upstream 500nt of exon2 TLR2 0.7 0.98 Downstream 500nt of exon 4 TLR2 0.75 0.99 Upstream 300nt of exon2 TLR2 0.55 0.99 Downstream 300nt of exon 4 TLR2 0.65 0.99 [118]Open in a new tab *Number listed represents average RF & SVM classifier scores calculated by RPISeq, the higher the score above 0.5, the higher the likelihood for interaction. Discussion In this study, we identified a circRNA transcriptome from RNA-seq data with total RNA, created a method of evaluating aligners that uses RNA-seq data from total RNA, and developed a more consistent and accurate quantitative method for measuring circRNA levels, to enable both circRNA and mRNA expression quantification at the same time from the same sample. Moreover, we, for the first time, characterized the circRNA transcriptome of PBMCs from sepsis patients at ICU-AD and ICU-DC and reported circASPH(2,4) and circASPH(2,3,4) are the most highly expressed circRNA in septic PBMCs, and the level of circASPH(2,3,4) at the admission to ICU correlated with the duration of a septic patient in the ICU. We identified differentially expressed circRNAs, explored potential roles of circASPH(2,4) and circASPH(2,3,4) in regulating RNAs and proteins, while predicting that up-regulated RPBs CPEB4 and TLR2 in the AD group facilitates the circRNA generation. RNA-seq with total RNA allows circRNA identification, expression quantification, and direct comparison to their linear counterparts in addition to the expression profiles of circRNA and linear RNA. However, it is very challenging to accurately identify circRNA from high abundance of mRNA when RNA-seq was conducted with total RNA, since circRNAs are generally expressed at low expression levels and have the same RNA sequences as their linear RNA, except the BSJ. Here, we generated and compared alignment data using multiple circRNA identification pipelines from total RNA-seq data. The three fusion-read aligners we used are well-established in circRNA research for circRNA-enriched RNA-seq data. TopHat-Fusion is a program that uses a single-end read alignment mechanism that is fairly equivalent to paired-end alignments, with a false-positive rate of 1 in every 32 million reads^[119]11. STAR first selects a point on each read to begin mapping to the genome and finds the longest alignment possible from that location. The process is repeated for the unmapped portion of the read. STAR then stitches the mapped segments together and selects the alignment with the highest quality as the final read sequence. STAR is much faster than the other two aligners, at the cost of requiring about 27 GB of free RAM for alignment to the human genome. STAR is optimized for rapid analysis of RNA sequencing results, and can produce alignments within hours^[120]6, but was found to have high false discovery rates and low F1 accuracy in abundance quantification^[121]12. MapSplice first divides all reads into segments of 20–25 bases in length and attempts to map them directly to the genome. Any unmapped segments are localized according to other segments of the read that are mapped, and splice junctions are identified and evaluated for confidence. Previous literature indicates that MapSplice is robust to background noise and does not experience reductions in sensitivity rate when noise is present^[122]13. MapSplice can detect large proportions of circRNAs in non-depleted RNA samples, making its use ideal when RNase R treatment is not applied to samples^[123]13. This algorithm boasts high sensitivity and robust identification of novel splice junctions, but it also requires more computational resources and time.^[124]13. Alignments from TopHat, STAR, and MapSplice were annotated by CircExplorer2, which extracts back-spliced reads and computes read counts. CircExplorer is extensively used in circRNA research, and its performance compared to other circRNA quantification software has been extensively investigated^[125]6,[126]13,[127]14. However, to our knowledge, we are the first to compare the performance of CircExplorer when different algorithms are used for its alignment step. The BBduk pipeline was used as an attempt to extract the “true positives” from the RNA-seq reads. Because RNA-seq reads can only be identified as circRNAs if they contain the BSJ, BBduk was used to search reads for BSJ sequences of known circRNAs, which was assembled using 15 bases from the 3’ and 5’ ends of mature circRNA sequences registered in circBase, or circAltas. Reads containing 25 consecutive bases from the 30-bp BSJ sequence with one mismatch or less were considered circRNA reads in this study. Overall, the BBduk pipeline identified twice the number of circRNA species and twice the circRNA abundance compared with the fusion-read aligners. Because there was a strong correlation between the BBduk and fusion aligner results for circRNA species jointly identified by BBduk and one of the fusion aligners, it is likely that some of the highly expressed “circRNAs” identified by BBduk are outliers and actually linear RNAs that contain a very similar sequence to a circRNA BSJ. The BBduk pipeline can be further improved to reduce false positives, for example, using 20 bases from the 3’ and 5’ ends of a mature circRNA sequence when generating a known circRNA library, but we may risk the missing of circRNAs expressed at low levels. The BBduk pipeline remains a fast and less resource-intensive method of screening known circRNA reads from RNA-seq data, although it cannot identify previously undocumented circRNAs. In our study, MapSplice identified the least number of unique circRNA species. CircRNA expression characterized by MapSplice is the most similar to BBduk. As BBduk is expected to select only reads that are no more than one mismatch from known circRNA junction sequences, MapSplice was regarded as the most accurate fusion-read aligner in the study. However, this does not mean that MapSplice is the gold standard of circRNA identification, as only about 60% of circRNA species identified by MapSplice were found by the other two fusion aligners and BBduk. Further investigation is required to reach a consensus on how conventional linear RNA-seq analysis should be modified for accurate circRNA identification and quantification. More recently, CirComPara was updated to integrate results from 9 existing aligners to maximize sensitivity and accuracy, though with a higher requirement for computational resources, time, and expertise^[128]15. Nonetheless, our results support recommendations that multiple circRNA aligners be used in conjunction to increase the sensitivity and precision of circRNA detection overall^[129]14,[130]16, and underscore the need for continued refinement of circRNA identification pipelines. Circexplore2 is a useful tool developed early and designed for short-read- based RNA-seq data. Circexplore2 has limitations in exon calling and is not able to identify some circRNA isoforms with the same BSJ and low read counts. In contrast, long-read-based sequencing or CIRI-full-like analysis showed significant differences between mRNA and circRNA sequences^[131]17,[132]18. The read length is the key determinant of circRNA isoform detection. Long read sequencing increases the ability to identify circRNAs, but it costs more. CIRI-full-like analysis can reconstruct full-length circRNA with paired-end data^[133]17, which helps identify circRNA isoforms and increases accuracy of exon calling. For single-end RNA-seq data or short-read-based RNA-seq data, CIRI-full-like analysis is not an option or does not improve circRNA calling. 75nt RNA-seq platform was conducted as single-end reads in our study, which is not ideal to CIRI-full-like analysis. Long-read-based sequencing and CIRI-full-like analysis will help improve accuracy and will be conducted in the future study. In our study, most identified circRNA species were found to originate from chromosome 1. This may be due to the sheer size of chromosome 1. We found that most of the circRNAs were back-spliced from exon 2 and that most genes expressed very few circRNA isoforms, while a few genes expressed several isoforms. Overall, our use of three fusion aligners, as well as documented BSJ sequences from circBase or circAtlas, highlights the limited overlap of circRNAs identified through different pipelines, and suggests that previous circRNA identification studies that used only one or two aligners likely yielded results with limited sensitivity and accuracy. qRT-PCR is the most commonly used quantitative method to measure circRNA and to validate RNA-seq data. A reverse transcription step is included to synthesize cDNA with hexamer random primers. The efficiency of the RNA to cDNA conversion varies across different transcripts^[134]19,[135]20. Extremely low levels of circRNA and short sequences relative to linear RNA (mRNA) in a cell make circRNAs less effectively transcribed than those abundant and longer sequences of mRNAs. As a result, much fewer cDNA copies of circRNA are used as a DNA template during PCR, which ends in a higher Ct value and big variations among technical repeats. Although circRNAs from the ASPH were the most abundant circRNA in septic PBMCs, it still needed 10 more cycles to get the same fluorescence value as linear ASPH. Moreover, it has a larger variation in Ct value among the technical replicates, which resulted in inconsistent data. RNase R digestion can help enrich circRNAs, but it loses linear RNAs and/or may lose those circRNAs sensitive to RNase R. In this study, we added a reverse primer specifically to a circRNA of interest to the reverse transcription reaction, to effectively and specifically reversely transcribe the circRNA. A relatively higher amount of circRNA cDNA was added to a PCR, which reduced Ct values by 8–10 cycles with < 0.5 Ct variances among technical replicates of PCR, increasing reproducibility and accuracy of qPCR. Of note, an addition of a reverse primer of circASPH only enhances reverse transcription of circASPH(2,4) and circASPH(2,3,4) and linear ASPH, without affecting other circRNAs. Theoretically, reliable and accurate qPCR results of each circRNA and its linear counterpart can be achieved with cDNA synthesized using a library of reverse primers for thousands of known different circRNAs and oligo dTs. RNA-seq is a high-throughput method, and it identifies a circRNA based on BSJ. Most RNA-seq needs to break down RNA into small fragments with size < 300nt. 50, 75, 100, and 150 bp platforms are commonly used, making it unable to discriminate different isoforms of circRNA with the same BSJ but with different components. In this study, we conducted RNA-seq with 75 bp platform as single-end read. Analyzed with CircExplorer2 and DEseq2 between the results of samples from ICU-AD and ICU-DC, circASPH(2,4) was identified as significantly altered. However, we observed the existence of two isoforms with the same BSJ sequence: circASPH(2,4) with 2 exons at length of 219 nt and circASPH(2,3,4) with 3 exons at the length of 264nt, identified by qRT-PCR according to dissociation curve, gel electrophoresis, and DNA sequencing results of PCR products. As such, gel electrophoresis following qRT-PCR helps to identify new isoforms of a circRNA and is suggested as a standard protocol for studying a circRNA especially for those RNA-seq data with short-length read or single end read. To design a couple of pairs of primers in different locations is helpful and important to detect any possible isoforms. As discussed earlier, long-read-based sequencing or CIRI-full-like analysis can reconstruct full-length circRNA and increase determination of circRNA isoforms. In fact, 150 nt platform has been used in high-throughput platforms, such as Hiseq, Novaseq, or MiSeq allows 300 nt paired-ends with limited depth. Furthermore, Nanopore/PacBIo has already allowed full-length circRNA sequencing of much longer variants without rRNA/linear RNA depletion^[136]18,[137]21–[138]24. These new platforms and analysis tool kits improve the identification capacity and accuracy. There were no significant changes in expression of circASPH(4,5,6,7,8,9,10,11,12,13,14), circFAM13B(8,9,10), circPCMTD1(3), and circVPS13C(10,11,12,13,14,15) between ICU-AD and ICU-DC measured by qRT-PCR, which is consistent with the results of RNA-seq. In this study, we used PBMCs isolated from sepsis patients at both ICU-AD and ICU-DC points. We first reported the circRNA transcriptome of septic PBMCs before and after intensive care and demonstrated that septic PBMCs at the time of ICU-AD express less diverse circRNA species, but at a higher total abundance, compared with those at ICU-DC. More importantly, our study is the first to show that the levels of circASPH(2,3,4) at ICU admission positively correlates with the length of stay at ICU, implying that circASPH(2,3,4) is a potential indicator of illness severity and that circASPH(2,3,4) level could be used as a stratification biomarker and potentially help physicians provide personalized care to septic patients. We also found that circRNA changes more dramatically than linear mRNAs, indicating that circRNA is more sensitive in response to infection than mRNA and can be a novel biomarker. RBPs regulate circRNA biogenesis and functions. In this study, we found three differentially expressed RBP TLR2, ZC3H12D, and CPEB4 mRNA with more than 1.5 fold-change between ICU-AD and ICU-DC from our RNA-seq data. TLR2 is a member of the TLR family and recognizes various foreign pathogens, including bacteria, fungi, virus, and endogenous materials^[139]25. ZC3H12D attenuates inflammation by reducing proinflammatory cytokine mRNA stability, playing an important role in immunity and inflammatory diseases^[140]26. The literature shows that there is a feedback loop between ZC3H12D and TLR signaling^[141]27. CPEB4 is associated with inflammation^[142]28. RPISeq analysis, which calculates the interaction possibility between a protein and an RNA sequence, found that both TLR2 and ZC3H12D had a high probability of interaction with circASPH(2,4), or circASPH(2,3,4), suggesting the circRNAs could function through interacting with TLR2 and ZC3H12D. Moreover, RPISeq analysis of TLR2 and CPEB4 and the flanking intronic regions of circASPH(2,4) and circASPH(2,3,4) shows that there is a high likelihood of both CPEB4 and TLR2 binding to the circRNAs’ flank intron sequences, to facilitate back-splicing events. This suggests that TLR2 and CPEB4 could promote the biogenesis of circASPH(2,4) and circASPH(2,3,4). Further studies should be performed to determine the exact role of CPEB4 and TLR2 in the biogenesis of septic circRNAs. Apart from RBPs, Our in silico analyses show circASPH(2,4) and circASPH(2,3,4) could bind to hsa-miR-6818-3p miR-7975, miR-384, and miR-670-3p, which are involved in multiple signaling pathways associated with the immune response, metabolism, cell cycle, vasopressin-regulated water reabsorption, and signal transduction associated with the immune response, metabolism, cell cycle, vasopressin-regulated water reabsorption, and signal transduction, autophagy and mitophagy, and tight junctions^[143]29–[144]33. In conclusion, we demonstrated that it is feasible to identify and acquire a circRNA transcriptome from RNA-seq with total RNA, and MapSplice is a better aligner with fewer false positive annotations. Adding reverse primers to cDNA synthesis could increase the accuracy and reproducibility of qRT-PCR for circRNA quantification. CircASPHs with the BSJ of exons 2 and 4 are the most upregulated circRNA, and circASPH(2,3,4) level is associated with the duration of ICU, having a clinic potential as an indicative factor for disease severity and patient stratification. The study offers insights into methodologies for the accurate and specific identification of a new circRNA. Materials and methods Sample collection Peripheral blood samples were collected from sepsis patients at ICU-AD and ICU-DC at London Health Sciences Centre. Written consent forms were collected from all participants. The study was approved by Lawson Health Research Institute and the Western University Research Ethics Board. All methods were performed in accordance with the relevant guidelines. Total RNA extraction Four mL of blood in EDTA tubes were diluted 1:1 with 10 × phosphate-buffered saline (Wisent Inc., Saint-Jean-Baptiste, QC) and overlaid on Ficoll-Paque PLUS (GE Healthcare, Mississauga, ON). After centrifugation at 630 g for 25 min and without breaking, PBMCs between the plasma and Ficoll layers were carefully collected and washed twice with phosphate-buffered saline. Total RNA from PBMC’s was isolated using the standard TRIzol method. Briefly, 1 mL TRIzol (Fisher Scientific, Toronto, ON) was added to isolated PBMCs and incubated at room temperature (RT) for 5 min. A 1:5 volume (200 µL) of chloroform (Fisher Scientific, Toronto, ON) was added, and the mixture was incubated again at RT for 5 min. After centrifugation for 15 min (14,000 g, 4 °C), the aqueous top layer was carefully extracted and mixed with an equal volume of isopropanol (Fisher Scientific, Toronto, ON). The mixture was incubated at RT for 5 min prior to centrifugation for 10 min (14,000 g, 4 °C). The pellet was collected, washed with 1 mL cold 75% ethanol, and centrifuged for 5 min (14,000 g, 4 °C). The pellet was dried for 10 min at RT prior to resuspension in 20 µL RNase Free H[2]O and heated for 15 min at 55 °C. RNA concentration was measured using a nanodrop (NanoDrop, Wilmington, DE, USA). Total RNA sequencing RNA-seq was conducted at the London Regional Genomics Centre (Robarts Research Institute, London, ON; [145]http://www.lrgc.ca). RNA quality was assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies Inc., Palo Alto, CA, USA) and the RNA 6000 Nano kit (Caliper Life Sciences, Mountain View, CA, USA). Total RNA samples were fragmented and rRNA depleted, and cDNA was synthesized, indexed, and cleaned-up using the Vazyme VAHTS Total RNA-seq (H/M/R) Library Prep Kit for Illumina (Vazyme, Nanjing, China) before PCR amplification. Libraries were then equimolar-pooled, and size distribution was assessed on a High Sensitivity DNA Bioanalyzer chip (Agilent, Santa Clara, CA, USA) and quantitated using the Qubit 2.0 Fluorimeter (Thermo Fisher Scientific, Waltham, MA, USA). PhiX spike-in controls were included. Sequencing was performed on an Illumina NextSeq 500 (Illumina Inc., San Diego, CA, USA) as single-end reads (1 × 76 bp) using a High Output v2 kit (75 cycles). Linear RNA expression was quantified through mapping to the Genome Reference Consortium Human Genome Build 38 patch release 13 (GRCh38.p13)^[146]34 using STAR v2.7.3a^[147]35 in Partek Flow (St. Louis, MO, USA). The sequencing data sets have been deposited to NCBI Sequence Read Archive (SRA) with the project ID PRJNA1191529 ([148]https://www.ncbi.nlm.nih.gov/bioproject/1191529). Circular RNA bioinformatic pipelines and alignment Circular RNAs were annotated using four pipelines: BBduk, TopHat-Fusion, STAR, and MapSplice2. A flowchart of bioinformatic pipelines is shown in Fig. [149]1. In both the BBduk and the TopHat pipelines, RNA-seq reads were first mapped to the hg38 reference genome (GRCh38.p13) in Partek Flow using STAR v2.7.8a. In the BBduk pipeline, RNA-seq reads were analyzed by the k-mer search tool BBduk^[150]36 (parameters: k = 25, mm = f, and hdist = 1) using a reference of known circRNA BSJs, built by retrieving the last 15 and first 15 bases of all circRNAs from circBase^[151]37 or circAtlas^[152]38. In the TopHat-Fusion pipeline, reads were first mapped to the hg38 reference genome in Partek Flow using STAR v2.7.8a, then unmapped reads were mapped to the reference genome (GRCh38.p13) and analyzed for fusion junctions using TopHat2 v2.1.1^[153]39 (parameters: –fusion-search –keep-fasta-order –no-coverage-search –fusion-min-dist 200). In the STAR and MapSplice2 pipelines, sequence reads were mapped to GRCh38.p13 and analyzed for fusion junctions using STAR (parameters: –chimSegmentMin 10) and MapSplice v2.2.1^[154]40 (parameters: –non-canonical-single-anchor –non-canonical-double-anchor –fusion-non-canonical, –min-fusion-distance 200, –min-map-len 25), respectively. Fusion junctions identified by TopHat-Fusion, STAR, and MapSplice2 were annotated using CircExplorer2^[155]41. Circular RNA differential expression analysis Differential expression analysis of circRNA between ICU-AD and ICU-DC was performed using the DESeq2 package^[156]42 in the R statistical computing environment (R Core Team 2020. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria). CircRNAs with a read count ≥ 4 in at least two samples, a fold-change ≥ 1.5 or ≤ − 1.5, and adjusted p < 0.05 were defined as differentially expressed, and their home genes were analyzed with Gene Ontology biological process terms. CircRNAs of interest were selected based on differential expression results for further validation. RNase R digestion 4 µg total RNA was digested with 20 U RNase R (ABM, Canada) in a 50 µl reaction for 2 h at 37 °C. 5 µl of RNase R digested RNA was used for cDNA synthesis. Ratios of linear transcript: circRNA RNA-seq results were analyzed using the circ-quant function from CIRCexplorer3-CLEAR to determine the ratios of linear transcript: circRNA from canonical splicing vs. back-spliced counts in fragments-per-billion. Quantitative real time PCR validation of circRNA expression Complementary DNA (cDNA) was synthesized from total RNA using OneScript® Plus cDNA Synthesis Kit (Applied Biological Materials Inc., Richmond, BC, Canada) according to the manufacturer’s protocol. Briefly, 2 µg total RNA was loaded into a reaction of 15 μL with 0.5 μM random primers, 0.5 mM dNTPs, in the presence or absence of 0.1 μM reverse primer R1 (Supplementary Table [157]S1), and incubated at 65 °C for 5 min before being cooled for 2 min on ice. Reverse transcription (RT) buffer and 1 µL Reverse transcriptase were then added up to 20 µL, and the mixture was incubated at 25 °C for 10 min, 55 °C for 15 min, then 85 °C for 5 min in a T100 thermal cycler (Bio-Rad Laboratories Ltd., Mississauga, Canada). Primers for aspartate β-hydroxylase (ASPH) mRNA, circASPH(2,4), circASPH(2,34), circASPH(4,5,6,7,8,9,10,11,12,13,14), circFAM13B(8,9,10), circPCMTD1(3), and circVPS13C (10,11,12,13,14,15) were designed and ordered from Invitrogen (Burlington, ON). 18S was used as a loading control. Primer information can be found in Table [158]S1. Quantitative real-time PCR (qRT-PCR) was performed using BlastTaq 2 × qPCR master mix (ABM Industries) with 200 nM primers and 2 µL of 1:10 diluted cDNA. Samples were heat-shocked at 95 °C for 3 min, followed by 40 cycles of 95 °C for 5 s, 60 °C for 10 s, and 72 °C for 10 s in a CFX Connect Real-Time PCR Detection System (Bio-Rad Laboratories Ltd.) or QuantStudio 3 (ThermoFisher Scientific, Mississauga, Canada). Relative circRNA expression was calculated using the 2^−∆∆Ct method. Gel electrophoresis PCR products were separated by 2% agarose gel electrophoresis, visualized with SYBR Gold Nucleic Acid Gel Stain (ThermoFisher Scientific), and imaged by a GelDoc Go Imaging System (BioRad). DNA sequencing PCR products were sequenced by the Genome Centre at Robarts Research Institute (Ontario, Canada) using PCR forward or reverse primers. CircRNA signal pathway analysis To predict potential miRNAs bound by circRNA, miRDB^[159]43 and mirBase^[160]44 were used, both of which allow “custom prediction” through inputting custom RNA sequences. CircRNA sequence with linear mRNA orientation, or circRNA sequence with 5’ end of the first exon joining 3’ end of the last exon was inputted to search potential miRNAs which were able to bind to circRNA. CircInteractome^[161]45 was also used to search potential miRNAs bound by circRNA. The following hierarchy was used to select miRNAs for further analysis and circRNA-miRNA-mRNA network construction: (1) miRNA binding site crossed the back-splice junction, (2) miRNA had more than one binding sites, and (3) miRNA was predicted by all three programs. TargetScan^[162]46, miRDB, microT-CDS from DIANA tools^[163]47, and miRTarbase^[164]48 were used to predict targets of miRNAs. Only genes predicted by all four programs were selected for pathway enrichment analysis and signal pathway network construction. Over-representation analysis using the WEB-based Gene Set Analysis Toolkit^[165]49 was conducted with KEGG pathways terms. The circRNA-miRNA-mRNA network was visualized using Cytoscape^[166]50. For circRNA-RNA binding protein (RBP) interactions, CircInteractome was first used to predict potential RBPs in circASPH. The list of human RBPs from the RNA-binding Proteins Database (RBPDB)^[167]51 was used to extract RBPs from our unpublished RNA-seq data. RBPs whose mRNA reads show a fold change > 1.5 or < − 1.5 and p -adj value < 0.05 between admission and discharge were selected. Likelihood of interaction between selected RBPs and circASPH was evaluated through RPISeq^[168]52, allowing for construction of the final circRNA-RBP interaction map. RPISeq was also used to determine the probability of CPEB4 and TLR2 interaction with the intronic regions flanking circASPH, to assess their potential contribution to circASPH biogenesis. Statistical analysis Volcano plots were produced in Microsoft Excel. Heatmaps were produced using the webtool Heatmapper^[169]53, and Venn diagrams were produced with the Venn Diagram Webtool by the VIB-UGent Center for Plant Systems Biology^[170]54. The student’s t-test, one- and two-way ANOVA followed by Tukey’s multiple comparisons test were conducted for comparison between two groups and three or more groups, respectively. Correlation between circASPH level and ICU stay time was conducted using GraphPad Prism 8. Statistical significance was defined as p < 0.05. Electronic supplementary material Below is the link to the electronic supplementary material. [171]Supplementary Material 1^ (300.4KB, pdf) Acknowledgements