Uncovering potential novel biomarkers and immune infiltration characteristics in persistent atrial fibrillation using integrated bioinformatics analysis

1.   Introduction

Atrial fibrillation (AF) is the commonest cardiac arrhythmia in the general population, associated with high mortality and morbidity. Based on AF duration, AF can be divided into paroxysmal AF, persistent AF (pAF), long-standing pAF, and permanent AF. AF is estimated to affect 34 million people worldwide and is increasing as the population ages ^[1]. AF can induce major adverse cardiovascular events, including higher risk of stroke, heart failure ^[2]. However, the mechanisms of AF are complex and incompletely understood. Inflammation and atrial fibrosis were reported to be the fundamental mechanisms of the pathophysiological processes in the development of AF ^[3]. Soluble biomarkers of inflammation, including CRP and IL-6, and fibrosis, containing TGF-1 and matrix metalloproteinases, were used for predicting AF ^[4].

MicroRNA (miRNA), a small non-coding RNA molecule, plays a crucial role in RNA silencing and post-transcriptional gene expression regulation. Recent studies have elucidated the involvement of miRNAs in the pathological process of cardiovascular disease ^[5], including regulating arrhythmogenesis and cardiac excitability ^[6]. Moreover, miRNAs can be involved in AF-induced ion channel remodelling and fibrosis ^[7]. Therefore, miRNAs can provide useful insights into the pathophysiology of AF. Additionally, regulation of immune cells is critical in the pathogenesis and recurrence of AF ^[8]. Previous studies proved that increased infiltration of immune cells could be related to the enlargement of the left atrium of AF patients, a risk factor for AF recurrence ^[9]. Although some researchers also focused on the function of various immune cells, such as macrophages, monocytes, neutrophils ^[10], how these cells change in specific subtype of AF still remains blank.

At present, microarray technology, together with integrated bioinformatics analysis, have been extensively applied to identify potential novel biomarkers and their roles in various diseases to further explore the pathogenesis and develop potential treatments ^[11,12]. Additionally, Weighted gene co‑expression network analysis (WGCNA) is gradually becoming a mature strategy in bioinformatics applications and can be used to identify underlying mechanisms, potential biomarkers or therapeutic targets in different types of diseases ^[13].

To our knowledge, few studies have investigated miRNAs, their target genes and immune cells in the heart tissue of patients with pAF. In the present study, three datasets (GSE31821, GSE41177 and GSE79768) were merged into an integrated dataset by the SVA method to remove batch differences. To identify the pivotal genes, miRNAs, pathways associated with pAF, differential expression analysis, WGCNA, enrichment analysis, as well as protein‐protein interaction (PPI) network, miRNA-TF-mRNA network were constructed successively. Finally, immune infiltration was investigated in pAF by CIBERSORT. The study focused on novel findings in pAF mechanisms and may help to the identification of potential biomarkers for diagnosis and treatment of pAF.

2.   Materials and methods

2.1.   Microarray data

The study flowchart is shown in Figure 1. Three pAF datasets (GSE31821, GSE41177, GSE79768) were downloaded from the NCBI Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) ^[14]. The above three datasets were all gene expression array and generated using GPL570 (HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array (Affymetrix, Santa Clara, CA) ^[15]. The dataset of GSE31821 included the auricle tissues of 4 patients with pAF and 2 normal samples; GSE41177 contained 38 atrial tissue samples, including 32 patients with pAF and 6 normal atrial tissue samples; GSE79768 consisted of 26 heart tissue samples including 14 pAF atrium heart tissue samples and 12 normal heart tissue samples.

Figure 1.  Study flowchart. Abbreviations: pAF, persistent atrial fibrillation; LIMMA, Linear models for microarray data.

DownLoad: Full-Size Img PowerPoint

2.2.   Data processing and identification of differentially expressed genes

The three raw datasets were pre-processed by affy in R, including background calibration, normalization, and log2 transformation ^[16]. When multiple probes corresponded to one common gene, the average value was taken as its expression value. Besides, the surrogate variable analysis (SVA) package from Bioconductor was used to remove batch effects and other unwanted variation between three datasets ^[17]. Finally, the LIMMA package was then applied to screen DEGs and a p-value < 0.05, and |log₂ Fold change (FC)| > 0.5 was set as a cut-off point for selecting DEGs ^[18].

2.3.   Weighted gene co-expression network analysis (WGCNA) construction and identification of modules

To explore the interactions between genes, a system biology approach, WGCNA, was applied for gene co-expression network construction. First, genes having upper 25% variation across samples in the integrated dataset were imported to the WGCNA. Second, to ensure that the results of network construction were reliable, outlier samples were removed. Third, Adjacency was calculated from the soft thresholding power β, which was derived by co-expression similarity, using pick-Soft-Threshold function. Then, the adjacency was transformed into a topological overlap matrix (TOM), and the corresponding dissimilarity (1-TOM) was calculated. Fourth, modules were detected through hierarchical clustering and dynamic tree cut function. To classify genes with similar expression profiles into gene modules, average linkage hierarchical clustering was conducted according to the TOM‐based dissimilarity measure with a minimum size (gene group) of 50 for the genes dendrogram ^[19]. Fifth, for the modules correlated to the clinical attributes, module membership (MM) and gene significance (GS) were calculated. Finally, the eigengene network was visualized. Information on genes in the modules were used for further analysis. The intersection of the DEGs screeded from the integrated dataset and the genes in significant module is common genes (CGs).

2.4.   Gene ontology (GO) annotation and kyoto encyclopaedia of genes and genomes (KEGG) pathway enrichment analyses of the CGs

DAVID (Database for Annotation, Visualization and Integrated Discovery) 6.8 (https://david.ncifcrf.gov/) is a gene functional classification tool that was applied for enrichment analysis to perform the Gene Ontology (GO) for the CGs ^[20]. The enriched GO terms were divided into Biological Process (BP), Cellular Component (CC), and Molecular Function (MF) ontologies. To understand the biofunction of the CGs, we uploaded the CGs to the DAVID database (Gene Functional Classification Tool, http://david.abcc.ncifcrf.gov/) ^[21] and KOBAS database (http://kobas.cbi.pku.edu.cn/kobas3) ^[22]. The significant enrichment for GO and KEGG analyses threshold was p‐value < 0.05 and count ≥ 2.

2.5.   Protein-protein interaction (PPI) network construction

The PPI interaction network of CGs was constructed with the Search Tool for the Retrieval of Interacting Genes (String) database (version 11.0; www.string-db.org) ^[23] and minimum interaction score ≥ 0.4 was defined as the cut-off value. Finally, the PPI network was visualized by Cytoscape ^[24].

2.6.   MiRNA (microRNA)-TF (transcription factor)-mRNA regulatory network construction

Target miRNAs of the CGs were predicted with the miRTarBase ^[25], Starbase ^[26] and Targetscan ^[27] database. To improve the accuracy of the prediction, only miRNAs predicted by all three databases were retained. Combined with the Enrichr database (http://amp.pharm.mssm.edu/Enrichr/), the (transcription factors) (TFs) targeted the CGs were predicted. The results with the p-value ≤ 0.05 were chosen as the cut-off values. After miRNA‐TF‐mRNA regulatory relationships were obtained, the miRNA‐TF‐mRNA regulatory network was visualised using Cytoscape.

2.7.   Immune cell infiltration analysis

To quantify the relative proportions of infiltrating immune cells in pAF, CIBERSORT, a method of analysing the different immune cell types of tissues, was adopted to analyse the merged expression data and calculate immune cell infiltrations ^[28]. p < 0.05 was used to filter the samples. And the percentage of each immune cell type in the samples was calculated and displayed in a bar plot. The heat map of the 22 immune cells was made using the "pheatmap" package ^[29]. The "vioplot" package was used to compare and visualize the levels of 22 immune cells between the pAF and control samples ^[30]. A correlation heatmap, which revealed the correlation of 22 types of infiltrating immune cells, was performed using the "corrplot" package.

2.8.   Statistical analysis

All data were expressed as the mean ± standard deviation (SD). Statistical analyses were performed using GraphPad Prism (version 8.3.4). A p‐value < 0.05 was considered significant for screening DEGs, module genes, enrichment analysis, PPI network, miRNA-TF-mRNA network and immune infiltration analysis.

3.   Results

3.1.   Identification of DEGs

A total of 264 genes were differentially expressed, of which 179 genes were up-regulated and 85 were down-regulated. The top 5 up-regulated genes involved S100 Calcium Binding Protein A12 (S100A12), C-X-C Motif Chemokine Receptor 2 (CXCR2), Myeloid Cell Nuclear Differentiation Antigen (MNDA), Nebulin Related Anchoring Protein (NRAP), and Four And A Half LIM Domains 2 (FHL2) while the top 5 down-regulated genes were Proteoglycan 4 (PRG4), Delta/Notch Like EGF Repeat Containing (DNER), Transcription Elongation Factor A Like 2 (TCEAL2), Family With Sequence Similarity 110 Member C (FAM110C), Mesothelin (MSLN) (Table S1). The heat map of DEGs and volcano plot of DEGs are shown in Figure 2A, B, respectively.

Figure 2.  Heat map and Volcano plot for the DEGs identified from the integrated dataset. (A) Each row of the heat map represents one DEG, and each column represents one sample, either normal or pAF. The red and blue colors represent upregulated and downregulated DEGs, respectively. (B) Red plot points represent upregulated DEGs, and blue plot points show downregulated DEGs. Abbreviations: DEG, differentially expressed gene; pAF, persistent atrial fibrillation.

DownLoad: Full-Size Img PowerPoint

3.2.   Identification of gene co-expression networks and modules

Firstly, the genes were ranked by variance from large to small, and the top 25% variance (5044) genes were selected for analysis. Secondly, the flashClust tools package was applied to perform the cluster analysis, by setting the threshold value to 80 and 1 outlier samples were subsequently detected and removed. Cluster 1 contains the 69 samples we want to keep. Thirdly, the power parameter ranging from 1–20 was screened out using the 'pickSoftThreshold' function in WGCNA package. In our research, we selected the power of β = 8 (scale free R² = 0.85) as the soft threshold to ensure a scale-free network. Set the threshold to 0.25 to merge similar modules in the cluster tree. As showed in Figure 3A, a total of 12 modules were obtained, in which genes had similar co-expression traits. The colours of all the modules were selected at random to distinguish between modules. The module eigengene (ME) in the magenta module (r = 0.48; p = 3E−5) exhibited the highest positive correlation and had significant correlation with pAF compared with other modules (Figure 3B). Therefore, the magenta module, including 616 genes, was identified as key modules for pAF and retained for further analysis. Additionally, we calculated correlation between gene's module memberships (MMs, correlation between a specific gene and module's eigengene) and gene significances (GSs, correlation between a specific gene and clinical variable) in the magenta module. As expected, significant positive correlations between gene's MM and GS were obtained in the magenta module (cor = 0.44, P = 1.5E-30) as shown in Figure 3C. The intersection of the DEGs screened from the integrated dataset and the genes in magenta module is common genes (CGs) (Figure 3D). The genes in the magenta module were clearly most significantly associated with pAF.

Figure 3.  Weighted gene co-expression network analysis and Venn diagram. (A) GCMs, represented by different colors under the gene tree. (B) Heat map of the association between modules and pAF. The magenta module was significantly correlated with pAF. The numbers inside and outside of the brackets represent p-values and correlation coefficients, respectively. (C) Correlation plot between MM (X-axis) and GS (Y-axis) of genes contained in the magenta module. (D) Venn diagram showing overlapping genes between the DEGs and the genes in the magenta module. Abbreviations: GCMs, gene co-expression modules; GS, gene significance; MM, module membership. Color images are available online; DEGs, differentially expressed genes; pAF, persistent atrial fibrillation.

DownLoad: Full-Size Img PowerPoint

3.3.   Enrichment analyses of the common genes

Functional enrichment analysis indicated that the CGs were mainly involved in biological process (BP) terms, including "angiogenesis", "osteoblast differentiation", "cellular protein modification process", "calcineurin-NFAT signalling cascade" and "endothelial cell morphogenesis". In the cell component (CC) ontology, the CGs were significantly enriched in "extracellular exosome", "Golgi membrane" and "Golgi apparatus". Molecular function (MF) analysis showed that the CGs were mainly enriched in "TPR domain binding", "protein binding", "integrin binding" (Table S2). Additionally, the KEGG pathway of the CGs was found to be enriched in "metabolic pathways", "protein processing in endoplasmic reticulum", "Oxytocin signalling pathway", "Hippo signalling pathway", "NOD-like receptor signalling pathway", "Rap1 signalling pathway" and "MAPK signalling pathway". (Figure 4) (Table S3). The results of BP revealed that two upregulated DEGs (RCAN1 and PPP3R1) were enriched in calcineurin-NFAT signalling cascade (p value = 0.023112).

Figure 4.  Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of common genes. The green bars represent biological processes, the red bars represent cellular components, the brown lines represent molecular functions, and the purple lines represent KEGG pathways. Abbreviations: BP, biological process; CC, cellular component; MF, molecular function; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

DownLoad: Full-Size Img PowerPoint

3.4.   PPI network

After removing the isolated nodes, the PPI network was constructed with 26 nodes and 27 edges, including 25 up-regulated genes and 1 down-regulated genes. The CGs in the PPI network were defined as hub genes, including interleukin enhancer binding factor 2 (ILF2), eukaryotic translation initiation factor 5A (EIF5A), ribosomal protein S11 (RPS11), snail family transcriptional repressor 2 (SNAI2), protein phosphatase 3 regulatory subunit B, alpha (PPP3R1), regulator of calcineurin 1 (RCAN1) and so on (Figure 5). In the network, the ILF2 and EIF5A had the highest degrees of connectivity; both connecting 5 nodes.

Figure 5.  Protein-protein interaction network. Red nodes represent upregulated genes, and the green node represents a downregulated gene.

DownLoad: Full-Size Img PowerPoint

3.5.   MiRNA‐TF‐mRNA regulatory network analysis

After miRNA‐gene and TF‐gene pairs were predicted, 97 miRNA‐TF‐mRNA regulatory relationships were obtained. The regulatory network (including 65 miRNAs, 14 transcription factors, 18 co‐upregulated genes, and 39 co-downregulated genes) was established by Cytoscape (Figure 6). Within the network, we discovered that PPP3R1 was targeted by Nuclear Factor of Activated T Cells 1 (NFATC1), whilst RCAN1 was targeted by Signal Transducer And Activator Of Transcription 2 (STAT2). We also found that miR-34a-5p could regulate RCAN1 and PPP3R1 at the same time.

Figure 6.  The miRNA‐TF‐mRNA regulatory network. The red balls represent upregulated genes, the blue triangles are predicted miRNAs, and the brown rhombuses represent TFs. Abbreviations: mRNAs, messenger RNAs; miRNAs, microRNAs; TFs, transcription factors.

DownLoad: Full-Size Img PowerPoint

3.6.   Immune cell infiltration analysis

We further performed the CIBERSORT algorithm to predict the immune cells infiltration between patients with pAF and the control group. The percentage of each of the 22 types of immune cells in each sample was shown in the bar plot and heat map (Figure 7A, B). The vioplot of the immune cell infiltration difference demonstrated that pAF patients had a higher level of monocytes, dendritic cells and neutrophils, and a lower level of CD8+ T cells and T cells regulatory (Tregs) compared with the control group (Figure 8A). The correlation of 22 types of immune cells revealed that activated Mast cells were positively related with T cells CD4 naive (r = 0.59) and Plasma cells were positively related with B cells naive (r = 0.41), whereas activated Mast cells were negatively related to resting Mast cells (r = −0.61), T cells CD8 were negatively related to resting T cells CD4 memory and Macrophages M2 negatively correlated with T cells gamma delta (r = −0.45) (Figure 8B). There is a significant difference in immune cell infiltration between heart tissue in patients with pAF and the normal rhythm group. Therefore, monocytes, dendritic cells, neutrophils, CD8+ T cells and T cells regulatory (Tregs) may be potential core immune cells, involved in stimulating the progression of pAF.

Figure 7.  The landscape of immune infiltration between pAF and normal controls. (A) The relative percentage of 22 types of immune cells. (B) The heat map of 22 types of immune cells. Abbreviations: pAF, persistent atrial fibrillation.

DownLoad: Full-Size Img PowerPoint

Figure 8.  Distribution and visualization of immune cell infiltration. (A) Comparison of 22 immune cell subtypes between patients in pAF and controls. Blue and red colors represent normal and pAF samples, respectively. (B) Correlation matrix of all 22 immune cell subtype compositions. Both horizontal and vertical axes demonstrate immune cell subtypes. Immune cell subtype compositions (higher, lower, and same correlation levels are displayed in red, blue, and white, respectively).

DownLoad: Full-Size Img PowerPoint

4.   Discussion

AF is the most common cardiac arrhythmia worldwide and the mechanisms of AF has not been fully excavated. Previous studies of gene expression files were limited by single method or small samples. In our study, we integrated three persistent AF datasets, used two different methods to identify common genes (WGCNA and LIMMA method). With the help of the LIMMA method, we identified a total of 264 DEGs between pAF samples and normal controls, amongst which 179 genes were up-regulated and 85 were down-regulated. Genes such as S100A12 and CXCR2 were significantly upregulated whereas PRG4 and DNER were significantly downregulated. The mechanisms for the occurrence and maintenance of AF are associated with atrial fibrosis, inflammation and oxidative stress ^[31]. S100A12 belongs to the S100/calgranulin subfamily of Ca2+-binding proteins, which induces an inflammatory response leading to the accumulation of immune cells in damaged tissue. The chemokine receptor CXCR2 is a critical regulator in hypertension and cardiac remodeling. Zhang, Y.L. identified that blocking CXCR2 prevents and reverses the development of AF in spontaneously hypertensive rats. Consequently, CXCR2 may be a potential therapeutic target for hypertensive AF. ^[32]. Proteoglycan is a structural constituent of the aortic valve and accumulates in areas of fibrosis and moderate calcification. PRG4 has been found in fibrotic parts of the aortic valves. Herder, C. et al. identified that DNER might operate as a biomarker of inflammation associated with the pathogenesis of the distal sensorimotor polyneuropathy ^[33]. In the WGCNA, we identified 12 gene modules, amongst which the magenta module, with 616 genes, has the highest correlation with pAF.

To make the result more accurate and reliable, we made the intersection of the genes in the most significant module and the differentially expressed genes in LIMMA method. 61 common genes were identified, enrichment and pathway analyses were further performed. Moreover, PPI network and ceRNA network were constructed. The biological process analysis revealed pAF was mainly related to calcineurin-NFAT signaling cascade.

Calcineurin-NFAT (nuclear factor of activated T-cells) signaling is critical for numerous biological processes, He RL and his colleagues ^[34] proved Calcineurin-NFAT could regulate pulmonary arterial hypertension rats through modulating the proliferation, migration and apoptosis of pulmonary artery smooth muscle cell. Zhou X et al ^[35] demonstrated that the activation of calcineurin-NFAT signaling pathway could lead to cardiac hypertrophy. In a rabbit AF model, the expression level of calcineurin and NFAT was up-regulated ^[36], speculating its role in regulating AF. However, the relationship of this signaling pathway and AF has not been fully reported. With the consideration of PPI network and gene enrichment, RCAN1 and PPP3R1 were enriched in the signaling pathway. RCAN1 (regulator of calcineurin 1) belongs to RCAN family, of which, RCAN2 and RCAN3 are homologous ^[37]. RCAN was reported to be associated with various diseases including Down Syndrome ^[38], Alzheimer's disease ^[39] and diabetes ^[40]. RCAN could bind to calcineurin, a Ca2+/calmodulin-dependent phosphatase, inhibits calcineurin activity. PPP3R1 is an important subunit of calcineurin ^[41]. Therefore, we conjectured that RCAN1 and PPP3R1 could regulate pAF through calcineurin-NFAT signaling pathway.

Atrial electrical and structural remodeling is the fundamental pathogenesis of AF ^[42]. MicroRNAs (miRNAs) are essential post-transcriptional modulators of gene expression. Lv, X., et al. found overexpression of miR-29b-3p reduced atrial structural and electrical remodeling in a rat fibrosis model, through the targeted regulation of the PDGF-B signaling pathway ^[43]. Taken miRNA-TF-mRNA network into account, miR-34a-5p functions as the bridge to connect RCAN1 and PPP3R1. Broad evidence proved the vital role of miR34a in cardiac disorders. Researchers found silencing miR-34 could attenuate cardiac dysfunction in moderate hypertrophic cardiomyopathy ^[44]. Zhu Y et al ^[45] demonstrated that miR34a may serve an important role in early electrophysiological remodeling and the development of AF via the regulation of AnkB expression. Furthermore, Diener C et al proved that PPP3R1 could target miR-34a-5p in regulating intracellular calcium signaling ^[46], one of the basic mechanisms of AF. While RCAN1 was identified as a novel direct target of miR-34a in regulating endothelial inflammation after cardiac bypass in goat placenta ^[47]. We speculated that miR-34a-5p might target RCAN1 and PPP3R1 in regulating AF through Calcineurin-NFAT signaling pathway. Therefore, miR-34a-5p/RCAN1 and miR-34a-5p/PPP3R1 may be viewed as potential therapeutic targets for pAF.

Additionally, the critical role of immune response in the pathogenesis of AF has attracted more attention, various immune cells, including innate and adaptive immune cells, infiltrated the atrium and interacted with atrial cardiomyocytes or secreting several chemokines and cytokines to regulate the microenvironment of the heart ^[48]. However, regarding the specific subtype pAF, immune infiltration analyses are not abundant, thus with the help of CIBERSORT, we evaluated the expression level of 22 kinds of immune cells in pAF. The results showed that monocytes, dendritic cells (DCs) and neutrophils were upregulated and CD8+ T cells and T cells regulatory were downregulated.

Recent studies have confirmed the function of monocyte in atrial structural remodelling after catheter ablation of AF. The activation of monocyte especially the enhanced migration ability could lead to the enlargement of atrial structure ^[49]. DCs are essential for regulating innate and adaptive immunity owing to its antigen-presenting role, the activation of DCs were related to rheumatic heart disease ^[50] and systematic inflammation response after coronary artery bypass grafting ^[51]. A prospective clinical study found the elevated ratio of neutrophils to lymphocyte (NTL) had an intimate positive corelation with the incidence of AF. The role of CD8+ T cell remained controversial recently, Kazem N et al. ^[52] constructed a prospective clinical study and found the fraction of CD8+ lymphocytes was significantly higher in individuals developing postoperative atrial fibrillation (POAF) comparing with individuals free of POAF. While with the bioinformatic analysis, Li and his colleagues ^[53] found the expression of CD8+ T cell decreased in AF patients, which is consistent with our findings. We speculated that the contrary results might be determined by different timepoints of the AF onset. Whether the patient catch AF initially or after surgery might lead to different findings, which needs to be further explored. Moreover, the expression of T reg cells were reported to be prominently decreased in peripheral whole blood and right atrial tissue of AF patients ^[54].

In this study, we found that pAF patients had an elevated level of CXCR2 and monocytes. Zhang, Y.L. also found that CXCR2 plays a critical role in driving monocyte infiltration of the atria, causing atrial remodelling and AF after hypertension. This is consistent with our study ^[55]. In addition, S100A12 and neutrophils are both highly expressed in pAF patients. Therefore, we surmise that S100A12, serving as a biomarker of inflammation, may induce neutrophil infiltration, leading to atrial remodelling.

However, there were some limitations to our study. Firstly, we only concentrated on several important enrichment results and genes associated with the enrichment, which might neglect potential genes related to pAF. Secondly, for immune infiltration, further interaction between immune cells and miRNAs should be focused. Besides, the conclusions obtained from this paper should be investigated in vivo and vitro.

5.   Conclusions

We conjectured that miR-34a-5p might target RCAN1 and PPP3R1 in regulating pAF through Calcineurin-NFAT signaling pathway. S100A12 may drive neutrophil infiltration through induction of inflammatory response, leading to atrial remodeling. In addition, Monocytes, dendritic cells, neutrophils, CD8+ T cells and T cells regulatory (Tregs) may be involved in the pathogenesis of pAF.

6.   Highlight

1) We combined the LIMMA differential gene expression analysis with WGCNA in order to identify the hub genes and pathways related to the pathogenesis of pAF.

2) MiR-34a-5p may regulate RCAN1 and PPP3R1 in controlling pAF through the Calcineurin-NFAT signaling pathway, which provides a molecular mechanism and potential therapeutic targets for pAF.

3) S100A12 may lead to atrial remodeling through driving neutrophil infiltration.

4) Monocytes, dendritic cells, neutrophils, CD8+ T cells and T cells regulatory (Tregs) could potentially be developed as targets of immunotherapy in patients with pAF.

Conflicts of interest

The authors declare that they have no competing interests.

Acknowledgments

We would like to express our gratitude to the providers of the datasets GSE31821, GSE41177 and GSE79768 for sharing the data online. This work was supported by the grants from Xuzhou Science and Technology Bureau to Defeng Pan [grant number: KC20097].