In Silico Transcriptomic Analysis of the Chloride Intracellular Channels (CLIC) Interactome Identifies a Molecular Panel of Seven Prognostic Markers in Patients with Pancreatic Ductal Adenocarcinoma.

BACKGROUND: Pancreatic ductal adenocarcinoma (PDAC) is associated with poor prognosis. In this context, the identification of biomarkers regarding the PDAC diagnosis, monitoring, and prognosis is crucial.
OBJECTIVES: The purpose of the current study was to investigate the differential gene expression profile of the chloride intracellular channel (CLIC) gene family network in patients with PDAC, in order to suggest novel biomarkers.
METHODS: In silico techniques were used to construct the interactome of the CLIC gene family, identify the differentially expressed genes (DEGs) in PDAC as compared to healthy controls, and evaluate their potential prognostic role.
RESULTS: Transcriptomic data of three microarray datasets were included, incorporating 114 tumor and 59 normal pancreatic samples. Twenty DEGs were identified; eight were up-regulated and twelve were downregulated. A molecular signature of seven genes (Chloride Intracellular Channel 1 - CLIC1; Chloride Intracellular Channel 3 - CLIC3; Chloride Intracellular Channel 4 - CLIC4; Ganglioside Induced Differentiation Associated Protein 1 - GDAP1; Ganglioside Induced Differentiation Associated Protein 1 Like 1 - GDAP1L1; Glutathione S-Transferase Pi 1 - GSTP1; Prostaglandin E Synthase 2 - PTGES2) were identified as prognostic markers associated with overall survival. Positive correlations were reported regarding the expression of CLIC1-CLIC3, CLIC4-CLIC5, and CLIC5-CLIC6. Finally, gene set enrichment analysis demonstrated the molecular functions and miRNA families (hsa-miR-122, hsa-miR-618, hsa-miR-425, and hsa-miR-518) relevant to the seven prognostic markers.
CONCLUSION: These outcomes demonstrate a seven-gene molecular panel that predicts the patients' prospective survival following pancreatic resection for PDAC.

INTRODUCTION

Pancreatic cancer (PC) is one of the leading causes of cancer-related death and the fourth cause of cancer mortality in the US [1, 2]. The majority of the cases diagnosed with PC are ductal adenocarcinomas (PDAC), most frequently located in the head of the pancreas [3, 4] and associated with poor prognosis [5]. Depending on the degree of differentiation and the tumor microenvironment, the malignancy may present poorly to well-formed glands or infiltrating cells forming sheets [3, 4]. Besides the great research efforts, the mortality rate regarding PC is increasing steadily, thus being projected that by 2030, it will represent the second cancer-related cause of mortality [6].

In the same context, it is well accepted that tumor development and growth depend on the tumor microenvironment and metabolism [7]. Chloride intracellular channels (CLICs) are a family of ion channels that have been found overexpressed in several malignancies [8, 9]. Although the available data on the potential expression or role of CLIC members in PDAC is limited [10], it has been demonstrated that especially CLIC1, CLIC3 and CLIC4 have significant roles in cancer and more specifically in tumor metastasis and aggressiveness, cell proliferation, along with epithelial to mesenchymal transition [8, 9]. This evidence prompted us to study the role of CLIC gene family members and their interactome in PDAC.

The purpose of the current study was to evaluate the differential expression of CLIC1-6 genes in PDAC compared with healthy tissue using microarray data of three independent PDAC datasets. Besides, we investigated the CLIC protein interactors and performed Kaplan-Meier analysis in order to identify novel candidate genes that, in conjunction with CLIC genes, may be used as biomarkers. Finally, we performed gene set enrichment analysis (GSEA) in order to identify the relative molecular functions and regulating miRNA families.

MATERIALS AND METHODS

Construction of the CLIC Interactome

In silico analysis of the network of interactors related to the CLIC gene family was performed to identify associated partners of the CLIC genes involved in PDAC. The CLIC gene network was produced using the GeneMANIA platform (http://genemania.org/) [11]. GeneMANIA is a software that predicts the function of one or more genes to construct a connectivity network based on gene ontology algorithms. The functions of the proteins related to the identified genes were extracted from the portal GeneCards (http://www.genecards.org/), which is a database of systematic information on the human genome.

In Silico Evaluation of the Transcriptomic Profile of the CLIC Gene Interactome in PDAC

We used the PubMed Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/gds) to evaluate the expression profile of the CLIC gene family network in patients with PDAC compared with healthy controls. PubMed GEO is a repository of publicly available curated gene expression datasets, along with original series and platform records. All microarray data used in this study were downloaded during July 2019. The differentially expressed genes (DEGs) were identified using three independent PDAC microarray datasets (GSE16515, GSE15471, GSE32676), incorporating 173 samples (36+36+42=114 primary tumor samples and 16+36+7=59 normal control samples). The PDAC samples were excised from the primary tumor and the normal control samples were provided by excision of the adjacent normal pancreatic tissue. All three datasets were generated using the [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array platform. The gene expression data were log-transformed. The null hypothesis suggested that the gene expression profile between tumor and normal tissue samples was comparable. Genes were considered as significantly overexpressed or underexpressed when the p-value corresponded to p < 0.05. The available baseline characteristics of the patients included in each dataset are demonstrated in Table .

Evaluation of the DEGs for Potential Prognostic Value in PDAC

We constructed survival curves, using the publicly available survival data of The Cancer Genome Atlas (TCGA), in order to evaluate the potential prognostic significance of the identified DEGs. To perform the analysis, we used the PROGgeneV2 Prognostic Database portal (Indiana University) [12, 13] and we calculated the median value of gene expression level to group the patients in either high expressing (above median) or low expressing (below median) each significantly differentially expressed gene.

Gene Set Enrichment Analysis (GSEA) Regarding the Molecular Functions and Regulating miRNAs of the Identified Prognostic Markers

The GSEA of Gene Ontologies (GO) was performed using the electronic tool ToppFun of the ToppGene platform (https://toppgene.cchmc.org/). ToppFun uncovers the significant enrichments regarding the molecular functions and regulating miRNAs of input genes based on the transcriptome, regulome, proteome, along with phenotype data [14]. The significant enrichments were further evaluated by calculating their false discovery rates (FDR). The analyses were performed during August 2019.

Statistical Analysis

All analyses were performed using GraphPad Prism 8.0 for Mac (GraphPad Software, San Diego, CA). The normal distribution of the data was performed by calculating the D’Agostino and Pearson Omnibus normality test. Comparisons of gene expression levels were performed with a two-tailed unpaired t-test regarding parametric data and Mann-Whitney U-test regarding nonparametric data. All p values were corrected for multiple comparisons by calculating the Q statistic (Benjamini-Hochberg). Differences demonstrating a Q < 0.05 were considered significant. Correlations were assessed by calculating the Pearson or the Spearman’s rank (ρ) correlation coefficients for parametric or non-parametric data, respectively. Deming regression analysis was performed in order to evaluate cause-and-effect relationships among significant genes. Kaplan-Meier survival curves were constructed using the PROGgeneV2 Prognostic Database (Indiana University) software. Differences were considered significant (rejection of the null hypothesis) with a p ≤ 0.05.

RESULTS

A trial flow of the current study is demonstrated in (Fig. ). The members of the CLIC gene family network that were extracted from the GeneMania platform are demonstrated in Table and Fig. (). A total of 26 interacting proteins were revealed through the construction of the CLIC interactome in homo sapiens.

Out of the 26 proteins, sufficient PubMed GEO data for their gene expression level was available regarding 23 genes (88.5%). A total of 20 DEGs were identified. No data was available regarding the GSTA2, GSTA5 and GSTT2B genes. Subsequently, we investigated the PubMed GEO database regarding the gene expression profile of the above genes. The cumulative outcomes are presented in Table . We identified eight overexpressed genes and twelve underexpressed genes (Table ; Figs. ( and ) and Supplementary Fig. ). No significant difference was reported regarding the expression level of three genes (Table ).

Significant Correlations Among DEGs

CLIC1 expression was found to be positively correlated with CLIC3 expression (p<0.0001 and Spearman’s r=0.5319); (Fig. ). The expression level of CLIC5 was also found positively correlated with CLIC4 (p<0.0001 and Spearman’s r=0.4015); (Fig. ) and CLIC6 (p<0.0001 and Spearman’s r=0.2808); (Fig. ). Deming regression analysis revealed the equations describing these correlations. The positive correlation between CLIC1 and CLIC3 is described by the equation: CLIC3=2.044*CLIC1-16.22. Moreover, the positive correlation between CLIC4 and CLIC5 is demonstrated by the equation: CLIC5=2.044*CLIC4-16.22. In addition, the correlation between CLIC5 and CLIC6 is described by the following equation: CLIC6=0.4278*CLIC5+2.883.

CLIC1, CLIC3, CLIC4, GDAP1, GDAP1L1, GSTP1, PTGES2 Expression were Prognostic Indicators of PDAC Patients

We evaluated whether the dysregulation of DEGs in PDAC could affect patient survival outcomes. PDAC data, along with gene expression and clinical information from The Cancer Genome Atlas (TCGA), were facilitated to assess their prognostic significance. A total of 173 pancreatic cancer patients were incorporated in this analysis. Kaplan-Meier survival curves were constructed for all DEGs. The median survival was 684 days for the CLIC1 low expression group, and it dropped to 498 days in CLIC1 high expression group (HR: 1.56 [95% CI: 1.16-2.1]; p=0.003; (Fig. ). CLIC3 low expression group had an increased median survival period, 666 days, as compared with the median survival of the high expression group, which was 498 days (HR: 1.16 [95% CI: 1.05-1.28]; p=0.0037); (Fig. ). CLIC4 low expression group was also associated with enhanced survival of 684 days compared to 593 days of the high expression patients (HR: 1.44 [95% CI: 1.12, 1.85]; p=0.005); (Fig. ). Furthermore, the median survival for the GDAP1 high expression group was 684 days, and for the GDAP1 low expression group was 532 days (HR: 0.71 [95% CI: 0.57-0.90]; p=0.004); (Fig. ). In addition, the median survival was 634 days regarding the GDAP1L1 high expression group and 603 days for the GDAP1L1 low expression group (HR: 0.95 [95% CI: 0.92, 0.98]; p=0.003); (Fig. ). Patients with a high level of PTGES2 expression were associated with increased median survival compared to the PTGES2 low expression group (684 days versus 592 days, respectively; HR: 0.6 [95% CI: 0.42, 0.87]; p=0.007); (Fig. ). Finally, the median survival was increased in the GSTP1 low expression group compared to the high expression group (614 versus 607 days; HR: 1.41 [95% CI: 1.12, 1.78]; p=0.004); (Fig. ). These outcomes demonstrated that CLIC1, CLIC3, CLIC4, GSTP1 were adverse factors while GDAP1, GDAP1L1, PTGES2 were beneficial factors regarding the median survival of PDAC patients.

The seven prognostic markers underwent GSEA. The top five enriched Gene Ontology terms for molecular functions are presented in Table , along with the miRNAs that

regulate the seven prognostic markers. Regulation of glutathione transferase, transferase activity (transferring alkyl or aryl groups), chloride channel, chloride transmembrane transporter, and anion channel activity represented the most significant molecular functions relevant to the seven markers. Finally, the GSEA indicated that the members of the hsa‐miR‐122, hsa‐miR‐618, hsa‐miR‐425, and the hsa‐miR‐518 miRNA families are important regulators of the seven prognostic markers.

DISCUSSION

The present in silico study evaluated the gene expression profile of the CLIC interactome in PDAC by incorporating data provided by three pancreatic cancer microarray datasets, thus providing enhanced accuracy at the molecular level, as compared with studies based on a single dataset. In fact, a total of 173 samples were included and analyzed in the present study. Furthermore, we have constructed the interactome of the CLIC gene family through which we report 20 DEGs related to PDAC. These novel gene candidates are CLIC1, CLIC3, CLIC4, CLIC5, CLIC6, GSTO1, GSTP1, GDAP1, GDAP1L1, HPGDS, GSTA1, GSTA3, GSTM2, GSTM4, GSTM5, GSTT1, GSTT2, GSTZ1, PTGES2, and TPRN. The current evidence should be further investigated to unveil the role of the DEGs in the biology of PDAC, along with their potential use as drug targets.

Our analysis also demonstrated that PDAC was marked by dysfunctions of chloride transport, glutathione derivative metabolic process, and regulation of the anion channel activity. CLIC1 is a protein that localizes primarily to the cell nucleus. It has been implicated in the pathogenesis of different types of cancer and has been found overexpressed in various types of cancer, including gastric, colorectal, and hepatocellular cancer [9, 15]. Our findings regarding the overexpression of CLIC1 in PDAC and its prognostic significance are in accordance with the outcomes of a previous study by Jia et al. [10]. CLIC3 also localizes in the membrane of the cell nucleus and has been implicated in different metastatic processes [16]. It has been reported that CLIC3 facilitates the translocation of α5β1 integrin from the late endosomes/lysosomes to the cell surface, thus promoting cell motility [17]. In the same context, the high expression level of CLIC3 was associated with decreased survival and a higher rate of metastasis [17]. In the present study, CLIC3 was overexpressed in PDAC and was associated with poor survival. Additionally, the expression level of CLIC4 was higher in PDAC compared to normal samples, while the higher CLIC4 expression was associated with poorer survival. These findings were in accordance with the only previous study [18] assessing the expression profile of CLIC4 in PDAC, which also demonstrated that the expression of CLIC4 associated with poor survival. Furthermore, our outcomes suggest positive correlations between the expression level of CLIC1-CLIC3, CLIC4-CLIC5, and CLIC5-CLIC6 that are described by certain equations.

According to the literature, GSTP1 facilitates xenobiotic metabolism and tumorigenesis [19, 20]. The drug metabolism through cytochrome P450 is also included in its relevant pathways. In the present study, we found that GSTP1 was upregulated in PDAC tissue samples. Moreover, the overexpression of GSTP1 was associated with decreased survival. GDAP1 and GDAP1L1 are paralogue protein-encoding genes regulating the mitochondrial network by promoting mitochondrial fission [21]. Both GDAP1 and GDAP1L1 were downregulated in PDAC and were associated with higher survival. Moreover, PTGES2 has been implicated in the pathogenesis of different types of cancer [22, 23]. According to the present outcomes, PTGES2 was a beneficial factor regarding the median survival in PDAC patients. The present study provided the first evidence, to the best of our knowledge, supporting the prognostic value of GSTP1, GDAP1, GDAP1L1, and PTGES2 in patients with PDAC.

Finally, the molecular functions mediated by the seven prognostic markers as revealed by the GSEA were associated with glutathione transferase, transferase activity (transferring alkyl or aryl groups), chloride channel, chloride transmembrane transporter, and anion channel activity. Furthermore, the GSEA predicted that the seven prognostic markers might be targets of the hsa‐miR‐122, hsa‐miR‐618, hsa‐miR‐425, and the hsa‐miR‐518 miRNA family members. MicroRNAs have been proposed as important factors implicated in the biology of PDAC, while the miRNA‐based therapy of PDAC has been suggested as a novel area of research [24]. Overall, the present study demonstrates, for the first time, to the best of our knowledge, the potential role of these miRNA families in PDAC biology.

The present study presents a seven-gene signature that can provide important prognostic information regarding PDAC patients, along with the related molecular functions and regulating miRNA families. This panel can be possibly used in conjunction with the current staging systems in order to provide enhanced prognostic information. In addition, the current study provides evidence regarding the expression profile of the CLIC interactome in patients with PDAC. This valuable information should be further investigated to enhance our level of knowledge regarding the PDAC biology and ameliorate our treatment options.

The limitations of the current study are (1) the processing of the tissue samples in three different laboratories, (2) the lack of available data regarding the baseline characteristics, PDAC TNM status and neoadjuvant treatment of all the included patients to perform further multivariate analyses, and (3) the lack of mutation/alteration data that provides stronger evidence compared to gene expression data. Nonetheless, the strengths of the present study are (1) the clear protocol, (2) the inclusion of three datasets and a large number of the incorporated tissue samples, thus increasing the level of evidence, (3) the use of the same Affymetrix chip in all datasets, (4) the analysis of survival data related to the gene expression profiles, and (5) the performance of GSEA which demonstrated the molecular functions and miRNA families related to the identified prognostic markers.

CONCLUSION

In the present study, we identified for the first time 20 novel genes differentially expressed in PDAC, along with the significant correlations between DEGs. Moreover, we identified seven genes with prognostic value in the context of PDAC. We also demonstrated the predicted molecular functions and miRNA families relevant to these prognostic markers. Given the in silico nature of the present study, further translational research is necessary to fully unveil the potential benefit from our outcomes regarding the PDAC diagnosis and treatment.

Table 1

Gene symbols and description of the CLIC interactome members.

Gene Symbol	Gene Description
		CLIC1	Chloride Intracellular Channel 1
CLIC2	Chloride Intracellular Channel 2
CLIC3	Chloride Intracellular Channel 3
CLIC4	Chloride Intracellular Channel 4
CLIC5	Chloride Intracellular Channel 5
CLIC6	Chloride Intracellular Channel 6
GDAP1	Ganglioside Induced Differentiation Associated Protein 1
GDAP1L1	Ganglioside Induced Differentiation Associated Protein 1 Like 1
GSTA1	Glutathione S-Transferase Alpha 1
GSTA2	Glutathione S-Transferase Alpha 2
GSTA3	Glutathione S-Transferase Alpha 3
GSTA4	Glutathione S-Transferase Alpha 4
GSTA5	Glutathione S-Transferase Alpha 5
GSTM2	Glutathione S-Transferase Mu 2
GSTM4	Glutathione S-Transferase Mu 4
GSTM5	Glutathione S-Transferase Mu 5
GSTO1	Glutathione S-Transferase Omega 1
GSTO2	Glutathione S-Transferase Omega 2
GSTP1	Glutathione S-Transferase Pi 1
GSTT1	Glutathione S-Transferase Theta 1
GSTT2	Glutathione S-Transferase Theta 2
GSTT2B	Glutathione S-Transferase Theta 2B
GSTZ1	Glutathione S-Transferase Zeta 1
HPGDS	Hematopoietic Prostaglandin D Synthase
PTGES2	Prostaglandin E Synthase 2
TPRN	Taperin

Table 2

Differential gene expression of the CLIC gene network in PDAC as compared to healthy controls.

Gene Symbol	Fold Changes (Actual)	Fold Changes (Hodges-Lehmann)	P values	Q values
Upregulated
CLIC1	1.640	1.375	< 0.0001	2.8636e-005
CLIC3	NA	NA	< 0.0001	2.8636e-005
CLIC4	0.5384	0.4133	0.0400	0.0065
CLIC6	0.3288	0.2223	0.0137	0.0027
GSTO1	0.7311	0.6437	< 0.0001	2.8636e-005
GSTP1	NA	NA	< 0.0001	2.8636e-005
HPGDS	0.8110	0.7026	0.0013	0.0003
TPRN	NA	NA	0.0013	0.0003
Downregulated
CLIC5	-0.4764	-0.2068	0.0153	0.0027
GDAP1	-0.4695	-0.4457	0.0124	0.0026
GDAP1L1	-2.278	-0.9061	< 0.0001	2.8636e-005
GSTA1	-3.436	-2.010	< 0.0001	2.8636e-005
GSTA3	-3.106	-1.994	< 0.0001	2.8636e-005
GSTM2	NA	NA	< 0.0001	2.8636e-005
GSTM4	-0.2976	-0.2804	0.0157	0.0027
GSTM5	-1.353	-0.9681	< 0.0001	2.8636e-005
GSTT1	-0.6718	-0.4820	0.0413	0.0065
GSTT2	-0.9199	-0.7902	0.0004	0.0001
GSTZ1	-0.7449	-0.6223	< 0.0001	2.8636e-005
PTGES2	-0.5889	-0.4797	< 0.0001	2.8636e-005
Not significantly different
CLIC2	0.4399	0.2214	0.3191	0.0437
GSTA4	-0.02907	-0.1858	0.2420	0.0347
GSTO2	-0.3611	-0.3582	0.0770	0.0116

Table 3

Enrichment analysis of gene ontologies for the prognostic factors; Top five relevant molecular functions and regulating miRNA families are presented.

Molecular Functions
-	ID	Name	Source	p Value	FDR B&H	FDR B&Y	Bonferroni	Genes from Input	Genes in Annotation
1	GO:0004364	glutathione transferase activity	NA	1.604e-16	7.380e-15	3.260e-14	7.380e-15	6	34
2	GO:0016765	transferase activity, transferring alkyl or aryl (other than methyl) groups	NA	8.095e-15	1.862e-13	8.223e-13	3.723e-13	6	63
3	GO:0005254	chloride channel activity	NA	2.723e-6	3.803e-5	1.680e-4	1.253e-4	3	81
4	GO:0015108	chloride transmembrane transporter activity	NA	4.001e-6	3.803e-5	1.680e-4	1.841e-4	3	92
5	GO:0005253	anion channel activity	NA	4.134e-6	3.803e-5	1.680e-4	1.902e-4	3	93
Regulating miRNAs
-	ID	Name	Source	p Value	FDR B&H	FDR B&Y	Bonferroni	Genes from Input	Genes in Annotation
1	hsa-miR-122:mirSVR highEffct	hsa-miR-122:mirSVR conserved highEffect-0.5	Micro RNA.org	2.015e-5	7.938e-3	5.753e-2	1.588e-2	3	607
2	hsa-miR-618:mirSVR highEffct	hsa-miR-618:mirSVR nonconserved highEffect-0.5	Micro RNA.org	1.124e-4	2.210e-2	1.602e-1	8.859e-2	3	1083
3	hsa-miR-425*:mirSVR highEffct	hsa-miR-425*:mirSVR nonconserved highEffect-0.5	Micro RNA.org	1.403e-4	2.210e-2	1.602e-1	1.105e-1	2	188
4	hsa-miR-518f:mirSVR highEffct	hsa-miR-518f:mirSVR nonconserved highEffect-0.5	Micro RNA.org	2.264e-4	2.739e-2	1.985e-1	1.784e-1	2	239
5	hsa-miR-518c:mirSVR highEffct	hsa-miR-518c:mirSVR nonconserved highEffect-0.5	Micro RNA.org	2.994e-4	2.949e-2	2.137e-1	2.359e-1	2	275

Abbreviations: GO=Gene Ontologies; FDR=False Discovery Rate; B&H: Benjamini–Hochberg; B&Y: Benjamini–Yekutieli.