Abstract
The expression of long noncoding RNAs (lncRNAs) in human heart failure (HF) has not been widely studied. Using RNA sequencing (RNA-Seq), we compared lncRNA expression in 22 explanted human HF hearts with lncRNA expression in 5 unused donor human hearts. We used Cufflinks to identify isoforms and DESeq to identify differentially expressed genes. We identified the noncoding RNAs by cross-reference to Ensembl release 73 (Genome Reference Consortium human genome build 37) and explored possible functional roles using a variety of online tools. In HF hearts, RNA-Seq identified 84,793 total messenger RNA coding and noncoding different transcripts, including 13,019 protein-coding genes, 2,085 total lncRNA genes, and 1,064 pseudogenes. By Ensembl noncoding RNA categories, there were 48 lncRNAs, 27 pseudogenes, and 30 antisense RNAs for a total of 105 differentially expressed lncRNAs in HF hearts. Compared with donor hearts, HF hearts exhibited differential expression of 7.7% of protein-coding genes, 3.7% of lncRNAs (including pseudogenes), and 2.5% of pseudogenes. There were not consistent correlations between antisense lncRNAs and parent genes and between pseudogenes and parent genes, implying differential regulation of expression. Exploratory in silico functional analyses using online tools suggested a variety of possible lncRNA regulatory roles. By providing a comprehensive profile of right ventricular polyadenylated messenger RNA transcriptome in HF, RNA-Seq provides an inventory of differentially expressed lncRNAs, including antisense transcripts and pseudogenes, for future mechanistic study.
Keywords
The results reported from the Encyclopedia of DNA Elements (ENCODE) project are transforming our nascent understanding of integrated human genomics. 1 Fully 75% of the human genome is transcribed into some type of RNA, 2 although only 1.22% of the genome encodes the exons that comprise the 20,687 known protein-coding genes. 1 The vast majority of the human genome, therefore, is transcribed into non–protein-coding RNAs, including the 8,801 small RNAs, 9,640 long noncoding RNAs (lncRNAs; defined as nontranslated RNA <200 base pairs [bp] in length) and 11,224 pseudogenes identified to date. 1 Evolving research has identified diverse epigenetic regulatory roles for lncRNAs in development, homeostasis, and disease. 3 lncRNAs are likely to play important regulatory roles in human heart failure as well. 4
Via high-throughput sequencing, RNA sequencing (RNA-Seq) provides a comprehensive profile of polyadenylated protein-coding messenger RNA (mRNA) and non–protein-coding lncRNA transcripts. 5 We used RNA-Seq to generate the lncRNA transcriptional profile of the right ventricle in end-stage human heart failure. Human right ventricular (RV) tissue may be obtained at either endomyocardial biopsy or surgery. Diagnostic endomyocardial biopsy, given its attendant risks, is rarely performed for patients with pulmonary arterial hypertension or isolated RV heart failure and yields small amounts of tissue. Ex-planted RV human myocardium, by contrast, is readily obtainable at cardiac transplantation and yields abundant tissue collectable under strictly controlled conditions to prevent RNA degradation. We thus examined the human RV transcriptome in hearts explanted at the time of transplantation for end-stage left ventricular (LV) failure.
METHODS
Clinical characterization of study subjects
De-identified clinical data were abstracted for all study subjects. RV contractile function was classified as normal or dysfunctional on the basis of the presence of mild or greater RV hypokinesis by visual inspection on the echocardiogram performed most closely to explantation. RV size and function were also quantified by tricuspid annular plane systolic excursion, the index of myocardial performance, and fractional area change according to American Society of Echocardiography guidelines. Pulmonary hypertension was defined as a mean pulmonary artery pressure greater than 25 mmHg by the right heart catheterization performed most closely to explantation.
RNA preparation
The myocardial samples were free of macroscopic infarction or fibrosis and were collected from approximately the same cross-sectional region of the unused donor and heart failure RVs. Ninety to 120 mg of tissue per study subject was used for RNA extraction. Total RNA was isolated using miRNeasy Mini kit (Qiagen) following the manufacturer's instructions.
RNA-Seq
We first performed a quality check of the input total RNA by running an aliquot on the Agilent Bioanalyzer to confirm integrity, and a Qubit RNA fluorometry assay was used to measure concentration. For each library, 100 ng of total RNA underwent enrichment for poly(A)-containing mRNA using poly(T) oligo-attached magnetic beads. Following purification, eluted poly(A) RNA was cleaved into small fragments of 120–210 bp (divalent cations, elevated temperature). The cleaved RNA fragments were copied into first-strand complementary DNA (cDNA) using SuperScript II reverse transcriptase and random primers. This was followed by second-strand cDNA synthesis using DNA polymerase I and RNase H. The cDNA fragments then underwent an end-repair process, the addition of a single “A” base, and ligation of the Illumina multiplexing adapters. The products were purified and enriched with polymerase chain reaction (PCR) to create the final cDNA sequencing library. The cDNA library underwent quality control check via an Agilent Bioanalyzer HS DNA assay to confirm final library size and via Agilent Mx3005P quantitative PCR (qPCR) machine using the KAPA Illumina library quantification kit to determine concentration. From a 2-nM stock, samples were pooled by molarity for multiplexing. From the pool, 12 pM was loaded into each well for the flow cell on the Illumina cBot for cluster generation. The flow cell was loaded onto the Illumina HiSeq 2500 using v3 chemistry and HTA 1.8 and sequenced at paired-end 50 bp with a target of 30 million pass filter reads per library. The raw sequencing reads in BCL format were processed through CASAVA-1.8.2 for FASTQ conversion and demultiplexing. The RTA chastity filter was used, and only the pass filter reads were retained for any further analysis.
Analysis of gene and transcript expression
We first performed quality control on raw sequence-derived data to identify potential outliers before doing any advanced analysis by applying tools such as Fastx Toolkit and FastQC. RNA read alignment and mapping were performed by Bowtie/TopHat, and transcriptome reconstruction was performed by Cufflinks for both mRNA and lncRNA. For Cufflinks, a minimum RPKM (reads per kilobase of exon model per million mapped reads) value of ≥1 was required for further analysis. 7 Cufflinks was also used to detect and quantitate alterative spliced transcripts and isoforms.
For differential expression quantification of mRNA and lncRNA genes, 9 we used three read-count-based negative binomial methods, DESeq, 10 edgeR, 11 and baySeq 12 to explore the consistency of differential expression quantification across methods. Each method aligned then adjusted the unadjusted read-count tables before analysis. For each platform, inflated type I error due to multiple comparisons was adjusted, and all adjusted tests were significant at the two-sided significance level of P ≤ 0.05. To examine parsimony and exclude potential bias across the three methods, we also compiled a sum of ranks from the three read-count-based methods to identify differentially expressed genes (DEGs) consistently selected across the three platforms. Following inspection of this ranking across the three platforms, for final selection of DEGs, we elected to use the results afforded by DESeq. DESeq afforded the best compromise between permissive (edgeR) and restrictive (baySeq) identification of DEGs and afforded the best reconciliation with the sum-of-ranks approach. In reporting adjusted total read counts for individual genes, we divided the unadjusted total gene read counts in a given library by the total read counts for all genes in same library and then multiplied this by 1 × 10 6 .
Exploration of lncRNA, pseudogene, and antisense function
Exploratory characterization of selected lncRNAs and antisense RNAs was provided by the online lncRNA search engine and bioinformatics and annotation tool LNCipedia (http://www.lncipedia.org), 13 a database for human lncRNA transcripts and genes that provides basic transcript information, gene structure and sequence from the genome sequence assembly tool Ensembl, secondary structure from RNAfold image, protein-coding potential from Coding Potential Calculator (CPC; used to assess protein-coding potential of a given RNA transcript) and HMMER, Proteomics Identification (PRIDE) database search, and microRNA (miRNA; nontranslated RNA usually 18–24 bp in length) binding site predictions from MiRTarget2 (an online miRNA target search engine and bioinformatics tool). Exploratory functional annotation of selected pseudo-genes was provided by the on-line resource tool pseudoMap (http://pseudomap.mbc.nctu.edu.tw/). 14 The pseudoMap tool maps a known pseudogene to both (1) the predicted endogenous short interfering RNAs (esiRNAs, which may silence mRNA translation via RNA interference) produced from that pseudogene and (2) the predicted miRNA targets for that pseudogene (by binding miRNA, a pseudogene may act as an “miRNA decoy,” providing competitive endogenous RNA-mediated regulation of mRNA translation). In selected instances, either the DIANA-LncBase (http://www.microrna.gr/LncBase) 15 or Segal Lab of Computational Biology at the Weismann Institute of Science online tools (http://genie.weizmann.ac.il) 16 was used for computational prediction of miRNA targets on lncRNAs or parent genes. For exploring noncoding RNA-protein interactions, RNA-Protein Interaction Prediction (RPISeq), which uses both random forest (RF) and support vector machine (SVM) classifiers to predict RNA-protein interactions using only sequence information (http://pridb.gdcb.iastate.edu/RPISeq), 17 and catRAPID, an online in silico RNA and protein binding prediction tool that uses sequence information to evaluate the interaction propensities of polypeptide and nucleotide chains using physicochemical properties (http://s.tartagilalab.com/catrapid/omics), 18 were used.
RESULTS
Subjects
The clinical and echocardiographic characteristics of the subjects who underwent transplant appear in Table 1.
Clinical and echocardiographic characteristics of study subjects
Note: Data are mean value ± standard deviation, unless otherwise indicated. Asterisk indicates P < 0.05. ACEI: angiotensin-converting enzyme inhibitor; ARB: angiotensin receptor blocker; FAC: right ventricular fractional area change; HF: heart failure; PDE V: phosphodiesterase type V; RAP, right atrial pressure; RIMP: RV index of myocardial performance; RV: right ventricular; TAPSE: tricuspid annular plane systolic excursion.
Heart failure RV transcriptome
By the read count methods, there were 791,427,072 total paired-end, Ensembl-annotated unadjusted reads summed across the 32 cDNA libraries. The average number of aligned, Ensembl-annotated unadjusted reads per library was 24,732,096, with a range of reads per library from 18,607,371 to 35,024,568. There were 29,501,480, 25,357,404, and 23,506,029 average unadjusted reads in the unused donor LV, unused donor RV, and heart failure RV libraries, respectively.
There were 77,768 and 83,793 total mRNA coding and noncoding transcripts identified in the unused donor RVs and heart failure RVs, respectively (Table 2). The heart failure RVs expressed 13,019 protein-coding genes, 2,085 lncRNAs, and 1,064 pseudogenes, and all types of transcripts were more abundant in the heart failure RVs than unused donor RVs.
Cumulative number of genes and transcripts expressed compared with GENCODE
Note: For this study, a gene was considered expressed if the value for fragments per kilobase of exon was ≥1. ENCODE: Encyclopedia of DNA Elements; lncRNA: long noncoding RNA. LV: left ventricle; RV: right ventricle.
lncRNAs
There were 2,085 and 1,474 lncRNAs expressed in heart failure RVs and unused RVs, respectively (Table 3). In heart failure RVs versus unused donor RVs, 78 lncRNAs (48 lncRNAs and 30 natural antisense transcripts), representing 3.7% of the total 2,085 lncRNAs, were differentially expressed genes. The 48 lncRNA differentially expressed genes included intergenic, intronic, and processed transcripts. Among the 48 lncRNA heart failure RV differentially expressed genes, 35 were decreased in expression and 13 were increased in expression compared with unused donor RVs (Tables 3, 4). As a rule, lncRNAs exhibited a wide range of expression as denoted by the mean (± standard deviation) and maximum-minimum range (Table 4). The most abundantly expressed lncRNA in the heart failure RV group was RP11–206l10.11 (ENSG00000228794.3), a novel lncRNA with no predicted miRNA targets (Table 3).
Differentially expressed long noncoding RNAs (lncRNAs)
Note: Chrom: chromosome; Ensembl: Ensembl release 74; inf: value cannot be calculated; L2F: log2-fold change; Trans: transcripts.
Long noncoding RNAs (lncRNAs): unused donor (DON) right ventricle (RV) versus heart failure (HF) RV read counts and LNCipedia notes
Note: bp: base pairs; lnc: long noncoding; MAX: read count maximum; MEAN: read count mean; MIN: read count minimum; miRNA: microRNA; SD: read count standard deviation; L2F: log2-fold change.
LNCipedia is available at http://www.lncipedia.org. 13
Exploratory in silico functional analyses using online tools suggested a range of possible lncRNA functions. For example, the most differentially expressed lncRNA AP000783.2 (ENSG00000255414) is a 2,471-bp, 2-exon novel lncRNA transcribed from an intergenic region on ch11:123,325,106–123,331,118 forward strand. Although its function is unknown, AP000783.2 possesses 81 DIANA-predicted miRNA target sites with miRNA targeted gene (miTG) scores >0.7 (the higher the score, the higher the probability of miRNA binding), including target sites for miRNA-942 (score, 0.981), miRNA-580 (score 0.961), and miRNA-4760–3p (score, 0.946), which suggests that AP000783.2 may serve as an miRNA “decoy” for these or other miRNAs. The second most differentially expressed lncRNA, RP11–403B2.6, is a 748-bp, 2-exon transcript from ch15:20,964,640–20,980,123 reverse strand. Interestingly, it was not expressed in any of the five unused donor RVs or LVs and was expressed in only 4 of the 22 heart failure RVs with up to 60 and 292 read counts when expressed. It has no DIANA-predicted miRNA target sites, and its function is unknown. The third most differentially expressed lncRNA, RP11–60A24.3 (ENSG00000265542), is a novel lncRNA transcribed from ch17:55,849,307–55,912,110 reverse strand with 2 transcripts (a 2-exon, 472-bp transcript and a 3-exon, 505-bp transcript) due to alternative splicing. Interestingly, the lncRNA RP11–60A24.3 intronic region includes the sequence for a separate lncRNA gene, RN7SKP94 (ENSG00000222976, lncRNA 7SK small nuclear pseudogene 94). RP11–60A24.3 was decreased 2.5 log2-fold (P value 1.9E-07) in heart failure RVs, expressed in all unused donor RVs, and not expressed in one heart failure RV. Its function is unknown, but it possesses 13 DIANA-predicted miRNA binding sites with miTG scores >0.7, including target sites for miRNA-654–5p (score, 0.937), miRNA-1207–3p (score, 0.936), and miRNA-541–3p (0.864), which suggests that RP11–60A24.3 may play a role as miRNA “decoy.”
As an example of exploratory in silico functional analysis using online tools of a somewhat better characterized lncRNA, RP11–298021.5 (ENSG00000244277.1) is a 7-exon, 2,350-bp transcript with 7 alternatively spliced isoforms transcribed from chr3:167,613,736–167,645,799, an intergenic region conserved in mice but not fish between the SERPINI1 and GOLIM4 protein-coding genes. Although its function is unknown, based upon a predicted CPC score of 2.1, RP11–298021.5 possesses modest protein-coding potential. There are no known RP11–298021.5 associations with proteins based upon the absence of hits in the PRIDE database. MirTarget2 prediction identified 8 miR targets with scores >80, denoting the possibility of interaction with these miRs. These miR targets do not overlap with the MirTarget2-predicted miRs interacting with SERPINI1 or GOLIM4, suggesting that, if lnc-SERPINI1–4 acts as a miR decoy, it may not serve as a decoy for nearby SERPINI1 or GOLIM4 transcripts.
Antisense transcripts
There were 30 natural antisense transcript differentially expressed genes in heart failure RVs versus unused donor RVs (Table 5). Of these 30, 12 were increased in expression and 18 were decreased in expression in heart failure RVs. Interestingly, RP11–524H19.s (parent gene unknown) and RP11–597M17.1 (parent gene SCARA5) natural antisense transcripts were not expressed in any heart failure RV, and RP3–467K16.7 (parent gene FHAD1) natural antisense transcript was not expressed in any unused donor RV.
Differentially expressed natural antisense transcripts
Note: adj: adjusted; Chrom: chromosome; EHZ2: predicted binding to EZH2 component of polycomb 2 repressive complex; Ensembl: Ensembl release 74; L2F change: log2-fold change; IncRNA: long noncoding RNA; Trans: transcripts.
LNCipedia is available at http://www.lncipedia.org. 13
The relationship between parent gene and natural antisense transcription pair expression in heart failure RV versus unused donor RVs was complex (Table 6). Parent and natural antisense transcript were both decreased in 8 pairs, parent and natural antisense transcript were both increased in 6 pairs, parent increased and natural antisense transcript decreased in 7 pairs, parent decreased and natural antisense transcript increased in 3 pairs, parent unchanged with natural antisense transcript increased in 2 pairs, parent unchanged with natural antisense transcript decreased in 1 pair, and status unknown in 3 pairs (two parent genes unknown, KIAA1370 parent of RP11–23N2.4 not expressed) in heart failure RVs. Unused donor versus heart failure RV parent/natural antisense transcript ratios were significantly increased for 7 parent : natural antisense transcript pairs and significantly decreased in 4 parent : natural antisense transcript pairs (Table 6).
Differentially expressed natural antisense transcripts and parent sense genes
Note: CORR: read count correlation; DON RV: unused donor right ventricle; Ensembl: Ensembl release 73; HF RV: heart failure right ventricle; MEAN: read count mean; N/C: not calculated due to absence of sufficient read counts; P value: sense:antisense ratio P value by t test; SD: read count standard deviation.
Exploratory in silico analyses using online tools suggested multiple possible functional roles of natural antisense transcripts, as previously reported. 19 Of the 27 transcripts with entries in LNCipedia, 13 possessed protein-coding potential (no identified or known protein transcripts), and 5 possessed high-likelihood miRNA target sequences. For example, the GATA2 parent gene and its natural antisense transcript were both targets of miR 4722 (MiRTarget scores 0.77 and 0.82, respectively), and CCDC85A parent gene and its natural antisense transcript were both targets of miR 1184 (MiRTarget scores 0.76 and 0.64, respectively), suggesting that some natural antisense transcripts may serve as miR decoys. By RPISeq, 7 natural antisense transcripts were predicted with probability ≥0.80 by both RF and SVM classifiers to be binding partners for the polycomb repressive complex 2 (PRC2) EZH2 subunit. By the catRAPID algorithm, the FHAD1 natural antisense transcript had a greater strength of RNA : protein interaction with EZH2 than did Xist (Fig. 1), a natural antisense transcript previously reported to bind to and regulate PRC2 EZH2. 19 Interestingly, an atrial natriuretic peptide (NPPA) natural antisense transcript was among the natural antisense transcript differentially expressed genes, and NPPA natural antisense transcripts have been previously reported to modulate the expression of NPPA splicing isoforms. 20
Interaction strength with EZH2. By the catRAPID binding prediction algorithm, the FHAD1 NAT (A) has a greater strength of RNA: protein interaction with the PRC2 component EZH2 than does Xist (B), a long noncoding RNA previously reported to bind to and regulate the H3K27me3 trimethylation activity of EZH2.
Pseudogenes
There were 1,064 and 793 transcribed pseudogenes expressed in heart failure RVs and unused donor RVs, respectively (Table 2). There were 27 pseudogenes, with differentially expressed genes representing 2.5% of all pseudogenes in heart failure RVs versus unused donor RVs (Table 2). These included 14 processed, 10 unprocessed, 2 unitary, and 1 unknown type pseudogenes. Of the 27 total pseudogenes, 10 were transcribed in antisense, and 17 were transcribed in sense. The parental coding genes included 2 unknown genes, 2 zing finger protein pseudogenes (i.e., 2 pseudogenes had parental pseudogenes), and known parental protein-coding genes with a variety of disparate, functional proteins.
As previously reported, in general, the expression of the parental coding gene was more abundant than the expression of the corresponding pseudogenes (Table 7). 21 There were exceptions to this pattern, however, because the expression of the FBX043 pseudogenes (ENSG00000 236155, RP11–231P20.2), for example, was considerably more abundant than the expression of the parental FBX043 gene. The ratio of parental gene read count number to pseudogene read count number was highly variable between individual subjects and across different parental gene : pseudogene pairings and significantly different for several parental gene : pseudogene pairings between unused donor and heart failure RVs (Table 8).
Differentially expressed pseudogenes
Note: Chrom: chromosome; Ensembl: Ensemble release 73; esiRNA: endogenous silencing RNA; L2F adj: log2-fold change adjusted; miRNA: microRNA; Trans: transcripts.
Differentially expressed pseudogene and parent gene read counts
Note: CORR: read count correlation; DON RV: unused donor right ventricle; Ensembl: Ensembl release 73; MAX: read count maximum; MEAN: read count mean; MIN: read count minimum; SD: read count standard deviation; HF RV: heart failure right ventricle; N/C: not calculated due to absence of sufficient read counts; P value: parent : pseudogene (TPG) ratio P value by t test.
Functional annotation provided by the online resource tool “pseudoMap” yielded a variety of both predicted pseudogene esiRNA products and pseudogene miRNA binding targets (Table 7). From the available results in pseudoMap, the predicted pseudogene esiRNA products did not typically target the parent gene mRNAs, and the pseudogene miRNA binding partners did not typically overlap with the parent gene mRNA or miRNA binding partners. There were exceptions, however. For example, both the RPL27A pseudogene and RPL27A parent gene mRNA 3′ end uracil-rich terminal repeat (3′-UTR; the most common mRNA binding site for micro-RNAs) were predicted targets for miRNA 612, suggesting that the RPL27A pseudogene could serve as a miRNA decoy for miRNA 612, and both the parent gene NMRAL1 and the NMRAL1 pseudogenes were targets for miR 342–3p and miR 1913. The esiRNAs produced from the FBXO43 pseudogenes were not predicted to target FBXO43 transcripts but may target FBXO43-associated transcripts, including MAPK1, MDM2, OSMR, and IRAK3 transcripts (for each pseudoMap target score >200, minimum fold energy for RNA sequence ≥50).
Patterns of parent gene and pseudogene expression also suggested possible regulatory interactions. For example, the parent gene NMRAL1 was not differentially expressed in heart failure versus unused donor RVs (374 ± 132 vs. 351 ± 79 read counts). However, the NMRAL1 pseudogene was significantly reduced in heart failure versus unused donor RVs (2 ± 2 vs. 13 ± 15 read counts; Table 8). pseudoMap predicted miR 342–3p and miR 1913 as regulators of both the parent gene and pseudogene. If the NMRAL1 pseudogene served as an miRNA decoy, then reduction in NMRAL1 pseudogene expression without parallel reduction in NMRAL1 expression might permit miR 324–3p and/or miR 1913 repression of NMRAL1 translation. NMRAL1 (HSCARG) is a redox sensor protein. 22 Reduced levels of NMRAL1 by RNA interference increases nitric oxide production and reduces cell viability, whereas overexpression of NMRAL1 increases cell viability.
DISCUSSION
This study is notable for the high quality of isolated RNA, robust RNA-Seq count numbers, and abundant expression of both protein-coding and non-protein-coding poly(A) mRNAs. Transcripts in the failing RVs (as annotated in GENCODE, version 16, a tabulation of all known genes based on known human genome sequence) included 17,247 total genes, including 13,019 protein-coding genes (65% of all known protein-coding genes), 2,085 lncRNAs (16% of all known lncRNAs), and 1,064 pseudogenes (8% of all known pseudogenes). These findings are in general agreement with earlier RNA-Seq studies in mice, albeit with slightly larger numbers of transcripts.23–25 The heart failure RVs expressed a larger number of RNA transcripts in all categories, including lncRNAs and pseudogenes, compared with the unused donor RVs. Overall, 3.7% of lnc RNAs and 2.5% of pseudogenes were differentially expressed in the heart failure RV versus unused donor RVs. Exploratory in silico analyses using online tools suggested multiple roles for lncRNAs, particularly roles as (1) potential miRNA decoys and (2) high-affinity binding partners for chromatin modifiers, such the EZH2 component of PRC2.
lncRNA regulatory roles
lncRNAs play a variety of defined critical regulatory roles, including X-inactivation and imprinting via interaction with DNA and/or DNA binding proteins, 26 gene repression or enhancement via interaction with PRC1 and PRC2, 27 and gene expression enhancement during development and differentiation in part via interaction with the mediator complex. 28 Additional myriad lncRNA regulatory mechanisms reported to date include (1) transcriptional silencing or augmentation by interaction with transcriptional machinery or chromatin-altering binding partners, (2) interaction with actively transcribing mRNA, (3) direct interaction with DNA itself, (4) nuclear esiRNA generation via lncRNA “parent” fragmentation, (5) cytoplasmic miRNA decoy via the competing endogenous RNA hypothesis of miRNA regulation, and (6) interference with cytoplasmic mRNA stability by competition for mRNA stabilization factors. 29 This functional diversity devolves from the wide array of protein, RNA and DNA binding domains, flexible scaffolding capabilities, modular secondary structures (Fig. 2), embedded esiRNAs, alternative transcripts, and alternatively spliced isoforms characteristic of lncRNAs. 30
Three representative examples of long noncoding RNA (lncRNA) secondary structure. Examples of the diverse lncRNA RNAfold proposed secondary structures for three of the differentially expressed lncRNAs. Notable are the unique secondary conformations coupled with preserved modular features, including loops. The RNAfold proposed secondary structures were downloaded from http://www.lncipedia.org. 13
lncRNAs in cardiac disease: evidence to date
Although the role of lncRNAs in epigenetic regulation in development, differentiation, and cancer is well established 3 and will likely provide opportunities for novel targeted therapies, 31 the role of lncRNAs in cardiac development and disease has not been widely studied. 32 In a recent study, the lncRNA “Braveheart” was required for both progression of nascent mesoderm toward a cardiac fate and maintenance of cardiac fate in neonatal mouse cardiomyocytes in part via interaction with the SUZ12 component of PRC2. 33 lncRNA NAT transcripts 34 have been shown to regulate expression of troponin I, 35 myosin heavy and light chains, 36 and atrial natriuretic peptide. 20 A higher antisense/sense mRNA ALC-1 ratio was associated with lower ALC-1 protein expression in hypertrophied human myocardium. 37 The DMPK 3′-UTR CTG triplet repeat, when expanded 50–2,000 times, as in myotonic dystrophy, leads to “RNA toxicity” by sequestration of RNA-binding proteins, downregulation of connexins 40 and 43, induction of Nkx2–5 expression, and promulgation of heart block. 38 The SRA1 gene yields both (1) lncRNAs termed “steroid receptor RNA activators” (SRAs), which are able to coactivate steroid nuclear receptors, and (2) mRNAs translated into steroid receptor RNA activator proteins (SRAPs). 39 SRAPs may bind SRAs, and the SRAP/SRA transcripts serve as a “bi-faceted” genetic system with important roles in myogenic differentiation 40 and human dilated cardiomyopathies. 41 In separate genome-wide association studies, single-nucleotide polymorphism (SNP) variants in the novel lncRNA MIAT conferred increased risk of myocardial infarction, 42 and SNPs in the lncRNA ANRIL locus carried the strongest genetic susceptibility for coronary artery disease. 43 As more mechanistic examples, the lncRNA ANRIL, a natural antisense transcript of the protein-coding INK4b/ARF/INK4a complex, binds chromobox 7 of PRC1, thereby regulating the PRC1-mediated methylation of histone H3 lysine 27 and effecting the transcriptional silencing of INK4a in cis. 44 The lncRNA cardiac hypertrophy related factor in mice acts as an endogenous “sponge” for mi489, thereby relieving the mi489-mediated repression of hypertrophy-inducing Myd88. 45 esiRNA silencing of the lncRNA MALAT1, which under normal conditions is highly expressed in endothelial cells of different origin, resulted in endothelial-dependent reduced vascular growth both in vitro and in vivo. 46
Speculative roles of lncRNAs in the heart
Little is presently known regarding the targets of the lncRNAs expressed in the heart. 47 We could not find published literature identifying the targets of the differentially expressed lncRNAs identified in this study. Thus, we used available online tools for lncRNA in silico preliminary functional analyses. Future efforts, including our own, will be directed toward elucidating lncRNA targets using RNA silencing (esiRNA) approaches in animal and cellular models. Elucidation of the roles of specific lncRNAs in the heart will prove challenging given the large number of lncRNAs, the large number of potential lncRNA targets and interaction partners, and the known mechanistic diversity of lncRNAs themselves. 47 Four broad functional roles as signals, decoys, guide, and scaffolds have been proposed for lncRNAs (Fig. 3). 48 First, lncRNAs may serve as signals themselves once transcribed and/or modulate existing signals by interacting with chromatin-binding protein complexes or other RNA species. 49 Second, lcnRNAs may serve as decoys for both RNAs and proteins. The “endogenous competing RNA hypothesis” proposes that multiple RNA species (mRNAs, lncRNAs, small nucleolar RNAs, and miRNAs) all “compete” for specific regulatory miRNA binding partners in a process of combinatorial regulation based on relative binding affinities and expression levels. 50 Third, lncRNAs often possess high-affinity binding sites for proteins and may serve as decoys. By binding to DNA and forming heteroduplexes or heterotriplexes, lcnRNAs may also guide chromatin-binding proteinstoregulatory cis and transgenomic regions. 48 Fourth, lncRNAs, due to their modular structure, may serve as scaffolds for protein complex and RNA-binding partner and assist or direct the assembly of larger chromatin regulatory complexes. 51 lncRNAs thus likely play key roles in the epigenetic regulation of gene expression during adaptive and maladaptive ventricular remodeling and the transition to failure.47,49
Domain architecture and function of long noncoding RNAs (lncRNAs). A, lncRNAs structural domains can sense or bind other RNAs via complementary base pair interactions, proteins, and DNA that can induce allosteric conformations changes to other structures in the lncRNA. B, Alternative splicing can combine these structural domains into lncRNA modular architecture. C, Coupling sensing and actuator domains permit lncRNAs to function as “control” devices. In the example on the left, binding of RNA (gray) induced a conformational change that prevents protein binding. On the right, the protein can bind in the absence of RNA, inducing the formation of a stem-loop secondary structure that can be processed and cleaved to generate RNA output. D, lncRNAs such as HOTAIR act as molecular scaffolds, binding multiple proteins to form complex ribonucleoprotein structures. E, lncRNAs such as Cist can target the catalytic function of proteins to specific genome sites. lncRNAs can recruit chromatin-modifying proteins (purple, blue) to target sites by association with a DNA-binding protein such as YY1 (red). The chromatin modifiers then modify local histones to influence the expression of adjacent genes.
lncRNAs in the RV
Animal models have demonstrated significant differences in RV mRNA expression transcriptional profiles in response to experimental pressure52–56 and volume overload.57,58 Most of these studies, however, examined only the expression of mRNA by earlier-generation microarrays, which did not include probes for lncRNAs and thus provided an incomplete readout of the comprehensive RV transcriptome.5,6 lncRNA differential expression by RNA-Seq in RV pressure versus volume overload experimental models has not been reported to date. Based upon the dynamic nature of RV mRNA expression in pressure versus volume overload models, however, it is likely that lncRNA expression differences exist in such models and play key regulatory roles in the epigenetic regulation of the differential mRNA transcription profiles reported. Differential expression of lncRNAs and associated regulatory pathways is likely a conserved feature of RV remodeling independent of the type of RV remodeling stimulus. We would anticipate that at least some of the differentially expressed lncRNAs and associated pathways identified in this study would play roles in RV adaptive and maladaptive remodeling and the transition to RV failure in pulmonary arterial hypertension.
lncRNA in the RV versus LV
In the mouse, the normal RV and normal LV are epigenetically distinguishable based on differential expression patterns of epigenetic regulators of histone acetylation and methylation. 59 Differences in RV versus LV mRNA transcription have been reported in the normal mammalian heart, 60 in some61,62 but not all 63 pressure overload experimental models, and in human congenital heart disease. 64 Although the RV versus LV may invoke similar regulatory pathways in response to pressure overload, the RV mRNA transcriptional response appears generally more vigorous than the LV response. 65 Based upon the aggregate microarray data to date, however, it is not yet clear whether targeting pathways specific to the stressed or failing RV, targeting pathways common to both the RV and LV, or targeting the mechanical interactions between the RV and LV will afford novel opportunities to attenuate RV adverse remodeling and the transition to RV failure.65,66
Pathobiology of RV failure in LV failure and future potential applications of lncRNA transcriptional profiling
The pathobiology of RV failure in the setting of LV failure is likely distinct based upon the etiology of LV failure. 65 Given the relative resistance of the RV to irreversible ischemic injury, RV failure in ischemic LV failure most likely devolves from the chronic elevation of pulmonary vascular impedance resulting from progressive LV diastolic hypertension. 67 In impedance overload, the RV remodels by concentric hypertrophy and adopts over time with unrelieved overload many of the geometric and cardiodynamic features of the LV. In nonischemic LV failure, by contrast, RV failure may result from the same heritable structural or contractile protein mutations or infectious/toxic injuries affecting the LV. 67 Heritable protein mutations or injuries delimit the RV remodeling repertoire to any additional insult or overload imposed by LV failure. Hence, RV failure pathobiology due to heritable protein mutations or injuries diverges sharply from the RV failure pathobiology of pulmonary vascular impedance overload alone. Beyond these broad etiologic factors, other modifying genetic and environmental factors—including common comorbid conditions that alter pulmonary hemodynamics independent of the heart (e.g., obstructive sleep apnea and chronic obstructive pulmonary disease), general medical comorbidities (e.g., diabetes mellitus), and specific heart failure pharmacologic therapies—introduce higher orders of heterogeneity into the pathobiology of RV failure in the setting of LV failure.65,66 To some degree, the RV transcriptome in end-stage LV failure would be expected to reflect the different etiologies of LV failure, the variable pulmonary vascular impedance conditions dependent on LV function, and the integration of these and other diverse genetic and environmental modifying factors. Collectively, these pathobiologic considerations in end-stage human LV heart failure suggest that (1) the RV transcriptome is not likely to mirror directly the LV transcriptome and (2) the RV transcriptome may exhibit substantial patient-to-patient heterogeneity. Such RV transcriptome heterogeneity was supported by the striking variability in subject-to-subject expression of many differentially expressed lncRNAs identified in this study and may afford a potential future opportunity to personalize therapies for RV failure.
Beyond providing insight into fundamental mechanisms of RV myocardial development and remodeling as well as affording an opportunity to personalize therapies, lncRNA transcriptional profiling may yield diagnostic clinical assays and prognostic biomarkers in human heart failure. A recent study of explanted human LVs reported that the LV lncRNA transcriptional profile discriminated the etiology of ischemic versus nonischemic heart failure with greater precision than did the LV mRNA or miRNA transcriptional profiling. 68 This study also showed that lncRNA expression profiles, but not mRNA or miRNAs, were markedly altered in response to the reverse remodeling observed with LV mechanical support. 68 lncRNA profiling and therapeutic targeting may complement the evolving roles of miRNA profiling and therapeutic targeting in the LV 69 and in RV dysfunction. 70 In a recent study, the mitochondrial lncRNA LIPCAR was upregulated in the plasma of patients who developed adverse ventricular remodeling and heart failure after myocardial infarction and was shown to be an independent predictor of survival in patients with heart failure. 71
Findings of additional speculative interest
Several of the findings in the present study are of additional speculative interest. First, we found expression of approximately 3,000 lncRNAs (including pseudogenes) representing 16% of known lncRNAs and 8% of known pseudogenes in RV myocardium. Approximately 6% overall of these lncRNAs were differentially expressed in heart failure versus unused donor RVs. These percentages of overall lncRNA expression and differential lncRNA are both modest. Most lncRNA expression, like mRNA expression, appears conserved in the failing RV, which is perhaps not surprising given the dependence of working myocardium (failing or not) on the highly expressed, highly conserved mitochondrial energetic and sarcomeric contractile protein components that dominate the myocardial transcriptome. 68 Second, there were a larger number of transcripts expressed overall in the heart failure versus unused donor RVs. This may have arisen from a variety of factors. The heart failure subjects had a greater number of comorbidities and medications, both of which may have induced different and richer expression patterns. Although the explantation techniques, cardio-preservation protocols, and tissue-preparation processes were identical in the unused donor and heart failure subjects, the unused donors and heart failure subjects were in markedly different physiologic states at the time of explanation. The unused donor subjects were brain dead, hormonally resuscitated, anesthetized without induction, were not supported on cardiopulmonary bypass, and had multiple vital organs ex-planted before cardiac explantation. Conversely, the heart failure subjects had intact brain function, were anesthetized with careful induction, and were supported by cardiopulmonary bypass with fully intact vital organ function at the time of explantation. Collectively, these differences in physiologic state may have contributed to the gene expression differences found. 72
Limitations
Our study was limited to explanted end-stage human myocardium. Although the RNA submitted for sequencing from the homogenized myocardial tissue was of high quality, we cannot exclude selective enrichment of RNAs due to sample-to-sample differences in tissue characteristics or composition, myocardial perfusion, or blood admixture. Although our tissue collection and preservation protocol was designed to minimize RNA degradation, the sequential processes of harvesting, transport, frozen storage, thawing, tissue homogenization, and RNA extraction cannot preserve all tissue RNA.
By design, we analyzed only RV tissues. Thus, we cannot compare differential expression of LV lncRNA versus RV lncRNA. The earlier study of LV lncRNA differential expression, cited above, did not provide a tabulation of differentially expressed LV lncRNAs to permit comparison with the tabulation of differentially expressed RV lcnRNAs in our study. 68 All patients were treated with medications, including inotropes, before transplantation, which likely altered gene expression in the RV. RNA-Seq has inherent potential limitations owing to 5′ bias and dependence on library-preparation methodology. 7 As demonstrated, substantial biological variation was found in read counts between subjects. The pooling of samples with such read count variation challenges statistical inference based on an assumed distribution of read counts. Statistical significance, as in all such studies, is accepted as a “proxy” for biological significance, realizing that meaningful biological variation may or may not result in statistical differences in expression if there are threshold, graduated, or nonlinear expression effects on translation.
We did not perform confirmatory qPCR or other assays for the differentially expressed genes reported. Multiple earlier studies have demonstrated equivalent precision of RNA-Seq and qPCR.5,6,73,74 Given the limitations of coverage, copy readout, and specificity of microarray probes, confirmatory qPCR is required to validate microarray findings. RNA-Seq, however, is not limited by coverage, copy readout, and specificity, and confirmatory qPCR is not currently recommended by experts for well-performed RNA-Seq.5,6,73,74 For this study, RNA-Seq was performed in the Vanderbilt Genomics Core laboratory, a facility with highly trained and experienced personnel and state-of-the-art equipment.
We compared RVs from subjects with end-stage LV failure with the RVs of unused human donors. Unused human donors do not have “normal” RV myocardium because of the state of brain death before harvest. Earlier studies have reported heterogeneity in gene expression in unused donor hearts attributable to clinical factors and harvest/storage factors. 72 Finally, all such studies are inherently limited by the lack of understanding of fundamental aspects of gene expression, including the relationship of transcribed protein-coding or non–protein-coding mRNA to translated protein. 75
Conclusion
RNA-Seq providesa comprehensive profile of the polyadenylated lncRNA transcriptome in human heart failure. Over 100 polyadenylated lncRNAs are differentially expressed in heart failure RVs versus unused donor RVs. Differentially expressed lncRNAs, including natural antisense transcripts and pseudogenes, likely possess a variety of regulatory functions.
Despite the likely critical roles of lncRNAs in myocardial pathobiology, myocardial lncRNA research is in its infancy. The current study was descriptive and exploratory by design. It describes a reasonably large number of differentially expressed lncRNAs in human heart failure and explores their possible functional roles via accessible, existing online resources and databases. It is hoped that these results provide investigators with an inventory of differentially expressed candidate lncRNAs, natural antisense transcripts, and pseudogenes for future mechanistic studies in human heart failure.