Identification of a hippocampal lncRNA-regulating network in cognitive dysfunction caused by chronic cerebral hypoperfusion

CCH-induced dysfunction of spatial learning and memory in the Morris water maze test

The results of the Morris water maze test showed that CCH rats had a much longer latency time to reach the platform from the third to seventh training day than the sham control rats (P < 0.01) (Figure 1A). CCH rats spent noticeably less time in the platform area from the fourth to seventh training day than the sham control rats (P < 0.01) (Figure 1B). After 7 days of training and 1 day of rest, the short-term memory test showed that the CCH rats had a dramatically longer latency time to reach the platform than sham rats (P < 0.01) (Figure 1C). After removing the platform, CCH rats spent less time in the platform area and less crossing time around the platform area than the sham control rats (P < 0.01) (Figure 1D–1E).

No obvious number changes of neurons in hippocampus and striatum after chronic cerebral hypoperfusion

In order to explore the potential causes of spatial cognitive dysfunction caused by chronic cerebral hypoperfusion, brain sections were stained by HE staining and the changes of morphology and number of neurons in the sections were observed. The results showed that there was no significant change in the number of neurons in CA3, CA1, DG and striatum (Figure 2A, 2B). This suggested that the cognitive dysfunction caused by chronic cerebral hypoperfusion may not be due to the neurons decrease of in hippocampus and striatum. Scale bar=50μm.

Expression profile of lncRNA and mRNA profiles in the CCH rat hippocampus

To investigate the underlying mechanisms of cognitive deficit caused by chronic cerebral hypoperfusion, we employed the RNA sequencing to explore the mRNA and lncRNA changes of hippocampus closely related with cognition after chronic cerebral hypoperfusion. Through the assessment of total RNA quantity (≥13μg), the sample concentration (≥433μg/μL), the purity of samples(OD260/280 ≥1.82) and the sample integrity(RIN≥8.5), all extracted RNAs from hippocampus tissue completely meet the requirement of library construct. After high throughput RNA sequencing, 654.12M raw reads and 98.11G raw bases were acquired, containing 340.81M reads and 51.12G for CCH models, and 315.31M reads and 46.99G for the sham controls. Through low quality sequencing data cleaning of the original data, we finally obtained 632.98 M clean reads and 89.54G clean bases in total, containing 329.81M reads and 46.64G bases in the CCH models, and 303.17M reads and 42.9 bases in the sham controls. The proportions of valid bases were 90.75%-91.69% and the Q30s were 94.19~94.68%. The average guanine-cytosine (GC) content was 48.07%.

A total of 22601 protein-encoding transcripts and 15460 lncRNAs were identified and subsequently analyzed in-depth. For lncRNA, 1371 novel lncRNAs were identified and subsequently analyzed. These have not been reported in the past. The sequence length of lncRNA transcripts ranged from 200–4300 bp. A total of 236 mRNAs and 189 lncRNAs were significantly altered in the CCH rat hippocampus compared to sham controls. Among these, 128 mRNAs and 91 lncRNAs were up-regulated, while 108 mRNAs and 98 lncRNAs were down-regulated.

The most up-regulated mRNAs and lncRNAs were Sftpa1 and TCONS_00019175, with FCs of 79.54 and 689.55, respectively compared to sham controls. The most down-regulated mRNAs and lncRNAs were Cdh17 and ENSRNOT00000092488, with FCs of 13.38 and 182.18, respectively compared to sham controls. The top 20 up-regulated and 20 down-regulated lncRNAs and mRNAs in the CCH group are listed in Tables 1 and 2. The clustering analysis and volcano plot visualization showed dramatically different expression levels of mRNAs and lncRNAs in the CCH and sham groups (Figures 3 and 4).

Expression profile validation

In order to verify the validity of RNA sequencing, we randomly selected two differentially up-regulated mRNAs of Katnal1 and Trpv4 and down-regulated mRNAs of Mgam and Cpg1 and two differentially up-regulated lncRNAs of XR_362089.3 and TCONS_00007808 and down-regulated lncRNAs of XR_001838876.1 and XR_595897.2 for detection by qRT-PCR. The qRT-PCR results showed that the change trends of the selected mRNA and lncRNA levels determined by qRT-PCR were consistent with these by RNA sequencing. This showed that the RNA sequencing data was valid (Figure 5A). The levels of Katnal1 and Trpv4 in hippocampus of CCH group were significantly up-regulated than in that of sham group (P<0.01), whereas Mgam and Cpg1 in hippocampus of CCH group were significantly down-regulated than in that of sham group (P<0.01) (Figure 5B). Similarly, up-regulation lncRNAs of XR_362089.3 and TCONS_00007808 and down-regulation lncRNAs of XR_001838876.1 and XR_595897.2 were noticeable in CCH group than in sham group (P<0.01) (Figure 5B). The primers used in qRT-PCR are listed in Table 3.

GO function and KEGG pathway enrichment analyses of DE mRNAs

The GO analysis predicted, and the KEGG analysis showed, that DE RNAs were associated with the probable functions and pathways. GO analysis indicated that the most enriched mRNAs related to monocarboxylic acid transport (GO: 0015718) in biological processes (Figure 6A), the extracellular space (GO: 0005615) in the cellular component (Figure 6B), and the peptide antigen binding (GO: 0042605) in the molecular functions (Figure 6C). KEGG pathway analysis showed that the top 10 deferentially-enriched KEGG pathways related to dysregulated mRNAs were involved in tryptophan metabolism, graft-versus-host disease, allograft rejection, thyroid hormone synthesis, type I diabetes mellitus, autoimmune thyroid disease, the intestinal immune network for IgA production, staphylococcus aureus infection, cell adhesion molecules (CAMs), and viral myocarditis (Figure 6D).

LncRNA-mRNA co-expression analyses

In total, 161 DE lncRNAs that were correlated with DE mRNAs were identified (P < 0.05 and correlation coefficient (COR) > 0.8). Among of all the lncRNA-mRNA coexpressions, lncRNA XR_589889.2 and Tgm2 showed p = 0.0083 and COR = 0.92 and lncRNA XR_595556.2 and Cdhr3 showed P= 0.018 and COR = 0.89. Ten mRNAs and their coexpressed lncRNAs were selected to construct the expression network (Figure 7A–7C).

Construction of lncRNA-miRNA-mRNA regulatory network

186 lncRNA-miRNA regulatory pairs were identified including 108 lncRNAs and mRNA. 368 miRNA-mRNA target pairs were identified including 147 miRNAs and 130 mRNAs. The 67 common miRNAs were identified (Figure 8A). Because of competitively binding miRNA as a miRNA sponge, lncRNAs could form a ceRNA network of lncRNA-miRNA-mRNA to booster miRNA target genes. Based on the regulatory pairs of miRNA-mRNA and lncRNA-miRNA, an lncRNA-miRNA-mRNA network was constructed. A total of 559 lncRNA-miRNA-mRNA target pairs were identified, including 83 lncRNAs, 67 miRNAs, and 108 mRNAs (Figure 8B).

Validation of lncRNA-miRNA-mRNA regulatory network

To validate the ceRNAs levels of lncRNA-miRNA-mRNA regulatory network, four random lncRNA-miRNA-mRNA target pairs, Itgb6-miR-1247-XR_001836753.1, Npsr1-miR-3594-XR_597317.2, Tsks-miR-667-3p-TCONS_00009790, Cyb5rl-miR-330-XR_353865.3, were selected for qRT-PCR. The results showed that Itgb6 and Npsr1 levels were significantly upregulated and Tsks and Cyb5rl levels were significantly downregulated in hippocampus of CCH animals than that of sham animals (P<0.01) (Figure 9A). As predicted by lncRNA-miRNA-mRNA regulatory network construction, microRNAs rno-miR-1247-5p and rno-miR-3594-5p were significantly down-regulated, and rno-miR-667-3p and rno-miR-330-5p were significantly upregulated in hippocampus of CCH animals than that of sham animals (P<0.01) (Figure 9B). Similarly, lncRNA XR_001836753.1 and XR_597317.2 were significantly upregulated, TCONS_00009790 and XR_353865.3 were significantly down-regulated in hippocampus of CCH animals than that of sham animals (P<0.01) (Figure 9C). The primers used in qRT-PCR are listed in Table 4.

Enrichment analyses of the DE lncRNA-miRNA-DE mRNA network

To further predict and analyze the probable functions and pathway of DE mRNAs regulated by DE lncRNA-miRNA, GO and KEGG analysis were performed. The GO analysis showed that the most enriched mRNAs related to lipid homeostasis (GO: 0055088) in biological processes (Figure 10A), the multivesicular body (GO: 0005771) in the cellular component, and hormone activity (GO: 0005179) in the molecular functions (Figure 10B, 10C). KEGG pathway analysis showed that the top 10 deferentially-enriched KEGG pathways related to dysregulated mRNAs were involved in tryptophan metabolism, thyroid hormone synthesis, pertussis, cholinergic synapse, osteoclast differentiation, the Jak-STAT signaling pathway, phagosome, the PI3K-Akt signaling pathway, proteoglycans in cancer, and cytokine-cytokine receptor interaction (Figure 10D).

Animals and CCH model surgery

Adult male Sprague-Dawley rats (200–220 g) were obtained from Hunan SJA Laboratory Animal Co., Ltd., and were housed with accessible food and water ad libitum. Rats were kept on a 12-h light/dark cycle with the light on from 7:00 am to 7:00 pm. Animal welfare and all experiments were approved by Ethics Committee of Renmin Hospital of Wuhan University.

The CCH model surgery was performed as previously described [12]. Briefly, the rats were anesthetized with intraperitoneal chloral hydrate (0.4 g/kg) and placed with a heating pad to maintain the body temperature at 37°C. After a ventral midline incision, both common carotid arteries were gently separated from the carotid sheath and vagus nerve [41]. Bilateral common carotid arteries were doubly ligated with 4-0 silk sutures just below the carotid bifurcation. In control rats, similar surgery was performed, but without vessel ligation. After the surgery was finished, the temperature of the rats was maintained under 37°C until recovery. During surgery, the brains of the rats did not soften, wholly or in part.

Morris water maze

After 30 days of CCH, all rats completed spatial memory training in the Morris water maze. The experiment was conducted as previously described [39]. The rats were trained in the water maze to find a hidden platform. This training comprised four trials per day with a 30-s inter-trial interval between 2:00 and 8:00 pm for seven consecutive days. Each trial started with the rat placed in the middle of the outer edge of one quadrant and facing the wall of the pool, and ended when the animal climbed onto the platform. Rats that could not find the platform in 60 s were guided to the platform. The Morris water maze video tracking analysis system (Shanghai, China) was used to record the activity trajectory of the rats. The swimming paths of the rats and latencies of the rats to find the hidden platform were recorded [42]. The time the rat spent before arriving at the platform during the first trial on each day over a 7-day period was recorded as the latency time. Upon removal of the platform, which occurred during the fourth trial on each day over the 7-day period, the time the rats stayed in the platform area was recorded. The latency time and the number of times the rat crossed the platform area was used to evaluate learning ability. After 1 day of rest, the short-term memory retention test was performed. The platform was either present or absent and rats were put into the first quadrant of the maze. The latency time to reach the platform area, the number of times the rats crossed the platform area, and the total time spent in the platform quadrant, were recorded.

HE staining

After finishing Morris water maze test, rats were anaesthetized by an overdose of chloral hydrate (1 g/kg), perfused, fixed and embedded with paraffins. Brains were cut into sections (5 μm) on slides. Paraffin section was dewaxed at 65°C after further dewaxing in xylene, the slices were hydrated in gradient alcohol. Then the sections were stained in hematoxylin staining solution for 5 minutes. After rinsing, the sections were separated color in alcohol hydrochloric acid for 45 seconds. After re-rinsing again, the sections were stained with eosin for 10 seconds and rinsed again, the slices were dehydrated with gradient alcohol and transparent in xylene. Finally, the sections were mounted with neutral gum and observed under microscope (Olympus BX51, Japan). The neuron morphology and number in hippocampus and striatum were observation and counted for analysis (number per mm²) to evaluate the neuronal pathological changes.

RNA extraction and library preparation

Hippocampi tissue was separated from the rats’ brains, frozen in liquid nitrogen, and stored at -80°C until use. Total RNAs were extracted using the TRIzol reagent (Invitrogen, Singapore) according to the manufacturer’s protocol. RNA purity and quantification were evaluated using the NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). RNA integrity was assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The qualified extracted hippocampal RNA samples must meet the following conditions: the sample concentration is not less than 100 μg/μL; total RNA quantity is greater than or equal to 1μg; the purity of samples requires OD260/280 value between 1.8-2.2; the sample integrity agilent2100 score is greater than or equal to 7 (RIN≥7). Then, the libraries were constructed using TruSeq Stranded Total RNA with Ribo-Zero Gold (Cat. No. RS-122-2301; Illumina, San Diego, CA, USA) according to the manufacturer’s instructions.

RNA sequencing and differentially-expressed RNA (DE RNA) analysis

The libraries were sequenced on an Illumina HiSeq X Ten platform. After removing the adapter, ploy-N, low quality reads, and reads with length less 50 by Trimmomatic software, clean data (clean reads) were obtained. The original sequencing quantity, effective sequencing quantity, Q30 and GC content were counted and evaluated comprehensively. If the quality value is Q30, the probability of error recognition is 0.1%, that is, the accuracy rate is 99.9%. Qualified sequencing results require a quality value of at least 85%.

Sequencing reads were mapped to the human genome (GRCh38) using HISAT2 [43]. For mRNAs, the fragments per kilobase per million (FPKM) [44] of each gene was calculated using Cufflinks 2.0 [45], and the read counts of each gene were obtained by HTSeq-count [46]. Differential expression analysis was performed using the DESeq2 R package [47]. A P value < 0.05 was set as the threshold for significantly differential expression. For lncRNAs, the transcriptome from each dataset was assembled independently using the Cufflinks 2.0 program [45]. All transcriptomes were pooled and merged to generate a final transcriptome using Cuffmerge (Cufflinks 2.0). All transcripts that overlapped with known mRNAs, other non-coding RNA, and non-lncRNA were discarded. Next, the transcripts longer than 200 bp and with a number of exons > 2 were picked out, and CPC (v. 0.9-r2) [48], PLEK (v. 1.2) [49], CNCI (v. 1.0) [50], Pfam (v. 30) [51] were used to predict transcripts with coding potential. The novel predicted lncRNAs were obtained through these processes. The characteristics (including length, type, number of exons) of lncRNA were analyzed after screening. Then, the novel predicted lncRNAs and known lncRNAs (from the NCBI and Ensemble databases) were both used for expression calculation and differential screening. Then, differential expression analysis was performed using the DESeq (2012) R package. All sequencing process and analyses were performed using OE Biotech Co., Ltd. (Shanghai, China).

Validation of quantitative real-time polymerase chain reaction (qRT-PCR)

To verify the validity and accuracy of RNA sequencing, we carried out qRT-PCR assessment to assess the data consistency between RNA sequencing and qRT-PCR. The RNA amplification by qRT-PCR has been described previously [12]. Total RNA of the 50 mg hippocampi tissue was extracted using the TRIzol Reagent according to the manufacturer’s protocol (Invitrogen). First strand complementary DNA (cDNA) was synthesized from total RNA using a First-strand cDNA Synthesis Kit (Thermo Fisher Scientific). The SYBR GREEN Mix (Invitrogen, Waltham, MA USA) reaction system was used for RT-PCR along with the forward primer, the reverse primer, and cDNA. The reaction process included: 1) a preincubation step at 95ºC for 3 min, 2) an amplification step of 45 cycles of 94ºC for 30 s, 3) different annealing temperature and 72ºC for 30 s, and 4) an elongation step of 72ºC for 10 min. A melting curve was recorded to verify the absence of primer dimers. Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was the endogenous control. RNA levels were assayed using the “ΔΔ Ct method” for relative expression.

GO function and KEGG pathway annotation analyses

Gene Ontology (GO) annotated the gene function in three levels: molecular function, biological process, and cellular component. The main functions of the DE genes could be analyzed by comparisons with genes with annotated functions; therefore, the functions that may be enriched by DE genes were deduced through hypergeometric algorithms and tests to acquire the P value. Similarly, the probable pathway of DE genes was acquired by comparison with genes with an annotated pathway in the database of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and by performing hypergeometric algorithms and other tests. A P value < 0.05 was considered statistically significant. The P value determined the degree of the function and pathway enrichment of the DE genes.

LncRNA-mRNA co-expression analyses

For co-expression analysis of lncRNA and the gene, according to the expression levels of DE lncRNA and genes, the Pearson correlation analysis test was used to calculate the correlation between the two expressions. A pairing of a correlation coefficient of > 0.8 and a P value of < 0.05 was considered to indicate a co-expression relationship.

ceRNA construction of the lncRNA-miRNA-mRNA network

The miRbase database and the miRanda program (v. 3.3a) [52] were used to predict the binding between these miRNA-differential mRNA/different lncRNA sequences, using the default parameters of miRanda v. 3.3a (S ≥ 150, ΔG ≤ −30 kcal/mol and demand strict 5' seed pairing), which predicts the miRNA bound to lncRNA and mRNA, respectively, and then takes the intersection to find the miRNA bound to both. S referred to the single residue pair match scores of the matching area, and ΔG referred to free energy of the double chains binding. Then, the intersecting miRNA and the corresponding lncRNA and mRNA were used to construct a ceRNA regulatory network of lnc-miRNA-mRNA using Cytoscape software [53].

Statistical analysis

Data were expressed as means ± standard error of the mean (SEM) and analyzed using SPSS 20.0 statistical software (SPSS Inc., Chicago, IL, USA). The repeated-measures analysis of variance procedure was used to determine the statistical significance of RNAs among the three groups. The one-way analysis of variance procedure followed by Dunnett's t-test was used to determine the statistical significance of differences of the means for comparisons between the two groups. P < 0.05 was considered a statistically significant difference. Fold changes (FCs) and P value of tests were used to determine the statistical significance of the RNA sequence data. A FC ≥ 1.5 and P < 0.05 were used as thresholds for DE lncRNAs and mRNAs.