PDE3A and GSK3B as Atrial Fibrillation Susceptibility Genes in the Chinese Population via Bioinformatics and Genome-Wide Association Analysis

Background: Atrial fibrillation (AF) is the most common cardiac arrhythmia, with uncovered genetic etiology and pathogenesis. We aimed to screen out AF susceptibility genes with potential pathogenesis significance in the Chinese population. Methods: Differentially expressed genes (DEGs) were screened by the Limma package in three GEO data sets of atrial tissue. AF-related genes were identified by combination of DEGs and public GWAS susceptibility genes. Potential drug target genes were selected using the DrugBank, STITCH and TCMSP databases. Pathway enrichment analyses of AF-related genes were performed using the databases GO and KEGG databases. The pathway gene network was visualized by Cytoscape software to identify gene–gene interactions and hub genes. GWAS analysis of 110 cases of AF and 1201 controls was carried out through a genome-wide efficient mixed model in the Fangshan population to verify the results of bioinformatic analysis. Results: A total of 3173 DEGs were identified, 57 of which were found to be significantly associated with of AF in public GWAS results. A total of 75 AF-related genes were found to be potential therapeutic targets. Pathway enrichment analysis selected 79 significant pathways and classified them into 7 major pathway networks. A total of 35 hub genes were selected from the pathway networks. GWAS analysis identified 126 AF-associated loci. PDE3A and GSK3B were found to be overlapping genes between bioinformatic analysis and GWAS analysis. Conclusions: We screened out several pivotal genes and pathways involved in AF pathogenesis. Among them, PDE3A and GSK3B were significantly associated with the risk of AF in the Chinese population. Our study provided new insights into the mechanisms of action of AF.


Introduction
Atrial fibrillation (AF) is the most common nonbenign cardiac arrhythmia in clinical practice and is one of the leading causes of stroke, heart failure, cardiovascular disease, and sudden death [1,2]. With the increasing age of populations, AF incidence is increasing rapidly worldwide [3]. The strongest risk factor for AF is old age, along with gender, smoking, alcohol consumption, body mass index, hypertension, left ventricular hypertrophy, significant heart murmur, heart failure, and myocardial infarction [4,5]. However, the etiology and pathophysiologic mechanisms of AF are incompletely understood. In addition, new drugs specially designed for the therapy of AF remain suboptimal, and patients have to depend on antique antiarrhythmic drugs, such as amiodarone, sotalol, propafenone, and flecainide, which have limited efficacy and significant side effects [6].
Recently, numerous studies have suggested that AF cases in the general population have a significant genetic component, even beyond traditional risk factors [7]. Genomewide association studies (GWAS) have been applied to more than 30 million individuals worldwide and have discovered more than 100 distinct genetic loci associated with AF [8,9]. However, previous studies have had a variety of limitations. Combined with the results of the genetic studies thus far, genetic variation only accounts for approximately 40-50% of the heritability for AF. The sample size and ethnic diversity of GWASs have become increasingly large, but reproducibility of the identified association signals has become a serious challenge. The effects of many SNPs are difficult to verify in different ethnical populations. In addition, more than 95% of these GWAS variants are localized in noncoding regions, which may act through the effect of gene expression or pathway regulation. The main reason for these limitations is that conventional GWASs have only been able to identify SNPs that are associated with AF rather than causative, and there is always a lack of explanation for the biological significance behind the association.
Common approaches to further investigate the potential mechanisms of AF in a specific population are to combine the results of GWAS analysis with multiomics analysis of human atrial tissue by bioinformatics mining, including differentially expressed genes (DEGs) analysis, drug-gene interactions (DGIs) analysis, and pathway enrichment analysis. Bioinformatic analysis and studies using microarrays to measure gene expression can be employed to screen molecular markers in patients and healthy individuals. Microarray studies are commonly used to obtain gene expression profiles to uncover the pathogenesis of complicated diseases and for biomarker identification. In the last decade, a number of studies have tried to determine the transcriptomic changes in both AF patients and animal models using microarray technologies [10,11]. Several key pathways related to microRNAs (miRNAs), such as Ca 2+ -dependent signaling pathways, inflammatory and immune pathways, and apoptotic and cycle pathways, have been found [12]. Nevertheless, the results of bioinformatic analysis often lack the verification of real-world people.
In this study, we first integrated the multilevel biological information resources of AF and screened the set of candidate AF genes of biological significance based on the differentially expressed genes, the reported GWAS susceptibility genes, and the potential drug target genes. At the preliminary stage, the union of genes is extracted from the results of multiple data sets at the same molecular level, that is, as long as the results have been reported at least once, they are considered candidate AF-related genes. Then, we screened out the hub genes of AF through pathway enrichment analysis and the construction of the gene-gene interaction network of the candidate genes. We further carried out a GWAS analysis to verify the association between loci discovered by bioinformatics and AF in the Chinese population. Our study may be helpful in revealing the genetic etiology and pathogenesis underlying AF.

Identification of DEGs and Susceptibility Genes
Gene-expression profiles of AF were collected from the GEO database (www.ncbi. nlm.nih.gov/geo (accessed on 7 August 2022)) [13]. The GSE2240 data set includes the atrial myocardium tissues of samples from 10 patients with AF and 20 controls with sinus rhythm. The GSE128188 data set includes left and right atrial appendage tissues of samples from 5 AF patients and 5 controls with sinus rhythm. The GSE115574 data set includes the atrial tissues of samples from 15 patients with AF and 15 controls with sinus rhythm. All data in the present study were collected from public databases, so ethical approval from our institution was not needed.
The  [14]. p-values < 0.05 after Bonferroni correction were chosen as cut-off criteria. Gene expression values of | log 2 (fold change, FC) | > 0 were labeled as upregulated genes, and values < 0 were labeled as downregulated genes. Then, we obtained the list of AF genetic susceptibility sites and their Biomedicines 2023, 11, 908 3 of 17 mapping genes from the GWAS-catalog database and Open Targets database (accessed on 10 August 2022). Differentially expressed atrial fibrillation genes and susceptibility genes were combined as potential atrial fibrillation-related genes for subsequent analysis.

Identification of Potential Therapeutic Targets
The DrugBank database (https://go.drugbank.com/ (accessed on 13 August 2022)) [15] was used to search clinical drugs with indications for atrial fibrillation. The STITCH database (http://stitch.embl.de/ (accessed on 13 August 2022)) [16] was used to search the drug-gene interaction network of corresponding drugs and summarize the direct targets of drug action in the interaction network. The TCMSP database (https://old.tcmsp-e. com/index.php (accessed on 13 August 2022)) [17] was used to search the Chinese herbal medicines related to the treatment of atrial fibrillation and their pharmaceutical active ingredients. In order to include the effective components of traditional Chinese medicine for atrial fibrillation, only the active chemical components with set oral bioavailability (OB) greater than or equal to 30% and drug likeness (DL) value greater than or equal to 0.18 were incorporated. The collected active ingredients were used to find their targets in the TCMSP, and the targets were converted into the corresponding gene name in the UniProt database (https://www.uniprot.org (accessed on 13 August 2022)) [18].

Pathway Enrichment Analysis
GO enrichment analysis and KEGG enrichment analysis were performed on genes related to atrial fibrillation that could be potential therapeutic targets to explore the core functional pathways related to the pathogenesis and treatment of atrial fibrillation. p-value < 0.01 after Benjamini-Hochberg correction was used as the threshold to determine the enrichment effect of a gene subset in GO or KEGG entries. GO and KEGG enrichment analyses were performed with the R package clusterProfiler and enrichplot. The Cytoscape3.8.2 software (Cytoscape Consortium, California, USA) (ClueGO package verision 2.5.9, Laboratory of Integrative Cancer Immunology, Paris, France) [19] was used to construct pathway enrichment network map based on KEGG enrichment results to further screen the core functions or pathways related to the pathogenesis of atrial fibrillation.

Construction of Pathway Gene Network
The interaction relationships among gene targets under the pathway module were predicted through the String database (https://cn.string-db.org/ (accessed on 14 August 2022)) [20]. The score between nodes was set to 0.4, so that only the node interaction relationships greater than 0.4 can be included in the gene interaction network. Based on Cytoscape 3.8.2 software MCODE package (Verision 2.0.2, Gary Bader & Christian Lopes &Vuk Pavlovic, Toronto, Canada), the gene interaction network was constructed for each KEGG pathway network gene to mine the hub gene sets under each pathway network.

GWAS Study Design and Subjects
The GWAS analysis consisted of 110 cases of AF and 1201 controls free of AF from the Fangshan Family-based Ischemic Stroke Study in China (FISSIC) [21]. FISSIC is an ongoing community-based case-control genetic epidemiological study that started in June 2005, enrolling families in Fangshan District, a rural area located southwest of Beijing, China. The inclusion criteria for the subjects were as follows: (1) age older than 18 years at enrollment; (2) variables of sex, age, or AF condition not missing; and (3) subjects without single-gene hereditary disease or cancer. The diagnosis of atrial fibrillation was confirmed by a second-class or higher-class hospital. Atrial fibrillation (ICD-10 code I48) was defined as any event with a date of occurrence before the participant's first visit for recruitment into the study. This study was approved by the Ethics Committee of the Peking University Health Science Center (Approval number: IRB00001052-13027), and written informed consent was provided by all participants.

Genotyping
DNA was extracted using a LabTurbo 496-Standard System (TAIGEN Bioscience Corporation, Taiwan, China). In addition, the purity and concentration of DNA were measured using ultraviolet spectrophotometry. Samples within the study were genotyped at the Capitalbio Technology Corporation using the Illumina ASA Chip. They were genotyped in 5 batches, grouped by origin of the samples, and with a balanced case-control mix on each array. Quality control (QC) [22] was performed on each sample, including >95% variant call rate, consistency between genotyped sex and the investigated sex, <3 SDs heterozygosity, consistency between IBDs and the investigated kinship, <5% Mendel errors, and no significant deviation from PCAs of ancestral background. QC was performed on each call set, including >95% sample call rate, Hardy-Weinberg equilibrium p > 1 × 10 −6 , minimum allele frequency (MAF) > 1%, and <10% Mendel errors. All QC was conducted using PLINK 1.9 software (https://www.cog-genomics.org/plink/ (accessed on 14 August 2022).

Statistical Analysis
Genome-wide association testing was performed using GEMMA (genome-wide efficient mixed model association) [23] software based on mixed effects model. The locus effect was decomposed into fixed effects on families and random effects on individuals by constructing the phylogenetic matrix, and the interindividual correlation within families was adjusted. Sex, age, and the first ten principal components were adjusted as covariates.
To correct for multiple testing, a genome-wide significance threshold of p < 1 × 10 −8 was performed. We inspected Manhattan plots and Q-Q plots for spurious associations and quantile-quantile plots to identify genomic inflation.
Analyses in addition to GEMMA were conducted using R version 4.2.1 (R Foundation for Statistical Computing, Vienna, Austria).

Identification of DEGs and Susceptibility Genes
A total of 3173 DEGs were screened from the three GEO data sets: expression of 1757 genes was upregulated and the expression of 1432 genes was downregulated between AF patients and controls (Table 1). Among them, 22 genes were differentially expressed in all three data sets ( Figure 1a). There were 16 genes with different expression regulatory directions in at least two data sets. The clusters of all DEGs are displayed in Figure 1a. A total of 356 lists of AF genetic susceptibility sites and their mapping genes were obtained from the GWAS-Catalog database and Open Targets database. Among all the DEGs in atrial fibrillation (n = 3173), variants in 57 genes were found to be significantly associated with the risk of AF (p < 5 × 10 −8 ). By combining GWAS and RNA expression information, 3472 genes were ultimately identified as potential AF-related genes, which are shown in Figure 1b.

Identification of Potential Therapeutic Targets
With atrial fibrillation as key words, 22 kinds of clinical drugs conforming to indications were searched based on the DrugBank database. A total of 137 direct targets were found in drug-gene interaction networks based on the STITCH database. A total of 439 herbal medicines related to the treatment of atrial fibrillation were found in the TCMSP database. Among them, 406 active ingredients were collected, which can be mapped to 278 targets. We intersected the targets with the potential AF-related genes. As a result, a total of 75 AF-related genes can be used as potential therapeutic targets, including 61 Chinese medicine targets and 26 chemical drug targets, as shown in Table 2 and Figure 2.  A total of 3472 genes were identified as potential AF-related genes, and 57 of them were overlapped in both sets.

Identification of Potential Therapeutic Targets
With atrial fibrillation as key words, 22 kinds of clinical drugs conforming to indications were searched based on the DrugBank database. A total of 137 direct targets were found in drug-gene interaction networks based on the STITCH database. A total of 439 herbal medicines related to the treatment of atrial fibrillation were found in the TCMSP database. Among them, 406 active ingredients were collected, which can be mapped to 278 targets. We intersected the targets with the potential AF-related genes. As a result, a total of 75 AF-related genes can be used as potential therapeutic targets, including 61 Chinese medicine targets and 26 chemical drug targets, as shown in Table 2 and Figure 2.

Pathway Enrichment Analysis
Pathway enrichment analysis was performed on the 75 AF-related genes above. GO enrichment analysis showed that 919 gene subsets were significantly enriched under specific biological processes. The first 20 biological processes with significant enrichment were selected to generate bubble maps (Figure 3a). In addition, KEGG enrichment analysis revealed significant enrichment of 79 gene subsets in specific pathways. The first 20 significantly enriched pathways were selected to generate the bubble map ( Figure 3b).

Pathway Enrichment Analysis
Pathway enrichment analysis was performed on the 75 AF-related genes above. GO enrichment analysis showed that 919 gene subsets were significantly enriched under specific biological processes. The first 20 biological processes with significant enrichment were selected to generate bubble maps (Figure 3a). In addition, KEGG enrichment analysis revealed significant enrichment of 79 gene subsets in specific pathways. The first 20 significantly enriched pathways were selected to generate the bubble map ( Figure 3b).

Construction of Pathway Gene Network
GO enrichment network maps suggest clustering of 919 distinct biological processes into 80 functional categories. The main biological functions of AF-related genes mainly include the regulation of polysaccharides and transmembrane transporters, muscle con-

Construction of Pathway Gene Network
GO enrichment network maps suggest clustering of 919 distinct biological processes into 80 functional categories. The main biological functions of AF-related genes mainly include the regulation of polysaccharides and transmembrane transporters, muscle contraction, and the negative feedback regulation process of catecholamine, glutaminergic compounds, and synaptic transmission. The pathway enrichment network constructed based on the KEGG enrichment results is shown in Figure 4. The nodes represent enriched biological processes or pathways. The node lines represent the number of common genes between biological processes or pathways, and the color indicates which functional group the node belongs to. According to shared genes among the pathways, 79 significant pathways were classified into 7 major pathway networks. The gene interaction network was constructed for each KEGG pathway network gene to mine the hub gene set under each pathway network. The hub genes under the seven subnetworks included the oxidative stress and cellular signaling pathway, hormone regulation pathway, cell adhesion pathway, tumor inhibition pathway, cell cycle regulation pathway, lipid metabolism and inflammatory pathway, and mental illness-associated pathway (Table 3 and Figure 5). A total of 35 common hub susceptibility genes were obtained by constructing the pathway-gene network. Twenty-six of them belonged to the oxidative stress and cellular signaling pathways.  The gene interaction network was constructed for each KEGG pathway network gene to mine the hub gene set under each pathway network. The hub genes under the seven subnetworks included the oxidative stress and cellular signaling pathway, hormone regulation pathway, cell adhesion pathway, tumor inhibition pathway, cell cycle regulation pathway, lipid metabolism and inflammatory pathway, and mental illness-associated pathway (Table 3 and Figure 5). A total of 35 common hub susceptibility genes were obtained by constructing the pathway-gene network. Twenty-six of them belonged to the oxidative stress and cellular signaling pathways.  TP53, ESR2, EGR1, HSPA5,  CREB1, CAV1, AKT1, TGFB1, GSK-3β,  MAPK1, NFE2L2, CCNA2, IGF2,  BCL2L1, CDKN1A, PTGS2, PTEN, CDK4,  ERBB2, MMP9, AR, IGF1R   A  2  3  3  3 AGT, SPP1, AGTR1   Table 3.

GWAS Analysis
A total of 496,798 genetic variants were tested after quality controls. Principal component analysis (PCA) revealed that participants in the present study are genetically East Asian, and there are no individuals who deviate significantly from their ancestral genetic background ( Figure S1). The GWAS association analysis of AF in the Fangshan population revealed 126 AF-associated loci at genome-wide significance (p < 1 × 10 −8 ) ( Figure 6 and Table 4). The significance level accounts for multiple testing of independent variants with MAF ≥ 0.1% using a Bonferroni correction. p values (two-sided) were derived from a genome-wide efficient mixed model with the least-squares method.  Table 3.

GWAS Analysis
A total of 496,798 genetic variants were tested after quality controls. Principal component analysis (PCA) revealed that participants in the present study are genetically East Asian, and there are no individuals who deviate significantly from their ancestral genetic background ( Figure S1). The GWAS association analysis of AF in the Fangshan population revealed 126 AF-associated loci at genome-wide significance (p < 1 × 10 −8 ) ( Figure 6 and Table 4). The significance level accounts for multiple testing of independent variants with MAF ≥ 0.1% using a Bonferroni correction. p values (two-sided) were derived from a genome-wide efficient mixed model with the least-squares method.

Bioinformatic Results of PDE3A and GSK-3β
We then sought to link risk variants to candidate genes by their effect on gene expression levels or potential drug targets based on the previous bioinformatic analysis to further enhance the biological understanding of the atrial fibrillation-associated loci. We found two genes overlapping in two approaches. PDE3A (p = 4.98 × 10 −5 ) and GSK-3β (p = 0.031) were first identified as upregulated DEGs between AF and sinus rhythm people in atrium tissue in the GSE2240 data set. Neither PDE3A nor GSK-3β were identified as GWAS susceptibility genes, but were included in subsequent analyses as potential AF-related genes. PDE3A was then identified as a potential drug target for AF. The drug active ingredients that interact with PDE3A are CGMP-inhibited 3 ,5 -cyclic phosphodiesterase A.
In pathway enrichment analysis, GSK-3β was enriched into pathways related to oxidative stress and hormone regulation, and was identified as a hub gene in the gene-gene network. The most significantly enriched pathway was the PI3K-Akt signaling pathway ( Figure 7). The PI3K-Akt pathway is an intracellular signal transduction pathway that promotes metabolism, proliferation, cell survival, growth, and angiogenesis in response to extracellular signals. This is mediated through serine and/or threonine phosphorylation of a range of downstream substrates. In the pathway, AKT1 regulates the phosphorylation of the Ser9 site of GSK-3β, which leads to its inactivation. Activated PI3K converts PIP2 to PIP3. These PIPs then mop up PDK1 and Akt to the cell membrane. When PDK1 and Akt are taken to the cell membrane, Akt gets activated and phosphorylated. Overexpression of phosphatase and tensin homolog deleted on chromosome 10 (PTEN) can inhibit AKT1 phosphorylation and further activate GSK-3β. GSK-3β activity inhibits the binding of GSK-3β to BCL2 and then promotes the activation of autophagy. GSK-3β can also phosphorylate MAPK1 kinases, which is implicated in fibrogenesis.

Discussion
In this study, we applied a complementary strategy to combine the results of GWAS analysis with bioinformatics data mining in multiomics, and found that variants of potential therapeutic target PDE3A and key mediator gene GSK-3β of AF were significantly susceptible to AF in the Chinese population.

Discussion
In this study, we applied a complementary strategy to combine the results of GWAS analysis with bioinformatics data mining in multiomics, and found that variants of potential therapeutic target PDE3A and key mediator gene GSK-3β of AF were significantly susceptible to AF in the Chinese population.
Phosphodiesterase 3A (PDE3A) encodes a member of the cGMP-inhibited cyclic nucleotide phosphodiesterase (cGI-PDE) family [24]. cGI-PDE enzymes hydrolyze both cAMP and cGMP, and play critical roles in many cellular processes by regulating the amplitude and duration of intracellular cyclic nucleotide signals [25]. The encoded protein mediates platelet aggregation and also plays important roles in the cardiac β-AR/AC/cAMP/PKA axis by regulating vascular smooth muscle contraction and relaxation. Our study first explored whether PDE3A was significantly upregulated in AF patients, which was consistent with previous research. Bernardo Dolce et al. demonstrated that in patients with persistent atrial fibrillation, the force responses to 5-HT are blunted, but they can be recovered after inhibition of PDE3 [26]. This suggests that the change of PDE3A expression may cause systolic dysfunction of atrial muscle and thus participate in cardiac remodeling. Combined with the DrugBank, STITCH, and TCMSP databases, we then identified PDE3A as a potential drug target for AF. It was also reported that inhibitors of the encoded protein of PDE3A may be effective in treating congestive heart failure [27,28]. In GWAS analysis, we found that minor allele A of the intron variant PDE3A rs11613698 was observed to increase the risk of atrial fibrillation in the Chinese population. Although previous GWAS studies in the Chinese population did not identify PDE3A as a risk gene for AF, some indirect associations have been reported. Zun Wang et al. applied a novel metaCCA method on the GWAS summary statistics data of stroke and AF and found that PDE3A was a potential pleiotropic gene, which may affect ischemic or hemorrhage stroke through multiple intermediate factors such as MAPK family [29]. It was also reported by Carmen Sucharov et al. that 29 polymorphisms in the human PDE3A gene promoter regulated PDE3A gene expression, and had further effects on the electrophysiological activity of the myocardium [30]. Therefore, mutation of PDE3A may change the expression level of its downstream proteins, and affect the risk of atrial fibrillation by mediating cardiac remodeling.
The protein encoded by glycogen synthase kinase 3β (GSK-3β) is a serine-threonine kinase, belonging to the glycogen synthase kinase subfamily. It is involved in energy metabolism, neuronal cell development, and body pattern formation [31]. Numerous studies have indicated that GSK-3β can be phosphorylated and inhibited by protein kinase C (PKC) and then regulate a wide variety of cardiac transcription factors [32]. Recent studies suggest that calcium channel modification is a possible mediator of the association between GSK-3β and AF. Yan Wang et al. developed an animal model and found that ethanol can inhibit GSK-3β through enhanced phosphorylation, thereby leading to upregulation of T-type calcium channels (TCCs) and increased AF susceptibility [33]. In addition, it was reported that GSK-3 can negatively regulate the sarco (endo)plasmic reticulum Ca 2+ -ATPase (SERCA) pump, a key regulator of Ca 2+ uptake in the heart [34]. In a prospective observational study, SERCA levels in peripheral lymphocytes were reported to be associated with the outcome of pericardial ablation in patients with persistent AF, and lower levels of SERCA expression could predict the recurrence of AF after pericardial ablation [35]. Taken together, overexpression of GSK-3β can cause abnormal regulation of calcium ions in cardiomyocytes by inhibiting SERCA levels, leading to myocardial electrical remodeling, and thus the occurrence and maintenance of AF.
In this study, we first found that GSK-3β had substantially higher levels of expression in patients with AF than in controls. Then, GSK-3β was identified as a hub gene in oxidative stress and cellular signaling pathways. There are several mechanisms that produce ROS in cardiac myocytes, including mitochondria, NADPH oxidase, uncoupled NO synthase, and xanthine oxidase [36,37], which increase oxidative stress and promote cardiac fibrosis. Many studies have reported that patients with AF have a decrease in antioxidant-related gene expression and an increase in ROS-related gene expression [38]. Consistent with our pathway enrichment results, previous studies have reported that regulation of oxidative stress-related gene expression was functionally associated with PI3K/AKT signaling, which is a key profibrotic element in various tissues and was reported to be capable of activating atrial fibroblasts to differentiate into myofibroblasts [39]. Several experimental studies have sought targets to inhibit atrial remodeling by affecting oxidative stress, inflammation, and the PI3K/Akt signaling pathway [40,41].
Our identification of the GSK-3β as a hub gene may help to provide a new target for the etiological study of oxidative stress pathway in AF. The minor allele G of the intron variant GSK-3β rs796944992 was observed to increase the risk of atrial fibrillation. Although there have been no previous reports on the GSK-3β variation and the risk of AF in other populations, many studies have reported that genetic alterations in GSK-3β are associated with various diseases mediated by oxidative stress pathways, including myocardial ischemia [42], myocardial infarction [43], and Alzheimer's disease [44]. Our results suggest that variants of GSK-3β directly confer a risk of AF only in the Chinese population, indicating strong population heterogeneity in the genetics of AF. The mechanism by which GSK-3β mediates atrial fibrillation through the oxidative stress pathway can be further explored in future studies.
There are complex mechanisms for the occurrence and maintenance of AF, including several pathways not identified in this study, such as overinflammation and epigenetic modification. Inflammation can alter the atrial electrophysiology and structure to increase the vulnerability to AF. In a study on patients with new-onset AF, the early recurrence of AF was related to inflammatory markers, and inflammatory markers were associated with development of permanent AF [45]. However, it was reported that although oral antioxidant treatment (α-lipoic acid, ALA) reduced serum levels of common markers of inflammation in ablated patients, ALA does not prevent AF recurrence after an ablative treatment [46]. This suggests that treatment targeting inflammatory biomarkers alone may not able to revert cardiac remodeling. This may also be the reason why few genes related to inflammation were screened out as hub genes in this study. As epigenetic regulators, miRNA plays an important role in cardiac development, and the dysregulation of miRNA expression is related to cardiac remodeling. It was reported that catheter ablation was related to miRNA modulation. Several miRNAs have been reported to be able to assess and predict the risk of recurrence in patients with AF after ablation [47]. In the future, the idea of this study could be applied to miRNA bioinformatics mining, which may help to identify miRNA drug targets with clinical value after AF ablation.
We acknowledge some limitations in our study. First, three GEO data sets were included to detect DEGs in AF. There were considerable differences in their study designs hindering straightforward comparison and merging with the studies [48]. The origin of mRNAs in tissue is inconsistent, and mRNA regulation in different tissues may be contradictory. Second, although bioinformatics data mining and GWAS analysis yielded 75 hub genes and 126 loci, there were only two overlapping genes between the two methods. This may be mainly due to the fact multiomics data on AF in the Chinese population are difficult to obtain, so the genetic background of the included open data is different from our GWAS population. These factors need to be carefully considered to avoid misinterpretation of the findings. Finally, the sample size of the GWAS analysis was limited and the ratio of cases to controls was imbalanced in this study, so the main purpose of GWAS analysis was to validate AF-related genes obtained from biological information analysis in the Chinese population. This is a strategy similar to that of candidate gene association studies, and GWAS analysis itself serves as a validation function. We hope to obtain a more independent AF cohort for further analysis in future studies.

Conclusions
In conclusion, this study is the first systematic report on the screening and verification of the association of PDE3A and GSK-3β with the risk of atrial fibrillation in the Chinese population, showing that PDE3A is a potential drug target for AF and GSK-3β is a hub gene in the gene-gene network of pathways related to oxidative stress and hormone regulation. Our study provides new insights into AF mechanisms in the Chinese population.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/biomedicines11030908/s1, Figure S1. Principal component analysis of Fangshan population and various ethnic groups. The figure shows the first two principal components to reveal population structure. The distribution of PCAs in Fangshan population and East Asian population is basically the same. Figure S2. LocusZoom plots PDE3A rs11613698 and GSK3B rs796944992. SNPs are colored based on their correlation (r 2) with the labeled top SNP. Arrows on the horizontal blue lines show the direction of transcription.
Author Contributions: Conceptualization, D.C.; methodology, X.L.; software, Z.Z. and X.L.; formal analysis, Z.Z. and X.L.; investigation, Y.W., Y.Z. and L.Y.; data curation, Y.W., Y.Z. and L.Y.; writingoriginal draft, Z.Z.; supervision, D.C. and X.W.; project administration, D.C. and X.W.; funding acquisition, X.W. All authors have read and agreed to the published version of the manuscript. Informed Consent Statement: Informed consent was obtained from all subjects involved in the study. Written informed consent has been obtained from the patient(s) to publish this paper.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to the policy of the Ethics Committee of the Peking University Health Science Center.

Conflicts of Interest:
The authors declare no conflict of interest.