Evaluation of lncRNA Expression Pattern and Potential Role in Heart Failure Pathology

During the last few decades, the morbidity and mortality of heart failure (HF) have remained on an upward trend. Despite the advances in therapeutic and diagnostic measures, there are still many aspects requiring further research. This study is aimed at finding potential long noncoding RNAs (lncRNAs) that could aid with the diagnosis and treatment of HF. We performed RNA sequencing on the peripheral blood of healthy controls as well as HF patients. The expression of lncRNAs was validated by RT-qPCR. Bioinformatic analysis was performed to investigate the possible mechanism of differentially expressed lncRNAs and mRNAs. The diagnostic value of lncRNAs was analysed by ROC analysis. Finally, a total of 207 mRNAs and 422 lncRNAs were identified. GO and KEGG pathway analyses revealed that biological pathways such as immune response, regulation of cell membrane, and transcriptional regulatory process were associated with the pathological progress of HF. The lncRNA-mRNA coexpression network was conducted, and several mRNAs were identified as key potential pathological targets, while lncRNA CHST11, MIR29B2CHG, CR381653.1, and FP236383.2 presented a potential diagnostic value for HF. These findings provide novel insights for the underlying mechanisms and possible therapeutic targets for HF.


Introduction
At present, statistics show that cardiovascular diseases (CVDs) are the leading cause of death around the world. Among them, it is estimated that about 40 million people suffer from heart failure (HF) in the whole world, and this number continues to grow due to the rapidly aging population as well as the global decline in birthrate [1,2]. In addition, the aggravating changes in people's lifestyles, sedentarism and the high-sugar and high-fat diets, are starting to show an increase in the incidence of HF in the younger populations, particularly in developing countries. In older individuals, HF tends to occur at the final stage of a series of CVDs, such as in the case of myocar-dial infarction and hypertension [3]. Left ventricular ejection fraction (LVEF) has always been regarded as its pathological phenotype, but it cannot be widely used due to its inconvenience [4,5]. Although there are novel and more effective measures to intervene in HF, the overall prognosis remains not optimistic [6,7]. Therefore, a deeper understanding of the pathological mechanism behind HF is critical for effective diagnosis, treatment, and prognosis.
Long noncoding RNA (lncRNA) is a class of noncoding RNA which does not encode any proteins but plays a role in regulating mRNA transcription, protein location, and other cellular biological process [8,9]. Previous studies revealed the role of dysregulated lncRNAs in the development and progression of CVDs [10][11][12][13]. Overexpression of the lncRNA cardiac physiological hypertrophy-associated regulator (CPhar) could prevent ischemia-reperfusion-induced apoptosis and cardiac dysfunction [14]. LncRNA noncoding repressor of nuclear factor of activated T cells (NRON) and myosin heavy-chain-associated RNA transcripts (MHRT) were both overexpressed in the plasma of HF patients [15]. The exploration of the correlation between lncRNAs and CVD has increased our understanding of CVD pathology and provided new therapeutic targets for CVD.
In this study, we performed RNA sequencing on the peripheral blood collected from 100 volunteers (50 HF patients and 50 healthy controls). RT-qPCR was used to identify the differentially expressed genes. Functional enrichment analysis was used on both the GO and KEGG databases. We conducted the lncRNA-mRNA coexpression network and reported several potential target genes by cisand transregulation prediction. We then hypothesized that the mRNA-lncRNA expression profiles can provide important insights to understand the underlying mechanism of HF, and the differentially expressed lncRNAs may facilitate the finding of new diagnostic targets as well as aid drug development.

Subjects and Blood
Sample. In the present study, peripheral whole blood samples were collected from a total of fifty HF patients and fifty normal controls at the Affiliated Hospital of Qingdao University between December 2018 and January 2020. The selection criteria for HF patients were as follows: (1) The patient must have been admitted to the hospital without any history of related disorders or treatments prior to being diagnosed; (2) NT-proBNP (N-teminal pro-B type natriuretic peptide)≥400 ng/L [16]; and (3) LVEF <50% [16,17]. For the control group, we selected healthy individuals (in a similar age group) without a history of heart-related diseases and with normal echocardiograms. Among all participants, a total of four HF patients as well as healthy controls were selected for sequencing. This study was approved by the ethics committee of the Affiliated Hospital of Qingdao University. Each volunteer agreed and signed an informed consent form. The detailed workflow of this study is shown in Figure 1.

RNA Isolation.
Total RNA was extracted from the blood samples using RNA isolate Total RNA Extraction Reagent (Vazyme, China) according to the manufacturer's instructions. In brief, 100 μl of blood was combined with 400 μl of RNA extraction reagent. After mixing, add a fifth volume of chloroform. The mix was shaken vigorously by hand for 15 s and incubated for 3 min at room temperature. Centrifuged the mix samples at 12000 g for 15 min at 4°C. Removed the aqueous phase of the sample and added the equal volume of 100% isopropanol, incubating for 10 min at room temperature. The mixed samples were centrifuged at 12000 g for 10 min at 4°C. And then, the RNA pellet was washed two times with 75% ethanol, and the dried RNA pellet was resuspended in RNase-free water. RNA quality was verified using a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and the ND-2000 (NanoDrop Technologies). RNA samples (OD260/280 = 1:8 − 2:2, OD260/ 230 ≥ 2:0, RIN ≥ 8, 28S : 18S ≥ 1:0, and total RNA > 10 μg) were used to construct a sequencing library.   The raw paired-end reads were trimmed, and the quality was controlled by SeqPrep (https://github.com/jstjohn/ SeqPrep) and Sickle (https://github.com/najoshi/sickle) using the default parameters. Then, the clean reads were separately aligned to the reference genome with orientation mode using HIASAT software (https://ccb.jhu.edu/software/hisat2/index .shtml) [18]. The mapped reads of each sample were assembled by StringTie (https://ccb.jhu.edu/software/stringtie/ index.shtml?t=example) in a reference-based approach [19]. The sequencing depth ranged from 15 to 30 M; approximately 98% of reads achieve the sequencing quality score of Q20, and 95% of the reads have the sequencing quality score of Q30. The sequencing error rates of sequencing data were from 0.015 to 0.017 (Table S1).

Differential Expression Analysis and Functional Enrichment.
To identify differential expression genes between two groups, the expression level of each transcript was calculated according to the transcripts per million reads (TPM) method. RSEM (http://deweylab.biostat.wisc.edu/rsem/) was used to quantify gene abundances. The differential expression analysis was performed using the DESeq2, DEGseq, and EdgeR software. Volcano plots and scatter diagrams showed the signal intensity of differentially expressed genes. Hierarchical clustering was performed to display the expression patterns of differentially expressed genes in each sample. For the function analysis, we performed the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) to the profile of differently expressed genes. The GO analysis was composed of three parts, including biological process (BP), cellular component (CC), and molecular function (MF). All differently expressed genes were mapped to the GO and KEGG databases. P value <0.05 indicates that there is significant enrichment in the pathway.

Construction of Target Gene Prediction and Regulatory
Network. The correlation between differentially expressed mRNAs and differentially expressed lncRNAs was evaluated using Pearson's correlation coefficient from matched mRNA and ncRNA expression profile data. We set Pearson's correlation coefficients greater than 0.9 as significant. The potential target genes of lncRNAs were predicted via cisregulation or transregulation patterns. Cisregulation prediction suggests that the function of lncRNA is related to the proteincoding genes adjacent to the coordinates, and lncRNA located upstream and downstream of coding proteins may intersect with promoters or other cisacting elements of coexpressed genes, thus regulating gene expression at the transcriptional or posttranscriptional level. We analysed and screened the protein-coding genes nearest 10 kb upstream or downstream of lncRNA transcription start site. Transregulation prediction was performed using the RNAplex (University of Vienna, Vienna, Austria) software to find the mRNAs that had complementary sequences to the lncRNAs [20]. The interaction network was built and visually displayed using Cytoscape software.
2.6. Quantitative Real-Time PCR (RT-qPCR). 1 μg of the total RNA was converted to cDNA using the HiScript III RT SuperMix for qPCR kit (Catalog number R323, Vazyme, China) according to the manufacturer's instructions. The reaction was performed in 0.2 ml PCR tubes using thermocyclers. The expression of lncRNAs was detected by the Bio-Rad CFX96 system and analysed by the CFX Manager software. SYBR Green mix kit was used for progress tracking (Catalog number Q711, Vazyme, China). Briefly, reactions were composed of 1X ChamQ Universal SYBR qPCR mix, 200 nM forward primer, 200 nM reverse primer, 10 ng cDNA, and ddH 2 O (up to 20 μl). The thermal cycling conditions were as follows: 95°C for 30 s for activation, and then 40 cycles of 95°C for 10 s, 60°C for 30 s, 95°C for 15 s, 60°C for 60 s, and 95°C for 15 s. All experiments were performed in triplicate. GAPDH was set as a reference gene [21,22]. Relative expression of the gene was calculated using the 2 −ΔΔCt method. The primers were designed using Primer 3 software, and the efficiency of all primer pairs was better than 95% [23]. The sequences of primers are listed in Table S2. 2.8. Statistical Analysis. All statistics were analysed using GraphPad Prism software (version 6.02). The data were presented as mean ± SD, and the Student's t-test was used to determine differences between the two groups. P < 0:05 was considered a significant difference.  Table S3 and Figure S1. There is no significant differences in physical characteristics between the two groups. The level of blood pressure, heart rate, haemoglobin, CK, HDL, apolipoprotein, creatinine, uric acid, glucose, and MYO showed no difference between the two groups. However, HF patients presented higher LDH, CHO, TG, LDL, and NEFA contents. Moreover, HF patients also had higher CKMB, hsTNT, and NT-proBNP levels and lower LVEF levels. Among the HF subjects, the number of NYHA grade II patients accounted for 20%, grade III for 58%, and grade IV for 22%.

Disease Markers
3.2. The Profile of Differentially Expressed lncRNAs and mRNAs between HF and Control Groups. In order to recognize the dysregulated lncRNAs and their potential pathological role in HF, RNA sequencing was performed and the lncRNA and mRNA expression profiles were analysed. We found a total of 207 mRNAs were differentially expressed between the two groups (Figures 2(a) and 2(c)). Among them, 93 were upregulated and 114 were downregulated. In addition, a total of 422 lncRNAs showed significantly different expression levels between the two groups, of which 80 were upregulated and 342 were downregulated. Volcano and scatter plots are shown in Figures 2(b) and 2(d). As shown in Figures 3(a) and 3(b), the hierarchical clustering analysis revealed the regulatory profiles of mRNA and lncRNAs between the CT and HF groups. The top-ten dysregulated mRNAs and lncRNAs are, respectively, shown in Tables 2 and 3. The chromosome distribution of dysregulated lncRNAs is displayed in Figure 2(e). Chromosome 1 presented the greatest number of lncRNAs, followed by chromosome 17 and chromosome 19. According to the position of lncRNAs relative to protein-coding genes on the genome, lncRNAs can be divided into five types: exon-sense overlap, intergenic, antisense, intron sense overlap, and bidirectional. Antisense and intergenic lncRNAs accounted for a large proportion (Figure 2(f)).

Functional Enrichment Analysis.
We performed GO and KEGG analyses to determine the coexpression of mRNAs. The top 15 significantly enriched GO terms are shown in Figure 4. The aberrantly expressed mRNA at BP levels was associated with defence responses to other organisms, response to biotic stimulus, immune response, response to external biotic stimulus, response to external stimulus, immune system process, response to oxygen-containing compound, signal transduction, response to interferonalpha, and cell surface receptor signalling pathway (Figure 4(a)). In terms of CCs, the pathways were found to be related to the protein-DNA complex, DNA packaging complex, external side of plasma membrane, side of membrane, transcription factor AP-1 complex, plasma membrane part, plasma membrane, immunoglobulin complex, new growing cell tip, nuclear nucleosome, basolateral plasma membrane, interleukin-23 complex, and cyclin K-CDK13 complex (Figure 4(b)). In the category of MF, the mRNAs were enriched in the RNA polymerase II-specific, DNAbinding transcription activator activity, cysteine-type endopeptidase inhibitor activity, antigen binding, doublestranded DNA binding, lipopolysaccharide-binding, RNA polymerase II proximal promoter sequence-specific DNA binding, sterol transporter activity, proximal promoter sequence-specific DNA binding, transmembrane signalling receptor activity, RNA polymerase II activating transcription factor binding, HMG box domain binding, cargo receptor activity, cAMP response element binding, and chemokine (C-C motif) ligand 12 binding (Figure 4(c)).

Construction of Coexpression and Target Prediction
Network. To explore the interaction of aberrantly expressed genes in HF, a coexpression network was constructed based on correlation analysis. A total of 328 pairs of lncRNA-mRNA interactions (including 125 lncRNAs and 70 mRNAs) were selected, and the network was constructed with Cytoscape software ( Figure 5). For example, lncRNA ENSG0000 0224789 was an upregulated novel transcript which could interact with GBP1, BATF2, IFITM3, RTP4, IFIT3, and ISG15. HIST1H2AK was associated with ENSG00000260708, ENSG00000273338, ENSG00000231246, ENSG00000235919, ENSG00000271869, ENSG00000273117, and ENSG0000 0272426. Other lncRNAs (such as ENSG00000204283, ENSG00000275092, and ENSG00000272843) and mRNAs (DIABLO, HIST1H3A, ANKRD42, C17ORF107, FOSB, FOS, and HLA-DOA) also demonstrated interactions. These data suggest that these genes may play critical roles in the coexpression network. All the interactions of lncRNA and mRNA are listed in Table S4.
We performed cisregulation and transregulation to analyse the potential targets of differentially expressed lncRNAs. As shown in Figure 6, several lncRNAs were found to be associated with multiple target genes, including HISTH2BN, GTF2H2, HES2, DOCK4, and STX11. These associations may provide important references for the further study of lncRNAs. The detailed prediction information is listed in Table S5.    Figure 7, the expression of those lncRNAs was indeed significantly increased in HF patients. These results were consistent with the sequencing data. A ROC curve analysis was carried out to evaluate the diagnosis ability of HF ( Figure 8, Table S6). The AUC values of NT-proBNP, CHST11, MIR29B2CHG, AP000873.3, CR381653.1, FP236383.2, and DLEU2 were 0.844, 0.84, 0.732, 0.613, 0.744, 0.699, and 0.551, respectively. However, the P value of AP000873.3 and DLEU2 was 0.054 and 0.443, which were greater than 0.05.

Discussion
lncRNAs account for the majority of ncRNAs, which were previously thought to be the by-product of RNA polymerase II transcription and have no biological functions. However, a large number of lncRNAs have been confirmed to regulate cellular processes and be related to the occurrence and development of many diseases [24][25][26]. HF is considered the common outcome of various myocardial diseases, such as cardiomyopathy, myocardial infarction, and hypertension [27]. Despite advances in pharmacology and medical technology, a large proportion of patients with heart disease progress to advanced HF [28,29]. The prevalence of HF has reached 20% at the age of 40, and the incidence increases with age [30]. Previous studies have reported that many aberrantly expressed lncRNAs play a critical role in the pathogenesis and progression of cardiac hypertrophy and myocardial injury [31][32][33]. Yang et al. reported that the lncRNA MIAT could affect the handling of calcium and contractile function, resulting in cardiac myocyte remodelling and hypertrophy [11].
Blood is an easily accessible tissue and is associated with the occurrence of CVDs and the expression of their risk factors [34]. It can also provide information about a patient's status and can be expanded to very large sample sizes for screening biomarkers [21,35,36]. In this study, we collected 100 blood samples from 100 volunteers, including 50 HF patients and 50 healthy controls. A total of 207 mRNAs and 422 lncRNAs were differentially expressed in the HF group compared with the controls. Many new differential genes were found in our study, and these genes have not been comprehensively verified or studied; hence, our study provides a comprehensive differential gene expression map that may help understand part of the regulatory mechanism of HF.
In order to understand the function of differentially expressed lncRNAs in HF, we performed GO and KEGG analyses on its related mRNAs. And then, the functional enrichment analysis showed that the differentially expressed genes were mainly involved in immune response, DNA 9 Disease Markers packaging complex, DNA-binding transcription activator activity, RNA polymerase II-specific, and related to the biological processes of the plasma membrane and cell membrane. These pathways seemed to be associated with myocardial hypertrophy, cardiac aging, and HF. HF is the end-stage of heart disease, and its pathological factors involve cardiac hypertrophy, apoptosis, functional impairment, myocardial fibrosis, and cardiac remodelling [37][38][39][40]. Similar pathways were analysed by KEGG pathway analysis. The top 15 enrichment pathways included the NOD-like receptor signalling pathway, human T-cell leukaemia virus 1 infection, NF-kappa B signalling pathway, TNF signalling pathway, and IL-17 signalling pathway. These data demonstrated that the dysregulated genes were involved in immune and inflammatory responses. Our results were consistent with previous studies, which reported that inflammatory responses and apoptosis were associated with the pathogenesis of HF [41][42][43].
lncRNAs usually act as sponges of the other genes, including miRNAs and mRNAs, to exert their regulatory functions. Based on the results of RNA-sequencing, we performed correlation analysis of lncRNAs and mRNAs and constructed the coexpression network with bioinformatic tools. Our results found a connection between several mRNAs and the dysregulated lncRNAs. To cite one, we found DIABLO possessing a high degree of connection. Additionally, the expression of DIABLO was associated with cardiomyocyte apoptosis, which is induced by cardiac ischemia [44]. Posttranslational protein modification plays an emerging and vital role in the occurrence and development of cardiovascular diseases [45]. HIST1H3A has been reported to regulate protein acetylation in cardiovascular pathologies [46]. Moreover, ischemia-reperfusion injury could affect the localization of FOSB in the cardiac myocyte [47]. In the cisregulation and transregulation networks, we found several mRNAs that could be potential targets of  Disease Markers lncRNAs. The high degree of centrality of the mRNA included FOSB, DOX4, and CDK19. Among those genes, CDK19 was a component of the transcriptional mediator kinase module regulating gene transcription [48]. Recent studies also found dysregulated CDK19 to be associated with coronary artery diseases [49]. EGR1 is a key gene in the process of cell differentiation and mitogenesis. Upregulated EGR1 can trigger inflammatory responses and cell death [50]. The high expression of EGR1 was found in the HF group, and this result was consistent with a previous study [35]. Although we found that several lncRNAs were related to specific mRNAs, further investigations are still needed to verify their relationship. We selected the top six upregulated lncRNAs (CHST11, MIR29B2CHG, AP000873.3, CR381653.1, FP236383.2, and DLEU2) and detected their expression in all samples. The val-idation results were consistent with those of the sequencing data. lncRNA CHST11 is a transcript variant of carbohydrate sulfotransferase 11, and the results of sequencing and RT-qPCR were both upregulated in the HF group. Previous studies reported that MIR29B2CHG dysregulation is associated with colorectal cancer, pulmonary adenocarcinoma, breast cancer, and adrenocortical carcinoma [51][52][53][54]. Interestingly, in this study, we found the expression levels of MIR29B2CHG to be upregulated in the HF group. DLEU2 has been reported to play a crucial role in the progression of cancers such as the case of liver cancer and endometrial cancer [55,56]. Some researchers have also reported that DLEU2 could regulate the differentiation of the skeletal muscle [57]. The lncRNAs CR381653.1, AP000873.3, and FP236383.2 are novel upregulated transcripts being reported for the first time in the present study. lncRNA AP000873.3 is an exon-sense overlap We found the expression of these lncRNAs to be highly different between the HF and control groups. According to previous studies, the highly dysregulated lncRNAs in cardiovascular diseases could be considered biomarker of these diseases [42,58]. We hypothesized that the expression of these lncRNAs was related to the pathogenesis of HF. So, we further performed ROC curve analyses to assess their diagnostic value. Moreover, we also assessed the diagnostic value of NT-proBNP. The results showed that the AUC value of lncRNAs CHST11, CR381653.1, MIR29B2CHG, and FP236383.2 were all above than 0.6, with the exception of CHST11 whose value was close to that of NT-proBNP. Conversely, the AUC values of AP000873.3 and DLEU2 were 0.613 and 0.551, respectively, and their P values were greater than 0.5. Although the expression differences between the two lncRNAs were significant, there was no significant potential relationship between their expression levels and HF in terms of diagnostic value, which we believe is due to the insufficient sample size. NT-proBNP levels are known to be affected by many factors, such as age, drug intervention, obesity, tachycardia, renal insufficiency, sepsis, chronic obstructive pulmonary disease, and cirrhosis [59,60]. Altogether, these results suggest that these lncRNAs may have a potential value as a diagnostic indicator for HF.

Disease Markers
In an earlier study, Liu et al. found that RNA sequencing technology had better resolution and lower error rates compared to other analysis methods while requiring a smaller sample size [61]. Zheng et al. reported lncRNA GAS5 to be downregulated in HF based on bioinformatic analysis. In addition, GAS5 may functions as a competing endogenous RNA through modulating miRNAs and mRNAs in the progression of HF [62]. Chen et al. revealed the coexpression network of mRNAs and lncRNAs in pressure overloadinduced HF by sequencing rat models [42]. At present, we selected 8 samples (4 HF and 4 controls) for sequencing analysis, while the differently expressed lncRNAs were detected and identified in 100 samples. We found that selected lncRNAs have high discrimination between HF and controls. These results suggest that these novel lncRNAs may function as potential biomarkers for HF.
However, this study also had several limitations. Even though the sample size was large compared to many previous studies, these lncRNAs still need to be further analysed in even larger-scale studies. All the participants in this study were Chinese, and thus, these findings should be confirmed in other populations. During the sample inclusion process, we included HF patients based on the combined clinical and lab test diagnosis, without being able to exclude the heterogeneity in their etiology, NYHA (New York Heart Association Class) class and/or HF type. In this current study, NYHA grade III and HFmrEF (HF with mildly reduced EF) types account for the majority of the HF participants. In addition, the underlying molecular mechanisms of how these dysregulated lncRNAs regulate HF still need to be explored.

Conclusions
In conclusion, our work revealed the expression profiles of differentially expressed lncRNAs and mRNAs between the HF and control groups by transcriptome sequencing. Several differentially expressed mRNAs are reported to be involve in various biological pathways related to the pathogenesis of HF. The lncRNA-mRNA coexpression network was conducted to find the hub genes related to the pathology of HF. Moreover, combined sequencing and validation analyses report several lncRNAs with high expression levels that present a confident diagnostic value for HF. These findings provide new insights for understanding the underlying mechanisms and therapeutic targets of HF.

Data Availability
The data used to support the findings of this study have been deposited in the NCBI's Sequence Read Archive (SRA), reference number (BioProject ID: PRJNA859795).

Ethical Approval
The study was conducted in accordance with the Declaration of Helsinki and approved by the ethics committee of the Affiliated Hospital of Qingdao University.

Conflicts of Interest
The authors declare that no conflict of interest regarding the publication of this paper.

Authors' Contributions
LHHA, PL, and HY designed the project. XC and LHHA analysed the data. XC wrote the manuscript, and LHHA edited the manuscript. CS, JG, YW, XL, and CZ supervised the work and reviewed the manuscript. JCCJ revised and polished the manuscript. All authors read and approved the final manuscript. 14 Disease Markers