Differential Circular RNA Expression Analysis for Screening the Potential Biomarkers in Atrial Fibrillation Patients

Shenjie Sun Tsinghua University Tingting Lv Beijing Tsinghua Changgung Hospital Siyuan Li Beijing Tsinghua Changgung Hospital Peng Liu Tsinghua University Ying Yang Tsinghua University Yuanwei Liu Beijing Tsinghua Changgung Hospital Fei She Beijing Tsinghua Changgung Hospital Rong He Beijing Tsinghua Changgung Hospital Ping Zhang (  zhpdoc@126.com ) Beijing Tsinghua Changgung Hospital https://orcid.org/0000-0002-5843-1108


Background
Atrial brillation (AF) is the most common arrhythmia worldwide [1], and its global prevalence is between 2% and 4% [2,3]. As a kind of chronic cardiac disease, AF seriously damages health by increasing the risks of stroke, heart failure, and mortality [4]. At present, rhythm or rate control and stroke prevention are major strategies for treating AF [5,6]. Exploring the etiology of AF is crucial in searching for improved treatment.
As a stable circular transcript, circular RNA (circRNA) is prior to become a kind of potential biomarker for the diagnosis and therapy of disease [7]. One previous study reported that circRNA exacerbated cardiac hypertrophy [6], increased cardiomyocyte proliferation [8], and suppressed myocardial brosis [9].
To clarify the potential role of circRNA/miRNA/mRNA network in AF, we adopted a more stringent adjusted P value and multi-omics analysis method to explore the complex interactions of circRNA/miRNA/mRNA regulatory network in AF. Finally, we screened for differentially expressed genes (DEGs) and differentially expressed circRNAs (DE circRNAs) in patients with AF and those with a sinus rhythm (SR) from public circRNA and mRNA expression data in the Gene Expression Omnibus (GEO) database. The DE circRNAs and DEGs were identi ed between the groups by bioinformatic analysis. The circRNA/miRNA/mRNA regulatory network was constructed to identify the key circRNAs associated with AF that could be further explored as candidate biomarkers.

Results
Differential expression analysis of circRNAs and mRNAs in AF 13,295 probes and 10,029 circRNAs in GSE97455 were analyzed in the present study. If the t-test adjusted P-value was less than 0.05 and the fold change was more than 2, further analysis using the limma package revealed 46 DE circRNAs, including 29 downregulated and 17 upregulated ones, which are shown in detail in Supplementary Table S1. The DE circRNA expression levels are presented in heat maps drawn by the limma package of R software (Fig. 1A). The DE circRNA distribution based on logFC and adjusted P value is illustrated in a volcano map drawn using the pheatmap package of R software ( Fig. 2A) Table S4). The results of the top 20 GO analysis were illustrated in Fig. 2C, ordered by GeneRatio and P value. The top ve GO terms by count were the biological process terms, including muscle system process, leukocyte migration, regulation of in ammatory response, regulation of lymphocyte activation, and positive regulation of response to external stimulus. The cellular component terms were external side of plasma membrane, collagencontaining extracellular matrix, and MHC class II protein complex. In addition, the top ve KEGG terms were cytokine-cytokine receptor interaction, in uenza A, hematopoietic cell lineage, leishmaniasis, and in ammatory bowel disease (Fig. 2D) (Supplementary Table S5).
Construction of circRNA/miRNA/mRNA regulatory network in AF We predicted 298 miRNAs targeted by 46 DE circRNAs and obtained 1,181 circRNA-miRNA pairs using Circular RNA Interactome online tools (Supplementary Table S6-7). Subsequently, we predicted that 13,821 genes appeared in any of at least two prediction tools of miRDB, miRTarBase, and TargetScan and obtained 68,526 miRNA-mRNA pairs, including 189 putative miRNAs (Supplementary Table S8-9).
Thirty-two overlapping genes were obtained by taking the intersection between 13,821 putative genes and the top 41 hub genes in GSE79768 (Supplementary Table S10). From this we constructed the circRNA/miRNA/mRNA regulatory network including 42 circRNAs, 90 miRNAs, and 32 mRNAs forming 376 circRNA-miRNA pairs and 134 miRNA-mRNA pairs in the network (  The relationship among circRNA, miRNA, and mRNA in AF The GO and KEGG functional enrichment analysis results were used to predict the potential functions using the circRNA/miRNA/mRNA regulatory network in AF. In the potential AF-associated signaling pathway, GO analysis showed that the calcium ion transport-related pathway and the second-messengermediated signaling pathway contained the hub genes ADRA2A (adrenoceptor alpha 2A), CXCL11 (C-X-C motif chemokine ligand 11), P2RY12 (purinergic receptor P2Y12), and PTGS2 (prostaglandinendoperoxide synthase 2); the programmed cell death pathway contained the hub genes IL1B, KIT (tyrosine-protein kinase Kit), and KITLG (KIT ligand); and the renin-angiotensin pathway included AGTR2 (angiotensin II receptor type 2). The potential relationship among the hub genes, the circRNA/miRNA/mRNA regulatory network, and functional enrichment analysis in AF (Supplementary  Table S20).

Discussion
In recent years, although cardiologists have made great progress for the diagnosis and treatment of AF, the high incidence and high disability rate induced by AF still constitute a heavy medical burden globally. And some points regarding the pathogenesis of AF are still unclear. In this study, we constructed a regulatory network of 7 circRNAs, 42 miRNAs, and 25 mRNAs using GEO datasets to explore the potential circRNA involved in AF. Combining the results of enrichment analysis, we determined the occurrence mechanism in calcium signaling-related pathway, second messenger-mediated signaling-related pathway, programmed cell death pathway, and renin-angiotensin pathway (Supplementary Table S20).
There are four hub genes, namely, ADRA2A, P2RY12, CXCL11, and PTGS2, involved in the regulation of calcium ion transport pathway and second messenger-mediated signaling pathway. The hsa_circ_0001666/hsa_miR-1827/ADRA2A pathway involved in AF through the second messenger signaling pathway. Previous study also showed that lacking ADRA2A elevated plasma noradrenaline concentrations and participated in the development of cardiac hypertrophy and dysfunction [11]. In addition, hsa_circ_0000288, hsa_circ_0021652/hsa_miR-548m/P2RY12 network also played an important role in AF through a similar mechanism. The latest research illustrated that the crosstalk between P2Y12 and IGF-I receptors can increase PKB phosphorylation through calcium-dependent activation of the Pyk2/Src pathway [12]. The hub gene PTGS2 involved in more complex circRNA/miRNA/mRNA regulatory networks in AF including 5 circRNA, such as hsa_circ_0001666, hsa_circ_0006168, hsa_circ_0021652, hsa_circ_0000288, and hsa_circ_0030162. Previous research also demonstrated that the combination of PTGS2 and miRNA-26 can inhibit the MAPK pathway and reduce the in ammatory response and myocardial remodeling of mice with myocardial infarction [13].
Moreover, there are three hub genes, namely, IL1B, KIT, and KITLG, involved in the programmed cell death pathway. Our research data revealed that the hsa_circ_0001666/hsa_miR-590-5p/IL1B involved in AF via programmed cell death pathway. And hsa_circ_0021652, hsa_circ_0000288/hsa_miR-583/KIT networks also participated in AF via similar pathway. Previous research results suggested that IL1B genetic variation is involved in mediating CRP expression [14]. In reviewing the literature, no data was found about the association between CXCL11, KIT, and KITLG in cardiac diseases. The hsa_circ_0001666/miR-1827/ADRA2A pathway, the hsa_circ_0000288 or hsa_circ_0021652/miR-548m/P2RY12 pathway, circ_0001666/miR-590-5p/IL1B pathway, will become the target regulation network of our research.
Furthermore, a total of four miRNAs, namely, miR-1263, miR-1299, miR-1322, and miR-330-3p, were involved in the renin-angiotensin system and second messenger-mediated signaling pathway with the AGTR2 gene. Previous studies revealed that AGTR2 is involved in regulating collagen accumulation and MMP expression [15]. The network showed that hsa_circ_0001666 can adsorb three miRNAs, namely, miR-1299, miR-1322, and miR-330-3p, while hsa_circ_0000367 interacts with miR-1263. The other three circRNAs, hsa_circ_0030162, hsa_circ_0021652, and hsa_circ_0000288, can also interact with miR-1322 to in uence AGTR2 expression in AF. The regulatory network involved in AGTR2 can become the most valuable target in our future research work.
The main limitation of this study is the small sample size, which may have resulted in selection bias. Therefore, further genomics research on more specimen sources is required to validate these results. We are currently preparing to perform knockdown and overexpression experiments to explore the effects of changes in the expression of circRNAs and associated key genes using cell and animal models. Despite these limitations, the ndings of this study suggest that circRNAs are also potential biomarkers or therapeutic targets for AF, offering a guide for further research to better understand the mechanism of AF.

Conclusions
In conclusion, we constructed a circRNA/miRNA/mRNA regulatory network including seven highly expressed circRNAs (hsa_circ_0000367, hsa_circ_0000691, hsa_circ_0000288, hsa_circ_0006168, hsa_circ_0001666, hsa_circ_0021652, and hsa_circ_0030162) through GEO datasets and performed GO and KEGG analyses to explore the potential mechanisms of action of circRNAs in AF. These circRNAs are worthy of further exploration for use as novel biomarkers in AF.

Data Retrieval
We downloaded original les on circRNA and mRNA expression from the GEO database (http://www.ncbi.nlm.nih.gov/geo/). The circRNA expression data were derived from the GSE97455 datasets, including 15 pairs of plasma specimens of AF and a control group. The GPL21825 platforms were generated using Arraystar Human CircRNA microarray V2. The mRNA expression data were obtained from GSE79768 datasets including seven left atrial specimens with AF and six specimens with SR. The miRNAs targeted by circRNAs were predicted using Circular RNA Interactome online tools (https://circinteractome.nia.nih.gov/) [16].

Differential expression analysis of circRNAs and mRNAs in AF
The probe name was converted to the circRNA international standard symbol by using the GPL21825 annotation le in the GSE97455 series matrix with Perl software 5.30.1. Differentially expressed circular RNAs (DE circRNAs) were analyzed using the limma package of R software. DE circRNAs were de ned based on the adjusted p-value was less than 0.05 and the fold change was more than 2. The CEL les of GSE79768 were cleaned including background correction, data normalization, and expression level calculation to form a gene expression level matrix. We converted the probe name to the gene symbol using the GPL570 platform annotation le by R software. DEGs were identi ed according to the adjusted p-value was less than 0.05 and the fold change was more than 2 using the limma package. The missing values were lled with the weighted k-nearest-neighbor technique [17]. The P value was adjusted by the Benjamini-Hochberg method.

PPI and functional enrichment analysis of mRNAs in AF
We constructed the protein-protein interaction (PPI) network to obtain the interactions of the identi ed DEGs in the GSE79768 by multiple protein function tools in the STRING database (STRING database, Version 11, https://string-db.org) [18]. The con dence limit of the interaction score median was set to 0.40. Cytoscape 3.7.2 was used to draw the network of the DEGs and the cytoHubba plugin components were to identify the hub genes in the PPI network sorted by the degree ranking table (https://cytoscape.org/) [19]. We performed functional enrichment analysis on DEGs in GSE79768 by use of the functional annotation tool DAVID 6.8 (https://david.ncifcrf.gov/) [20]and the R clusterPro ler package.
Construction of the circRNA/miRNA/mRNA regulatory network in AF We predicted miRNAs and established circRNA-miRNA pairs using the DE circRNAs in Circular RNA Interactome. Putative miRNAs targeted by DE circRNAs were used to predict mRNAs by utilizing miRDB [21], miRTarBase [22] and TargetScan [23]. The intersection of predicted mRNAs and above-mentioned hub genes from GSE79768 were used to establish miRNA-mRNA pairs. Then, the circRNA/miRNA/mRNA regulatory network was constructed using Cytoscape 3.7.2.

Declarations
Ethics approval and consent to participate Not necessary. The raw data of this study are derived from the GEO data portal, which are publicly available databases.