Identification of molecular mechanisms of glutamine in pancreatic cancer.

The aim of the present study was to explore the critical genes and molecular mechanisms in pancreatic cancer (PC) cells with glutamine. By analyzing microarray data GSE17632 from the Gene Expression Omnibus database, the DEGs between PC cells treated with glutamine and without glutamine were evaluated. Additionally, function enrichment analyses and protein-protein interaction (PPI) network construction of DEGs were performed. Network module and literature mining analyses were performed to analyze the critical DEGs in PC cells. In total, 495 genes were selected as DEGs between control and glutamine cells in PC. These DEGs were mainly enriched in several Gene Ontology (GO) terms in biological process, cellular components and molecular function. Additionally, they were also enriched in certain pathways, including metabolic pathways and the Janus kinase-signal transducer and activator of transcription (JAK-STAT) signaling pathway. MYC, heat shock 70kDa protein 5 (HSPA5), interleukin 8 (IL8), and chemokine (C-X-C motif) receptor 4 (CXCR4) were hub genes in the PPI network. Furthermore, two sub-network modules of PPI network and two co-occurrence networks were obtained. The DEGs of MYC, HSPA5, IL18 and CXCR4 may exert important roles in molecular mechanisms of PC cells with glutamine.

Introduction

Pancreatic cancer (PC) is a deadly malignant disease, with an overall 5-year survival rate of <5% (1). It is the 7th leading cause of cancer-associated mortalities worldwide (2) and poses a great threat to the health of individuals. Annually, the mortality rate of PC patients is almost identical to its incidence rate (3). Currently, pancreatectomy remains the most effective therapy modality for PC patients, which offers the only potential for successful treatment. However, the patients undergoing resection treatment have a median survival of only 12–22 months (4). Therefore, it is urgent to understand the molecular pathophysiology of PC to promote the development of effective therapeutic strategies.

In recent decades, considerable efforts have been made to investigate the pathogenesis and therapeutic strategies for the treatment of PC. It is well known that cell growth is controlled by a coordinated response to nutrients and growth factors. Alterations in nutrient sensing and growth factors may lead to cancer incidence (5). Glutamine, as the necessary nutrient in nucleic acid synthesis and cell proliferation, plays an important role in the process of tumor anabolic processes (6,7). Pancreatic ductal adenocarcinoma cells have been found profoundly sensitive to glutamine deprivation, indicating that glutamine is critical for pancreatic ductal adenocarcinoma growth (8). In particular, one study has found that oncogenes could regulate nutrient metabolism in the development of malignancy. MYC, for example, can drive glutamine uptake and catabolism by activating the expression of genes, including glutaminase and solute carrier family 1 (neutral amino acid transporter), member 5 (9). Although several genes associated with glutamine metabolism in PC have been studied, it is far from sufficient to fully understand the molecular mechanisms of PC.

Therefore, in the present study, the expression profile data GSE17632 (5) was assessed to identify the differentially-expressed genes (DEGs) between PC cells treated with glutamine and without glutamine. With these selected DEGs, Gene Ontology (GO) functional and pathway enrichment analyses were performed, and the protein-protein interaction (PPI) network was constructed. Additionally, network module and literature mining analyses were also performed to further study the functions of DEGs. The present study explored the critical genes and molecular mechanisms in PC cells with glutamine by bioinformatics methods.

Materials and methods

Microarray data source

The mRNA expression profile data of GSE17632 were downloaded from the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/) database of National Center for Biotechnology Information. In the dataset, 4 samples (GSM440132, GSM440133, GSM440134 and GSM440135) of PC cells treated with glutamine (glutamine group) compared with PC cells without glutamine treatment (control group) were selected for analysis. The PC cells were BxPC-3 pancreatic cancer cells. The platform was GPL4133 Agilent-014850 Whole Human Genome Microarray 4×44K G4112F (Feature Number version).

Data preprocessing

The dataset were dual-channel chips, including the Cy3 channel and Cy5 channel. The Cy3 channel consisted of control samples, and the growth condition of the PC cells was glucose and glutamine depleted. The Cy5 channel consisted of experimental samples, of which the growth condition was glutamine. The original data was pre-processed using locally weighted scatterplot smoothing (LOWESS) (10) and the pre-processed data were used to analyze the DEGs.

DEG analysis

The DEGs between the glutamine and control groups were analyzed using the Bioconductor limma package (11). The unpaired t-test was used to calculate the P-value and false discovery rate (FDR). Additionally, the fold-change (FC) among the sample groups was also calculated. Only genes with FDR<0.01 and |log2FC|>1 were selected as the DEGs.

GO functional and pathway enrichment analyses

In the present study, the function enrichment analysis of the DEGs were analyzed using GO (http://www.geneontology.org) (12) database, which provided the function annotations of DEGs in biological process (BP), molecular function (MF) and cellular component (CC), respectively. In addition, the pathway enrichment analysis was performed through Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.ad.jp/kegg/) (13) database. During the process of enrichment analyses, the significant threshold of the hypergeometric test was set as 0.05.

PPI network construction and network module analysis

The Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org/) database (14) is a precomputed global resource for evaluating PPI information. In the present study, the STRING database was used to predict the PPI for the DEGs. With a PPI score of 0.7, the PPI network was constructed and was visualized using Cytoscape (15) which was a general bioinformatics package used for visualizing biological network and integrating data.

Additionally, in the present study, the Cluster One (16) plugin in Cytoscape was used for mining modules in the PPI network. The network modules with P<0.001 were selected.

Literature mining of key genes in the network module

Subsequent to analysis of the modules, the DEGs in the modules were analyzed using literature mining to explore their relevance in previous studies. The Gene Cluster with Literature Profiles 2.0 (GenCLiP 2.0; http://ci.smu.edu.cn/) (17) online tool was used for literature mining of human genes and network. The input gene set was the key gene set in the PPI network. The Literature Mining Gene Networks module of GenCLiP was used to construct gene co-occurrence networks of the input genes and to analyze the hotspot-associated genes in the literature. The biological function of hotspot genes were then analyzed by the Gene Cluster with Literature Profiles module, with the parameters of P≤1×10−10 and Hit≥4 (Hit represents the number of articles mentioning the corresponding gene that also contain the search term used).

Results

DEGs analysis

In the present study, a total of 495 genes were selected as DEGs in the glutamine group, including 329 upregulated DEGs and 166 downregulated DEGs.

GO functional and KEGG enrichment analyses of DEGs

The significant enrichment result of DEGs in BP, CC and MF was shown in Table I. The most significant terms of BP, CC and MF enriched by upregulated DEGs were, respectively, GO:0010033 response to organic substance, GO:0044421 extracellular region, and GO:0051787 misfolded protein binding. The downregulated DEGs were mainly enriched in BP terms associated with anatomical structure development, CC terms associated with cell junction and MF terms associated with region sequence-specific DNA binding transcription factor activity.

Table I.

The most significant enrichment results of differentially expressed genes in GO terms.

A, Upregulated
Go ID	Term	Count	P-value
BP
GO:0010033	Response to organic substance	80	2.97×10−11
GO:0006950	Response to stress	99	2.02×10−10
GO:0035966	Response to topologically incorrect protein	17	6.35×10−10
GO:0006986	Response to unfolded protein	16	1.96×10−9
GO:0033993	Response to lipid	34	7.07×10−9
CC
GO:0044421	Extracellular region part	101	6.71×10−11
GO:0005615	Extracellular space	50	1.31×10−10
GO:0005576	Extracellular region	107	2.72×10−8
GO:0070062	Extracellular vesicular exosome	73	3.71×10−7
GO:1903561	Extracellular vesicle	73	3.71×10−7
MF
GO:0051787	Misfolded protein binding	4	1.64×10−5
GO:0005509	Calcium ion binding	25	3.86×10−5
GO:0045236	CXCR chemokine receptor binding	4	1.103×10−4
GO:0051087	Chaperone binding	7	1.165×10−4
GO:0051082	Unfolded protein binding	8	1.522×10−4
B, Downregulated
Go ID	Term	Count	P-value
BP
GO:0007259	JAK-STAT cascade	7	3.34×10−5
GO:0002573	Myeloid leukocyte differentiation	8	9.61×10−5
GO:0046427	Positive regulation of JAK-STAT cascade	5	1.33×10−4
GO:0046425	Regulation of JAK-STAT cascade	6	1.68×10−4
GO:0048856	Anatomical structure development	56	1.78×10−4
CC
GO:0036056	Filtration diaphragm	2	1.51×10−3
GO:0036057	Slit diaphragm	2	1.51×10−3
GO:0030054	Cell junction	18	1.55×10−3
GO:0032809	Neuronal cell body membrane	2	3.48×10−3
GO:0044298	Cell body membrane	2	3.49×10−3
MF
GO:0034046	Poly(G) binding	2	3.50×10−4
GO:0004402	Histone acetyltransferase activity	4	8.47×10−4
GO:0090595	Acetyl-CoA:L-lysine N6-acetyltransferase	4	8.47×10−4
GO:0000982	RNA polymerase II core promoter proximal region sequence-specific DNA binding transcription factor activity	8	1.14×10−3
GO:0000983	RNA polymerase II core promoter sequence-specific DNA binding transcription factor activity	2	3.10×10−3

Count represents the number of significantly enriched genes. GO, Gene Ontology; CXCR, chemokine (C-X-C motif) receptor; JAK-STAT, Janus kinase-signal transducer and activator of transcription; Acetyl-CoA, acetyl-coenzyme A; BP, biological process; CC, cellular component; MF, molecular function.

In addition, the result of KEGG pathway enrichment analysis was shown in Table II. The upregulated DEGs were mainly enriched in 15 pathways, including protein processing in endoplasmic reticulum, metabolic pathways and cytokine-cytokine receptor interaction. By contrast, the downregulated DEGs were mainly enriched in 9 pathways, including the Janus kinase-signal transducer and activator of transcription (JAK-STAT) signaling pathway and mitogen-activated protein kinase (MAPK) signaling pathway.

Table II.

KEGG pathway enrichment analysis results of differentially expressed genes.

A, Upregulated
KEGG ID	Term	Count	P-value
04141	Protein processing in endoplasmic reticulum	14	4.13×10−6
00250	Alanine, aspartate and glutamate metabolism	5	3.78×10−4
01100	Metabolic pathways	37	8.71×10−4
00520	Amino sugar and nucleotide sugar metabolism	5	2.50×10−3
00900	Terpenoid backbone biosynthesis	3	2.95×10−3
00670	One carbon pool by folate	3	5.06×10−3
03430	Mismatch repair	3	1.02×10−2
04610	Complement and coagulation cascades	5	1.18×10−2
00072	Synthesis and degradation of ketone bodies	2	1.29×10−2
04060	Cytokine-cytokine receptor interaction	11	1.61×10−2
00790	Folate biosynthesis	2	1.93×10−2
05150	Staphylococcus aureus infection	4	2.34×10−2
05020	Prion diseases	3	3.18×10−2
00051	Fructose and mannose metabolism	3	3.42×10−2
03030	DNA replication	3	3.42×10−2
B, Downregulated
KEGG ID	Term	Count	P-value
04630	JAK-STAT signaling pathway	6	4.88×10−4
04010	MAPK signaling pathway	6	7.85×10−3
04920	Adipocytokine signaling pathway	3	1.01×10−2
04512	ECM-receptor interaction	3	1.84×10−2
05219	Bladder cancer	2	3.13×10−2
04930	Type II diabetes mellitus	2	4.01×10−2
04670	Leukocyte transendothelial migration	3	4.11×10−2
04510	Focal adhesion	4	4.25×10−2
05213	Endometrial cancer	2	4.63×10−2

Count represents the number of significantly enriched genes. KEGG, Kyoto encyclopedia of genes and genomes; JAK-STAT, Janus kinase-signal transducer and activator of transcription; MAPK, mitogen-activated protein kinase; ECM, extracellular matrix.

PPI network and network module analyses

The constructed PPI network of DEGs is shown in Fig. 1A. In total, 173 nodes and 290 interacting protein pairs were contained in the PPI network. There were 12 DEGs with degree ≥10, such as early growth response 1 (EGR1; degree, 20), MYC (degree, 19), heat shock 70kDa protein 5 (HSPA5; degree, 16), interleukin 8 (IL8; degree, 15), and chemokine (C-X-C motif) receptor 4 (CXCR4; degree, 10).

Figure 1.

PPI network of DEGs. (A) The PPI network of DEGs in pancreatic cancer cells with glutamine. (B) The sub-network module 1 in PPI. (C) The sub-network module 2 in PPI. The red nodes represent upregulated DEGs and the green nodes represent downregulated DEGs. PPI, protein-protein interaction; DEGs, differentially expressed genes.

From the constructed PPI network, two sub-network modules were obtained. The genes in the two modules were all upregulated DEGs. In total, 10 DEGs and 38 interacting pairs were contained in module 1 (Fig. 1B), including IL8, CXCR4 and CXCR3. Additionally, 11 DEGs and 28 interacting pairs were contained in module 2 (Fig. 1C), including HSPA6 and HSPA5.

Literature mining of the network module

The co-occurrence network of module 1 is shown in Fig. 2A. In total, 8 genes were contained in the network. In addition, according to the enrichment score, the DEGs of module 1 were significantly enriched in 5 clusters and 1 single function (Table III). A heat map based on the genes and functions in Table IV was constructed (Fig. 2B).

Figure 2.

Literature mining results for network modules. (A) The co-occurrence network of DEGs in module 1 from literature mining. (B) The functional heat map of DEGs in module 1. (C) The co-occurrence network of DEGs in module 2 from literature mining. (D) The functional heat map of DEGs in module 2. PPI, protein-protein interaction; DEGs, differentially expressed genes.

Table III.

Functional enrichment results of differentially expressed genes in module 1.

A, Cluster1 (enrichment score, 13.92)
Keywords	P-value	Gene list
Chemokine receptor	7.36×10−15	C3; CCL20; CXCL1; CXCL14; CXCL2; CXCL3; CXCR3; CXCR4; IL8
Quantitative real time	5.29×10−14	C3; CCL20; CXCL1; CXCL14; CXCL2; CXCL3; CXCR3; CXCR4; IL8
Extracellular matrix	8.83×10−15	C3; CCL20; CXCL1; CXCL14; CXCL2; CXCL3; CXCR3; CXCR4; IL8
Cell migration	6.27×10−15	C3; CCL20; CXCL1; CXCL14; CXCL2; CXCL3; CXCR3; CXCR4; IL8
B, Cluster2 (enrichment score, 13.49)
Keywords	P-value	Gene list
Monocyte chemotactic protein	5.60×10−13	CCL20; CXCL1; CXCL2; CXCL3; CXCR3; CXCR4; IL8
Macrophage inflammatory protein	2.70×10−12	CCL20; CXCL1; CXCL2; CXCL3; CXCR3; CXCR4; IL8
Chemokine	2.26×10−17	CCL20; CXCL1; CXCL14; CXCL2; CXCL3; CXCR3; CXCR4; IL8
C, Cluster3 (enrichment score, 11.32)
Keywords	P-value	Gene list
Tumor necrosis factor	2.98×10−11	C3; CCL20; CXCL1; CXCL2; CXCL3; CXCR3; CXCR4; IL8
Monocyte chemoattractant protein	1.61×10−13	C3; CCL20; CXCL1; CXCL2; CXCL3; CXCR3; CXCR4; IL8
Innate immune response	2.35×10−11	C3; CCL20; CXCL1; CXCL2; CXCL3; CXCR3; CXCR4; IL8
D, Cluster4 (enrichment score, 10.42)
Keywords	P-value	Gene list
Innate immune system	2.34×10−11	C3; CCL20; CXCL1; CXCL2; CXCR3; CXCR4; IL8
Adaptive immune response	6.24×10−11	C3; CCL20; CXCL1; CXCL2; CXCR3; CXCR4; IL8
E, Single1 (enrichment score, 10.20)
Keywords	P-value	Gene list
Central nervous system	6.25×10−11	C3; CCL20; CXCL1; CXCL14; CXCL2; CXCR3; CXCR4; IL8; RGS4
F, Cluster5 (enrichment score, 10.11)
Keywords	P-value	Gene list
CC chemokine ligand	9.19×10−11	CCL20; CXCL2; CXCR3; CXCR4; IL8
Draining lymph node	6.59×10−11	CCL20; CXCL1; CXCL2; CXCR3; CXCR4; IL8

CXCR, chemokine (C-X-C motif) receptor; CCL, chemokine (C-C motif) ligand; IL, interleukin; RGS4, regulator of G protein signaling 4; C3, complement C3.

Table IV.

Functional enrichment results of differentially expressed genes in module 2.

A, Single1 (enrichment score, 12.14)
Keywords	P-value	Gene list
Heat shock protein	7.17×10−13	DNAJB1; DNAJB11; DNAJB4; DNAJB5; DNAJC12; DNAJC3; HSPA1A; HSPA5; HSPA6; PDIA4
B, Cluster1 (enrichment score, 11.25)
Keywords	P-value	Gene list
Unfolded protein response	1.47×10−9	DNAJB11; DNAJB9; DNAJC3; HSPA1A; HSPA5; PDIA4
Glucose regulated protein	2.18×10−14	DNAJB11; DNAJB9; DNAJC12; DNAJC3; HSPA1A; HSPA5; PDIA4
C, Cluster2 (enrichment score, 10.23)
Keywords	P-value	Gene list
Chaperones	1.20×10−9	DNAJB1; DNAJB11; HSPA5; HSPA6; PDIA4
Folding	2.34×10−9	DNAJB1; DNAJB11; DNAJC3; HSPA5; PDIA4
Endoplasmic reticulum stress	1.98×10−11	DNAJB1; DNAJB11; DNAJB9; DNAJC12; DNAJC3; HSPA5; PDIA4
Protein folding	2.18×10−13	DNAJB1; DNAJB11; DNAJB9; DNAJC12; DNAJC3; HSPA1A; HSPA5; PDIA4
Endoplasmic reticulum	6.08×10−11	DNAJB1; DNAJB11; DNAJB4; DNAJB9; DNAJC12; DNAJC3; HSPA1A; HSPA5; PDIA4
D, Single2 (enrichment score, 8.64)
Keywords	P-value	Gene list
ATPase activity	2.30×10−9	DNAJB1; DNAJB11; DNAJB5; DNAJC12; HSPA1A; HSPA5; HSPA6

DNAJ, DnaJ heat shock protein family; HSPA, heat shock protein family A; PDIA4, protein disulfide isomerase family a member 4.

Additionally, the co-occurrence network of module 2 was revealed in Fig. 2C. In total, 9 genes were included in the network. The DEGs of module 2 were significantly enriched in 2 clusters and 2 single functions (Table IV). The constructed heat map based on the genes and functions in Table IV was shown in Fig. 2D.

Discussion

In the present study, a total of 495 genes were identified as DEGs between the glutamine and control groups. These DEGs were mainly enriched in functions associated with response to organic substance, and metabolic pathway and JAK-STAT signaling pathway. Additionally, in the PPI network, MYC, HSPA5, IL18 and CXCR4 had high connectivitydegree. The majority of the DEGs were found to be hotspot genes based on literature mining.

In the PPI network, MYC had a high connectivity degree and was considered as a hub gene. MYC encodes a multifunctional, nuclear phosphoprotein that plays an important role in cell cycle progression, apoptosis and cellular transformation. The overexpression of MYC can promote cell transformation between G1 and S phase and lead to cell proliferation and formation of cancer (18). Studies have shown that knockdown of MYC results in inhibited growth of PC cells (19,20). Notably, MYC has been documented to induce the expression of mitochondrial glutaminase to stimulate glutamine catabolism, which plays an important role in cancer cell metabolism (21).

MYC can be regulated by the JAK-STAT signaling pathway, which was a significant pathway in the present study (22). The JAK-STAT signaling pathway participates in immune function and cell growth and differentiation (23). Additionally, components of the pathway, such as STAT3, have been shown to promote uncontrolled cell growth through dysregulation of gene expression involved in apoptosis, and cell-cycle regulation (24). As a result, it was hypothesized that glutaminase may have important roles in PC cell metabolism by regulating the JAK-STAT signaling pathway.

In particular, module analysis of the PPI network showed that two modules were obtained in the present study. By combining with literature mining, CXCR4 and IL8 were found to be key DEGs in module1. CXCR4 encodes the 7 trans-membrane G-protein-coupled receptor and a chemokine receptor specific for stromal cell-derived factor 1 (SDF1) (25). In cancer, CXCR4 is associated with metastasis to tissues that have a high concentration of SDF1 (26). The expression of CXCR4 has been suggested to play an important role in tumor cell invasion and metastasis in PC (27). In addition, IL8 encodes a chemokine that has pro-inflammatory effects (28). The association between inflammation and cancer has been well established. IL8 can also promote cancer stem-like cell invasion and metastasis, as well as treatment resistance (29). Therefore, glutaminase may increase the expression of CXCR4 and IL8 to promote the invasion and metastasis of PC cells.

In addition, in module 2, HSPA5 was a key DEG. HSPA5 is a regulator of endoplasmic reticulum (ER) function (30). The expression of HSPA5 is induced by ER stress and its overexpression has been reported in numerous types of cancer cells (31). Studies have shown that HSPA5 can inhibit the etoposide-mediated apoptosis by inhibiting activation of caspase-7 in cancer cells (32). Additionally, HSPA5 contributes to the growth of tumor and can induce drug resistance of cancer cells (31). Therefore, the expression of HSPA5 plays an important role in the progression of PC cells.

In conclusion, analysis of the gene expression profiles, significant differences in gene expression were found between glutamine and control group. Through analysis of DEGs, it was found that MYC, IL18, CXCR4 and HSPA5 may exert important roles in molecular mechanisms of PC cells with glutamine. However, additional experiments with larger samples are required to verify the present results.