The AF dataset, GSE2240, was obtained from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/). GSE2240 contained data on 30 right atrial appendages from 10 patients with permanent AF, defined as AF duration longer than 3 months, and 20 patients with sinus rhythm (SR), with no history of AF.

Raw data were pre-processed identically R for background correction and normalization. R package annotation was conducted to match probes and gene symbols, and the probes matching several genes were removed. The median was regarded as the final expression value for the gene matched by multiple probes. We calculated the standard deviation (SD) for each gene and ranked them from the largest to the smallest, and the top 5000 genes were chosen for WGCNA refer to the official tutorials.

A gene co-expression network was constructed using the R package WGCNA (10, 17). In this analysis, nodes represented genes, and edges represented the degree of co-expression. Firstly, the hclust function was used to cluster the samples to determine whether there were outliers samples. Then, weighted correlation coefficients were introduced to calculate the adjacency matrix of the expression profile genes. The co-expression similarity between genes i and j was defined as S_ij = |cor(i,j)|. The correlated adjacency of the genes was further analyzed using the power function: a_ij = |S_ij|^β. After that, soft threshold β was chosen within a specific range to satisfy the scale-free network and better network connectivity. Fit index R² > 0.85 was set to make the connections between genes obey the approximate scale-free network distribution. And the pickSoftThreshold function was applied to automatically filter the appropriate soft threshold β. Finally, the blockwiseModules function was used for network construction and module detection to generate the topological overlap matrix (TOM) with a minimum module size of 50 and a merge cut height of 0.25 for co-expressed gene modules. Based on the heterogeneity of TOM, the genes with similar expression patterns were divided into the same module by means of average link hierarchical clustering, and the gene modules were identified by dynamic shearing tree method.

Gene modules are clusters of genes that are closely related to co-expression. WGCNA uses a hierarchical clustering approach to identify gene modules and represents each cluster with a different color. Genes not assigned to any module are placed in gray modules. Principal component analysis was performed for each module, and the module eigengenes (MEs) were calculated using the first principal component representing the overall expression level of the module. Module-trait correlation coefficients were calculated to provide a heat map of the correlation coefficients between modules and traits (18). Then, we screened the gene modules significantly associated with the traits using the eigenvectors of modules, correlation coefficients of traits, and module significance (P < 0.05). The values of gene significance (GS) and module membership (MM) were calculated. GS represents the correlation between gene expression and traits within the calculated module, and MM represents the correlation coefficient between the expression of a gene and the expression of the main component of the gene within the module. Finally, the list of genes calculated by the networkScreening function was screened by setting the value range of GS, MM, and q.weighted to identify and characterize the key hub genes.

Gene modules significantly associated with traits were selected, and network maps were drawn by weighted co-expression relationships between genes using the Cytoscape software. Functional annotation analysis of genes in the core module based on the Gene Ontology (GO) database and signaling pathway enrichment analysis of genes in the core module based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database were performed to determine the biological functions and potential biological pathways of genes in the trait-related modules.

All protocols and the use of human blood samples were in accordance with the Declaration of Helsinki and approved by the Human Ethics Review Committee of the Chinese PLA General Hospital. All leukocyte samples were obtained from the PLA General Hospital Biosample Bank, including 9 elderly AF patients and 9 elderly non-AF controls. Blood samples were collected in EDTA anticoagulation tubes. Leukocytes were isolated from human peripheral blood by Ficoll density gradient centrifugation according to the manufacturer's instructions and stored at −80°C until further analysis. Detailed patient characteristics are summarized in Supplementary Table S1.

Total RNA was extracted from leukocytes using an RNAprep Pure Blood Kit (TIANGEN Biotech Corporation, China), according to the standard protocol. Absorption spectrophotometry using NanoDrop-1,000 (Thermo Fisher Scientific, Yokohama, Japan) was used to determine RNA concentrations and purity (260/280 ratio >1.8). The RNA (1 μg per sample) was reverse transcribed using the High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems, CA, USA) in a total reaction volume of 20 μl, following the manufacturer's instructions. Then, RT-qPCR was performed as follows: pre-denatured at 95°C for 30 s; denatured at 95°C for 5 s, and 40 cycles at 60°C for 30–34 s. All experiments were performed in triplicates. The housekeeping gene GAPDH was used as the endogenous reference. The primer sequences of PTGDS, COLQ, ASTN2, VASH1, RCAN1, AMIGO2, RBP1, MFAP4, and ALDH1A1 used in this study are listed in Supplementary Table S2. The relative mRNA expression levels were calculated using the 2^−ΔΔCt method. A flow chart of this study is shown in Figure 1.

Statistical analyses were performed using SPSS statistical package (version 21.0; SPSS Inc., Chicago, USA). Continuous variables were presented as mean ± standard error (SE), and categorical variables were expressed as percentages. Student's t-test or the chi-square test was used to determine the differences between AF patients and controls. The sample size was analyzed by PASS (version 15.0; NCSS, USA). Differences were considered statistically significant at P < 0.05.

Gene expression values of all samples were subjected to sample clustering and phenotypic heat map analysis. As shown in Figure 2A, there were no significant outlier samples. Therefore, all the samples were included in the subsequent data analysis. Using the WGCNA, we calculated and selected β = 8 as the soft threshold for the dataset based on the scale-free network fit index and average connectivity (Figures 2B,C). The adjacency and TOM matrices between genes were calculated. A hierarchical clustering tree of the genes was constructed based on the TOM matrix. Then, the genes were divided into 19 modules using the dynamic shearing tree method. Each module was represented by a colored rectangle, and the vertical coordinates were the gene occupancy ratio. The tree branches correspond to 10 different gene modules, each leaf on the tree corresponds to a gene, and similar genes are clustered into modules of the same color (Figures 3A,B). The module name and the number of genes included in each module were as follows: turquoise (1,617), blue (513), brown (489), yellow (430), green-yellow (105), black (188), magenta (163), pink (177), red (199), purple (135), yellow (332), grey (21), tan (95), salmon (92), cyan (91), midnight blue (79), light cyan (68), grey60 (68), and light green (58). Genes from each module are listed in Supplementary Table S3. The genes in the gray module could not be clustered to any other module and were removed from the subsequent analysis.

The correlation heat map between co-expression modules and clinical traits was drawn by calculating the relationship between each gene module and clinical traits. Hierarchical clustering and heat map analyses were performed for each module. The correlations between the modules are shown in Figure 4A. The pink and green modules are distributed in different clustering subtrees. In this study, the pink module exhibited the highest positive correlation (r = 0.76, P < 0.0001) with the SR phenotype, and the green module exhibited the highest negative correlation with the AF phenotype (r = −0.8, P < 0.0001) (Figure 4B). Therefore, we focused on these two modules in the subsequent analysis.

The pink and green modules were used as key modules for GS and MM analyses, respectively. The correlation coefficient between GS and MM was r = 0.75, P < 0.0001 (Figure 5A) and r = 0.65, P < 0.0001 (Figure 5B) for the green and pink modules, respectively. Figures 5C,D represent the heat map of eigengene expression for the two modules. Figures 5E,F show the heat and clustering map between the genes of the two modules and each clinical trait.

Finally, three criteria were used to screen for key hub genes in the pink and green modules: GS > 0.65, MM > 0.8, and q.weighted <0.01. After excluding non-coding genes, 16 hub genes were screened in the green module (Table 1) and 23 hub genes in the pink module (Table 2).

The size of the node shape represents the number of edges connected to that node, and the width of the edge represents the weight of the connection between two nodes. GO enrichment analysis was performed on the hub genes in the two modules separately; the respective module genes are shown in Figure 6. The gene functions of the green module were enriched in cyclic adenosine monophosphate (cAMP)-mediated signaling (GO:0019933), activation of protein kinase A activity (GO:0034199), regulation of inflammatory response (GO:0050727), cellular calcium ion homeostasis (GO:0006874), and other biological processes. In addition, the gene functions of the pink module were mainly enriched in the reactive oxygen species (ROS) metabolic process (GO:0072593), Wnt signaling pathway (GO:0060071), aging (GO:0007568), calcium ion transport (GO:0006816), and regulation of mitochondrial membrane potential (GO:0051881).

The KEGG pathway analysis showed that the green module was mainly enriched in the calcium signaling pathway (hsa04020), cAMP signaling pathway (hsa04024), and peroxisome proliferator-activated receptors (PPAR) signaling pathway (hsa03320). The pink module gene was mainly enriched in the transforming growth factor beta (TGF-β) signaling pathway (hsa04350). The corresponding data are presented in Supplementary Tables S4, S5. The screened hub gene list was inputted into the Cytoscape software. The interaction relationship graph between the co-expression network genes was made according to the weights of the genes, as shown in Figures 7A,B.

We performed RT-qPCR to detect gene expression levels in peripheral blood leukocytes collected from elderly AF patients and controls to validate the hub genes further. As shown in Figure 8A, we found that the expression of PTGDS, COLQ, ASTN2, VASH1, and RCAN1 genes from the green module was significantly increased in the AF group compared with those in the control (SR) group. In contrast, the expression of AMIGO2, RBP1, MFAP4, and ALDH1A1 genes in the pink module was decreased in the AF group compared with those in the SR group, suggesting that these genes can potentially be involved in the molecular mechanisms of AF.

Furthermore, based on gene expression levels, receiver operating characteristic (ROC) analysis was used to investigate whether the nine genes could predict AF in the elderly. As shown in Figure 8B, the area under the curve (AUC) and P value of PTGDS, COLQ, ASTN2, VASH1, and RCAN1 from the green module were 0.8272 (0.0193), 0.7840 (0.0423), 0.8765 (0.0071), 0.7901 (0.0380), and 0.8519 (0.0118), respectively. The AUC and P value of AMIGO2, RBP1, MFAP4, and ALDH1A1 in the pink module was 0.7654 (0.0576), 0.8519 (0.0118), 0.7531 (0.0703), and 0.8765 (0.0071), respectively (Figure 8C). These results indicate that the seven genes PTGDS, COLQ, ASTN2, VASH1, RCAN1, RBP1, and ALDH1A1 could be gene markers to differentiate AF patients and controls.