Downregulation of AC 061961 . 2 , LINGO 1-AS 1 , and RP 11-13 E 1 . 5 is associated with dilated cardiomyopathy progression

Purpose: This study aimed to explore lncRNAs implicated in dilated cardiomyopathy (DCM) using high-throughput sequencing. Methods: From February 2016 to May 2017, ten samples of failing hearts collected from the left ventricles of DCM patients undergoing heart transplants, and ten control samples obtained from normal heart donors were included in this study. After sequencing, raw data were processed, and differentially expressed genes (DEGs) and lncRNAs between samples from patients with DCM and controls were screened, followed by functional enrichment analysis and Weighted Gene Coexpression Network Analysis (WGCNA). Five key lncNRAs were validated through real-time PCR. Results: After processing sequence data, 1398 DEGs were identified between DCM and control groups, including 267 lncRNAs. WGCNA identified seven modules that were significantly correlated with DCM. The top 50 genes in the three modules (black, dark-green, and green-yellow) were significantly correlated with the DCM disease state. Four core enrichment lncRNAs, such as AC061961.2, LINGO1-AS1, and RP11-557H15.4, in the green-yellow module were associated with functions of neurotransmitter secretion. Five core enrichment lncRNAs, such as KB-1299A7.2 and RP11-13E1.5, in the black module co-regulated functions associated with the regulation of blood circulation and heart contraction. AC061961.2, LINGO1-AS1, and RP11-13E1.5 were confirmed to be downregulated in DCM tissues by real-time PCR. Conclusion: The present study suggests that downregulation of AC061961.2, LINGO1-AS1, and RP11-13E1.5 is associated with DCM progression, and these may serve as key diagnostic biomarkers and therapeutic targets for DCM.


INTRODUCTION
Cardiomyopathies are a group of heart muscle diseases characterized by abnormal chamber size and wall thickness or abnormal contractile function [1].Dilated cardiomyopathy (DCM) is the most common cardiomyopathy worldwide, characterized by enlargement in left ventricular cavity, as well as impaired contractility in the absence of significant coronary artery disease [2].Additionally, some abnormal loading conditions, including primary valvular diseases, congenital heart diseases and arterial hypertension may also result in DCM.DCM may lead to sudden cardiac death and progressive heart failure and thus, is related to serious morbidity and premature mortality [3].Despite therapeutic advances, the 5-year mortality for DCM remains up to 20% [4].Studies have found that genetic etiologies are associated with the pathogenesis of DCM [5,6].Therefore, understanding the genetic basis of DCM and identifying causative genes in individuals with DCM may help in improving the clinical management of this disease.
Lakdawala et al. [6] recently studied a diverse population of patients with DCM for clinical genetic testing.They identified relevant genetic mutations in 17.4% of those DCM patients, commonly involving the myosin heavy chain 7, lamin A/C, and troponin T2 genes.The completion of the Human Genome Project revealed that 98% of the genome consists of non-coding sequences, including microRNAs and long non-coding RNAs (lncRNAs).LncRNAs are transcribed RNA molecules longer than 200 nucleotides that actively participate in various physiological and pathological processes [7].Several recent studies have identified many lncRNAs associated with diseased hearts [8][9][10].However, there are still not enough studies on the functions of lncRNAs in DCM [11,12].
In this study, we aimed to explore more lncRNAs implicated in DCM using high-throughput sequencing of lncRNAs in myocardial tissues collected from DCM patients and normal heart donors.Based on sequencing data, we screened differentially expressed genes (DEGs) and lncRNAs (DELs), followed by a series of bioinformatic analyses for identifying key lncRNAs involved in DCM.The findings were then validated through RT-PCR.

Raw read quality control and reference genome alignment
The number of raw reads in paired ends of each sample was more than 30,000,000 (data not shown).After filtration, numbers of clean reads ranged from 28,248,054 to 35,511,833 in each sample.The rates of clean reads were higher than 86.2%.The rates of read alignment to the human reference genome (hg19) were between 83.1% and 85.8% in the DCM group and between 69.3% and 85.3% in the control group (data not shown).

Differential gene and lncRNA analysis
A total of 21,065 genes were expressed in all samples of at least one group (DCM or control).Among these genes, 4328 were lncRNAs.With thresholds of |log (fold change)| > 1 and adjusted p-value < 0.05, 1398 DEGs (447 upregulated and 951 downregulated ones) were identified.Additionally, these DEGs included 267 DELs (93 upregulated and 174 downregulated ones).Heat maps of DEGs and DELs are shown in Figure 1.DEGs and DELs could clearly distinguish between the DCM and control groups.

Functional enrichment analysis
Upregulated DEGs were significantly associated with biological processes linked to heart and muscle contraction as well as the pathway of cardiac muscle contraction (hsa04260) (Figure 2A).Downregulated DEGs were significantly connected to the following GO terms: regulation of membrane potential (GO:0042391), synapse organization (GO:0050808), and signal release (GO:0023061) as well as pathways of protein digestion and absorption (hsa04974), circadian entrainment (hsa04713), and calcium signaling pathway (hsa04020) (Figure 2B).

PPI network and subnetwork module analysis
The PPI network consisted of 3565 edges and 799 nodes (252 upregulated and 547 downregulated genes).The top 10 gene nodes with high degrees included ALB, TNF, SRC, JUN, RIPK4, FOS, MYC, ACTA1, ACTA2, and EDN1.In addition, 21 modules were extracted from the PPI network.The two modules containing the most nodes are shown in Figure 3A and 3B.Module 1 contained 35 nodes and 287 edges, while module 2 included 36 nodes and 279 edges.DEGs in module 1 were significantly involved in neuroactive ligand-receptor interaction, phospholipase D signaling pathway, and calcium signaling pathway and DEGs in module 2 were particularly associated with the AGE-RAGE, TNF, and IL-17 signaling pathways (Figure 3C).

WGCNA analysis
DEGs were subjected to WGCNA module analysis, and each module contained at least 10 genes.A total of 13 modules were identified.Correlations between WGCNA modules and DCM are shown in Table 1.Upon applying the thresholds of p-value < 0.01 and correlation coefficient > 0.7, seven modules (indicated by different colors) were correlated with DCM, among which the magenta, blue, and purple ones were positively correlated with DCM, while the green-yellow, salon, dark-green, and black ones were negatively correlated with it.In this study, we focused on the top 50 genes in these seven modules.LncRNAs were included among the top 50 genes in each module.

Correlation analysis of WGCNA modules and the DCM disease state
Using GSAASeqSP, correlations between the top 50 genes in the seven modules and the DCM disease state were analyzed.The correlation p-values between three modules (black, dark-green, and green-yellow) and the DCM disease state were less than 0.05.These three modules are shown in Figure 4.After GASS analysis, the core enrichment genes in these three modules were obtained.The black module had 35 core enrichment genes, which included five lncRNAs (AC009878.2,CTB-174D11.1/KB-1299A7.2,RP11-557H15.4,and RP4-534N18.2); the dark-green module had 12 core enrichment genes without any lncRNAs; and the green-yellow module had 31 core enrichment genes, including four lncRNAs (AC061961.2,LINGO1-AS1, RP11-13E1.5, and RP11-218M11.3).

lncRNA regulatory network analysis
The mRNAs involved in the three above-mentioned modules were regarded as lncRNA correlated genes.Based on the nine lncRNAs identified in the three modules, we constructed lncRNA regulatory networks, as shown in Figure 5.

Functional enrichment analysis of lncRNAs
These lncRNA correlated genes were subjected to functional enrichment analysis.As shown in Figure 6, the five core enrichment lncRNAs in the black module co-regulated functions associated with regulation of blood circulation, heart contraction, membrane potential, heart process, heart contraction, and muscle contraction.Meanwhile, in the green-yellow module, the four core enrichment lncRNAs co-regulated functions of synaptic vesicle exocytosis, regulation of calcium ion-dependent

Real-time PCR verification of the expression of key lncRNAs
Expression levels of the top five key lncRNAs with low p-values, namely, AC061961.2,LINGO1-AS1, RP11-557H15.4,KB-1299A7.2,and RP11-13E1.5,were determined using real-time PCR.As shown in Figure 7, AC061961.2,LINGO1-AS1, and RP11-13E1.5 were     significantly downregulated in samples from DCM patients compared with that in control samples (p < 0.01).RP11-557H15.4was downregulated in the DCM group, albeit without a significant difference.However, KB-1299A7.2was not detected in samples from patients with DCM, which may be due to its expression level being too low in them.

DISCUSSION
In the present study, we obtained 1398 DEGs from the comparison between DCM and control groups; these DEGs included 267 DELs.WGCNA of these DEGs identified seven modules significantly correlated with DCM.Moreover, the top 50 genes in the black, dark-green, and green-yellow modules were significantly correlated with the DCM disease state.Additionally, nine core enrichment lncRNAs, including AC061961.2,LINGO1-AS1, RP11-557H15.4,KB-1299A7.2,and RP11-13E1.5,were obtained from the black and green-yellow modules.Expression levels of AC061961.2,LINGO1-AS1, and RP11-13E1.5 in tissue samples from patients with DCM were verified through real-time PCR.
these lncRNAs were not identified in our study, which may be due to the sample heterogeneity.Cardiac neurotransmitter systems are impaired in heart diseases.In patients with heart failure, abnormalities in the cardiac neurotransmitter systems contribute to arrhythmogenesis and cardiac dysfunction [17].Qin et al. [18] also found that the heightened sympathetic state during congestive heart failure was associated with decreased cardiac sympathetic neurotransmitters.The role of the cardiac neurotransmitter system in DCM has not been fully investigated.In the present study, AC061961.2,LINGO1-AS1, and RP11-13E1.5 were significantly associated with the functions of neurotransmitter secretion and transport.Additionally, they were downregulated in samples from patients with DCM, which suggest their roles in DCM progression.
LncRNA regulatory networks were constructed based on WGCNA modules, and SOX5 was predicted to be a correlated gene of AC061961.2,LINGO1-AS1, and RP11-13E1.5.SOX5 belongs to the SOX family of transcription factors, which is expressed in many human tissues, including the heart [19,20].Genome-wide association studies have implicated SOX5 as a candidate gene for susceptibility to the electrocardiographic PR interval [21], higher resting heart rate [22], left ventricular mass [23], and atrial fibrillation [24], suggesting its importance in the heart.A recent study also revealed that SOX5 regulated cardiac function by negatively regulating Wnt signaling, which had been implicated as an important event regulating cardiac function and development [25].Taking these findings together, these three lncRNAs may be involved in the progression of DCM by interacting with SOX5.
The correlated genes of lncRNAs in the black module co-regulated functions associated with the regulation of blood circulation, heart contraction, heart process, and membrane potential.For instance, ADORA1 was a correlated gene of RP11-557H15.4 and KB-1299A7.2,which was also involved in the PPI subnetwork module 1. ADORA1 was significantly enriched in the above GO functions.ADORA1 is a member of the adenosine receptor family of G-proteincoupled receptors with adenosine as an endogenous ligand [26].Adenosine is implicated in cardioprotective events, particularly in cardiac hypertrophy [27].For example, Xu et al. [28] reported that reduced adenosine levels results in

GAPDH-hR AGGCAGGGATGATGTTCTGGAGAG
obvious cardiac hypertrophy in rat model.A recent study also suggested that ADORA1 prevents transition from compensated myocardial hypertrophy to decompensated heart failure [29].Therefore, RP11-557H15.4 and KB-1299A7.2may involve in DCM pregression through the downregulation of ADORA1.Establishment of the relationship among diseases, physiological processes, and the action of smallmolecule therapeutics is a fundamental challenge arising throughout biomedicine.Connectivity maps are used to search for connections among small molecules sharing a mechanism of action, physiological processes and chemicals, and diseases and drugs [30].In this study, we used a connectivity map for predicting potential drug small molecules for treating DCM.Many small molecules, such as N-[4-(3-oxo-3-phenylprop-1-en-1-yl)phenyl]-2,2-diphenylacetamide, 4,5-dianilinophthalimide, and 1,4-chrysenequinone, were predicted, but their roles in DCM require further validation.
In conclusion, the present study suggests that AC061961.2,LINGO1-AS1, and RP11-13E1.5 are downregulated in myocardial tissues of patients with DCM.Additionally, downregulation of RP11-557H15.4 and KB-1299A7.2may also be associated with DCM, although this requires further confirmation.These lncRNAs may serve as key diagnostic biomarkers and therapeutic targets for DCM.

Samples
Ten human failing heart tissue samples were collected from the left ventricles of DCM patients who accepted heart transplants.These patients aged 38-59 (average age 47.4 ± 2.4) years and were enrolled from February 2016 to May 2017.The inclusion criterion was that they should be idiopathic DCM.The patients with secondary cardiomyopathy were excluded.Additionally, ten healthy heart samples (as control) were obtained from the left ventricles of normal heart donors (average age 51.2 ± 2.3 years) that were unsuitable to transplant for several technical reasons.Detailed patient information is shown in Table 3.All patients provided informed consent before the study.All procedures in this study were approved by our hospital's protection of human ethics committee.

Total RNA isolation and quality control
Total RNAs were isolated from 10 myocardial tissue samples (5 samples from DCM group and 5 from control group) using a TRIzol reagent (Invitrogen, CA, USA).RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA).NanoPhotometer ® spectrophotometer (IMPLEN, CA, USA) was used to check RNA purity, and Qubit ® RNA Assay Kit in Qubit ® 2.0 Flurometer (Life Technologies, CA, USA) was used to measure the RNA concentration.

Library preparation
For each sample, 3 μg of RNA was used for the preparation of RNA samples.Ribosomal RNA was removed using the Epicentre Ribo-zeroTMrRNA Removal Kit (Epicentre, USA).rRNA-free residue was removed by ethanol precipitation.Sequencing libraries were then generated with the NEBNext ® Ultra TM Directional RNA Library Prep Kit for Illumina ® (NEB, USA).Briefly, fragmentation was performed using divalent cations in NEBNext First Strand Synthesis Reaction Buffer (5×).First-strand cDNA was synthesized using M-MuLV Reverse Transcriptase and random hexamer primers, and second-strand cDNA was synthesized using RNase H and DNA Polymerase I.In the reaction buffer, dTTP of dNTPs was replaced by dUTP.After adenylation of the 3′ ends of the DNA fragments, NEBNext Adaptor with a hairpin loop structure was ligated to prepare for hybridization.To select cDNA fragments preferentially, 250-300 bp in length, the library fragments were purified with the AMPure XP system (Beckman Coulter, MA, USA).Then, 3 μL of USER Enzyme (NEB, USA) was used with sizeselected, and adaptor-ligated cDNA at 37°C for 15 min followed by 5 min at 95°C.Next, PCR was performed based on Universal PCR primers, Index (X) Primer, and Phusion High-Fidelity DNA polymerase.The products were purified (AMPure XP system), and the library quality was assessed on the Agilent Bioanalyzer 2100 system.

Clustering and sequencing
Index-coded samples were clustered on a cBot Cluster Generation System using HiSeq4000 PE Cluster Kit (Illumina).Then libraries were sequenced on an Illumina Hiseq4000 platform and 150-bp paired-end reads were generated.Obtained sequencing data were deposited in the NCBI SRA database (http://www.ncbi.nlm.nih.gov/sra) under series accession number SRP108571.

Raw read quality control and reference genome alignment
Obtained raw reads were quality-filtered to remove low-quality reads and adaptor sequences using Fastx toolkit version 0.0.13 and Prinseq-lite version 0.20.4[31].Obtained clean data were aligned to the human reference genome (hg19) using TopHat version 2.0.8 [32].

Quantitative and differential gene and lncRNA analyses
Based on the gene and lncRNA annotation information in the GENCODE version 24 (mapped to GRCh37) database [33], we obtained the fpkm value and read count of genes using the StringTie version 1.2.2 [34].DEGs including DELs were identified between the disease and normal control groups using the edgeR [35] package in R by applying the following criteria: |log (fold change)| > 1 and adjusted p-value < 0.05.

Prediction of lncRNA correlated genes
Weighted Gene Coexpression Network Analysis (WGCNA) is used for identifying clusters (modules) of highly correlated genes [36].In this study, all DEGs including DELs were subjected to WGCNA using the WGCNA package [36] in R. mRNAs in the module that contained DELs were regarded as correlated genes of that lncRNA.

Functional enrichment analysis
DEGs and DELs correlated genes were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses using clusterprofiler [37] in R.

Correlation analysis of WGCNA modules and DCM disease state
The correlation between WGCNA modules and DCM was analyzed using one-way ANOVA and the module with a p-value of < 0.05 was selected.Additionally, based on the expression values of module genes, we analyzed the correlation between module genes and the DCM disease state using a toolset of Gene Set Association Analysis for RNA-Seq with Sample Permutation (GSAASeqSP, http://gsaa.unc.edu)[38], and the module with a nominal p-value of <0.05 was considered to be significantly correlated with the DCM disease state.Furthermore, according to the core enrichment gene obtained from GSAASeqSP, we screened key lncRNAs and constructed the lncRNA regulatory network on the basis of the WGCNA module.The lncRNA regulatory network was then visualized through Cytoscape [39].

Protein-protein interaction (PPI) network analysis
PPI networks of DEGs and DELs correlated genes were predicted using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING version 10.0, http:// string-db.org)[40] and were visualized using Cytoscape [39].In addition, PPI subnetwork modules were identified using the MCODE [41] plugin in Cytoscape.Parameters of MCODE were set as follows: degree cut-off = 5, node score cut-off = 0.2, and K-Core = 2.

Prediction of potential drug small molecules for treating DCM
For the selected key lncRNAs in WGCNA modules, we predicted their potential drug molecules from the connectivity map [42] database according to their expression levels.

Real-time PCR verification of the expression of key lncRNAs
A total of 20 samples, including 10 samples from patients with DCM and 10 from controls, were used for real-time PCR verification of the expression of key lncRNAs.Briefly, total RNAs were isolated using a TRIzol reagent (Invitrogen, CA, USA).RNA concentration and quality were determined on a TECAN infinite M100 PRO Biotek microplate reader (TECAN, CA, USA).Then 0.5 μg of the total RNA was used from cDNA synthesis using the PrimeScript RT Master Mix (RR036A; Takara, Dalian, China).PCR was performed using the SYBR GREEN kit (4367659; Thermo, USA) in the Viia7 Real-Time PCR System (Applied Biosystems, USA).The primers used in this study are listed in Table 4.

Statistical analysis
Data are presented as mean ± standard deviation.Statistical analysis was performed using SPSS 22.0 (IBM, Armonk, NY, USA).Differences in miRNA expression levels among different groups were analyzed by one-way analysis of variance and were considered significant if p < 0.05.

Figure 1 :
Figure 1: The heat maps of (A) differentially expressed genes and (B) differentially expressed lncRNAs between dilated cardiomyopathy (DCM) and control samples.

Figure 2 :
Figure 2: Gene Ontology (GO) functions and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched by (A) upregualted and (B) downregulated differentially expressed genes.The x-axis represents GO and pathway terms, and the y-axis represents the number of genes enriched in terms.BP: biological process; CC: cellular component; MF: molecular function.

Figure 3 :
Figure 3: The top 2 protein-protein interaction modules containing more nodes.(A) Module 1 contains 35 nodes and 287 edges.(B) Module 2 includes 36 nodes and 279 edges.Red hexagon represents upregualted genes and green circle represents downregualted genes.(C) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched by differentially expressed genes in module 1 and module 2.

Figure 5 :
Figure 5: The lncRNA regulatory networks based on the nine lncRNAs and their correlated genes.Red hexagon represents upregualted genes; green circle represents downregualted genes; red triangle represents upregualted lncRNA; green triangle represents downregualted lncRNA.

Figure 6 :
Figure 6: The Gene Ontology (GO) functions of (A) biological process, (B) cellular component and (C) molecular function enriched by the correlated genes of the nine lncRNAs.