Authors: Yihui Yang, Liwen Ren, Sha Li, Xiangjin Zheng, Jinyi Liu, Wan Li, Weiqi Fu, Jinhua Wang, Guanhua Du
Categories: Original Article, immune infiltration, pancreatic cancer, tumor microenvironment
Source: Translational Cancer Research
Doi: 10.21037/tcr-21-2021
Pancreatic cancer is one of the most commonly diagnosed and lethal malignancies worldwide and has few good biomarkers and therapeutic targets. GABRP is the π subunit of the gamma-aminobutyric acid (GABA) A receptor, which is expressed in a number of non-neuronal tissues. GABRP is significantly upregulated in pancreatic cancer, but its biological and immunological role as well as its clinical diagnostic and prognostic value in pancreatic cancer is still incompletely known.
In this study, pancreatic adenocarcinoma (PAAD) cohorts from TCGA and GEO datasets were used to compare GABRP mRNA levels in cancerous and normal tissues and protein expression was evaluated using immunohistochemistry. The Kaplan-Meier plotter and GEPIA2 database were used to analyze the correlation between GABRP expression, overall survival, and disease-free survival in pancreatic cancer patients. Gene set enrichment analysis (GSEA) was performed with the Linked Omics database to explore the molecular mechanisms of GABRP in pancreatic cancer. And the correlation between GABRP expression and immune infiltration was explored using the TIMER database, CIBERSORT database and ESTIMATE algorithm.
GABRP mRNA was significantly overexpressed in TCGA-PAAD cohorts (P<0.0001) and enhanced GABRP expression predicted poorer overall survival according to Kaplan-Meier plotter database (P=0.0024) and GEPIA2 (P=0.038). Hypomethylation of promoter (P<0.01) and the regulation of hsa-miR-3655 may contribute to the overexpression of GABRP in pancreatic cancer. GSEA analysis revealed that GABRP played an important role in the immune response. GABRP expression was also correlated with immune infiltration and immune cell markers. Higher GABRP expression was significantly associated with greater infiltration of immune cells and stromal cells into pancreatic cancer microenvironments as well as higher expression of six important immune check point genes including PDCD1 (P<0.05), CD274 (P<0.05),* CTLA4* (P<0.01), PDCD1LG2 (P<0.01), TIGHT (P<0.01) and TIM3 (P<0.01).
GABRP is a potential prognostic biomarker and is correlated with immune infiltration and tumor microenvironment in pancreatic cancer. This suggests that GABRP may serve as a potential prognostic biomarker and therapeutic target in pancreatic cancer as well as a possible regulator of tumor microenvironment affecting the efficacy of immunotherapy. Further studies are needed to elucidate the molecular mechanism of the immunoregulatory role of GABRP.
**Keywords: ** GABRP, pancreatic cancer, immune infiltration, tumor microenvironment
Pancreatic cancer, primarily pancreatic ductal adenocarcinoma (PDAC), is one of the most commonly diagnosed malignancies which ranks 14^th^ among all cancers and the 7^th^ leading cause of cancer-related death worldwide. Global cancer statistics reveal that there were 495,773 new cases and 466,003 deaths of pancreatic cancer in 2020 (1). Pancreatic cancer is a fatal malignancy with a 5-year survival of about 5% when diagnosed and a median survival of <11 months (2). Surgical resection is the only potential curative treatment for pancreatic cancer at present and adjuvant chemotherapies have improved patient outcomes, but at the cost of a higher incidence of adverse effect (3,4). However, just a minority of pancreatic cancer patients are eligible for surgical therapy since the disease lacks early symptoms until developed into advanced stages (5). In recent years, targeted therapy and immunotherapy have brought therapeutic benefit in many solid tumors. Several targeted therapies such as* VEGF* inhibitors and KRAS inhibitors have been evaluated in pancreatic cancer, but received minimal effect. The same can be seen in immunotherapy, which shows promising efficacy in several solid tumor especially melanoma including PD-1 inhibitors and CTLA-4 inhibitors while the response rate is unsatisfactory in pancreatic cancer (6-8). Thus, it remains imperative to identify potential prognostic biomarkers and therapeutic targets to facilitate the clinical treatment and discovery of clinical candidate compounds and to exploit possible strategies of immunotherapy for pancreatic cancer.
GABRP is the π subunit of the gamma-aminobutyric acid (GABA) A receptor, which is a hetero-pentameric and ligand-gated chloride channel composed of five subunits. These subunits are encoded by 19 genes in 8 subclasses based on sequence homology including α1–6, β1–3, γ1–3, δ, ε, θ, π and ρ1–3. GABAA receptors mainly distributed in the central nervous system which frequently consists of at least two α, two β subunits and one γ subunit. The γ subunit can be replaced by a δ, ε, π subunit, but the function of the π subunit is still unclear (9-11). Unlike other *GABAA *receptors, GABRP is detectable in multiple non-neuronal normal human tissues including uterus, prostate, ovaries, small intestine and lung. It is reported that GABRP may play an important role in endometrial receptivity and control receptor sensitivity to endogenous steroids and pregnenolone (12-14).
Although the function of GABRP in normal tissues has not been fully elucidated, several studies have revealed dysregulation of GABRP expression in solid tumors including breast, ovarian, cervical and pancreatic cancer (15-18). Jiang et al. proposed that GABRP promoted pancreatic tumorigenesis via KCNN4-mediated Ca^2+^ signaling to regulate chemokine and macrophage recruitment, which indicated the significance and potential immunomodulatory role of GABRP in pancreatic cancer (18). However, there has been no study to comprehensively explore the expression, prognostic value and regulatory role of GABRP in pancreatic cancer and whether GABRP was associated with immune infiltration or could affect the tumor microenvironment. Here, we integrated multiple ‘omics’ data about pancreatic cancer from multiple authoritative databases and used bioinformatics analysis to comprehensively and systematically analyze the expression profile, clinicopathological characteristics and potential molecular mechanism of GABRP in pancreatic cancer. For the first time, we elucidated the potential mechanism of GABRP dysregulation and showed how GABRP might participate in regulating pancreatic tumor infiltrating cells and tumor microenvironment, which further confirmed the important role of GABRP as a prognostic biomarker and therapeutic target in pancreatic cancer as well as a possible regulator of tumor microenvironment affecting the efficacy of immunotherapy (Figure 1A). We present the following article in accordance with the REMARK reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-21-2021/rc).
Figure 1 GABRP mRNA expression was upregulated in pancreatic cancer. (A) The flowchart of study procedures; (B) the expression level of GABRP was measured in different types of tumor tissues and compared with normal tissues in the GEPIA2 database; (C) the expression level of GABRP was measured in different types of tumor tissues and normal tissues using the Oncomine database (P value is 0.001, fold change is 2, and gene rankings apply to all); (D) GABRP expression in tumors and matching normal tissues using two independent cohorts (GSE15471, n=78; GSE16515, n=52) derived from Gene Expression Omnibus datasets; (E) immunohistochemical analysis of GABRP expression in a human PDAC tissue microarray. Representative GABRP images are shown in the left panel. The percentages of tissues displaying low and high staining in normal pancreas and PDAC tissues are shown in the middle panel. ****, P<0.0001. TCGA, The Cancer Genome Atlas; GEO, Gene expression omnibus; DEGs, differential expressed genes; IHC, immunohistochemistry; ACC, adrenocortical carcinoma; BLCA, bladder urothelial carcinoma; BRCA, breast invasive carcinoma; CESC, cervical squamous cell carcinoma and endocervical adenocarcinoma; CHOL, cholangiocarcinoma; COAD, colon adenocarcinoma; DLBC, lymphoid neoplasm diffuse large B-cell lymphoma; ESCA, esophageal carcinoma; GBM, glioblastoma multiforme; HNSC, head and neck squamous cell carcinoma; KICH, kidney chromophobe; KIRC, kidney renal clear carcinoma; KIRP, kidney renal papillary cell carcinoma; LAML, acute myeloid leukemia; LGG, brain lower grade glioma; LIHC, liver hepatocellular carcinoma; LUAD, lung adenocarcinoma; LUSC, lung squamous cell carcinoma; MESO, mesothelioma; OV, ovarian serous cystadenocarcinoma; PAAD, pancreatic adenocarcinoma; PCPG, pheochromocytoma and paraganglioma; PRAD, prostate adenocarcinoma; READ, rectum adenocarcinoma; SARC, sarcoma; SKCM, skin cutaneous melanoma; STAD, stomach adenocarcinoma; STES, stomach and esophageal carcinoma; TGCT, testicular germ cell tumors; THCA, thyroid carcinoma; THYM, thymoma; UCEC, uterine corpus endometrial carcinoma; UCS, uterine carcinosarcoma; UVM, uveal melanoma; NP, normal pancreas; PDAC, pancreatic ductal adenocarcinoma.
The data on GABRP mRNA expression in pancreatic adenocarcinoma (PAAD) tissues and normal tissues and corresponding clinical data were retrieved from the TCGA database (https://portal.gdc.cancer.gov/). A total of 178 tumor tissues and 4 normal tissues obtained from The Cancer Genome Atlas (TCGA) database were included in the data analysis. The GABRP mRNA expression profiles of two pancreatic cancer datasets from the gene expression omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) were also included in the data analysis, which consisted of microarrays with accession number GSE15471 (GEO platform GPL570) which contained 39 tumor samples and 39 normal samples and the GSE16515 (GEO platform GPL570) which contained 36 tumor samples and 16 normal samples (19,20). The GABRP mRNA expression profiles of three pancreatitis datasets from the GEO database were also included in the data analysis including microarrays with the accession number, GSE15471 (GEO platform GPL6246), which contained six pancreatitis samples and five normal samples, the GSE65146 (GEO platform GPL6246), which contained 64 pancreatitis samples and 9 normal samples, and the GSE41418 (GEO platform GPL1261) which contained 6 pancreatitis samples and 6 normal samples.
The Oncomine database (http://www.oncomine.org) incorporated 86,733 samples and 715 gene expression datasets into a single comprehensive database allowing various analytical methods (21). We used the Oncomine database to evaluate the mRNA expression level of the GABRP gene among various cancer types. The threshold significance was determined by a P value of 0.001, fold-change of 2, gene rank of ‘all’ and data type of ‘mRNA’.
The gene expression profiling interactive analysis 2 (GEPIA2) database (http://GEPIA2.cancer-pku.cn/index.html) was utilized to evaluate the expression of GABRP in different cancer types. Featuring 198,619 isoforms and 84 cancer subtypes, GEPIA2 (22) is an updated and enhanced version of GEPIA, a web resource for gene expression analysis based on tumor and normal samples from the TCGA and the GTEx databases. The gene expression profiles of GABRP among all tumor samples paired normal tissues were shown in a dot-plot.
Two datasets, GSE15471 and GSE16515, were used to analyze the expression of GABRP in pancreatic cancer between normal tissues and tumor tissues. The data was processed by R software v4.0.3 and the plot was constructed using the R software package, ggplot2.
The Kaplan-Meier plotter (www.kmplot.com) is a web-based tool for exploring the impact of 54,000 genes on patient survival covering 21 cancer types (23). The database is based on gene expression data, relapse-free and overall survival information downloaded from TCGA, GEO and the European Genome-Phenome Archive (EGA) database. The Kaplan-Meier plotter was used to perform prognostic analysis of GABRP and hsa-miR-3655 expression levels in PAAD, and the association between GABRP expression and survival outcome based on different immune cell contents with auto-selected cutoff values. GEPIA2 was used to perform prognostic analysis of GABRP expression levels in PAAD based on GABRP expression with a specific cutoff value using the R software package ‘survival’. The log rank P value, the hazard ratio (HR), and survival curves were calculated and annotated for each plot.
Univariate and multivariate Cox regression analysis was performed using the prognosis model from the Assistant for Clinical Bioinformatics (https://www.aclbi.com/static/index.html#/prognosis) to identify the appropriate terms to build the nomogram. The Forest was used to show the P value, HR and 95% CI of each variable through ‘forestplot’ in the R package. The variable that shows significance in multivariate Cox proportional analysis shows that it should be an independent indicator of prognosis. Based on the results of multivariate Cox proportional hazards analysis, a nomogram was developed to predict the X-year overall recurrence through the ‘rms’ R package.
The cBioPortal for cancer genomics was used to explore the relationship between GABRP DNA methylation and GABRP expression in TCGA-PAAD (24). TargetScan, miRDB and miRwalk were used to predict the regulatory miRNAs of GABRP (25-27). The common miRNAs predicted by the above three databases were included in the prediction cohort. The common miRNAs significantly downregulated in PAAD in the three cohorts were then regarded as potential regulatory miRNAs of PARPBP in HCC. UALCAN (28) (http://ualcan.path.uab.edu/index.html) is an interactive web-portal web resource to perform comprehensive analysis of TCGA gene expression data. We used UALCAN to assess the methylation levels of *GABRP *and the expression of miR-3655 in PAAD and its effect on prognosis.
LinkedOmics is a publicly available portal that includes multi-omics data from all 32 TCGA cancer types and provides a unique platform for molecular biologists and clinicians to access, analyze and compare multi-omics cancer data within and across tumor types (29). We searched for GABRP co-expressed genes in the PAAD dataset (n=178) based on the Pearson correlation coefficient. In the LinkInterpreter module, gene set enrichment analysis (GSEA) module was used to assess the biological processes (BPs), cell components (CCs), molecular function (MFs), and KEGG pathways based on co-expressed genes. The GSEA program was run with 500 simulations. If the P value and false discovery rate (FDR) were both less than 0.05, the gene set was considered significantly enriched.
The tumor immune estimation resource (TIMER, https://cistrome.shinyapps.io/timer/) is a web server designed for comprehensive analysis of tumor-infiltrating immune cells in 32 cancer types from TCGA (30). The TIMER database was used to analyze the relationship between the level of GABRP expression and six tumor-infiltrating immune cells (TIICs) and tumor purity in PAAD. The six TIICs include B cells, CD4^+ ^T cells, CD8^+ ^T cells, macrophages, neutrophils, and dendritic cells. The SCNA (somatic copy number alteration) module was used to compare tumor infiltration levels in PAAD with different somatic copy number alterations for GABRP. Then the correlation between GABRP and several immune cell markers were assessed by using the correlation module (31).
Cell type identification by estimating relative subsets of RNA transcripts (CIBERSORT, https://cibersort.stanford.edu/) is a deconvolution algorithm developed to characterize 22 types of immune cells composition in complex tissues (32). By integrating TCGA-PAAD datasets, GSE15471 and GSE16515 datasets, we analyzed the integrated mRNA expression matrix using the CIBERSORT R script acquired from the CIBERSORT website. After filtering samples with P≥0.05, 59 normal samples and 205 tumor samples were included in the study.
The ESTIMATE algorithm (33) uses gene expression signatures to estimate the fraction of stromal cells and immune cells and calculates the immune score and stromal score of PAAD tumor samples. The algorithm was implemented using the ‘estimate’ package by R software. Two subtypes of PAAD were classified based on GABRP expression (high vs. low) and compared using the Kruskal-Wallis test.
The PAAD chip (Cat. No. HPanA060CS03-M-109) was purchased from SuperChip in Shanghai, China. The chip contained two normal pancreatic tissues and 29 PAAD tumor and adjacent tumor tissues, which included phases I (3 cases), II (23 cases) and III (3 cases) PAAD tumor tissues. After routinely dewaxing and hydration, antigen retrieval was performed in 10 mM saline sodium citrate at sub-boiling temperature for 90 s. Then, the activity of endogenous peroxidases in specimens was blocked with 3% hydrogen peroxide. After permeabilization with 0.1% Triton X-100 and blocking with 5% bovine serum albumin (BSA), the chip was incubated with primary antibody against rabbit anti-GABRP antibody (1:200 dilution, PA5-46830, Invitrogen) at 4 ℃ overnight. The chip was balanced at room temperature for 30 min, washed with PBS and incubated with horseradish peroxidase-conjugated goat anti-mouse IgG (DAKO, Wuhan, China) for 60 min at room temperature. GABRP expression was visualized by DAB (DAKO, Wuhan, China) staining. The chip was scanned by Pannoramic Scan (3DHISTECH Ltd., Hungary).
The Wilcoxon rank sum test was employed to explore the different expression pattern between PAAD tumor samples and normal samples. The PAAD cohorts were divided into two groups using the optimal threshold evaluated by X-tile 3.6.1 on the basis of GABRP expression levels (34). The relationship between GABRP expression and clinicopathological parameters was analyzed using Chi-squared test. The correlation between GABRP and TIICs was calculated by Spearman correlation. The Cox proportional hazards model was used to identify the independent indicators related to OS, and the hazard ratios (HRs) with 95% confidence intervals (CI) were also calculated. Nomograms were constructed based on the results of multivariate Cox analysis on the basis of OS to predict the X-year overall recurrence using the ‘rms’ R package. All statistical analyses were performed using R software (Version 4.0.2), and the threshold of statistical significance was set at P<0.05.
The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). Institutional ethical approval and informed consent were waived.
The mRNA expression of GABRP in different cancers and normal tissues was evaluated by analyzing the GEPIA2 and Oncomine databases (Figure 1). In GEPIA2, the expression of GABRP was obviously elevated in PAAD compared to normal tissues (Figure 1B). In the Oncomine database, the expression of GABRP was significantly higher in PAAD than normal tissues (Figure 1C). Next, the expression of GABRP in pancreatic cancer was further validated by analyzing two GEO cohorts (GSE15471 and GSE16515) (19,20,35,36). In GEO cohort, GSE15471 and GSE16515 (Figure 1D), GABRP was significantly overexpressed in tumor samples compared to normal pancreas tissues. To further check the expression pattern of GABRP protein in PDAC, immunohistochemical (IHC) analysis was carried out in tissue microarrays of human PDAC. It was shown in Figure 1E that GABRP expression is higher in PDAC than in normal pancreas, while commonly detected in the PDAC cells (Figure 1E). However, there was no significant change in the expression of GABRP in pancreatitis (Table S1). In summary, GABRP was highly expressed in pancreatic cancer.
The clinicopathological data of 178 PAAD patients were obtained from the TCGA database and separated into low or high GABRP expression based on the optimal threshold of OS. As shown in Table 1, there was no significance in age, gender, TNM stage and living status while grade was significantly associated with the expression level of GABRP (P<0.05). The correlation between GABRP expression and survival outcome of PAAD was next investigated using the Kaplan-Meier plotter. Results from Kaplan-Meier (KM) survival analysis showed that elevated expression of GABRP was linked to poor prognosis in PAAD. Patients with higher GABRP expression had shorter overall survival than those with lower GABRP expression (Figure 2A, HR =1.91, P=0.0024). Similarly, statistically significant lower disease-free survival (DFS) was observed in patients with higher GABRP expression (Figure 2B, HR =5, P=0.00073). GEPIA2 was used to confirm the effect of GABRP expression on pancreatic cancer prognosis (Figure 2C,2D). Similar results from GEPIA2 disclosed that patients with high GABRP expression had worse prognosis in TCGA PAAD cohorts.
Figure 2 Overexpressed GABRP independently predicts poor OS in PAAD. (A,B) Correlation between GABRP and prognosis of pancreatic cancer in Kaplan-Meier plotter database; (C,D) correlation between *GABRP *and prognosis of pancreatic cancer in GEPIA2 database. OS, overall survival; PAAD, pancreatic adenocarcinoma.
The prognostic value of GABRP in PAAD was first assessed by using the Kaplan-Meier plotter and GEPIA2. In the next step, Cox regression analysis was performed to explore independent risk factors of OS by using data from the TCGA PAAD cohort. In the univariate model, GABRP expression, age and grade were significantly related to OS in PAAD (Figure 3A, all P<0.05). Multivariate analysis further confirmed that GABRP expression was an independent risk factor of OS in PAAD [Figure 3B, HR =1.09879 (1.00705, 1.19888), P=0.03419]. Finally, a prognostic nomogram for PAAD was constructed to further predict the prognostic value of GABRP in PAAD based on mRNA expression and GABRP expression was identified as an independent indicator of OS. As shown in Figure 3C, a nomogram was generated to assess the performance of the risk score in combination with the GABRP expression of PAAD patients in predicting the OS. The concordance index (c-index) value was 0.541 for PAAD, which measures the level of agreement between prediction by nomogram and the observed survival status. The calibration curves for 1-, 2-, 3- and 5-year survival were also plotted for PAAD (Figure 3D).
Figure 3 Overexpressed GABRP independently predicts poor OS in PAAD. (A) Univariate Cox proportional hazards regression model analysis based on overall survival in PAAD; (B) multivariate Cox proportional hazards regression model analysis based on overall survival in PAAD; (C) postoperative prognostic nomogram for patients with PAAD; (D) the calibration curve of the nomogram for predicting overall survival at 1, 2, 3 and 5 years. Observed OS was plotted on the y-axis and nomogram predicted probability of OS was plotted on the x-axis. OS, overall survival; PAAD, pancreatic adenocarcinoma.
The potential mechanism of GABRP overexpression in PAAD was explored from the perspective of epigenetic regulation. Increasingly, studies have revealed abnormal DNA methylation in cancers and many non-neoplasm diseases. Promoters with cancer-associated hypomethylation are always linked to the activation of genes that were silent in normal tissues (37). The analysis via cBioportal revealed a significant negative correlation between *GABRP *methylation and GABRP expression (Figure 4A, Spearman’s r=−0.77, P=6.52e-37, Pearson’s r=−0.66, P=1.40e-23). In addition, UALCAN was used to verify the promoter methylation of GABRP in PAAD, and the level of *GABRP *promoter methylation in PAAD was significantly lower than in normal tissues (Figure 4B). Patients with PAAD were stratified based on the stage, grade, nodal metastasis status and TP53 mutant status. The results showed that lower levels of GABRP promoter methylation were associated with higher-grade tumors, advanced individual cancer stages, lymph node metastasis and TP53 mutant status (Figure S1A-S1D).
Figure 4 Promoter hypomethylation and miR-3655 downregulation contributed to GABRP upregulation in PAAD (A) Correlation analysis between GABRP expression and GABRP methylation; (B) comparison of GABRP promoter methylation level in tumor tissues compared with normal tissues; (C) common miRNAs containing GABRP-binding sequences using miRDB, TargetScan and MirWalk databases; (D) expression of miR-3655 in PAAD tissues compared with normal tissues; (E) low miR-3655 predicts worse survival in PAAD; (F) putative binding site of GABRP 3'-UTR by miR-3655. **, P<0.01; ****, P<0.0001. PAAD, pancreatic adenocarcinoma.
We next explored the potential regulatory miRNAs of GABRP in PAAD using miRDB, TargetScan and MirWalk. Seventeen common miRNAs were selected as candidates containing GABRP-binding sequences (Figure 4; Table S2). After evaluation with UALCAN, hsa-miR-3655 was identified as a potential regulatory miRNA of GABRP. Hsa-miR-3655 was significantly downregulated in PAAD compared with normal tissues in the TCGA-PAAD dataset (Figure 4D, P=0.0067). KM analysis revealed that low expression of hsa-miR-3655 was associated with worse survival (Figure 4E, P=0.015).
LinkedOmics was used to explore the co-expressed genes of GABRP in PAAD which should shed light on the biological meaning of GABRP in the process of tumorigenesis and development. A total of 6,647 genes had co-occurrences in PAAD (false discovery rate, FDR <0.01). A total of 4,100 genes (red dots) were significantly positively correlated with GABRP expression and 2,547 genes (green dots) were significantly negatively correlated with GABRP expression in PAAD (Figure S2A). The top 50 ranked genes significantly correlated with GABRP expression were visualized in a heatmap (Figure S2B,S2C). GABRP expression showed a strong negative association with expression of PEMT (r=−0.5772, P=3.34E-17), ATP6AP1 (r=−0.5402, P=7.08E-15), HDAC5 (r=−0.5217, P=8.17E-14), etc. (Figure S2B). GABRP expression also showed a strong positive association with expression of GTB4 (r=0.6364, P=1.33E-21), MMP7 (r=0.6193, P=3.14E-20), LOC284578 (r=0.6085, P=2.06E-19), etc. (Figure S2C). GSEA of co-expressed genes of GABRP was performed to explore the underlying biological functions of GABRP overexpression in PAAD. Gene ontology (GO) term annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis by GSEA are presented in Figure 5. In biological processes (BPs), GABRP co-expressed genes were correlated with regulation of leukocyte activation, immune response-regulating signaling pathway, adaptive immune response, etc. (Figure 5A). In cellular components (CCs), GABRP co-expressed genes were related to immunological synapse, cell substrate junction, MHC protein complex, etc. (Figure 5B). In molecular functions (MFs), GABRP co-expressed genes were associated with collagen binding, cytokine receptor activity, cytokine binding, etc. (Figure 5C). KEGG pathway analysis showed enrichment in the NOD-like receptor signaling pathway, NF-kappa B signaling pathway, cytokine-cytokine receptor interaction, etc. (Figure 5D). The GO analysis and KEGG analysis results revealed that GABRP co-expressed genes had a widespread impact on immune responses.
Figure 5 Analysis of genes co-expressed with GABRP in PAAD. The bubble plot for (A) GO-BP, (B) GO-CC and (C) GO-MF and (D) KEGG pathway analysis by GSEA. NES, normalized enrichment score; PAAD, pancreatic adenocarcinoma; GO, Gene ontology; BP, biological process; CC, cellular component; MF, molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes; GSEA, gene set enrichment analysis; FDR, false discovery rate.
We systematically depicted the pattern of immune cells in pancreatic cancer by extracting and processing the signature gene expression profiles from TCGA-PAAD and GEO cohorts (GSE15471 and GSE16515) performed by CIBERSORT. After filtering the samples for P≥0.05, 59 normal samples and 205 tumor samples were included in the analysis. The differential immune infiltration landscape of pancreatic cancer and normal tissues was displayed in Figure 6. There was a different pattern of infiltrating immune cells between normal tissues and pancreatic cancer tissues (Figure 6A). Moreover, different correlation patterns among 22 immune cells were found in pancreatic cancer and normal tissues (Figure 6B). In the study cohort, the fractions of 22 immune cells were correlated with each other in different degrees suggesting a complex tumor infiltrating environment in pancreatic cancer. The violin plot of 22 immune cells showed that the proportion of different infiltrating immune cells significantly varied between tumor and normal tissues except for macrophages M1 and M2 (Figure 6C). Compared to normal tissues, the fractions of naïve B cells, CD4^+^ T cells and macrophage M0 were relatively higher in tumor tissues while CD8^+^ T cells and NK cells were significantly lower in tumor tissues.
Figure 6 Differential immune infiltration landscape of pancreatic cancer and normal tissues. (A) Heatmap and (C) violin plot indicate differences in immune cell distribution between tumor (red) and normal (blue) tissues in PAAD. (B) Correlation patterns among 22 immune cell subsets in PAAD. PAAD, pancreatic adenocarcinoma.
Based on the results of GSEA analysis, we first assessed the relationship between GABRP and six types of infiltrating immune cells including B cells, CD8^+^T cells, CD4^+^T cells, macrophages, neutrophils and dendritic cells by accessing the TIMER database (Figure 7A). The results suggested that the expression level of GABRP was positively correlated with the infiltration degree of B cells (r=0.22, P=3.86E-03), CD8^+^ T cells (r=0.298, P=7.73E-05), macrophages (r=0.194, P=1.10E-02), neutrophils (r=0.277, P=2.41E-04) and dendritic cells (r=0.264, P=4.94E-04) in PAAD. In addition, a higher GABRP expression level was associated with lower tumor purity (r=−0.295, P=8.63E-05), which may contribute to worse clinical features and outcomes (38). Genes with somatic copy number alterations (sCNAs) could influence tumorigenesis, progression and immunotherapy response (39). GABRP copy number alterations were related to immune cell infiltration by B cells, CD4^+^ T cells, macrophages, and neutrophils in PAAD (Figure 7B).
Figure 7 Correlation between GABRP expression and immune infiltration in pancreatic cancer. (A) Association between *GABRP *expression and immune cell infiltration in PAAD; (B) associations between GABRP gene copy numbers and immune cell infiltration in PAAD. *, P<0.05; **, P<0.01; ***, P<0.001. PAAD, pancreatic adenocarcinoma.
In order to clarify the reason why GABRP was associated with immune infiltration, we further explored the correlations between GABRP expression and a series of immune markers related to various types of immune cells using the TIMER database. Several immune cells were characterized by immune markers including B-cells, CD8^+^ T cells, neutrophils, macrophages, dendritic cells, Th1 cells, Th2 cells and Treg cells, and their correlation with GABRP expression was determined using the TIMER database (40,41). As shown in Table 2, GABRP expression had a significantly positive correlation with immune makers of B cells, CD8^+^ T cells, neutrophils, dendritic cells, Th1 cells, Th2 cells and Treg cells, without regard to age or tumor purity. The macrophage markers were significantly correlated with GABRP expression, except for CD84 and CD163 in the condition adjusted by tumor purity optionally. The above analysis laid the foundation for the important role of GABRP in immune infiltration.
Given the effect of GABRP on prognosis and immune cell infiltration in PAAD, we determined whether GABRP expression could affect the prognosis of pancreatic cancer due to immune cell infiltration. We performed a prognostic analysis of GABRP based on different immune cell contents in PAAD using the Kaplan-Meier plotter. The results indicated that higher GABRP expression related to worse prognosis in terms of enriched B cells, decreased CD4^+^ T cells, enriched CD8^+^ T cells, decreased macrophages, enriched NK cells, and decreased Treg cells compared to decreased B cells, enriched CD4^+^ T cells, decreased CD8^+^ T cells, enriched macrophages, decreased NK cells, and enriched Treg cells (Figure 8). In addition, the different profiles of CD4^+^ T cells and macrophages could considerably alter the effect of GABRP expression on survival outcome of PAAD (Figure 8C,8D,8G,8H). The above results suggested that higher expression of GABRP may affect survival of pancreatic cancer partly due to immune infiltration.
Figure 8 Prognostic analysis of GABRP based on different immune cell levels in pancreatic cancer. (A-L) Relationships between GABRP of different immune cell subgroups and PAAD prognosis. PAAD, pancreatic adenocarcinoma.
As one of the most lethal solid tumors worldwide, pancreatic cancer is characterized by an immunosuppressive and obstructive tumor microenvironment as well as a highly malignant phenotype (42). Thus, gaining insight into the tumor environment is critical to improving immunotherapy. Given the impact of GABRP expression on immune infiltration in pancreatic cancer, we speculated whether the expression level of GABRP was associated with the status of the tumor microenvironment. ESTIMATE was used to evaluate the composition and infiltration of immune cells and stromal cells in the tumor microenvironment using gene expression signatures based on GABRP expression. Compared to the cohorts with lower GABRP expression, higher GABRP expression was significantly correlated with higher immune score, stromal score and lower tumor purity (Figure 9A-9D). Furthermore, we examined the expression of six immune checkpoint genes including PDCD1,* CD274*,* CTLA4*,* PDCD1LG2*,* TIGHT* and TIM3), which participated in immune escape and affected immunotherapy response. The results showed that expression of the above immune checkpoint genes was significantly upregulated in the cohort with higher GABRP expression (Figure 9E-9J, P<0.05).
Figure 9 Evaluation of tumor microenvironment in pancreatic cancer based on* GABRP* expression. (A) Heatmap of two different immune subtypes based on GABRP expression; (B-D) evaluation of immune scores, tumor purity and stromal scores for the two immune subtypes. Correlation between expression of (E) PDCD1, (F) CD274, (G)* CTLA4*, (H) PDCD1LG2, (I) TIGHT and (G) TIM3 evaluated by Mann-Whitney U-test; bars indicate medians. *, P<0.05; **, P<0.01.
Pancreatic cancer is one of the most fatal malignancies worldwide with the worst survival of all carcinomas (43). In recent years, the incorporation of immunotherapy has greatly revitalized and achieved great therapeutic advances in the treatment of many solid tumors (44). However, pancreatic cancer remains almost completely resistant to immunotherapy and the majority of clinical trials related to immunotherapy failed to prove clinical efficacy in pancreatic cancer (7). Thus, it is worthwhile to identify novel prognostic biomarker and therapeutic targets as well as to understand the factors affecting the efficacy of immunotherapy for pancreatic cancer. In this study, we performed a systematic bioinformatics analysis of GABRP in pancreatic cancer integrating the mRNA expression profiles from TCGA-PAAD and GEO cohorts. In pancreatic cancer, enhanced GABRP expression was correlated with a worse prognosis and was an independent predictor of survival. We found that the expression level of GABRP was associated with the number of infiltrating immune cells and several immune cell markers indicating the role of GABRP in immunological interactions. We found a potential link between GABRP expression and the tumor microenvironment, highlighting the need for further research to elucidate the immunological role of GABRP in pancreatic cancer.
We assessed GABRP expression in pancreatic cancer compared with normal tissues and found that GABRP was significantly upregulated in pancreatic cancer using the independent Oncomine and GEPIA2 databases, which was further validated in two GEO datasets (GSE15471 and GSE16515). We observed a remarkable correlation between GABRP and relatively poor prognosis in pancreatic cancer related to OS and DFS. Multivariate analysis revealed that the expression level of GABRP was an independent risk factor for unfavorable survival in pancreatic cancer. However, GABRP expression was not correlated with TNM stage, but grade, in the TCGA cohort, which may be due to the limited number of cases with metastasis. To further validate the prognostic value of GABRP in pancreatic cancer, nomograms were constructed based on GABRP expression. Even if the nomograms reached a moderate prognostic accuracy (C index =0.541), they may not be precise enough due to the lack of enough independent risk factors identified in multivariate Cox analysis. However, the above results suggested that GABRP holds promise as a valuable prognostic biomarker for pancreatic cancer to improve clinical decisions.
We further explored the potential mechanism of GABRP overexpression in pancreatic cancer from the perspective of epigenetics. DNA methylation plays an important role in regulating gene expression and aberrant DNA methylation is involved in cancer progression (45). DNA hypomethylation is frequently observed in human tumors and contributes to the aberrant activation of specific genes related to cancer progression (46). In pancreatic cancer, there exists a significantly negative correlation between GABRP methylation and *GABRP *mRNA expression. Moreover, methylation of the GABRP promoter was significantly upregulated in pancreatic tumor tissues. The above analysis implied that the enhanced expression of GABRP might be partially attributed to DNA hypomethylation. Micro-RNAs (miRNAs) can regulate the expression of target genes in many processes and have been identified as significant participants in the progression of pancreatic cancer (47). In the present study, hsa-miR-3655 was identified as a candidate regulatory miRNA for GABRP, and we showed that the expression of hsa-miR-3655 was significantly downregulated in pancreatic cancer compared with normal controls, and its downregulation was related to poor prognosis. The putative binding between the GABRP 3’-UTR site and miR-3655 further confirmed that miR-3655 could be a potential upstream regulator of *GABRP *in pancreatic cancer. The above findings suggested that epigenetic regulation might contribute to the dysregulation of GABRP.
To further clarify the underlying biological function of GABRP in pancreatic cancer, GSEA was performed to determine GABRP co-expressed genes. GO analysis showed consistent enrichment of genes in several immunologically relevant biological processes and pathways, which indicated that GABRP might participate in immune responses. An inflammatory microenvironment is now recognized to be an integral factor contributing to tumorigenesis (48). KEGG pathway analysis showed that several inflammation-related pathways were enriched such as the NOD-like receptor pathway, NF kappa B pathway, and TNF pathway, all of which are crucial inflammatory regulators involved in tumor initiation and progression (49-51). The results of GSEA analysis suggested that the role of GABRP in pancreatic cancer was closely connected with an inflammatory immune response.
Immune checkpoint therapy (ICT) has demonstrated great clinical benefit in several cancer treatments in recent years, although pancreatic cancer studies reported limited clinical response to ICT (52-54). Thus, a comprehensive and in depth understanding of tumor immune environment and immune infiltration in pancreatic cancer could greatly assist in the development of effective immunotherapeutic strategies for pancreatic cancer. We compared the immune infiltration landscape between pancreatic cancer and normal controls and found an obvious increase in macrophages, especially M0 macrophages, and CD4^+ ^T cells, and a decrease in CD8^+^ T cells. The infiltration of macrophages into solid tumors was associated with poor survival outcomes while the infiltration of CD8^+^ T cells was correlated with better prognosis (55,56). We further explored this subject and found a relationship between GABRP expression and immune infiltration. The expression of GABRP was positively correlated with the infiltration of B cells, CD8^+^ T cells, macrophages, neutrophils and dendritic cells, and GABRP copy number alteration was associated with B cells, CD4^+^ T cells, macrophages, and neutrophil infiltration in PAAD. In addition, GABRP expression was associated with different immune cell markers, and the differences in numbers of infiltrating immune cells could significantly affect the prognosis of pancreatic cancer patients. All the above results indicated that GABRP was correlated with immune infiltration which affected the prognosis of pancreatic cancer.
The immune cells in the pancreatic cancer microenvironment play a significant role in tumorigenesis and progression (57). Given the impact of GABRP on immune infiltration, we further explored the relationship between GABRP expression and the components of the pancreatic cancer microenvironment using ESTIMATE to evaluate 22 types of immune cells and stromal cells in tumor samples. We found that higher GABRP expression was significantly correlated with higher immune score and stromal score and lower tumor purity, which was in accordance with the analysis of immune infiltration via TIMER. Inflammation is intimately associated with immune responses, and inflammatory infiltration into the pancreatic cancer microenvironment contributed to tumor growth and metastasis (58). Pancreatic cancer is characterized by a pathological hallmark of stroma desmoplasia, which creates a mechanical barrier that reduces chemotherapy efficacy and results in low immune cell infiltration (57). According to our results, higher GABRP expression was correlated with increased density of immune cells and stromal cells, which resulted in lower tumor purity and complexity of the tumor microenvironment. The above results indicated an underlying connection between GABRP and tumor microenvironment. Higher GABRP expression was markedly associated with higher expression of six immune checkpoint genes, indicating that GABRP expression might potentially affect the efficacy of immunotherapy (59). The above correlation might provide new insights into immunological treatment for pancreatic cancer patients who received limited benefit from immune checkpoint inhibitors.
In summary, GABRP is a potential prognostic marker for pancreatic cancer and correlated with immune infiltration and tumor microenvironment. However, the present study still has some limitations. While the conclusions of this study were based on GABRP mRNA and protein expression, our results only indicated an overall relationship between GABRP and immune infiltration and tumor microenvironment. Further studies are needed to verify the epigenetics regulation of GABRP expression level and the underlying mechanism of the specific interactions between GABRP, immune infiltration and tumor environment.
In the present study, we demonstrated that GABRP was overexpressed in pancreatic cancer and that enhanced GABRP expression correlated with poor prognosis. Promoter hypomethylation and the regulation of hsa-miR-3655 might contribute to the overexpression of GABRP in pancreatic cancer. GABRP expression was correlated with immune infiltration and microenvironment. Further studies are needed to elucidate the molecular mechanism of the regulatory role of GABRP on the immune response.
We would like to thank Gary Bentley for his help in polishing our paper.
Funding: This research was funded by Beijing Natural Science Foundation (7212157), CAMS Innovation Fund for Medical Sciences (2016-I2M-3-007), National Natural Science Foundation of China (81803584, 81703536), Technology Major Projects for “Major New Drugs Innovation and Development” (2018ZX09711001-005-025, 2018ZX09711001-012).
*Ethical Statement: *The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). Institutional ethical approval and informed consent were waived.