A focal cerebral ischemia model induced by MCAO, prepared as previously described (18), was prepared using 20 7-week-old male Sprague-Dawley rats of a specific pathogen-free grade (weighing 200±20 g), purchased from the experimental animal center of Anhui Medical University (Anhui, China). The study protocol was ethically approved by the Committee on the Ethics of Animal Experiments of Anhui University of Chinese Medicine (approval no. 2012AH-036-03). In brief, the animals were fasted overnight but allowed ad libitum access to water. They were then anesthetized with chloral hydrate (350 mg/kg, intraperitoneal injection). A 4-0 silicon-coated monofilament nylon suture with a round tip was inserted through an arteriectomy in the common carotid artery just below the carotid bifurcation and then advanced into the internal carotid artery ~18 mm distal to the carotid bifurcation until a mild resistance was felt. Following 2 h of MCAO, the filament was removed to allow reperfusion. As a control, control-operated rats underwent identical surgery but did not have the suture inserted. Four days subsequent to MCAO, the left hemispheres were collected and immediately frozen in liquid nitrogen.

RNA-seq was performed by Ao-Ji Bio-Tech (Shanghai, China). Briefly, total RNA was extracted using an RNeasy Mini kit (Qiagen GmbH, Hilden, Germany), according to the manufacturer's protocol. The RNA quality control was performed using Nanodrop 2000 and Agilent 2100, and mainly depended on the concentration, purity and integrity of the RNA. Ribosomal RNA was removed from total RNA using Ribo-Zero rRNA removal beads (Illumina, Inc., San Diego, CA, USA). Libraries were constructed according to the standard TruSeq protocol (19). Purified cDNA libraries were prepared for cluster generation and sequencing on an Illumina HiSeq 2500 (Illumina, Inc.) according to the manufacturer's protocol. Subsequently, data analyses were performed in silico.

Quality control of the RNA-Seq reads was conducted using FastQC (v0.11.3) (The Babraham Institute, Cambridge, UK). Reads were trimmed using the software seqtk (github.com/lh3/seqtk) for known Illumina TruSeq adapter sequences, poor reads and ribosome RNA reads. Trimmed reads were aligned to the rat genome (Rn6) using Hisat2 (version 2.0.4) (20). Transcripts were assembled using Stringtie (v1.3.0) (20,21). Transcripts constructed from Stringtie were compiled together by gffcompare (v0.9.8) (20,21). Transcripts detected in at least five samples (half of the total number) were considered to be bona fide transcripts. Transcripts, with the exception of those with just one exon and shorter than 200 base pairs, were further analyzed for the identification of lncRNAs. Transcripts with class codes ‘i,’ ‘u,’ and ‘x,’ were considered to be potential novel long transcripts. Pfam (22), Coding Potential Calculator (CPC) (23) and Coding-Non-Coding Index (CNCI) (24) were used to estimate the coding potential of each novel transcript. Transcripts with a Pfam score <0, CNCI <0 and CPC non-significant were considered to lack coding potential. Transcripts were compared with annotation databases, including NONCODE (v4) (http://www.noncode.org) and Ensembl (25). The matched transcripts were considered to be known lncRNAs, and others were considered to be novel lncRNAs. All lncRNAs were quantified using Stringtie. According to the positional association between lncRNA and mRNA in the genome, lncRNA may be classified into six types: Bidirectional, exonic_antisense, exonic_sense, intergenic, intronic_antisense and intronic_sense (26).

Initially, the DELs and differentially expressed genes (DEGs) were analyzed using EdgeR (27). For DEGs, log2| [fold change (FC)] |>1 and P<0.05 were used as the cutoff values. Meanwhile, log2| (FC) |>1 and false discovery rate (FDR) <0.05 were used as the threshold for DELs. Hierarchical clustering of DELs was performed based on mean signals using a Euclidean distance function. In addition, a volcano plot was generated. The Pearson's correlation coefficient (PCC) between lncRNAs and mRNAs was calculated (cutoff value, PCC>0.9, P<0.05) and the lncRNA-mRNA regulatory network was structured using Cytoscape 2.8.3 (28).

cis- or trans-acting algorithms were used to predict the potential targets of lncRNAs. The first algorithm predicted potential target genes of cis-acting lncRNAs that were physically located within 10 kb upstream or 20 kb downstream of lncRNAs using liftOver genome browser (genome.ucsc.edu/cgi-bin/hgLiftOver). The second algorithm predicted potential target genes of trans-acting lncRNAs based on the lncRNA-mRNA complementary sequences, and predicted lncRNA-mRNA duplex energy. First, BLASTN (29) was performed to detect potential target mRNA sequences with >95% identity and E value <1×10⁻⁵ (https://blast.ncbi.nlm.nih.gov/Blast.cgi). Then, RNAplex (30) was used to calculate the complementary energy between lncRNAs and their potential trans-regulated target genes with RNAplex-10⁻³⁰. Gene Ontology (GO) (31) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (32) enrichment analyses of the identified potential target genes were performed using the Database for Annotation, Visualization and Integrated Discovery (33); and P<0.05 was considered to indicate a statistically significant difference.

Total RNA was extracted from left hemisphere samples using TRIzol^® reagent (Invitrogen; Thermo Fisher Scientific, Inc., Waltham, MA, USA) and reverse-transcribed using a Thermo Fisher Scientific RevertAid First Strand cDNA Synthesis kit (cat. no. K1622; Thermo Fisher Scientific, Inc.) according to the manufacturer's protocol at 42°C for 60 min. To further confirm the expression data from RNA-seq, a cutoff value (FC>2, P<0.05) was randomly selected for qPCR verification. The expression levels of six randomly DELs (NONRATT027551.2, MSTRG.1836.1, MSTRG.4344.10, MSTRG.7720.11, NONRATT005132.2 and MSTRG.20633.3) were assayed using a SYBRGreen flurophore (Applied Biosystems; Thermo Fisher Scientific, Inc.) using the PikoReal real-time PCR system (Thermo Fisher Scientific, Inc.) under the following conditions: Initial denaturation at 95°C for 30 sec, followed by 40 cycles at 95°C for 30 sec and 60°C for 30 sec, and a final extension step at 4°C for 20 min. FC was determined using the 2^−ΔΔCq method (34). GAPDH mRNA was used as an internal control. The primers used are listed in Table I.

The comparisons between the MCAO group and the control group were determined using a Student's t-test for the RT-qPCR results by SPSS 22.0 statistical software (IBM Corp., Armonk, NY, USA). P<0.05 was considered to indicate a statistically significant difference. The PCC between lncRNAs and mRNAs was calculated using the Hmisc package in R based on the expression determined using RNA-seq (PCC>0.9, P<0.05). The correlation analysis between the RT-qPCR results and RNA-seq results was calculated in Excel 2013 (Microsoft Corporation, Redmond, WA, USA) with the function of CORR.

The present study characterized the lncRNA landscape and expression by performing deep RNA-seq experiments on three control and three MCAO tissue samples. Subsequent to the seqtk quality assessment of sequencing, >33 million total original reads for each sample were obtained, and the proportion of bases with quality values >20 was >94%. These results indicated that the quality of the sequencing results was acceptable (Table II). Subsequent to filtering out the adaptor sequence and low quality reads, the percentage of clean reads within the raw reads accounted for 94% of the total sequences in two groups. Hisat2 software was used to map the obtained clean reads to the Rattus norvegicus reference genome. As presented in Table II, ~97% of the trimmed reads were mapped onto the reference genome. In total, 24,304 lncRNAs were screened from six samples, and there were 23,255 shared lncRNAs detected in the MCAO and control groups (Fig. 1A). The majority of the identified lncRNAs were transcribed from protein-coding exons; others were from introns and intergenic regions (Fig. 1B). In addition, the present study analyzed the distribution of the identified lncRNAs on the rat chromosomes; 24,304 lncRNA transcripts were identified in all chromosomes, and chromosome 1 included the most lncRNAs (Fig. 1C).

EdgeR was used to filter DEGs and DELs and differentiate their expression between the control and MCAO groups. A total of 1,007 DEGs (|FC|>2, P<0.05) were identified, including 785 upregulated genes and 222 downregulated genes. Similarly, as presented in Fig. 2, 134 DELs (|FC|>2, FDR<0.05) were identified in the MCAO group (Fig. 2A and B), including 77 upregulated and 57 downregulated DELs (Fig. 2C and D). In the present study, it was revealed that the FC values of certain DELs were equal to positive infinity and negative infinity, meaning that these lncRNAs are switched-on or off with MCAO. Essentially, positive or negative infinity indicates zero expression of the lncRNA in normal or MCAO groups. It was speculated that this may be associated with the abundance of lncRNAs and the sensitivity to RNA-seq. The top five upregulated DELs were NONRATT027551.2, MSTRG.1836.1, MSTRG.4344.10, NONRATT028102.2 and MSTRG.31500.2; the top five downregulated DELs were MSTRG.7720.11, NONRATT005132.2, MSTRG.20633.3, NONRATT020232.2 and MSTRG.1836.3.

The cutoff correlation r-values (|PCC|>0.9) and P-values (P<0.05) were selected to structure a lncRNA-mRNA co-expression network between DEGs and DELs. As Fig. 3 presents, 46 DEGs, 104 DELs and 664 edges were filtered out using Cytoscape to construct the co-expression network. The co-expression-associated top 30 GO terms and pathway terms enrichment analyses presented in Figs. 4 and 5 suggest that these DELs were associated with the cellular response to inorganic substances, alkaloids, estradiols, reactive oxygen species, metal ions and oxidative stress. In particular, the heme oxygenase 1 (HO-1) gene participates in many of these functions. A multitude of pathways were implicated, including metabolic pathways, chemokine signaling pathways, malaria, Parkinson's disease and the cell cycle. Notably, the BUB1 mitotic checkpoint serine/threonine kinase B (BUB1B) and C-C motif chemokine ligand 2 (CCL2) genes were associated with the cell cycle.

A total of 91 cis-regulatory genes of 94 DELs, including 55 upregulated lncRNAs in the MCAO group were identified; 14 of the 91 cis-regulatory genes exhibited differential expression. A total of 13 of the DEL/cis-regulatory gene pairs had positive correlations as follows: NONRATT021925.2 (Rho GDP dissociation inhibitor β), NONRATT004791.2 (G protein subunit γ transducing 2), NONRATT015286.2 (periostin), MSTRG.30235.10 (LRR binding FLII interacting protein 1), NONRATT015403.2 (IQ motif containing GTPase activating protein 3), NONRATT008267.2 (kinesin family member 14), NONRATT009960.2 (non-SMC condensin I complex subunit G), NONRATT011312.2 (retinol binding protein 3), NONRATT005985.2 (DNA topoisomerase IIα), NONRATT016680.2 (ENSRNOG00000053081), NONRATT013960.2 (galanin receptor 1), NONRATT016022.2 (CART prepropeptide) and NONRATT024616.2 (δ like non-canonical Notch ligand 1) (Table III). Additionally, 90 trans-regulatory genes of lncRNAs were filtered by BLASTN and RNAplex, with a negative correlation identified between ENSRNOT00000092040 and Ccl9 (Table IV).

From the data in Fig. 6, NONRATT027551.2, MSTRG.1836.1 and MSTRG.4344.10 were identified to be significantly upregulated in the MCAO group compared with the control (P<0.05), consistent with the RNA-seq data, while MSTRG.7720.11, NONRATT005132.2 and MSTRG.20633.3 were significantly downregulated in the MCAO group compared with the control (P<0.01), also consistent with the RNA-seq data. These results, revealing that the RNA-seq results were consistent with the RT-qPCR results, verified that the RNA-seq results were reliable (Fig. 6).