Gene expression datasets (GSE14975, GSE41177, and GSE79768) for AF were retrieved from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The transcriptional profiling in GSE14975 (Platform: GPL570) comprises left atrial myocardium from 5 AF patients and 5 matched samples of patients in sinus rhythm (controls) (Adam et al., 2010). GSE41177 (Platform: GPL570) consists of 19 left atrial appendage samples from 16 AF patients and 3 controls (Yeh et al., 2013). These two datasets were merged (the clinical information of these two datasets was displayed in Supplementary Table 1), yielding a training set of 21 AF cases and 8 controls, after which the “sva” package (Leek et al., 2012) was employed to eliminate batch effects, hence ensuring data uniformity. GSE79768 (Platform: GPL570) was the validation set in this study, which included 13 left atrial samples from 7 AF patients and 6 controls (Tsai et al., 2016). Besides, a total of 68 mitochondria-associated ER membrane-related genes (MAMRGs) were extracted from the early research (Li H. M. et al., 2024).

To identify AF-related genes, the “limma” R package (version 3.54.0) was employed to discover the differentially expressed genes (DEGs) between AF and controls in the training set (Ritchie et al., 2015). The |logFoldChange| > 1 and adjusted P-value <0.05 were screening thresholds. Moreover, results were displayed as a volcano plot [using the “ggplot” R package (version 3.4.1) (Wickham, 2016)] and a heatmap [using the “ComplexHeatmap” R package (version 2.14.0) (Gu et al., 2016)]. Subsequently, AF-MAMRGs were determined through utilizing a Venn diagram to link the AF-related genes and MAMRGs. The STRING database (https://cn.string-db.org/) was applied to construct the protein-protein interaction (PPI) network to explore the direct and indirect associations among the AF-MAMRGs (parameter setting: Low confidence = 0.4).

To comprehend the potential biological functions of AF-MRMRGs, the functional enrichment analysis was conducted with the “clusterProfiler” R package (version 4.2.2), including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses (P < 0.05) (Wu et al., 2021). GO analysis comprised three components: Biological Process (BP), Cellular Component (CC), and Molecular Function (MF).

The support vector machine-recursive feature elimination (SVM-RFE) is a technique for feature selection that recursively eliminates the least important features using SVM. To evaluate the significance of AF-MRMRGs in the diagnosis of AF, SVM-RFE was first utilized to screen characteristic genes by the “e1071” R package. Simultaneously, a least absolute shrinkage and selection operator (LASSO) regression model with 10-fold cross-validation was constructed using the glmnet R package (v4.1.7) to enhance confidence in the identified genes (Huang et al., 2024). After that, the expression of characteristic genes between AF and controls was analyzed in both training and validation sets, which with the significant alterations and consistent expression trends in two sets were denoted as biomarkers for AF.

A nomogram is a graphical tool widely used to predict the probability of a particular outcome based on a set of variables. Therefore, based on the expression levels of the identified biomarkers, a nomogram was developed to predict the probability of AF. Moreover, the calibration curve, receiver operating characteristic (ROC) curve, and decision curve analysis (DCA) curve were plotted to assess the predictive performance of the nomogram.

The GSEA (Guo et al., 2024a; Xu et al., 2025a; Yang et al., 2025) and GSVA analyses were conducted to investigate the potential mechanisms of biomarkers in the occurrence of AF. For GSEA, the “c2. cp.kegg.v7.4. symbols.gmt” was extracted from the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb) as the reference set. The correlations between biomarkers and all genes were calculated, and correlation coefficients were ranked. GSEA was carried out by the “clusterProfiler” R package (version 4.2.2) with a P-value <0.05. Additionally, the samples in the training set were categorized into high- and low-expression groups according to the expression of each biomarker. With “h.all.v7.4. symbols.gmt” downloaded from the MSigDB database as the reference set, the GSVA score of each pathway in two expression groups was calculated by the “GSVA” R package (Li J. et al., 2024), and the differences in GSVA scores between the two groups were assessed via the “limma” R package with |t| > 2 and P-value <0.05.

Previous studies have demonstrated that the immune system undergoes significant changes during AF and interacts with the environment and cells involved in the initiation and maintenance of AF (Yao et al., 2022). Consequently, the infiltration abundance of 28 immune cells in AF and controls was evaluated by single-sample gene set enrichment analysis (ssGSEA) (Cui et al., 2025), and the differences between AF and controls were analyzed via the Wilcoxon test (P < 0.05). To further explore the correlations between biomarkers and immune infiltration, the Spearman algorithm was applied.

The DGIdb database (https://www.dgidb.org/) (Cannon et al., 2024) is an online database of drug-gene interactions, with data sourced from multiple drug databases (DrugBank, PharmGKB, ChEMBL), clinical trial databases, and PubMed literature. The potential drugs were predicted for AF using biomarkers as the keywords. To further investigate the specific mechanism, molecular docking was conducted by AutoDock Vina (Eberhardt et al., 2021). The 3D structures of drugs were downloaded from the PubChem database (SDF file) (Kim et al., 2025), and the SDF files were transferred into PDB files via Babel GUI(O'Boyle et al., 2011). Moreover, the 3D structures of proteins were extracted from the Protein Data Bank Database (Burley et al., 2025). Finally, the PyMol software (Mooers and Brown, 2021) was employed to view and visualize the results.

The scRNA-seq dataset of AF (GSE224959, GPL18573) was also retrieved from the GEO database, collecting left atrial tissue from 7 AF samples and 5 controls (Hulsmans et al., 2023). The processing of data was performed via the “Seurat” R package (Satija et al., 2015) according to the following criteria: (1) cells were removed with less than 200 genes expressed; (2) genes were excluded with less than 3 cells covered; (3) cells were excluded with more than 20% mitochondrial genes; (4) cells were removed with 200 ≤ gene numbers ≤3,000; and (5) genes were excluded with 200 ≤ count numbers ≤10,000. After identification of the top 2,000 highly variable genes and PCA dimensionality reduction analysis, the cluster analysis was conducted for remaining cells using the FindNeighbors and FindCluster functions (Di et al., 2022) in the “Seurat” R package. Subsequently, cell clusters were annotated through comparing the differentially expressed genes of each cluster with the marker genes of each cell type in the CellMarker database (Zhang et al., 2019). Afterwards, the percentage of each cell type was counted, and the expression of each biomarker in each cell type between AF cases and controls was measured. Moreover, the pseudotime analysis of the cell type where the expression of biomarkers between AF and controls was at significant levels was conducted to investigate the change of biomarker expression during cell differentiation. Finally, to reveal interactions among these cell types, CellPhoneDB (Efremova et al., 2020) was employed.

To further investigate the causal effect of biomarkers on AF, MR analysis was carried out using TwoSampleMR R package (Hemani et al., 2018). The data of the exposure and outcome were sourced from the IEU OpenGWAS database, including eqtl-a-ENSG00000089022 (MAPKAPK5) (31,470 samples and 15,599 single nucleotide polymorphisms (SNPs)), eqtl-a-ENSG00000141510 (TP53) (31,684 samples and 18,457 SNPs), eqtl-a-ENSG00000204632 (HLA-G) (25,690 samples and 28,777 SNPs), and bbj-a-71 (AF) (36,792 samples (8,180 AF and 28,612 control) and 5,018,048 SNPs). To satisfy the three assumptions of MR analysis, SNPs employed as IVs must meet the following criteria: (1) IVs must be significantly relevant to the exposure with P < 5 × 10⁻⁶; (2) the SNPs with linkage disequilibrium (LD) were removed (parameter settings: r ² = 0.01, kb = 100); (3) the F statistic value was calculated to evaluate the strength of SNPs, and those with F less than 10 were excluded. Afterwards, the two-sample MR analysis was conducted through combining the harmonise_data and mr functions with five methods (Weighted median, MR Egger, Simple mode, Inverse variance weighted (IVW), and Weighted mode), and among these methods, the IVW was the main method. Causal effects were represented by odds ratios (OR) and their 95% confidence intervals (CIs), and statistical significance was defined as P less than 0.05. Furthermore, to measure the robustness of results, the sensitivity analysis was conducted through the heterogeneity test, horizontal pleiotropy test, and leave-one-out (LOO) test. The choice between IVW-fixed effects (IVW-FE) and IVW-random effects (IVW-RE) was determined by the presence of heterogeneity. Specifically, IVW-RE was selected when heterogeneity existed, whereas IVW-FE was adopted in the absence of such heterogeneity. The mr_heterogeneity () function was employed to detect heterogeneity. Horizontal pleiotropy was identified using the MR-Egger intercept method: a significant difference between the intercept term of MR-Egger regression and 0 (P < 0.05) indicated the existence of horizontal pleiotropy. Additionally, the LOO method was used to evaluate the impact on the MR results after excluding each individual SNP.

The AF model was established by rapid atrial pacing using canines. The sham group received pacemaker implantation procedure under sterile conditions without atrial pacing. The pacing group underwent pacemaker implantation with continuous rapid atrial pacing (450 beats/min) for 2 months before the electrophysiological measurements. All electrophysiological measurements were recorded in a computerized electrophysiology system (Lead 7,000, Jinjiang Inc., China). An S1S1 programmed stimulus method (with cycle length of 120 ms, 100 ms, and 75 ms; 5 s each, performed 3 times per frequency) was used to evaluate AF durations. AF was characterized as an irregular atrial rate exceeding 500 bpm and that lasts for more than 5 s.

All trial protocols were approved by the Laboratory Animal Welfare & Ethics Committee of Dongguan Songshan Lake Central Hospital. The experiments involving live vertebrates were carried out in strict accordance with the relevant guidelines and regulations. All methods were reported in accordance with ARRIVE guidelines.

Total RNA was isolated from canine atrial tissue samples utilizing RNA extraction solution (G3013, Saiwei Biotechnology, CN) in accordance with the procedure. Tissue homogenization was conducted using a three-dimensional cryogenic grinder (KZ-5F-3D, Saiwei Biotechnology, CN), followed by purification processes involving chloroform extraction, isopropanol precipitation, and washing with 75% ethanol solution. The purity and concentration of RNA were evaluated using a Nanodrop 2000, and the concentration was standardized to 200 ng/μL cDNA synthesis was performed utilizing SweScript All-in-One RT SuperMix for qPCR (One-Step gDNA Remover, G3337, Saiwei Biotechnology, CN) to reverse transcribe RNA into cDNA, with the reaction protocol established at 25 °C for 5 min, 42 °C for 30 min, and 85 °C for 5 s. The 2×Universal Blue SYBR Green qPCR Master Mix (G3326, Saiwei Biotechnology, CN) was utilized for qPCR, with each 15 μL reaction comprising the master mix, primer mixs for HLA-G, TP53, and MAPKAPK5, cDNA template, and nuclease-free water. Primer sequences are presented in Table 1. Each sample was run in triplicate. The amplification was performed on a CFX Connect Real-Time PCR System (CFX Connect, Bio-Rad, CN) with the thermal cycling conditions: initial denaturation at 95 °C for 30 s, followed by 40 cycles of 95 °C for 15 s and 60 °C for 30 s, and a melting curve stage from 65 °C to 95 °C. Relative gene expression levels were calculated using the 2^(−ΔΔCT) method.

The expression levels of TP53, MAPKAPK5, and HLA-G in atria were detected by Western blotting. After atrial tissues were washed with PBS, lysis buffer prepared with proteinase inhibitors and RIPA strong lysis buffer was added. Tissues were homogenized on ice for 30 min using a liquid nitrogen grinder. Samples were centrifuged at 4 °C, 12,000 rpm for 20 min using a refrigerated high-speed centrifuge. After centrifugation, collect the supernatant and measure the protein concentration using the BSA kit. An appropriate amount of 5× loading buffer was added to the remaining supernatant, and heat at 95 °C in a metal bath for 10 min. Electrophoresis of protein samples was performed using sodium dodecyl sulfate polyacrylamide gel electrophoresis (SDS-PAGE). After electrophoresis, proteins were transferred onto a PVDF membrane and blocked with 5% non-fat dry milk for 2 h. After blocking, the membrane was washed 5 times with TBST solution, each for 6 min. The membrane was then incubated with primary antibodies: anti-TP53 (IPD-ANP1072, IPODIX, China, dilution 1:500), anti-MAPKAPK5 (IPD-ANP9498, IPODIX, China, dilution 1:500), anti-HLA-G (GB115645, Servicebio, China, dilution 1:500), and anti-GAPDH (GB15004, Servicebio, China, dilution 1:10,000) at 4 °C for at least 12 h. After the blocking period, wash the membrane 5 times with TBST, each for 6 min. Next, add the corresponding secondary antibody working solution and incubate for 2 h at room temperature. After incubation, wash the membrane 3 times with TBST, each for 6 min. To ensure that the loaded samples were of equal concentration, the ratio of band intensity to GAPDH was calculated to quantify the relative expression levels of these proteins.

All statistical analyses were conducted employing R software (version 4.2.2). The disparities were analyzed via the Wilcoxon test (n = 2). The data obeyed normal distribution characteristics and are shown as the mean ± standard deviation (SD). Two-sample independent Student’s t-test was used to compare the means of two groups. A statistically significant P-value is less than 0.05.

After removing the batch effects, PCA demonstrated that AF samples could be clearly separated from control samples (Figure 1A). Through differentially expressed analysis, a total of 4,384 DEGs between AF and controls were obtained, including 4,380 upregulated and 4 downregulated genes (Figures 1B,C). Subsequently, through combination of MAMRGs and DEGs, there were 18 AF-MAMRGs obtained (Figure 1D). The PPI network of AF-MAMRGs included 16 nodes and 23 edges (Figure 1E), with TP53 had interactions with most proteins, such as MIF, MAPKAPK5, and SPI1. Additionally, the enrichment analysis illustrated that a total of 958 GO items and 75 KEGG pathways were markedly enriched (Supplementary Table 2). The top 10 GO items and the top 5 KEGG pathways were displayed, including “cellular senescence (GO-BP item),” “integral component of lumenal side of endoplasmic reticulum membrane (GO-CC item),” “DNA-binding transcription factor binding (GO-MF item),” and “Human T-cell leukemia virus 1 infection (KEGG)” (Figures 1F,G).

To discover the genes with diagnostic value for AF among AF-MAMRGs, we first identified 15 characteristic genes through SVM-RFE and these 15 genes were incorporated in the LASSO results (17 genes) to ensure robustness (Figures 2A,B). For the expression examination, in both training and validation sets, TP53, MAPKAPK5, and HLA-G were significantly upregulated in AF compared to controls (Figures 2C,D).

Consequently, these three genes were identified as biomarkers for AF. To further analyze the causal effect of biomarkers on AF, a two-sample MR analysis was conducted. Since the number of SNPs of TP53 and HLA-G were less than 3 after IV selection, only the causal effect of MAPKAPK5 on AF was measured. Results illustrated that MAPKAPK5 had a significant causal association with AF (P = 0.022) and was a risk factor for AF onset (OR = 1.065, 95%CI = 1.009–1.125) (Table 2). Moreover, the scatter plot, forest plot, and funnel plot also demonstrated the consistency and reliability of the findings (Figures 2E–G). The Q statistic P-values for MR in this analysis were all below 0.05, which prompted the selection of a fixed-effect model. No heterogeneity and horizontal pleiotropy was detected (Q value = 0.415; Egger_intercept = 0.074; P = 0.143), and removal of any individual SNP had no effect on results, suggesting that the results were reliable and robust (Table 3; Figure 2H).

We enhanced the nomogram in both training and validation cohorts to forecast the occurrence of AF based on the expression of three biomarkers (Figures 3A,E). Both calibration curves indicated that the slopes were approximately 1 (Figures 3B,F). The area under the ROC curves were 0.821 and 1, respectively (Figures 3C,G). Furthermore, the DCA indicated that the nomogram exhibited the highest net benefit (Figures 3D,H). These results indicated that the nomogram had high accuracy and robustness for predicting AF.

Next, we analyzed the specific signaling pathways involved in the three biomarkers and explored the impact of the biomarkers on the signaling pathways associated with AF progression. GSEA results demonstrated that three biomarkers were involved in immune-related pathways, such as “cytokine-cytokine receptor interaction,” “chemokine signaling pathway,” “TOLL like receptor signaling pathway,” “MAPK signaling pathway,” “natural killer cell mediated cytotoxicity” (Figures 4A–C) (Supplementary Table 3). Moreover, these biomarkers were also involved in metabolism-related pathways, including “glycerophospholipid metabolism” and “propanoate metabolism” (Supplementary Table 3). Results of GSVA illustrated that the HALLMARK pathways enriched of the two expression groups of TP53 and HLA-G were identical. The high-expression group was enriched in “mitotic spindle,” “oxidative phosphorylation,” “UV response DN” and other 13 pathways (Figure 4D; Supplementary Table 3). The low-expression group was enriched in “allograft rejection,” “inflammatory response,” “IL6-JAK-STAT3 signaling” and other 6 pathways (Figure 4D; Supplementary Table 3).

Previous studies have demonstrated that the immune response may play an important role in the development and maintenance of AF. Moreover, results of GSEA revealed that three biomarkers were relevant to immune-related pathways in our study. Therefore, we further measured the infiltration abundance of 28 immune cells in AF and controls, and group differences were analyzed (Figures 4E,F). Between AF and control groups, the infiltration abundance of 18 immune cells was considerably increased in the AF group compared to the control group, including macrophages, activated dendritic cells, neutrophils, monocytes, etc.(Figure 4F). Both activated CD4⁺/CD8⁺ T cells and central memory CD4⁺/CD8⁺ T cells exhibited significant differences in infiltration levels between groups, highlighting the involvement of diverse T cell subsets. Moreover, all biomarkers had marked positive relevance to these immune cells, and the highest correlation was between TP53 and monocytes (cor = 0.91, P = 7.495e-12), meanwhile, myeloid-derived suppressor cells (MDSC) was highly correlated with all biomarkers (cor >0.95) (Figure 4G).

To find drugs for AF therapy, we searched the DGIdb database. Results demonstrated that two potential drugs interacted with HLA-G, two potential drugs interacted with MAPKAPK5, and 49 potential drugs interacted with TP53 (Figure 5A). Moreover, according to the ranking of interaction scores, the drug with the highest interaction score for each biomarker was selected for molecular docking. Results demonstrated that the docking energy between MAPKAPK5 and GLPG-0259 was −9.5 kcal/mol, and SIMVASTATIN had a strong binding affinity with residues LYS-6 and GLU-232 of HLA-G through hydrogen bonding, with a docking energy of −5.74 kcal/mol (Figures 5B,C). These results indicate a stable binding capacity between MAPKAPK5 and HLA-G and their respective targeted drugs, suggesting their potential as therapeutic candidates for further investigation in related pathways. However, the docking energy between TP53 and THIOUREIDOBUTYRONITRILE was −3.84 kcal/mol, indicating that the bond between them was not stable.

Following quality control of scRNA-seq data, the top 2,000 highly variable genes were determined, and the top 30 principal components (PCs) were selected (Supplementary Figure S1). Through cluster analysis, a total of 16 cell clusters were determined, which were annotated as 6 cell types: mononuclear phagocytes and dendritic cells (MP/DCs), lymphoid cells, mural cells, neutrophils, endothelial cells, and fibroblasts (Figures 6A–C). Among these cell types, MP/DCs accounted for the highest proportion in both AF and controls (Figure 6D). Expression analysis illustrated that between AF and controls, the expression of TP53 in lymphoid cells, neutrophils, and endothelial cells and the expression of MAPKAPK5 in lymphoid cells and neutrophils were significantly different (Figure 6E). However, between AF and control groups, the expression of HLA-G in each cell type did not differ significantly (Figure 6E). Afterwards, the pseudotime analysis was employed on lymphoid cells and neutrophils. Results demonstrated that TP53 and MAPKAPK5 were expressed throughout cell differentiation (Figures 7A,B). Additionally, cell-cell interaction analysis illustrated communication patterns among 6 cell types in the AF microenvironment. Compared to the control group, the bidirectional interaction between neutrophils and endothelial cells in AF was increased (Figure 7C), which thus drew our focus. Further analysis of ligand-receptor pairs revealed an increased signaling intensity of the APP-CD74 axis, a key communication pathway between neutrophils and endothelial cells. This suggests that the pathway may exacerbate immune activation and local atrial inflammatory responses by enhancing the “recruitment signals” from endothelial cells to neutrophils. In contrast, the ITGB2-ICAM2 axis between neutrophils and endothelial cells was diminished in the AF group. Unlike controls, the interactions between lymphoid cells and mural cells as well as between lymphoid cells and neutrophils in AF were decreased (Figure 7C).

AF model was established by rapid atrial pacing in canines to further investigate the expression of TP53, MAPKAPK5 and HLA-G in atria. As shown in Figures 8A,B, electrophysiological measurements demonstrated that both the induction times and duration of AF were significantly higher in the AF canine model. Consistent with these findings, electrocardiographic (ECG) recordings revealed characteristic changes in the AF model (Figure 8C). Subsequently, Western blotting (Figures 8D–G) and qRT-PCR (Figures 8H–J) were employed to ascertain the expression levels of TP53, MAPKAPK5, and HLA-G in both the AF group and the sham group. The findings demonstrated that TP53, MAPKAPK5, and HLA-G were markedly overexpressed in the canine model of AF. The alignment between transcriptional and protein levels enhances our assurance in the differential expression of these genes in the etiology of AF.