Derivation and Validation of the Potential Core Genes in Pancreatic Cancer for Tumor-Stroma Crosstalk

Background Pancreatic cancer is a fatal malignancy with a poor prognosis. The interactions between tumor cells and stromal cells contribute to cancer progression. Pancreatic stellate cells (PSCs) play a key role in tumor-stroma crosstalk of pancreatic cancer. The in-depth exploration for tumor-stroma crosstalk is helpful to develop novel therapeutic strategies. Our aim was to identify the potential core genes and pathways in tumor-stroma crosstalk. Methods 3 microarray datasets were from Gene Expression Omnibus (GEO). Differentially expressed genes (DEGs) were screened through bioinformatics analysis. Gene Ontology (GO) enrichment, Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment, and protein-protein interaction (PPI) network were used to obtain the biological roles of DEGs. The top 15 DEGs were explored by principal component analysis. We validated the top 15 DEGs expression in the tumor-stroma crosstalk model in which PSCs were treated with the mixture of Aspc-1 and Panc-1 supernatant. Results A total of 221 genes were filtered as DEGs for tumor-stroma crosstalk. The results of principal component analysis for the top 15 DEGs can distinguish three groups. According to the KEGG enrichment, there were 8, 7, and 7 DEGs enriched in cancer related pathway, PI3K-Akt signaling pathway, and microRNAs, respectively. In the tumor-stroma crosstalk model, significant differences can be validated in the AKAP12, CLDN1, CP, FKBP1A, LAMB3, LSM4, MTMR3, PRKARIA, YWHAZ, and JUND expressions. Conclusions These results identified the potential core genes and pathways in pancreatic cancer for tumor-stroma crosstalk, which could provide potential targets for the treatment of pancreatic cancer.


Background
Accompanied with nearly 100% of 5-year mortality rate, pancreatic cancer is one of the most quickly fatal cancers around the world [1]. Although in recent year we have some amazing improvements in the surgery, radiation therapy, and chemotherapy, pancreatic cancer still has a desperate prognosis [2]. It is one of the main causes for clinical treatment difficulties that pathogenesis and development of pancreatic cancer are not fully understood [3]. Thus, an in-depth exploration into the molecular mechanism of pancreatic cancer biology is urgently needed to develop effective therapeutic approaches.
Cancer is not only actuated by the accumulation of variety of somatic aberrations, but also accelerated by the interaction between cancer cells and the ambient microenvironment [4]. The tumor microenvironment consists of a variety of cell types, such as immune cells, pericytes, fibroblasts, bone-marrow-derived cells, and vascular endothelial cells, embedded in the extracellular matrix (ECM). In recent years, the opinion that stromal cells contribute a great effort to tumor initiation and progression was extensively accepted [5]. Cancer-associated fibroblasts (CAFs) can induce the tumorigenesis through ECM remodeling, angiogenesis, and the secretion of soluble factors.
Remarkable desmoplasia is the pathological feature of pancreatic cancer and leads to its malignant potential. Desmoplasia includes an excessive amount of ECM, which inhibits drug delivery to tumor cells, resulting in chemoresistance [6]. Now, several therapeutic agents have been developed to decrease excessive ECM, such as ECM protein with inhibitor of hyaluronic acid (HA), pegylated recombinant human hyaluronidase (PEGPH20), a novel agent that degrades HA to enhance the delivery of cytotoxic agents, which has demonstrated promising preclinical results and early clinical evidence of efficacy in the first-line treatment 2 BioMed Research International of metastatic PDCA with acceptable tolerability. However, for existing therapeutic agents, the potential to augment therapies remains rational next steps and more results of upcoming clinical trials will provide critical guidance [7,8].
As a main member of CAFs, pancreatic stellate cells (PSCs) are the main sources of ECM production in desmoplasia of pancreatic cancer. In the normal pancreas, PSCs are located between pancreatic lobules and acinar [7]. Under pathological conditions, resident PSCs are activated and secrete excessive amounts of ECM proteins, leading to desmoplasia. There are profound interactions between PSCs and pancreatic cancer cells [8]. Pancreatic cancer cells can stimulate PSCs to produce excessive ECM proteins and cytokines [9]. On the other hand, PSCs also can enhance pancreatic cancer cell proliferation, motility, invasion, and chemoresistance [10]. Based on the evidences above, it is suggested that the tumor-stroma crosstalk of pancreatic cancer is a complicated process. Despite the progress that has been done, the molecular mechanisms of tumor-stroma crosstalk remain unclear. In this study, we used bioinformatics methods to analyze the mRNA expression data during coexposure of pancreatic tumor and PSCs, to identify potential core genes and pathways in tumor-stroma crosstalk, aiming to provide valuable information for further pathogenesis mechanism elucidation and supply the groundwork for therapeutic target identification for pancreatic cancer.

Materials and Methods
. . Affymetrix Microarray Data. The source of our data, GSE49583, GSE49584, and GSE49586 transcriptional profile were provided by Giese NA et al. They were in the GEO database (http://www.ncbi.nlm.nih.gov/geo/) from the National Center for Biotechnology Information (NCBI).

. . Identification of Differentially Expressed Genes (DEGs).
The raw data of the mRNA expression profiles were downloaded and analyzed by R language software. Background correction, probe summarization, and quartile data normalization were applied to the original data. The limma method [11] in Bioconductor (http://www.bioconductor.org/) was utilized to identify DEGs; the significance of DEGs was calculated by t-test and was represented by P value. To reduce the risk of false positives, the Benjamini-Hochberg False Discovery Rate (FDR) method was used for adjusted P values for multiple testing. The corrected P value was represented by FDR. The | log 2 FC| > 1.5 and FDR < 0.01 were used as the cut-off criteria. We were interested in the intersection of these sets.
. . Hierarchical Clustering Analysis of Expression Profiling in PDAC of Tumor-Stroma Crosstalk. To reveal samples in which the most similar groups were neighbors, a two-way hierarchical clustering analysis [12] was applied to genes using the "heatmap" package in the R language. The results were displayed using a heat map.
The GO-BP terms, GO-CC terms, and GO-MF terms were filtrated by the standard of P value smaller than 0.01 and the significant enriched KEGG pathways were discerned with a P value smaller than 0.05. (TFs). FunRich (Functional Enrichment analysis tool) from ExoCarta (http://www.exocarta.org/) was used to perform this analysis. As for TFs, the DEGs were chosen and annotated to judge the related TFs. And the cut-off standard was P<0.01.

. . Principal Component Analysis of the Top
DEGs. To analyze the tumor-stroma crosstalk in pancreatic cancer, 15 DEGs with the smallest significant P value were chosen to implement the principal component analysis through R language package pca2d. Principal component analysis is up to determine the major variants in a multidimensional data series [16].
. . PPI Network Construction and Identification. The context involving molecular mechanism of cellular processing is provided by the functional interactions between proteins. In the current study, DEGs protein-protein interaction (PPI) network was constructed by the Search Tool from the Retrieval of Interacting Genes (STRING, http://string.embl.de/) database and then was visualized through Cytoscape [17].
. . Isolation and Identification of PSCs and the Construction of Tumor-Stroma Crosstalk Model. Healthy male C57BL/6 mice weighing 13 to 25 g and aged 3 to 11 weeks were used in our tests. Normal C57BL/6 mouse PSCs were primary isolated and cultured using a published method [9]. In detail, pancreatic tissues in C57BL/6 mice were washed by fetal bovine serum, minced into 0.4-1.0 mm3, and then cultured in sterile culture flasks. When they reached 70%-85% confluence after 4-5 days, primary PSCs were collected and then passed on. Primary mouse PSCs were used in the second to sixth passage of this study.
The changes in intracellular lipid droplets were detected using oil red O staining (Sigma-Aldrich); -SMA and desmin expression were tested by immunocytofluorescent staining. Intracellular lipid droplets, -SMA, and desmin expression were used to identify primary PSCs. The list of primary antibodies was shown in Supplementary Table S1 online.
Isolated primary C57BL/6 mice PSCs were treated by DMEM with 15%FBS containing the mixture of Aspc-1 and Panc-1 supernatant for 72 hours as the tumor-stroma crosstalk model.
. . Western Blotting Analysis. Western blotting analysis was performed as we described above [9]. Protein concentrations were measured by the bicinchoninic acid (BCA) Protein Assay Kit (P0010, Beyotime Biotechnology, China). Samples of 15 g total protein were used for western blotting. The list of primary antibodies was shown in Supplementary Table S1 online.
. . e Top DEGs Traits in Human Pancreatic Cancer Tissue. The top 15 DEGs protein expressions in pancreatic cancer tissues were determined from the human protein atlas (www.proteinatlas.org).
. . Statistical Analysis. Data are shown as means ± SD. All statistical analyses were performed via SPSS 16.0 for Windows (SPSS Inc., IL, USA). A two-tailed Student's t-test or oneway ANOVA was used to contrast intergroup variance. The P value of < 0.05 was considered statistically important.

. . Identification of DEGs and Hierarchical Clustering Analysis of Expression Profiling in Tumor-Stroma
Crosstalk. Total 28 samples are the expression data from primary human PSCs treated with 8 tumor cell lines, MiaPaCa2 cell treated with human naïve PSCs and MiaPaCa2 treated with human prestimulated PSCs, respectively. The samples (named from GSM1202262 to GSM1202279 and from GSM1202282 to GSM1202291) were found gathered in pancreatic cancer in tumor-stroma crosstalk clusters (Supplementary Figure S1 online). 28 samples were grouped into three clusters by hierarchical clustering. There were 221 DEGs in the intersection. The top 15 DEGs with the minimum significant P value genes were listed in Table 1.

. . GO Analysis of Identified DEGs.
To gain insight into the functional characteristics of identified 221 DEGs, we conducted GO enrichment analyses. As shown in Table 2 . . KEGG Pathway Enrichment of Identified DEGs. According to the KEGG pathway-related database, there were 8, 7, 7, 6, 6, 5, 5, 4, and 4 DEGs enriched in pathways in cancer, PI3K-Akt signaling pathway, microRNAs in cancer, lysosome, IL-17 signaling pathway, pancreatic cancer, MAPK signaling pathway, P53 signaling pathway, and TNF signaling pathway, respectively. Meanwhile, the results of pathway enrichment analysis via DAIVD were shown in Table 3.     . . PPI Network Construction. PPI network analysis can identify the key hub members among a cluster of molecules. Based on STRING database, a PPI network of the top 15 DEGs was constructed. The network consisted of 75 nodes and 675 edges. Average node degree was 18. Clustering coefficient was 0.772. As shown in Figure 1, no interaction of CP, LAMB3, CLDN1, ZBTB33, MTMR3, S100Z, and KLRC3 with other proteins was showed, while the rest 8 DEGs constituted 60 PPI pairs. PPI enrichment P value was less than 0.01.
A PPI network of 221 DEGs was also constructed. The network consisted of 274 nodes and 1261 edges. Average node degree was 9.2. Clustering coefficient was 0.736. PPI enrichment P value was 6.76E-12 (Supplementary Figure S4 online).
. . Characterization of PSCs. For 5 days, the primary PSCs cells crawled out. The cells contained lipid droplets by oil red O staining and stained positively for desmin but negatively for -SMA by IF. When passed on, PSCs were auto-activated, becoming myofibroblast-like cells. Cells were positive for -SMA by IF (Figure 2).

. . Validation of the Top
DEGs Expressions in Tumor-Stroma Crosstalk. After PSCs were treated with the mixture of Aspc-1 and Panc-1 supernatant for 72 hours, the mRNA of PSCs were extracted to validate the top 15 DEGs expressions. The mRNA levels of the top 15 DEGs were shown in Figure 3. CLDN1, CP, FKBP1A, LAMB3, LSM4, MTMR3, YWHAZ, and JUND mRNA levels were much higher in the tumor-stroma crosstalk. Compared with the control group, PRKARIA and AKAP12 were lower in the tumor-stroma crosstalk. There were no significant differences in SAR1B, ZBTB33, STYX, KLRC3, and S100Z mRNA levels. The significant genes were also identified by western blotting analysis. Our result showed that mRNA and protein data were consistent (Figure 4). the top 15 proteins expression in clinical specimens from the human protein atlas (www.proteinatlas.org). Except for S100Z, the other 14 DEGs expression profiling in pancreatic cancer specimens was shown in Figure 5. No data of S100Z were reported in human atlas.

Discussion
Despite advances in medical and surgical therapy, pancreatic cancer still has a poor prognosis. An in-depth exploration into the malignant essence of pancreatic cancer is urgently needed to improve survival rate. There has been increasingly accumulating evidences that support substantial two-way interactions between the stromal components and cancer cells [10,19,20]. As a member of stromal cells, PSC plays a vital role in pancreatic fibrosis correlated with pancreatic cancer. Our study identified the potential core genes and pathways in pancreatic cancer for tumor-stroma crosstalk, which is pivotal to the development of new treatment strategies.
Cancer is a complex structure including malignant cells and multifarious surrounding cells. Remarkable desmoplasia is the pathological feature of pancreatic cancer. PSCs are the main sources of ECM production in desmoplasia of pancreatic cancer. Evidence is emerging that there exists a symbiotic relationship between pancreatic cancer cells and PSCs, which lead to an overall increase in the rate of growth of the tumor [21]. Pancreatic cancer cells and PSCs mutually promote each other's differentiation and proliferation [22].
Lately, microarray technology has been widely applied to reveal the genetic alteration in the development of diseases. So, in this study, we used bioinformatics methods to analyze the mRNA expression data during coexposure of pancreatic tumor and PSCs, which were available on the GEO database. In order to gain insight into the functional characteristics of the identified DEGs, we conducted GO, KEGG pathway enrichment analyses, PPI network analysis, and principal component analysis to analyze key genes and pathways in the tumor-stroma crosstalk. Then, we validated the top 15 DEGs expression in the tumor-stroma crosstalk model. Apart from SAR1B, ZBTB33, STYX, KLRC3, and S100Z, there were significant differences in the other 10 DEGs levels in the tumor-stroma crosstalk. It is supported that they may play a key role in pancreatic cancer for tumor-stroma crosstalk.
CLDN1 is shown to function as a tumor suppressor or oncogene based on some specific cellular context [23]. In pancreatic cancer, CLDN1 can lead to TNF-alpha-dependent proliferation. CLDN1 is also correlated to EMT of pancreatic cancer [24]. It is implied that CLDN1 immunophenotype is closely relevant to the malignant behavior of pancreatic cancer.
Serum levels of Ceruloplasmin (CP) were increased in different tumors. The Capan-1 cells produce CP which contains sLex. Pancreatic cancer cells can result in the synthesis of increased sLex on CP found in patients with pancreatic cancer [25]. Since CP is involved in angiogenesis and neovascularisation of cancer, CP should be a key gene in tumor-stroma crosstalk of pancreatic cancer.
PRKAR1A is a tumor suppressor in the pancreas and points to the PKA pathway as a possible therapeutic target for these lesions. Loss of this gene leads to nonfunctional neuroendocrine tumors with an acinar component. Since loss of PRKAR1A is sufficient to cause endocrine (as well as exocrine) tumorigenesis in the pancreas, PRKAR1A may play a key role in tumor-stroma crosstalk.
LAMB3 can be found in many different kinds of epithelial tissues and tumor microenvironment. It has been reported that LAMB3 may play a significant role in the prognosis and progression of pancreatic cancer. Over-expression of LAMB3 is correlated to clinicopathologic features and reduced survival in pancreatic adenocarcinoma patients [26].
In conclusion, we selected the DEGs and explored the underlying molecular mechanism in pancreatic cancer for tumor-stroma crosstalk by bioinformatics methods. The top 15 DEGs may play a significant role during pancreatic cancer pathogenesis and development. These results will provide valuable information for further molecular mechanism elucidation of tumor-stroma crosstalk in pancreatic cancer. Nevertheless, further experiments are still required.