The arrhythmogenic cardiomyopathy-specific coding and non-coding transcriptome in human cardiac stromal cells
BMC Genomics volume 19, Article number: 491 (2018)
Arrhythmogenic cardiomyopathy (ACM) is a genetic autosomal disease characterized by abnormal cell-cell adhesion, cardiomyocyte death, progressive fibro-adipose replacement of the myocardium, arrhythmias and sudden death. Several different cell types contribute to the pathogenesis of ACM, including, as recently described, cardiac stromal cells (CStCs). In the present study, we aim to identify ACM-specific expression profiles of human CStCs derived from endomyocardial biopsies of ACM patients and healthy individuals employing TaqMan Low Density Arrays for miRNA expression profiling, and high throughput sequencing for gene expression quantification.
We identified 3 miRNAs and 272 genes as significantly differentially expressed at a 5% false discovery rate. Both the differentially expressed genes as well as the target genes of the ACM-specific miRNAs were found to be enriched in cell adhesion-related biological processes. Functional similarity and protein interaction-based network analyses performed on the identified deregulated genes, miRNA targets and known ACM-causative genes revealed clusters of highly related genes involved in cell adhesion, extracellular matrix organization, lipid transport and ephrin receptor signaling.
We determined for the first time the coding and non-coding transcriptome characteristic of ACM cardiac stromal cells, finding evidence for a potential contribution of miRNAs, specifically miR-29b-3p, to ACM pathogenesis or phenotype maintenance.
Arrhythmogenic cardiomyopathy (ACM) is a genetic autosomal disease that affects about 1 in 2000–5000 young people, in particular athletes [1, 2]. It is mainly caused by mutations in genes encoding desmosomal proteins even if other genes have been associated with ACM too . ACM is characterized by cardiomyocyte death and progressive accumulation of fibro-adipose tissue, but the pathogenic mechanisms remain still largely unknown. Abnormal cell-cell adhesion and intracellular signaling, caused by the mutations, are thought to alter Wnt/β-catenin  and Hippo  pathways resulting in myocyte death, fibro-adipogenesis, gap junction and ion channel remodeling. This leads to myocardium wall thinning, progressive heart failure, inefficient current conduction, and life threatening arrhythmias . For a long time, ACM has been considered a cardiomyocyte-specific disease, but this concept has been questioned when other cell types were found to be implicated in ACM pathogenesis. In particular, cardiac stromal cells (CStCs), an abundant cardiac cell population with differentiation potency  and a role in cardiac homeostasis and remodeling during pathological conditions , have recently been characterized as a source of adipocytes in ACM . Human primary CStCs from ACM patients express desmosomal proteins and show, compared to healthy controls, a greater propensity to undergo adipogenic differentiation in vitro in a PKP2- and Wnt pathway-dependent fashion . Importantly, CStCs can contribute to the formation of adipose tissue in the ACM affected heart.
While extensive characterizations have been performed on the genetic basis of this disease, only few studies have been conducted investigating differences at transcriptional levels [9,10,11]. The present work addresses, for the first time, the characterization of the ACM-specific coding and non-coding transcriptome in CStCs. To this end we compared expression of genes and miRNAs in cells derived from diagnostic biopsies of ACM patients or controls. In order to reveal biological pathways potentially involved in ACM we performed category enrichment and network analyses on the differentially expressed genes and validated target genes of the deregulated miRNAs.
Biological samples and cell culture
Biopsy samples were obtained from ACM patients’ right ventricle during diagnostic procedures in the area adjacent to the electroanatomical scar. Control samples were obtained from the right ventricle of explanted hearts of cadaveric donors (healthy hearts of organ donors). To limit the variability due to regional differences, we selected small portions of the endocardium in the right ventricle free wall of the donors’ hearts, as this matches the most likely source of ventricular biopsies in ACM patients. All enrolled individuals are males aged between 37 and 64 (mean age for ACM patients: 47.6, standard deviation 8.6, mean age for controls: 48.2, standard deviation 7.8). All ACM patients fulfill definite ACM diagnosis as defined by the Task Force 2010 criteria (either two major criteria or one major and two minor or four minor). PKP2 mutations in ACM patients were determined by Sanger sequencing. Samples, both from ACM patients and control individuals, were processed, as previously described, to obtain CStCs  on which the expression profiling was performed. Adherent cells were cultured, propagated and used between passage 3 and 9. As previously described , mesenchymal lineage was confirmed and endothelial and hematopoietic origin was excluded through flow cytometry analysis of surface markers. See Table 1 for an overview of samples, clinical data and experiments in which they were included.
RNA was TRIZOL-extracted and quantified with a NanoDrop ND-1000 Spectrophotometer. RNA quality was assessed by an Agilent 2100 Bioanalyzer (Small RNA Kit) for miRNA expression and by an Experion electrophoresis system (RNA High Sense Analysis Kit, Bio-Rad) for gene expression profiling. Only high quality RNAs (A260/280 and A260/230 ratios > 1.8 and an RQI > 9.5) were used for subsequent analyses.
miRNA expression analysis
miRNA expression was measured by TaqMan Low Density Arrays (TLDA, Human MicroRNA Panel A and B Card Sets v3.0, Applied Biosystems) on 480 ng of total RNA from 3 ACM patients and 3 healthy subjects using a 7900TH System (Applied Biosystems) following manufacturer’s instructions. The complete data set consisted of 24 TLDAs: 2 replicates for each sample and two cards, A and B, for each. Primary data analysis was conducted with SDS RQ manager (Applied Biosystems) applying a manual threshold (0.2) and baseline (based on cycles 10–16). Further analysis was performed in R (version 3.4.2). The data was normalized with the global geometric mean scaling method . miRNA assays were annotated to miRBase v20 using the LifeTech online search tool. Replicated measurements per individual were averaged prior to differential expression analysis and all assays that did not yield a valid measurement in at least 2 out of the 3 individuals in one of the two sample groups (because their measurements were flagged by the SDS software as “Undetermined” or “No amplification in well”) were excluded. This reduced the data set from 768 assayed to 208 detected miRNAs. Significance of differential expression was assessed with the moderated t-test , p-values were adjusted for multiple hypothesis testing with the method from Benjamini and Hochberg controlling the false discovery rate (FDR) . miRNAs with an adjusted p-value padj < 0.05 were considered significant.
High throughput sequencing and data analysis
Paired-end mRNA-sequencing libraries were prepared using the Illumina TruSeq RNA Sample Prep Kit v2-Set A starting from 1 μg of total RNA obtained from CStCs of 3 ACM patients and 3 healthy controls, following manufacturer’s instructions. Samples were multiplexed and sequenced on 4 lanes of an Illumina HiSeq2000 resulting in between 53 and 85 million read pairs per sample and a high average base call quality (Qphred score > 35). Subsequent analysis was performed in R (https://www.r-project.org, version 3.4.2, Bioconductor version 3.5). The 100 nt-long raw reads were filtered removing all reads with 7 or more undetermined nucleotides and subsequently trimmed as soon as 3 of 5 consecutive nucleotides in a read’s sequence had a base-call quality Qphred < 20. Finally, all reads shorter than 41 nt were removed. Alignment of the reads was performed against genome GRCh38 (Ensembl 81) using GSNAP .
Gene quantification was performed using GenomicAlignments  counting the number of reads falling completely within the exons of a gene, not considering multi-mapping reads or reads with a low alignment quality (phred score < 30). Genes with a normalized read count < 10 in both sample groups were removed, reducing the data set from 60,411 human genes defined in Ensembl to 18,076. Differential gene expression analysis was performed using DESeq2  disabling the automatic pre-filtering. Raw p-values were adjusted for multiple hypothesis testing using the method from Benjamini and Hochberg . All genes with an adjusted p-value padj < 0.05 were considered significant.
Quantitative reverse transcriptase–polymerase chain reaction
One μg of total RNA was reverse transcribed using the SuperScript® VILO™ cDNA Synthesis Kit (Invitrogen), according to manufacturer’s instructions. cDNA was amplified using TaqMan® Universal PCR Master Mix (Applied Biosystems). For evaluation of JUP expression, cDNA was amplified by SYBR-GREEN qPCR. miRNA expression was measured with TaqMan assays for the individual miRNAs. Real-time reactions were performed following manufacturer’s instructions on a CFX96™ RT-qPCR Detection System (Bio-Rad, USA). See extended Material and Methods in the Supplement for the list of primers and assays.
RT-qPCR data analysis was performed in R. Two technical replicates per sample were averaged and normalized against house-keeping control TBP (for genes) or miR-28-3p (for miRNAs) to yield the ΔCq values . miR-28-3p was selected due to its stable expression across all samples in the whole genome miRNA expression analysis. Where present, measurements in RNA from the same individual but from different experiments/passages were averaged prior to differential expression analysis. Significance of differential expression was assessed using the Student’s t-test.
Category enrichment analysis
The GOstats package  was used for category enrichment analyses, with categories representing either sets of genes annotated to pathways in the Reactome database, or to terms from one of the 3 Gene Ontologies Biological Process, Cellular Component or Molecular Function. Genes considered to be expressed in the present data set were used as background gene set. The p-values from the hypergeometric tests were adjusted for multiple hypothesis testing using the Bonferroni method and pathways/GO terms with an adjusted p-value < 0.05 were considered significant.
miRNA target gene and host gene analysis
Validated miRNA target gene information was taken from the miRTarbase (version 6.1)  and the analysis conducted using our “mirtarbase” R-package (https://github.com/jotsetung/mirtarbase). For identification of potential miRNA host genes/primary transcription units  we used our “mirhostgenes” R-package (https://github.com/jotsetung/mirhostgenes) with host gene predictions for miRBase v21 and genes defined in Ensembl version 81.
Systematic evaluation of the downregulation of ACM-miRNA target genes was performed as follows: for each miRNA, the average M-value (representing the differential expression of the gene) was calculated across all of its target genes and compared to average M-values of randomly selected, same-sized, sets of miRNA target genes in a permutation test (1000 permutations). The p-value represents the fraction of average M-values of the random gene sets that were < = the average M-value of the miRNA’s target genes (one-tailed test). Background genes were selected randomly among validated (and expressed) target genes of expressed miRNAs. We additionally used a Student’s t-test to test the alternative hypothesis of an average M-value < 0.
Gene network analysis
Gene interaction networks were defined using information from BioPlex  (version 4), mentha  (human data release 2015.11.27), STRING  (version 10, requiring confidence score of at least 0.9), and a functional similarity network  based on Gene Ontology (GO) biological process annotations (GO version 2015.09). Only genes that yielded a read count in the present transcriptome data set or that are expressed in the heart tissue  were included in the networks. Clusters within each network were identified using ClusterOne  with default parameters, requiring at least four genes per clister in the BioPlex/mentha/STRING network and five genes in the functional similarity network. Full network analysis results are presented on https://gemex.eurac.edu/bioinf/acm/. Network visualizations were created using python-igraph  and Cytoscape.js  and Cytoscape, GO annotations and enrichments (FDR < 5%, backgrounds comprising all genes in the respective networks) were computed using Dintor .
To delineate the miRNA expression profile characteristic for ACM cells, we extracted RNA from CStCs of 3 ACM patients and 3 (sex and age matched) control individuals and subjected it to whole genome miRNA expression profiling. We identified the 3 miRNAs miR-520c-3p, miR-29b-3p and miR-1183 to be significantly differentially expressed at a 5% FDR (i.e. with a padj < 0.05; Fig. 1b), all of them being more than 4-fold higher expressed in ACM than in control samples (Table 2; the complete results are provided in Additional file 1). For validation, we performed RT-qPCR measurements of miR-520c-3p and miR-29b-3p in RNA samples from the same 3 ACM and 3 control individuals, but from different passages in two separate experiments, and, to further strengthen the biological significance of the results, also in cells from 5 additional ACM patients and 3 additional controls. Validation was thus performed in RNA samples from 8 ACM patients and 5 control individuals (see Table 1 for clinical information of the individuals).
Differential expression was confirmed for miR-29b-3p but failed for miR-520c-3p, whose expression was below the detection levels in all tested samples (Table 2).
To identify biological pathways being potentially influenced by the differential expression of the identified miRNAs, we next performed Gene Ontology and Reactome pathway enrichment analyses on their target genes. We defined target genes using miRTarBase  and, by restricting to only high confidence target genes, identified 84 genes (77 target genes for miR-29b-3p and 7 for miR-520c-3p, no target gene was reported for miR-1183). These were tested for enrichment in Reactome pathways against a background gene set consisting of the 1486 validated target genes of all 208 expressed miRNAs and pathways related to extracellular matrix organization, cell-matrix adhesion and collagen formation were identified (Table 3). The results from the GO enrichment analysis for biological process, cellular component and molecular function supported this finding (Tables S1-S3 in Additional file 2).
Since these results based mainly on target genes of miR-29b-3p, we repeated the analysis considering also target genes with lower supporting evidence and thus increased the number of target genes from the other two miRNAs in the analysis (Additional file 2: Table S4). The results from this analysis were comparable (Additional file 2: Tables S6-S9).
Summarizing, we identified 3 miRNAs (miR-520c-3p, miR-29b-3p and miR-1183) to be differentially expressed between ACM patient-derived and control CStCs and confirmed the de-regulation of miR-29b-3p in an extended data set including samples from in total 8 ACM patients and 5 controls. We also found evidence for a potential contribution of miR-29b-3p to ACM development by targeting proteins involved in extracellular matrix organization.
Genes differentially expressed between ACM and control samples
In order to delineate differences in the gene expression profiles between ACM and control CStCs, we subjected poly-A enriched RNA from 3 ACM and 3 control samples from different individuals to high throughput sequencing for whole genome expression analysis. Gene expression quantification was performed by counting the reads completely aligned within the exon boundaries of any transcript. Most of these reads were assigned to protein coding genes but a considerable number was also assigned to genes from other biotypes such as lincRNAs or snRNAs, allowing quantification of some non-coding genes in the present experiments too (see Figure S1 in Additional file 2 for an overview of read counts per gene biotype). A principal component analysis (PCA) performed on the gene expression profiles showed the expected separation of ACM and CTRL samples (Fig. 2a), but did also suggest the presence of an unknown confounding factor represented by PC1. Including a variable describing the above factor into the linear model and thus alleviating for its potential influence, we performed gene-wise tests for differential expression and identified 272 genes more than two-fold differentially expressed at a 5% FDR (i.e. with an padj < 0.05; Fig. 2b and Table S10 in Additional file 2; Additional file 3 for the complete results).
For validation, we selected 4 significantly differentially expressed genes (PCDH1, SAXO2, SEMA3D, NEDD9). These genes were selected because they were among the most significant differentially expressed genes with a previously described relation to cardiac diseases. For validation of the RNA-seq results we assessed their differential expression in RNA samples from the same individuals, but from cells from different passages and experiments (to ensure the results being independent of the passage/experiment). Regulation of all 4 genes was reconfirmed (Additional file 2: Table S11). To test and validate the biological significance of their regulation in the context of ACM, we evaluated differential expression of the genes in an extended data set consisting of RNA from CStCs derived from 7 ACM patients and 7 control individuals (Table 4).
Regulation of SAXO2 and NEDD9 was reconfirmed also on the larger data set, while PCDH1 and SEMA3D, despite being on average more than two-fold, did not reach significance levels.
Subsequent Gene Ontology enrichment analyses of the differentially expressed genes revealed significant enrichment in cellular component terms “integral component of plasma membrane” and “intrinsic component of membrane” (Additional file 2: Table S12) and biological process terms “dopamine transport” and “catecholamine secretion” (Additional file 2: Table S13). Together, these results indicate an enrichment of membrane associated genes among the identified differentially expressed genes.
We next evaluated the expression of genes with ACM-associated variants  and, with the exception of JUP, did not find any of these genes being differentially expressed (Additional file 2: Table S14). Interestingly, five of them, CTNNA3, DES, DSC2, PLN and RYR2, were not expressed, or expressed below the detection limit, in any of the analyzed samples. From the remaining genes, 4 were highly expressed (LMNA, DSP, TMEM43 and JUP), and 4 moderately (TTN, TGFB3, DSG2 and PKP2). Expression of JUP was highly variable between within-group samples and was on average higher in ACM than in CTRL samples. This was reconfirmed by real-time RT-qPCR measurements in the same RNA (Additional file 2: Figure S2). The observed differential expression did however disappear if additional samples were considered (Additional file 2: Figure S2). Thus, in line with previous reports , no biologically significant differential expression of ACM associated genes could be observed on mRNA levels.
Taken together, the results from the gene expression analysis are in agreement with the miRNA target gene analysis, indicating that an integral part of the deregulated genes are functional in the outer cell membrane.
Combined miRNA and gene expression analysis
Given that some miRNAs are known to repress their target genes mainly by destabilizing their mRNA , we combined the miRNA and gene expression data set to evaluate whether a downregulation of target genes of the identified ACM-specific miRNAs is already observable on mRNA levels. Indeed, the average log2 fold-change values (M-values) of target genes for most miRNAs were slightly negative, representing an on average lower expression of these genes in ACM than in control samples (Additional file 2: Tables S16-S17). To test whether this observed repression was larger than expected by chance, we compared it to average M-values of same-sized sets of randomly selected target genes of expressed miRNAs. The resulting p-values from these tests are listed in Table 5. Of these, the average M-values of miR-29b-3p target genes and the combination of all of the miRNAs’ target genes were significant suggesting that the identified miRNAs are functional and that a tendency of repression is already detectable on mRNA levels.
Most human miRNAs are encoded in introns of protein coding genes and share in many instances the primary transcription unit with these . Using our mirhostgenes R-package we identified potential host genes for 2 of the ACM-specific miRNAs: the two lincRNAs C1orf132 and AC058791.1 for miR-29b-3p and the protein coding gene SP4 for miR-1183. While none of these were significantly differentially expressed, they all showed the expected trend of a higher expression in ACM samples (Additional file 2: Table S18). According to the protein expression atlas , SP4 is expressed in heart muscle tissue. No protein expression data is available for the other two as they are non-coding. According to the Ensembl Expression Atlas however, also C1orf132 is highly expressed in heart tissue. These miRNAs might therefore, in addition to CStCs, also be expressed in other cell types of the heart.
In conclusion, our analysis on the expression of the miRNAs’ targets genes on mRNA levels suggests the identified miRNAs being functional in the analyzed cells which strengthens the evidence for their potential involvement in ACM.
Gene network analysis
We next performed gene network analyses to evaluate relations between regulated genes, target genes of differentially expressed miRNAs and genes with known ACM-related genetic variants. Such networks allow to identify groups of input gene products that physically interact, that are part of common protein complexes, or that are functionally related. We built a first network based on functional annotations (functional network) of the genes mentioned above. Connected nodes in this network represent thus genes implicated in common biological processes. Clusters of highly connected components were identified and biological processes characteristic for each cluster were determined. Overlapping clusters were further combined resulting in a total of 6 functional modules (FM; resulting network model shown in Fig. 3; see Additional file 4 for module details). The largest, FM_1, consisted of a large number of miRNA targets that are implicated in extracellular matrix organization. Module FM_2 contained multiple transcription factors that are up- or down-regulated in ACM. Module FM_4 included the largest number of differentially expressed genes (9), four of which encoding lipid transport proteins (ABCC3, ATP8A1, ESYT3 and ABCG1). Genes in FM_5 are related to cardiac muscle cell adhesion and most of them encode desmosomal proteins. Modules FM_3 and FM_6 contained mostly genes involved in signaling processes (small GTPase mediated signal transduction for FM_3 and G-protein coupled receptor signaling for FM_6).
A second gene network was generated using functional relationships from STRING and protein interactions from BioPlex and mentha (combined network). For this analysis we included all genes with a valid read count in the present dataset or which, according to the protein atlas , are expressed in the heart tissue. Highly connected clusters of genes were identified as above and those with more than one differentially expressed gene, ACM-miRNA target gene or gene with ACM-related variants were further investigated. The full results are shown online at https://gemex.eurac.edu/bioinf/acm/ (“Clusters in full networks”). Some of the identified clusters showed considerable overlap with the previously identified functional modules, such as c0515 that overlaps FM_5, or the collagen genes containing cluster c0502 resembling FM_1. Additional promising clusters were identified too, like the cell adhesion related cluster c0997 that contained two genes from the protocadherin family (PCDH12 and PCDH17), both being lower expressed in ACM than in control samples, or cluster c0065 that links the differentially expressed gene ITGA10 with the ACM-miRNA target gene ITGB1 and with all of its 10 genes being also involved in extracellular matrix organization. Another interesting cluster (c0703) consisted of 16 genes, 6 of which were related to cell adhesion. From these genes, two, the ephrin receptors EPHA4 and EPHB2, were up-regulated and one, the EFNB class ephrin ligand EFNB2, was down-regulated in ACM cells (Additional file 2: Figure S3).
Summarizing, the functional networks provide an overview of the main biological processes in which both the ACM differentially expressed genes and the ACM-specific miRNA targets play a role. The results indicate that, in addition to cell adhesion related processes, deregulated genes are also involved in lipid transport, inflammation and fibrosis related processes.
In this study, we performed for the first time a comprehensive characterization of the coding and non-coding transcriptome of cardiac stromal cells derived from ACM patients and control individuals. Expression studies in ACM have been performed in the past [9,10,11], but these were either focused on protein coding genes or on miRNAs and were conducted on RNA extracted directly from heart tissue biopsies or even from explanted hearts. Of the 199 deregulated genes from Gaertner et al.  we could detect 191 in our samples, but none of them were differentially expressed. Similarly, none of the 3 genes reported by Akdis et al.  were deregulated in our samples. Gaertner et al. performed expression profiling in RNA extracted directly from explanted (final stage) ACM- or non-failing hearts, Akdis et al. in RNA from biopsies, that were even reported to be inhomogeneous in the original article . Missing concordance between studies can thus most likely be attributed to the different sample types being studied, different cell types, or to inhomogeneous cell populations. By performing all analyses in cells of the same type we ensured the results being independent of cell type or tissue-specific differences, limiting the validity of our conclusions however to the cardiac stromal cells analyzed. Recently, 28 miRNAs have been described to be differentially expressed in mouse cardiac muscle cells (HL-1) after shRNA-mediated PKP2 knock-down . We evaluated the expression for the human equivalents for 12 of these miRNAs in our CStCs. Only two of them, miR-671-3p and miR-487b-3p showed a tendency of differential expression, albeit not reaching significance levels (Additional file 2: Table S15). Also in this case, differences can be explained by different cell types being studied, different species and also the fact that PKP2 knock-down was artificially induced by shRNA silencing .
In our analysis, we identified 3 miRNAs to be differentially expressed between CStCs from ACM and control samples. Two of them (miR-29b and miR-1183) have previously been related to a pathway associated with ACM or other cardiac diseases. miRNAs from the miR-29 family have been shown to play a role in cardiac fibrosis after myocardial infarction targeting collagen, fibrillin and elastin genes . Recently, also desmocollin-2 (DSC2), was identified as a direct target of miR-29b in mouse keratinocytes , suggesting a potential direct involvement of this miRNA in ACM. miR-1183 on the other hand was identified as a putative serum biomarker for rheumatic heart disease . Analyses on the validated target genes of the ACM-specific miRNAs revealed their involvement in ACM related pathways, such as “extracellular matrix organization”. These results, along with the literature mentioned above, provide evidence for a potential prominent role of miRNAs in the ACM pathogenesis or phenotype maintenance driven by CStCs.
Among the significantly differentially expressed genes were also some potentially interesting candidates: PCDH1, the strongest differentially expressed gene, highly expressed in control and very low in ACM samples, encodes the protocadherin 1 membrane protein. Elevated expression of PCDH1 was found to be related to higher cell aggregation activity . Among the strongest up-regulated genes in ACM was also SEMA3D which has been implicated in vascular development and which can inhibit cell motility and tubulogenesis . PCDH1 and SEMA3D were however not de-regulated in all tested cases of the extended validation set, possibly due to a different genetic background of the analyzed ACM cases, and suggesting that different pathways might contribute to ACM in its different subtypes. Another interesting differentially expressed gene was NEDD9, a focal adhesion protein that acts as a scaffold to regulate signaling complexes important in cell attachment, migration and invasion . Fibroblasts from NEDD9 −/− mice showed significantly decreased adhesion strength  while over-expression of NEDD9 was shown to enhance cell speed and haptotaxis towards fibronectin . NEDD9 is also a validated target gene of the ACM-specific miRNA miR-29-3p and of miR-18a-5p, a miRNA also more than two-fold up-regulated in ACM which regulation did however not reach significance levels (adjusted p-value = 0.08). Its more than two-fold down-regulation might hence be mediated or consolidated by these miRNAs.
In analyzing the expression of genes with ACM-associated variants we found that five of the 13 analyzed genes, namely CTNNA3, DES, DSC2, PLN and RYR2, were not expressed. Therefore, potential pathogenic effects of mutations in those genes is unlikely to be mediated by CStCs.
By evaluating the expression of target genes from the ACM-specific miRNAs we found them to be significantly, albeit only slightly, lower expressed in ACM samples. This moderate downregulation was not completely unexpected given the comparably low expression of the miRNAs and the current opinion of miRNAs acting predominantly through fine-tuning gene expression [43, 44].
The functional network analysis resulted in some interesting clusters of related genes. Among these was the functional module FM_4 containing 9 differentially expressed genes, four of which encode lipid transport proteins. These proteins might play a role in the lipid accumulation described in ACM CStCs. Genes from two other modules, FM_1 containing collagen encoding genes and FM_6 inflammatory response related genes, might be involved in the molecular events leading to fibrosis in ACM. Regarding the analysis of the combined network, an interesting cluster was c0703, a cluster consisting of 16 genes with 13 being involved in ephrin receptor signaling. The membrane bound Eph receptors and ligands are involved in a wide array of developmental processes including cardiovascular, with one of their most prominent biological functions being the modulation of cell adhesion . EphB receptors and ephrins are part of the genetic program activated by β–catenin/Wnt signaling . Eph-ephrin signaling has also been implicated in cardiac stem cell migration into injured tissue after myocardial infarction  and has even been proposed as new therapeutic target for that disease . Also an involvement in the modulation of electrical coupling between cardiomyocytes has been reported .
Limitations of the study
At last some limitations of the present study have to be acknowledged. The results presented here base on cultured cells and may hence not be directly expanded to CStCs in their original microenvironment. Different mutations in desmosomal genes might also have an impact on the expression profiles of the cells. Also, it is currently unclear if ACM pathogenesis and disease progression is variant-specific or caused by a common canonical pathway. The small numbers of available ACM cases without any or with common mutations precluded however such mutation-specific differential expression analyses in the present data set. In addition, individuals in the present study do not share a similar genotype and thus any changes possibly identified between samples can not be attributed exclusively to specific variants or sets of variants in specific ACM-related genes. This interesting aspect could however be analyzed in future studies based on samples from individuals with similar genetic background, such as affected and non-affected individuals from the same family.
The small number of biological replicates used in the genome wide screens represent another limitation of the study and hence not all of the results might apply to ACM in general. Nevertheless, we could validate the main findings in a larger set of ACM patients and controls. Still, due to the limited number of samples, the present study should be considered as initial, hypothesis-generating only and further studies evaluating the interconnection between the networks regarding the progression of the disease are needed to fully characterize ACM pathogenesis. Regarding the low number of replicates, it is however also important to highlight that ACM is a rare disease with low incidence in the population [1, 2] and that it is difficult to obtain CStCs from ACM patients due to the small rate of patients that undergo cardiac biopsies.
Selection of appropriate control samples is crucial for any study. All control samples used for transcriptome analysis and miRNA profiling derive from deceased individuals and might thus not be completely identical to biopsies from living hearts. By taking the samples from deceased organ donors without heart disease, which ensures special care to the explanted organs, we tried to minimize any potential biases. CStCs from control samples were also generated using the same protocols and cultured under the same conditions as CStCs from ACM individuals.
Due to the phenotypic overlap of ACM with those of other cardiomyopathies, particularly with the arrhythmic form of idiopathic dilated cardiomyopathy (DCM) , we can not completely exclude the miRNA and gene expression patterns presented here being shared across some cardiomyopathies. Recently, genetic overlaps between ACM and DCM have been described: mutations in genes encoding desmosomal proteins have been reported in DCM cases  and mutations in TTN , PLN  and FLNC  have been associate with both ACM and DCM. Abnormalities such as epsilon waves and right ventricular dilation, considered classical hallmarks of ACM, have been described recently in patients with Brugada syndrome [54, 55], along with mutations in the PKP2 gene . We can therefore not exclude that the identified molecular signatures are specific for ACM and do not overlap, at least to a certain degree, with signatures in other primary cardiomyopathies. Still, for the potential overlap with DCM, it has to be noted that the shared phenotypic and structural overlap of both diseases is seen in end-stage heart failures, while the results in this study have been performed in early-stage ACM patients.
Finally, all of the analyzed individuals in this study are male, and some conclusions or results might thus not be applicable to female ACM patients. Despite the lower prevalence of ACM in females, it would however be of utmost interest, also in the light of gender medicine, to dissect in future studies gender-specific differences in ACM disease development and progression.
We determined for the first time the coding and non-coding transcriptome of ACM in cardiac stromal cells. We found evidence for a potential prominent role of miRNAs in ACM pathogenesis or phenotype maintenance. Due to the small sample sizes used for whole genome expression profiling, the presented results might however not be extended to all cell types involved in ACM or to all different subtypes of ACM. Taken together, the comprehensive data set and the results presented here might serve as a resource for future studies elucidating the role of the identified genes, miRNAs or pathways in ACM.
Cardiac stromal cells
False discovery rate
Principal component analysis
TaqMan Low Density Arrays;
Sen-Chowdhry S, Morgan RD, Chambers JC, McKenna WJ. Arrhythmogenic cardiomyopathy: etiology, diagnosis, and treatment. Annu Rev Med. 2010;61:233–53.
Thiene G, Nava A, Corrado D, Rossi L, Pennelli N. Right ventricular cardiomyopathy and sudden death in young people. N Engl J Med. 1988;318:129–33.
Lazzarini E, Jongbloed JDH, Pilichou K, Thiene G, Basso C, Bikker H, et al. The ARVD/C genetic variants database: 2014 update. Hum Mutat. 2015;36:403–10.
Garcia-Gras E, Lombardi R, Giocondo MJ, Willerson JT, Schneider MD, Khoury DS, et al. Suppression of canonical Wnt/beta-catenin signaling by nuclear plakoglobin recapitulates phenotype of arrhythmogenic right ventricular cardiomyopathy. J Clin Invest. 2006;116:2012–21.
Chen SN, Gurha P, Lombardi R, Ruggiero A, Willerson JT, Marian AJ. The hippo pathway is activated and is a causal mechanism for adipogenesis in arrhythmogenic cardiomyopathy. Circ Res. 2014;114:454–68.
Pilichou K, Thiene G, Bauce B, Rigato I, Lazzarini E, Migliore F, et al. Arrhythmogenic cardiomyopathy. Orphanet J Rare Dis. 2016;11:33.
Rossini A, Frati C, Lagrasta C, Graiani G, Scopece A, Cavalli S, et al. Human cardiac and bone marrow stromal cells exhibit distinctive properties related to their origin. Cardiovasc Res. 2011;89:650–60.
Sommariva E, Brambilla S, Carbucicchio C, Gambini E, Meraviglia V, Russo Dello A, et al. Cardiac mesenchymal stromal cells are a source of adipocytes in arrhythmogenic cardiomyopathy. Eur Heart J. 2016;37:1835–46.
Gaertner A, Schwientek P, Ellinghaus P, Summer H, Golz S, Kassner A, et al. Myocardial transcriptome analysis of human arrhythmogenic right ventricular cardiomyopathy. Physiol Genomics. 2012;44:99–109.
Akdis D, Medeiros-Domingo A, Gaertner-Rommel A, Kast JI, Enseleit F, Bode P, et al. Myocardial expression profiles of candidate molecules in patients with arrhythmogenic right ventricular cardiomyopathy/dysplasia compared to those with dilated cardiomyopathy and healthy controls. Heart Rhythm. 2016;13:731–41.
Zhang H, Liu S, Dong T, Yang J, Xie Y, Wu Y, et al. Profiling of differentially expressed microRNAs in arrhythmogenic right ventricular cardiomyopathy. Sci Rep. 2016;6:28101.
Mestdagh P, Van Vlierberghe P, De Weer A, Muth D, Westermann F, Speleman F, et al. A novel and universal method for microRNA RT-qPCR data normalization. Genome Biol. 2009;10:R64.
Smyth GK. Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol. 2004;3(Art 3):29. (electronic)
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J Roy Statist Soc Ser B. 1995;57:289–300.
Wu TD, Nacu S. Fast and SNP-tolerant detection of complex variants and splicing in short reads. Bioinformatics. 2010;26:873–81.
Lawrence M, Huber W, Pagès H, Aboyoun P, Carlson M, Gentleman R, et al. Software for computing and annotating genomic ranges. PLoS Comput Biol. 2013;9:e1003118.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.
Schmittgen TD, Livak KJ. Analyzing real-time PCR data by the comparative C(T) method. Nat Protoc. 2008;3:1101–8.
Falcon S, Gentleman R. Using GOstats to test gene lists for GO term association. Bioinformatics. 2007;23:257–8.
Chou C-H, Chang N-W, Shrestha S, Hsu S-D, Lin Y-L, Lee W-H, et al. miRTarBase 2016: updates to the experimentally validated miRNA-target interactions database. Nucleic Acids Res. 2016;44:D239–47.
Rainer J, Ploner C, Jesacher S, Ploner A, Eduardoff M, Mansha M, et al. Glucocorticoid-regulated microRNAs and mirtrons in acute lymphoblastic leukemia. Leukemia. 2009;23:746–52.
Huttlin EL, Ting L, Bruckner RJ, Gebreab F, Gygi MP, Szpyt J, et al. The BioPlex network: a systematic exploration of the human Interactome. Cell. 2015;162:425–40.
Calderone A, Castagnoli L, Cesareni G. Mentha: a resource for browsing integrated protein-interaction networks. Nat Meth. 2013;10:690–1.
Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, et al. STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 2015;43:D447–52.
Weichenberger CX, Palermo A, Pramstaller PP, Domingues FS. Exploring approaches for detecting protein functional similarity within an Orthology-based framework. Sci Rep. 2017;7:381.
Uhlén M, Fagerberg L, Hallström BM, Lindskog C, Oksvold P, Mardinoglu A, et al. Proteomics. Tissue-based map of the human proteome. Science. 2015;347:1260419.
Nepusz T, Yu H, Paccanaro A. Detecting overlapping protein complexes in protein-protein interaction networks. Nat Meth. 2012;9:471–2.
Csardi G, Nepusz T. The igraph software package for complex network research. InterJ Complex Syst. 2006;1695:1-9.
Franz M, Lopes CT, Huck G, Dong Y, Sumer O, Bader GD. Cytoscape.Js: a graph theory library for visualisation and analysis. Bioinformatics. 2016;32:309–11.
Weichenberger CX, Blankenburg H, Palermo A, D'Elia Y, König E, Bernstein E, et al. Dintor: functional annotation of genomic and proteomic data. BMC Genomics. 2015;16:1081.
Rampazzo A, Calore M, van Hengel J, van Roy F. Intercalated discs and arrhythmogenic cardiomyopathy. Circ Cardiovasc Genet. 2014;7:930–40.
Guo H, Ingolia NT, Weissman JS, Bartel DP. Mammalian microRNAs predominantly act to decrease target mRNA levels. Nature. 2010;466:835–40.
Ha M, Kim VN. Regulation of microRNA biogenesis. Nat Rev Mol Cell Biol. 2014;15:509–24.
Gurha P, Chen X, Lombardi R, Willerson JT, Marian AJ. Knockdown of Plakophilin 2 downregulates miR-184 through CpG Hypermethylation and suppression of the E2F1 pathway and leads to enhanced Adipogenesis in vitro. Circ Res. 2016;119:731–50.
van Rooij E, Sutherland LB, Thatcher JE, DiMaio JM, Naseem RH, Marshall WS, et al. Dysregulation of microRNAs after myocardial infarction reveals a role of miR-29 in cardiac fibrosis. PNAS. 2008;105:13027–32.
Kurinna S, Schäfer M, Ostano P, Karouzakis E, Chiorino G, Bloch W, et al. A novel Nrf2-miR-29-desmocollin-2 axis regulates desmosome function in keratinocytes. Nat Commun. 2014;5:5099.
Li N, Lian J, Zhao S, Zheng D, Yang X, Huang X, et al. Detection of differentially expressed MicroRNAs in rheumatic heart disease: miR-1183 and miR-1299 as potential diagnostic biomarkers. Biomed Res Int 2015;2015:524519–11.
Sano K, Tanihara H, Heimark RL, Obata S, Davidson M, St John T, et al. Protocadherins: a large family of cadherin-related molecules in central nervous system. EMBO J. 1993;12:2249–56.
Aghajanian H, Choi C, Ho VC, Gupta M, Singh MK, Epstein JA. Semaphorin 3d and semaphorin 3e direct endothelial motility through distinct molecular signaling pathways. J Biol Chem. 2014;289:17971–9.
Tikhmyanova N, Little JL, Golemis EA. CAS proteins in normal and pathological cell growth control. Cell Mol Life Sci. 2010;67:1025–48.
Zhong J, Baquiran JB, Bonakdar N, Lees J, Ching YW, Pugacheva E, et al. NEDD9 stabilizes focal adhesions, increases binding to the extra-cellular matrix and differentially effects 2D versus 3D cell migration. PLoS One. 2012;7:e35058.
Fashena SJ, Einarson MB, O'Neill GM, Patriotis C, Golemis EA. Dissection of HEF1-dependent functions in motility and transcriptional regulation. J Cell Sci. 2002;115:99–111.
Grosshans H, Filipowicz W. Molecular biology: The expanding world of small RNAs. Nature. 2008;451:414–6.
Nazarov PV, Reinsbach SE, Muller A, Nicot N, Philippidou D, Vallar L, et al. Interplay of microRNAs, transcription factors and target genes: linking dynamic expression changes to function. Nucleic Acids Res. 2013;41:2817–31.
Arvanitis D, Davy A. Eph/ephrin signaling: networks. Genes Dev. 2008;22:416–29.
Palmer A, Klein R. Multiple roles of ephrins in morphogenesis, neuronal networking, and brain function. Genes Dev. 2003;17:1429–50.
Goichberg P, Bai Y, D'Amario D, Ferreira-Martins J, Fiorini C, Zheng H, et al. The ephrin A1-EphA2 system promotes cardiac stem cell migration after infarction. Circ Res. 2011;108:1071–83.
O'Neal WT, Griffin WF, Dries-Devlin JL, Kent SD, Chen J, Willis MS, et al. Ephrin-Eph signaling as a potential therapeutic target for the treatment of myocardial infarction. Med Hypotheses. 2013;80:738–44.
Ishii M, Mueller I, Nakajima T, Pasquale EB, Ogawa K. EphB signaling inhibits gap junctional intercellular communication and synchronized contraction in cultured cardiomyocytes. Basic Res Cardiol. 2011;106:1057–68.
Akdis D, Brunckhorst C, Duru F, Saguner AM. Arrhythmogenic Cardiomyopathy: Electrical and structural phenotypes. Arrhythm Electrophysiol Rev 2016;5:90–101.
Elliott P, O'Mahony C, Syrris P, Evans A, Rivera Sorensen C, Sheppard MN, et al. Prevalence of desmosomal protein gene mutations in patients with dilated cardiomyopathy. Circ Cardiovasc Genet. 2010;3:314–22.
van der Zwaag PA, van Rijsingen IAW, Asimaki A, Jongbloed JDH, van Veldhuisen DJ, Wiesfeld ACP, et al. Phospholamban R14del mutation in patients diagnosed with dilated cardiomyopathy or arrhythmogenic right ventricular cardiomyopathy: evidence supporting the concept of arrhythmogenic cardiomyopathy. Eur J Heart Fail. 2012;14:1199–207.
Ortiz-Genga MF, Cuenca S, Dal Ferro M, Zorio E, Salgado-Aranda R, Climent V, et al. Truncating FLNC mutations are associated with high-risk dilated and Arrhythmogenic cardiomyopathies. J Am Coll Cardiol. 2016;68:2440–51.
Catalano O, Antonaci S, Moro G, Mussida M, Frascaroli M, Baldi M, et al. Magnetic resonance investigations in Brugada syndrome reveal unexpectedly high rate of structural abnormalities. Eur Heart J. 2009;30:2241–8.
Letsas KP, Efremidis M, Weber R, Korantzopoulos P, Protonotarios N, Prappa E, et al. Epsilon-like waves and ventricular conduction abnormalities in subjects with type 1 ECG pattern of Brugada syndrome. Heart Rhythm. 2011;8:874–8.
Cerrone M, Lin X, Zhang M, Agullo-Pascual E, Pfenniger A, Chkourko Gusky H, et al. Missense mutations in Plakophilin-2 can cause Brugada syndrome phenotype by decreasing sodium current and Nav1.5 membrane localization. Heart Rhythm. 2013;10:1743.
We thank all patients consenting to participate in this study. We thank Dr. Deborah Mascalzoni for ethical and legal advice. We thank Drs. Veronica De Sanctis and Roberto Bertorelli (NGS Facility, CIBIO-LaBSSAH, University of Trento) for NGS sequencing and helpful discussions. We thank the Department of Innovation, Research and University of the Autonomous Province of Bozen/Bolzano for covering the Open Access publication costs.
This work was supported by the Department of Innovation, Research and Universities of the Autonomous Province of Bolzano-South Tyrol and by 5 × 1000 from Italian Ministry of Health [2012 to Centro Cardiologico Monzino IRCCS]. The funding bodies had no role in the design of the study and collection, analysis and interpretation of data and in writing the manuscript.
Availability of data and materials
The RNA-seq dataset supporting the conclusions of this article is available in the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/, accession number GSE95120), with raw read data available in the Short Read Archive (https://www.ncbi.nlm.nih.gov/sra, accession number SRP100426). The complete results from the miRNA and gene expression analyses are included within the article (and its additional files). The results from the gene network analyses are available at https://gemex.eurac.edu/bioinf/acm/.
The authors declare that they have no competing interests.
Ethics approval and consent to participate
The study complies with the Declaration of Helsinki and was approved by the Ethic Committees of Centro Cardiologico Monzino (on 04/11/2015). Control samples have been obtained from cadaveric donors at Fondazione Banca dei Tessuti di Treviso, according to the recommendations of the Italian Transplant Network. Written informed consent was obtained from all participants.
Consent for publication
A “consent for publication” was given to each participant. Clinical and personal information is only provided for participants who signed this consent.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
miRNA expression results. Contains the complete data and results from the miRNA expression profiling. (XLSX 82 kb)
Supplementary information. Contains supplementary information, tables and figures. (PDF 1224 kb)
Gene expression results. Contains gene quantification data and results from the transcriptome profiling. (XLSX 4484 kb)
Network analysis results. Contains functional module definitions and results from the functional network analysis. (XLSX 17 kb)
About this article
Cite this article
Rainer, J., Meraviglia, V., Blankenburg, H. et al. The arrhythmogenic cardiomyopathy-specific coding and non-coding transcriptome in human cardiac stromal cells. BMC Genomics 19, 491 (2018). https://doi.org/10.1186/s12864-018-4876-6
- Arrhythmogenic cardiomyopathy
- Cardiac stromal cells
- microRNA expression
- Gene expression
- Transcriptome profiling