Analysis of NFKB2‑mediated regulation of mechanisms underlying the development of Hodgkin's lymphoma.

Nuclear factor‑κB (NF‑κB) is widely involved in various lymphoid malignancies. However, its exact functional role and potential regulatory mechanisms in Hodgkin's lymphoma (HL) remains unclear. The present study aimed to investigate the regulatory mechanism of NF‑κB in HL by analysis of a gene expression profile that was obtained from HL cells with or without NF‑κB subunit 2 (NFKB2) knockdown. The GSE64234 dataset containing 6 HL cell line specimens transfected with small interfering (si)RNA against NFKB2 and 6 control specimens transfected with non‑targeting siRNA sequences was downloaded from the Gene Expression Omnibus database. Based on these data, differentially expressed genes (DEGs) were screened for following data preprocessing. Functional enrichment analysis was subsequently conducted among the identified upregulated and downregulated DEGs. Additionally, a protein‑protein interaction (PPI) network was constructed and module analyses were performed. Finally, microRNAs (miRNAs/miRs) targeting the identified DEGs were predicted for the construction of a miRNA‑target regulatory network. A total of 253 DEGs were identified, consisting of 109 upregulated and 144 downregulated DEGs. Pathway enrichment analysis revealed that B‑cell lymphoma 2‑like 1 (BCL2L1) was significantly enriched in the NF‑κB signaling pathway, and colony‑stimulating factor 2 (CSF2) and BCL2L1 were enriched in the Jak‑signal transducer and activator of transcription (STAT) signaling pathway. BCL2L1 and CSF2 were determined to be hub genes in the PPI network. A total of 6 miRNAs, including let‑7a‑5p, miR‑9‑5p, miR‑155‑5p, miR‑135a‑5p, miR‑17‑5p and miR‑375, were identified in the miRNA‑target regulatory network. The results of the present study indicated that NFKB2 may be involved in HL development through regulation of BCL2L1, CSF2, miR‑135a‑5p, miR‑155‑5p and miR‑9‑5p expression, as well as the modulation of Jak‑STAT and NF‑κB signaling pathways.

Introduction

Hodgkin's lymphoma (HL) is a common malignant tumor among adolescents (1) that is divided into classical and nodular lymphocyte predominant variants based on the 2008 World Health Organization classification (2). Using chemotherapy and radiotherapy, alone or in combination, the durable remission rate for HL is ~60–80% (3,4). However, survivors have an increased risk of long-term sequelae, including cardiovascular disease and secondary malignancy (5,6). Therefore, understanding the underlying molecular mechanisms of HL development is essential for the discovery of novel therapeutic targets for HL.

Previous studies have demonstrated that the development and maintenance of various solid tumors and lymphoid malignancies is driven by nuclear factor-κB (NF-κB) (7,8). Constitutive NF-κB activation has been identified in classical HL (9). A recent study by de Oliveira et al (10) investigated the functional roles and specific transcriptomes of NF-κB dimers with distinct subunit compositions by knockdown of NF-kB subunit 2 (NFKB2) in a HL cell line; it was reported that NFKB2 knockdown upregulated genes associated with hematopoietic and lymphoid organ development. NFKB2 is a member of the transcription factor complex NF-κB family, which promotes tumor cell proliferation and inhibits apoptosis (11). Truncations in the C-terminal region of NKFB2 are closely associated with tumorigenesis in various hematopoietic tumors, including multiple myeloma, chronic lymphocytic leukemia and cutaneous T-cell lymphoma (12,13). A clear understanding of the regulatory role of NFKB2 in HL may contribute to the development of HL therapies.

de Oliveira et al (10) only analyzed the differentially expressed genes (DEGs) associated with NFKB2 knockdown in a HL cell line, and the regulatory mechanism of NFKB2 in HL remains poorly understood. Therefore, in the present study, the regulatory mechanism of NFKB2 in the tumorigenesis of HL was investigated by downloading the gene expression data (GSE64234) deposited by de Oliveira et al (10). DEGs between the NFKB2 knockdown and control groups were screened for prior to functional enrichment analysis of the identified DEGs, followed by protein-protein interaction (PPI) network and module analyses. Furthermore, microRNA (miRNA/miR) targets were predicted for the identified DEGs and a miRNA-target regulatory network was constructed.

Materials and methods

Affymetrix microarray data

The GSE64234 gene expression profile dataset deposited by de Oliveira et al (10) was downloaded from the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/) database, which was based on the GPL6244 Affymetrix Human Gene 1.0 ST Array (Thermo Fisher Scientific, Inc., Waltham, MA, USA). Data on 12 HL cell line specimens were available in GSE64234, including 6 treated with non-targeting small interfering (si)RNA (control group) and 6 with NFKB2-targeted siRNA (knockdown group).

Data preprocessing

Probe data in CEL files were transformed into expression profiles. If numerous probes corresponded to the same gene, the mean expression value of the probes was used as the gene expression value. The oligo package (version 1.34.0) (14) in R (version 3.3.2) (15) was used to perform raw data format conversion, missing value imputation and background correction normalization for expression profile data. Quartiles were used for data normalization.

DEG screening

The gene expression matrix derived from the knockdown and control groups was analyzed to identify the DEGs. P-values for genes were calculated with the non-paired t-test using the limma package (version 3.10.3) (16) in R and subsequently adjusted using the Benjamin and Hochberg method. An adjusted P-value of <0.05 was set as the threshold value. To directly assess DEG expression in these samples and confirm the classification, a DEG heat map was drawn using the pheatmap package (version 1.0.8) (17) in R.

Functional enrichment analyses

Upregulated and downregulated DEGs were evaluated based on Gene Ontology (GO) (18) and Kyoto Encyclopedia of Genes Genomes (KEGG) (19) pathway enrichment analyses, using the ‘gene set function-functional enrichment-mRNA enrichment’ module in the online Multifaceted Analysis Tool for Human Transcriptome (www.biocloudservice.com). Functional enrichment analyses were based on Fisher's exact test and P<0.05 was selected as the threshold value.

PPI network and subnetwork module analyses

PPI pairs were predicted by the Search Tool for the Retrieval of Interacting Genes database (version 10.0; http://www.string-db.org/) (20) with a confidence score of >0.4. Using the obtained PPI pairs, the PPI network was constructed with Cytoscape version 3.2.0 (21). All nodes in the PPI network were DEGs. In the PPI network, the node score was calculated based on degree centrality. Nodes with a high score, which likely have an important role in the PPI network, were referred to as hub proteins. Proteins typically interact with other proteins in a module. Therefore, following PPI network analysis, significant modules from the PPI network were analyzed using the MultiContrast Delayed Enhancement (MCODE) version 1.4.2 (22) plugin in Cytoscape (23).

miRNA-target regulatory network analysis

miRNAs have a major role in the regulation of biological processes in organisms (24). Therefore, miRNAs that may target the identified DEGs were predicted using the ‘gene-miRNA interaction information retrieval system’ in miRWalk 2.0 (25) (http://mirwalk.uni-hd.de/), with a minimum seed length of seven and statistical significance set at P<0.05. miRWalk is a comprehensive database that not only documents miRNA binding sites within the complete sequence of a gene, but also combines this information with a comparison of binding sites resulting from 12 existing miRNA-target prediction programs (25). In this study, the identified miRNA-target gene pairs required prediction by miRWalk, which included miRanda, miRDB 4.0, miRmap and TargetScan 6.2 databases. Predicted miRNAs were reviewed using the ‘holistic view of validated disease-miRNA interactions’ system in miRWalk 2.0 to select miRNAs that were confirmed to be associated with HL. These verified miRNA-target pairs were used to construct the regulatory network in Cytoscape (21).

Results

A total of 253 DEGs were identified following data preprocessing, which included 109 upregulated DEGs and 144 downregulated DEGs. The heat map clearly indicated the distinct separation of the knockdown and control groups (Fig. 1).

Figure 1.

Heat map analysis of differentially expressed genes in the siNT and siNFKB2 groups. Red indicates high expression values and green indicates low expression values. si, small interfering RNA; siNT, non-targeting siRNA; NFKB2, nuclear factor-κB subunit 2; siNFKB2, NFKB2-targeting siRNA.

The upregulated DEGs were significantly enriched in ‘platelet activation’, ‘osteoclast differentiation’ and ‘long-term depression’ KEGG pathways, while downregulated DEGs were predominantly enriched in ‘intestinal immune network for IgA production’, ‘HTLV-I infection’ and ‘NF-κB signaling pathway’ KEGG pathways (Table I).

Table I.

Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis for upregulated and downregulated DEGs.

A, Upregulated DEGs
ID	Name	P-value	Genes
hsa04611	Platelet activation	0.03878	PLA2G4F, FCGR2A, PRKG2, PRKG1
hsa04380	Osteoclast differentiation	0.03953	TYK2, GAB2, FCGR2A, MAPK10
hsa04730	Long-term depression	0.04684	PLA2G4F, PRKG2, PRKG1
B, Downregulated DEGs
ID	Name	P-value	Genes
hsa04940	Type I diabetes mellitus	0.0003212	CD80, HLA-DRB5, IL12B, HLA-DPB1, LTA
hsa04672	Intestinal immune network for IgA production	0.0004971	CD80, HLA-DRB5, TNFRSF17, AICDA, HLA-DPB1
hsa05330	Allograft rejection	0.003062	CD80, HLA-DRB5, IL12B, HLA-DPB1
hsa05166	HTLV-I infection	0.003913	WNT5A, CSF2, TLN1, HLA-DRB5, NFKB2, BCL2L1, HLA-DPB1, LTA
hsa04064	NF-κB signaling pathway	0.004905	NFKB2, BCL2L1, BIRC2, PLAU, LTA
hsa04060	Cytokine-cytokine receptor interaction	0.009454	CSF2, CCR4, TNFRSF17, TNFRSF8, TNFRSF14, IL12B, LTA
hsa05145	Toxoplasmosis	0.01411	HLA-DRB5, BCL2L1, IL12B, HLA-DPB1, BIRC2
hsa05168	Herpes simplex infection	0.01472	NECTIN1, HLA-DRB5, TNFRSF14, IL12B, HLA-DPB1, LTA
hsa04062	Chemokine signaling pathway	0.0157	CCL22, CCR4, PTK2B, HCK, GRK2, PRKCD
hsa04514	Cell adhesion molecules (CAMs)	0.02597	ALCAM, CD80, NECTIN1, HLA-DRB5, HLA-DPB1
hsa05332	Graft-versus-host disease	0.02806	CD80, HLA-DRB5, HLA-DPB1
hsa05323	Rheumatoid arthritis	0.03255	CSF2, CD80, HLA-DRB5, HLA-DPB1
hsa05202	Transcriptional misregulation in cancer	0.04407	CSF2, REL, RARA, BCL2L1, PLAU

DEGs, differentially expressed genes.

Additionally, upregulated DEGs were significantly enriched in 10 GO terms, including ‘integral component of membrane’, ‘plasma membrane’ and ‘Golgi apparatus’ terms (Fig. 2A). Downregulated DEGs were significantly enriched in 53 GO terms, including ‘protein binding’, ‘cytoplasm’, ‘transcription factor activity’ and ‘sequence-specific DNA binding’. The top 15 downregulated GO terms are presented in Fig. 2B.

Figure 2.

GO enrichment analysis for DEGs in BP, CC and MF terms. The number on the x-axis indicates the enriched count of differentially expressed genes. (A) Enriched BP, CC and MF GO terms in the upregulated DEGs. (B) Top 15 enriched GO terms (BP, CC and MF) for downregulated DEGs. GO, Gene Ontology; DEGs, differentially expressed genes; BP, biological process; CC, cellular component; MF, molecular function.

PPI network analysis

The PPI network included 142 nodes and 220 PPI pairs (Fig. 3). Based on degree centrality analysis, 8 upregulated and 12 downregulated DEGs were determined to be hub genes (Table II). KEGG pathway enrichment analysis was subsequently performed to investigate the function of the 20 identified hub genes. As presented in Table III, these hub genes were significantly enriched in 14 KEGG pathways. For instance, B-cell lymphoma 2-like 1 (BCL2L1) was significantly enriched in ‘NF-κB signaling pathway’ and ‘Jak-STAT signaling pathway’. Colony-stimulating factor 2 (CSF2) was significantly enriched in ‘transcriptional misregulation in cancer’, ‘TNF signaling pathway’ and ‘Jak-STAT signaling pathway’. Tyrosine kinase 2 was significantly enriched in ‘Jak-STAT signaling pathway’, ‘toxoplasmosis’, ‘Epstein-Barr virus infection’ and ‘osteoclast differentiation’. CD44 antigen was enriched in ‘Epstein-Barr virus infection’.

Figure 3.

Protein-protein interaction network analysis of DEGs. Red circles represent upregulated DEGs and green diamonds represent downregulated DEGs. DEGs, differentially expressed genes.

Table II.

Degree centrality analysis for 20 hub genes.

Node	Description	Degree
BCL2L1	Downregulated	19
CDH1	Upregulated	13
CAT	Upregulated	12
TYK2	Upregulated	12
REL	Downregulated	11
PTK2	Upregulated	10
CSF2	Downregulated	9
HCK	Downregulated	8
BIRC2	Downregulated	8
NFKB2	Downregulated	8
WDFY2	Upregulated	7
MAPK10	Upregulated	7
PRKCD	Downregulated	7
ATM	Upregulated	6
TRIP10	Downregulated	6
CD44	Upregulated	6
PTK2B	Downregulated	6
FMOD	Downregulated	5
FNBP1	Downregulated	5
IER3	Downregulated	5

Table III.

Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis for the 20 hub genes.

ID	Name	P-value	Genes
hsa05202	Transcriptional misregulation in cancer	0.0003731	CSF2, PTK2, REL, BCL2L1, ATM
hsa04064	NF-κB signaling pathway	0.0007863	NFKB2, BCL2L1, BIRC2, ATM
hsa05200	Pathways in cancer	0.001083	PTK2, CDH1, NFKB2, BCL2L1, MAPK10, BIRC2
hsa05145	Toxoplasmosis	0.001901	TYK2, BCL2L1, MAPK10, BIRC2
hsa04062	Chemokine signaling pathway	0.006877	PTK2, PTK2B, HCK, PRKCD
hsa05169	Epstein-Barr virus infection	0.007294	TYK2, CD44, NFKB2, MAPK10
hsa04210	Apoptosis	0.007715	BCL2L1, BIRC2, ATM
hsa05222	Small cell lung cancer	0.01415	PTK2, BCL2L1, BIRC2
hsa04912	GnRH signaling pathway	0.01611	PTK2B, MAPK10, PRKCD
hsa05166	HTLV-I infection	0.01643	CSF2, NFKB2, BCL2L1, ATM
hsa04668	TNF signaling pathway	0.02149	CSF2, MAPK10, BIRC2
hsa04380	Osteoclast differentiation	0.03187	TYK2, NFKB2, MAPK10
hsa04068	FoxO signaling pathway	0.03323	CAT, MAPK10, ATM
hsa04630	Jak-STAT signaling pathway	0.03840	TYK2, CSF2, BCL2L1

Subnetwork module analysis

Module analysis identified three modules (termed A, B and C; Fig. 4). Module A was comprised of 5 nodes, 10 interaction pairs and was enriched in ‘amyotrophic lateral sclerosis (ALS)’ and ‘NF-κB signaling pathway’ KEGG pathways. There were 4 nodes and 6 interaction pairs in module B, without KEGG pathway enrichment. Module C contained 4 nodes and 6 interaction pairs and was enriched in the ‘hematopoietic cell lineage’ KEGG pathway.

Figure 4.

Subnetwork module analysis of DEGs identified three modules. Red circles represent upregulated DEGs and green diamonds represent DEGs. DEGs, differentially expressed genes.

The miRNA-target regulatory network was constructed with 32 nodes and 29 PPI pairs (Fig. 5). Among the 32 nodes, there were 16 upregulated DEGs, 10 downregulated DEGs and 6 miRNAs, which included let-7a-5p, miR-9-5p, miR-155-5p, miR-135a-5p, miR-17-5p and miR-375.

Figure 5.

miRNA-target regulatory network analysis. Blue triangles represent miRNA, red circles represent upregulated DEGs, green diamonds represent downregulated DEGs and the arrows indicate the direction of regulation. miRNA, microRNA; DEGs, differentially expressed genes.

Discussion

In the present study, a total of 109 upregulated and 144 downregulated DEGs were identified between the NFKB2 knockdown and control groups. Pathway enrichment analysis revealed that BCL2L1 was significantly enriched in the NF-κB signaling pathway, and CSF2 and BCL2L1 were enriched in the Jak-STAT signaling pathway. Additionally, BCL2L1 and CSF2 were hub genes in the PPI network. Prediction of miRNA-target interactions identified 29 miRNA-target pairs involving in 6 miRNAs (let-7a-5p, miR-9-5p, miR-155-5p, miR-135a-5p, miR-17-5p and miR-375). The findings of the present study may provide guidelines for understanding the regulatory mechanism of NFKB2 and contribute to the identification of therapeutic targets for HL.

BCL2L1 was demonstrated to be downregulated in HL cell lines following NFKB2 knockdown. This gene is a member of the Bcl-2 family and is a vital apoptosis-mediating gene that encodes both an antiapoptotic and proapoptotic splice variant (26). Additionally, in the present study, module analysis revealed that BCL2L1 interacts with NFKB2 in module A, which indicated that knockdown of NFKB2 may suppress BCL2L1 expression and subsequently regulate HL cell apoptosis. Notably, this result was consistent with the findings of de Oliveira et al (10), who reported that NFKB2 knockdown may downregulate BCL2L1 expression. Furthermore, BCL2L1 was enriched in the NF-κB signaling pathway in the present study. The NF-κB signaling pathway has been associated with the progression of several tumors, including HL (7,11). Therefore, it was speculated that BCL2L1 may be involved in NF-κB-dependent control of HL cell survival.

In addition to BCL2L1, CSF2 was also downregulated in the NFKB2 knockdown group. CSF2 encodes a cytokine that has been reported to be involved in the pathogenesis of HL (27) and is highly expressed in HL cells (28). In accordance with these findings, the present study demonstrated that CSF2 was upregulated in HL control cells compared with HL cells with NFKB2 knockdown, indicating that NFKB2 may be involved in HL tumorigenesis through the upregulation of CSF2. KEGG pathway enrichment analysis revealed that CSF2 was enriched in the Jak-STAT signaling pathway. The Jak-STAT signaling pathway is regulated by a large array of cytokines and growth factors that induce proliferation or differentiation (29). Activation of the Jak-STAT signaling pathway is a hallmark of various B-cell lymphomas, including classical HL (30). Taken together, these data indicate that NFKB2 may be implicated HL development by regulating CSF2 expression, which in turn activates the Jak-STAT signaling pathway.

miRNAs regulate gene expression and cellular processes (24). They are also potential diagnostic and prognostic molecular markers in cancer, as well as targets for the development of highly specific therapies (31). Based on the DEGs identified in the present study, several miRNAs were predicted. Among these, miR-135a-5p has been previously reported to have critical involvement in HL cell survival by targeting Jak2 (32). Additionally, miR-155-5p has been recommended as a diagnostic marker in hematological tumors (33). Notably, dysregulated expression of miR-155 has been reported in various types of lymphoma (33). Differential expression of miR-9-5p has been demonstrated in hematological cell lines compared with normal lymphocyte populations (34). Furthermore, a target gene of miR-9-5p, PR/SET domain 1, was reported to be downregulated in HL cells (35), further indicating a potential role for miR-9-5p in HL. Therefore, NFKB2 may be implicated in HL development through interaction with these miRNAs.

In conclusion, the present study demonstrated that NFKB2 may be involved in the development of HL by interacting with several genes and miRNAs, including BCL2L1, CSF2, miR-135a-5p, miR-155-5p and miR-9-5p. Additionally, NFKB2-mediated modulation of the Jak-STAT and NF-κB signaling pathways may have a role in HL tumorigenesis. The genes, miRNAs and signaling pathways identified in the present study may serve as promising therapeutic targets for HL. Future experimental verification in vivo is required to confirm the findings of the present study.