Characteristics of circRNA expression pattern in stenosed AVF tissues
circRNA sequencing to explore the profile was initiated on Illumina Hiseq platform by analyzing four stenosed AVF tissues and four control vascular tissues from the ESRD patients. 17,620 circRNAs were identified, and their distribution characteristics on chromosomes are illustrated in Supplementary Fig. 1A. These circRNAs were transcribed from all human chromosomes including 22 autosomes and 2 gonosomes, most of which were distributed on chromosome 1, chromosome 2, chromosome 4, and chromosome 6. The least distribution is on chromosome Y. Moreover, the top 10 upregulated or downregulated circRNAs were mainly scattered on chromosome 1. Based on the properties of parental genes from which circRNAs are derived, the high-quality circRNAs identified by find_circ were categorized as annot_exons, antisense, intergenic, intron_exon, intronic, and one_exon. The proportion of annot_exons was found to be the largest (Supplementary Fig. 1B).
Identification of differentially expressed circRNAs between AG and CG
The correlation of circRNA expression pattern between samples was tested before the DE circRNA analysis, and we show the correlation coefficients in Fig. 1 (the coefficient is closer to 1, the expression pattern between samples is more similar). Hierarchical cluster analysis was used to display the expression levels of circRNAs in the two group samples. The heatmap exhibited that there was significant difference in circRNA expression in the two groups (Fig. 2). Based on the filter criteria of fold change > 2 or < -2 and P value < 0.05, we identified 208 DE circRNAs between AG and CG. Among these, 92 DE circRNAs were significantly over-expressed and 116 circRNAs decreased more than two-fold in AG comparing to CG. Volcano plot authenticated the fold change of circRNAs after log2 transformations were clearly differentiated, which showed that the distribution was roughly symmetrical in the two groups (Supplementary Fig. 2). Supplementary Table 3 lists the top 10 DE circRNAs that were significantly upregulated or downregulated.
Functional enrichment analysis of DE circRNA parental genes
We observed that the parental genes of 208 DE circRNAs were assigned to 1,959 GO terms, which included 1,450 GO terms enriched for BP, 269 GO terms enriched for MF, and 240 GO terms enriched for CC. In Fig. 3, we see that there are 6 GO terms assessable with P value < 0.05. The majority of the hosting protein-coding genes of DE circRNAs were significantly enriched in BP such as single-multicellular organism process, anatomical structure development, and multicellular organism development. For MF, there are four parental genes related to alpha-actinin binding (P value < 0.05). Other genes may be associated with metal ion binding and cytoskeletal protein binding though the results were not significant (P value > 0.05). Some genes enriched in CC correlated with voltage-gated calcium channel complex, cell leading edge, and basal lamina. Unfortunately, there was no significant correlation (P value > 0.05).
The top 20 KEGG pathways that significantly enriched are presented in Supplementary Fig. 3, which include focal adhesion (FA), regulation of actin cytoskeleton, vascular smooth muscle contraction, and ErbB signaling pathway. FA was the most significantly enriched pathway, and most of DE circRNAs were involved in the FA also (Transcript count = 11). Results indicated that DE circRNAs may be involved in AVF stenosis development and progression through interference of the motility and metabolism of VSMC from the AVF outflow vein.
Validation of candidate circRNAs in AVF stenotic tissue
Seven candidate circRNAs (downregulated: hsa_circ_015036, hsa_circ_005718, hsa_circ_017814, hsa_circ_019184, upregulated: hsa_circ_007293, hsa_circ_026394, and hsa_circ_035577) were filtered out to validate the reliability of previous RNA-sequencing results, and their expression levels between AG and CG were detected using qRT-PCR. The candidate circRNAs satisfied following conditions: (1) fold change > 2 or < -2; (2) P value < 0.05; (3) the host protein-coding genes of candidate circRNAs were involved in FA signal pathway, which may be highly correlated with AVF dysfunction. Based on qRT-PCR results, the hsa_circ_026394 and hsa_circ_035577 expression trends in the two groups were consistent with the RNA-sequencing results, but there was no significant statistical difference (P value > 0.05). Moreover, hsa_circ_019184 expression abundance in the two samples was too low to be detected. Therefore hsa_circ_019184, hsa_circ_026394 and hsa_circ_035577 were excluded for further study. The hsa_circ_007293 expression levels in AG were significantly higher in CG, and hsa_circ_015036, hsa_circ_005718, and hsa_circ_017814 were significantly under-expressed in stenosed AVF tissues comparing to normal vascular tissues, which mostly followed the circRNA transcript expression trends that detected by RNA-sequencing (Fig. 4).
Establishing interaction network between target miRNAs and candidate circRNAs
We further predicted miRNA binding sites on four candidate circRNAs (hsa_circ_015036, hsa_circ_005718, hsa_circ_017814, and hsa_circ_007293) through miRanda software, which included three downregulated circRNAs and one upregulated circRNA. Based on the criteria of max score ≥ 140 and max energy ≤ -25 [14], a total of 134 miRNA were paired with four DE circRNAs after discarding the repeated values. For each candidate circRNA, top 10 target miRNAs with the highest degree of correlation were considered (a lower max energy is indicative of a stronger correlation) [14], then we acquired 38 miRNAs after discarding repeating target miRNAs. Figure 5 illustrates the miRNA-circRNA interaction network regarding the four candidate circRNAs.