Transcriptomic Profiling of Canine Atrial Fibrillation Models After One Week of Sustained Arrhythmia.

[Figure: see text].

Atrial fibrillation progression is associated with electrical and structural (fibrotic) remodeling of the atria.

Patients with atrial fibrillation show substantial gene expression changes in the atria, but most available data are from patients with advanced phases of the disease.

A portion of the atrial remodeling associated with atrial fibrillation is due to changes in ventricular rhythm and can be controlled in animal models by producing atrioventricular block and pacing the ventricles at a fixed, physiological rate.

After one week of tachypacing, extensive gene-regulation changes are identified by RNA-sequencing, with greater dysregulation in dogs with atrioventricular block compared to those without.

The imprinted DLK1-DIO3 locus is transcriptionally activated after one week of sustained arrhythmia.

Differentially expressed miRNAs of the DLK1-DIO3 locus identify target genes implicated in glutamate signaling, highlighting the potential importance of this system that has been shown to affect cardiomyocyte excitability.

Atrial fibrillation (AF) is the most common sustained arrhythmia, with an estimated lifetime risk of 22% to 26% and association with increased morbidity and mortality.[1] Despite advances in antiarrhythmic therapies, their suboptimal efficacy and adverse effects have limited their use.[2] Therefore, there is a need to further characterize fundamental arrhythmia mechanisms to discover new therapeutic targets.[2] Although AF is known to be a final common end point of atrial remodeling resulting from a variety of heart diseases, it can also be, in turn, a cause of remodeling. This vicious cycle is called AF begets AF[3] and explains the progressive nature of this arrhythmia and the complexity of its management.

Atrial remodeling is characterized by ion channel dysfunction, Ca2+ handling abnormalities, and structural changes, which result in AF induction and maintenance.[4,5] Heart disease, and even rapid atrial activity itself, cause the development of atrial fibrosis, which is a hallmark of structural remodeling. The degree of fibrosis is positively correlated with the persistence of AF.[6] Atrial cardiomyocytes subjected to rapid activation release factors that induce fibroblast-to-myofibroblast differentiation that leads to increased collagen synthesis.[7]

Any arrhythmia causing a rapid ventricular rate, including AF, is a well-recognized inducer of ventricular dysfunction, so-called arrhythmia-induced cardiomyopathy.[8] Heart failure enhances atrial stretch and sympathetic tone, making AF more resistant to rate- or rhythm-control treatments.[9] AF promotion results from the rapid atrial rate, but rapid ventricular rates due to inadequate rate control also promote AF-related atrial remodeling with a different profile from the remodeling produced by rapid atrial rate alone.[10] Radiofrequency atrioventricular node ablation (AVB) with right ventricular pacing is a nonpharmacological strategy for rate control that can improve symptoms and outcomes.[11]

Our previous work in canine AF models showed that maintaining AF for one week by rapid atrial pacing activates fibroblasts, collagen gene expression, and cardiomyocyte ion channel changes, without yet causing fibrosis.[12] Continued electrical maintenance of AF for 3 weeks produces fibrosis but electrically maintained AF with ventricular rate control through AVB produces less profibrillatory remodeling than 3 weeks of AF alone[10]. However, how AF with and without AVB impact atrial remodeling at the molecular level has not yet been assessed comprehensively. To answer this question, we took advantage of our well-characterized AF dog models and performed RNA-sequencing (RNA-seq) of cardiomyocyte-enriched atrial samples after one week to capture the molecular actors of atrial remodeling. In comparison with control (CTL) dogs, we found thousands of mRNAs, long noncoding RNA (lncRNA), and microRNA (miRNA) that are differentially expressed (DE) in the atria of the canine AF models. Pathway analyses of the transcriptomic data highlighted known biological processes but also potential novel modulators of arrhythmia initiation which may shed new light on our understanding of AF in humans.

Methods

All results and R code are available at https://github.com/lebf3/DogAF.

Canine AF Model

A total of 18 adult mongrel dogs of either sex, weighing 18 to 32 kg, were obtained from LAKA Inc and randomly assigned to CTL group (n=6) and 2 canine AF models (n=6/group; Table I in the Data Supplement). We elected to study an N of 6 animals per group based on previously published RNA-seq studies and our previous experiments with these models. Further, the number of DE genes identified in our analyses (post hoc) indicates that this sample size is sufficient to capture the main transcriptional changes that occur in the atrium of these dog models. Animals were handled in accordance with the Guide for the Care and Use of Laboratory Animals established by the National Institutes of Health as approved by the Montreal Heart Institute Ethics Committee (2016-47-01, 2019-47-03 for control dogs, 2015,47-01, 2018.47.12 for AF dogs).

To induce AF, animals were subjected to atrial tachypacing without (AF)[13] and with (AF+AVB)[14] atrioventricular node ablation under 0.07 mg/kg acepromazine (IM), 5.3 mg/kg ketamine (IV), and 0.25 mg/kg diazepam (IV), and 1.5% isoflurane anesthesia. In the AF group, a bipolar pacing lead was placed in the right atrial appendage under fluoroscopic guidance. In the AF+AVB group, pacing leads were inserted into the right atrial appendage and right ventricular apex. Pacing leads were connected to a subcutaneous pacemaker implanted in the neck (right side). In the AF+AVB group, radiofrequency catheter ablation was used to create AF+AVB. For this purpose, a quadripolar catheter with fluoroscopic guidance was placed across the tricuspid valve via the right femoral vein. Radiofrequency energy was then used to perform ablation when action potential at the His bundle was detected. Twenty-four to seventy-two hours after surgery, dogs in the AF group were subjected to AF-maintaining atrial tachypacing at 600 bpm for 7 days. In the AF+AVB group, RA and right ventricle were paced at 600 and 80 bpm, respectively. In animals of the CTL group, no pacemaker was inserted. No adverse event was recorded, and no dog was excluded.

Enrichment of Dog Atrial Cardiomyocytes

Cardiomyocytes were enriched from the left atrium with enzymatic digestion through the coronary artery-perfused Langendorff system, as previously described.[15] Briefly, dogs were anesthetized with 2 mg/kg morphine (IV) and 120 mg/kg alpha-chloralose and mechanically ventilated. Hearts were aseptically and quickly removed after intraatrial injection of 10 000 U heparin and placed in Tyrode solution containing 136 mmol/L NaCl, 5.4 mmol/L KCl, 2 mmol/L CaCl2, 1 mmol/L MgCl2, 10 mmol/L dextrose, 5 mmol/L HEPES, 0.33 NaH2PO4 (pH was adjusted to 7.3 with NaOH). The left coronary artery of the isolated heart was cannulated, and the left atrium was dissected free and perfused with 100% oxygenated Tyrode solution (37 °C, 1.8 mmol/L Ca2+). The arterial branches were ligated to have a leak-free system, and left atrium tissues were perfused with Ca2+-free Tyrode solution for ≈10 minutes, followed by ≈1-hour perfusion with ≈0.45 mg/mL CLSII (collagenase, Worthington, Lakewood, NJ) and 0.1% bovine serum albumin (Sigma–Aldrich, Oakville, ON) in Ca2+-free Tyrode solution for enzyme digestion. Digested tissue was removed from the cannula and cut into small pieces, and atrial cardiomyocytes were harvested.

RNA-Seq/miRNA-Seq

Library Preparation and Sequencing

mRNA and miRNA libraries were prepared at Genome Québec. mRNA libraries were made with the NEBNext_dual kit (rRNA-depleted stranded) and sequenced on NovaSeq 6000 S2 PE100 Illumina platform generating 32-123M Paired-end reads per sample. miRNA libraries were prepared with TruSeq smRNA and sequenced on the HiSeq 4000 SR50 Illumina platform generating 10 to 12 mol/L reads per sample.

Bioinformatic Processing and DE Analysis

The complete analysis can be found at https://github.com/lebf3/Dog_AF_transcriptomic. Briefly, mRNA reads were pseudomapped on reference transcriptome CanFam3.1.98 with Kallisto[16] with the options quant -t 5 -b 100 and the rest as default. We aggregated transcripts by genes with tximport[17] and quantified with DESeq2.[18] Genes with 0 reads in >12 samples were removed. Shrunken log2 transformed expression corrected for library size, with and without fibroblast fraction as a covariate (within DEseq2’s model) were then analyzed for DE with Wald test for all pairwise comparisons of CTL, AF, and AF+AVB and likelihood ratio test for total assay DE. We did not adjust our differential gene expression analyses for biological sex because no genes where DE between female and male dogs in our experiment. We plotted the principal component analysis with fibroblast fraction as a covariate from log2 transformed expression values corrected with Limma’s removeBatchEffect() function[19] for visualization of fibroblast effect on the top 1000 most variable genes. We then compared sets of GENEIDs found to be up (log-2 fold change >0 and P<0.01) or down (log-2 fold change <0 and P<0.01) in all possible contrasts.

For miRNAs, we trimmed reads using fastp[20] with default settings and aligned them to CanFam3.1.98 genome with STAR v2.7.1a[21] according to Encyclopedia of DNA Elements protocol.[22] DEseq2 DE analysis was then conducted with the same parameters as described above for mRNAs.

Deconvolution of RNA-Seq Data

To account for potential tissue heterogeneity, we used a murine atrial gene signature matrix described in Donovan et al[23] and our matrix of gene expression in Fragments Per Kilobase of exon model per Million reads mapped in CIBERSORTx online tool.[24] We then performed nonparametric Wilcoxon test on all possible comparisons for fibroblast fraction with a statistical significance threshold of P<0.05.

Gene Set Enrichment Analyses

For each gene sets described above, we performed hypergeometric testing against the human Gene Ontology Biological Processes from Molecular Signatures Database v7.1 with the HypeR package.

miRNA Target Prediction

For DE miRNA present in the 5 most cited miRNA databases (DIANA, Miranda, PicTar, TargetScan, and miRDB), we defined genes as targets if they were (1) annotated with a human homolog in the ensemble database, (2) predicted targets by at least 3 out of 5 databases queried with the MiRNAtap package, (3) DE (mRNA false discovery rate [FDR] <0.01), (4) inversely correlated (Pearson r<−0.5) log2 expression, corrected for the fibroblast fraction (expression values corrected with Limma’s removeBatchEffect() function). We then performed a gene set enrichment analyses (GSEA) with the remaining 82 predicted targets of the miRNA located on the syntenic region of the Dlk1-Dio3 locus (CanFam3.1 Chr8:68961744-69696779) as described above.

RNA-Seq and miRNA-Seq DE Genes Comparison Between Human AF Patients and Canine AF Models

DE genes in our canine AF models with annotated human orthologues in the ENSEMBL database were compared with a meta-analysis of miRNA DE in human AF and a large RNAseq study on left atrial appendages obtained from 261 patients undergoing valve surgery.[25,26] Precursors of the human miRNAs listed in Shen et al Table S8 (n=53, 21 upregulated, and 32 downregulated)[25] were retrieved with the R package miRBaseConverter and then compared across species. Because only 5 miRNA were found to overlap with human DE miRNA, only mRNA genes found to be DE in human and our canine AF models are represented as an Upset plot. The counts for human mRNA data were downloaded from GEO database (GSE69890). DE testing was conducted as described above for the 3 groups; no AF (CTL, n=50), AF in AF rhythm (AF, n=130), and AF in sinus rhythm (AF.SR, n=81), with inclusion of sex as a covariate.

Mitochondrial Genes DE in Canine AF Models

The human MitoCarta3.0[27] database was queried for genes with mitochondrial localization in the heart (n=539). We represented DE genes with likelihood ratio test FDR <0.01 from that list as volcano plots.

Proteomics

Dog cardiomyocytes were lysed by sonication, reduced, and alkylated. Protein was precipitated, resuspended, quantified, and subjected to tryptic digest. Peptides (500 ng) were analyzed by reverse phase nano-HPLC coupled to a Bruker maXis II mass spectrometer (positive mode, mass range 150–2200 m/z, collision induced dissociation of top 20 precursors). LC-mass spectrometry/mass spectrometry data were analyzed for protein identification and label-free quantification using MaxQuant[28] (1.6.1.0) against the public database UniProt with taxonomy Canis lupus familiaris and common contaminants (downloaded on 01.08.2019, 29809 sequences) with carbamidomethylation on Cys as fixed and oxidation on Met as variable modification with decoy database search included (mass tolerance 0.006 Da for precursor, 80 ppm for product ions; 1% PSM and protein FDR, match between runs enabled, minimum of 2 ratio counts of quantified razor and unique peptides).

DE Analysis and Correlation

Proteins with >3 missing values per treatment were removed. The remaining missing intensities were replaced with random values taken from the Gaussian distribution centered around a minimal value from the 10th quantile with the DEP package’s Minprob function, to simulate a relative label-free quantification value for those low abundant proteins. Two-sample t tests with subsequent multiple testing correction by FDR were used to identify DE proteins (P<0.01) with the fibroblast fraction as covariate using the Limma package.

Because proteomic processing does not always converge to a single protein, only 755 genes out of the 1029 in the proteomic matrix were correlated to their corresponding RNA-seq data. We compared overlapping genes’ mean log2 transformed expression in proteomic and RNA-seq. The distribution of mean RNA-seq expression of the 755 overlapping genes was then compared with the full mean RNA-seq gene expression values.

Results

RNA-Sequencing of Cardiomyocyte-Enriched Atrial Samples From Canine AF Models

We analyzed data from 3 groups of 6 dogs. The first group (CTL) was the control group without atrioventricular ablation (AVB) nor pacemaker, in the second group (AF), right atrial tachypacing at 600 beats per minute was used to maintain AF electrically for one week, and the third group (AF+AVB) included dogs with electrically maintained AF for one week in the presence of AVB and ventricular pacing at 80 beats per minute to control the ventricular rate. We reasoned that transcriptomic profiling of atria from these animals should allow us to discover the molecular changes that occur over the first week after the onset of AF and play a role in the development of the tissue remodeling accompanying the transition from paroxysmal to persistent AF.

Initial analysis of bulk RNA-seq data hinted at some heterogeneity of cellular composition across samples. Therefore, we estimated the fraction of the major cell types in each sample using an in silico deconvolution technique implemented in CIBERSORTx (Figure 1A).[24] Because of the induced tissue remodeling due to the AF treatments, we found that both AF and AF+AVB dogs had more fibroblasts in their atria than CTL animals (Figure 1B). To emphasize the transcriptional differences between conditions that are not a result of variable cellular composition, we included the fibroblast fraction as a covariate in all subsequent DE analyses. Correction for this confounding variable reduced intergroup variability (Figure 1C and 1D).

Figure 1.

Deconvolution of canine atria cell composition using bulk RNA-sequencing.A, We inferred cell fractions with CIBERSORTx and an atrial-specific gene signature matrix obtained using orthologous murine genes.[23] We present cell fractions for each dog sample that we analyzed in this study. B, When we group animals per treatment arm, we observed a significantly higher fraction of fibroblasts in the atrial fibrillation dog models (AF and AF+AVB) than in the control animals (AFvsCTL Wilcoxon test, P=0.0087 and AF+AVBvsCTL, P=0.015). Principal component analysis of the top 1000 most variable genes expressed in canine atria before (C) and after (D) correction for fibroblast fraction show treatment-dependent clustering after correction for cell composition. AF indicates Atrial fibrillation; AF+AVB, AF with atrioventricular block; and CTL, control.

Proteomic Analysis Largely Confirms the Transcriptomic Results

To validate our RNA-seq results, we took advantage of mass spectrometry-based protein quantification results from the same 18 dog atrial cardiomyocyte-enriched cell extractions that were generated in a parallel study (detailed proteomic results will be presented elsewhere). After stringent quality control, we obtained relative quantification for 755 proteins. For these genes, the relative RNA and protein levels were strongly correlated (Pearson r=0.49, P=1.57×10−46; Figure 2A). Many of the genes that are well-correlated encode abundant cardiomyocyte proteins, such as titin (TTN), myosin light chain-4 (MYL4), desmin (DES), and tropomyosin-1 (TPM1). We found that RNA-seq could profile transcripts with a wider range of expression profiles, whereas mass spectrometry-based proteomics preferentially captured proteins whose genes are expressed at high levels (Figure 2B).

Figure 2.

Validation of highly expressed RNA by proteomics.A, In 18 atrial samples, 755 genes (NS=619, RNA.sign=122, both.sign=14) found in both data sets are highly correlated at the protein (x axis) and RNA (y axis) levels (Pearson, r=0.49; P=1.57×10−46). For reference, we annotated 15 genes that are differentially expressed in the RNA-seq experiment and have high protein expression levels. NS, not differentially expressed in the RNA-seq or proteomic experiment; RNA.sign, genes that are differentially expressed in the RNA-seq assay only; Both.sign, differentially expressed genes in both the RNA-seq and proteomic experiments. DE genes in the RNAseq data set have an false discovery rate (FDR) <0.01 (likelihood ratio test) and proteomics data set an FDR <0.05 (F test). The gray area around the line corresponds to the 95% CI. B, Relative expression level of all transcripts measured in the RNA-seq experiment. The histogram shows that genes that are present in both the RNA-seq and proteomic experiment are highly expressed (Common, dark gray) in comparison to the expression levels of all transcripts measured (all_RNA, light gray).

Transcriptomic Changes in Cardiomyocyte-Enriched Atrial Samples

Pairwise comparisons of gene expression levels between the three groups of dogs identified 434, 5971, and 7867 genes that are DE (FDR <0.01) in atrial cardiomyocyte-rich fractions in atrial fibrillation versus control samples (AFvsCTL), atrial fibrillation+AV block versus control samples (AF+AVBvsCTL), and AFvsAF+AVB, respectively (Figure 3A and 3B). All differential gene expression level results are available in Table II in the Data Supplement and https://github.com/lebf3/DogAF). Many genes previously implicated in AF are dysregulated in both AF and AF+AVB dogs when compared with controls, thus validating the experimental design. This includes FHL1 involved in myofilament regulation,[29] SORBS2 involved in intercalated disc gap junction regulation,[30] and KCNA5, which regulates atrial action potential repolarization.[31] Previous studies have established an important role for mitochondrial dysfunction in the cause of AF.[27] Accordingly, we identified 54 genes that encode mitochondrial proteins that are DE in our AF canine models (Figure I in the Data Supplement). In particular in the AF+AVB group, we noted the upregulation of 2 key beta-oxidation genes (CPT1A and ACADL) and the downregulation of the electron transport chain genes COX17 and NDUFA8. However, our data also implicates genes not previously recognized to be involved in AF, such as leukocyte receptor cluster member-8 (LENG8), transcription elongation regulator-1 (TCERG1), ligand dependent nuclear receptor corepressor (LCOR), formin-binding protein-4 (FNBP4), and ENSCAFG00000049959 (orthologue of the lncRNA MEG3; Figure 3A, Table III in the Data Supplement, and https://github.com/lebf3/DogAF).

Figure 3.

Analyses of differentially expressed atrial genes identify many biological pathways that are dysregulated in atrial fibrillation (AF) dog models.A, Volcano plots of all transcripts that we analyzed in this study. Transcripts in red have a false discovery rate (FDR) <0.01. We found 434, 5971, and 7867 genes that were differential expression (DE) in the AFvsCTL, AF+AVBvsCTL, and AFvsAF+AVB analyses, respectively. The full DE results are available in Table II in the Data Supplement. B, Upset plot showing the intersection of upregulated and downregulated DE genes (FDR <0.01) in each analysis. C, The 5 most significant biological pathways identified using gene set enrichment analyses (GSEA) for each set of DE genes (FDR <0.01). Full results are available in Table IV in the Data Supplement. AF+AVBvsCTL indicates atrial fibrillation+AV block versus control samples; AFvsCTL, atrial fibrillation versus control samples; and AVB, atrioventricular block.

To understand what pathways are modulated in the atria of these canine AF models, we performed GSEAs on the DE genes (Figure 3C and Table IV in the Data Supplement). In AFvsCTL, we noted an upregulation of genes associated with profibrotic pathways (eg, extracellular structure organization, biological adhesion, response to wounding) and a downregulation of genes implicated in angiogenesis, such as blood vessel morphogenesis. Genes implicated in muscle biology were upregulated in the AF+AVBvsCTL analysis (eg, muscle structure development, striated muscle cell differentiation) whereas the same comparison implicated downregulated genes involved in ion transport and signaling pathways (eg, sensory perception). We confirmed that this enrichment was not due to a smaller fraction of cardiac neurons found in the atria of AF+AVB dogs (Kruskal-Wallis, P=0.32). Because of the large overlap in genes that are downregulated in AF+AVBvsCTL and upregulated in AFvsAF+AVB (Figure 3B), we identified similar pathways in the GSEA for these 2 comparisons (in Figure 3C, compare AF+AVBvsCTL and AFvsAF+AVB). Finally, genes that were downregulated in the AFvsAF+AVB analysis implicated genes with more generic functions in gene expression and chromatin modifications, such as the histone-lysine N-methyltransferase SETD5 and the DNA methyltransferase TET2. Dysregulation of the expression of these chromatin-related genes and pathways is consistent with the extensive transcriptomic changes observed in the atria of AF+AVB dogs, in sheep models of AF (Figure II in the Data Supplement) as well as in patients with AF.[32,33]

Dysregulation of miRNA Expression

Because miRNA play important roles in AF biology[34] but are not detected in standard RNA-seq protocols, we performed in parallel miRNA-seq on the same dog samples. We found 31, 19, and 21 miRNA that are DE (FDR <0.01) in AFvsCTL, AF+AVBvsCTL, and AFvsAF+AVB, respectively (Figure 4A, Table II in the Data Supplement). When comparing miRNA expression in the 2 AF models, MIR185 on the dog chromosome 26 was the most DE miRNA with strong upregulation in the atria of AF animals. We also noted that 11 of the most strongly DE miRNA in the AF+AVBvsCTL and AFvsCTL analyses (MIR136, MIR411, MIR370, MIR127, MIR493, MIR494, MIR485, ENSCAFG00000025655 [96.20% identity to hsa-mir-379], MIR758, MIR543, MIR889) mapped to the chr8:68 900 000 to 69 700 000 region in the dog reference genome CanFam3.1 (Figure 4B). This region, highly conserved in mammals, is syntenic to the imprinted 14q32 region in humans (also known as the DLK1-DIO3 locus).[35] The lncRNA MEG3, which we described above as being over-expressed in the AF canine models is also located in the same DLK1-DIO3 syntenic dog locus.

Figure 4.

Eleven differentially expressed microRNAA, Volcano plots of all miRNA that we measured in our experiments. We identified 31, 19, and 20 miRNA that are differentially expressed (false discovery rate [FDR] <0.01) in the AFvsCTL, AF+AVBvsCTL and AFvsAF+AVB analyses, respectively. B, Miami plots of miRNA and their corresponding statistical significance (y axis) for the AF+AVBvsCTL (top) and ATvsCTL (bottom) analyses. An arrow indicates the miRNA cluster located on the canine chromosome 8 region that is syntenic to human DLK1-DIO3. The odd and even chromosomes FDR values are in blue and red, respectively. C, Upset plot showing the DE miRNA targets located in the syntenic DLK1-DIO3locus and their corresponding number of potential target RNA. We identified potential targets with the MiRNAtap package (predicted by ≥3 databases) from DE miRNA (FDR <0.01) and DE mRNA (FDR <0.01). D, Gene set enrichment analyses (GSEA) with the potential gene targets (x axis) of the DE miRNA located at the syntenic DLK1-DIO3 locus. We only present the top 5 pathways enriched in this analysis. A red square in the heatmap indicates membership of a given target gene to the biological pathways located on the left (empty columns were removed for clarity). GSEA FDR and AF+AVBvsCTL DE FDR are on the right and top of the heatmap, respectively. AF indicates atrial fibrillation; AF+AVBvsCTL, atrial fibrillation+AV block versus control samples; AFvsCTL, atrial fibrillation versus control samples.

The dysregulation of the expression of lncRNA and miRNA at the same locus suggested that they might co-regulate the expression of genes implicated in the same biological pathway(s). To address this possibility, we used in silico predictions to infer the DE mRNA that are possible direct targets of these DE miRNA located at the syntenic DLK1-DIO3 locus. For this analysis, we focused on DE miRNAs and DE mRNAs that are predicted to physically interact by at least 3 out of 5 databases and that have expression levels that are negatively correlated in the RNA-seq/miRNA-seq experiments (Pearson r<−0.5). Using these filters, we identified 82 potential target genes for the DE miRNAs at this locus, with most genes targeted by a single miRNA (Figure 4C). GSEA with these 82 genes indicated a common role in synaptic signaling involving glutamate signaling (Figure 4D and Table V in the Data Supplement). Some of the key genes within these pathways are metabotropic glutamate receptor-1 and -8 (GRM1, GRM8), glutamate ionotropic receptor delta type subunit-1 (GRID1), glutamate ionotropic receptor AMPA (α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid) type subunit-1 (GRIA1), and corticotropin-releasing factor-binding protein (CRHBP).

Partial Differential Transcriptomic Overlap Between Human and Dog AF Atrial Samples

To assess the ability of our canine AF models to capture early transcriptomic changes which might be missed by profiling the atria of human AF patients who have developed the disease over years, we compared the DE genes identified in AF dogs with results from human left atrial transcriptomic profiling experiments.[25,26] Hsu et al performed bulk RNA-seq experiments on left atrial appendages from 261 patients who underwent cardiac surgery to treat AF, valve disease, or other cardiac disorders. For differential gene expression analyses, these patients were divided between no AF (n=50), AF in AF rhythm (n=130), and AF in sinus rhythm (n=81). When we intersected this list of human DE genes with the list of DE genes in our AF dog models (with clear human orthologs), we identified 668 genes (Figure 5 and Table VI in the Data Supplement). We found the strongest overlap between dog AF+AVB and human AF in AF rhythm. Of note, most of the strongest signals in our dog study are also present in this human study (LENG8, SORBS2, BMP10, FNBP4, and glutamate receptor-related genes [GRM1, GRM8, GRIA1, GRIK2, GRID1]). For miRNA, we compared our results with data from a large meta-analysis involving 40 articles and 283 DE miRNA in AF (in different tissues and species).[25] Of the 53 AF-associated miRNAs that were previously identified in human cardiac tissues and showed consistent results in the meta-analysis, we found 5 miRNAs in our analyses of the dog transcriptomic data sets (MIR144, MIR142, MIR146B, MIR223, and MIR451). These miRNAs have not been characterized functionally yet for a role in AF. Generally, the canine AF group matched better the directionality of change reported in the human AF DE miRNA meta-analysis.

Figure 5.

Overlaps in genes differentially expressed in canine AF models and human AF patients. We compared differentially expressed genes in our canine AF models with annotated human orthologues that are DE in human AF left atrial appendages.[26] Homo Sapiens; hsa, hsa_AF; AF in AF rhythm, hsa_AF.SR; AF in sinus rhythm, hsa_CTL; no AF.

Discussion

In this study, we used a transcriptomic approach to comprehensively assess the molecular architecture of AF-induced remodeling with and without AVB from atrial cardiomyocyte-enriched samples. We validated the robustness of our RNAseq data by correlating it with a proteomic analysis, which showed a strong correlation in gene expression, with well-defined cardiomyocyte genes being most highly expressed in both datasets (eg, TTN, MYL4). We confirmed the involvement of known AF factors like the reactivation of developmental pathways, but also found a strong and novel association with microRNAs and lncRNA from the DLK1-DIO3 locus, including the MEG3 canine orthologue. This finding is concordant with the many chromatin remodeling genes dysregulated in our models, which is an emerging phenotype of AF both in human and sheep models.[36]

Molecular Remodeling in AF With Versus Without AVB

We did not expect to find a smaller number of DE genes in the AFvsCTL analysis than in the AF+AVBvsCTL analysis, given our previous observation that AF treatment alone without AVB results in more important tissue remodeling.[10] One possible explanation is that our prior histological studies were done in AF animals treated for 3 weeks,[10] whereas the results presented here reflect RNA changes after 1 week of AF. The transcriptomic changes in the AF+AVB group show that cells are under active chromatin modification, indicating ongoing adaptation to the stimulus. This is not observed in the AF group (lacking AVB), which may indicate that this adaptation has already occurred. This idea would be consistent with the downregulation of chromatin-related genes recently noted in the atria of sheep AF models[36] and may be a result of earlier establishment of profibrotic transdifferentiation in AF compared with AF+AVB canine models.

Our analyses highlighted many genes not previously implicated in AF. While the functions of some of these genes remain uncharacterized (eg, LENG8 in AF+AVB), we can speculate on the activities of others. For instance, TCERG1 and FNBP4, which are upregulated in the atria of AF+AVB dogs, encode coregulated proteins that are involved, respectively, in RNA splicing and translation.[37] The upregulation of these genes, with general actions on gene transcripts, may (partly) explain why more genes are dysregulated in AF+AVB animals when compared with the CTL or AF groups (Figure 3B). Another interesting candidate is LCOR, which is upregulated in the AF+AVB model and encodes a transcriptional cofactor that interacts with PPARγ and RXRα to control gene expression.[38] While further experiments are needed to determine the extent by which LCOR modifies gene expression in AF+AVB and contributes to the pathology, our RNA-seq experiments detected the dysregulation of 2 of its likely targets based on the literature: CPT1A (see above) and the cell cycle regulator CDKN1A.[39]

Potential Role of Noncoding Genes at the DLK1-DIO3 Locus in Early AF

The highly conserved DLK1-DIO3 locus hosts 2 differentially DNA-methylated regions modulating the expression of its noncoding RNA clusters, where in humans the maternal allele is hypomethylated with concomitant expression on the hypermethylated paternal allele of noncoding RNA and other protein-coding genes (DLK1, RLT1, and DIO3).[35] In both AF+AVBvsCTL and AFvsCTL, we found a large proportion (58% and 23%, respectively) of DE miRNA at this locus, underlying its importance in AF-related adaptation. We also found dysregulation of the MEG3 lncRNA canine orthologue at this locus. MEG3 is a highly expressed lncRNA that has been studied in various pathologies, including cancer[40] and more recently cardiovascular diseases.[41] Noncoding RNAs at this locus have been shown to mediate various cardiac developmental programs.[35] More specifically, MEG3 can contribute to the recruitment of the polycomb repressive complex-2 (PRC2),[42] a key chromatin modulating factor. Of particular interest, Mondal et al[43] showed that through interaction with the H3-Lys-27 methyltransferase EZH2, MEG3 can repress TGF-beta target genes, which are known to promote a profibrotic response. Data have been presented that suggest an important role of EZH2 and EZH2-regulated genes in AF.[44]

Glutamate Receptor Regulation by miRNAs From the DLK1-DIO3 Locus

Our GSEA analysis-predicted gene targets of DE miRNA at the DLK1-DIO3 locus suggest a role for glutamate signaling in AF. Immunostaining has confirmed the presence of glutamate receptors on cardiomyocytes.[45] Glutamate was also found to be significantly increased in AF patient left atrial appendages.[46] Glutamate signaling is important in vagal afferent neurons,[47] and remodeling of the glutamate system in AF may relate to the extensive previous evidence of autonomic dysfunction in patients with AF.[48] Moreover, a recent study has shown fundamental roles for glutamatergic receptors in rat atrial cardiomyocytes and induced pluripotent stem cell-derived atrial cardiomyocytes, including a reduction in cardiomyocyte excitability after GRIA3 knockdown.[49] Therefore, the DLK1-DIO3 miRNA cluster may be an adaptative regulator of cardiomyocyte excitability or of neural cells in the presence of AF.

Limitations

We used cardiomyocyte-enriched samples in an attempt to obtain clearer results from the transcriptomic analysis by excluding extrinsic variability due to changes in cell composition. However, while our samples are enriched in cardiomyocytes, they do not constitute a pure cardiomyocyte preparation. A disadvantage is that variability due to changes in cell composition is not eliminated. On the other hand, our cardiomyocyte-enriched (but not pure) samples allow us to detect potential features of AF related to noncardiomyocyte cells, such as autonomic dysregulation mediated by neural cells; however, we cannot unambiguously attribute DE genes to transcriptomic changes in a specific cell type. In part, we were able to control for fibroblast composition by adjustment through analysis for expression of fibroblast-related RNA-expression patterns. Nevertheless, features underlined here should be confirmed in pure cell cultures or single cell transcriptomic assays. A second limitation is the difficulty in extrapolating our findings to gene expression changes in humans. We found only modest overlap of DE genes in our model compared with reported DE gene patterns in human; several factors could explain this (eg, differences in biospecimen preparation, tissue heterogeneity, fundamental differences between dog and human AF pathology). It is also possible that different transcriptomic programs may be involved at the initiation of arrhythmia and tissue remodeling (AF and AF+AVB dog models) when compared with those dysregulated in the atria once the pathology has been present for years.

We compared our transcriptomic results with proteomic data on the same samples and found a very high level of correlation (Figure 2). These results are consistent with a large measure of transcriptomic control over protein expression and validate the relevance of transcriptomic analysis of these data. A further in-depth look at the proteomic signature in these models would be of interest but is beyond the scope of the present article.

Conclusions

Understanding the pathophysiology of chronic human diseases such as AF is challenging because they develop over many years and initially present with only unremarkable preclinical symptoms. In this study, we took advantage of 2 well-characterized canine AF models to chart the transcriptomic changes that occur at the earlier phases of arrhythmia. Despite the inherent limitations in relating dog models to human AF, our results offer interesting new hypotheses for future testing, including in man. In particular, the upregulation of miRNAs at the DLK1-DIO3 locus after 1 week of AF suggests that they may be early biomarkers of tissue remodeling and adaptation in the atria.

Sources of Funding

This work was funded by the Fonds de Recherche en Santé du Québec (FRQS), the Canada Research Chair Program, the Montreal Heart Institute Foundation (MHIF), Heart and Stroke Foundation of Canada (grant No. 18-0022032), Canadian Institutes of Health Research (CIHR; grant No. 148401), the Austrian Science Fund (FWF; projects KLI645, W1226, and F73 to RBG). F.J.A. Leblanc received scholarships from the CIHR, FRQS, MHIF, and Université de Montréal.

Disclosures

None.

Supplemental Materials

Supplemental Table I–VI

Supplemental Figures I–II