Analysis of Genes Involved in Persistent Atrial Fibrillation: Comparisons of ‘Trigger’ and ‘Substrate’ Differences Open Access

Atrial fibrillation (AF), the most common dysrhythmia, is a major risk factor for cardiovascular (CV) mortality and morbidity resulting from thromboembolism and hemodynamic changes owing to the high-frequency fibrillation of the atrial wall and stagnation of blood flow [1]. Currently, there are 30 million cases of AF worldwide and approximately one-third of cases of CV-related mortality and morbidity are caused by AF [2, 3]. Recently there has been an increasing focus on the relationships between subclinical episodes of AF and the risk of silent cerebral infarct (SCI) and stroke. Marfella’s group reported findings related to silent episodes of AF and observed that subclinical episodes of AF may increase the risk of SCI and stroke, especially in patients with type 2 diabetes [4]. Emerging evidence has linked heart failure and comorbid AF to increased CV-related mortality [5]. Notably, AF also carries a greater risk of poor cardiovascular prognosis in patients with HF [5], indicating the importance of early AF detection, particularly to improve the quality of life of HF patients.

However, symptomatic patients have limited treatment options and currently available antiarrhythmic drugs are generally associated with potential adverse effects or toxicity. Furthermore, catheter ablation and surgical MAZE have limited efficacy in patients with persistent AF [6]. Surgical ablation is associated with a recurrence rate of ∼30% and potential complications in 6% [6]. Similarly, in the synergy ablation system, the ABLATE trial demonstrated that the Cox Maze-IV procedure is effective in 76% and the rate of major adverse events is greater than 9% [7]. Preventing AF recurrence after catheter ablation is a major challenge, and a clinical trial conducted by Sardu’s group illustrated that antioxidant treatment after ablation could reduce the serum levels of inflammatory factors, but regrettably it did not improve the poor prognosis in the first year after catheter ablation [8]. That study underscores the fact that the pathogenesis of AF and its underlying molecular mechanisms are not completely understood. Recently, a series of preclinical studies suggested that gene-based approaches could be used for the diagnosis and treatment of AF [9, 10]. Most studies have had encouraging results that indicate the potential of gene-based therapies to combat AF. Amit’s group utilized an adenoviral-mediated delivery of KCNH2 (negative mutant of the IKr potassium channel) into a porcine AF model to enable epicardial atrial gene transfer and pacemaker implantation, as well as electrophysiological and molecular studies [9]. The authors found that KCNH2 transfer prolonged atrial action potential duration and prevented AF [9]. Similarly, Zhang’s group found that miR-206 modulated intrinsic cardiac autonomic nerve remodeling in AF via right atrial tachypacing via the regulation of superoxide dismutase 1 (SOD1) in vivo, and this may be a potential therapeutic target for AF because it shortens the atrial effective refractory period [10]. Furthermore, many loci or target genes associated with ion channels, autonomic innervation, gap junctions, and substrate modifiers have been identified in preclinical studies to attenuate fibrosis, autonomic nerve sprouting, and AF duration [11-13]. However, inflammation and the immune response, oxidative stress, and genetic mechanisms do not fully explain the development of AF, and autonomic system dysfunction may exert direct effects on the initiation and maintenance of AF due to heterogeneous changes in the atrial electrophysiology [14, 15]. In a large clinical study, Rizzo’s group offered compelling evidence of the relationships between autonomic dysfunction and silent AF in type 2 diabetes, thus revealing a potential prognostic factor and principle for target therapies that potentially offer greater specificity and efficacy in the control of AF [15]. In sum, AF may be a final common endpoint of atrial remodeling or electrical, structural, cardiac autonomic nerves, and ion channel changes that cause irregular atrial rhythms and unique spatial heterogeneous conduction [16]; thus it is essential to develop a clinical scheme for the early and accurate detection of AF and appropriate interventions [13]. At present, gene-based transformation therapy and tailored treatments may be promising solutions because of the advent of genome-wide association studies.

We downloaded the expression microarray data (CEL data) from the GSE41177 dataset from Gene Expression Omnibus (GEO) (http://ncbi.nlm.nih.gov/geo/) [17]. The GSE41177 dataset included 16 pairs of pulmonary veins and surrounding left atrial (LA-PV) junctions and left atrial appendage (LAA) specimens obtained from persistent AF patients after valvular surgery. We measured gene expression using GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array (Affymetrix Inc., Santa Clara, CA). In addition, platform information, probe annotation, and clinical information for all samples were pooled for further analysis.

We processed RAW data (CEL file) using the Affy R package (http://www.bioconductor.org/packages/release/bioc/html/affy.html), provided by the Bioconductor project. We used linear models for microarray data (LIMMA) package in R language, Robust multi-array analysis method, and k-nearest neighbors algorithm to identify differentially expressed (DE) genes between the LA-PV and LAA specimens [18]. We also assessed background correction, quantile normalization, and probe summarization. Averages were considered as gene expression values if the multiple probes were mapped to the same gene symbol. We used a false discovery rate to calculate fold-changes (FC) and also the Benjamini-Hochberg method [19]. Gene expression values of log2FC > 1 and adjusted p < 0.05 were chosen as thresholds. Hierarchical clustering analysis was used to build a hierarchy of clusters of the DE genes [20]. Subsequently, we performed a visualization analysis of chromosome loci and applied it to the initiation site information to map gene regions and obtain biological information on the DE genes. The database of hg19_ refGene was referenced in the gene mapping analysis [21].

Weighted correlation network analysis (WGCNA), using the Eigengene network methodology, was applied to study the biological networks and select clinical features based on weighted correlations between quantitative variables (normalized gene expression data) [22, 23]. Here, we chose the first 25% least mean square deviation of the expression variable to guarantee the stability and accuracy of the bioinformatics statistics to a certain extent to perform further co-expression network analysis. Soft thresholding, based on the scale-free topology criterion, was employed for further analysis to preserve the continuous nature of co-expression information [22, 23]. The dynamic branch cutting approach was applied to identify the gene modules, and the Eigengenes method was used to draw characteristic values based on more complex predictive models, including decision trees and Bayesian networks [23]. Here, the smallest value of the soft thresholding parameter was determined by a condition of high-scale independence. Heat maps were generated to evaluate the connection strengths and module Eigengene of the correlation network. The analysis strategy is presented in Fig. 1.

We used the online Search Tool for the Retrieval of Interacting Genes (STRING database, Version 10.5; http://string-db.org/) to construct a PPI network based on uniquely comprehensive coverage and predictive function of genome-wide data [24]. We used a reliability threshold of a combined score of > 0.4 to construct the PPI network and used Cytoscape software (Version 3.5.1; http://cytoscape.org/) to visualize and analyze the biological networks [24]. Thus, the degree of each node was used to screen “hub” genes that interacted with other genes or proteins.

The Database for Annotation, Visualization, and Integrated Discovery (DAVID) is a web-based analytical tool used to extract comprehensive biological information about candidate genes that can then provide functional interpretation of large lists of genes derived from genomic studies [25]. We used DAVID resources to perform Gene Ontology (GO) functions and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. We assessed GO terms in biological processes, molecular functions, and cellular components. GO terms and KEGG maps of biological functions were enriched with a threshold p-value of < 0.05.

After normalization, 152 genes were found to be differentially expressed between LA-PVs samples and LAA samples. Among these DE genes, 104 were down-regulated and 48 were up-regulated. A heat map was generated to assess the expressed values of the DE genes to distinguish LA-PV samples and LAA samples in AF patients. Interestingly, the DE genes included more down-regulated genes in LA-PV than in LAA (Fig. 2A). Additionally, genetic mapping information showed that more DE genes were associated with the negative log of the corresponding p-value between 10 × 10^–8 and 15 × 10^–8, and 7 hub genes were at more than 10 × 10^-8 (Fig. 2C).

To confirm the WGCNA study, we performed cluster analysis among the top 5121 genes (first 25%), and thus GSM1005428, GSM1005444, GSM1006247, and GSM1006249 were excluded from the sample cluster analysis (Fig. 4A). In addition, the power of β = 12 was chosen to ensure that the scale independence was up to 0.78 with higher mean connectivity (Fig. 4B). As such, we identified 16 gene clusters/modules of co-expressed genes based on the dynamic branch cutting approach (the merge cut height was 0.25) (Fig. 4D). In addition, following the module-trait relationships analysis, the grey module was primarily associated with AF duration (p = 0.005), the cyan module was mainly associated with left ventricular ejection fraction (LVEF) (p = 0.002), the blue module was significantly associated with left atrium diameter (LAD) (p = 0.001), and the magenta module was particularly important for ischemic heart disease, hypertension, and diabetes comorbidities (the p-values were 0.005, 0.02, and 0.05, respectively) (Fig. 4C). Interaction analysis of the co-expression module demonstrated non-statistically significant results with a brightness area among different modules and a higher scale of independence among these clusters (Fig. 4F). Furthermore, the Eigengene adjacency plot also illustrated a significant difference in connectivity effects (Fig. 4E). The sample clinical information is depicted in Fig. 4A.

KEGG pathway maps of biological functions associated with up- and down-regulated DE genes were obtained from the DAVID web-based search tool, respectively. Table 1 shows the up- and down-regulated genes and the pathways in which they interact. We assessed the co-pathways of the DE genes and identified hsa04512: extracellular matrix (ECM)-receptor interaction (up-regulated: COMP, THBS2, THBS4; down-regulated: COL6A6, RELN, TNN) as possible points of similarity between the LA-PV and LAA samples, reflecting the special functional roles of ‘trigger’ and ‘substrate’ (these data are shown in Table 1). The expression of hsa04512-related DE genes is depicted in Fig. 5A.

In the WGCNA study, the grey module genes were mainly enriched in hsa04512: ECM-receptor interaction (including 10 genes with p-values of 1.59E-04) and further revealed the biological mechanisms in relation to AF duration. The blue module genes were mainly enriched in hsa04062: chemokine signaling pathway (including 27 genes with p-values of 2.13E-07) and closely associated with LAD. In addition, the cyan module genes were significantly associated with cfa04110: cell cycle (including 9 genes with p-values of 3.51E-10) related to LVEF, and the magenta module genes were primarily associated with ssc04668: TNF signaling pathway (including 9 genes with p-values of 2.68E-07) may be significantly associated with genetic traits of ischemic heart disease, hypertension, and diabetes (Fig. 5C).

We performed GO enrichment analysis of the up-regulated and down-regulated DE genes following the LIMMA study. For the up-regulated genes, we obtained several clusters of interest based on official gene symbol information. Among the biological processes, the first cluster was significantly associated with the Wnt signaling pathway, including the non-canonical Wnt signaling pathway, positive regulation of the canonical Wnt signaling pathway, and negative regulation of the canonical pathway. Among the molecular functions, our results showed that the up-regulated DE genes were significantly associated with the binding of heparin. In cellular components enrichment analysis, the terms of the extracellular region and extracellular space were also enriched. Additionally, different terms were significantly enriched among the down-regulated genes, including regulation of blood pressure, notch signaling path, and PI3K signaling. However, following molecular function analysis, calcium ion binding was significantly enriched in the down-regulated DE genes. Furthermore, the cellular component terms of the extracellular regions, extracellular space, and integral components of plasma membranes were significantly enriched in our study (Fig. 3).

In the WGCNA study, the genes in the grey module, which were associated with AF duration, were mainly enriched in association with the positive regulation of inflammatory responses. Similarly, the genes in the blue module were significantly associated with immune responses, illustrating the potential molecular mechanism in relation to LAD. In addition, we also found that the term mitotic cytokinesis was significantly enriched in the cyan module genes, which were closely associated with LVEF. Nevertheless, regarding the genetic traits of comorbid ischemic heart disease, hypertension, and diabetes, the terms immune response and fat cell differentiation were primarily identified among the genes in the magenta module (Fig. 5B).

Based on the STRING database, we constructed a PPI network containing 152 DE genes with a combined score of > 0.4 (Fig. 2D). There were 66 nodes with network density and heterogeneity of 0.04 and 0.88, respectively. After calculating the degree of each DE gene, 7 genes were considered central to the maintenance of AF with a degree of > 6. In addition, leptin (LEP), FBJ murine osteosarcoma viral oncogene homolog (FOS), endothelin 1 (EDN1), neuromedin U (NMU), calbindin 2 (CALB2), tachykinin (TAC1), and pro-platelet basic protein (chemokine (C-X-C motif) ligand 7) (PPBP) could interact with each other in the PPI network (as shown in Fig. 2D). The expression levels of these hub genes are shown in Fig. 2B.

The LA-PV junction is thought to be an important ‘driver’ or ‘trigger’ region for the maintenance and immediate recurrence of AF and short bursts of LA-PV tachycardias. LAA may act as a “substrate” or “tract” in AF. Pison’s group reviewed the potential mechanisms of sustained AF triggered by ectopic beats arising from muscle sleeves of the LA-PV to drive focal activity or local re-entry emanating from spiral waves from a curved central core. The mechanisms behind the maintenance of AF are increasingly being attributed to these re-entry substrates [26]. Romero’s group reported that the LA-PV is the key to AF maintenance—they used PV isolation to treat AF [27]—and further alterations of atrial substrates related to atrial inflammation and fibrosis may reduce the atrial potential duration (APD) and cause electrical heterogeneity and greater vulnerability of single or multiple re-entry beats and random-wave front propagation [28-30]. However, LAA is considered the primary source of thromboemboli in patients with persistent or permanent AF [31].

To develop novel strategies or individualized treatment, it is essential to discriminate region-specific gene expression and understand the pathophysiological mechanisms associated with AF maintenance. Recently, Yeh’s group used a genome-wide approach to identify region-specific gene expression in paired LA-PV junction and LAA specimens of patients with AF [32], and found that the LA-PV junction had a higher expression of arrhythmia-related genes; conversely, several thrombosis-related genes were up-regulated in LAA samples [32]. In our LIMMA study, 152 DE genes (48 up-regulated and 104 down-regulated genes) and 7 hub genes were identified, and the ECM-receptor interaction and chemokine signaling pathway were both significantly enriched in the up-regulated DE genes, as were both the ECM-receptor interaction and PI3K-Akt signaling pathway in the down-regulated DE genes. Moreover, to assess the additional relationships between the gene modules and external clinical information, we carried out a WGCNA study and found that ECM-receptor interactions, chemokine signaling pathways, and the cell cycle may be potentially targeted for AF duration, LAD, and LVEF, respectively, at least with a soft threshold cutoff of 12.

In this study, both differential gene expression and co-expression network analyses powerfully revealed an important regulatory role for the ECM-receptor interaction pathway in the comparison of trigger and substrate, which turned out to be closely related to AF duration. Recently, Husser’s group performed clinical and genetic predictor analysis based on pathway enrichment analysis of GWAS (Genome-wide Association Study) data, and demonstrated that the ECM-receptor interaction pathways were significantly associated with left atrial enlargement and AF type and response to AF catheter ablation [33]. Presumably, ECM remodeling revealed complex relationships between cell-to-cell interactions and communication, and was validated as a central mediator that tightly links ion channels with the development of an ectopic pacemaker [33]. Although the previous studies and our results have provided evidence to show that ECM interaction is a potential target, the precise molecular mechanisms of the ECM in the progressive nature of AF are still unknown.

LEP is a hormone synthesized by adipocytes that has signaling ability and regulates food intake and energy metabolism. Abd El-Aziz’s group compared circulating LEP and genotyped the LEP gene and leptin receptor (LEPR) gene and their gene polymorphisms among 100 coronary artery disease (CAD) patients and 100 healthy subjects [34]. The group reported greater LEP expression in CAD patients, and gene polymorphisms were more frequent in the LEP gene (AA genotype) and the leptin receptor gene (RR genotype) in CAD patients who experienced heart failure [34]. In addition, Li’s group demonstrated that LEP-related obesity may be related to cardiac contraction defects resulting from damage to nicotinamide adenine dinucleotide phosphate (NADPH) oxidase activation, oxidative modification of sarco(endo) plasmic reticulum Ca²⁺-ATPase, and myosin heavy chains [35]. Similarly, based on the rapidly growing data from experimental studies, the overexpression of FOS may play a central role in heart disease because of its association with inflammation and cardiac remodeling. Zhang’s group illustrated that the expression of c-fos in cardiomyocytes and serum is rapidly increased at days 1, 3, 7, and 10 after coronary ligation [36]. Drosos’s group assessed elevations in the expression of LEP protein in aortic root and coronary arterial perivascular adipose tissue (C-PVAT) and internal mammary arterial perivascular adipose tissue (IMA-PVAT) and found higher expression of hypoxia-inducible factor-1α (HIF-1α) and Fos-like antigen (FOSL)2 in hypoxia C-PVAT compared with IMA-PVAT via regulation of the hypoxia-inflammation-fibrosis axis in atherosclerosis, which reflected an increased risk of atherosclerotic plaque and CV-related death [37]. In addition, EDN1 is an important intermediate for guiding the intrinsic ability of these progenitors to differentiate and self-organize into a functional vasculature [38], and has a pivotal role in the maintenance of normal contractile function via the superoxide-metalloproteinase 9-(Mmp9) cascade. Moreover, less than 35% decrease in EDN1 expression is sufficient to cause significant malignant cardiovascular changes [39]. NMU produces a brain-gut peptide with potent contractile effects that is significantly expressed in human plasma, epicardial adipose tissue, and the vasculature, including the left ventricle, coronary artery, and the saphenous vein. Several studies have shown that NMU precursor mRNA is significantly up-regulated in the left ventricle in dilated cardiomyopathy and ischemic heart disease patients [40, 41].

In conclusion, we combined the powerful LIMMA differential gene expression analysis with WGCNA to identify the hub genes and pathways related to AF perpetuation, recurrence, and maintenance. Of note, 152 DE genes and 7 hub genes, including LEP, FOS, EDN1, NMU, CALB2, TAC1, and PPBP were identified and may contribute to the pathogenic and pathophysiological mechanisms. Moreover, the pathway of ECM-receptor interaction was primarily enriched and closely linked to the initiation and maintenance of AF.

Although gene-based analysis results are encouraging, we currently do not understand how candidate genes contribute to AF or whether these pathways will translate into improved clinical outcomes to prevent AF. It is worth noting that while the hub genes and pathways in the comparison between AF ‘triggers’ and ‘substrates’ were found, these findings need to be verified in animal models. Moreover, the PPI network information, GO terms, and KEGG pathway maps from the STRING and DAVID databases are not available in similar databases.

This work was supported by the National Natural Science Foundation of China (NSFC No.81771165), the Special Program for Applied Research on Super Computation from the NSFC-Guangdong Joint Fund (the second phase), the Natural Science Foundation Project in Guangdong province, China (Grant No. 2016A030313295), the Major Project Development and Emerging, Interdisciplinary Funding Projects of Sun Yat-sen University (Grant No. 15ykjc17b) and the Guangzhou science and technology project of Major Special Research Topics on International Collaborative Innovation (Grant No. 201704030032).

The authors have no conflict of interests to disclose.

R. Zou, M. Yang and W. Shi contributed equally to this work.