Pathological Grade-Associated Transcriptome Profiling of lncRNAs and mRNAs in Gliomas

The aim of the present study was to explore the expression profiles of lncRNAs and mRNAs in glioma patients and to elucidate any potential relationship between lncRNAs and mRNAs in glioma. High-throughput transcriptome sequencing of mRNAs and lncRNAs from six normal tissues and 16 glioma tissues (grade II, six cases; grade III, four cases; and grade IV, six cases) was performed. Series test of cluster (STC) analysis was used to screen significant trending models associated with glioma. Gene co-expression networks were constructed for the differentially expressed lncRNAs and mRNAs, and gene-ontology (GO) and pathway-enrichment analyses were further performed. Quantitative real-time PCR was performed to validate the five most differentially expressed lncRNAs and mRNAs. After filtering the raw sequencing data, we found 578 lncRNAs and 3,216 mRNAs that were significantly dysregulated in glioma (fold change ≥ 2, p < 0.05). Twenty model profiles of lncRNA and 10 model profiles of mRNA were summarized, and three patterns of lncRNAs and two patterns of mRNAs were of clinical significance. Three gene co-expression networks between mRNAs and lncRNAs were built to clarify the relationship between lncRNAs and mRNAs in glioma. GO and pathway analyses indicated that the differentially expressed lncRNAs and mRNAs were enriched in several biological processes and signaling pathways associated with tumorigenesis. Both lncRNAs and mRNAs exhibited dynamic differential expression profiles that indicated their potential roles in different degrees of glioma malignancy. A series of bioinformatics analyses indicated that most of these lncRNAs and mRNAs are involved in important biological processes and pathways associated with the pathogenesis of glioma. These results provide potential directions and valuable resources for future investigations via the comprehensive integration of these lncRNAs and mRNAs.


INTRODUCTION
Gliomas are the most common type of primary brain tumor (1) representing 75% of all malignant primary central-nervoussystem (CNS) tumors in adults (2). In the updated 2016 version of the World Health Organization (WHO) classification of CNS tumors, gliomas are divided into circumscribed gliomas (WHO grade I) and diffusely-infiltrating gliomas (whether astrocytic or oligodendroglial, WHO grades II-IV) (3). Compared to circumscribed gliomas, diffusely-infiltrating gliomas exhibit a more relentless malignant progression, a reduced efficacy to various therapeutic approaches, and a higher risk of recurrence (4). Despite efforts to promote various new therapies and advances in the research of tumor biology, the prognosis for patients with gliomas, especially diffuselyinfiltrating gliomas, is still bleak (2,5). This is mainly due to a lack of accurate biomarkers and a poor understanding of the pathogenesis of gliomas, which leads to delayed diagnoses and ineffective therapeutic outcomes. Therefore, there is an urgent need to better understand the mechanisms underlying glioma and to find potential biomarkers and accurate therapeutic targets.
Long non-coding RNAs (lncRNAs) account for a major class of non-coding RNAs (ncRNAs) and measure a length >200 nucleotides (6,7). Recent studies have demonstrated that lncRNAs may be involved in gene expression via four different processes: epigenetic regulation, translational regulation, transcriptional regulation, and post-transcriptional regulation (8)(9)(10)(11)(12). Increasing evidence has suggested that lncRNAs are vital epigenetic regulators of mRNA expression and constitute an important fraction of the human transcriptome. Furthermore, there is a growing number of studies that have reported that lncRNAs play important roles in tumor genesis, progression, and metastasis, as well as many other cellular processes (13)(14)(15)(16). Aberrant expression of lncRNAs may contribute to glioma pathogenesis, including cellular proliferation, apoptosis, and metastasis (17)(18)(19)(20). The dysregulation of lncRNAs may serve as diagnostic biomarkers of early stages of glioma and could be exploited as therapeutic targets (21)(22)(23). However, the potential pathological and biological roles of lncRNAs and mRNAs in different degrees of glioma malignancy have yet to be elucidated.
Recently, the deep sequencing of transcriptomes is being utilized with a higher sensitivity for the identification of differential expression. Advances in next-generation, deepsequencing technology have identified a number of ncRNAs. Here, we performed high-throughput transcriptome sequencing of glioma tissues and normal tissues to determine lncRNA and mRNA profiles, to investigate novel tumor-related lncRNAs and mRNAs in glioma, and to generate model profiles for future studies. We constructed a gene co-expression network for the differentially expressed lncRNAs and mRNAs in glioma tissues to further investigate the relationship between lncRNAs and mRNAs. We also conducted gene-ontology (GO) enrichment analysis and pathway-enrichment analysis for the differentially expressed lncRNAs and mRNAs. In addition, the five most differentially expressed lncRNAs and mRNAs were verified by quantitative real-time PCR (qRT-PCR).

Patients and Samples
We recruited 50 patients diagnosed with glioma or epilepsy and collected their tumor or normal tissues from March 2014 to December 2018. All the tumor patients were explicitly diagnosed with glioma by histopathological examination after surgery and were classified as grade II, grade III, and grade IV according to the CNS tumor-classification criteria (fourth edition) published by the WHO in 2016. Six normal tissues and 16 glioma tissues (grade II, six cases; grade III, four cases; and grade IV, six cases) were selected at random for high-throughput transcriptome sequencing, and qRT-PCR analysis was performed in the other samples. All the patients had no prior chemotherapy or radiotherapy and did not have any other serious diseases. The brain tissues for RNAsequencing were transferred to −80 • C storage within 60 min of resection. This experiment was approved by the Ethics committee of Nantong University, and all patients provided informed consent.

RNA Library Construction and RNA Sequencing
Total RNA was extracted from the brain tissue samples using Trizol reagent (Invitrogen, Cat no.15596-026, USA), following the manufacturer's protocol. Ribosomal RNA was removed from the total RNA samples using Ribo-Zero rRNA Removal Kits (Illumina, USA), as per the manufacturer's instructions. RNA libraries were constructed using rRNAdepleted RNAs with TruSeq Stranded Total RNA Library Prep Kit (Illumina, USA), according to the manufacturer's instructions. The libraries were controlled for quality and quantified using the BioAnalyzer 2,100 system (Agilent Technologies, USA). Then, the 10-pM libraries were denatured as single-stranded DNA molecules, captured on Illumina flow cells, amplified in situ as clusters, and finally amplified (150 cycles) and sequenced on an Illumina HiSeq Sequencer, according to the manufacturer's instructions.
The high-quality reads generated were aligned to the human reference genome (UCSC hg19) with hisat2 software. Then, guided by the Ensembl gene-annotation file, cuffdiff software (part of cufflinks) was used to reveal the expression profile of the lncRNAs and mRNAs in terms of Fragments Per Kilobase of transcript per Million mapped reads (FPKM) values, from which the fold change between groups and the corresponding p-values were calculated. Subsequently, differentially expressed lncRNAs and mRNAs were identified and lncRNA target genes were predicted by their locations to nearby genes.

Series Test of Clusters (STC)
STC analysis was used to screen the significant trending models associated with glioma and their corresponding differentially expressed mRNAs and lncRNAs. Fisher's exact test was used to identify significant profiles, and p ≤ 0.05 was used as a threshold of significance. Transcript_ID, the transcript ID; Gene_ID, the Ensembl gene identifier; Log 2 (fold change), fold change between the two groups; p_value, p-value between the two groups; regulation, "up" indicates upregulation, and "down" indicates downregulation.
FIGURE 1 | Differentially expressed lncRNAs and mRNAs in glioma tissues and normal tissues. The hierarchical clustering and heat map show differential lncRNA (A) and mRNA (B) expression profiles of all targets among samples: red represents high relative expression, and green represents low relative expression. The variation of differential lncRNAs and mRNAs between the glioma and normal groups is shown in a scatter plot (C,D). The red dots and green dots denote a fold change >2.

Gene Co-expression Analyses
To explore the interactions between the DEGs and differentially expressed lncRNAs, gene co-expression networks were build based on their co-expression patterns. The lncRNAs with a related coefficient of R ≥ 0.95 or R ≤ −0.95 were screened for functional analysis.

Gene Ontology (GO) Analysis
All differentially expressed genes (DEGs) were mapped to GO terms in the GO database (http://www.geneontology. org/). A hypergeometric test was applied to find significantly enriched GO terms in the input list of DEGs, based on "GO::TermFinder" (http://smd.stanford.edu/help/GO-TermFinder/GO_TermFinder_help.shtml).
A Bonferroni correction was applied to adjust the p-value. The false discovery rate (FDR)adjusted p ≤ 0.05 was used as a threshold and GO terms fulfilling this condition were defined as significantly enriched.

Pathway Analysis
Pathway analysis was used to identify pathways involving the DEGs, according to the Kyoto Encyclopedia of Genes and Genomics (KEGG). Pathways with FDR-adjusted p ≤ 0.05 were defined as significantly enriched. Cytoscape was used to generate graphical representations of the pathways.

Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)
Reverse transcription was performed with the High-Capacity cDNA Reverse Transcription Kits (Applied Biosystems, Foster City, USA), according to the manufacturer's instructions. qRT-PCR was performed on an ABI 7500 thermocycler (Applied Biosystems, Foster City, CA, USA) by using SYBR Green Real-Time PCR Master Mix (Toyobo, Japan). GAPDH was used for normalization. All qPCR reactions were performed in biological triplicates. Primer sequences are listed in Supplemental Table 1.

Statistical Analysis
Student's t-tests and one-way analyses of variance (ANOVA) with Bonferroni corrections for multiple comparisons were performed to determine significant differences between different groups. A false discovery rate (FDR)-adjusted p ≤ 0.05 was regarded as statistically significant. All statistical details are specified in the figure legends.

Differentially Expressed lncRNAs and mRNAs in Glioma Tissues Compared With Normal Tissues
The lncRNA and mRNA expression levels were compared in glioma tissues and normal tissues. We found 578 lncRNAs and 3,216 mRNAs that were significantly dysregulated in glioma tissues. Among these, 509 lncRNAs and 2,282 mRNAs were upregulated and 69 lncRNAs and 934 mRNAs were downregulated (fold change ≥ 2.0, p < 0.05). We then used hierarchical-clustering analysis to reveal betweengroup comparisons of lncRNA and mRNA expression levels (Figures 1A,B). In addition, the variation of differential lncRNAs and mRNAs between the glioma and normal groups is shown in a scatter plot (Figures 1C,D). The

Model Profile Analysis of lncRNAs and mRNAs in Glioma Tissues and Normal Tissues
To narrow down the number of highly significant differentially expressed lncRNAs and mRNAs, we further analyzed their specific expression patterns. Twenty model profiles of lncRNAs and 10 model profiles of mRNAs were summarized. Among the 20 patterns, we identified nine patterns of lncRNAs that exhibited significant p-values (Figure 2A; p-values in red boxes).
Among the nine significant patterns, the expression of lncRNAs in profile No. 1, 2, and 22 were of clinical significance (Figures 2B-D). The lncRNA model profile No.1 and No. 2 contained 78 lncRNAs and 58 lncRNAs, respectively, the expressions of which were decreased consistently in glioma tissues (grades II-IV). Additionally, lncRNA profile No. 22 was constructed with 457 lncRNAs, which exhibited consistently up-regulated expression in glioma tissues (grades II-IV).

Establishment of the Gene Co-expression Network for lncRNAs and mRNAs in Glioma Tissues
To clarify the relationship between lncRNAs and mRNAs in glioma, we performed correlation analyses for lncRNAs and mRNAs in terms of their expression values in glioma tissues. Additionally, a gene co-expression network between mRNAs and lncRNAs was constructed (Figures 3-5). In profile No. 22, there were 83 square nodes and 287 circular nodes that represented lncRNAs and mRNAs, respectively. Moreover, the edges showed the interaction between the lncRNAs and mRNAs. These results indicated that lncRNAs may play vital roles in the pathogenesis of glioma.

GO-and Pathway-Enrichment Analyses of Differentially Expressed lncRNAs
To further identify the functional roles of these differential lncRNAs found to be dysregulated in the tumor group, we conducted GO-enrichment and pathway-enrichment analyses. The GO analysis returned terms associated with three categories: molecular function (MF), cellular component (CC), and biological process (BP). The number of lncRNAs in profile No. 22 found associated with each GO term was counted and are shown in a pie chart (Figure 6).
The 10 most enriched GO terms (in descending order of enrichment score) within the three categories are shown in Figures 6D-F. Furthermore, KEGG pathway analysis demonstrated that the differentially expressed lncRNAs were significantly enriched in various important pathways. The dot plot in Figure 7A shows the eight highest enrichment scores (lowest log 10 pvalues) of the significant pathways (Figure 7A), and Figure 7B shows the regulatory roles of the lncRNAs involved in cancer pathways.

GO-and Pathway-Enrichment Analyses of Differentially Expressed mRNAs
Since the functions of the differentially expressed mRNAs are different from those of the lncRNAs, we conducted independent GO-enrichment and pathway-enrichment analyses for the mRNAs. Figure 8 shows the number of mRNAs associated with each GO term (Figures 8A-C), and the 10 most enriched GO terms (in descending order of enrichment score) in each of the three categories are shown in Figures 8D-F. As observed with the lncRNAs, KEGG pathway analysis demonstrated that the differentially expressed mRNAs were significantly enriched in various important pathways. The dot plot shows the eight most significantly enriched pathways (enrichment score = -log 10 p-value) (Figure 9A), and the regulatory roles of the mRNAs involved in systemic lupus erythematosus and the staphylococcus aureus infections are shown in Figure 9B.

DISCUSSION
As new research techniques have developed, a growing number of high-throughput platforms have been used to analyze gene expression in glioma. However, previous studies mostly focused on mRNA, DNA, or protein levels in glioblastoma or highgrade glioma using proteomics or microarrays (24)(25)(26). In the present study, we assessed the significantly differentially expressed lncRNAs and mRNAs between normal tissues and glioma tissues using high-throughput sequencing. To our knowledge, this study is the first to identify the pathological grade-associated transcriptome profiles of lncRNAs and mRNAs  in glioma. We identified 578 lncRNAs and 3,216 mRNAs that were significantly dysregulated in glioma tissues. Among these, 509 lncRNAs and 2,282 mRNAs were upregulated, whereas 69 lncRNAs and 934 mRNAs were downregulated (fold change ≥ 2, p < 0.05). This result indicated that these lncRNAs and mRNAs may be involved in glioma initiation and/or progression. However, further research is required to elucidate the detailed mechanisms of the involvement of these lncRNAs and mRNAs in glioma.
The lncRNAs and mRNAs that dynamically changed with differing degrees of glioma malignancy may play crucial biological roles in the disease process (27, 28). Through model-profile analysis, the dynamic expression profiles of lncRNAs and mRNAs were obtained, and the nine significant dynamic expression profiles of lncRNAs were then screened. The three profiles of clinical significance were profiles No. 1,No. 2,and No. 22. Each of these profiles contained a large number of lncRNAs that were consistently downor up-regulated in the different grades of glioma tissues. In the most significant profile (No. 22), seven lncRNAs were identified to be involved in viral carcinogenesis. The lncRNAs with significant dynamic-expression changes may be more correlated with the malignancy of glioma and, hence, may play vital roles in the regulation of glioma initiation and progression.
To further discern the key lncRNAs associated with glioma, we integrated lncRNA and mRNA co-expression networks in profiles No. 1,No. 2,and No. 22. In total, 83 lncRNAs and 287 mRNAs were identified to play vital regulatory roles in model profile No. 22. Many studies have reported that a large number of lncRNAs, such as ZEB1-AS1 and PCNA-AS1 (29,30) played important roles in the accurate and complicated co-expression networks. We screened a number of lncRNAs and mRNAs that included key genes that are closely associated with the pathogenesis of glioma. For example, SOX21 has been reported to be closely related to the tumorigenesis of glioblastoma, hepatocellular carcinoma, and colorectal cancer (31)(32)(33). This suggests that these co-expressed lncRNAs and mRNAs may participate in cancer-related pathways. It was also found that most lncRNAs can be co-expressed with various mRNAs, indicating that each lncRNA may regulate many mRNAs.
A growing number of lncRNAs have been identified in tumors (34,35); however, the various functions of lncRNAs still remain poorly characterized. To infer the functional roles of the lncRNAs in glioma, GO and pathway analyses were performed for these differentially expressed lncRNAs. The GO analysis indicated that these lncRNAs are most significantly enriched in domains of "biological regulation, " "metabolic process, " and "regulation of biological process, " which are strongly associated with glioma. In addition, the pathway analysis revealed that "pathway in cancer, " "calcium signaling pathway, " "cAMP signaling pathway, " and "Ras signaling pathway" are related with the pathogenic process of glioma.
In addition, there were also some limitations of our present study. Regardless of the technology performed to detect expression levels and the sample sizes that are used, the actual gene expression levels vary among individuals because expression is a random process. Consequently, the analysis data may not be powerful enough to reflect these actual expression levels across individuals. However, biological variability will decrease with the increase of the scale of samples. Therefore, we will perform further studies with more samples in the future.

CONCLUSIONS
In summary, we characterized the expression profiles of lncRNAs and mRNAs in normal tissues and glioma tissues (WHO grades II-IV). Then, we analyzed the dynamic differentially expressed profiles of lncRNAs and mRNAs, which indicated their potential vital roles in gliomas with different degrees of malignancy. A series of bioinformatics analyses indicated that most of these lncRNAs and mRNAs are involved in important biological processes and pathways associated with the pathogenesis of glioma. These results provide potential directions and valuable resources for future studies via the comprehensive integration of these lncRNAs and mRNAs.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in GenBank (SRA accession: PRJNA604108).

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethics committee of Nantong University. The patients/participants provided their written informed consent to participate in this study.