The study flowchart was shown in Figure 1. The datasets (GSE41177 and GSE79768) were obtained from the Gene Expression Omnibus database (GEO, http://www.ncbi.nlm.nih.gov/geo/), There are 14 AF samples and 12 sinus rhythm samples from 7 patients with AF and 6 patients with sinus rhythm in the dataset GSE79768 (GPL570 HG-U133_Plus_2 Affymetrix Human Genome U133 Plus 2.0 Array) (10), and 32 AF samples and 6 sinus rhythm samples from 16 patients with AF and 3 patients with sinus rhythm in the dataset GSE41177 (GPL570 HG-U133_Plus_2 Affymetrix Human Genome U133 Plus 2.0 Array) (11). Therefore, 46 AF samples and 18 sinus rhythm samples from 23 patients with AF and 9 patients with sinus rhythm were included in our analysis. Of note, left atrial tissue sample and right atrial tissue sample were respectively obtained from every patient in the dataset GSE79768, and tissue sample in left atrial appendage and tissue sample in left atrium-pulmonary vein junction were respectively obtained from every patient in the dataset GSE41177.

Subsequently, the normalization of the datasets was performed using the “SVA” package. We first read the genes information of two datasets to obtain each gene symbol. When different expression values were presented in a same gene, the mean values were calculated. Logarithmic transformation was performed for the non-logarithmic value of genes expression. Then, two datasets were combined and internally normalized. Finally, the normalization of the dataset was output.

The protocols of 2 studies analyzed in our study were conducted according to the Helsinki Declaration, and obtained the ethical approval. All patients enrolled in the 2 studies were given written informed consent (10, 11).

The genes were analyzed by using “limma” package in R software, thresholded at |logFC| > 0.5 and adjusted p value <0.05. When different expression values were presented in a same gene, the mean values were calculated. When the genes expression values were 0, the genes were excluded. Then, the differential expression genes (DEGs) between samples with AF and sinus rhythm were obtained, when adjusted p values <0.05 correcting for the false discovery rate method (12). Finally, the DEGs were output and visualized by heatmap.

We then performed GO, KEGG, and GSEA enrichment analyses of DEGs by using “clusterProfiler and enrichplot” packages, and the immunological signature genome was used as a reference for GSEA. A p value <0.05 was considered statistically significant.

We constructed a weighted co-expression network using the WGCNA package to normalize the expression profiles of GSE41177 and GSE79768 batches. We checked the data by using “goodSampleGenes” and used the “pickSoftThreshold” function to obtain the appropriate soft threshold (β), then created a scale-free topology network while converting the matrix data into an adjacency matrix. After calculating module eigengene (ME) and merging similar modules in the clustering tree according to ME, we obtained a hierarchical clustering dendrogram. Gene significance (GS) and module significance (MS), which were used to measure genetic and clinical information, were calculated, and the correlations between modules and models were analyzed. Finally, we calculated the module membership (MM) for each gene to analyze the GS in the modules.

We chose the LASSO regression algorithm to identify important variables as well as to improve the accuracy of prediction. The central genes with the lowest p value, high inter-module connectivity with |GS| > 0.50, and |MM| > 0.80, were firstly selected for screening of hub genes. Then, “Venn” package was used to intersect the central genes with DEGs. Finally, hub genes were obtained via “glmnet” analysis.

We used the “ggpubr” package to construct box plots to assess hub genes expression levels between patients with AF and sinus rhythm, and plotted receiver operating characteristic curves (ROC) were used to analyze the accuracy of the hub gene for diagnosis by the “pROC” package.

We constructed an animal AF model using the rats with transverse aortic constriction (TAC) (13–15). Meanwhile, the rats without TAC were used as control group. Eight rats were included in AF group, and 7 rats were included in control group. The temperature, appetite, and weight were detected in rats in each group in everyday, to exclude infectious diseases and heart failure. One month after TAC, the rats in each group were anesthetized and fixed supinely on an animal test bench. Then, vulnerability to AF of the rats in each group were measured by 3 cyclically transesophageal burst atrial pacing at a frequency of 1,500 stimuli/min for 30 s, with 5 min interval between cycles.

The experimental protocols were performed in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals, and were approved by the Ethics Committee of Guangxi Medical University Laboratory Animal Center (Nọ: 202207002).

The atrial tissue of rats in each group were collected after measurement of vulnerability to AF, and were snap frozen in liquid nitrogen for RNA isolation. Total RNA was extracted from the tissue using RNA extraction reagent (Servicebio, Wuhan, China). The concentration and purity of the extracted RNA were detected using an ultra-microspectrophotometer (Thermo Fisher Scientific, Massachusetts, USA). 1 μg of total RNA was reverse transcribed into cDNA using the Prime Script™ RT reagent Kit with gDNA Eraser (Takara, Tokoyo, Japan). The quantitative polymerase chain reaction (qRCR) procedure was carried out as follows: 95°C for 30 s, 95°C for 5 s, and 60°C for 30 s for 40 cycles. Melt-curve analysis was performed at 65–95°C. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as the internal control for normalizing gene expression. The data obtained were calculated by the 2^–ΔΔCt method (16). The sequences of primers were shown in Table 1.

The ssGSEA algorithm was used to analyze the proportion of 28 types of immune cells in the immune infiltrative microenvironment of AF, the “vioplot” package was used to assess the immune cells infiltration between AF and control group, and Spearman correlation analysis was used to calculate the correlation between Hub genes expression and the proportion of immune cells infiltration.

qPCR data were analyzed using SPSS 23.0 software (SPSS Inc., Chicago, IL, USA). Two unpaired Student's t-test was performed to compare the difference between AF group and control group. All the results of qPCR in the study were visualized by graph pad Prism 8 (GraphPad Software, Inc., San Diego CA, USA). A p value <0.05 was considered statistically significant. All other statistical analyses were performed using R software version 4.1.0 (www.r-project.org/).

After the datasets were merged and normalized, a total of 298 DEGs were obtained. The results were shown in Figure 2.

To understand the biological functions and signaling pathways associated with AF, we performed the differential enrichment analysis. GO enrichment analysis revealed that DEGs were associated with neutrophil degranulation, neutrophil activation involved in immune response, immune receptor activity, RAGE receptor binding, and collagen-containing extracellular matrix (Figure 3A). KEGG enrichment analysis showed that DEGs were associated with cytokine-cytokine receptor interactions, chemokine signaling pathways and B cell receptor signaling pathways, etc. (Figure 3B). The MsigDB database of the immune signature gene set was referenced for the GSEA to explore the potential mechanisms of immune function in AF. GSEA revealed that DEGs were significantly enriched in B cells, CD⁴⁺ T cells, CD⁸⁺ T cells (Figure 3C).

We used the WGCNA package to construct a gene co-expression network. After the samples were clustered to handle missing values and remove outliers, a scale-independent topological network was built (soft threshold β = 7, scale-free R² = 0.83; slope = −1.59). A co-expression matrix was also constructed using a one-step method. Mixed cuts were made to construct a hierarchical clustering tree, and 10 gene modules were subsequently obtained by fusion of similar modules.

The correlation of the above modules with AF or sinus rhythm was clarified using heatmap. The results indicated that the module (including CLEC4A, COTL1, EVI2B, FCER1G, GAPT, HCST, NCF2, PILRA, TLR8, and TYROBP) had the highest correlation with AF (ρ = 0.48, p = 5 × 10⁻⁵) (Figure 4). Furthermore, there was a good correlation between GS and MM in this module (ρ = 0.42, p = 3.4 × 10⁻²⁴). Therefore, we finally chose this module as the candidate module.

Ten genes from the black module were obtained using |GS| > 0.50 and |MM| > 0.80 as a cutoff value, and analyzed using the LASSO algorithm. The 4 hub genes including PILRA, NCF2, EVI2B, and GAPT were obtained finally by LASSO regression analysis after intersecting with DEGs (Figure 5).

We used box plots to verify the expression levels of the 4 Hub genes, the results presented in Figure 6 suggested that the expression of PILRA, NCF2, EVI2B, and GAPT were higher in atrial tissue of patients with AF, compared to those in atrial tissue of patients with SR. To verify whether the 4 hub genes have a good diagnostic value, we computed the values of the area under the curve (AUC) of the Hub genes. The AUC values of PILRA, NCF2, EVI2B, and GAPT were 0.856, 0.890, 0.861, and 0.862, respectively (Figure 7). These results indicated that these genes had a good diagnostic value for AF.

The expression levels of 4 hub genes in rats with AF and sinus rhythm were measured by qPCR. There were 7 rats in Control group and 8 rats in AF group. Two rats in AF group with TAC were dead of heart failure. Finally, 7 rats in control group and 6 rats in AF group were performed transesophageal burst atrial pacing to measure the vulnerability of AF. The electrocardiograms of rats in control and AF groups were shown in Figure 8. The results have been shown in Table 2. One rat in Control group was dead after transesophageal burst pacing. The expression levels of PILRA were measured in 6 rats in each group, and expression of NCF2, EVI2B, and GAPT were measured in 5 rats in each group. The results shown in Figure 9 revealed that the expression level of PILRA was significantly elevated in rats with AF, compared to those in rats with SR, but the expression levels of NCF2, EVI2B, and GAPT were no obvious difference between the rats with AF and the rats with SR. These results further proved that PILRA may be a potential target for intervention of AF.

We used the ssGSEA to explore the different types of immune cells between AF and SR, and evaluate their relationship. The distribution of immune cells in AF and SR, shown in Figures 10A,B revealed that immunity had a correlation with AF. Activated CD⁴ T cell (p = 0.019), Activated CD⁸ T cell (p = 0.024), Activated dendritic cell (p = 0.002), Gamm delta T cell (p = 0.016), Immature B cell (p < 0.001), myeloid-derived suppressor cell (MDSC) (p < 0.001), Macrophage (p = 0.001), Mast cell (p = 0.001), Monocyte (p = 0.002), Neutrophil (p < 0.001), Plasmacytoid dendritic cell (p = 0.003), Regulatory T cell (p < 0.001), Type 17 T helper cell (p = 0.001), Effector memory CD⁴ T cell (p = 0.005), Central memory CD⁴ T cell (p = 0.008), Central memory CD⁸ T cell (p = 0.001) and Effector Memory CD⁸ T cell (p < 0.001) were positively correlated with AF (Figure 10C).

The Spearman correlation analysis showed that PILRA was significantly and positively correlated with Immature B cell (p < 0.01), dendritic cell (p < 0.001), Myeloid–derived suppressor cell (p < 0.001), mast cell (p < 0.01), macrophage (p < 0.001), and T cells and their subpopulations (including Regulatory T cell, Effector memory CD⁸ T cell, Effector memory CD⁴ T cell, Central memory CD⁴ T cell, Activated CD⁸ T cell, Activated CD⁴ T cell (p < 0.001), Gamma delta T cell (p < 0.01), Central memory CD⁸ T cell (p < 0.05). These results suggested that PILRA may play a crucial role in regulation of immune cells infiltration that enhanced inflammatory response to promote AF, and PILRA may be a novel target for intervention of AF.