The circular transcriptome significantly differs between frail and robust individuals
CircExplorer2 and CIRI2 detected a total of 52,835 and 122,727 circRNAs, respectively in all the samples. A good detection overlap was observed between algorithms with 41,108 circRNAs in common, out of which 15,189 were detected with a minimum of 2 reads and, therefore, were considered bona fide circRNAs. Most of those circRNAs are poorly expressed, with 94.7% detected with less than 5 reads/sample (14,380 circRNAs). Of the highly expressed circRNAs, 766 (5.0%) are detected with 5-40 reads and only a few (43 circRNAs, 0.3%) are detected with more than 40 reads (Figure 1A).
When the origin of identified circRNAs was investigated, we found that 64.3% of the circRNAs (9,775 circRNAs) are located in a known gene locus, while 5,414 are intergenic. Remarkably, although we have detected 9,775 genic circRNAs, only 3,490 unique genes are identified, indicating that some of the circRNAs originated from the same genic region. Indeed, out of the 3,490 genic regions, 1,917 regions produce more than one circRNA, and 134 regions could be considered circRNA hotspots for producing more than 10 circRNAs (Figure 1B). Regarding the comparison between robust and frail individuals, out of the 15,189 circRNAs, 44.1% (6,692 circRNAs) are detected both in frail and robust participants while the rest are exclusive to one of the groups. There are 4,042 circRNAs detected in at least one of the pools from frail individuals and absent in all the robust ones. Similarly, we detected 4,455 circRNAs which are exclusively detected in some of the pools corresponding to robust individuals (Figure 1C).
For the circRNAs detected in both groups, we performed a differential expression analysis to report the expression differences that could be associated with frailty. Interestingly, we identified 89 significantly differentially expressed (DE) circRNAs, with a similar proportion of upregulated and downregulated circRNAs - 43 circRNAs were found to be upregulated in frail individuals whereas 46 were downregulated (Figure 2 and Supplementary Table 1).
In order to select circRNA candidates for their subsequent validation by qPCR, we included some additional criteria. Among the circRNAs found to be DE we filtered out those detected with a low number of reads (Base Mean < 5 reads/sample) and selected the top 5 differentially expressed circRNAs based on their Fold Change (FC > |1.8|). Additionally, among the group exclusive circRNAs, the two circRNAs with the highest Base Mean were selected for validation.
The validation was performed on all 70 individuals enrolled. Samples were studied independently, without pooling them. qPCRs were conducted using divergent primers spanning the respective BSJ for each of the 7 circRNA candidates. The correct amplification of the BSJs was confirmed by agarose gel electrophoresis and Sanger sequencing (Supplementary Figure 1). 4 of the circRNAs selected from the RNA-Seq analysis were validated by qPCR (Table 1). The circRNA hsa_circ_0007817 was found to be 2.11 times more expressed in frails when compared to robust individuals (p<0.0001) whereas for hsa_circ_0101802 we validated an upregulation of 1.24 times (p-value=0.03). Remarkably, the two circRNAs selected for being group exclusive could be detected in all the 70 individuals, indicating that their presence/absence cannot be used as a clear criterion for fail/robust classification. Hsa_circ_0060527 and hsa_circ_0075737, were significantly upregulated in frail individuals (FC=1.39; p-value<0.0001 and FC=1.26; p-value=0.005, respectively) (Figure 3 and Table 1).
We also evaluated whether the expression of the 7 candidate circRNAs correlates with the loss of functional abilities. The obtained results show that the expression of hsa_circ_0007817, hsa_circ_0101802, hsa_circ_0058514, hsa_circ_0060527, and hsa_circ_0075737 gradually increase with the time required by the elderly to complete the TUG test. In line with this, as the GS score decreases, the levels of hsa_circ_0007817, hsa_circ_0101802, hsa_circ_0060527 and hsa_circ_0075737 increase. Similarly, the levels of hsa_circ_0007817 and hsa_circ_0060527 increase as the SPPB score decreases. Therefore, a higher expression of the mentioned circRNAs is associated with a worse score in the frailty tests (Supplementary Table 2).
circRNAs are potential biomarkers of frailty
To explore the diagnostic potential of the aforementioned circRNA candidates, we performed a ROC curve analysis. As shown in Figure 4, the 4 validated circRNAs showed a good performance with an area under the curve (AUC) higher than 0.65, indicating that based on the leukocyte expression of any of those circRNAs the probability of correctly classifying them as frail or robust is higher than 65%. Remarkably, the AUC calculated for hsa_circ_0007817 is 0.85 which shows a great potential to function as a frailty biomarker.
Additionally, we investigated whether the use of several circRNAs in combination could improve the potential to correctly discriminate between frails and robust. To do so, we created a Linear Discriminant Analysis to detect the best variables or combination of variables that better discriminate the frailty status. Based on this model, we reported that the combination of the expression levels for hsa_circ_0007817, hsa_circ_0075737 and hsa_circ_0079284 (Combination 1) performs with an AUC of 0.959, indicating that it has a great diagnostic value, with a 95.9% probability of correctly classifying the frail and robust individuals (Figure 4). Moreover, the combinations of the four significant and validated circRNAs (Combination 2: hsa_circ_0007817, hsa_circ_0101802, hsa_circ_0060527, hsa_circ_0075737) and the combination of hsa_circ_0007817 + hsa_circ_0101802 + hsa_circ_0060527 (Combination 3) also outperform the discriminating capacity of individual circRNAs (AUC 0.900 and 0.882, respectively) (Figure 4 and Supplementary Table 3).
The expression of hsa_circ_0079284 is reduced after a physical intervention
Then, we investigated circRNA expression in elder donors after a physical intervention. The expression of the five circRNAs with more promising results regarding their potential to be frailty biomarkers (hsa_circ_0007817, hsa_circ_0101802, hsa_circ_0060527, hsa_circ_0075737, and hsa_circ_0079284) was determined in these samples by qPCR. Only three out of the five selected circRNAs were sufficiently amplified in these samples. No consistent data was obtained for hsa_circ_0007817 and hsa_circ_0060527 due to their low expression.
We found that, after the 3 months of continued exercise, hsa_circ_0079284 levels were reduced in 10 out of the 12 donors with a p-value of 0.02 (Figure 5A). In the case of hsa_circ_0101802 and hsa_circ_0075737, they also show the same reduction trend after the intervention (8/12 and 10/12 donors with reduced levels correspondingly) although these differences were not statistically significant (Figure 5B and C) (Supplementary Table 4). In addition, for the 5 donors from whom samples after 3 months of finishing the intervention were available, we observed a consistent (p-value=0.05) increase in hsa_circ_0079284 levels (Figure 5A) (Supplementary Table 4). When the expression of hsa_circ_0079284 and TUG scores were compared right before and after the intervention, 7 out of the 12 participants obtained concordant results: reduction in hsa_circ_0079284 and TUG performance time or increased hsa_circ_0079284 and TUG performance time. Similarly, for 8 of the 12 patients, there is a concordance between the reduction in hsa_circ_0079284 levels and an increase in the SPPB score (Supplementary Figure 2 and Supplementary Table 4).
Insights into circRNA function: miRNA sponging and secondary structure are not predicted to be leading the way
Furthermore, we attempted to explore the biological function of the circRNAs found to be altered in frailty. The miRNA sponging function was the first described (Hansen et al., 2013) and it is the most studied regulatory mechanism of circRNAs. Increasing evidence supports miRNA binding site (BS) density (number of BS/nucleotide) as a reliable criteria for assessing the miRNA sponging potential (Guo, Agarwal, Guo, & Bartel, 2014). Through the prediction of TargetScan recorded in the CircInteractome database, out of the 15,189 circRNAs detected in our study, 160 circRNAs present more than 20 miRNA BS. Nevertheless, only three of those circRNAs were detected in several samples and present a miRNA BS density with at least 1 miRNA BS/100nts. (hsa_circ_0001946, hsa_circ_0106898 and hsa_circ_0106897), but none of them was found to be significantly differentially expressed with frailty.
On the other hand, recent publications have drawn attention to circRNA secondary structure as a determinant feature that can define their interaction with different proteins, and thus, their function (Fischer, Busa, Shao, & Leung, 2020; Liu et al., 2019). It has been reported that circRNAs can form intramolecular double-stranded regions resulting in highly structured circRNAs, which preferentially interact with proteins such as the double-stranded RNA activated protein kinase R (PKR) (Liu et al., 2019) or the UPF1 RNA binding protein and its associated protein G3BP1 (Fischer et al., 2020). Among all the 15,189 circRNAs detected in our samples, 24.7% are predicted to be highly structured (according to Fischer et al), but no differences could be observed linked to frailty status.