Circular RNAs in peripheral blood mononuclear cells are more stable than linear RNAs upon sample processing delay

Abstract Circular RNAs (circRNAs) are a novel class of RNAs with closed loop structure. Blood circRNAs are widely acknowledged to be more stable than linear mRNAs, which show promising prospect to be liquid biopsy biomarkers for clinical applications. However, accumulating studies have demonstrated that sample processing delays have profound effects on blood transcriptome expression profiles, wherein knowledge remains elusive about the impacts of prolonged sample processing on blood expression profiles of circRNAs. We collected whole blood samples from three donors and isolated peripheral blood mononuclear cells (PBMCs) at six different incubation time points. We measured total RNA expression profiles using RNA sequencing (RNA‐seq) and investigated the differentially expressed circRNAs, mRNAs and lncRNAs upon blood processing delay. Meanwhile, we explored the underlying inducement of aberrant expression of circRNAs against their corresponding mRNA transcripts. Finally, we utilized rMATS‐turbo and CIRI‐AS, respectively, to screen out differential alternative splicing (AS) events in linear mRNAs and circRNAs. Sample incubation at 4°C lasting to 48 hours (h) led to minimal effects to circRNAs' expression. However, it induced extensive alterations for mRNAs and lncRNAs when the incubation time was beyond 12 h. Additionally, only 2 h processing delays may result in profound impacts on AS events of linear mRNAs, while less impact on the equivalence of circRNAs. Our results suggested that PBMC circRNAs are stable upon sample processing delay, which are more suitable to be liquid biopsy biomarkers.

portable devices and automated workstations. [13][14][15] The standardization of these steps has greatly improved the whole procedure of blood sample processing, which may reduce technical noises and increase data reproducibility. 12 However, the incubation time from blood draw to RNA extraction is hardly to be standardized due to the limitation of working time, location and other situations. 10 The prolonged storage of blood samples has been observed to make substantial changes on the measured blood transcriptomes. [8][9][10][11] For example, Dvinge et al.
performed bulk RNA-seq of blood samples from four healthy donors, and found rapid transcriptional and post-transcriptional changes upon different blood incubation times at room temperature or cryopreservation. 8 In addition, Massoni-Badosa et al. performed single-cell RNA-seq and single-cell ATAC-seq of human blood samples from two healthy donors and three leukaemia patients. They also concluded that ex vivo prolonged blood storage induced marked alterations of transcriptional profiles and chromatin accessibility at the single-cell level. 9 Similarly, Savage et al. performed multi-omic profiling of human peripheral blood samples at different handling time points and investigated the effect of delayed blood processing on the multi-omic datasets, including targeted bulk PBMCs transcriptomics, PBMC single-cell transcriptomics, cell numbers and plasma proteomes. 10 They found extensive changes of single-cell transcriptome and plasma proteome after 4 h incubation, while accumulating differences were observed for the targeted bulk transcriptomes and the number of immune cells during an overnight rest (18 h) after blood draw. 10 All these studies indicated that blood mRNA transcripts were sensitive to handing delays, which may confound the biological findings from blood mRNA expression experiments and the translation of blood mRNA signatures for disease management. 8,9 To overcome this limitation, many recent studies have proposed blood circRNAs as a new kind of blood RNA biomarkers for human diseases. 16 Unlike linear mRNAs, circRNAs have a unique circular structure that lacks free ends, poly(A) tails and 5′ caps. 17 This makes circRNAs resistant to de-adenylation, decapping and exonucleases. 17 Therefore, circRNAs are observed to be more stable than linear mRNA transcripts. 18,19 Specifically, the median half-life of cir-cRNAs was at least 2.5 times longer than that of their linear mRNA counterparts in mammary cells. 18 Additionally, the expression levels of circRNAs in serum exosomes were stable for serum samples that incubated at room temperature for up to 24 h. 19 These results hinted that blood circRNA expressions should be more robust than linear mRNAs under different blood incubation times. In contrast, Rochow et al. found gradually reduced expression values of circRNAs in kidney cancer tissue and cell lines with RIN (RNA integrity number) value reduction. 20 They suggested that circRNAs were subjected to degradative processes in clinical samples, which was similar to linear mRNAs. 20 The different conclusions made in previous studies make it difficult to predict how circRNAs will respond to blood incubation times in sample processing. Therefore, there is still an urgent need to investigate the effect of blood sample processing delays on circRNA expression profiles.
In this study, we measured the expression profiles of linear mRNAs, long non-coding RNAs (lncRNAs) and circRNAs, in human

| Expression quantification of PBMC RNA transcripts
For each RNA-seq dataset, we identified the expressed circRNA transcripts using CIRI-full 21 with GRCh38 reference genome, Ensembl 94 gene annotation and the default parameters. Next, we constructed a reference library of expressed PBMC circular transcripts by combining the de novo constructed circular transcripts in annotated human genes from CIRI-full output and the known blood circRNA transcripts from isoCirc catalogue. 22 Then, we quantified the expression values of both circular and linear RNA transcripts using AQUARIUM 23 with the compiled reference library of circular transcripts, Ensembl 94 gene annotation and the default parameters. We chose to use AQUARIUM for RNA expression quantification, since we have observed its superior performance in estimating the expression values of both circular and linear RNAs at the transcript level. 23 After calculating the transcripts per million (TPM) values for all linear and circular RNA transcripts in each RNA-seq dataset, we integrated all the expressed transcripts in 18 RNA-seq datasets. Finally, those lowly expressed transcripts (a transcript that has a TPM value smaller than one in more than four samples) were excluded for further analysis.
For circRNAs, the transcripts whose biotypes of parental genes were not protein coding or lncRNAs based on Ensembl 94 gene annotation were further filtered out.

| Differential expression analysis
To investigate the transcriptome changes between different blood incubation time points, we imported the transcript expression profiles from AQUARIUM output using tximport 24 and calculated the expression differences of both circular and linear RNA transcripts using DESeq2 with likelihood ratio test and apeglm shrinkage method. 25 We chose DESeq2 for differential expression analysis, since it has better performance for alignment-free isoform quantification tools. 26 Transcripts or genes with |log 2 (fold change)| > 0.5 and adjusted p-value < .05 were considered as significantly differential expression. Among them, some were further classified into newborn or degraded transcripts, which suggests the dynamic gain or loss of RNA transcripts upon incubation. We defined the transcripts with a TPM value larger than one in at least two replicated samples at a time

| Differential alternative splicing analysis
Other than expression abundance, alternative splicing (AS) event can produce various RNA transcripts from one gene. To identify the changes of AS events between PBMC samples at different incubation time points, we first detected the AS events of linear mRNAs and circRNAs using rMATS-turbo 28 and CIRI-AS, 29 respectively. Both methods used a ratio value (Ψ) to estimate the inclusion possibility of a targeted exon. Next, we calculated the difference of Ψ values (△Ψ) between replicates at different incubation time points for each exon.
Then, we used paired t-test to calculate the p-value of differential splicing events with △Ψ, of which a stringent threshold, p-value < .05 and |△Ψ| ≥ .05, was adopted to define significantly differential AS. Finally, the number of abnormal AS at each time course was normalized by dividing the total number of identified AS exons of CIRI-AS or rMATS-turbo.

| Functional enrichment analysis
To explore the biological functions of blood incubation-related transcripts, we performed the Gene Set Enrichment Analysis (GSEA) 30 using gseGO() function and visualized the enriched gene sets using enrichplot() function in clusterProfiler package. 31 Biological pathways (BP) that have p-value less than .05 were considered as significantly enriched.

| Expression landscape of circular and linear RNA transcripts in human PBMC samples
Whole blood samples were collected from three healthy donors in anticoagulant blood collection tubes and were immediately incubated at 4°C for six scheduled time intervals, including 0 h, 2 h, 6 h, 12 h, 24 h and 48 h. PBMCs were subsequently isolated from these 18 samples, and total RNAs were extracted for transcriptome profiling. Although the RIN value decreases with the incubation time interval (Figure 1A), the quality of extracted RNAs was consistently good in all samples ( Figure S1). The average RIN value of all these samples was at 9.3, and the sample with the lowest RNA quality had a RIN value at 7.8 ( Figure 1A). These RNA samples were used for RNA-seq library construction and transcriptome profiling.
For each RNA-seq data, circRNA transcripts were identified using a home-built computational pipeline (see Materials and Methods).
Then, the expression of circRNAs, linear mRNAs and lncRNAs were quantified. RNA transcripts with low expression abundance in these blood samples were further filtered. Additionally, 62 circular transcripts were excluded due to biotypes of corresponding genes not being protein coding or lncRNA, of which 58 were pseudogenes, 1 misc_RNA, 1 snoRNA and 2 TR genes (see Materials and Methods).
Finally, a repertoire of expressed RNA transcripts in all these PBMC samples was constructed. In this PBMC transcriptome repertoire, 41,936 RNA transcripts were expressed in total. Among them, 5007 (11.9%) were circular transcripts, and 2709 (6.5%) were linear lncRNAs. The remaining 34,222 (81.6%) transcripts are coding mRNA transcripts ( Figure 1B). As expected, the coding mRNA has the highest number of expressed transcripts ( Figure 1B) and expression values ( Figure 1D). Although lncRNAs have the smallest number of expressed transcripts, their expressed abundance (22.8%) is higher than that of circRNAs, which account only 1.8% of the total expressed RNA transcripts ( Figure 1D). This can be explained by the lowest mean expression value of circRNA transcripts ( Figure 1C), since most circRNAs were lowly expressed in mammalian samples. 32 For exonic circRNAs (ecircRNAs), most were composed of no more than five exons ( Figure 1E). Additionally, the length of most ecircR-NAs is less than 1000 base pairs ( Figure 1F). In general, this de novo constructed PBMC transcriptome has similar characteristics that were observed in previous studies. 33,34

| Dynamic expression changes of circular and linear RNA transcripts upon incubation
To quantify the effects of sample processing delays from blood collection to PBMC isolation on blood transcriptome, we first in-  Table S1). The number of dysregulated transcripts gradually increased with the time interval of sample processing (Figure 2A-C).
This is consistent to the findings observed in previous studies of mRNA transcripts. [8][9][10] Comparing different types of RNA transcripts, we found that circRNAs had the least number of transcripts that were dysregulated in prolonged handling procedures, while mRNAs had the largest number of induced changes (Figure 2A-C). This is still the case when the changes were normalized to the proportion of dysregulated transcripts by dividing the total number of transcripts in the class ( Figure 2D). Taken together, our results sug-  Figure 3D). Furthermore, we found that the overlap of differential transcripts at certain time point with those in the previous time point was gradually increased as well ( Figure S4). This indicates that the dynamic gain or loss of RNA transcripts is accumulated along the course of blood incubation. and immune response-related pathways ( Figure 4B). This suggests that circRNAs and mRNAs may perform different biological functions in response to external stimulus during sample handling delays. CircRNAs are likely to trigger signal cascade by acting as indirect regulators, while mRNAs tend to transport functional proteins to cell membrane and induce immune response by acting as direct executors. The co-operation of circRNAs and mRNAs may mediate the structure morphogenesis, and even the apoptotic or death of blood cells.

| Potential causes that induce circRNA dysregulation during incubation
CircRNA expression is the product of the transcribed expression level of circRNA host gene and its proportion of back-splicing. 17 Figure 5A,B, red and green dots). Notably, there were still some splice-derived dysregulated circRNAs whose parental genes showed no significant expression changes ( Figure 5A,B, purple and blue dots). The expression trends of representative circRNAs and their parental genes during incubation using TPM values, rather than fold change also confirmed the conclusions ( Figure S5). GSEA analysis showed that biological functions of the host genes of splice-derived dysregulated circRNAs ( Figure 5C) had clear differences with those of transcription-derived dysregulated circRNAs ( Figure 5D). Interestingly, the splice-derived dysregulated

| Alternative splicing of RNA transcripts during blood incubation
AS is one of the key processes of multi-exonic gene expression during pre-mRNA maturation, including skipped exon (ES), alternative 5′ splice site (A5SS), alternative 3′ splice site (A3SS) and retained intron (RI). Previous studies have found these events are common in cir-cRNA formation as well. 28,35 To gain insights of incubation-induced AS, we identified AS events of both circRNAs and mRNAs at each incubation time point and then performed differential alternative splicing analysis compared to the original time point (Figure 6). For linear mRNAs, we saw a gradual increase of all four AS events along the incubation course ( Figure 6A). Comparing to linear mRNAs, cir-cRNAs had a far smaller number of AS events upon blood incubation ( Figure 6A). When normalized by the number of identified AS exons of each RNA type, the ratio of AS events occurred in circRNAs was still far lower than linear mRNAs ( Figure 6B). This suggests that AS F I G U R E 5 Correlation of log2(fold_change) of circRNAs versus log2(fold_change) of corresponding linear RNA expression at 24 h (A) and 48 h (B). Red and green dots represent transcription-derived circRNAs that were upregulated or downregulated because of consistent upregulation or downregulation of their parental genes. Purple and blue dots represent splice-derived circRNAs that were upregulated or downregulated whose parental genes showed no significant expression changes. Grey dots represent circRNAs that had no differential expression as well as their parental genes. The top 20 enriched biological pathways of the host genes of splice-derived circRNAs (C) and transcription-derived circRNAs (D) by GSEA (p-value < .05) were also shown F I G U R E 6 Sample handling delay induced differential AS events: (A) The number of differential AS events of circRNAs and mRNAs at different incubation time points; (B) The percentage of differential alternative splicing of circRNAs and mRNAs among the total identified AS circRNA or mRNA exons by CIRI-AS or rMATS-turbo events of circRNA transcripts were more tolerated to handing delays than linear mRNAs. We further compared the distribution of parental genes that experienced differential AS events and observed significantly different size between circRNAs and mRNAs ( Figure S6A), which might be explained by the longer circRNAs than mRNAs ( Figure S6C). However, the expression abundance of circRNAs with differential AS events is statistically similar to that of mRNAs ( Figure S6B), although circRNAs with AS events are more likely to have higher expression abundances than mRNAs ( Figure S6D).
These results suggested that the differences in detecting AS events were less likely to be relevant to the size and abundance of linear RNAs and circRNAs.

RNA expression in PBMCs or whole blood is important indicators
of the host's immune status, and their aberrant expression is closely related to many disease conditions, creating favourable prospect for liquid biopsy. 36  Our conclusion is consistent to the results proposed in several previous studies that circRNAs are more stable 19,32 and have longer halflives 18 than linear mRNAs. In this study, we focused on the effect of sample processing delay on circRNAs expression. To minimize the sample-inherent biases, such as disease duration or severity, we used peripheral blood samples from healthy donors rather than patients. Although our study contains a relatively smaller number of 18 samples from 3 participants, our emphasis is on the variation between groups rather than inter-individual variation. Therefore, it is appropriate to dissect expression changes of transcripts over time using three biological replicates at each time point. 37 20 This difference may be explained by the reduced RNA quality of clinical samples in their study, which may not be the case for the blood samples upon incubation ( Figures 1A and S1). In another study, Savage et al. suggested that single cells were more active during sample incubation, and transcriptome alterations appeared earlier at the single-cell level. 10 Therefore, it is interesting to further analyze the dynamic changes of circRNA expression and evaluate its robustness at the single-cell level.
In addition, we found that circRNA dysregulation was mainly derived from the dysregulated expression of their parental genes ( Figure 5). However, GSEA analysis indicated that handling delay induced different changes between circRNA and mRNA transcripts ( Figure 4), underscoring that circRNAs were not simple by-products of their linear counterparts. Specifically, splice-derived dysregulated circRNAs tended to perform their biological functions by interacting with Pol II ( Figure 5C). Interestingly, some circRNAs have been experimentally validated to act as Pol II interactors. For example, circEIF3J and circBPTF can interact with U1 snRNP to form an RNAprotein complex, and then bind to Pol II at the promoter region to enhance the transcription of their parental genes. 38 Ci-ankrd52 can also associate with Pol II to regulate the expression of its parental gene by modulating the elongation of Pol II. 39 These suggest blood incubation can cause distinct circRNA changes, although these changes may be neglectable even for samples with 48 h incubation.
Finally, AS events are the post-transcriptional process to diversify transcriptome and proteome by adjusting incorporated exons or introns, which are ubiquitous in the formation of mRNAs and cir-cRNAs. 29,35 Particularly, dysregulation of AS has been highly associated with human diseases, and is potential diagnostic biomarkers or therapeutics targets. 40 Moreover, Dving et al. have proposed that sample incubation would induce isoform switch. 8 Herein, we systematically investigated the effects of sample delays on four types of AS events for both circular RNAs and linear mRNAs. We found blood incubation resulted in profound impacts on AS of mRNAs, but not circRNAs ( Figure 6). Therefore, it is imperative to take this technical bias into consideration when interpreting the results of differential AS analysis of blood mRNA transcripts. Meanwhile, we observed an enrichment of NMD-related genes in differentially expressed mRNAs upon blood incubation ( Figure 4B). Couple with the extensive mRNA AS events, NMD could be a post-transcriptional mechanism in regulating gene expression upon incubation. 35,41 Specifically, we proposed that dysregulated AS may create isoforms with premature termination codons and truncated proteins under environmental stress, and indirectly participate NMD-related pathways to accelerate cell death.
In summary, PBMC circRNAs have smaller transcriptome changes than mRNAs and lncRNAs upon sample processing delays, no matter the expression level or AS events. Therefore, circRNAs are superior to linear transcripts as the blood biomarker candidates, especially when the sample handling process of clinical blood samples cannot be normalized.

ACK N OWLED G EM ENTS
This work was funded by grants from National Natural Science Foundation of China to WG (61571109) and from Innovation Team and Talents Cultivation Program of National Administration of Traditional Chinese Medicine (ZYYCXTD-C-202208).

CO N FLI C T O F I NTE R E S T
The authors have no conflicts of interest to declare.

DATA AVA I L A B I L I T Y S TAT E M E N T
The raw data of this study are available from the corresponding author upon reasonable request.