Comprehensive Analysis of Immune Cell Infiltration and M2-Like Macrophage Biomarker Expression Patterns in Atrial Fibrillation

Introduction

Atrial fibrillation (AF) is a common arrhythmia worldwide, and its epidemiological characteristics and clinical manifestations have been extensively researched.^1,2 However, the treatment of AF still has some limitations, such as the potential adverse effects of antiarrhythmic drugs, the risk of bleeding associated with anticoagulation therapy, and the long-term efficacy of catheter ablation.³ These limitations are due to an incomplete understanding of the mechanisms of AF. Previous studies have indicated that electrical and structural remodeling of the atria underlies AF.⁴ Factors like ion channel abnormalities, cardiomyocyte remodeling, and cardiac neuromodulation disorders^5,6 contribute to reduced atrial contractile function and blood pooling, increasing the risk of thrombosis and stroke.⁷ The function of immune cells in atrial tissue is closely related to the progression of AF.⁸ To better comprehend the development of AF, further clarification of the remodeling mechanisms in AF is necessary at the cellular and molecular levels.^9,10

Chronic inflammation is considered a pathological pathway in several cardiac diseases, including AF.¹¹ Fibrosis, another consequence of inflammation, is closely linked to the structural remodeling of the atria.^12,13 Studies have shown that AF patients have higher levels of inflammatory markers compared to those in sinus rhythm (SR),^14,15 and the differentially expressed genes (DEGs) between AF and SR are primarily associated with inflammatory responses.^16,17 However, the molecular mechanisms driving fibrosis through inflammation are not fully understood. Increasing evidence suggests that immune cells play a crucial role in AF development,⁸ and the heterogeneity of monocytes and macrophages is closely linked to myocardial fibrosis and cardiac remodeling.^18,19 Macrophages, the most abundant leukocytes in the heart, are a major source of cytokines involved in various immune, infectious, endo-environmental, homeostatic, and metabolic processes. They also play a role in stimulating angiogenesis and fibrosis progression.^8,20–22 Macrophages are classified into two main phenotypes based on their polarization state: classical M1 (pro-inflammatory) and M2 (anti-inflammatory, tissue repair). M2-type macrophages contain multiple functionally distinct subpopulations, some of which exhibit features²³ that overlap with the M1 subset. M2-type macrophages have been suggested to be beneficial for AF due to their anti-inflammatory, apoptosis-promoting, and tissue-repairing properties.²⁴ However, their atypical pro-inflammatory,²⁵ pro-fibrotic^26,27 roles are often overlooked. A recent study using single-cell RNA sequencing (scRNA-seq) analysis found a higher infiltration of pro-fibrotic M2-like macrophages and CD8+ T cells in left atrial appendage (LAA) tissue of patients with AF compared to SR.²⁸ These results highlight the central role of inflammation and fibrosis in AF and suggest that macrophages recruited during atrial remodeling may serve as an immunotherapeutic target in patients with AF.

The main challenges in macrophage research involve identifying different populations and exploring the functions of marker genes in various phenotypes.²⁹ Studies at the Bulk RNA sequencing (Bulk RNA-seq) level have struggled to characterize intra- and inter-tissue heterogeneity, but recent advancements in ScRNA-seq have offered an unprecedented opportunity to systematically characterize specific subpopulations of macrophages and patterns of gene expression in the AF.³⁰ Macrophage localization and function in the heart are distinct, as demonstrated by Dick et al’s single-cell RNA sequencing, which identified four functionally distinct cardiac macrophage clusters.³¹ Maarten et al also identified five macrophage clusters through single-cell sequencing of left atrial tissue from patients with heart valve disease, and found that pro-inflammatory and pro-fibrotic SPP1+³² macrophage clusters were significantly expanded in atrial tissue from AF mice.³³ Integration of bulk RNA-seq and scRNA-seq analyses improved the accuracy of identifying diagnostic biomarkers and revealed the important role of macrophage subpopulations in the pathophysiologic process of AF at the cellular and molecular transcriptional levels.

Materials and Methods

Data Collection and Processing

AF and SR gene expression profiling data were downloaded from the Gene Expression Omnibus (GEO) platform (https://www.ncbi.nlm.nih.gov/geo/). Five bulk RNA-seq datasets containing RNA content of human LAA tissue (GSE14975, GSE31821, GSE79768, GSE41177, and GSE115574) were included for further evaluation. All of the above five datasets were based on the GPL570 platform and in total contained LAA tissue samples from 31 patients in SR and 44 patients with AF (Table S1). Using R (version 4.3.1) software, we merged the above five datasets, normalized them, performed batch effect removal using the “combat” function of the “sva” R package, and visualized them by Principal Component Analysis (PCA) plots. The scRNA-seq dataset was obtained from GSE224959 and contains LAA samples from 5 AF and 5 SR patients. The flow of data analysis in our study is shown in Figure 1.

Figure 1 Flowchart of the study. Abbreviations: AF, atrial fibrillation; SR, sinus rhythm; &, and; scRNA-seq, single-cell RNA sequencing; WGCNA, weighted gene correlation network analysis; CIBERSORT, cell type identification by estimating relative subsets of RNA transcripts; hdWGCNA, high-dimensional WGCNA; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes. ROC, receiver operating characteristic; LASSO, least absolute shrinkage and selection operator; qPCR, quantitative real-time PCR.

Identification of DEGs and the Functional Enrichment Analysis

The “limma” R package was used to perform differential gene expression analysis on the bulk RNA-seq data of the AF and SR groups, with a threshold of P < 0.05 and |log2 Fold change (FC)|> 1. The expression of DEGs was visualized using the “pheatmap” R package. The “ClusterProfiler” R package was employed for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of the DEGs, with a threshold set at <0.05.

Evaluation of Immune Cell Infiltration

The abundance of immune cell infiltration in AF patients in the bulk RNA-seq dataset was quantified in R using the cell-type identification by estimating relative subsets of RNA transcript (CIBERSORT) algorithm.³⁴ We filtered the samples using P<0.05 to obtain infiltration ratios and correlations for 22 immune cell subtypes. The percentage of each immune cell type in the samples was quantified using the “ggplot2” R package, and the “brewpub” R packages are shown in the bar graphs. “Corrplot” package was used to generate heatmaps related to 22 infiltrating immune cells. The “vioplot” package was used to visualize the infiltration level of 22 immune cells between AF and SR groups, and P<0.05 was considered to be different.

Constructing the Gene Co-Expression Network

Using the “WGCNA” R package to identify gene modules associated with M2 macrophages. The top 10,000 genes with the highest absolute median difference in the merged bulk RNA-seq data were chosen, converted into a neighbor-joining matrix, and clustered. Outlier samples were removed to ensure the reliability of the network construction results. The “pickSoftThreshold” function was utilized to select and validate the best soft threshold (β) value. Hierarchical clustering dendrograms were then constructed (with a minimum number of genes per module set to 30 and a mergeCutHeight of 0.25), and similar modules in the clustering tree were merged using module eigengene (ME) values. For regulators associated with immune infiltrating cells, module membership (MM) and gene significance (GS) were recalculated and marked with different colors. The module with the smallest P-value and highest relevance was chosen as the key module, and the genes within the key modules were further analyzed.

Cell Clustering Analysis

The “Seurat” R package was used for downstream analysis of scRNA-seq data, with quality control (QC) requirements of >200 total unique molecular identifiers (UMIs) in a single cell, 250 < nfeatureRNA < 2500, and mitochondrial genes < 10%. Data were normalized using LogNormalize, the first 2000 highly variable genes were screened using the vst method, data were downscaled using PCA combined with highly variable genes with Principal Components (PCs) of 20, and batch effects were removed using harmony. The Louvain algorithm in FindNeighbors and FindClusters was used for intercellular clustering of cells, with a resolution of 0.8 for large classes of cells and 0.4 after subclustering, and visualized by the Uniform Manifold Approximation and Projection (UMAP) method. The selection criteria for DEGs were that the logFC was greater than 0.25 after subfamily clustering of medullary-like cells, and only positive values were selected. Subsequently, we used the CellMarker database (http://bio-bigdata.hrbmu.edu.cn/CellMarker/index.html), the “SingleR” R package, and manual hand annotation to identify cell types. The “ggplot2” package was used to show the compositional ratio of different cell types between the AF and SR groups.

Analysis of Cell-Cell Communication

To further explore the interaction and function of macrophage subtypes with other immune cells in AF, we conducted cell communication analysis. We loaded Seurat objects into the Cell Chat program and utilized the “aggregateNet” function to analyze the interaction counts and weights among cell clusters. The “netAnalysis_signalingRole_heatmap” function was employed to identify the outgoing or incoming communication signals of the aggregated cell-cell communication network. Significant cell-cell interactions mediating cell-cell communication (ligand-receptor pairs in the bubble diagram) were identified by assessing the cell-cell communication networks. Identify cell-specific signaling interactions using “network centrality analysis”.

hdWGCNA Analysis of Macrophages in AF

hdWGNCA can construct co-expression networks across spatial hierarchies and is more suitable for scRNA-seq data.³⁵ To further identify the characteristic gene modules of M2 macrophages in the AF group, the “hdWGCNA” R software package was utilized to conduct weighted gene co-expression network analysis of the scRNA-seq dataset. The gene expression data were first normalized, and then Pearson correlation coefficients were calculated to capture the co-expression relationship between each pair of genes. Parameters were scanned for various soft power thresholds (from 1 to 30) using the “TestSoftPowers” function. The “PlotSoftPowers” function was then used to visualize the optimal soft threshold for the scale-free network. Based on these thresholds, an adjacency matrix and a topological overlap matrix (TOM) were generated. The “dynamicTreeCut” algorithm was employed to identify co-expression modules from the TOM and assign a unique color to each module. The connectivity values (kME) of the feature genes were calculated using the “ModuleConnectivity” function to assess their core gene connectivity within the respective modules. The top 25 feature gene scores in each module were obtained using Seurat’s Ucell function, and a feature plot was constructed for each module using the “ModuleFeaturePlot” function. The “GetModuleTraitCorrelation” function was used to determine the importance of the modules, and the “corrplot” R package was used to visualize the correlation between each module. The average expression of module-specific core genes in different macrophage clusters was visualized by Seurat’s “DotPlot” function. Finally, the top 50 genes in each module were extracted as core genes for further analysis.

Identify Key Genes as AF Biomarkers and Evaluate Their Diagnostic Value

To identify genes highly correlated with M2 macrophage infiltration in AF, the Venn diagram was used to take intersecting genes, including DEGs in the AF and SR groups of the bulk RNA-seq dataset, genes highly associated with M2 macrophages in WGCNA, differentially expressed genes in M2 macrophages in scRNA-seq, and hdWGCNA Genes. Two machine learning algorithms, LASSO and random forest, based on the “glmnet” and “randomForest” R packages, were used to screen the genes that could identify AF. Binary logistic regression models were constructed based on the characterized genes using the “glm” R package and the “caret” R package. Validation of model stability using external dataset GSE2240. Additionally, to more accurately assess the diagnostic validity, an automated search for the optimal parameter settings was performed based on the best predicted results using a five-fold cross-validation method. The `pROC` R package was used to plot ROC curves to assess the diagnostic effectiveness of the model by calculating AUC values.

Clinical Sample Collection and Reverse Transcription Quantitative Polymerase Chain Reaction (RT-qPCR) Analyses

To further validate the reliability of the Hub genes, we recruited 6 AF and 6 SR patients between 01 November 2023 and 15 November 2023 Information on the basic characteristics of the patients is shown in Table S2. All subjects had their blood collected in the morning on an empty stomach, and venous blood was collected immediately after using EDTA tubes to centrifuged to separate plasma and stored at −80°C for subsequent experiments. The study complied with the Declaration of Helsinki, and this clinical blood collection protocol was approved by the Ethics Committee of Yunnan University Hospital [No. 2022200], and informed consent was provided in writing to the recruited subjects, (Appendix 1).

An online tool (https://www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi?link_loc=blasthome) was used to design primers and assess primer specificity. Primer lengths were set between 20bp and 25bp and β-actin was selected as the internal reference gene. Total RNA was extracted from fresh plasma obtained from the AF (n=6) and SR (n=6) groups using TRIzol (ambion) according to the manufacturer’s protocols. The Hifair^® III Ist Strand cDNA Synthesis SuperMix for qPCR (gDNA digester plus) kit was then used to The RNA was reverse transcribed into cDNA. The total volume of fluorescent quantitative PCR amplification was 20 ul. The experimental rows were replicated three times in parallel, and the relative quantitative levels of mRNA for each subgroup were calculated using the 2-∆∆CT method. β-actin was used as an internal control. The primers used in this study are shown in Table S3.

Construction of a Multifactor Regulatory Network Based on Key Signatures

To further explore the expression patterns of hub genes, experimentally supported databases of functional target regulators, miRTarBase (https://mirtarbase.cuhk.edu.cn/~miRTarBase/miRTarBase_2022/php/index.php) and Regnetwork (http://www.regnetworkweb.org), were used to predict microRNAs (miRNAs) and transcription factors (TFs) targeting hub genes. Results with a p-value ≤ 0.05 were chosen as critical values. After obtaining the TF-miRNA-mRNA regulatory relationships, Cytoscape was used to visualize the regulatory network.

Single-Gene GSEA Enrichment Analysis

To further explore the potential signaling pathways of Hub genes in AF, we conducted single-Gene Set Enrichment Analysis (GSEA) enrichment analysis using the merged BulkRNA-seq dataset, based on the “c2.cp.kegg.v7.0.symbols” gene set from the MSigDB database. In order to obtain the standardized enrichment scores for each analysis, gene set permutations were set to 1000-fold, and a p-value < 0.05 was considered as a significantly enriched pathway.

Statistics

All analyses in this study were conducted using R software 4.3.1 and GraphPad Prism 9.0. Group comparisons were carried out using t-tests. The accuracy of the prediction model was evaluated using ROC curves. Fisher’s exact test was employed to analyze scRNA-seq cell proportions, and correlation analysis was used to examine the correlation coefficients between the two variables. Statistically significant differences were defined as p-values less than 0.05.

Result

DEGs Identification and Functional Enrichment Analysis

After merging the five bulk RNA-seq datasets, the Results of removing the batch effect are presented by a two-dimensional PCA clustering diagram, as shown in Figure 2A, which shows the samples of the AF and SR groups clustered together separately. We identified a total of 593 DEGs in the AF group, of which 387 were up-regulated and 216 were down-regulated. The volcano plot demonstrated the top 10 genes that were up- or down-regulated (Figure 2B). Figure 2C depicts the results of GO enrichment analysis based on DEGs, this includes processes of cellular ketone metabolic process, transmembrane receptor protein serine/threonine kinase signaling pathway, hormone metabolic process, collagen-containing extracellular matrix, platelet alpha granule, extracellular matrix structural constituent, and heparin binding in relation to metabolism, oxidation, coagulation, matrix related biological functions. KEGG pathway analysis of DEGs showed that DEGs were significantly enriched in AF-related immune, inflammatory and fibrotic pathways, such as PI3K-Akt signaling pathway, TGF-β signaling pathway, rheumatoid arthritis, phagolysosomes, PPAR signaling pathway, ECM-receptor interactions, and Th1 and Th2 cell differentiation (Figure 2D). These biological functions are intricately linked to the primary mechanisms of AF, such as triggered activity and matrix remodeling.

Figure 2 Identification of DEGs between AF and SR samples and enrichment analysis of DEGs using the GO and KEGG databases. (A) Principal component analysis clustering of gene expression in AF and SR tissues. (B) Volcano plots for the genetic differential analysis. The red points represent upregulated DEGs, and the blue points denote downregulated DEGs. (C) The circled graph on the left shows the relationship between key DEGs and the most enriched biological processes, while the graph on the right illustrates the top 6 enriched BPs, CCs, and MFs. (D) Functional enrichment analysis of the KEGG database. The words on the left indicate enriched KEGG, the size of the balls indicates the number of genes enriched, and the color indicates the level of enrichment. Abbreviations: DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; AF, atrial fibrillation; SR, sinus rhythm; BPs, biological process; CCs, cellular component; MFs, molecular function.

The Difference of Immune Cell Infiltration

Based on the enrichment analysis results suggesting that DEGs were concentrated in immune-related signaling pathways and responses, we further applied the CIBERSORT algorithm to assess the immune cell infiltration in the two groups of specimens. The bar graphs demonstrated the relative percentages of 22 immune cell subtypes in 31 SR samples and 44 AF samples (Figure 3A), and we found that there was a significant difference in immune cell infiltration between SR samples and AF samples. The correlation heatmap (Figure 3B) showed that M2 macrophages showed the strongest positive correlation with monocytes (r=0.51), while M2 macrophages showed a significant negative correlation with M0 macrophages (r=−0.82). Figure 3C shows that the infiltration of activated CD4 memory T cells, M1 macrophages, M2 macrophages, and activated mast cells was higher in the AF group compared with the SR group, while naive B cells, resting dendritic cells, and resting mast cells were relatively lower. Our findings suggest that these 7 immune cells are involved in the progression of AF, in particular, M2 macrophages show a markedly active expression state in AF and the strongest correlation with monocytes and M0 macrophages. This finding is consistent with previous findings²⁸ that activated M2 macrophages play an important role in the onset and progression of AF. This finding deepens our understanding of proinflammatory, profibrotic M2-type macrophages in AF and provides key clues for future studies.

Figure 3 Visualization and assessment of immune cell infiltration. (A) The relative percentages of the 22 immune cell types in 75 left atrial appendage samples. Each color represents one cell type. (B) Heatmap showing the correlation of the 22 immune cell types. Red signifies a positive correlation, whereas blue indicates a negative correlation. (C) The fraction of infiltrating immune cells in the AF and SR groups. The SR group is denoted by the color blue, while the AF group is denoted by the color red. A significance level of P<0.05 was used to indicate statistical significance. Abbreviation: NK, natural killer.

Identification of Gene Co-Expression Networks and Modules

We went to further explore gene modules highly associated with M2 macrophage subtypes by WGCNA methods, The optimal β-value is 10 (Figure 4A), and the resulting scale-free topology fit index was calculated to be 0.81, meeting the criterion for scale-free networks. Using dynamic shear tree analysis, we created co-expression modules, which were then merged to form a total of 4 modules (Figure 4B). The sample trait heatmap displayed the differential distribution of the seven immune infiltrating cell types associated with AF across the samples (Figure 4C). We then ranked the module feature matrices and calculated the correlation between the model feature matrices and the infiltrating immune cells to generate a heat map of the module-trait relationships (Figure 4D). The results showed that the black module was significantly positively correlated with M2 macrophages (r = 0.64, p < 0.001) and contained 3288 genes associated with M2 macrophages (Figure 4E).

Figure 4 WGCNA screening of M2 macrophage-associated gene modules. (A) Evaluation of the scale-free fit index and mean connectivity across different soft-thresholding powers. (B) Associations between consensus modules and samples, each comprising a group of closely connected genes. Each branch represents an individual gene, with each color indicating a co-expression module. (C) Trait heatmap displaying the distribution of the seven immune cell types in each sample. (D) Heatmaps showing the correlations between 4 modules and 7 types of immune cells. (E) Scatter plots describing the relationship between gene significance and gene module membership in the black module. Each dot represents a gene in the black module. Abbreviation: WGCNA, weighted gene coexpression network analysis.

Single-Cell RNA Sequencing Data Preprocessing and Analysis

The results of the scRNA-seq raw data QC are shown in Figure S1A, where harmony further integrated and eliminated batch effects (Figure S1B). Based on the Seurat standard steps, we finally identified 21 clusters, which were then embedded into the UMAP graph as shown in Figure 5A. A dataset containing 18,856 cells in 5 AF samples and 5 SR samples was further analyzed. Cell annotation was performed by manual annotation combined with the CellMarker database and the “SingleR” R package, and we identified eight cell types, including Tissue stem cells (SOX2, OCT4), T cells (CD3, CD4), Natural Killer cells (NK cells) (CD56, CD16), Endothelial cells (CD31, CD34), Myeloid cells (CD11b, CD14), Smooth muscle cells (SMA, calponin), Neutrophils (CD15, CD16), and B cells (CD19, CD20). The UMAP plot shows the distribution difference of these cell types between AF and SR groups (Figure 5B), and the results of cell types annotated by classical genes are shown in Figure 5C, which demonstrates that the results of our cellular annotation are clear and accurate. The cell ratio stacking diagram (Figure 5D) further illustrates that the number of Myeloid cells was significantly increased in the AF group compared with the SR group, while the number of Smooth muscle cells, Endothelial cells, NK cells, and T cells was significantly decreased. We found that Myeloid cells were the cell type that focused on the changes in the AF group, which is consistent with our previous expected results. We further annotated myeloid cell subsets, generating UMAP plots (Figure 6A) and bar graphs (Figure 6B). The differentially expressed genes among myeloid cell subpopulations are shown in Figure 6C, CTSB, MSR1, CTSA, CD36, CD68, FCGR2A, MAFB, and CD14 marker genes were highly expressed in the M2 Macrophage cell population, Up-regulation of these marker genes indicates the polarization status of M2 macrophages and their immunomodulatory role in AF. Figure 6D shows that the AF group had higher proportions of DC cells, macrophage M2, and macrophage M0, and lower proportions of macrophage M1 and monocytes compared with the SR group. Consistent with the findings in our bulk data of the active expression status of M2 macrophages in AF, we hypothesized that the over infiltration of M2 macrophages might be an important factor in the progression of AF.

Figure 5 Heterogeneity of immune cells in the AF and SR groups. (A) UMAP plot showing the cell distribution of the AF and SR groups. (B) UMAP plot showing the results of the cell annotation between the AF and SR groups. (C) Bubble dot graphs showing the expression of the top three different marker genes for each cell type. The colors of the dots represent the average expression, and the sizes of the dots represent the average percentage of cells that expressed the DEGs. (D) Stacked plot of cell proportions between different groups and each sample. Abbreviations: AF, atrial fibrillation; SR, sinus rhythm; UMAP, Uniform Manifold Approximation and Projection; DEGs, differentially expressed genes.

Figure 6 Cell annotation for myeloid cell subsets between groups in the AF and SR cohorts. (A) UMAP plot showing the distribution of myeloid cell subsets in the AF and SR groups. (B) Bar graphs showing the proportions of cell clusters in the different groups. (C) Heatmap showing the DEGs in the myeloid cell subsets. Yellow for high expression, purple for low expression. (D) Distribution of myeloid cell subsets in each group was determined using a UMAP plot and a stacked plot of myeloid cell subset proportions between the different groups. Abbreviations: DEGs, differentially expressed genes; AF, atrial fibrillation; SR, sinus rhythm; DEGs, differentially expressed genes; UMAP, Uniform Manifold Approximation and Projection.

Cell-Cell Communication Analysis

In the AF group, we used the“CellChat”R package to explore intercellular communication and to describe the signaling pathways that interacted in the single-cell data. Figure 7A shows a significant increase in interactions between macrophage M1, macrophageM2 and DC cells in AF. As shown in Figure 7B, macrophage subtypes, DC cells, and smooth muscle cells represented the cell types with the highest number of incoming and outgoing signaling patterns in AF. The strengths of cell output and afferent interactions are plotted in Figure 7C. To further elucidate the potential cell-cell communication between macrophage M1 and macrophage M2 and other cell types of cells, based on the relative expression of ligand-receptor (L-R) pairs, we analyzed the output signaling of macrophage M2 and macrophage M1 (Figure 7D). At the level of cell-cell contact, macrophage M0 and DC was most active in the TNF-TNFRSF1A ligand-receptor pair, affecting macrophage M1, macrophage M2 and monocytes by activating the TNF signaling pathway (Figure 7E). Smooth muscle cells were the only senders and receivers in the ANGPTL signaling pathway network (Figure 7F).

Figure 7 Cell‒cell communications in AF. (A) Integrated cell‒cell communication network plotted by interaction and weight. The circle sizes are proportional to the number of cells in each cell group, and the edge width represents the communication probability. (B) Heatmaps showing the relative strength of the signaling pathways among the 11 cell populations in the outgoing and incoming signaling patterns. A gradual change in color from green to white indicates a change in the relative interaction strength from high to low. (C) Dot plot of outgoing/incoming interaction strength for 11 cell clusters. Dot size is positively correlated with strength. (D) Bubble plot of macrophage subsets outgoing interaction signaling pathways. The dot color reflects the communication probabilities, and the dot size represents the computed p value. (E and F) The relative contribution of TNF-α and ANGPTL signaling to each cell group. The term “sender” refers to a source of signaling, “receiver” refers to the target of signaling, “mediator” refers to the gatekeeper of a cell, and “influencer” refers to the ability to influence the flow of information within a signaling network. Abbreviation: AF, atrial fibrillation.

hdWGCNA Analysis of Macrophage Subsets in AF

To identify genes characterizing macrophage subpopulations in the LAA tissue of patients with AF, we performed hdWGCNA. The optimal soft threshold for constructing the scale-free network was 5 (Figure 8A). Dynamic shear dendrograms were visualized for hierarchical clustering based on gene expression patterns, and nine gene modules were generated (Figure 8B). Figure S2A and B shows the top 10 hub genes of these 9 modules and illustrates the relationships between the modules. In addition, the expression distribution of the hub genes from these nine modules was mapped onto macrophage subsets (Figure 8C). Figure 8D shows that blue, yellow and pink modules are positively correlated with macrophage M2, macrophage M0 and monocyte cluster.

Figure 8 hdWGCNA identifies the modular signature genes of macrophage subsets in AF. (A) Soft power = 5 was selected to construct the scale-free network. (B) Hypervariable genes were clustered by hdWGCNA into 9 modules. Each leaf in the dendrogram represents an individual gene, while the color assigned at the bottom indicates its membership in a specific co‐expression module. The “gray” module consists of genes that were not grouped into any co‐expression module. (C) t-SNE plot showing the expression distribution of hub genes for each module across the 9 clusters. Darker colors represent higher module gene expression. (D) Dot plot presenting the average expression of module-specific hub genes in different macrophage subsets. The size of the dots represents the percentage of cell subtypes and the color of the dots represents the average expression of the module genes. Abbreviations: hdWGCNA, high-dimensional weighted gene coexpression network analysis; AF, atrial fibrillation; t-SNE, t-stochastic neighbor embedding.

Identification of M2 Macrophages Related Biomarkers in AF

Using the Venn diagram in the multidimensional data, we found five DEGs associated with M2 macrophages in AF (Figure 9A). We further screened for signature genes with the ability to diagnose AF by applying machine learning algorithms. After selecting the best λ, the LASSO algorithm screened four genes (IGF1, PDK4, RAB13 and TMEM176B) (Figure 9B and C). The relative importance of the five candidate genes in the random forest algorithm is shown in Figure 9D. Eventually, the LASSO and random forest algorithms shared the Four genes, IGF1, PDK4, RAB13, and TMEM176B, were identified as biomarkers of AF. We developed a binary logistic regression model using the four mentioned biomarkers. The AUC curves in Figure 9E showed that the four biomarkers had high predictive validity (AUC > 0.75). We validated the stability of the model using the external dataset GSE2240 from the GPL96 and GPL97 platforms, respectively. As shown in Figure S3A and B, the data validation results (average AUC = 0.654) from the two platforms complement each other. Although the results of the ROC curves did not fully meet our expectations, these data still showed certain trends that supported our initial findings. Additionally, we assessed the reliability of the model’s efficacy using 5-fold cross-validation. The results indicated that the model demonstrated good and reliable predictive efficacy (AUC ≥ 0.8) with a mean value of 0.8814 (Figure 9F).

Figure 9 Machine learning algorithms to identify Hub genes and ROC curve validation. (A) The feature gene selection from different dimensions. Set1 includes the genes highly associated with M2 macrophages in WGCNA; Set2 consists of the DEGs from the bulk RNA-seq dataset; Set3 comprises the differentially expressed genes of M2 macrophages in the scRNA-seq; Set4 contains the genes highly correlated with the subpopulation of M2 macrophages in the hdWGCNA. (B and C) Adjustment for feature selection by LASSO logistic regression analysis. (D) Genes are displayed in descending order of importance by random forests. (E) The diagnostic efficacy of the 2 algorithms for the 4 crossover genes is represented by ROC curves. (F) Model stability was assessed using fivefold cross-validation. Abbreviations: WGCNA, weighted correlation network analysis; DEGs, Differentially Expressed Genes; hdWGCNA, high-dimensional weighted gene co‐expression network analysis; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic.

Validation of the Key Genes

To further validate the reliability of our diagnostic genes, we collected serum samples from 6 pairs of AF-SR patients in the clinic and examined the expression levels of each of the 4 key genes by qRT-PCR to verify their differential expression between AF and SR. Similar to the results of previous analyses (Figure 10A), the expression levels of IGF1, RAB13 and TMEM176B were higher in the AF group than in the SR group (Figure 10B–D), whereas the expression of PDK4 was not significantly different between the 2 groups (Figure 10E).

Figure 10 Validation of the key genes in AF patients through qPCR (A) Expression levels of the Hub genes in the merged Bulk-seq dataset. (B–E) The relative expression levels of key genes in whole blood between the SR and AF groups. (B) IGF1; (C) RAB13; (D) TMEM176B; (E) PDK4. *P < 0.05, ***P < 0.001. Abbreviations: AF, atrial fibrillation; qPCR, quantitative real-time PCR; SR, sinus rhythm; ns, not significant.

TF-miRNA-mRNA Regulatory Network Analysis

To further explore the potential regulatory mechanisms of Hub genes, 152 TF-miRNA-mRNA regulation were obtained after predicting miRNA-gene pairs and TF-gen pairs gained relationships. Regulatory network (including 117 miRNAs, 41 transcription factors and 4 target mRNAs) visualized by Cytoscape (Figure 11). In this network, we found that miRNAs and TFs targeting IGF1 were the most abundant, and IGF1 and PDK4 were targeted by activated TF factor (STAT5A) and hsa-miR-152. In addition, RAB13 and IGF1 were targeted by activated TF factors (JUN and CTCF), while hsa-miR-19a and hsa-miR-19b were bridges connecting RAB13 and IGF1.

Figure 11 The TF-miRNA‒mRNA regulatory network. The red balls depict mRNAs, the blue balls represent miRNAs, and the brown balls represent TFs. Abbreviations: mRNAs, messenger RNAs; miRNAs, microRNAs; and TFs, transcription factors.

Signaling Pathways Analysis of Macrophage-Related Biomarker

In order to explore the biological pathways of Hub genes in AF pathology, we performed single gene GSEA enrichment analysis. The results showed that IGF1, RAB13, TMEM176B and PDK4 were enriched in Circadian rhythm, Fatty acid metabolism, and PPAR signaling pathway, Phagosome, ECM-receptor interaction, and Endocytosis. (Figure 12A–D).

Figure 12 Identification of signaling pathways affected by target genes via GSEA. The major signaling pathways were enriched based on the expression of target genes in the merged Bulk RNA-seq dataset. (A) IGF1; (B) RAB13; (C) TMEM176B; (D) PDK4. Abbreviation: GSEA, Gene Set Enrichment Analysis.

Discussion

Trigger and atrial matrix remodeling are thought to be the main causes of initiating and sustaining AF, increasing the likelihood of automaticity and re-entry, and are closely related to immune remodelling.⁹ Inflammation and fibrosis are considered significant elements in the immune remodeling process of AF.^36,37 Monocytes/macrophages are key inflammatory cells in the heart, and their contribution to repair and fibrosis in the early stages of cardiac injury is well documented.^38,39 Macrophages are heterogeneous populations of immune cells recruited and activated during different cardiac processes, and exhibit diverse functions.^40–42 Recent studies have shown that macrophage heterogeneity is closely associated with myocardial fibrosis, cardiac remodeling, and may serve as a potential therapeutic target for atrial fibrosis.¹⁸ However, their specificity in the heart remains controversial. M2 macrophages have been implicated in myocardial fibrosis because of their anti-inflammatory, pro tissue repair, remodeling, angiogenesis, and maintenance of homeostasis, but their deleterious aspects are not fully understood. IL-4 (M2 activator) has been reported to induce M2 macrophages to exhibit a nonclassical proinflammatory phenotype.²⁵ A study by Westermann et al reported an increased number of TGFβ-expressing leukocytes with characteristics of M2 macrophages in patients with heart failure in whom the ejection fraction remained unchanged.⁴³ During ischemia and reperfusion, the transition from M1 to M2 stimulates a pro-fibrotic environment, including maturation of monocyte precursors into fibroblasts.⁴⁴ These results illustrate the evidence for activation of M2 macrophages in disease pathogenesis. Our study comprehensively analyzed a large amount of RNA-seq and single-cell RNA-seq data to characterize immune cell infiltration in LAA tissue from patients with AF. Using multiple machine learning algorithms and combining them with experimental validation, we identified biomarkers associated with M2 macrophages in AF and speculated on the potential molecular regulatory mechanisms of these biomarkers in the pathophysiological process of AF.

Our enrichment analysis of differential gene lines in AF revealed that these DEGs were significantly enriched in the TGF-β signaling pathway, PPAR signaling pathway, and ECM-receptor interactions.⁴³ It is worth noting that TGF-β is an important pro-fibrotic factor secreted by macrophage M2,⁴⁴ and that macrophage migration inhibitory factor (MIF) promotes cardiac fibrosis in patients with AF through activation of the TGF-β signaling pathway.⁴⁵ Excessive M2 macrophage activity may promote TGF-β1 and IL-10 secretion, thereby inhibiting ECM uptake and degradation and cardiac fibrosis in patients with AF.⁴⁵ Excess M2 macrophage activity may promote the secretion of TGF-β1 and IL-10, which inhibits ECM uptake and degradation and increases susceptibility to myocardial fibrosis and AF.⁴⁶ PPAR regulates macrophage metabolism, and multiple signaling proteins such as PPARα, PPARc, and SIRT3 have been shown to be Ang II–induced atrial an important target of inflammation and fibrosis.^47,48 These results suggest that in AF, an imbalance of M2 macrophages may lead to the activation of multiple inflammatory, fibrotic, and metabolic-related signaling pathways. Understanding pro-inflammatory and pro-fibrotic macrophage subtypes and their gene functions may help further elucidate the complex mechanisms of electrical, structural and neural remodeling in heterogeneous atrial conduction. By the CIBERSORT algorithm, we found that T cells and macrophage subpopulations showed high variability between AF and SR groups. In addition, consistent with previous findings,⁴⁹ M2 macrophages and T cells were increased in the AF group. The new study demonstrated a higher infiltration abundance of the M2 typical marker CD163 in the LAA of patients with AF by quantitative analysis.²⁸ These results characterize the immune cell infiltration in the LAA tissue of patients with AF. Second, by WGCNA analysis, we found that M2 macrophages correlated most significantly with the black module. From this, we hypothesized that the black module gene is a M2 macrophages-related gene in AF progression.

By scRNA-seq analysis, we found that the number of Myeloid cells was significantly increased in the AF group, while the number of Smooth muscle cells, Endothelial cells, NK cells, and T cells was decreased. To further clarify the heterogeneity of cell types, we performed subclustering analysis of myeloid cells. The results showed that DC cells and M2 macrophages had high variability between AF and SR groups. Similar to the results of the immune infiltration analysis in the bulk data, the proportion of M2 macrophages was increased in the AF group. The results of the study by Chris J. Watson showed that pro-fibrotic CD163+ M2 macrophages, collagen amount, and BNP expression were higher in atrial tissue of AF patients compared with SR patients.⁵⁰ In addition, a study by Sheng et al found that M2-like macrophages with high pro-inflammatory scores and CD8+ T cells were recruited in the left atrium of the AF left atrium.²⁸ The above results suggest that that the immune function and related molecular features of macrophage M2 in AF progression deserve further exploration. Intercellular communication can not only elucidate the underlying intercellular communication and cellular states, but also uncover the intracellular conduction and amplification of specific signaling pathways.^51,52 We found that the pattern of intercellular interactions in LAA tissues of patients with AF was predominantly characterized by cell-substrate interactions. At the level of cellular ligand-receptor-mediated signaling, we found that TNF-TNFRSF1A signaling was most active in macrophage M0 and DC cells, which interacted with macrophage M1, macrophage M2, and monocytes through the TNF signaling pathway. Notably, TNF, an inflammatory, fibrosis-related signaling pathway involved in the progression of AF, was found to be activated to varying degrees in different subtypes of macrophages, and alterations in TNF/TNFRS1A signaling may increase proinflammatory signaling.^53,54 In addition, we found that smooth muscle cells are the sole senders and receivers in the ANGPTL signaling network. The TGF-β1/Smad pathway in smooth muscle cells was found to be an important pathway affecting atrial remodeling and fibrosis.⁵⁵ Animal experiments have shown that ANGPTL4 attenuates Angii-induced atrial fibrillation and fibrosis through the PPAR pathway.⁴⁸ Our results emphasize the important pathway of M2 macrophage-modulated immune-inflammatory responses in AF.

We identified IGF1, PDK4, RAB13, and TMEM176B as biomarkers associated with M2 macrophages in AF using two machine-learning algorithms, LASSO and Random Forest. These biomarkers were experimentally validated through qPCR. We found that the relative expression levels of IGF1, TMEM176B and RAB13 were consistent with the expression trend of the analyzed results, while PDK4 was not significantly differentially expressed in AF. Upon literature review, we found evidence supporting PDK4 as a valid biomarker for distinguishing AF from SR by analyzing whole blood gene expression before and after electrical cardioversion (ECV) in AF patients.⁵⁶ IGF-1 is synthesized and secreted primarily in the liver, and it has been reported to be produced by cardiac fibroblasts that are also exposed to pressure overload stimulation in the murine heart The expression of IGF-1R mRNA in cardiac fibroblasts is higher than that in myocardial fine, and the effect of IGF-1 on cardiac fibrosis may vary depending on the duration and dose of treatment.¹⁹ PDK4 is a key enzyme in the regulation of glucose homeostasis.⁵⁷ The activation of NF-kappaB in cardiac cells inhibits PPAR activity, which leads to decreased PDK4 expression and increased glycolysis.⁵⁸ PDK4 is a regulator of fibrogenesis⁵⁹ and increased glycolysis induces myocardial fibrotic remodeling.⁶⁰ RAB13 is a protein involved in organelle-to-organelle communication and plays a key role in intracellular processes such as vesicle transport and phagocytosis.⁶¹ The researchers analyzed the mRNA expression dataset in the GEO database to identified 15 biomarkers of septic shock including RAB13.⁶² These results suggest that RAB13 is closely associated with phagocytosis and inflammation. TMEM176B is a member of the membrane-spanning 4-domain (MS4) family of transmembrane proteins, and TMEM176 and the enzyme-like membrane-bound protein CD39 may be potential targets for unleashing the activation of the inflammasome.⁶³ Marcelo Hill’s study showed that TMEM176B acts as a dual immunomodulator that can inhibit effector immune responses in some settings while promoting immunity through antigen presentation in others.⁶⁴ However, the roles of RAB13 and TMEM176B in atrial fibrillation have never been investigated. We hypothesize that IGF1, PDK4, RAB13, and TMEM176B may influence the progression of AF through M2 macrophage infiltration, and more in-depth studies are needed to confirm the important roles and pathways of these genes.

The expression of TFs and their regulatory activities are important for resolving complex life activities, and dysregulation of gene expression has been associated with disease.⁶⁵ Alterations in the expression levels of miRNAs cause adverse effects such as cardiomyocyte damage, cardiac fibrosis, angiogenesis, and inflammatory responses.⁶⁶ Dysregulation of miRNA and TF expression in early stages underlies the development of AF.⁶⁷ Among the 4 target genes forming the TF-miRNA-mRNA regulatory network, we found that miRNAs and TFs targeting IGF1 were most abundant, and IGF1 and PDK4 were targeted by activated TF factor (STAT5A) and hsa-miR-152. In addition, RAB13 and IGF1 were targeted by activated TF factors (JUN and CTCF), whereas hsa-miR-19a and hsa-miR-19b were bridges connecting RAB13 and IGF1. Lu et al demonstrated that up-regulated hsa-miR-9, and hsa-miR-152 prevented intracellular Ca(2+) concentration by regulating the concentration of AF.⁶⁸ Evidence that miR-19a and miR-19b of the miR-19 family regulate cardiomyocyte proliferation and inhibit apoptosis.⁶⁹ We hypothesized that hsa-miR-152/IGF1 and hsa-miR-152/PDK4 might be potential therapeutic targets for AF; the activated expression of RAB13 and IGF1 might be caused by the activated expression of JUN,⁷⁰ which mediates the transcription of TGF-β, and CTCF,⁷¹ which regulates the state transition of immune cells; and the silencing of hsa-miR-9 and hsa-miR-152 might contribute to the expression of RAB13 and IGF1, whereas blockade of this regulation might delay the progression of AF. Many of experimental studies will be needed in the future to validate the rules and pathways proposed here.

To further illustrate the biological pathways of the core genes in the pathological process of AF, we performed single-gene GSEA enrichment analysis. The results showed that four core genes were enriched in Circadian rhythm, Fatty acid metabolism, PPAR signaling pathway, Phagosome, ECM-receptor interaction, and Endocytosis. The results of the analysis drew our attention to fatty acid metabolism and PPAR. Previous studies have found that plasma free fatty acid may be a potential biomarker for predicting outcome events in stroke patients with AF in conjunction with the CHADS 2 and CHA 2 DS 2 -VASc scoring systems.⁷² Pericardial fat is a rich source of inflammatory mediators, leading to structural remodeling of the atria and fibrosis by contributing to the onset and persistence of AF.⁷³ In this study, we examined the role of collagen remodeling and fatty acid metabolism in the development and persistence of AF. Correlated biomarkers between collagen remodeling and dysregulation of fatty acid metabolism are involved in structural remodeling in AF.^74,75 PPAR is an important signaling pathway that regulates macrophage metabolism and promotes M2-phenotyping.⁷⁶ PPAR regulates the expression of a variety of genes involved in myocardial energy metabolism to adapt to various pathophysiological conditions in the heart.⁷⁷ Some studies have revealed that AF is a major cause of cardiovascular disease by means of transcriptomics and proteomic histology analyses revealed that the expression of various genes and proteins in the atrial tissues of patients with AF with VHD was closely related to the PPAR signaling pathway.⁴⁷ Combined with our results, we hypothesized that IGF1, RAB13, TMEM176B, and PDK4 may regulate the progression of AF through the PPAR signaling pathway and fatty acid metabolism.

Our study comprehensively analyzed the cellular heterogeneity of AF and speculates on the expression patterns of AF gene markers associated with M2 macrophages. However, due to the limited number of human LAA specimens in public databases, increasing the sample size further may enhance the credibility and authenticity of the findings. Additionally, expanding the study into multi-omics fields such as proteomics, metabolomics, and epigenetics could contribute to a more comprehensive understanding of the disease process and uncover the underlying molecular mechanisms of AF. Finally, the putative regulatory networks and signaling pathways based on target genes need to be validated through additional experimental studies.

Conclusion

In Conclusion, our findings suggest that over infiltration of M2 macrophages may be an important factor in the progression of AF, while M2 macrophage-related genes IGF1, RAB13, TMEM176B, and PDK4 may regulate the progression of AF through the PPAR signaling pathway and fatty acid metabolism. Our study reveals the characteristics of immune cell infiltration and the expression pattern of M2-like macrophage-related genes in LAA tissues of AF patients at the cellular and molecular levels, providing new theoretical insights for the early diagnosis and treatment of AF.