Identification of key miRNA-gene pairs in chronic lymphocytic leukemia through integrated analysis of mRNA and miRNA microarray.

The aim of the present study was to explore the miRNA-Gene regulatory mechanism in chronic lymphocytic leukemia (CLL), and identify new targets for the therapy of CLL. The miRNA expression dataset GSE62137 and mRNA expression dataset GSE22529 were downloaded from National Center of Biotechnology Information Gene Expression Omnibus database. In CLL samples compared with normal B cell samples, differentially expressed miRNAs (DEMs) were identified via the GEO2R instrument of GEO and differentially expressed genes (DEGs) were obtained via the limma package of R. Functional enrichment analysis of the DEGs was performed via the Database for Annotation, Visualization and Integrated Discovery. The targets of the DEMs were identified based on the miRNAWalk platform. The overlaps between the DEGs and the targets of the DEMs were selected, and the miRNA-Gene regulatory network was constructed based on the overlaps and the corresponding DEMs. A total of 63 DEMs and 504 DEGs were identified in CLL samples compared with normal B cell samples. Eleven enriched functional clusters of the DEGs were obtained. 405 miRNA-Gene regulatory pairs were identified. The miRNA-Gene regulatory pairs contained 351 target genes of the DEMs, including 9 overlaps with the DEGs. A miRNA-Gene regulatory network was constructed. Bioinformatics methods could help us develop a better understanding of the molecular mechanism of CLL. MiRNAs may play a critical role in regulating the process of CLL. They may affect CLL by regulating the processes of immunoreactivity and protein degradation. Genes such as Neurogenic Locus Notch Homolog Protein 2, PR/SET domain 4 and A-kinase anchoring protein 12 may be their regulating targets in CLL.

Introduction

B-cell chronic lymphocytic leukemia (CLL) accounts for about 25% of all leukemia and is the most common form of lymphoid malignancy in western countries (1). CLL is characterized by the accumulation of the clonal expansion and resistance to apoptosis of immature CD5(+) B cells (2). The disease displays a high heterogeneity in its clinical course (3,4). Staging and prognostication of CLL is performed by two similar clinical staging systems (5,6). However, the clinical staging systems cannot fully reflect the high variability of CLL and predict survival of CLL patients, as well as the response to therapy (7). Molecular aberration plays an important role in the initiation and progression of CLL. Several factors that can predict the clinical course have been identified (8–10). The most frequent deletion of genomic DNA in CLL occurs in chromosome 13q13.4. This deletion occurs in ~50% of CLL cases and is associated with a long interval between diagnosis and the need for treatment (11). One study was reported that an increased DPF3 expression in CLL patients that was linked with a significant increase of STAT5 activation in myeloid lineage cells (granulocytes) and not in neoplastic B cells (12). However, the molecular mechanism of these associations is still largely unknown.

miRNAs are small noncoding RNA molecules that negatively regulate gene expression via degradation or translational repression of their targeted mRNAs (13). miRNAs played important roles in the pathogenesis of various human cancers and ~50% of the known human microRNAs were located at cancer-associated regions of the genome (14–16). In addition, miRNAs played critical roles in the development and oncogenesis of B-cell. They were involved in the pathogenesis of CLL (17). MiRNAs expressions were associated with clinical characteristics of CLL, and they were likely to be served as diagnostic and prognostic biomarkers, as well as potential therapeutic targets in CLL (18,19). Negrini et al revealed that miRNAs whose expression was distinctive of cases with mutated vs. unmutated IGHV genes or cases with 13q, 11q, and 17p deletions and trisomy 12 and miR26a, miR532-3p, miR146-5p, and miR29c* were strongly associated with progression-free survival in CLL (20). However, the specific regulatory mechanism of miRNA in CLL was still needed to further explore.

In the current studies usually analyzed gene expression by high throughput RNA sequenc, though miRNA is also known as important in the regulation of gene expression and cancer formation. In this study, we firstly identified the key miRNA-Gene pairs in CLL.

In this study, differentially expressed genes (DEGs) and miRNAs (DEMs) in CLL samples compared with normal B cell samples were identified via bioinformatics methods. The miRNA-Gene regulatory network in CLL was constructed. A better understanding of the regulatory mechanism of miRNAs in CLL was obtained. Our study may provide references for the diagnosis and therapy of CLL.

Materials and methods

Microarray data

The miRNA expression dataset GSE62137 (21) and mRNA expression dataset GSE22529 (22) were downloaded from National Center of Biotechnology Information (NCBI) Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) database. The miRNA dataset GSE62137 contained 44 samples, including 38 CLL cell samples and 6 normal B cell samples. Some of these samples were treated with IL-4. Our study was performed based on the untreated samples (23 CLL cell samples and 3 normal B cell samples cultured with nothing). MiRNA expression profile was detected via the Agilent-021827 Human miRNA Microarray G4470C platform. The mRNA dataset GSE22529 contained 52 samples, including 41 CLL cell samples and 11 normal B cell samples. The microarray data we used was detected based on GPL96 [HG-U133A] Affymetrix Human Genome U133A Array platform.

Identification of differently expressed miRNAs and genes

For the miRNA dataset, DEMs in CLL samples compared with normal B cell samples were identified via the GEO2R application of GEO. The screening threshold was FDR corrected P<0.05 and |log2 (fold-change)|>1. For the mRNA dataset, the raw data with CEL files were background corrected, normalized and log2 transformed using the affy package in R. (23). If multiple probes correspond to one gene, the mean expression value was defined as expression value. DEGs in CLL samples compared with normal B cell samples were identified via the limma package (24) of R. The DEGs were identified according to the criteria of |log2 (fold-change) |>1 and FDR corrected P<0.05.

Functional clustering analysis

The Database for Annotation, Visualization and Integrated Discovery (DAVID) (https://david.ncifcrf.gov/) was a widely used web-based tool for genomic functional annotations (25). In the present study, Gene Ontology (GO) terms enrichment analysis were performed via DAVID. Enrichment score >1 was used as the criteria to identify significant functional clusters.

Screening of target genes of the DEMs

The targets of the DEMs, as well as the miRNA-Gene pairs, were identified based on the miRNAWalk platform. Five databases: miRnada (26), miRDB (27), miRWalk (28), RNA22 (29) and TargetScan (30) were used in this study to predicted the targets of the DEMs. Only the targets that simultaneously appeared in the five databases were screened out.

Construction of the miRNA-Gene regulatory network

The overlaps between the DEGs and the targets of the DEMs were selected. The miRNA-Gene regulatory network was constructed based on the miRNA-Gene pairs, and then visualized via Cytoscape software (version 3.4.0; www.cytoscape.org).

Results

The DEMs and the DEGs

A total of 63 DEMs were identified in CLL samples compared with normal B cell samples, including 51 downregulated ones and 12 upregulated ones. The top 20 DEMs were listed in Table I. For the mRNA dataset, gene expression values of each sample after normalization were shown in Fig. 1. 504 DEGs were identified in CLL samples compared with normal B cell samples, including 316 downregulated ones and 188 upregulated ones. Table I showed the top 20 DEGs. Cluster analysis of CLL samples and normal B cell samples based on the DEGs expression values was shown in Fig. 2. From the heatmap (Fig. 2), we found that the gene expression of 10 normal B cell samples were distinguished from the 41 CLL samples. One normal B cell samples was assigned to CLL samples. Accuracy rate of prediction is 98.07%.

Table I.

The top 20 DEMs in CLL samples compared with normal B cell samples.

MiRNA ID	P-value	LogFC
hsa-miR-582-5p	5.32×10−17	−8.23
hsa-miR-181a	8.86×10−04	−8.56
hsa-miR-132	8.86×10−04	−6.42
hsa-miR-95	8.86×10−04	−4.34
hsa-miR-548c-3p	8.86×10−04	−3.65
hsa-miR-181c	8.86×10−04	−3.22
hsa-miR-150	8.86×10−04	1.54
hsa-miR-486-5p	8.86×10−04	5.65
hsa-miR-451	8.86×10−04	6.43
hsa-miR-144	8.86×10−04	8.56
hsa-miR-28-5p	1.81×10−03	1.57
hsa-miR-885-3p	1.42×10−02	−3.2
hsa-miR-199a-3p	1.71×10−02	−5.62
hsa-miR-155	1.74×10−02	2.00
hsa-miR-126	2.22×10−02	−6.29
hsa-miR-29a	2.26×10−02	1.27
hsa-miR-21	2.32×10−02	1.42
hsa-miR-202	2.47×10−02	−3.53
hsa-miR-29b	3.93×10−02	1.15
hsa-miR-199a-5p	3.95×10−02	−3.42

miRNA, microRNA; DEMs, differentially expressed miRNAs; FC, fold-change; CLL, chronic lymphocytic leukemia.

Figure 1.

Gene expression value of each sample in GSE22529 after normalization.

Figure 2.

Cluster analysis of chronic lymphocytic leukemia samples and normal B cell samples based on the differentially expressed genes. Lighter red in the heat-map represents high expression and darker green indicates low expression, black denotes medial expression.

Enriched functional clusters

11 enriched functional clusters of the DEGs were obtained. The enriched GO terms in the top 2 clusters (cluster 1 and cluster 2) were shown in Fig. 3. The GO terms of clusters 1 were significantly enriched in the process of immune response, including B cell receptor signaling pathway (GO:0050853), complement activation classical pathway (GO:0006958), antigen binding (GO:0003823) and so on. The GO terms of clusters 2 were mainly enriched in the process of immune recognition receptors includes complement activation classical pathway (GO:0006958), antigen binding (GO:0003823), receptor-mediated endocytosis (GO:0006898) and so on. The most enriched GO term in these 2 clusters involve an immune response.

Figure 3.

The enriched Gene Ontology terms in the top 2 functional clusters (cluster 1 and cluster 2).

The miRNA-gene regulatory network

A total of 405 miRNA-Gene regulatory pairs were identified via the miRNAWalk platform. The miRNA-Gene regulatory pairs contained 351 target genes of the DEMs, including 9 overlaps with the DEGs. These overlaps were regulated by 10 DEMs, and then 10 miRNA-Gene regulatory pairs were obtained. Afterwards, a regulatory network between these overlaps and the DEMs were constructed and visualized (Fig. 4). As seen in Fig. 4, 9 target genes differentially expressed (TRAM2, BASP1, PRDM4, DGKG, TOP1, JOSD1, NOTCH2, AKAP12, CNTN6) are regulated by 9 different miRNAs (hsa-miR-138b, hsa-miR-352, hsa-miR-1, hsa-miR-302b, hsa-miR-326, hsa-miR-136, hsa-miR-181c, hsa-miR-145, hsa-miR-150).

Figure 4.

The miRNA-Gene regulatory network in chronic lymphocytic leukemia. The rectangle represents the gene, and the rhombus represents the miRNA. Has, homo sapiens; miRNA, microRNA; Notch2, neurogenic locus notch homolog protein 2; PRDM4, PR/SET domain 4; AKAP12, A-kinase anchoring protein 12; BASP1, brain abundant membrane attached signal protein 1; DGKG, diacylglycerol kinase γ; CMTM6, CKLF like MARVEL transmembrane domain containing 6; TRAM2, translocation associated membrane protein 2; TOP1, topoisomerase (DNA) I; JOSD1, josephin domain containing 1.

Discussion

CLL is the most common form of leukemia in western countries with the incidence of approximately 1 out of 100,000 patients per year (31), but little is known about its initiation and progression (32). MiRNAs play critical roles in B-cell oncogenesis by regulating the expression of many genes, and they are involved in the pathogenesis of CLL. In this study, a miRNA-Gene regulatory network in CLL was constructed via bioinformatics methods, which would help us make a better understanding of the molecular mechanism of CLL.

In the enriched GO terms in cluster 1, immunoglobulin complex, circulating and positive regulation of B cell activation were most closely related to the initiation and progression of CLL. While, in the enriched GO terms in cluster 2, proteolysis and serine-type endopeptidase activity had the most closely relationship with CLL. The mutational status of immunoglobulin heavy-chain variable-region (IgVH) genes in the leukemic cells of CLL was an important prognostic factor in the disease. Approximately 50–70% of patients with CLL had somatic hypermutation in IgVH genes of the leukemic cells (9,33). In addition, immunoglobulin could also be used as a therapeutic agent for CLL, and intravenous immunoglobulin G may alter the response of CLL cells to chemotherapy (34). B cell activation was an important initiation factor of CLL. CLL was characterized by the progressive accumulation of clonal mature B cells in the blood, bone marrow, and secondary lymphoid organs (35). B cell receptor signaling represented one of the central pathways to enhance survival and proliferation in CLL (10). One research showed that differences existed in the activation of B-CLL cells in vivo, and these patterns were correlated with disease activity (36). Although the importance of protein proteolysis in biological functions was frequently unclear, new technologies had started to unravel the critical role of clipping in cellular homeostasis and disease (37,38). Specific histone H2A proteolysis as disruption of the histone code was reported to drive hematopoietic cells in lymphomagenesis and result in lymphoid malignancies (39). Furthermore, matriptase was reported to be highly upregulated in CLL, which promoted cancer cell invasion (40). Serine-type endopeptidase also played critical roles in CLL. A serine endopeptidase, kallikrein B1 (KLKB1) was reported to be overexpressed in CLL, and its expression could be served as a novel molecular biomarker for the diagnosis of CLL (41).

The miRNA-Gene regulatory network contained 9 DEGs; they were BASP1, NOTCH2, CMTM6, DGKG, TRAM2, TOP1, PRDM4, JOSD1 and AKAP12. DAVID is a comprehensive system of biological function annotation information for large-scale genes, so it is too little to perform the DAVID analysis on this 9 overlapping DEGs. According to the European Molecular Biology Laboratory (EMBL) database (http://www.ebi.ac.uk/), we found 9 DEGs were enriched the following GO terms. BASP1, NOTCH2, TOP1, PRDM4, AKAP12 and JOSD1 were included in the biological process of protein binding (GO:0005515) which child terms included antigen binding, immunoglobulin receptor binding. JOSD1 was aslo included in the biological process of proteolysis (GO:0006508). JOSD1, NOTCH2 and AKAP12 wwere included in the biological process of hydrolase activity (GO:0016787), receptor activity (GO:0004872), regulation of protein kinase C signaling (GO:0090036), respectively, which were closely related to complement activation, classical pathway, positive regulation of B cell activation. TRAM2 and DGKG were respectively included in the biological process of integrall component of menbrane (GO:001601), intracellular signal transduction (0035556) which were closely related to signal transduction. To sum up, the GO terms of 9 DEGs were associated with the GO terms of all DEGs. Many of these genes were associated with the pathobiology of CLL. For example, NOTCH2 was a member of the NOTCH gene family. The Notch gene family encodes transmembrane receptors that modulate differentiation, proliferation and apoptotic programs in response to extracellular ligands expressed on neighboring cells (42). Enforced expression of NOTCH1IC in bone marrow stem cells caused T-cell leukemia in mice, indicating a causative role for NOTCH family in T-cell oncogenesis (43). The NOTCH2 oncogene was reported to be overexpressed in B-CLL cells, and was also related to the failure of apoptosis characteristic for this disease. Deregulation of NOTCH2 signaling was involved in the aberrant expression of CD23 in B-CLL (44). The overexpression of CD23 was one of the hallmarks of B-CLL cells (45). PRDM4 was a highly conserved member of the PRDM family (46). One study reported that PRDM4 could control proliferation and differentiation, which played critical roles in tumorigenesis (47). Another study reported that PRDM4 protein mapped to a tumor suppressor locus on human chromosome and could affect the processes of ovarian, gastric, and pancreatic cancers (48). AKAP12, which functioned as a kinase scaffold protein and as a dynamic regulator of the b2-adrenergic receptor complex, was one of the A-kinase anchoring proteins (49). AKAP12 played an important role in regulating cell cycle, cytokinesis cell adhesion, signaling, and oncogenic suppression (50). AKAP12 expression was closely related to tumorigenesis. The downregulation of AKAP12 expression had been reported in human prostate cancers in vivo, suggesting that the inactivation of AKAP12 expression was associated with oncogenesis (51). In gastric cancer, AKAP12 may function as an important negative regulator of the survival pathway (49). In addition, AKAP12 expression could also be served as a predictor for survival in CLL (52).

In conclusion, miRNAs played a critical role in regulating the process of CLL. They could affect CLL by regulating the processes of immunoreaction and protein degradation. Genes such as NOTCH2, PRDM4 and AKAP12 proved be their regulating targets in CLL. These DEGs which were related to CLL, could potentially serve as biomarkers for detection, prognosis, monitoring and predicting therapeutic responses in CLL. However, further studies were still needed to confirm our results.

Table II.

The top 20 DEGs in CLL samples compared with normal B cell samples.

Gene name	P-value	LogFC
LMO2	2.84×10−19	−3.95
GLRX	1.59×10−18	−3.31
IGHV1-69	2.46×10−17	−5.23
DENND3	5.06×10−17	−1.87
CCDC88A	5.06×10−17	4.07
HOPX	9.31×10−16	−2.71
CST3	1.10×10−15	−4.84
IGHV3-48	1.10×10−15	−4.42
PIK3C2B	1.10×10−15	−3.27
CEBPD	1.10×10−15	−2.47
FPR1	1.10×10−15	−1.86
CD14	2.09×10−15	−5.12
CPVL	2.09×10−15	−3.74
IGKV3-20	2.09×10−15	−2.94
CSF1R	2.09×10−15	−2.62
SCN3A	3.03×10−15	−3.69
SLC2A3	6.60×10−15	−3.21
SORL1	1.51×10−14	−5.08
FOS	3.61×10−14	−6.24
SLC2A14	3.72×10−14	−3.32

DEGs, differentially expressed genes; FC, fold-change; CLL, chronic lymphocytic leukemia.