The Construction and Comprehensive Analysis of a ceRNA Immunoregulatory Network and Tissue-Infiltrating Immune Cells in Atrial Fibrillation.

BACKGROUND: At present, the mechanisms behind atrial fibrillation (AF) pathogenesis are still unclear. We construct a ceRNA immunoregulatory network to further understand the mechanism of AF.
METHODS: Four AF mRNA datasets from the Gene Expression Omnibus (GEO) database were integrated by SVA method. AF-related immune genes (AF-IRGs) were selected via combining ImmPort database with the genes in the module most associated with AF obtained by a weighted gene coexpression network analysis (WGCNA). Then, circRNA and miRNA expressions from the GEO database were extracted and mapped with related databases. Next, an immune-related circRNA-miRNA-mRNA ceRNA network was constructed and hub genes were filtered from a protein-protein interaction (PPI) network, and the differentially expressed (DE) hub genes in AF were further screened. Additionally, immune infiltration was investigated in AF by using CIBERSORT. Subsequently, the relationships between DE hub genes and AF-related infiltrating immune cells were performed by using Pearson correlation coefficients. Ulteriorly, the immune-cells-related ceRNA subnetwork in AF was built.
RESULTS: A total of 95 AF-IRGs were detected, and an immune-related ceRNA network in AF was constructed with 12 circRNAs, 7 miRNAs and 50 mRNAs. The immune infiltration analysis indicated that a higher level of neutrophils, as well as a lower level of T cells regulatory (Tregs) and NK cells activated in AF. Four DE hub genes (CXCL12, IL7R, TNFSF13B, CD8A) were associated with Tregs or NK cells activated immune cells (P < 0.05). Tregs or NK cells activated immune cells-related ceRNA subnetwork including 5 circRNAs (has_circ_0001190, has_circ_0006725, has_circ_0079284, has_circ_0005299, and has_circ_0002103), 4 miRNAs (has-miR-198, has-miR-623, has-miR-1246, and has-miR-339-3p) and 4 DE hub genes was eventually constructed in AF.
CONCLUSION: Our results provide new insights into the molecular mechanisms governing AF progression from the perspective of immune-related ceRNA network.

Introduction

Atrial fibrillation (AF) is the most frequently occurring type of arrhythmia and its incidence increases with age. In addition, AF is closely associated with poor quality of life, embolic stroke, heart failure, and elevated mortality.1 However, the underlying mechanisms of AF remain unclear. Understanding the mechanisms could help us to find new strategies for treating AF. Therefore, it is necessary to add our further insights into cellular and biological mechanisms of AF and to evolve more effective treatments for AF.

At present, the mechanisms behind AF pathogenesis are mainly reflected in atrial electrical and structural remodeling.2 Furthermore, research have shown that a link between inflammation and atrial fibrosis, and inflammation and its associated immune response may play a critical role in the occurrence and maintenance of AF.3,4 Study found that C-reactive protein (CRP) was independently associated with the occurrence of AF, and the incidence of AF increased with the increase of CRP expression.5 In addition, several immune-related genes (IRGs), such as C-X-C motif chemokine receptor 2 (CXCR2) and C-X-C motif chemokine 10 (CXCL10), may be implicated in the structural remodeling of AF.6,7 The infiltration of immune cells has also been shown to influence the development of AF. For example, monocyte activation especially the enhanced migration ability plays a vital role in the pathogenesis of progression in atrial remodeling in AF patients.8

The role of the competing endogenous RNA (ceRNA) network has attracted much attention in recent years.9 As a novel type molecule with special biological functions, circRNAs can be used as microRNA (miRNA) sponge to competitively bind miRNAs, modulate alternative splicing, and regulate the expression of parental genes.10,11 Emerging data have elucidated that circRNA-miRNA-mRNA regulatory axes are implicated in the pathological progression of cardiovascular disease. For example, ciRS-7 can be used as miRNA-7a sponge to inhibit the expression of PARP and SP1, resulting in promoting myocardial apoptosis.12 CircRNA_000203 can enhance the expression of fibrosis-associated genes (Col1a1 and CTGF) by suppressing targets of miR-26b-5p, leading to cardiac fibrosis, which is associated with AF.13 Therefore, circRNA-miRNA-mRNA regulatory network can deliver novel insights into the pathophysiology of AF. However, few studies have explored the relationship between immune infiltration and an immune-related circRNA-miRNA-mRNA ceRNA network in AF.

In the present study, we provided new insights into the molecular mechanism of ceRNA immunoregulation during the progression of AF. Firstly, we identified AF-related immune genes (AF-IRGs) via combining ImmPort database with the genes in the module most associated with AF obtained by a weighted gene coexpression network analysis (WGCNA). Then, we constructed an immune-related circRNA-miRNA-mRNA ceRNA network based on circRNA-miRNA pairs and miRNA-mRNA pairs, and filtered differentially expressed (DE) hub genes from a protein–protein interaction (PPI) network. In addition, we performed an immune cell infiltration analysis in AF by using CIBERSORT.14 Finally, we selected DE hub genes associated with AF-related infiltrating immune cells by Pearson correlation coefficients to construct an immune-cells-related ceRNA subnetwork. Our study flowchart is shown in Figure 1.

Figure 1

Flow diagram of study.

Methods

Data Acquisition and Processing

One circRNA microarray dataset (GSE12940915), one miRNA microarray dataset (GSE7088716) and four mRNA microarray datasets (GSE41177,17 GSE79768,18 GSE31821,19 and GSE1497520) were downloaded from gene expression omnibus (GEO, ).21 All samples are from atrial tissue and the samples are classified as persistent AF and SR control groups; GSE129409 and GSE70887 were obtained from GPL21825 Arraystar Human circRNA Array V2, and GPL19546 Agilent-021827 Human miRNA Microarray respectively. Four mRNA datasets were generated using GPL570 (HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array. The GSE129409 contained left atrial appendage tissue from 3 persistent AF patients and 3 sinus rhythm (SR) controls. The GSE70887 contained atrial appendage tissue from 4 persistent AF patients and 2 SR controls. The GSE41177, GSE79768, GSE31821 and GSE14975 contained left atrial appendage tissue from 16 persistent AF patients and 3 SR controls, 7 persistent AF patients and 6 SR controls, 2 persistent AF patients and 2 SR controls, and 5 persistent AF patients and 5 SR controls, respectively.

The probe IDs were converted to symbol codes by using ActivePerl 5.24.2 software () to process series matrix files. The GSE14975 was log base 2 transformed by affy in R.22 After merging mRNA microarray data, the combat function in the sva package was used to eliminate batch effects and other unwanted variation between four mRNA datasets.23 Finally, expression values were normalized using the “normalizeBetweenArrays” function of the package of “limma” in R software.24

The study protocol was approved by the Ethics Committee of Xiangtan central Hospital.

Identification of Differential Expression of circRNAs, miRNAs and mRNAs

The Limma package was utilized to screen DE circRNAs, DE miRNAs and DE mRNAs between SR samples and AF samples.25 The GSE129409 was analyzed with |log2 FC| > 3 and p.value <0.05 set as the cut-off point for the selection of DE circRNAs. To assess DE miRNAs in AF, GSE70887 was used, and transcripts with a cut-off point of |log2 FC| > 1 and p.value <0.05 were retrieved. The integrated dataset was analyzed with p.value <0.05 set as the cut-off point for the selection of DE mRNA. Heatmaps and volcano maps for DE circRNAs and DE miRNAs were created by using the “pheatmap” and “ggplot2” package of R software, respectively.

Weighted Gene Co-Expression Network Analysis (WGCNA) Construction and Identification of AF-IRGs

To identify candidate biomarkers and therapeutic targets, a system biology approach, WGCNA was used in gene co-expression network construction. In our study, the mRNA integrated dataset was constructed to gene co-expression network using the WGCNA package in R.26 The soft thresholding power was set as 4; then, analyses of module–trait correlations, gene significance (GS), and module membership (MM) were performed. Statistically significant modules were defined as P < 0.05.

AF-IRGs were screened by matching 1793 IRGs to genes in the most relevant module with AF using the IRG dataset from the ImmPort database ().

Construction of an Immune-Related circRNA-miRNA-mRNA ceRNA Network

Information about circRNAs can be obtained in circBase ().27 Target miRNAs were predicted by the cancer-specific circRNA database (CSCD, ).28 For each DE circRNA, all predicted miRNAs were obtained. Next, miRNAs that overlapped with predicted and DE miRNAs were gathered. Then, the TargetScan database () was used to predict the miRNA-targeted mRNAs.29 To this end, mRNAs recognized by TargetScan database were regarded as candidate targets and were intersected with AF-IRGs. Finally, an immune-related ceRNAs regulatory network in AF was constructed based on circRNA-miRNA pairs and miRNA-mRNA pairs. The ceRNA network was visualized by Cytoscape version 3.8.0 ().

Enrichment Analysis

GO and KEGG functional enrichment analysis of AF-IRGs, and genes in the immune-related ceRNA network, respectively, were performed using the “clusterprofiler”, “org.Hs.eg.db”, “enrichplot”, “ggplot2 Goplot” and “digest” packages of R/Bioconductor. p.value <0.05 was regarded as a statistically significant difference.

Construction of PPI Regulatory Network and Identification of DE Hub Genes

A PPI network of these genes in the immune-related ceRNA network was established by using the STRING database (),30 and minimum interaction score ≥0.4 was defined as the cut-off value. The PPI network was visualized by the CytoCope 3.8.0 software. In addition, top 10 hub genes in the PPI network were filtered by using the MCC algorithm of the Cytoscape plug-in cytoHubba, and DE hub genes were screened by intersection with the DE mRNAs.

Immune Cell Infiltration Analysis

To quantify the relative proportions of infiltrating immune cells in AF, CIBERSORT in R, was adopted to analyse the integrated mRNA dataset and calculate immune cell infiltrations. P < 0.05 was used to filter the samples. Principal component analysis (PCA) cluster plot of immune cell infiltration in AF and SR tissues was carried out by the “ggplot2”, “FactoMineR”, and “factoextra” package. Spearman correlation analysis using R software was applied to the 22 types of infiltrating immune cells, and a correlation heatmap was performed using the “Corrplot” package to visualize the result. The “vioplot” package was used to compare and visualize the levels of 22 immune cells between the AF and SR samples.

Correlation Analysis Between DE Hub Genes and AF-Related Infiltrating Immune Cells

The relationship between DE hub genes and AF-related infiltrating immune cells was performed by using Pearson correlation coefficients, which was visualized by using the “ggpubr” package of R.

Construction of the Immune Cells-Related ceRNA Subnetwork

After obtaining the DE hub genes related to AF-related immune cells by Pearson’s method (P < 0.05), we reconstructed an immune-cells-related ceRNA subnetwork from the immune-related circRNA-miRNA-mRNA ceRNA network and visualized it through the CytoCope 3.8.0 software. The boxplot maps of cicrRNAs, miRNAs and mRNAs expression levels in the subnetwork in their respective datasets were created by using the “ggpubr” and “reshape2” package of R software.

Diagnostic Analysis of mRNAs of ceRNA Subnetwork for AF

Receiver operator characteristic (ROC) curve analysis was performed using the “pROC” package of R to determine the effectiveness of these genes of ceRNA subnetwork in predicting AF.

Results

Identification of DE circRNAs and DE miRNAs in AF

In total, 103 DE circRNAs (59 up-regulated and 44 down-regulated circRNAs) were identified in the GSE129409 dataset (Figure 2A and B). In addition, 16 up-regulated and 13 down-regulated miRNAs were screened in the GSE70887 dataset (Figure 2C and D).

Figure 2

Identification of DE circRNAs, DE miRNAs in AF. Heat map (A) and volcano plots (B) of DE circRNAs in the GSE129409 dataset. Heat map (C) and volcano plots (D) of DE miRNAs in the GSE70887 dataset.

Coexpression Module Construction and AF-Related IRGs Screening

In order to find the functional clusters in AF patients, gene co-expression networks were constructed from the integrated mRNA dataset with the WGCNA package. After setting the power as 4, 9 modules were excavated (Figure 3A). The results of the module-trait relationships are presented in Figure 3B, revealing that the red module was the most relevant modules with AF (r = 0.45; P = 0.002), including 280 genes. As GS and MM exhibit a very significant correlation, this finding implies that the genes in the red module tend to be highly correlated with AF (cor = 0.5, P = 4.1e-19) (Figure 3C). In addition, the intersection of 1793 IRGs from the ImmPort database and the genes in red module was AF-related IRGs (AF-IRGs) (Figure 3D), including 95 genes.

Figure 3

Identification of modules associated with AF in the integrated mRNA dataset and AF-related IRGs screening. (A) A cluster dendrogram was built based on the dissimilarity of the topological overlap. (B) Module-trait relationships. The red module was significantly correlated with AF. (C) Correlation plot between MM (X-axis) and GS (Y-axis) of genes contained in the red module. (D) Venn diagram showing intersecting genes between IRGs and the genes in the red module.

Enrichment Analysis of AF-IRGs

The results of the GO analysis showed that the AF-IRGs were predominantly enriched in neutrophil chemotaxis (biological processes, BP), external side of plasma membrane (cellular components, CC), and immune receptor activity and chemokine activity (molecular functions, MFs) (Figure 4A). KEGG pathway analysis showed significant enrichment in the viral protein interaction with cytokine and cytokine receptor, cytokine and cytokine receptor interaction and chemokine signaling pathway (Figure 4B).

Figure 4

AF-IRGs enrichment analysis and Venn diagram. (A) GO enrichment analysis for AF-IRGs; (B) KEGG pathway enrichment analysis for AF-IRGs; (C) Venn diagram of overlapping miRNAs from predicted miRNAs and DE miRNAs in the GSE70887 dataset. (D) Venn diagram of overlapping mRNAs from AF-IRGs and predicted mRNAs.

Construction of an Immune-Related circRNA-miRNA-mRNA ceRNA Network in AF

Among the 103 DE circRNAs, 15 circRNAs were not searched in the CSCD database. We predicted that the remaining 88 circRNAs might target 1995 miRNAs by using the CSCD database. Then, 9 intersecting miRNAs were acquired by the intersection of target miRNAs and DE miRNAs (Figure 4C). Next, a total of 13,135 potential target genes were predicted for the 9 miRNAs by using the TargetScan database. Subsequently, 56 intersecting mRNAs were obtained by the intersection of target genes and AF-IRGs (Figure 4D). Finally, an immune-related circRNA-miRNA-mRNA ceRNA network in AF was constructed with 12 up-regulated circRNAs, 7 down-regulated miRNAs and 50 up-regulated mRNAs (Figure 5A and Table 1).

Figure 5

Construction and analysis of an immune-related ceRNA network. (A) The immune-related circRNA-miRNA-mRNA network in AF. Red diamonds indicate upregulated circRNAs, green triangles represent downregulated miRNAs, and blue ovals represent AF-IRGs; (B) GO enrichment analysis for mRNAs in the ceRNA network; (C) KEGG pathway enrichment analysis for mRNAs in the ceRNA network; (D) PPI network construction for mRNAs in the ceRNA network; (E) Hub gene selection from PPI network by using the MCC algorithm.

Table 1

Immune-Related ceRNA Network of circRNAs, miRNAs, and mRNAs in AF

circRNA	miRNA	mRNA
hsa_circ_0005299	hsa-miR-1246	MSR1, CYSLTR1, TNFSF10, HGF, C5AR1, IL7, CD8A, TNFSF13B, C3AR1, CCR7, CXCL13, PRKCQ, TLR8, AZGP1, CSF2RB
hsa_circ_0004771hsa_circ_0004390hsa_circ_0001615	hsa-miR-483-5p	C5AR1, CCR5, PLXNC1, RSAD2, CCR7, CD74, CD28, TLR7, NOD2, PNOC, CD3D
hsa_circ_0001190	hsa-miR-198	IL10RA, CCR7, ISG20, CXCL9, CTSS, HLA-DPB1, PIK3CG, CD74, CCR5, IL7R, HGF, RAC2
hsa_circ_0079284hsa_circ_0006725	hsa-miR-623	CMA1, CCR7, RSAD2, CXCL9, HLA-DPB1, CSF2RB, C3AR1, ITK, CD1E, SDC1, CCR5, RASGRP1, CD28, CX3CR1, IL7R, HGF, HCK, AZGP1, CD8A
hsa_circ_0002103	hsa-miR-339-3p	CD72, CSF2RB, CXCL16, CXCL12, CCR7, TNFSF13B, CCL26, SDC1
hsa_circ_0004979	hsa-miR-508-5p	RSAD2, CXCL16, PLXNC1, CX3CR1, CXCL10, CXCL13, LYZ, SDC1, CD3D, PIK3CG, CD28, TLR7, CTSS, HLA-DPB1, LCK, IL10RA, TLR8, PTGER2, CCR7, CXCL11, ITGAL, IL2RB
hsa_circ_0006324hsa_circ_0083092hsa_circ_0003638	hsa-miR-516a-5p	RSAD2

Enrichment Analyses for mRNAs in the Immune-Related ceRNA Network

The results of the GO analysis showed that the identified mRNAs in the network were predominantly enriched in T cell activation and cell chemotaxis (biological processes, BP), external side of plasma membrane (cellular components, CC), and receptor ligand activity (molecular functions, MFs) (Figure 5B). KEGG pathway analysis showed significant enrichment in the cytokine and cytokine receptor interaction and chemokine signaling pathway (Figure 5C).

PPI Network Analysis

A PPI network including 47 nodes and 594 edges were built for the 50 mRNAs of the ceRNA network by using the STRING database after removing unconnected nodes (Figure 5D). Hub genes in the PPI network were identified by using the MCC algorithm. According to the MCC sores, the top ten highest-scored genes (Figure 5E), including C-X-C motif chemokine 10 (CXCL10), C-X-C motif chemokine 9 (CXCL9), T-cell surface glycoprotein CD8 alpha chain (CD8A), Interleukin-7 (IL-7), Interleukin-7 receptor subunit alpha (IL7R), Tumor necrosis factor ligand superfamily member 13B (TNFSF13B), Stromal cell-derived factor 1 (CXCL12), CD28, C-C chemokine receptor type 5 (CCR5) and C-C chemokine receptor type 7 (CCR7) were selected as the hub genes, and 5 of them were upregulated in the AF patients compared to the SR controls, known as DE hub genes (CD8A, IL7R, TNFSF13B, CXCL12, CCR5) (Table 2).

Table 2

DE Hub Genes with Higher MCC Score in Protein–Protein Interaction Network

Gene Symbol	MCC Score	logFC	P value	Gene Title
CD8A	473738	0.431	0.036	T-cell surface glycoprotein CD8 alpha chain
IL7R	342360	0.440	0.031	Interleukin-7 receptor subunit alpha
TNFSF13B	339480	0.487	0.018	Tumor necrosis factor ligand superfamily member 13B
CXCL12	285,522	0.500	0.012	Stromal cell-derived factor 1
CCR5	238,786	0.504	0.005	C-C chemokine receptor type 5

Note: Positive logFC values correspond to up-regulated genes.

Abbreviations: DE, differentially expressed; FC, fold change; MCC, maximal clique centrality.

Immune Infiltration Analyses

The proportion of 22 immune cells detected by the CIBERSORT algorithm in each sample were displayed in the bar plot and heat map (Figure 6A and B). PCA was used to analyze the immune cells in AF and SR tissues, and it was found that there were significant group bias clustering and individual differences (Figure 6C). The correlation of 22 types of immune cells are shown in Figure 6D. Compared with SR, AF patients commonly contained a higher proportion of neutrophils, and a lower proportion of T cells regulatory (Tregs) and NK cells activated (P < 0.05) (Figure 6E).

Figure 6

Immune infiltration analyses for AF. The proportion (A) and heatmap (B) of 22 immune cells detected by CIBERSORT algorithm. (C) The principal component analysis for immune cell in AF and SR tissues. (D) Correlation heat map for immune cells. (E) The violin plot of immune cells.

A total of 4 DE hub genes were linearly associated with Tregs or NK cells activated. Tregs were negatively associated with CXCL12 (R = −0.4, P = 0.007), IL7R (R = −0.36, P = 0.013) and TNFSF13B (R = −0.38, P = 0.008). NK cells activated were negatively correlated with CD8A (R = −0.34, P = 0.022), IL7R (R = −0.36, P = 0.013) and TNFSF13B (R = −0.41, P = 0.004) (Figure 7).

Figure 7

The correlation between DE hub genes and AF-related immune cells. The relationships between DE hub genes and AF-related immune cells were calculated using Pearson correlation coefficients (A–F).

Construction of the Immune Cells-Related ceRNA Subnetwork and ROC Curve Analysis

Four DE hub genes (CD8A, IL7R, TNFSF13B, and CXCL12) associated with AF-related immune cells (Tregs and NK cells activated), 4 miRNAs and 6 circRNAs were eventually used to construct immune cells-related ceRNA subnetwork (Figure 8). The expression levels of cicrRNAs, miRNAs and mRNAs in the ceRNA subnetwork in their respective datasets are shown in Figure 9. The expression of two genes (TNFSF13B and CXCL12) were significantly associated with a diagnosis of AF (0.7

Figure 8

Construction of the immune cells-related circRNA-miRNA-mRNA subnetwork. Red diamonds indicate upregulated circRNAs, green triangles represent downregulated miRNAs, and yellow ovals represent upregulated mRNAs.

Figure 9

The boxplot maps of cicrRNAs (A) miRNAs (B) and mRNAs (C) expression levels in the immune cells-related ceRNA subnetwork in their respective datasets. *P < 0.05, **P < 0.01, ***P < 0.001.

Figure 10

The receiver operator characteristic curves of CXCL12 and TNFSF13B for AF.

Discussion

In recent years, more and more evidences have revealed that the pathogenesis of AF is related to immune and inflammatory responses. And emerging data have elucidated that circRNA-miRNA-mRNA ceRNA network is implicated in the pathological progression of cardiovascular disease. Therefore, we sought to perform a combined analysis of a ceRNA immunoregulatory network and tissue-infiltrating immune cells in AF to provide theoretical basis for finding new targets for the prevention and treatment of AF.

Firstly, we screened 95 AF-IRGs, and the enrichment function confirmed that these genes were primarily associated with immune and inflammatory responses. Together with AF-IRGs, we constructed an immune-related circRNA-miRNA-mRNA ceRNA network in AF, indicating that these circRNAs, miRNAs and mRNAs in the network may play significant roles in AF immune mechanisms. The results of the GO enrichment analysis showed that 50 mRNAs in the ceRNA immunoregulatory network were also mainly implicated in the inflammation and immune system. Further, KEGG pathway analysis results were predominantly enriched in cytokine-cytokine receptor interactions and the chemokine signaling pathway, which showed a close relationship with myocardial ischemia/reperfusion injury as reported in Liu et al’s study.31 Some studies have shown that chemokines and their receptors regulate immune cell recruitment and activation and participate in persisting of the local inflammatory response in atherosclerosis,32 as well as the pathophysiology of left ventricular remodeling and cardiac dysfunction caused by pressure overload.33 These findings indicated that cytokine-cytokine receptor interactions and the chemokine signaling pathway play an important role in the cardiovascular field, which may be associated with the development of AF. Subsequently, a PPI network was established by using mRNAs in the ceRNA immunoregulatory network and we filtered 5 DE hub genes (CD8A, IL7R, TNFSF13B, CXCL12, CCR5). Among them, 4 hub genes (CCR5, CXCL12, TNFSF13B, IL7R) were enriched in cytokine-cytokine receptor interactions and 2 hub genes (CCR5, CXCL12) were enriched in chemokine signaling pathway. Therefore, our results suggested that these DE hub genes may be primarily involved in the immune and inflammatory processes of AF.

Next, immune infiltration was investigated in AF by CIBERSORT. In this program, we found that a higher proportion of neutrophils, as well as a lower proportion of T cells regulatory (Tregs) and NK cells activated may be associated with the pathogenesis of AF. Neutrophils are the most abundant circulating leukocyte in humans, and study found that myocardial inflammatory infiltration may be a cause of AF, including neutrophil and inflammation markers.34 Moreover, Wu et al35 showed that neutrophil-to-lymphocyte ratio could be an independent predictor of long-term outcomes in patients with AF. Chen et al36 showed that the proportion of Tregs was obviously lower in AF patients, which was consistent with the results of our study. Another study confirmed that Tregs mitigate cardiac hypertrophy and myocardial fibrosis in Angiotensin II–induced hypertensive mice.37 NK cells account for more than 1% of cardiac lymphocytes and are involved in the regulation of cardiac diseases.38 Reduction in the number and function of NK cells has been found in patients with coronary artery disease (CAD),39 especially in patients with acute coronary syndrome.40 In addition, IFN-γ-producing NK cells have been demonstrated to negatively regulate liver fibrosis, suggesting that low expression of NK cells may be involved in atrial fibrosis,41 which requires further experimental verification. Therefore, neutrophils, Tregs and NK cells may be potential key subtypes implicated in driving the progression of AF disease.

Ulteriorly, our results found that only 4 DE hub genes (CD8A, IL7R, TNFSF13B, and CXCL12) were associated with Tregs or NK cells activated, but not with neutrophils by Pearson correlation coefficients, and we constructed the immune cells-related ceRNA subnetwork. And only CXCL12 and TNFSF13B could accurately predict AF (0.7

Human B cell activating factor (TNFSF13B, BAFF) is a tumor necrosis factor superfamily member. Study has shown that BAFF can inhibit Tregs, thereby promoting allergic inflammation of the airway,52 and another study has also shown that BAFF is involved in the process of myocardial fibrosis after acute myocardial infarction (MI).53 Our results found that Tregs expression was low in the left atrial tissue of AF patients and was negatively correlated with TNFSF13B. These reports suggested TNFSF13B may be involved in the immune inflammatory mechanism of AF. Additionally, the sub-network showed that downregulated miR-1246 and miR-339-3p, and target upregulated TNFSF13B were under the regulation of has_circ_0005299 and has_circ_0002103. Huang et al54 revealed that miR-1246 in exosomes derived from Endothelial progenitor cells can enhance angiogenesis in vitro and in vivo after MI, and these improvements may be related to amelioration of myocardial injury and myocardial fibrosis after MI. So, we speculated that the low expression of miR-1246 in the heart tissue of AF patients may be related to promoting cardiac fibrosis. Therefore, we speculated that has_circ_0005299-miR-1236-TNFSF13B axis and has_circ_0002103-miR-339-3p-TNFSF13B axis associated with Tregs and NK cells might be related to the immune mechanism of AF.

The remaining two genes (CD8A, IL7R) in the immune-cells-related ceRNA subnetwork have not been reported in AF studies, but animal study55 has suggested that CD8A+ immune cell infiltration may be related to cardiac fibrosis and left ventricular contractility decline in diabetic cardiomyopathy and another study56 has shown that IL7R is involved in the inflammatory response in atherosclerosis. As for the immune cells-related ceRNA subnetwork constructed by us, it was found that the upregulated circRNAs (has_circ_0001190, has_circ_0006725, has_circ_0079284, and has_circ_0005299), downregulated miRNAs (has-miR-198, has-miR-623, and has-miR-1246) and upregulated mRNAs (CD8A and IL7R) may be immunologically related to AF, but the specific mechanism still needs further study.

However, our study has some limitations. First, AF sample size for circRNA and miRNA microarray analysis was small in the present study. Second, the samples analyzed for miRNA were not from the left atrial appendage, and the difference between the left and right atrial in AF patients may bias the results. Thirdly, for better clinical significance, the biomarkers for the diagnosis of AF should be verified in blood samples of AF patients. Finally, the results of this study need to be further verified by future in vivo or in vitro experiments.

Conclusion

Our study found that an immune-related circRNA-miRNA-mRNA ceRNA network in AF which includes 5 circRNAs (has_circ_0001190, has_circ_0006725, has_circ_0079284, has_circ_0005299, and has_circ_0002103), 4 miRNAs (has-miR-198, has-miR-623, has-miR-1246, and has-miR-339-3p) and 4 mRNAs (CXCL12, IL7R, TNFSF13B, and CD8A) may depend on the interaction among tissue-infiltrating immune cells, including Tregs or NK cells activated. Our results may provide new insights into the molecular mechanisms governing AF progression, and provide theoretical basis for finding new targets for the prevention and treatment of AF.