A study of the differential expression profiles of Keshan disease lncRNA/mRNA genes based on RNA-seq

Introduction

Keshan disease (KSD) is an endemic cardiomyopathy unique to China, and patients are mainly found in the narrow and low-selenium zone from the northeast to the southwest of China. According to statistics (1), 16 provinces (regions), with a population of 60 million, are affected. KSD is mainly characterized by myocardial involvement, and the clinical manifestations are ventricular enlargement and a decrease in myocardial contractility, which often causes heart failure, arrhythmia and even sudden death. Thus, the prognosis is poor.

The cause of KSD is not yet fully understood, but research results confirm there exists a close association between low selenium levels and KSD (2,3). As an essential trace element, selenium realizes its biological function mainly through selenoproteins. At present, a large number of studies have revealed that selenium deficiency is closely related to the occurrence and development of KSD. The selenium content of the environment of the KSD area, as found in the soil, local food, and hair and blood of the people, and the plasma selenoprotein P are lower than those in non-affected areas (4,5). A systematic review and meta-analysis of selenium deficiency and KSD revealed that the supplementation of selenium could significantly reduce the incidence of KSD (2).

Current studies suggest that the mechanism of chronic KSD is mainly related to mitochondrial oxidative damage, apoptosis and fibrosis regulation (6-8). Low selenium can cause damage to myocardial mitochondria, and selenium deficiency can lead to a decrease in the synthesis of glutathione peroxidase. This then leads to the aggravation of oxidative stress damage, in turn causing myocardial damage and, ultimately, resulting in cardiomyopathy. Selenium deficiency may also amplify other pathogenic factors, including viral and other infections (3). Recent comparative studies on the gene expression of KSD and dilated cardiomyopathy have found that KSD and dilated cardiomyopathy have a different pathogenesis: the pathogenesis of KSD mainly manifests as viral infection, cell membrane damage and oxidative stress (9).

In recent years, the advances in RNA sequencing technology mean that it has developed into a powerful research platform indispensable for exploring diseases. At present, it is widely used for research into a variety of complex human diseases, such as malignant tumors, immune diseases and unknown infectious diseases, and it is crucial to the search for pathogens and biomarkers, the exploration of pathogenesis, and the formulation of targeted treatment programs (10). In this study, RNA-sequence (RNA-seq, a transcriptome sequencing technology) was used to construct differential expression profiles of lncRNA in the peripheral blood plasma of KSD patients. These profiles were then analyzed by means of bioinformatics, providing a theoretical basis for shedding light on the mechanism of disease at the level of gene expression regulation.

We present the following article in accordance with the MDAR reporting checklist (available at: http://dx.doi.org/10.21037/cdt-20-746).

Methods

Clinical data

Ten KSD patients in the severe KSD area of Shandong Province were enrolled as the KSD group. There were five male patients and five female patients, ranging from 38 to 45 years of age, with an average age of 42.6 years. All the patients had their medical histories taken, and underwent a physical examination, electrocardiography, X-ray photography, cardiac ultrasound and related biochemical examinations, and the examination results met the diagnostic criteria of KSD (GB17021-1997). In addition, 10 healthy people, five males and five females, were enrolled locally as the control group, having undergone a comprehensive physical examination to confirm that they had no heart disease. The clinical data of the KSD group and control group are presented in Table 1. The present study was conducted in accordance with the Declaration of Helsinki (as revised in 2013) and was approved by the ethics committee of Liaocheng People’s Hospital (No. 201510), and all patients signed written informed consent.

Table 1 Clinical data of the KSD group and control group
Full table

Main reagents and instruments

Main reagents: Trizol kit (Invitrogen life technologies, USA), DEPC-treated water (Invitrogen life technologies, USA), RNase inhibitor (Epicentre, USA), SuperScript^TM III reverse transcriptase (Invitrogen, USA), 5× RT buffer (Invitrogen, USA), 2.5 mM dNTP mixture (HyTest Ltd., Finland). Main instruments: NanoDrop ND-2000 instrument (Thermo Fisher Scientific, USA), Ribo-Zero rRNA Removal kits (Illumina, USA), TruSeq Stranded Total RNA Library Prep Kit (Illumina, USA), BioAnalyzer 2100 instrument (Agilent Technologies, USA), Illumina HiSeq4000 (Illumina, USA).

Experimental methods

Specimen acquisition

The subjects fasted for 12 hours, and blood was withdrawn from the elbow vein in the morning, after they had got up. Five mL of peripheral venous blood was withdrawn from each subject in both the KSD group and the control group, with an anticoagulation negative pressure tube, and stored at −80 °C.

RNA extraction

The blood sample was removed from the −80 °C refrigerator, thawed, and centrifuged at 12,000 g at 4 °C for 10 minutes to remove impurities. Next, 250 µL of the blood was taken and transferred to a 1.5-mL centrifuge tube, 750 µL of TRIzol LS reagent and 20 µL of glacial acetic acid were added, and the tube was shaken to mix the solution up well. After homogenization, the sample was incubated at 15 to 30 °C for five minutes, 0.2 mL of chloroform was added, the tube was vigorously shaken for 15 seconds, incubated at 30 °C for two minutes, and then centrifuged at 12,000 g and 4 °C for 15 minutes. After centrifugation, the blended liquid had separated out into layers, with the RNA in the colorless water layer. The upper water layer was transferred to another centrifuge tube, and 0.5 mL of isopropanol was added to it. The mixture was left to stand at room temperature for 10 minutes, then centrifuged at 12,000 g and 4 °C. After 10 minutes, RNA gel-like precipitates could be seen on the bottom and sides of the tube. The upper liquid layer was removed, and 1 mL of 75% ethanol was added to clean off the gelatinous sediment attached to the tube sides. It was fully oscillated, centrifuged at 75,000 g and 4 °C for five minutes. The upper liquid layer was removed again, and the sediments were left standing but not to completely dry out. About 5–10 minutes later, RNase-free water was added to dissolve the RNA, and the solution was incubated at 55–60 °C for 10 minutes and, finally, stored at −80 °C.

CDNA library preparation and sequencing

A Ribo-Zero rRNA Removal Kit (Illumina, USA) was used to remove rRNAs from the RNA, and a TruSeq Stranded Total RNA Library Prep Kit (Illumina, USA) was used to pre-process the RNA to construct a sequencing library. A BioAnalyzer 2100 instrument (Agilent Technologies, USA) was used for library quality control and quantification. In accordance with the Illumina sequencing instructions, the 10 pM library was denatured into single-stranded DNA molecules, which were captured on an Illumina flow cell and amplified in situ into clusters. They then underwent 150 cycles of sequencing on the Illumina HiSeq sequencer.

Statistical analysis

Original data was obtained after sequencing by an Illumina HiSeq4000 sequencer. The first quality control of the raw data was performed using the Q30 value. Cutadapt software (v1.9.3) was used to remove adapters and low-quality reads, and then Hisat2 software (v2.0.4) was used under the guidance of an Ensembl gtf gene annotation file to align high-quality adapter-free reads to the human reference genome (UCSC MM10). Cuffdiff software (v2.2.1) was used to obtain a gene-level lncRNA/mRNA FPKM (Fragments Per Kilobase of transcript per Million fragments mapped) value, as the lncRNA/mRNA expression profile. Agilent Gene Spring software was used for data analysis, and standardized processing was performed to obtain two sets of fold change values. After statistical analysis using a t-test, differential lncRNA/mRNA expression profiles were obtained. The standard was that the fold change in the two groups was more than double and P<0.05. Gene Spring Software 11.0 was used to do a cluster analysis of the differentially expressed lncRNA/mRNA. As a regulatory factor, lncRNA can regulate adjacent protein-coding genes, and target gene detection finds differentially expressed genes and marks them as target genes. Differential genes were subjected to Gene Ontology (GO) biological function analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. P<0.05 was considered statistically significant.

Results

Clinical data of the KSD group and the control group

The differences in gender composition, age, systolic blood pressure (SBP), diastolic blood pressure (DBP) and heart rate between the two groups were not statistically significant; compared with the control group, in the KSD group, the left ventricular end-diastolic diameter and left atrial diameter were significantly larger, the left ventricular ejection fraction was significantly lower (P<0.01), and the heart function was significantly lower (P<0.01, Table 1).

Quality control of RNA and library

The RNA concentration of each sample was measured by NanoDrop ND-2000 (Thermo Fisher Scientific, Waltham, MA, USA), and the OD260/OD280 value was used as an indicator of RNA purity. In this study, the OD260/OD280 range was 1.8–2.1, which is to say, the purity of the RNA was qualified. RNA integrity was tested using denatured agarose gel electrophoresis. It was observed that the RNA5s, 18S and 28S bands were clear and intact without any obvious degradation.

Differential expression of lncRNA/mRNA

By using the RNA-seq technique, a total of 89,905 lncRNAs and 20,315 mRNAs were detected (Figure 1). Data analysis revealed that the difference in the expression of lncRNA in plasma between the KSD group and the control group was statistically significant. Statistical analysis revealed that 921 lncRNAs had obvious differential expression (fold change >2, P<0.05), of which 36 were up-regulated and 885 were down-regulated; 2,771 mRNAs presented with obvious differential expression, of which 253 were up-regulated and 2,518 were down-regulated. A cluster analysis of the differential expression of lncRNA was conducted and the difference in lncRNA expression between the KSD group and the control group could clearly be seen. The same group of samples was gathered together, and the results indicated that the gene expression trend was consistent (Figure 2). The differentially expressed lncRNAs were tested for target genes, and 117 genes were found to be regulated by differential lncRNA.

Figure 1 Differential expression profiles of lncRNA/mRNA. KSD, Keshan disease; CTRL, control.

Figure 2 Cluster analysis of the differential expression of lncRNA. (A) Cluster analysis of differentially expressed lncRNA between the KSD group (KSD 1-10) and control group (ctrl1-10), the red part represents the up-regulation of lncRNA, and the green part represents the down-regulation of lncRNA. (B) Cluster analysis of mRNA differentially expressed between the KSD group (KSD 1-10) and control group (ctrl1-10), the red part represents the up-regulation of mRNA, and the green part represents the down-regulation of mRNA. KSD, Keshan disease.

Bioinformatic analysis

In order to determine the biological process (BP) and molecular mechanism of differentially expressed genes in KSD, an enrichment analysis of the function and pathway of differentially expressed genes was also carried out. An enrichment analysis of the function and pathway was carried out using the database for annotation, the visualization and integrated discovery (DAVID) online tool (11). The tool is based on GO and the KEGG database. GO analysis is a biological gene bank used to annotate genes, gene products and sequences. It is divided into three parts, and determines the molecular function (MF), BP and cell component (CC). In this study, a CO function analysis of the differential expression of lncRNA was conducted to learn more about the function of these lncRNAs, GO terms with a P value of <0.05 were considered to be statistically significant. GO function analysis revealed that the differentially expressed lncRNAs were mainly involved in the BPs of neuronal differentiation, cell morphogenesis and differentiation, actin filament binding, superoxide anion generation and hydrogen peroxide biosynthesis (Figure 3). A KEGG pathway analysis was done to map molecular data sets from genomics, transcriptomics, proteomics and metabonomics to a KEGG pathway map for biological function interpretation, with P≤0.05 as the threshold of significant enrichment. The KEGG pathway analysis showed that differentially expressed lncRNAs were enriched in six signaling pathways (Table 2), among which the FOXO signaling pathway ranked first (Figure 4).

Figure 3 Functional distribution of differentially expressed lncRNA by GO analysis. (A) GO enrichment analysis of differentially expressed lncRNAs in biological process (BP). (B) GO enrichment analysis of differentially expressed lncRNAs in cell components (CC). (C) GO enrichment analysis of differentially expressed lncRNAs in molecular functions (MF). GO, Gene Ontology.

Table 2 Signaling pathways with enriched differential genes
Full table

Figure 4 FoxO signaling pathway (hsa04068).

QRT-PCR validation

The names and the fold changes of the top 20 lncRNA/mRNA are in Table 3. Three significantly up-regulated lncRNAs and three significantly down-regulated lncRNAs were randomly selected (ENST00000421185, ENST00000410140, TCONS_l2_00004571; down-regulated: ENST00000355500, ENST00000359888, ENST00000376025). A quantitative real-time polymerase chain reaction (qRT-PCR) was performed on 10 patients with KSD and 10 healthy controls. It was revealed that, in the expanded samples, the PCR results were basically consistent with the RNA-seq results (Figure 5). lncRNA/mRNA network is in Figure 6.

Table 3 The names and the fold changes of the top 20 lncRNA/mRNA
Full table

Figure 5 lncRNA (up-regulation: ENST00000421185, ENST00000410140, TCONS_l2_00004571; down-regulation: ENST00000355500, ENST00000359888, ENST00000376025) the results in the expanded samples verified by RT-PCR (*, P<0.05; **, P<0.01; ***, P<0.001). KSD, Keshan disease; CTRL, control.

Figure 6 ncRNA/mRNA network for lncRNA-mRNA with Pearson correlation coefficient >0.99.

Discussion

KSD was first found in Keshan County, Heilongjiang Province, China, in 1935. After the Second World War, similar cases were also found in the Nagano Prefecture of Japan and the northern mountainous area of Korea (12,13). After more than half a century of prevention and control, the type of KSD also changed from acute KSD and subacute KSD to chronic KSD. In patients with chronic KSD, the heart is significantly enlarged and the cardiac function is decompensated. This is mainly characterized by severe cardiomyopathy and multifocal myocardial necrosis, which are the focus of prevention and treatment research at present. Therefore, in this study, patients with chronic KSD were chosen as the study subjects, and RNA-seq technology was used to construct the differential expression profiles of lncRNA/mRNA in KSD, in order to explore the molecular mechanism and target genes of KSD at the gene level.

At present, research into KSD mainly focuses on the association between selenium deficiency and KSD. However, studies on the mechanism have mainly focused on selenium-related mitochondrial oxidative damage, apoptosis and myocardial fibrosis. Selenium deficiency may also increase the infection of enterovirus and other pathogens, through the above-mentioned mechanism, leading to the occurrence and progress of KSD, and finally causing heart failure (3,6-8). In a previous study by this research team, miRNA sequencing was also conducted in patients with dilated cardiomyopathy, and bioinformatic analysis revealed that a neuronal differentiation process and a mitogen-activated protein kinase (MAPK) signaling pathway were involved in the occurrence and development of dilated cardiomyopathy (14). Therefore, there are some differences in miRNA between KSD and dilated cardiomyopathy.

RNA-seq can be used to study gene expression differences at the whole genome level, which has the advantages of more accurate quantification, higher repeatability, a wider detection range and more reliable analysis. It uses high-throughput sequencing technology to sequence all cDNA library from RNA reverse transcription in tissues or cells in order to calculate the expression levels of different RNAs by counting the number of related reads. Thus, more comprehensive genetic information, such as transcriptional location and splicing, can be determined (10).

In the present study, a total of 117 genes were found to be regulated by differential lncRNAs. GO function analysis of the differentially expressed genes revealed that the differentially expressed lncRNAs were mainly involved in the BPs of neuronal differentiation, cell morphogenesis and differentiation, actin filament binding, superoxide anion generation and hydrogen peroxide biosynthesis. The results showed that KSD patients differ from normal people in cardiomyocyte differentiation and oxidative stress. Pei et al. (7) found that the level of oxidative stress increased significantly and was positively correlated with the degree of myocardial damage in patients with KD. The present study confirmed the argument that oxidative stress is related to myocardial injury in KSD at the gene expression level.

KEGG pathway analysis demonstrated that the differentially expressed lncRNAs were enriched in six signaling pathways, in which the FoxO signaling pathway, which is related to apoptosis, ranked first. FOXO plays an important role in cell metabolism, cell survival and resistance to oxidative stress (15). It has been confirmed in one study that by increasing the resistance of mice to oxidative stress, FOXO could reduce myocardial injury. FOXOs play a critical role in promoting cardiomyocyte survival during conditions of oxidative stress through the induction of antioxidants and cell survival pathways (16). Another experiment verified that FOXO expression was down-regulated in mice with oxidative stress-induced myocardial injury (17). This study revealed that the differentially expressed lncRNAs were significantly enriched in the FoxO signaling pathway and noticeably down-regulated. It can be inferred that the occurrence and development of KSD may be related to the weakening of myocardial resistance to oxidative stress caused by the down-regulation of the FoxO signaling pathway. The intervention of this signaling pathway may provide new ideas for the treatment of KSD. The key genes affecting the occurrence and development of KSD were further screened using the RNA-seq technique, and the two candidate genes detected in the preliminary results were IGF1R and TGFB2. IGF1R participates in the process of myocardial hypertrophy by changing cell cycle, apoptosis and interaction with the surrounding environment. A previous study revealed that the IGF1 level in the myocardium of patients with idiopathic hypertrophic cardiomyopathy was significantly higher than that in patients with aortic stenosis. Immunohistochemistry and in situ hybridization revealed that IGF1 expression was up-regulated in patients with idiopathic hypertrophic cardiomyopathy. Therefore, it is suggested that the up-regulation of IGF1 expression is not the result of myocardial hypertrophy but the cause of myocardial hypertrophy (18). One study induced IGF1R silencing with RNA interference, and the results were that the severity of myocardial hypertrophy was significantly reduced, and the myocardial function was improved. TGFB2, another gene screened candidate, is a transforming growth factor. A further study revealed that, after RXRA gene knockout, TGFB2 gene expression and cardiomyocyte apoptosis increased, and the embryonic cardiac outflow tract development was abnormal (19). The investigators believed that the abnormal expression of a target gene is related to the occurrence and development of KSD. The next step is to construct a mouse model, with oxidative stress injury, selenium deficiency and vitamin deficiency, to verify this conclusion, using RNA interference induced gene silencing.

In summary, in this study, RNA-seq technology was used to construct the differential expression profiles of lncRNA/mRNA in KSD, a bioinformatics analysis of differentially expressed genes was performed, and IGF1R and TGFB2 were preliminarily screened out as candidate genes. The investigators will next conduct animal experiments to verify these candidate genes, and continue to study the differential gene expression and regulation mechanism of KSD by increasing the sample size, in order to reveal the mechanism of the occurrence and development of cardiomyopathy of KSD and, ultimately, facilitate drug treatment.

Acknowledgments

We thanked all people who make contributions to this work.

Funding: National Natural Science Foundation of China (81573095).

Footnote

Reporting Checklist: The authors have completed the MDAR reporting checklist. Available at: http://dx.doi.org/10.21037/cdt-20-746

Data Sharing Statement: Available at: http://dx.doi.org/10.21037/cdt-20-746

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at: http://dx.doi.org/10.21037/cdt-20-746). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The present study was conducted in accordance with the Declaration of Helsinki (as revised in 2013) and was approved by the ethics committee of Liaocheng People’s Hospital (No. 201510), and all patients signed written informed consent.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Yan C, Fang W. Status and etiology progress of Keshan disease. Advances in Cardiovascular Diseases 2017;38:225-9.
2. Zhou H, Wang T, Li Q, et al. Prevention of Keshan Disease by Selenium Supplementation: a Systematic Review and Meta-analysis. Biol Trace Elem Res 2018;186:98-105. [Crossref] [PubMed]
3. Loscalzo J. Keshan disease, selenium deficiency, and the selenoproteome. N Engl J Med 2014;370:1756-60. [Crossref] [PubMed]
4. Zhang X, Wang T, Li S, et al. A Spatial Ecology Study of Keshan Disease and Hair Selenium. Biol Trace Elem Res 2019;189:370-8. [Crossref] [PubMed]
5. Zhang X, Wang T, Li S, et al. A spatial ecological study of selenoprotein P and Keshan disease. J Trace Elem Med Biol 2019;51:150-8. [Crossref] [PubMed]
6. Jiang S, Hou J, Liu H, et al. The roles of PPARs and PPARγ coactivator-1 alpha in mitochondrial cardiomyopathy. Chin J Endemiol 2018;37:943-6.
7. Pei J, Fu W, Yang L, et al. Oxidative stress is involved in the pathogenesis of Keshan disease (an endemic dilated cardiomyopathy) in China. Oxid Med Cell Longev 2013;2013:474203 [Crossref] [PubMed]
8. Ma ZG, Yuan YP, Wu HM, et al. Cardiac fibrosis: new insights into the pathogenesis. Int J Biol Sci 2018;14:1645-57. [Crossref] [PubMed]
9. Huang GY, Xiang YZ, Liu JW, et al. Comparative analysis of lncRNA-mRNA co-expression between Keshan disease and dilated cardiomyopathy. Chin J Endemiol 2019;38:361-7.
10. Hawkins RD, Hon GC, Ren B. Next-generation genomics: an integrative approach. Nat Rev Genet 2010;11:476-86. [Crossref] [PubMed]
11. Huang da W. Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc 2009;4:44-57. [Crossref] [PubMed]
12. Tan J, Huang Y. Selenium in geo-ecosystem and its relation to endemic diseases in China. Water Air Soil Pollut 1991;57-58:59-68. [Crossref]
13. Kakehashi Y. Why did Keshan disease occur? The relationship between White muscle disease, Shinshu cardiomyopathy and Keshan disease. Chin J Endemiol 2000;19:150-2.
14. Huang G, Cao M, Huang Z, et al. Small RNA-sequencing identified the potential roles of neuron differentiation and MAPK signaling pathway in dilated cardiomyopathy. Biomed Pharmacother 2019;114:108826 [Crossref] [PubMed]
15. Nho RS, Hergert P. FoxO3a and disease progression. World J Biol Chem 2014;5:346-54. [Crossref] [PubMed]
16. Sengupta A, Molkentin JD, Paik JH, et al. FoxO transcription factors promote cardiomyocyte survival upon induction of oxidative stress. J Biol Chem 2011;286:7468-78. [Crossref] [PubMed]
17. Guo Y, Li Z, Shi C, et al. Trichostatin A attenuates oxidative stress-mediated myocardial injury through the FoxO3a signaling pathway. Int J Mol Med 2017;40:999-1008. [Crossref] [PubMed]
18. Li RK, Li G, Mickle DA, et al. Overexpression of transforming growth factor-beta1 and insulin-like growth factor-I in patients with idiopathic hypertrophic cardiomyopathy. Circulation 1997;96:874-81. [Crossref] [PubMed]
19. Wang L, Li X, Huang W, et al. TGFβ signaling regulates the choice between pluripotent and neural fates during reprogramming of human urine derived cells. Sci Rep 2016;6:22484. [Crossref] [PubMed]