Identification of Critical Pathways and Potential Key Genes in Poorly Differentiated Pancreatic Adenocarcinoma

Introduction

Pancreatic adenocarcinoma (PDAC) accounts for more than 95% of human pancreatic cancer, is a highly aggressive tumor whose mortality practically parallels incidence.¹ Five-year survival of PDAC patients is no more than 5% in the world.² By the year 2030, PDAC is predicted to be the second leading cause in tumor-related deaths.³ Less than 20% of patients can be completely resected, and most patients will get early progression even after potential curative resection, with a five-year survival rate of less than 25%.^4,5 Despite advances in surgery and adjuvant chemotherapy for pancreatic cancer in recent years, the survival rate was not significantly improved. Therefore, it is urgent to explore effective targets to improve the prognosis of PDAC.

Aberrant gene expression is a common theme of human disease, corresponding mRNA play pivotal roles in the occurrence and progress of tumor-related disease.⁶ The emergence of high-throughput sequencing technology allows investigators to detect a complete data of the mutational and transcriptional landscape of most tumors.⁷ To date, numerous investigators have performed comprehensive next-generation sequencing analysis on PDAC samples, aiming to determine novel biological markers whose mutations are involved in tumor formation. The RNA-seq has been considered as a promising tool to filtrate the hub genetic markers of PDAC progression.⁷ Bioinformatic analysis of the RNA-seq data and transcriptomes has found a group of altered genes in PDAC, which dramatically promoted our understanding of the biological processes of tumorigenesis as well as remedy in human tissue samples.⁸ It is demonstrated that SMAD3 take an active part in inducing epithelial-mesenchymal transition (EMT) and was considered as an effective biomarker indicating poor prognosis.⁹ Some study revealed that THBS1 enriched in the phosphoinositide 3-kinase-Akt pathway could promote the development and progression of PDAC.¹⁰ Some scholars^6,11 have considered KRAS, TP53, DPC4, and SMAD4 as tumor driver genes that may result tumor recurrence and antitumor immunity in PDAC, meanwhile, targeted drugs against these pathways can be developed.^12,13 However, despite enormous progress in understanding the tumor initiation and progression of PDAC has made in gene level, there exist no widely recognized prognostic biomarkers with high sensitivity and specificity which could predict prognosis precisely and can be applied for postoperative individualized treatment. At present, RNA-seq technology combined with bioinformatics analysis enable it a promising way to comprehensively explore the aberrations of mRNA expression across the formation and development of PDAC. Investigation into these genes could provide valuable information that can help to devise the optimal treatment for patients and even predict disease relapse.

In the present study, we carried out a complicated differential analysis of PDAC and paired normal samples by integrating RNA-seq data both from our cohort aiming to elucidate the differentially expressed genes (DEGs) and explore the potential biofunctions by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses. Furthermore, the protein-protein interaction (PPI) network was constructed, and five hub genes were identified from it. Furthermore, survival analysis revealed that the hub genes were obviously correlated with the prognosis of PDAC. The present study may provide new insight into the understanding of PDAC pathogenesis and the identified hub genes may serve as potential targets for diagnosis and treatment.

Materials and Methods

Sample Preparation and RNA-seq Analysis

Seven patients underwent radical resection of the primary tumor was included in this study and all of them were confirmed by postoperative pathology as poorly differentiated adenocarcinoma of the pancreas. PDAC tissues and their adjacent normal tissues were collected for transcriptome sequencing analysis and were stored in RNA later solution (Ambion, USA) immediately after surgery and stored at –80°C prior to RNA extraction. This study has been approved by the Ethics Committee of Henan Provincial People’s Hospital, and all patients have signed informed consent before the surgery.

Total RNA of each specimen was isolated using TRIzol reagent (the RNeasy Mini Kit; Qiagen, Germany) following the manufacturer’s instructions. An Agilent 2100 RNA Nano 6000 Assay Kit (Agilent Technologies, CA, USA) was used to qualify the integrity and concentration of total RNA. Later, all the specimens were sent to BGI (The Beijing Genomics Institute) Corporation (Wuhan, China) for further RNA-seq detection and analysis via BGISEQ-500 sequencer.

Identification of DEGs

DEGs between PDAC and non-PDAC samples from the generated data (RNA-seq dataset) were screened using R scripts, only the data with |log2FC| ≥ 1 and Q-value ≤0.05 were left for as DEGs for further investigation.

GO Analysis, Enrichment, KEGG Pathway Analyses

To uncover the functional roles of DEGs, we determined the transcriptional profiles acquired using RNA-seq and performed GO and KEGG enrichment analyses. Based on the DR.TOM system of BGI, we analyzed the functional enrichment of DEGs, concentrating on the pathways interrelated with these genes that were enriched in the terms biological process, molecular function and cellular component. The KEGG pathway analysis was applied to determine the significant pathways of the DEGs. The threshold of the hypergeometric distribution test for default enrichment results was 0.05.

Protein-protein Interaction (PPI) Network Analysis and Identification of Hub Genes

Search tool for the retrieval of interacting genes (STRING) database is online tool supplied to illustrating the PPI information. To further explore the potential interplay among those DEGs, we mapped the DEGs to STRING, and only interactions that enjoyed a minimum required combined score >0.4 were set as significant. Subsequently, the plug-in MCODE (version 1.5, http://apps.cytoscape.org/apps/mcode) was carried out to identify significant modules of the constructed network. Furthermore, the key genes were identified based on the PPI network by using cytoHubba (version 0.1, apps.cytoscape.org/apps/cytohubba), which is another plug-in of Cytoscape.

Survival Analysis and Correlation of DEGs

To study prognostic implications of hub genes, online server PROGgeneV2 was employed to conduct a survival analysis of data from TCGA (https://cancergenome.nih.gov/). In the public available TCGA datasets, PDAC samples were dichotomized into hub gene high- and low-expression groups by the median transcript per million (TPM) value as a cutoff, respectively. Cox regression analysis was carried out to unravel the correlation between hub gene expression and prognosis in PDAC and the log rank test was used for the hypothesis tests. The genes interrelated with OS were carried out with further analysis, incorporating Pearson correlation analysis as well as analysis of expression levels between tumor and normal tissues via GEPIA Online tools. A p-value of <0.05 is considered statistically significant.

Quantitative Reverse Transcriptional PCR (qRT-PCR)

Total RNA of the seven patients’ specimens were isolated using the E.Z.N.A. Total RNA Kit I (Invitrogen) following standard procedures. Complementary DNA (cDNA) was synthesized subsequently by using the Prime Script RT Reagent Kit (Takara Bio, Japan). Real-time PCR was performed on the ABI 7500 Touch RealTime PCR Detection System (ABI. USA) to determine the mRNA expression level of screened hub genes according to the comparative Ct method. GAPDH was set as a normalization control. The sequences of primers are presented in Supplementary material.

Results

DEGs Associated with PDAC

Paired differential gene expression analysis by RNA-seq revealed that 126 genes were found to be differentially expressed in seven paired PDAC samples from our cohort. Among the DEGs, 82 were upregulated genes and 44 were downregulated genes (Supplementary material). Heatmaps of the distributions of DEGs are presented in Figure 1.

Figure 1 Heat map of the differentially expressed genes in seven PDAC samples.

GO Enrichment Analyses of DEGs

To further elaborate the function of DEGs in human PDAC, GO annotations and KEGG pathway enrichment analysis was carried out. GO analysis determined 51 terms significantly enriched, incorporating terms correlated with molecular function (BP), biological process (CC) and cellular component (MF) (Figure 2A–C). Changes in BP were significantly enriched in cell adhesion, blood coagulation, proteolysis, keratinization, extracellular matrix organization. For CC, the DEGs were mainly concerned with extracellular region, extracellular space, integral component of membrane, plasma membrane, extracellular exosome. Moreover, for category MF, the main functions of proteins encoded by the DEGs were calcium ion binding, peptidase activity, metallopeptidase activity, extracellular matrix structural constituent, structural molecule activity.

Figure 2 GO and KEGG pathway enrichment analyses for DEGs. (A) Chord plot of GO terms under BP category. (B) Chord plot of GO terms under CC category. (C) Chord plot of GO terms under MF category. (D) Dot plot of KEGG pathway enrichment analyses.

KEGG Pathway Analysis

The pathway terms associated with the obtained DEGs were tested for significance using Fisher’s exact test (p<0.05). Pathway terms with prominent enrichment were figured out and DEGs were opt for the pathway category analysis. Interestingly, most of the enriched pathways were correlated to the signal transduction and specific types of cancer-related pathways, including chemical carcinogenesis, IL-17 signaling pathway. In addition, the DEGs were also enriched in protein digestion and absorption, drug metabolism-cytochrome p450 and metabolism of xenobiotics by cytochrome p450 (Figure 2D).

Protein-protein Interaction (PPI) Network Analysis and Functional Module Analysis

This network consists of convinced interactions from curated databases and those that were experimentally determined; involved interactions analysis with regards to gene neighborhood, gene fusions and gene co-occurrence; and text-mining, co-expression, respectively. The direct or indirect function interaction of proteins was visualized with the help of String database (Figure 3A). Furthermore, clustering analysis carried out by Cytoscape plug-in MCODE identified four functional modules from the whole network (Figure 3B–E), and the subsequently pathway enrichment analyses for genes of the four modules which could help to comprehend the potential biofunctions showed that these modules were mostly involved in tumorigenesis, cell adhesion and energy metabolism incorporating cell adhesion molecules (CAMs), pathways in cancer, cell cycle, PPAR signaling pathway, estrogen signaling pathway, cholesterol metabolism (Table 1).

Table 1 KEGG Enrichment Analysis of Top Four Modules Identified from PPI Network

Figure 3 PPI network construction and functional module analysis. The whole PPI network (A). Network of functional module 1 (B). Network of functional module 2 (C). Network of functional module 3 (D). Network of functional module 4 (E).

Survival Curves, Expression Levels and Correlation Analysis

CEACAM5, KRT6A, KRT6B, KRT7, KRT17 were suggested to have prognostic relevance with PDAC, patients with high expression of the five genes may have a worse outcome in survival condition (Figure 4). We further explored the relationship between mRNA expression level and human tumor samples. Expression levels of CEACAM5, KRT6A, KRT6B, KRT7, KRT17 are displayed in Figure 5. CEACAM5, KRT6A, KRT6B, KRT7, KRT17 demonstrate high expression levels in PDAC tissues. And then, we carried out Pearson correlation analyses among the expressed genes revealed in Figure 6. It was demonstrated that CEACAM5 was positively correlated with KRT7 (p=1.7e−08, R=0.3); KRT7 was also positively correlated with KRT6A, KRT6B, KRT17: KRT6A (p-value=0, R=0.43); KRT6B (p-value=2.7e−15, R=0.41); KRT17 (p-value=0, R=0.57).

Figure 4 Prognostic survival analysis of CEACAM5, KRT6A, KRT6B, KRT7, KRT17 in PDAC were analyzed by the online server PROGgeneV2 tool, data of survival analysis is from TCGA (https://cancergenome.nih.gov/).

Figure 5 Expression analysis of a CEACAM5, KRT6A, KRT6B, KRT7, KRT17 in normal and tumor pancreatic tissues.

Figure 6 Pearson correlation analyses of CEACAM5, KRT6A, KRT6B, KRT7 and KRT17. CEACAM5 was positively correlated with KRT7 (a); KRT7 was also positively correlated with KRT6A, KRT6B, KRT17 (b, c, d).

Expression Validation

To further validate the reliability of the RNA-seq data, we validated the expression of hub genes in cancers and noncancer samples using qRT-PCR. As a result, the relative expression level of CEACAM5, KRT6A, KRT6B, KRT7, KRT17 in PDAC samples were significantly higher compared to nonPDAC samples (p<0.05). This result was consistent with the RNA-seq analysis (Figure 7).

Figure 7 Validation of hub genes expression in PDAC and non-PDAC clinical samples using qRT-PCR. The relative expression level of the five hub genes.

Discussion

Tumorigenesis is sophisticated and usually accompanied by multiple genetic alterations at the transcription level.¹⁴ The high mortality rate of pancreatic cancer has made it a worldwide public health problem. Many clinicians are not able to make an early diagnosis, this might be caused by lack of identification of specific biomarkers, either in genes or proteins.¹⁵ Fortunately, the emergence of RNA-seq technology has provided a possibility for exploring mechanisms underlying tumor formation at the mRNA level as well as identification of effective molecular therapeutic targets. The combination of bioinformatics and high-throughput sequencing have produced a vast amount of gene expression profiles, which could perform statistical analysis in disease occurrence, progression, and recurrence, serving an important role in tumor analysis.

In the current study, 126 DEGs were identified based on transcriptome analysis. GO and KEGG enrichment analyses revealed that the DEGs were predominately associated with cell adhesion, integral component of membrane, signal transduction and chemical carcinogenesis, IL-17 signaling pathway. Furthermore, we performed an integrative analysis of potential prognostic association between PDAC and DEGs by PROGgeneV2. CEACAM5, KRT6A, KRT6B, KRT7, KRT17 have potential prognostic value for PDAC patients. The PPI network and module analysis were carried out to further explore functional sub-networks and CEACAM5, KRT6A, KRT6B, KRT7, KRT17 represented core hub genes in the network. Moreover, these genes are not only correlated with PDAC, but also to the signal transduction pathway.

CEACAM5 is a key member of the glycosylphosphatidylinositol (GPI)-linked immunoglobulin (Ig) superfamily, attached to the plasma membrane of cell.¹⁶ The CEACAM5 gene, encoding the most extensively studied tumor markers-the carcinoembryonic antigen (CEA), is notorious for it is vital role in the clinical work as a tumor marker for several malignancies.¹⁷ The encoded protein as an intercellular adhesion molecule, plays an irreplaceable role in several process incorporating cell-cell recognition, proliferation, differentiation and cancer invasion and metastasis.¹⁸ Furthermore, it is indispensable in forming tissue architecture and neovascularization. Expression of CEACAM5 has been revealed to have significant correlation with the PDAC patients’ lymph node status.¹⁹ It has been conclusively demonstrated that the overexpression of CEACAM5 will contribute to pancreatic cancer progression, particularly invasion and metastasis.²⁰ Knockdown of the CEACAM5 gene could effectively inhibit tumor proliferation of tumor-bearing mouse.²¹ Distant metastasis is an important sign of a poor prognosis in PDAC patients. Expression of CEACAM5 in PDAC tissues is higher than that in normal tissues, which may contribute to its invasion and it is consistent with previous studies. Our study found a significant difference of OS between overexpression of CEACAM5 and the rest, this may be a conceivable marker in predicting poor prognosis of patients with PDAC.

The keratins are the typical intermediate filament proteins of epithelia, revealing a remarkable degree of molecular diversity.²² As part structural basis of the epithelial, keratins play pivotal roles in the mechanical stability and integrity of epithelial cells and tissues, involving intracellular signaling pathways.²³ Keratins are abundantly expressed in hair and skin cells, its aberrant expression in other tissues may cause pathological process, involved in the progression of various tumors.²⁴ It has been extensively explored that the role of keratins in cancer development was related with epithelial-derived tumors such as bladder cancer,²⁵ lung cancer, colorectal cancer, urothelial cancer, and gastric cancer.²⁶ Here in this work, we focused on the roles of KRTA/6B, KRT7, KRT17 in PDAC.

KRT6 belongs to the type II intermediate filament protein, including three isoforms (KRT6A, KRT6B, and KRT6C).²⁷ KRT6A could be detected in ductal myoepithelial cells in human mammary duct.²⁸ KRT6A/KRT6B proteins have been unraveled to play roles in various cancers. Slusser-Nore et al found that the expression of KRT6 is regulated by the activation of the ERK1/2 pathway in human urothelial cells.²⁹ It is reported that the interaction of KRT6 with notch1 could aggravate the invasion of renal cell carcinoma.³⁰ However, the role of KRT in PDAC was unknown. Our RNA sequencing data demonstrated that mKRT6A/mKRT6B was hyperexpressed in PDAC. Based on the survival analysis, we found for the first time that patients with higher mKRT6A/mKRT6B expression showed shorter OS period. Thus, further studies are needed to investigate the function in detail and it may be the potential biomarker and novel therapeutic target for pancreatic cancer. In addition, we found that KRT6A/B is positively correlated with KRT7/KRT17 in tumor tissues; moreover, they were almost absent in the adjacent tissue counterparts.

Our recognition of KRT7 (keratin 7) in PDAC is concordant with previous findings. Actually, KRT7 is an intermediate filament protein, which is selectively expressed in simple epithelia.³¹ In the normal pancreas, only minimal expression of it could be detected in the duct, but are absent from acinar and endocrine cells,³² is considered as a characteristic marker for ductal cells of the hepato-pancreatic ductal tree.³³ It was demonstrated that KRT7 is significantly overexpressed in pancreatic cancer tissues, the coding protein was routinely detected by immunohistochemistry analysis for the diagnosis of PDAC.³⁴ KRT7 also participated in the context of differentiation of pancreatic cancer cells and it was investigated that KRT7 promoter is available for targeting PDAC cells either in vitro or in vivo. The AdK7Luc virus invasion assays further demonstrated that by incorporating the KRT7 promoter fragment in the adenovirus could increase cancer specificity.³² Studies have shown that overexpression of KRT7 in pancreatic cancer was markedly associated with liver metastasis and it was recommended as marker in distinguishing pancreatic and ampulla of Vater adenocarcinomas.³⁵ KRT7 has been extensively explored in cancer research, nevertheless, no studies have linked the expression of KRT7 to poor prognosis in pancreatic cancers. Our survival analysis showed for the first time that patients with high expression of KRT7 mRNA enjoyed significantly shorter OS than that of the low expression one, and the gene was positively correlated with KRT6A, KRT6B, and KRT17, respectively.

KRT17 belongs to type I intermediate filament, which is mainly involved in formation of basal epithelial cells. It was striking, however, that KRT17 was routinely regenerated and highly expressed in human tissue samples.³⁶ Mature pancreatic duct epithelium cell express no or little KRT17, nevertheless, it was determined that the KRT17 mRNA level was eminently upregulated in PDAC samples, comparing with normally differentiated pancreatic tissues, and an identical trend was detected in neoplastic cells, compared with normal pancreatic cells.³⁷ Furthermore, immunohistochemistry analysis revealed that KRT17 was overexpressed in tumor tissues compared with adjacent noncancerous tissues.³⁸ Due to it is function as structural constituents of the tissue, KRT17 has been reported to play a key role as a multifunctional promoter and oncogene to accelerate proliferation, metastasis, invasion of malignancies.³⁹ In addition, Chen et al discovered that inhibition of the proliferation of PANC‑1 human pancreatic cancer cells could be replicated by silencing of KRT17 through lentivirus‑mediated short hairpin RNA.⁴⁰ And it was found that upregulation of KRT17 promoted PDAC cell proliferation, invasion, and metastasis through EMT.⁴⁰ In addition, some studies raised concerns that overexpression of KRT17 mRNA was associated with a poor outcome for PDAC patients.⁴¹ In our study, RNA-seq analysis demonstrated that mRNA expression of KRT17 was markedly increased in cancer tissues compared to adjacent normal tissue. And the bioinformatics analysis predicted that overexpression is associated with tumor metastasis and poor prognosis, indicating that KRT17 could serve as a potential prognostic biomarker and a therapeutic target for PDAC.

Taken together, we identified 126 DEGs in patients with PDAC by employing RNA-seq method, and defined the potential novel therapeutic targets as well as evaluating their prognostic values in PDAC patients by bioinformatic analysis. GO-term and KEGG-pathway analysis were enriched in cellular component, cellular process related GO terms, and signal transduction as well as specific types of cancers related KEGG pathways. CEACAM5, KRT6A, KRT6B, KRT7, KRT17 may serve as novel biomarkers of poor survival in PDAC, and CEACAM5, KRT6A/B, KRT17 may even play a carcinogenic role which might provide effective anti-cancer target in PDAC therapy. Of course, our study has limitation, since TCGA database does not rovide sufficient information for us to identify poorly differentiated PDAC from the non-PDAC ones, the survival analysis of the recognized hub genes was carried out using samples of all pancreatic cancer subtypes. Thus, although this result contributes to the understanding of the biological relevance of the hub genes, further studies are needed to verify the function of hub genes.

Conclusion

In this work, 126 DEGs were identified between poorly differentiated pancreatic adenocarcinoma and the matched normal tissues. Among the DEGs, five genes were considered as hub genes which were related to many pathways of tumor progress and leading to poor prognosis. The result of this study may provide novel biomarkers for PDAC which is valuable for further investigation.