Characterization of novel lncRNAs in upper thoracic spinal cords of rats with myocardial ischemia-reperfusion injuries.

Myocardial ischemia-reperfusion injury (MIRI) is a significant problem in clinical cardiology, and refers to a more serious myocardial injury caused by blood recanalization after a period of myocardial ischemia, as compared with injury caused by vascular occlusion. The spinal cord, as the primary afferent and efferent center of cardiac sensory and sympathetic nerve fibres, has received increased attention in recent years with regards to the regulation of MIRIs. Previous studies have revealed that MIRI has a strong correlation with the abnormal expression of long non-coding (lnc)RNAs in the myocardium; however, there are limited reports on the effects of the altered expression of lncRNAs in the spinal cord following MIRI. To investigate the expression patterns of lncRNAs in the spinal cord after MIRI and their potential role in the early stage of reperfusion, a MIRI model was established in rats. After 30 min of myocardial ischemia and 2 h of reperfusion, the upper thoracic spinal cord tissues were immediately dissected and isolated. lncRNAs and mRNAs in spinal cord tissues were screened using transcriptome sequencing technology, and the expression of several highly deregulated mRNAs, including Frs3, Zfp52, Dnajc6, Nedd4l, Tep1, Myef2, Tgfbr1, Fgf12, Mef2c, Tfdp1 and lncRNA, including ENSRNOT00000080713, ENSRNOT00000090564, ENSRNOT00000082588, ENSRNOT00000091080, ENSRNOT00000091570, ENSRNOT00000087777, ENSRNOT00000082061, ENSRNOT00000091108, ENSRNOT00000087028, ENSRNOT00000086475, were further validated via reverse transcription-quantitative PCR. The number of altered expressed lncRNAs was 126, among which there were 41 upregulated probe sets and 85 downregulated sets. A total of 470 mRNAs were differentially expressed, in which 231 probe sets were upregulated and 239 were downregulated. Gene Ontology analysis indicated that dysregulated transcripts related to biological processes were mainly associated with 'cell-cell signaling'. Moreover, pathway analysis demonstrated significant changes in the 'PI3K/Akt signaling pathway' and the 'p53 signaling pathway'. Thus, the altered expression of lncRNAs in the spinal cord may be of considerable importance in the process of MIRI. The present results could provide an insight into the potential roles and mechanism of lncRNAs during the early stage of reperfusion. Copyright: © Li et al.

Introduction

Ischemic heart disease (IHD) is a leading cause of premature mortality worldwide (1,2). IHD affects ~126 million individuals globally (1,655 per 100,000), which is estimated to be 1.72% of the world's population (1). Restoration of coronary artery perfusion is currently one of the main treatment options but it can lead to myocardial ischemia-reperfusion injury (MIRI), which refers to a more serious myocardial injury caused by blood recanalization after a period of myocardial ischemia compared with that induced by vascular occlusion (3,4). MIRI is reported to increase the risk of cardiovascular disease, which imposes health, social and human burdens on societies worldwide, and thus is a public health issue (5-8). Therefore, treatment of MIRI is important for the prevention of cardiovascular disease, the end stage of cardiac dysfunction and mortality.

Previous studies have reported that the spinal cord is closely associated with the occurrence and development of MIRI (9-11). Following cardiac ischemia and reperfusion, chemical mediators such as dynorphin and suppresses substance P within the myocardium activate cardiac afferent fibres projecting from the dorsal root ganglion to the upper thoracic dorsal horn of the spinal cord. This results in the release of inflammatory factors and neuropeptides within the spinal cord (12-15), which triggers abnormal changes in gene transcription and translation within the spinal cord (16-18). Despite these recognized spinal responses, the genomic mechanisms that enable differences during these spinal processes of MIRI remain largely elusive.

It has been revealed that long non-coding (lnc)RNAs serve a key role in the processing of gene expression signals via the regulation of cell cycle distribution, stem cell differentiation and transcriptional and post-transcriptional control (19,20). Moreover, there are aberrant expression levels of spinal lncRNAs in peripheral diseases, such as pruritus (21,22), acute kidney injury (23), neuropathic pain and inflammation-related pain (24).

In the present study, the altered expression of mRNAs and lncRNAs in the thoracic (T)1-T4 spinal cord were analyzed using high-throughput RNA-sequencing (RNA-seq) between a group of rats established with MIRI and a sham group. Using microarrays, the current study identified 470 mRNAs and 126 lncRNAs that were significantly differentially expressed in rats with MIRI, and further analyzed the biologic roles of these mRNAs and lncRNAs by performing Gene Ontology (GO) and pathway analyses. Subsequently, the differentially expressed mRNAs and lncRNAs were validated using reverse transcription-quantitative (RT-q) PCR. The present results may facilitate the development of molecular-targeting therapies for MIRI.

Materials and methods

Animals

In total, 26 male Sprague Dawley rats (specific pathogen free grade; weight, 250-300 g; age, 8-10 weeks) from Hunan Slack Jingda Experimental Animal Company were used in the present study. The animals were housed at a controlled room temperature (25±0.5˚C) and relative humidity (40-60%) with a 12-h light/dark cycle and had unlimited access to food and water. The number of animals in each group was three, and each animal was used only once. All experiments were performed with the approval of the Institutional Animal Care and Use Committee of Tongji Hospital, Huazhong University of Science and Technology (IRB ID, TJ-A0804) and the National Institutes of Health Guide for the Care and Use of Laboratory Animals (25).

MIRI

Rats were randomly divided into MIRI (model, n=13) and sham (control, n=13) groups. MIRI experiments were performed according to previously published protocols (18,26-29). The rats were maintained in deep anesthetic state on a heating pad using sodium pentobarbital (50 mg/kg; intraperitoneal), which was indicated by disappearance of the foot withdrawal reflex.

In the model group, the left anterior descending coronary artery (LAD) was ligated 2-mm below the left atrial appendage for 30 min and then reperfused for 2 h, while sutures were placed without LAD ligation in control operations. The reperfused rats were considered as the animal model of MIRI. Sham-operated age-matched rats served as controls. During the experiment, the electrocardiogram (ECG) of animals was monitored using a electrocardiogram machine (model no. ECG-2303B; Guangzhou Sanrui Electronic Technology Co., Ltd.) at the following time points: Before the operation, ischemia for 0, 15 and 30 min, reperfusion for 0 min, 30 min and 2 h. After reperfusion for 2 h, 150 mg/kg sodium pentobarbital via intraperitoneal injection was used for euthanasia and the confirmation of euthanasia was decapitation. Then, the hearts were collected for 2,3,5-triphenyltetrazolium chloride (TTC) staining, hematoxylin and eosin (H&E) staining and Masson staining. At the end of reperfusion, the rats were deeply anesthetized, blood (~1 ml) samples were collected from the aorta for cardiac troponin I (cTnI) detection before euthanasia.

In the present experiment, ~20% of the animals died due to hemorrhage and arrhythmia caused by ligation of coronary artery at the early stage of reperfusion. When ventricular fibrillation and bleeding were uncontrollable, euthanasia was performed on animals as humane endpoints. Time points of 30 min for ischemia and 2 h for reperfusion were selected for this animal model is sufficient to result in significant histopathological and functional myocardial injury, which can be confirmed using the results of ECG, cTnI, TTC, H&E and Masson staining (30).

cTnI detection and myocardial tissue staining

cTnI level was determined using an Automated Blood or Urine Chemical Analyzer (Vitro 350; Ortho Clinical Diagnostics). H&E staining and Masson trichrome staining were conducted in the myocardial tissue samples in order to observe the myocardial pathology. Heart tissues were fixed in 10% formaldehyde for 24 h at room temperature and embedded with paraffin, before 4-µm thick sections were cut. After deparaffinization by washing with xylene followed by rehydration in a descending ethanol gradient, sections were immersed in hematoxylin for 5-7 min, differentiated in 1% acid alcohol for 2-5 sec and stained with 0.5% eosin for 2 min. After rinsing with distilled water for 30 sec, sections were dehydrated with an ascending ethanol gradient and were cleared in xylene. All procedures were carried out at room temperature. Using a light microscope (magnification x200; Leica DM 4000B; Leica Microsystems, Inc.), pathological changes of myocardial tissues were observed. The procedure of Masson's trichrome staining was similar to H&E staining, which was described previously (31).

The infarct size was assessed using TTC staining 2 h after reperfusion. After surgery, the hearts were removed and frozen for 20 min at -20˚C, then transversally cut into sections with a thickness of 1-2 mm. Tissue sections were incubated in 2% TTC for 10 min at 37˚C in dark conditions and then fixed in 10% formaldehyde at 4˚C overnight. Under light microscopic observation, the infarct area was a white color, and the healthy tissues had a red color. The degree of myocardial infarction was represented by the percentage of infarcted area accounting for the left ventricle.

Tissue extraction and RNA isolation

In total, three rats in the MIRI group and three from the sham group were anesthetized, and the T1-T4 spinal cord segments were quickly extracted. Total RNAs were collected from the dorsal horn of the spinal cord using TRIzol® reagent (Invitrogen; Thermo Fisher Scientific, Inc.) based on the manufacturer's instructions. The RNA concentration and purity were examined using a spectrophotometer (Invitrogen; Thermo Fisher Scientific, Inc.) (32,33). Equal amounts of mRNA (1 µg) from one rat with the same treatment was mixed as one sample after these tests. A total of six samples were selected for microarray analysis.

Gene expression analysis

High-throughput lncRNA-seq was used to identify novel genes in the spinal cord that regulate cardiac function in the rat MIRI model. Briefly, total RNA was isolated from the T1-T4 of spinal cord tissues from the three MIRI- and sham-group rats, and the concentration, purity and integrity of the RNA were assessed on an Agilent 2100 Bioanalyzer (Agilent Technologies, Inc.) and checked using RNase-free 1% agarose gel electrophoresis containing GelRed® (cat. no. 41003; Biotium, Inc.). Ribosomal RNA was then removed from the total RNA. The mRNA and ncRNA were mostly retained and the strand-specific library was then constructed using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (cat. no. E7530L; New England BioLabs, Inc.) according to the manufacturer's protocols. Briefly, the mRNA was fragmented and used as a template to synthesize cDNA with random primers. The first strand cDNA was synthesized using the following progress: 25˚C for 10 min, 42˚C for 15 min, 70˚C for 15 min and final holding at 4˚C. A Second Strand Marking Master Mix was used to synthesize second strand cDNA at 16˚C for 60 min. After using a ligate adapter, uracil-N-glycosylase treatment and PCR amplification were used to construct a sequencing library. The loading concentration of DNA was 12.182 nM. The PCR thermocycling parameters were set as follows: Initial denaturation at 98˚C for 30 sec; followed by 15 cycles of 98˚C for 10 sec, 60˚C for 30 sec and 72˚C for 30 sec; 1 cycle at 72˚C for 5 min and final holding at 10˚C. The library was then sequenced by Illumina HiSeqÔ 4000 with pair ends reads of 150 bp. The data were analyzed by Gene Denovo Biotechnology Co. (https://www.genedenovo.com; Guangzhou, China). The differentially expressed mRNAs and lncRNAs were identified with |log2 fold change (FC)| and P-values as calculated with a t-test. |log2FC| >1 and P<0.05 were set as the up- and downregulated gene thresholds. Hierarchical clustering and volcano plots determined the different expression patterns of mRNAs and lncRNAs among the samples.

Bioinformatics analysis

The differentially expressed lncRNAs and mRNAs with statistical significance were identified via volcano plot filtering. The threshold used to screen RNAs with |log2FC| >1 was P<0.05. |log2FC| >1 is to screen genes with multiples of difference >2 times (34,35). Gene Ontology (GO) analysis (http://www.geneontology.org) described the molecular function, cellular component and biological process of differentially expressed transcripts. GO analysis and a Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were applied to determine the biologic roles of these differentially expressed lncRNAs and mRNAs based on the latest KEGG data (http://www.genome.jp/kegg/).

RT-qPCR analysis for upper thoracic spinal cord

The total RNA, extracted from the upper T1-T4 segments of the spinal cord using the TRIzol® reagent (Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions, was used for the generation of cDNA (21,22,36-40). The primers were designed using the Primer Express 3.0 software (Applied Biosystems; Thermo Fisher Scientific, Inc.) and the specific forward and reverse primer sequences are listed in Table I. RNA samples were quantified using a spectrophotometer (BioPhotometer; Eppendorf AG) and then synthesized to cDNA by reverse transcription using the PrimeScript™ RT reagent kit (Takara Bio, Inc.). The temperature protocol was conducted using the following protocol: 15 min at 37˚C, 5 sec at 85˚C and holding at 4˚C cDNA was quantitated via RT-qPCR using a TB green® Premix Ex Taq (cat. no. RR420A; Takara Bio, Inc.). The thermocycling conditions for PCR were as follows: Initial denaturation for 30 sec at 95˚C, followed by 40 cycles of 15 sec at 95˚C, 15 sec at 60˚C and 45 sec at 72˚C. The threshold cycle (Cq) was used to estimate the amount of target mRNA. The comparative Cq method with a formula for relative FC (2-ΔΔCq) was used to quantify the amplified transcripts (32,33,41). The relative gene expression was determined via normalization to GAPDH. Experiments were evaluated in triplicate and repeated ≥3 times.

Table I

Primer sequences for reverse transcription-quantitative PCR.

A, mRNA
Primer name	Primer sequences (5'→3')
Frs3	F: GAGCCAGTCATCATCACACGGAAC
	R: GAAGCCATTGGAGAAGCTGGAGAC
Zfp523	F: CGCTGACTTCCTGTGAGTGTGAC
	R: ACCTTGGCTCTGGCTCTCGTC
Dnajc6	F: GACAAGCCTCATGGAGCCAAGAAG
	R: GGAGAAGTTCGTGCTCACATCGG
Nedd4l	F: TACGCAGTGGCACCGACCTAG
	R: GTGTGCGGCCTCCTGATTGATC
Tep1	F: GGATGGATACGGAGCTGCTGAATG
	R: GGACACTGACGAAGAGGCACTTG
Myef2	F: CCAGGTGGACAGCCAATTAGTGC
	R: GCCTGCCATGCTATTCATTGCTTC
Tgfbr1	F: CATTGCTGGTCCAGTCTGCTTCG
	R: TGGTGAATGACAGTGCGGTTATGG
Fgf12	F: TGCTGTGCGACTGTGAGATTGC
	R: CGTGCGGCTCGTTGTACTCG
Mef2c	F: TGCTGTGCGACTGTGAGATTGC
	R: CGTGCGGCTCGTTGTACTCG
Tfdp1	F: AAGAAGATGCCGCCGAGTTGTG
	R: CGCTGGCTTCACGACACATCC
GAPDH	F: CGCTAACATCAAATGGGGTG
	R: TTGCTGACAATCTTGAGGGAG
B, Long non-coding RNA
Primer name	Primer sequences (5'→3')
ENSRNOT00000080713	F: CACTTCCGCCCACATCACA
	R: AAACCCCTCAGTCTCAAAACCT
ENSRNOT00000090564	F: GACCACGCATCTACA
	R: GCGAACTGAGGTGAGAAGGGA
ENSRNOT00000082588	F: TGCTTCACTTCTCCCACTCTTC
	R: AGGCCTTCCATCAAACAAATC
ENSRNOT00000091080	F: GGTGGTGGACAAGGATGAG
	R: TGGACGGAGGTTGGGGAGAT
ENSRNOT00000091570	F: CAAGAAATCCAATACGACCC
	R: TCCCCGTGGAAATAGGTGTAA
ENSRNOT00000087777	F: GGAAGAGGGACGTTGGGTA
	R: AGCCAGGTCCAAGTCACAGG
ENSRNOT00000082061	F: CATTTCTTCCTCCCTTCCCTT
	R: ACGGTTCTCCAATCGCACA
ENSRNOT00000091108	F: CCGTGGAATATGGGAAGGAC
	R: CAAGATGGAAAGGCAACGAG
ENSRNOT00000087028	F: ACGGGGTGAAAGGGAAGA
	R: TGCTGGTTTGTTAAGGGGATG
ENSRNOT00000086475	F: CAGTTTGGTGTCTGTTTCCC
	R: TGCTTCTTCGAGCCGGTATT

F, forward; R, reverse.

Statistical analysis

Data are presented as the mean ± SEM, and were analyzed using GraphPad Prism software v5.0 (GraphPad Software, Inc.). RT-qPCR parameters were analyzed using unpaired t-test for repeated measures. P<0.05 was considered to indicate a statistically significant difference.

Results

Characteristics of ischemic myocardial tissues

Serum cTnI levels in the model group were significantly increased compared with the control group (Fig. 1B). In addition, the results of TTC staining (Fig. 1A), H&E staining and Masson staining (Fig. 1C) all confirmed the successful establishment of the rat MIRI model. In the control group, the myocardium had a physiological myocardial architecture, and the myocardial fibres and myocardial cells were complete and arranged in an orderly manner. In the model group, structural disorder of the cardiac tissue was observed, with different degrees of vacuolar degeneration and necrosis, as well as loose stroma. Moreover, the number of cardiomyocyte fibres were markedly increased after MIRI.

Figure 1

Evaluation of MIRI. (A) Representative images of 2,3,5-triphenyltetrazolium chloride staining (left) and quantitative analysis of infarct size (right) in hearts subjected to Sham and MIRI surgery. Scale bar, 2 mm. (B) Serum cTnI concentrations in rats with Sham or MIRI surgery. Data are presented as the mean ± SEM, n=4 rats per group. (C) Representative images of H&E staining and Masson trichrome staining from rat left ventricles 2 h after MIRI and in sham group. Scale bar, 100 µm. *P<0.05 vs. sham group. H&E, hematoxylin and eosin; MIRI, myocardial ischemia-reperfusion injury.

mRNAs and lncRNAs expression profiles in MIRI

The patterns of the mRNA and lncRNA expressions in the T1-T4 spinal cord at 2 h after MIRI or sham operations were then examined using microarray. From the volcano map and hierarchical clustering analysis results, a landscape of the expression characteristics of the mRNAs and lncRNAs was obtained. The volcano plots demonstrated that large numbers of the mRNAs and lncRNAs were differentially expressed between the two groups (Fig. 2A and C). The hierarchical cluster analysis of the lncRNAs and mRNAs indicated that the samples of two groups were clustered together and the signal intensity was consistent in the two groups (Fig. 2B and D). Furthermore, these differential alterations of the mRNAs and lncRNAs in the T1-T4 spinal cord were associated with MIRI.

Figure 2

MIRI results in expression profile changes of mRNA and lncRNA in the T1-T4 spinal cord. (A) Volcano plot comparing global mRNA gene expression profiles in the spinal cord between the two groups. (B) Hierarchical cluster analysis of the mRNAs. (C) Volcano plot comparing global lncRNA gene expression profiles in the spinal cord between the two groups. (D) Hierarchical cluster analysis of the lncRNAs. The cluster analysis demonstrated that the three MIRI or three sham samples were clustered together. Upregulated, downregulated and not statistically different genes were colored in red, green and black, respectively. MIRI, myocardial ischemia-reperfusion injury; lncRNA, long non-coding RNA.

Differentially expressed mRNAs and lncRNAs

To further analyze the differentially expressed mRNAs and lncRNAs, a significance analysis of the microarrays method with the criteria P<0.05 and FC >1 was performed.

In the differentially expressed mRNAs, there were 470 genes that exhibited an FC that was >1. The number of downregulated mRNAs was 239, whereas the number of upregulated mRNAs was 231. The most upregulated mRNAs were Tfdp1, Map2, Cep19, Aurkaip1, Tpbg, Myef2, Ap2m1, Tgfbr1, Fgf12 and Mef2c. The most downregulated mRNAs were Frs3, Zfp523, Tpm4, Apold1, Ctnnb1, Lss, Psen2, Dnajc6, Nedd4l and Tep1. Detailed information regarding the differentially expressed mRNAs is listed in Table II.

Table II

Details of the top 10 downregulated and upregulated mRNAs in the T1-T4 spinal cord between MIRI and sham groups.

A, Downregulation
ID	Log2FC (MIRI/sham)	Symbol	P-value	Description
ENSRNOT00000019404	-12.0846	Frs3	4.50x10-20	Fibroblast growth factor receptor substrate 3
ENSRNOT00000000595	-11.2761	Zfp523	1.51x10-14	Zinc finger protein 523
ENSRNOT00000089056	-10.8893	Tpm4	0.001717	Tropomyosin 4
ENSRNOT00000010289	-10.5635	Apold1	1.1610-9	Apolipoprotein L domain containing 1
ENSRNOT00000079085	-10.1716	Ctnnb1	8.32x10-5	Catenin beta 1
ENSRNOT00000081903	-9.81911	Lss	0.013918	Lanosterol synthase (2,3-oxido squalene-lanos terol cyclase)
ENSRNOT00000091715	-9.75933	Psen2	0.000372	Presenilin 2
ENSRNOT00000086539	-9.30682	Dnajc6	5.25x10-5	DnaJ heat shock protein family (Hsp40) member C6
ENSRNOT00000083363	-9.15482	Nedd4l	0.005334	Neural precursor cell expressed, developmentally down-regulated 4-like, E3 ubiquitin protein ligase
ENSRNOT00000088981	-8.70275	Tep1	0.001598	Telomerase associated protein 1
B, Upregulation
ID	Log2FC (MIRI/sham)	Symbol	P-value	Description
ENSRNOT00000087192	12.22179	Tfdp1	6.42x10-8	Transcription factor Dp-1
ENSRNOT00000045766	11.86754	Map2	5.07x10-5	Microtubule-associated protein 2
ENSRNOT00000037795	11.22882	Cep19	2.22x10-5	Centrosomal protein 19
ENSRNOT00000025531	11.16113	Aurkaip1	0.000724	Aurora kinase A interacting protein 1
ENSRNOT00000014326	10.76763	Tpbg	0.000962	Trophoblast glycoprotein
ENSRNOT00000081533	10.69407	Myef2	7.98x10-9	Myelin expression factor 2
ENSRNOT00000088821	10.64386	Ap2m1	0.004070	Adaptor-related protein complex 2, mu 1 subunit
ENSRNOT00000081090	9.813781	Tgfbr1	0.005577	Transforming growth factor, beta receptor 1
ENSRNOT00000071586	9.781360	Fgf12	0.004286	Fibroblast growth factor 12
ENSRNOT00000076230	8.748193	Mef2c	1.98x10-5	Myocyte enhancer factor 2C

FC, fold change; T, thoracic; MIRI, Myocardial ischemia-reperfusion injury.

The results identified that 126 lncRNAs, which included 41 upregulated and 85 downregulated lncRNAs, were significantly altered in the MIRI group compared with the sham group. The most upregulated lncRNAs were ENSRNOT00000086475, ENSRNOT00000087028, ENSRNOT00000091108, ENSRNOT00000031815, ENSRNOT00000082061 and ENSRNOT00000087777. The most downregulated lncRNAs were ENSRNOT00000080713, ENSRNOT00000090564, ENSRNOT00000082588, ENSRNOT00000091080, ENSRNOT00000091570, ENSRNOT00000077319, ENSRNOT00000077509, ENSRNOT00000089288, ENSRNOT00000083187 and ENSRNOT00000091435. Additional information regarding the differentially expressed lncRNAs is presented in Table III.

Table III

Details of the downregulated and upregulated long non-coding RNAs in the T1-T4 spinal cord between MIRI and sham group.

A, Downregulation
ID	Log2FC (MIRI/sham)	P-value	Symbol
ENSRNOT00000080713	-9.98726	0.000442	LOC102550026
ENSRNOT00000090564	-8.10329	0.040576	AABR07050135.1
ENSRNOT00000082588	-7.49185	4.01x10-5	AABR07000385.2
ENSRNOT00000091080	-6.9542	1.14x10-6	AABR07044980.1
ENSRNOT00000091570	-6.71425	0.025341	AABR07071198.1
ENSRNOT00000077319	-6.22882	0.005841	AABR07071383.1
ENSRNOT00000077509	-5.32193	0.048185	AABR07027371.1
ENSRNOT00000089288	-4.16993	0.01174	AABR07044001.4
ENSRNOT00000083187	-3.52356	0.024012	AABR07060560.4
ENSRNOT00000091435	-3.38082	0.016282	AABR07042668.1
B, Upregulation
ID	Log2FC (MIRI/sham)	P-value	Symbol
ENSRNOT00000086475	8.797662	0.013366	AABR07070801.2
ENSRNOT00000087028	6.129283	0.027308	LOC103692514
ENSRNOT00000091108	4.058894	0.020796	AABR07006724.1
ENSRNOT00000031815	3.297681	0.041965	Rn50_X_0744.6
ENSRNOT00000082061	2.966833	0.034871	AABR07071395.1
ENSRNOT00000087777	2.917538	0.024331	AABR07016845.1

FC, fold change; T, thoracic; MIRI, Myocardial ischemia-reperfusion injury.

Functional prediction of differentially expressed mRNAs in MIRI

To investigate the spinal molecular mechanisms in MIRI, GO and KEGG pathway analyses of the differentially expressed mRNAs in MIRI vs. sham were performed.

To identify major biochemical metabolic pathways and signal transduction pathways for differential gene enrichment, KEGG pathway analysis (Fig. 3) was conducted. Using the KEGG pathway analysis, it was found that differential genes were mainly enriched in the ‘PI3K/Akt signaling pathway’, ‘protein digestion and uptake’, ‘p53 signaling pathway’ and ‘neuroactive ligand-dependent body interaction’ (Fig. 3; Table IV).

Figure 3

Heatmap of Kyoto Encyclopedia of Genes and Genomes pathway analysis. Size of the bubble represents the number of enriched genes, and the color of the bubble from red to green indicates a gradual increase in the P-value. ECM, extracellular matrix; TRP, transient receptor potential.

Table IV

Differential gene enriched in Kyoto Encyclopedia of Genes and Genomes pathway.

Pathway ID	Pathway	Gene number	P-value
ko04151	PI3K/Akt signaling pathway	20	0.00039
ko04974	Protein digestion and absorption	8	0.001448
ko04115	p53 signaling pathway	7	0.001795
ko04510	Focal adhesion	13	0.005058
ko04080	Neuroactive ligand-receptor interaction	16	0.005635
ko04512	ECM-receptor interaction	7	0.008425
ko00534	Glycosaminoglycan biosynthesis-heparan sulfate/heparin	3	0.019966
ko04960	Aldosterone-regulated sodium reabsorption	4	0.023467
ko04015	Rap1 signaling pathway	12	0.027612
ko04721	Synaptic vesicle cycle	5	0.030579
ko00600	Sphingolipid metabolism	4	0.036733
ko04020	Calcium signaling pathway	11	0.03857

ECM, extracellular matrix.

The differential genes associated with cellular components were mainly enriched in the ‘extracellular fraction’, ‘cytoplasm’, ‘microtubule cytoskeleton’, ‘cytoskeleton’ and ‘intracellular’, whereas the genes participating in biological processes were mainly enriched in ‘cell signal transduction’ (Fig. 4; Table V).

Figure 4

Functional analysis of differential lncRNAs. (A) Gene Ontology analysis heat map of co-expressed genes of differential lncRNAs. (B) Kyoto Encyclopedia of Genes and Genomes pathway analysis heat map of co-expressed genes of differential lncRNAs. lncRNA, long non-coding RNA; PPAR, peroxisome proliferator activated receptor; NOD, nucleotide-binding oligomerization domain-containing protein.

Table V

Functional classification of differential genes via GO analysis.

A, Cellular component
GO ID	Description	Gene number	P-value
GO:0044424	Intracellular part	48	0.013126
GO:0005737	Cytoplasm	30	0.025862
GO:0015630	Microtubule cytoskeleton	3	0.026874
GO:0005856	Cytoskeleton	7	0.041896
GO:0005622	Intracellular	63	0.044992
B, Biological process
GO ID	Description	Gene number	P-value
GO:0030198	Extracellular matrix organization	3	0.010728
GO:0043062	Extracellular structure organization	3	0.010729
GO:0007267	Cell-cell signaling	4	0.028507
GO:0032196	Transposition	1	0.035381

GO, Gene Ontology.

Gene co-expression networks were presented to identify interactions between lncRNAs and their co-expression genes (Fig. 5). The co-expressed network showed that 10 highly-dysregulated lncRNAs were highly co-expressed with 277 target genes, of which 10 were positively correlated (colored with yellow) and 267 were negatively correlated (colored with blue). ENSRNOT00000080713 was at the core of the lncRNA-mRNA co-expression network and had a co-expression relationship with many mRNAs.

Figure 5

Sub-network of 10 highly dysregulated lncRNAs co-expressed with mRNAs. Genes colored in red are upregulated lncRNAs and those in green are downregulated lncRNAs. Upregulated mRNAs are colored with yellow and downregulated mRNAs are colored with blue. lncRNA, long non-coding RNA.

RT-qPCR confirmation of mRNA and lncRNA expression from microarray

To confirm the reliability of the microarray data and analyze the temporal alterations in mRNA and lncRNA expression after MIRI, the lncRNAs, ENSRNOT00000080713, ENSRNOT00000090564, ENSRNOT00000082588, ENSRNOT00000091080, ENSRNOT00000091570, ENSRNOT00000087777, ENSRNOT00000082061, ENSRNOT00000091108, ENSRNOT00000087028 and ENSRNOT00000086475, and the mRNAs, fibroblast growth factor receptor substrate 3 (Frs3), zinc finger protein 76 (Zfp523), DnaJ heat shock protein family (Hsp40) member C6 (Dnajc6), NEDD4 like E3 ubiquitin protein ligase (Nedd4l), telomerase associated protein 1 (Tep1), myelin expression factor 2 (Myef2), transforming growth factor β receptor 1 (Tgfbr1), fibroblast growth factor 12 (Fgf12), myocyte enhancer factor 2C (Mef2c) and transcription factor Dp-1 (Tfdp1), were randomly selected and detected via RT-qPCR.

The results demonstrated that expression of lncRNAs ENSRNOT00000080713, ENSRNOT00000090564, ENSRNOT00000082588, ENSRNOT00000091080 and ENSRNOT00000091570 were markedly decreased in MIRI group compared with those in the control group, whilst ENSRNOT00000087777, ENSRNOT00000082061, ENSRNOT00000091108, ENSRNOT00000087028 and ENSRNOT00000086475 were markedly increased. In addition, the mRNAs Frs3, Zfp523, Dnajc6, Nedd4l and Tep1 were markedly decreased in the MIRI group compared with those in the control group, but Myef2, Tgfbr1, Fgf12, Mef2c and Tfdp1 were markedly increased. These RT-qPCR results were consistent with the results from the transcriptome sequencing (Fig. 6), which further supported the reliability of the microarray data.

Figure 6

RT-qPCR confirmation of mRNA and lncRNA expression from the microarray. (A) mRNA expression levels in spinal cord tissue of MIRI rats. (B) lncRNA expression levels in spinal cord tissue of MIRI rats. The expression of mRNAs and lncRNAs in the spinal cord of rats with MIRI was consistent with the results of transcriptome sequencing, indicating the reliability of the microarray data. lncRNA, long non-coding RNA; Frs3, fibroblast growth factor receptor substrate 3; Zfp523, zinc finger protein 76; Dnajc6, DnaJ heat shock protein family (Hsp40) member C6; Nedd4l, NEDD4 like E3 ubiquitin protein ligase; Tep1, telomerase associated protein 1; Myef2, myelin expression factor 2; Tgfbr1, transforming growth factor β receptor 1; Fgf12, fibroblast growth factor 12; Mef2c, myocyte enhancer factor 2C; Tfdp1, transcription factor Dp-1; RT-qPCR, reverse transcription-quantitative PCR; seq, sequencing; MIRI, myocardial ischemia-reperfusion injury.

RT-qPCR validation of mRNAs and lncRNAs in T1-T4 spinal cord

The results of RT-qPCR demonstrated that the mRNA expression levels of Frs3, Nedd4l and Tep1 were significantly downregulated in the model group, while the expression of mRNA Myef2 was significantly upregulated in the model group. The expression levels of mRNA Zfp523, Dnajc6 and Myef2 were not significantly different between the two groups (Fig. 7). The expression levels of lncRNA ENSRNOT00000080713, ENSRNOT00000090564, ENSRNOT00000082588 and ENSRNOT00000091570 were significantly downregulated in the model group, but the expression levels of lncRNA ENSRNOT00000087777, ENSRNOT00000082061, ENSRNOT00000087028 and ENSRNOT00000086475 were significantly upregulated in the model group. However, the expression levels of lncRNA ENSRNOT00000091080 and ENSRNOT00000091108 were not significantly different between the two groups (Fig. 8).

Figure 7

Reverse transcription-quantitative PCR validation of the mRNAs in the T1-T4 spinal cord. The expression levels of mRNA Frs3, Nedd4l and Tep1 were significantly downregulated in model group, whereas the expression of mRNA Myef2 was significantly upregulated. The expression levels of mRNA Zfp523, Dnajc6 and Myef2 were not statistically different between the two groups. Data are presented as the mean ± SEM (n=6). *P<0.05 vs. control. Frs3, fibroblast growth factor receptor substrate 3; Zfp523, zinc finger protein 76; Dnajc6, DnaJ heat shock protein family (Hsp40) member C6; Nedd4l, NEDD4 like E3 ubiquitin protein ligase; Tep1, telomerase associated protein 1; Myef2, myelin expression factor 2; Tgfbr1, transforming growth factor β receptor 1; Fgf12, fibroblast growth factor 12; Mef2c, myocyte enhancer factor 2C; Tfdp1, transcription factor Dp-1; MIRI, myocardial ischemia-reperfusion injury.

Figure 8

Reverse transcription-quantitative PCR validation of 10 lncRNAs in the T1-T4 spinal cord. The expression levels of lncRNA ENSRNOT00000080713, ENSRNOT00000090564, ENSRNOT00000082588 and ENSRNOT00000091570 were significantly downregulated in model group, whereas the expression levels of lncRNA ENSRNOT00000087777, ENSRNOT00000082061, ENSRNOT00000087028 and ENSRNOT00000086475 were significantly upregulated. The expression levels of lncRNA ENSRNOT00000091080 and ENSRNOT00000091108 were not statistically different between the two groups. Data are presented as the mean ± SEM (n=6). *P<0.05, **P<0.01,***P<0.001 vs. control. MIRI, myocardial ischemia-reperfusion injury; lncRNA, long non-coding RNA.

Discussion

MIRI causes nociceptive chemical and electrical signals that affect dorsal root ganglia and the upper thoracic spinal cord (42). Previous studies have reported that spinal gene expression is associated with in the development of MIRI (9,27) and have identified the molecular mechanisms during MIRI (26,27,42); however, the impact of the spinal process during MIRI requires further investigation.

Emerging evidence indicates that the inflammatory and immune responses that occur in the spinal cord are vital contributors to the progression of MIRI (9,13,14,16). Moreover, autonomic dysregulation following MIRI is an important pathogenic event (43-47). Interactions between the heart and spinal autonomic nervous system serve a role in both the development and maintenance of MIRI, and are also crucial in myocardial ischemia-induced angina pectoris (48-52). Based on the importance of the spinal autonomic nervous system in the pathophysiology of MIRI, there has been a notable effort in this field to study changes in the genomics, transcriptomics, proteomics and metabolomics in the spinal cord after MIRI, which could facilitate both existing interventions and novel therapies directed at cardiac-autonomic targets (53).

In the present study, the differentially expressed lncRNAs in the T1-T4 spinal cord during MIRI, along with their characteristics and possible relationships with mRNAs, were identified. A total of 126 lncRNAs were observed in the T1-T4 spinal cord of the MIRI model rats. Among these lncRNAs, 41 lncRNAs were upregulated, and 85 lncRNAs were downregulated at 2 h after MIRI. Moreover, RT-qPCR was used to validate the microarray analysis results, for which the conclusion was consistent. Based on these findings, it was suggested that the aforementioned lncRNAs may serve a role in MIRI pathology. Furthermore, the present research contributes to understanding the molecular mechanism of MIRI.

To further investigate the roles of these differentially expressed lncRNAs in the development of MIRI, GO and KEGG pathway analyses were performed. Based on the GO and KEGG enrichment analyses of these lncRNAs, the results suggested that the significantly enriched biologic processes and molecular functions of the upregulated genes after MIRI were associated with gene sets termed as follows: ‘PI3K/Akt signaling pathway’, ‘protein digestion and absorption’, ‘p53 signaling pathway’ and ‘neuroactive ligand-dependent body interaction’.

The PI3K/Akt signaling pathway is associated with MIRI (54). p53 is a tumor suppressor, which exerts an important role in the process of apoptosis and is associated with the occurrence and development of cardiovascular disease (55-57). Previous studies have shown that apoptosis is a typical manifestation and primary pathological mechanism of MIRI (58). The loss of cardiomyocytes caused by apoptosis is an important reason for the development of various heart diseases (59). Preventing the process of apoptosis can prevent myocardial damage caused by reperfusion, thus slowing or preventing the occurrence of heart failure and even mortality (60).

The present study has several limitations. For example, the reperfusion time for tissue collection lasted only 2 h. Therefore, it remains unknown if the expression of lncRNA has the same regulation for a longer period after ischemia. Moreover, a functional analysis of differentially expressed genes was performed without further verification. LncRNAs are poorly conserved and only part of the sequence is conserved by multiple species (61,62). Previous studies have reported that there are 481 segments longer >200 bp (range, 200-779 bp) that are conserved between orthologous regions of the mouse, rat and human genomes (63), and these conserved lncRNAs can provide an insight into their key functional roles (64,65). As a result, further research is required to examine the functions of lncRNAs and the conservation of differentially expressed lncRNAs in rats and humans, which will identify their predicted targets and provide a detailed expression profile of lncRNA.

In conclusion, the present study demonstrated that lncRNA transcripts of the upper thoracic spinal cord were differentially expressed after MIRI. These spinal lncRNAs may regulate the expression of MIRI-related proteins and genes, and contribute to the pathogenesis of MIRI. Future investigations of these spinal lncRNAs found in the present study will mainly focus on their functions and their relationship with MIRI. Nonetheless, the present study may provide important evidence for further basic research and clinical investigations on improved therapeutic interventions of MIRI.