Bioinformatics Analysis of Competing Endogenous RNA Network and Immune Infiltration in Atrial Fibrillation

Background There is still no clear understanding of the pathogenesis of atrial fibrillation (AF). For this purpose, we used integrated analysis to uncover immune infiltration characteristics and investigated their relationship with competing endogenous RNA (ceRNA) network in AF. Methods Three AF mRNA data sets (GSE14975, GSE79768, and GSE41177) were integrated using the SVA method from Gene Expression Omnibus (GEO). Together with AF circRNA data set (GSE129409) and miRNA data set (GSE70887) from GEO database, we built a ceRNA network. Then hub genes were screened by the Cytoscape plug-in cytoHubba from a protein-protein interaction (PPI) network. As well, CIBERSORT was employed to investigate immune infiltration, followed by Pearson correlation coefficients to unravel the correlation between AF-related infiltrating immune cells and hub genes. Ulteriorly, circRNA-miRNA-mRNA regulatory axises that could be immunologically related to AF were obtained. Results Ten hub genes were identified from the constructing PPI network. The immune infiltration analysis revealed that the number of monocytes and neutrophils was higher, as well as the number of dendritic cells activated and T cells regulatory (Tregs) was lower in AF. Seven hub genes (C5AR1, CXCR4, HCK, LAPTM5, MPEG1, TLR8, and TNFSF13B) were associated with those 4 immune cells (P < 0.05). We found that the circ_0005299–miR-1246–C5AR1 and circRNA_0079284-miR-623-HCK/CXCR4 regulatory axises may be associated with the immune mechanism of AF. Conclusion The findings of our study provide insights into immuno-related ceRNA networks as potential molecular regulators of AF progression.


Introduction
Atrial fibrillation (AF) is a common arrhythmia, increasing with age, reaching 7.5% or more in people older than 80 [1]. In addition, the risk of embolic stroke, heart failure, and mortality increase with AF [1,2]. Although treatment strategies of AF have advanced dramatically in recent years, their efficacy is not ideal, especially in the radiofrequency ablation treatment of persistent AF (PAF) [3]. e reason is the incomplete knowledge of the AF mechanisms. To develop more effective treatments for AF, we therefore need to gain a deeper understanding of the molecular and cellular mechanisms involved in AF.
CircRNAs are new RNA molecules with unique biological functions that can act as sponges for microRNAs (miRNAs) to bind competitively and regulate parental genes expression [4,5]. A growing body of evidence suggests that circRNA-miRNA-mRNA regulatory axis are involved in cardiovascular disease's pathogenesis. CiRS-7, for example, is useful as a sponge for miRNA-7a, which promotes myocardial apoptosis via inhibiting PARP and SP1 expression [6]. e competitive binding of heart-related circRNA (HRCR) with endogenous miR-223 increases the expression of ARC gene, thereby inhibiting heart failure and cardiac hypertrophy [7]. Further, circRNA_000203 promote fibrosis-associated gene expression by inhibiting miR-26b-5p targets, contributing to myocardial fibrosis [8]. Hence, competing endogenous RNA (ceRNA) networks may shed light on AF pathophysiology. Studies have also shown that inflammation and the immune response it triggers are crucial to the development of AF [9,10]. However, rarely have studies examined the relationship between immune cells infiltration and ceRNA networks in atrial tissue of AF patients.
In our study, on the one hand, we integrated three PAF data sets by the SVA method and used weighted gene coexpression network analysis (WGCNA) and differential expression analysis to identify common genes (CGs). Differentially expressed (DE) miRNAs and DE circRNAs in PAF were identified from GSE70887 and GSE129409 data set, respectively. Together with CGs, ceRNA network was built based on circRNA-miRNA pairs and miRNA-mRNA pairs. Hub genes were then filtered using the Cytoscape plug-in cytoHubba by analyzing protein-protein interaction (PPI) networks. On another hand, CIBERSORT was used to study immune infiltration in AF [11]. In the following analysis, Pearson correlation coefficients were used to determine the correlation between AF-related infiltrating immune cells and hub genes. Finally, we gained novel insight into the mechanisms that govern the progression of AF by the analysis of immune-related ceRNA networks. e study flowchart is shown in Figure 1.

Data Acquisition.
ree AF mRNA data sets (GSE14975, GSE79768, and GSE41177), AF circRNA data set (GSE129404), and AF miRNA data set (GSE70887) were downloaded from gene expression omnibus (GEO) [12] database. ree left atrial appendage samples from PAF patients and three sinus rhythm (SR) controls were included in GSE129409, while four atrial appendage samples from PAF patients and two SR controls were contained in GSE70887. Among the GSE14975, GSE79768, and GSE41177, left atrial appendage tissue was obtained from five PAF patients and five SR controls, seven PAF patients and six SR controls, and sixteen PAF patients and three SR controls, respectively.
All data sets originated from a free open-access database on the Internet; thus, this study does not require ethical approval and patient consent. (GSEA). In order to transform gene probe IDs to gene symbol codes, the series matrix files were processed by ActivePerl 5.24.2 software (https://www.activestate.com/ products/perl/). In the three mRNA data sets, the data in GSE14975 data set were transformed into log base 2 data by affy package in R [13]. e sva package's combat function was applied to remove batch effects and other undesired variation between the three mRNA microarray data sets after merging all data [14]. In the end, R software's "limma" package contains the "normalizeBetweenArrays" function [15], which normalized expression values. Gene set enrichment analysis (GSEA) is a computational algorithm for determining whether a predefined set of genes exhibit consistently significant differences between two states [16]. In a GSEA, the sequenced genes of AF and SR samples are   analyzed after imputting gene annotation files, reference function sets, and all the gene data from both samples. e pathways enriched in each phenotype were analyzed based on nominal p value and normalized enrichment score (NES).

Weighted Gene Co-Expression Network Analysis (WGCNA) Construction and Identification of Modules.
Gene co-expression network was constructed using the integrated data set with the help of a system biology approach of WGCNA [17]. e soft thresholding power β was set as 5 and 20 and was selected using the function pick-Soft reshold. In order to classify genes with similar expression profiles into gene modules, average linkage hierarchies were clustered according to topological overlap matrix (TOM)-based difference measure, and the minimum size (gene group) of the genes dendrogram was 50 [18]. Finally, module membership (MM), gene significance (GS), and module-trait correlations analyses were conducted. P < 0.05 was defined as statistically significant module.   Common genes (CGs) are the intersection of DE mRNA identified from the integrated mRNA data set and the genes found in yellow module.

Construction of a circRNA-miRNA-mRNA Regulatory
Network. CircRNAs information can be found in CircBase (https://www.circb ase.org/) [19]. e cancer-specific circRNA database (CSCD, https://gb.whu.edu.cn/CSCD/) [20] was able to predict target miRNAs for each DE circRNA. en, we gathered miRNAs that overlapped both DE and predicted miRNAs and used TargetScan database [21] to predict targeted genes. Next, those targeted genes were considered as candidate targets and overlapped with CGs. Lastly, we constructed a ceRNA regulatory network of AF and visualized it using Cytoscape version 3.8.0.

Construction of PPI Regulatory Network and Identification of Hub Genes.
Using the STRING datebase (https:// string-db.org) [22], a PPI network was constructed for these genes in the ceRNA network, and a minimum interaction score of 0.4 was considered the cutoff point. CytoCope 3.8.0 software was utilized to visualize the PPI network, and the Maximal Clique Centrality (MCC) arithmetic of the Cytoscape plug-in cytoHubba was used to filter hub genes in the PPI network. Finally, a circRNA-miRNA-hub gene subnetwork was constructed. Boxplot maps representing differential expression of circRNAs, miRNAs, and mRNAs in their microarray data sets in the circRNA-miRNA-hub gene subnetwork were generated with the help of "reshape2" and "ggpubr" packages of R software.

Immune Cell Infiltration
Analysis. e integrated data set was analyzed using CIBERSORT in R software to compute the relative proportion of infiltrating immune cells in AF, and the samples were filtered using P < 0.05. A principal component analysis (PCA) was performed on immune cell infiltration using the "ggplot2" package. e 22 types of infiltrating immune cells were subjected to a Spearman correlation analysis using R software, and we generated a correlation heatmap with the help of "Corrplot" package in R for visualizing the results. e expression of 22 immune cells was compared and visualized using "vioplot" package between PAF and SR samples.

Correlation Analysis between Hub Genes and Infiltrating
Immune Cells Associated with AF. A Pearson correlation coefficient was applied to examine the relationship between hub genes and AF-related infiltrating immune cells, which was visualized with the "ggpubr" package of R.

Diagnostic Analysis of Hub Immune-Related Genes for AF.
For the purpose of determining the effectiveness of hub immune-related genes in predicting AF, receiver operator characteristic (ROC) curve analysis was conducted with the help of "pROC" package.

Enrichment Analysis of Merged Expression Data through
GSEA. GSEA was applied to analyze the significant difference between AF and SR groups for the integrated data set. e enrichments for upregulated gene sets in the significant order (size of NES) were related to physiological cardiac muscle hypertrophy (GO) (Figure 2

Identification of Gene Co-Expression Networks and
Modules. Using the WGCNA package, gene co-expression network was built from the integrated gene data set to identify functional clusters in AF patients. Eight modules were excavated after setting the power to 5 (Figures 3(a)). In Figure 3(b), the module-trait relationships are illustrated, showing that the yellow module has the greatest relationship with AF (r � 0.47, p � 0.002), encompassing 365 genes. e greatly significant correlation between GS and MM indicates that genes in the yellow module are greatly associated with AF (cor � 0.5, p � 1.7e − 24) (Figure 3(c)).

Identification of DE circRNAs, DE miRNAs, DE mRNAs, and CGs.
In the merged mRNA data set, a total of 439 DE mRNAs were screened in AF (Figures 4(a), 4(b)). en, CGs were defined as the intersection of DE mRNAs from the integrated mRNA data set and genes from yellow module (Figure 4(c)), including 110 genes. Totally, 103 DE circRNAs

Construction of ceRNAs Regulatory Networks in AF.
Fifteen DE circRNAs were not found in the CSCD database. Based on this database, 1,995 targeted miRNAs were predicted from the remaining 88 DE circRNAs. In the next step, 9 miRNAs were obtained through the intersection of DE miRNAs and predicted miRNAs ( Figure 5(e)). Using the TargetScan database, 9 miRNAs predicted 13,135 potential targets. en, 64 mRNAs were acquired by the intersection of CGs and predicted target genes ( Figure 5(f )). Finally, the ceRNA network associated with AF was constructed ( Figure 6(a)).

Functional Enrichment Analyses for mRNAs in the ceRNAs
Network. Go functional enrichment analysis revealed that those genes in the circRNA-miRNA-mRNA ceRNA network were primarily involved in biological process (BP) terms, including "T cell activation," "lymphocyte proliferation." In the cell component (CC) ontology, those genes were mainly enriched in "external side of plasma membrance," "endocytic vesicle." Molecular function (MF) analysis indicated that those genes were significantly enriched in "immune receptor activity," "coreceptor activity" (Figure 6(b)). e KEGG pathway of those genes were primarily involved in "hematopoietic cell lineage," "chemokine signaling pathway," "leukocyte transendothelial migration" (Figure 6(c)).
3.6. PPI Network Analysis. Using STRING database, 55 mRNAs in ceRNAs network were constructed into a PPI network consisting of 45 nodes and 160 edges after removing unconnected nodes (Figure 6(d)). To explore and construct the crucial circRNA-miRNA-hub genes regulatory axis in the progression of AF, the MCC algorithm was used to identify hub genes in the PPI network. rough the MCC scoring method, the top ten genes were defined as hub genes (Table 1), which were C-C chemokine receptor type 5 (CCR5), C-X-C chemokine receptor type 4 (CXCR4), Tolllike receptor 8 (TLR8), stromal cell-derived factor 1 (CXCL12), C5a anaphylatoxin chemotactic receptor 1 (C5AR1), hematopoietic cell kinase (HCK), tumor necrosis factor ligand superfamily member 13B (TNFSF13B), interferon regulatory factor 8 (IRF8), macrophage gene 1 protein (MPEG1), and lysosomal-associated transmembrane protein 5 (LAPTM5), respectively ( Figure 6(d)). And those hub genes were upregulated in AF. Subsequently, a circRNA-miRNA-hub gene subnetwork was showed in Figure 6(e). Figure 7 shows the differential expression of each molecule in the ceRNA subnetwork in the microarry data sets. e basic information of the 9 circRNAs in the ceRNA subnetwork are listed in Table 2.

Immune Infiltration Analyses.
Using the CIBERSORT algorithm, bar plots and heat maps display the relative proportion of the 22 types of immune cells detected in each sample (Figures 8(a), 8(b)). Using PCA, the immune cells in the atrial tissues of AF patients and SR revealed distinct group bias clustering and individual differences (Figure 8(c)). Correlation analysis between infiltrating immune cells and T cells regulatory (Tregs) was positively related to macrophages M0 (r � 0.69) and negatively correlated to neutrophils (r � 0.40) and the activated dendritic cells were positively related to activated NK cells activated (r � 0.54). In contrast, T cells CD8 were negatively related to T cells CD4 memory resting (r � −0.67), mast cells resting were negatively related to activated mast cells (r � −0.58), and T cells gamma delta negatively correlated to macrophages M2 (r � −0.52) (Figure 8(d)). e monocyte and neutrophil counts tend to be higher in AF than SR, while dendritic cells activated and T cells regulatory (Tregs) tend to be lower (P < 0.05) (Figure 8(e)).

Discussion
AF is the most frequent arrhythmia that presents in clinical practice and increases the risk for heart failure, stroke, and death. By further understanding the mechanisms, we may find new strategies to treat AF. We performed a combined analysis of ceRNA networks and immune infiltration associated with AF to explore the molecular mechanism. First, we performed GSEA to analyze the significant difference between AF and SR groups for the integrated data set. Gene sets related to cardiac muscle hypertrophy, collagen binding, innate immune response, chemokine signaling pathway, renin angiotensin system, and T cell receptor signaling pathway were differentially enriched with AF phenotype. e chemokine signaling pathway plays a key role in cardiovascular disease. For example, chemokines and their receptors are important in the recruitment and activation of immune cells and the persistence of the local inflammatory response in atherosclerosis [23]. Chemokines and their receptors have also been shown to be involved in the pathophysiology of cardiac remodeling and heart failure resulting from excessive pressure load [24]. e renin angiotensin system is known to be closely associated with mechanism of AF development. e atrial electrical and structural remodeling are the core part of AF, and the activation of the renin-angiotensin-aldosterone system (RAAS) contributes to atrial remodeling [25]. Conversely, studies have demonstrated that renin-angiotensin system inhibitors can delay heart remodeling in patients with AF [26] and prevent recurrence of AF after ablation [27]. So, GSEA indicated that those gene sets in AF patients were primarily involved in cardiac muscle, immune, and inflammatory responses.
Next, we performed two different methods (WGCNA and LIMMA method) to screen CGs related to AF from the integrated data set. Together with DE circRNAs and DE miRNAs related to AF, the ceRNA network was constructed, suggesting that these mRNAs, miRNAs, and circRNAs in the circRNA-miRNA-mRNA ceRNA network could play an important role in the pathogenesis of AF. GO functional enrichment analysis revealed that these genes in the ceRNA network were significantly involved in the regulation of lymphocyte and mononuclear cell proliferation and immune receptor activity. KEGG pathway analysis results were predominantly enriched in "Chemokine signaling pathway" and "Cytokine-cytokine receptor interaction." According to GO and KEGG results, ceRNA network about AF in this study was also mainly correlated with inflammation and immunity. Subsequently, we constructed a PPI network based on 55 mRNAs in the ceRNA network and used the MCC algorithm in the CytoHubba plug-in to filter 10 hub genes (CCR5, CXCR4, TLR8, CXCL12, C5AR1, HCK, TNFSF13B, IRF8, MPEG1, and LAPTM5). Among them, the highest scored gene was CCR5. Research has shown that CCR5 is involved in autoimmune and inflammatory diseases such as rheumatoid arthritis [28] by regulating the activation and migration of immune cells. And studies also have indicated that CCR5 may play a role in Ang II-induced hypertension and vascular dysfunction [29], as well as in the development of arthrosclerosis and cardiovascular disease [30]. Additionally, CCR5 inhibition protects against pressure overload-induced cardiac dysfunction through P38 and ERK1/2 signaling pathways [31]. erefore, we speculated that CCR5 may be implicated in the pathogenesis of AF through immune and inflammatory responses, which needs to be further explored.
en, CIBERSORT was used to study immune infiltration in AF. According to the study, we found that an increase in monocytes and neutrophils, along with a decrease in dendritic cells and regulatory T cells (Tregs), might be linked to AF pathogenesis. And our bioinformatics analysis showed 7 hub genes (C5AR1, CXCR4, HCK,       LAPTM5, MPEG1, TLR8, and TNFSF13B) were associated with those 4 AF-related immune cells, and these genes may be able to predict AF based on the ROC curve analysis (0.7 < AUC < 1). So, we conjectured that these 7 genes may be involved in AF pathogenesis by regulating these 4 immune cells.

B cells naive B cells memory Plasma cells T cells CD8 T cells CD4 naive T cells CD4 memory resting T cells CD4 memory activated T cells follicular helper T cells regulatory (Tregs) T cells gamma delta NK cells resting NK cells activated
Ulteriorly, we performed a combined analysis of circRNA-miRNA-hub gene subnetwork and immune infiltration associated with AF. It has been found that inflammatory infiltration of myocardium, including neutrophils and inflammation markers, may contribute to AF [32]. Neutrophils dominate the inflammatory cells in AF patients who undergo pericardiotomy, atriotomy, or catheter ablation according to some studies [33,34]. Furthermore, elevated neutrophil-to-lymphocyte ratio (NLR) has been shown to be independently associated with a higher risk of all-cause mortality and combined end point events in patients with AF and to be an independent predictor of longterm prognosis in AF patients [35]. Our study also found neutrophils were upregulated in AF tissues. We found that T cell regulatory (Tregs) proportion in atrial tissues of AF patients was significantly lower, which was similar to the results of Chen et al. [36]. And they validated IL-6-miR-210 inhibits Tregs function by targeting Foxp3 to promote atrial fibrosis. Recently, we downloaded immune genes from the database and constructed the immune cell-related ceRNA subnetwork through bioinformatics analysis, which the results showed that Tregs were also underexpressed in atrial auricular tissue of AF [37]. One more study revealed that Tregs alleviate myocardial fibrosis and cardiac hypertrophy in hypertensive mice caused by angiotensin II [38]. is suggests that neutrophils and Tregs might be important core cell subtypes in driving AF disease progression. Furthermore, our finding showed that C5AR1 was highly expressed in the auricle tissue of AF patients and was positively associated with neutrophils and negatively correlated with Tregs. Studies confirmed that high levels of C5a and its interaction with C5aR1 led to excessive activation of central neutrophil functions [39,40]. Meanwhile, C5AR1 was involved in cardiovascular disease. C5aR1 activation has been reported within atherosclerotic plaques [41,42], and C5AR1-deficient mice revealed obviously mitigate cardiac remodeling and inflammation after Ang II infusion [43]. From our constructed circRNA-miRNA-hub gene subnetwork, two downregulated miRNAs (miR-1246 and miR-483-5p) and their upregulated C5AR1 target gene were regulated by four upstream upregulated circRNAs. Study found that miR-1246 in endothelial progenitor cell (EPC)derived exosomes enhanced in vitro and in vivo angiogenesis in myocardial infarction (MI), and these improvements may be involved in the reduction of myocardial injury and cardiac fibrosis after MI [44]. erefore, we hypothesized that low expression of miR-1246 in patients with AF may ultimately lead to atrial fibrosis by targeting and enhancing fibrosis-related gene expression. ese reports suggested that circ_0005299-miR-1246-C5AR1 regulatory axis could be associated with the immune mechanism of AF.
Previous study has shown that dendritic cells are present in damaged heart tissue and play a significant role in cardiac remodeling after MI [45]. However, we found that activated dendritic cells were downregulated (P < 0.05), and resting   dendritic cells tended to be upregulated (P > 0.05) in the left atrial tissue of patients with AF. We speculated that different timepoints of the AF onset may determine the contrary results.
In AF patients, the proportion of intermediate monocytes was higher compared with the control group [46], especially in those with low-voltage zones [47]. Furthermore, the activation of monocytes, more specifically their enhanced migration ability, is crucial in the pathogenesis of atrial remodeling in AF patients [48,49]. Our analysis also found that monocytes were highly expressed in atrial tissues of AF patients, and GO enrichment analysis of genes in the ceRNA network was mainly enriched in regulating the proliferation of monocytes, which was positively correlated with HCK gene. Study showed that HCK gene expression was increased in LPS-stimulated human peripheral blood monocyte-derived macrophages [50]. And HCK as critical for regulating alternative activation of monocytes [51]. Our study found that Tregs were also negatively associated with other 5 hub genes (CXCR4, LAPTM5, MPEG1, TLR8, and TNFSF13B) in AF tissues. It has been reported that CXCR4 and LAPTM5 are involved in the negative regulation of Tregs [52,53]. Moreover, Wang et al. [54] showed that CXCR4 is overexpressed in AF patients, which may lead to the occurrence of AF by modulating atrial fibrosis and structural remodeling. Additionally, the ceRNA subnetwork found that downregulated miR-623 and its upregulated HCK and CXCR4 target genes were coregulated by circ_0006725 and circ_0079284. Ring Finger 216 (RNF216) is the host of has_circRNA_0079284, and study reported that it may be involved in innate immune signaling [55]. Study on intervertebral disc degeneration has shown that miR-623 directly bound CXCL12 to reduce levels of inflammatory factors in LPS-injured nucleus pulposus cells [56].
ose indicated that circRNA_0079284, miR-623, HCK, and CXCR4 may play roles in immune and inflammation. Our result suggested that the interaction between circRNA_0079284-miR-623-HCK/CXCR4 may be associated with the immune mechanism of AF. Other circRNA-miRNA-mRNA regulatory axes composed of these 7 hub genes in the ceRNA subnetwork may also be related to immunity (Figure 6(e)).
However, there are some limitations to this literature. First, only a few samples were used for miRNA and circRNA microarray analysis. Second, not all samples used for miRNA microarray analysis were from left atrial appendage, and tissue samples from different parts of atrial of patients with AF did not achieve homogeneity between samples, which may produce bias to the results. ird, in this study, the genes that predicted AF were derived from tissue samples, and the biomarkers for diagnosing persistent AF need to be tested in blood samples from those patients. Finally, further in vitro and in vivo experiments are needed in the future to validate the results deduced by bioinformatics analysis.

Conclusion
rough the analysis of immune-related ceRNA networks, our findings provide novel insight into the molecular mechanisms underlying the progression of AF. And we found that the circ_0005299-miR-1246-C5AR1 and circRNA_0079284-miR-623-HCK/CXCR4 regulatory axises may be associated with the immune mechanism of AF.
Data Availability e data sets analyzed in this study are available from the corresponding author upon reasonable request.