Comprehensive analysis of differentially expressed profiles of lncRNAs and mRNAs reveals ceRNA networks in the transformation of diffuse large B-cell lymphoma.

Diffuse large B-cell lymphoma (DLBCL) is one of the malignancies with a high mortality rate. The molecular mechanisms involved in transformation of DLBCL remain unclear. Therefore, it is critically important to investigate the biological mechanisms of DLBCL. Accumulating evidence indicates that long non-coding RNAs (lncRNAs) serve key functions in tumorigenesis, cancer progression and metastasis. Compared with follicular lymphoma (FL), a total of 123 upregulated lncRNAs and 192 downregulated lncRNAs in DLBCL were identified. Subsequently, a specific DLBCL-associated competing endogenous RNA (ceRNA) network and a specific FL-associated ceRNA network was constructed. Gene Oncology and Kyoto Encyclopedia of Genes and Genomes pathway analysis revealed that differentially expressed lncRNAs served key functions in regulating signal transduction, transcription, cell adhesion, development and protein amino acid phosphorylation. Furthermore, the molecular functions of PRKCQ antisense RNA 1, HLA complex P5, OIP5 antisense RNA 1, growth arrest specific 5 and taurine upregulated 1 were investigated, and it was revealed that these lncRNAs served important functions in regulating a series of biological processes, including anti-apoptosis, cell cycle, DNA repair, response to oxidative stress and transcription. The present study may provide a potential novel therapeutic and prognostic target for the treatment of DLBCL.

Introduction

Follicular lymphoma (FL) accounts for ~30% of all non-Hodgkin's lymphomas (1). FL is usually indolent, and patients have long survival rates (2). However, in 25–60% of all patients with FL, FL undergoes transformation into diffuse large B-cell lymphoma (DLBCL), which results in rapid progression, treatment resistance and mortality (2). DLBCL is a malignancy with a high mortality rate due to the lack of biomarkers for early diagnosis and efficient therapeutic strategies (3). Therefore, it is critically important to identify biomarkers for DLBCL and to investigate the biological mechanisms of DLBCL.

Accumulating evidence reveals that lncRNAs serve key functions in tumorigenesis, cancer progression and metastasis (4,5). Long non-coding RNAs (lncRNAs), a major class of non-coding RNAs, are RNA polymerase II transcripts that are >200 bp and do not encode proteins (6). Multiple reports have revealed that lncRNAs may regulate the expression of protein-coding genes through transcriptional, post-transcriptional, post-translational and/or epigenetic regulation (7,8). Previously, a number of studies had revealed that lncRNA expression may be deregulated in various types of human cancer (9,10). For example, prostate cancer associated 3 (non-protein coding) was significantly upregulated in prostate cancer, compared with health tissues (11). Additionally, it was indicated that H19, imprinted maternally expressed transcript (non-protein coding) was overexpressed in hepatocellular carcinoma and that this overexpression was disease-associated (12). According to the competing endogenous RNA (ceRNA) hypothesis (13), ceRNAs may compete for the same micro RNA (miRNA) response elements to regulate each other (14). Previously, studies revealed that the ceRNA network may serve prognostic or diagnostic functions in cancer. For example, Zhou et al (15,16) identified dysregulated lncRNA-associated ceRNA networks as biomarkers for pancreatic and ovarian cancer. Previously, a number of studies indicated that altered expression of certain lncRNAs may be an important mechanism of DLBCL progression. A number of lncRNAs, including HOX transcript antisense RNA (HOTAIR) (17), tumor protein p53 pathway corepressor 1 (lincRNA-p21) (18), paternally expressed 10 (PEG10) (19), MEF2C antisense RNA 2, SACS antisense RNA 1, RP11-25K19.1, MME antisense RNA 1, RP11-360F5.1 and CSMD2 antisense RNA 1 (20) which were significantly associated with the survival outcomes of DLBCL. Peng et al (21,22) reported that hepatocellular carcinoma upregulated long non-coding RNA and leukemia-associated non-coding IGF1R activator RNA 1 were associated with cell proliferation in DLBCL. Zhou et al (23) identified a 17-lncRNA signature for subtype classification and prognosis prediction by analyzing differentially expressed lncRNAs between germinal center B-cell-like and activated B-cell-like subtypes. However, the molecular mechanisms and functions underlying the involvement of lncRNAs in the transformation of DLBCL remain largely unknown.

In the present study, the aim was to identify differentially expressed lncRNAs and mRNAs involved in the transformation of DLBCL by analyzing a cohort of previously published datasets from the Gene Expression Omnibus (GEO). In an attempt to provide novel information on the molecular mechanisms and functions of lncRNAs, a bioinformatics analysis was conducted to identify the lncRNA-miRNA-mRNA regulatory axis in DLBCL. Subsequently, Gene Ontology and Kyoto Encyclopedia of Genes and Genomes analysis was performed in order to investigate the potential functions of dysregulated lncRNAs.

Materials and methods

Microarray data and data preprocessing

Microarray data was downloaded from a previous study by Brodtkorb et al (24), which was referenced in the GEO database (accession no. GSE53820; www.ncbi.nlm.nih.gov/geo/). In this dataset, preprocessed usingt the limma package in R (version, 3.34.2; www.r-project.org/), expression profiles were obtained from a total of 81 biopsies, which were taken from 41 patients diagnosed with FL using the Affymetrix HG U133 Plus 2.0 Gene Chip (Affymetrix; Thermo Fisher Scientific, Inc., Waltham, MA, USA). Of these, 49 biopsies (43 with a histological diagnosis of FL and 6 with DLBCL) were sourced from 24 patients with a subsequent transformation to a higher-grade lymphoma (namely, DLBCL) and 32 biopsies were sourced from 17 patients without any sign of transformation. The cut-off values used for selecting differentially expressed mRNAs were fold change ≥2 and P<0.05.

lncRNA classification pipeline

In order to evaluate the expression of lncRNAs in microarray data, a pipeline, previously described by Zhang et al (25), was employed to identify the probe sets uniquely mapped to lncRNAs from the Affymetrix array by using the following criteria: i) For the probe sets with Refseq IDs, those labeled as ‘NR’ (where NR indicates non-coding RNA in the Refseq database) were retained; ii) for the probe sets with Ensembl gene IDs, those annotated with ‘lncRNA’, ‘processed transcripts’, ‘non-coding’ or ‘misc_RNA’ in Ensembl annotations were retained; iii) the probe sets obtained were refined by filtering pseudogenes, ribosomal RNAs, microRNAs, transfer RNA (tRNA)s, small nuclear RNAs and small nucleolar RNAs. A total of 2,448 annotated lncRNA transcripts with corresponding Affymetrix probe IDs were obtained. The cut-off values used for selecting differentially expressed lncRNAs were fold change ≥2 and P<0.05.

Functional group analysis

GO analysis and KEGG analysis were employed to determine the biological functions of the identified differentially expressed mRNAs, based on the freely available online MAS 3.0 system from CapitalBio Corporation (http://bioinfo.capitalbio.com/mas3/; Beijing China). The P-value (hypergeometric P-value) denotes the significance of the pathway associated with the conditions. P<0.05 was considered to indicate a statistically significant difference.

Construction of the lncRNA-miRNA-mRNA network

To predict the functions of the differentially expressed lncRNAs, co-expression networks of differentially expressed lncRNAs were constructed for further bioinformatics analysis, as previously described by Guttman and Rinn (26) and Shen et al (27). The StarBase dataset (27) was used to identify potential dysregulated lncRNA-miRNA pairs. StarBase and TargetScan (28) databases were also used to identify miRNA-mRNA pairs. Finally, a co-expression network based on association analysis between the differentially expressed lncRNAs and mRNAs was constructed. The lncRNA-mRNA interaction was integrated into the co-expression networks according to positive regulation and only gene pairs with |R| >0.5 were selected.

Hierarchical clustering analysis

To generate an overview of lncRNA and mRNA expression profiles between FL and DLBCL, hierarchical clustering analysis was performed based on the expression values. Cytoscape 3.0 was applied to plot the lncRNA-miRNA-mRNA ceRNA networks involved in the transformation of DLBCL.

Statistical analysis

All numerical data (log-transformed) are presented as the mean ± standard deviation of at least 3 determinations. Statistical comparisons between groups of normalized data were performed using a Student's t-test (unpaired) or Mann-Whitney U-test according to the test conditions. P<0.05 was considered to indicate a statistically significant difference with a 95% confidence level. All of the above statistical analyses are analyzed with R software, version 3.2.4 (https://www.r-project.org/).

Results

Systematic comparison of differentially expressed mRNAs and lncRNAs between FL and DLBCL

In order to compare differentially expressed mRNAs and lncRNAs between FL and DLBCL, a publicly available gene expression database (accession no. GSE53820) was utilized. This database includes 75 FL samples and 6 DLBCL samples. Differentially expressed mRNAs in the GSE53820 database were analyzed, and it was identified that 1,884 genes were upregulated and 814 genes were downregulated in DLBCL compared with FL.

Based on the NetAffx annotation of the probe sets and the Refseq and Ensemble annotations of lncRNAs, a total of 2,448 lncRNA transcripts (corresponding to 1,970 lncRNA genes) were identified in the GSE58320 database. lncRNA expression patterns between FL and DLBCL were compared, and a total of 123 lncRNAs were significantly upregulated and 192 lncRNAs were significantly downregulated (P<0.05) in DLBCL compared with FL.

GO and KEGG analysis of differentially expressed mRNAs

To identify the potential functions of differentially expressed mRNAs, GO and KEGG analysis were performed using MAS 3.0 software. GO analysis indicated that the upregulated genes were primarily involved in the regulation of cell cycle, cell division, mitosis, DNA-dependent regulation of transcription and DNA replication, which are mainly associated with cell proliferation (Fig. 1A). KEGG pathway analysis revealed that upregulated genes were primarily enriched in pathways associated with cell cycle, pyrimidine metabolism, ubiquitin mediated proteolysis and Wnt signaling pathway (Fig. 1B).

Figure 1.

GO and KEGG pathway analysis of differentially expressed mRNAs between follicular lymphoma and diffuse large B-cell lymphoma. (A) GO and (B) KEGG pathway analysis of the upregulated mRNAs. (C) GO and (D) KEGG pathway analysis of the downregulated mRNAs. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MAPK, mitogen-activated protein kinase; MHC, major histocompatibility complex; CAMs, cell adhesion molecules.

Meanwhile, downregulated genes were mainly enriched in categories associated with transcription, immune response, interspecies interaction between organisms and signal transduction (Fig. 1C). These results suggest that these pathways may participate in regulating the transformation of FL. Downregulated genes were mainly associated with T cell receptor signaling pathway, cell adhesion molecules, cytokine-cytokine receptor interaction, MAPK signaling pathway and natural killer cell mediated cytotoxicity (Fig. 1D).

GO and KEGG analysis of differentially expressed lncRNAs

Co-expression networks were constructed to identify the association between differentially expressed mRNAs and lncRNAs using the GSE53820 database. The cut-off values used for selecting differentially expressed lncRNAs were a fold change ≥2 and P<0.05. GO and KEGG analyses were performed for each lncRNA using the set of co-expressed mRNAs.

In the present study, the top 500 differentially expressed lncRNAs and mRNAs were classified according to GO terms (Fig. 2A and C). GO analysis revealed that the upregulated lncRNAs were enriched in transcription, cell cycle, cell division, mitosis, and protein amino acid phosphorylation (Fig. 2A), while the downregulated lncRNAs were enriched in transcription, transport, cell cycle, interspecies interaction between organisms, and oxidation reduction (Fig. 2C).

Figure 2.

GO and KEGG pathway analysis of differentially expressed lncRNAs between follicular lymphoma and diffuse large B-cell lymphoma. (A) GO and (B) KEGG pathway analysis of the upregulated lncRNAs. (C) GO and (D) KEGG pathway analysis of the downregulated lncRNAs. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; lncRNA, long non-coding RNA; TCA, citrate acid cycle; tRNA, transfer RNA.

According to KEGG pathway analysis, upregulated lncRNAs were primarily enriched in pathways associated with cell cycle, T-cell receptor signaling pathway, pyrimidine metabolism and purine metabolism (Fig. 2B). Downregulated lncRNAs were enriched in pathways associated with cell cycle, purine metabolism, aminoacyl-tRNA biosynthesis and valine, leucine and isoleucine degradation (Fig. 2D).

Construction of the ceRNA networks

In order to investigate the molecular mechanisms of lncRNAs, the lncRNA-miRNA-mRNA axis was predicted in the present study. Firstly, the interactions between differentially expressed lncRNAs and their theoretical target miRNAs was predicted using the StarBase database (27). Then, TargetScan (28) and the StarBase database were employed to identify mRNAs targets that are suppressed by miRNAs. Finally, a co-expression network based on the association analysis between the differentially expressed lncRNAs and mRNAs was constructed. The lncRNA-mRNA interaction was integrated into the co-expression networks according to positive regulation and only gene pairs with |R| >0.5 were selected. lncRNA-miRNA-mRNA ceRNA networks involved in the transformation of DLBCL were constructed using Cytoscape 3.0 (http://www.cytoscape.org/).

The results of the present study revealed a specific DLBCL-associated and a specific FL-associated ceRNA network. As presented in Fig. 3, 14 lncRNAs, (including OIP5-AS1, SNHG16, HOXA11-AS and NUTM2A-AS1), 198 miRNAs, and >1,200 mRNAs were involved in the specific DLBCL-associated ceRNA network. It was revealed that the FL-associated ceRNA network included 8 lncRNAs, (including HCP5, COX10-AS1, PRKCQ-AS1 and LEMD1-AS1), 71 miRNAs, and >200 mRNAs (Fig. 4). The networks were constructed using Cytoscape 3.0.

Figure 3.

The network of competing endogenous RNA involved in the transformation of diffuse large B-cell lymphoma. Red triangle nodes represent lncRNAs, white square nodes represent miRNAs, and blue square nodes represent mRNAs.

Figure 4.

Construction of the competing endogenous RNA networks involved in the transformation of follicular lymphoma. Triangle nodes represent lncRNAs, circular nodes represent miRNAs, and square nodes represent mRNAs. miR, microRNA.

Investigating the molecular functions of PRKCQ-AS1, HCP5, OIP5-AS1, growth arrest specific 5 (GAS5) and taurine upregulated 1 (TUG1)

According to the ceRNA networks, it was revealed that PRKCQ-AS1, HCP5, OIP5-AS1, GAS5 and TUG1 functioned as key regulators (Figs. 5 and 6). However, the molecular functions of PRKCQ-AS1, HCP5, OIP5-AS1, GAS5 and TUG1 in the transformation of DLBCL remain unknown. By analyzing co-expressed mRNAs, it was revealed that TUG1, OIP5-AS1 and GAS5 were associated with anti-apoptosis, cell cycle, DNA repair, mitosis, transcription, mitosis, G2/M transition of mitotic cell cycle and protein amino acid phosphorylation functions (Fig. 5). PRKCQ-AS1 was associated with the response to oxidative stress, regulation of smooth muscle cell proliferation and acute-phase response functions (Fig. 6). HCP5 was associated with transcription, cell adhesion, lipid metabolism and immune response functions (Fig. 6).

Figure 5.

Analysis of the associated molecular functions of OIP5-AS1, GAS5 and TUG1 in the transformation of diffuse large B-cell lymphoma. Red triangle nodes represent lncRNAs, pink circular nodes represent biological processes and square nodes represent mRNAs. OIP5-AS1, OIP5 antisense RNA 1; GAS5, growth arrest specific 5 (non-protein coding); TUG1, taurine upregulated 1 (non-protein coding).

Figure 6.

Analysis of the associated molecular functions of PRKCQ-AS1 and HCP5 in the transformation of diffuse large B-cell lymphoma. Triangle nodes represent lncRNAs, circular nodes represent biological processes and square nodes represent mRNAs. PRKCQ-AS1, PRKCQ antisense RNA; HCP5, HLA complex P5 (non-protein coding).

Discussion

The molecular mechanisms involved in the transformation of DLBCL had previously been unclear. Therefore, it was critically important to investigate the biological mechanisms of DLBCL. In the present study, differentially expressed mRNAs and lncRNAs between FL and DLBCL were identified using the GEO database accession no. GSE53820. Subsequently, a specific DLBCL-associated ceRNA network and a specific FL-associated ceRNA network were constructed. GO and KEGG pathway analyses revealed that differentially expressed lncRNAs served key functions in regulating signal transduction, transcription, cell adhesion, development and protein amino acid phosphorylation.

DLBCL is a malignancy with a high mortality rate due to a lack of biomarkers for early diagnosis and efficient therapeutic strategies (3). Previously, studies had indicated that lncRNAs served key functions in tumorigenesis, cancer progression and metastasis (3,4). An increasing number of studies have additionally demonstrated that the expression of lncRNAs may be deregulated in various types of human cancer, including DLBCL (9,10,17). In prostate cancer, Crea et al (29) identified prostate cancer associated transcript 18 as a novel biomarker and potential therapeutic target for metastatic prostate cancer. Wan et al (30) also reported that androgen-responsive lncRNAs may function as biomarkers for prostate cancer. In the present study, differentially expressed mRNAs and lncRNAs between FL and DLBCL were identified using a publicly available gene expression database, GSE53820. From the microarray expression profiles, it was identified that 1,654 genes were upregulated and 1927 genes were downregulated in DLBCL compared with FL. It was also revealed that 152 lncRNAs were significantly upregulated, and 37 lncRNAs were significantly downregulated between the DLBCL and FL groups.

In order to predict the functions of the differentially expressed lncRNAs, co-expression networks were constructed and GO and KEGG analysis was performed for each lncRNA by using a set of co-expressed mRNAs. According to the KEGG pathway analysis, upregulated lncRNAs were primarily enriched in pathways associated with the cell cycle, T cell receptor signaling pathway, pyrimidine metabolism and purine metabolism. Downregulated lncRNAs were enriched in pathways associated with cell cycle, purine metabolism, aminoacyl-tRNA biosynthesis and degradation of valine, leucine and isoleucine. GO analysis revealed that the upregulated lncRNAs were enriched in transcription, cell cycle, cell division, mitosis, and protein amino acid phosphorylation, whilst the downregulated lncRNAs were enriched in transcription, transport, cell cycle, interspecies interaction between organisms and oxidation reduction.

Previously, a number of reports had revealed that the altered expression of certain lncRNAs may be an important mechanism of DLBCL progression. A number of lncRNAs, including HOTAIR (17), LincRNA-p21 (18) and PEG10 (19) were significantly associated with the progression of DLBCL. However, the molecular mechanisms and functions underlying the involvement of lncRNAs in the transformation of DLBCL remain largely unknown. In the present study, in order to investigate the molecular mechanisms involved in the regulation of DLBCL progression by lncRNAs, lncRNA-miRNA-mRNA ceRNA networks were constructed based on our analysis. From the present study, it was revealed that TUG1, PVT1, MALAT1 and HCP5 served key functions in lncRNA-mediated ceRNA networks. According to GO analysis, the molecular functions of TUG1, PVT1, MALAT1 and HCP5 in DLBCL were investigated. According to the ceRNA networks constructed in the present study, it was revealed that PRKCQ-AS1, HCP5, OIP5-AS1, GAS5 and TUG1 functioned as key regulators. However, the molecular functions of PRKCQ-AS1, HCP5, OIP5-AS1, GAS5 and TUG1 in the transformation of DLBCL remained unknown. By analyzing co-expressed mRNAs, it was revealed that TUG1, OIP5-AS1 and GAS5 were associated with anti-apoptosis, cell cycle, DNA repair, mitosis, transcription, mitosis, G2/M transition of mitotic cell cycle and protein amino acid phosphorylation. PRKCQ-AS1 was associated with the response to oxidative stress, regulation of smooth muscle cell proliferation and acute-phase response. HCP5 was associated with transcription, cell adhesion, lipid metabolism and immune response.

In conclusion, differently expressed lncRNAs between FL and DLBCL were identified for the first time, screened by using a microarray. Compared with FL, a total of 123 upregulated lncRNAs and 192 downregulated lncRNAs in DLBCL were identified. Subsequently, a specific DLBCL-associated ceRNA network and a specific FL-associated ceRNA network were constructed. GO and KEGG pathway analyses revealed that differentially expressed lncRNAs served key functions in regulating signal transduction, transcription, cell adhesion, development and protein amino acid phosphorylation. The present study would provide a potential novel therapeutic and prognostic target for the treatment of DLBCL.