Preserved correlation matrices pinpoint extracellular matrix organization as a critical factor in pancreatic ductal adenocarcinoma

Introduction

With one of the lowest survival rates among any cancer, pancreatic cancer (PC) is a gastrointestinal malignancy with exceptional lethality. Despite its relatively low incidence rates, PC is the seventh leading cause of cancer-related death worldwide and the third in the USA after lung and colorectal cancer. The incidence rates are higher in more developed countries¹ and display a rising tendency. By 2030, pancreatic cancer is predicted to be the second leading cause of cancer-related mortality in the USA.² Approximately 90% of pancreatic cancer cases belong to pancreatic ductal adenocarcinomas (PDAC), arising from precursor lesions termed pancreatic intraepithelial neoplasia. The cause of the disease is multifactorial, but cigarette smoking appears to be one of the primary risk factors. However, alcohol consumption, obesity, aging, diabetes mellitus, and family history are significant contributors as well.^3–5

Unfortunately, most patients seldom exhibit symptoms at the early stages of the disease. The wide range of non-specific symptoms, including jaundice, back- and abdominal pain, and weight loss, hampers precise diagnosis. In addition, existing screening biomarkers, including serum levels of CA19-9, display relatively low sensitivity and specificity.⁶ As a result, 80-85% of patients are diagnosed with unresectable tumors characterized by regional spread or metastatic disease. Only 24% of diagnosed patients survive their first year, while five-year survival is as low as 10%. Currently, surgery is the only available potentially curative option.⁷ Based on the extent of the disease, patients are allocated into one of four categories: resectable, borderline resectable, locally advanced, and metastatic. Altogether, only 15-20% of patients are good candidates for pancreatectomy, and the five-year survival after surgery remains low at 25%. Adjuvant therapy after surgery, including gemcitabine plus capecitabine, or modified FOLFIRINOX (oxaliplatin, irinotecan, leucovorin, and fluorouracil), is the preferred option, as both treatment regimens improve overall survival.^8,9 The combination therapies based on modified FOLFIRINOX, or gemcitabine with nab-paclitaxel for metastatic pancreatic cancer, became the gold standard.^10–12 In the recent PRODIGE24-trial, a modified FOLFIRINOX regimen resulted in prolonged disease-free (21.6 vs. 12.8 months) and overall survival (54.4 vs. 35.0 months) compared with gemcitabine monotherapy. Nearly 20% of FOLFIRINOX treated patients were alive at 18 months compared to only 6% of those treated with gemcitabine alone.⁹ Locally advanced disease with vascular involvement is regularly treated with systemic chemotherapies approved in the metastatic setting. Despite standard chemotherapy, the vast majority of resected patients experience recurrence and die within a short period.

Despite extensive molecular characterization, personalized treatment strategies of PDCA are still limited. About 5-9% of pancreatic cancer patients carry germline BRCA1 or BRCA2 mutations. In 2019, the FDA approved olaparib, a PARP-inhibitor, as maintenance therapy for metastatic patients with BRCA1/BRCA2 mutations who had not progressed during first-line platinum-based chemotherapy. Olaparib therapy improved progression-free survival in this setting compared to placebo and became the first approved targeted therapy for pancreatic cancer.¹³ Other approaches, including the adoption of immune-checkpoint inhibitors, did not yet translate into clinical benefits, as pancreatic tumors are characterized by highly immuno-suppressive microenvironments.¹⁴ The diverse response to treatment sparked initiatives for patient classifications, but classical histopathological features lack clinical utility.¹⁵ Based on the further genome-wide screening of transcriptional changes, PDCA has been classified into two,¹⁶ three,¹⁷ or four¹⁸ distinct molecular subtypes. Cross-comparison of these studies revealed reproducible biologic subgroups with some prognostic relevance, but the described molecular subtypes have not yet been validated or translated into successful therapeutic decisions.¹⁹ Nevertheless, principal drivers of pancreatic tumorigenesis have been identified, including KRAS, SMAD4, CDKN2A, and TP53.^1,20,21 Because so few patients can have a successful surgery, and the number of therapies is limited, identifying novel potential biomarkers for early diagnosis and targeted interventions are of utmost importance.

In this study, our aim was to exploit gene expression alterations between normal and tumor tissue samples. However, instead of searching for genes upregulated in tumors, we approached biomarker selection from a novel perspective. In normal, diploid tissue samples, gene co-expression correlations frequently signal shared biological functions with coordinated regulation. On the other hand, malignant tumors frequently undergo abnormal genetic transformations, but some correlations remain consistent across normal tissues and tumor samples. We speculated that maintained correlations might be vital for cellular proliferation and survival; thus, identifying preserved correlations in gene expression across normal and tumor tissue samples can provide new targets for developing alternate therapies.

Methods

Results

The combined microarray (GEO) dataset consists of 248 pancreatic tumors and 108 normal tissue samples, allowing the analysis of 12,210 genes. The RNA-seq dataset incorporates 177 tumors and four normal samples from TCGA, and 248 normal tissue samples from GTEx, enabling the analysis of 21,479 genes.

Each dataset underwent multistep filtering: gene expression differences were assessed between tumor and normal tissue samples based on log2FC values. Genes with log2FC > 0.5 or log2FC < -0.5 from the GEO dataset, and genes with log2FC > 1 or log2FC < -1 from the RNA-seq dataset were incorporated into subsequent analyses. Altogether 6088 genes had altered expression at p < 0.01 across tumor and normal tissue samples in the GEO dataset, while 17,967 genes at p < 0.01 in the RNA-seq dataset. Genes with median expression below the cutoff value of 500 were excluded from subsequent analyses to reduce the chance of noise discoveries; therefore, the number of genes was reduced to n = 2,205 in the GEO dataset and n = 3,773 in the RNA-seq dataset.

We identified 309 common genes in both RNA-seq and microarray datasets with altered gene expression across normal and tumor tissues (Figure 1A). In general, more overexpressed genes were identified in tumors compared to normal tissues. Out of the 309 genes, 267 had increased expression in primary tumors compared to 42 genes overexpressed in normal samples, both in GEO and RNA-seq datasets. Genes with altered expression across normal and tumor tissues are graphically represented by a VOLCANO plot (Figure 1B).

Figure 1. The process of data selection from the TCGA and GEO databases (A). Genes with altered expression across tumors and normal tissue samples, in RNA-seq (green) and gene chip (red) datasets (B). Representation of the gene ontology analysis, including the 262 genes participating in maintained correlations. The size of the dot indicates the number of genes associated with a particular biological process (C).

Gene overrepresentations

Seventy-nine genes from the identified 309 differentially expressed genes participated in preserved correlations, out of which 73 genes were upregulated in tumors compared to normal tissue samples (Supplemental Table 2).

Some genes overexpressed in tumor samples were strikingly overrepresented among the 371 preserved correlations, including SPARC, COL6A3, MMP2, HTRA1, FN1, PALLD, and COL3A1, some participating in as many as 58 correlations (Figure 2A). Over 70% of all preserved correlations (263 altogether) involved at least one of these seven genes.

Figure 2. Seven genes were heavily overrepresented in maintained correlations, among those SPARC participating in as many as 58 correlations (A). Increased expression of 28 genes from preserved correlations was associated with worse overall survival in pancreatic cancer patients, including ENO1, DSG2, CLIC1, ADAM10, ACTR3, FN1, MYL12A, and COL4A2 (B). A STRING protein network analysis including the 79 genes with differential expression from preserved correlations. FN1 has a central role in the associated network of proteins (C).

Gene expression alterations during tumor progression

The 371 preserved correlations incorporated 79 of the initially identified differentially expressed genes, out of which the majority, 73 genes, were overexpressed in tumors compared to normal tissues. With the help of TNMplot.com, we have also compared the expression pattern of these genes between primary tumors and metastases. The expression of MYL12A increased most significantly in metastatic samples compared to primary tumors (Figure 3A). Of note, the expression of both MYL12A and MYL12B gradually intensified across normal samples, primary tumors, and metastatic tissue samples (MYL12A, KW p = 2.57E-11; MYL12B KW p = 1.38E-14). We have observed a similar pattern in seven additional genes, including ENO1 and CLIC1, that were previously identified as robust prognostic biomarkers of survival (Figure 3B).

Figure 3. Comparing gene expression across tumors, normal and metastatic samples.

Expression differences are color-coded: red encodes increasing expression from normal to tumor and from tumors to metastatic samples, while green encodes decreasing expression (A). The expression of MYL12A, MYL12B, CLIC1 and ENO1 genes gradually increased from normal tissue samples to primary tumors and metastatic samples (B).

The close-knit network of 262 genes involved in preserved correlations is mainly associated with ECM organization and function (C). Gene ontology analysis involving the 79 differentially expressed genes from preserved correlations revealed similar top functions (D).

Discussion

Pancreatic cancer is one of the most lethal malignancies with limited options of therapies and biomarkers. In order to pinpoint survival-specific vulnerabilities, we aimed to identify upregulated genes with preserved correlations across normal and tumor samples. Previous research revealed profoundly different biomolecular networks underlying cancer conditions compared to healthy tissues.³² Notwithstanding, we speculated that sustained correlations might offer insight into fundamental biological functions, providing potential novel targets for therapeutic interventions. To increase robustness, we collected data from two independent platforms, including RNA-seq and microarray-based studies. Then, we identified genes with an altered expression between tumors and normal tissues and looked for correlations preserved across both data platforms.

We have identified 371 preserved correlations comprising 262 genes, including 79 genes with altered expression across tumors and normal samples. These preserved correlations were robust, as were independently identified in both microarray and RNA-seq datasets. Moreover, a few of these differentially expressed genes remarkably dominated preserved correlations, including SPARC, COL6A3, COL3A1, MMP2, HTRA1, FN1, and PALLD. Altogether, 16 genes upregulated in tumor tissues were linked to dismal disease outcomes, out of which the expression of MYL12A gradually increased during disease progression.

Upon close examination, most of these genes are associated with the tumor microenvironment, participating in extracellular matrix (ECM) organization. PDAC has a dense fibrotic stroma, composed of abundant extracellular matrix and other non-neoplastic cell types, responsible for chemo- and radiotherapeutic resistance.³³ The abundant ECM deposition, called desmoplastic reaction, plays a crucial role in PDAC. Beyond structural support, the ECM functions as a reservoir of growth factors and a mediator of complex signaling pathways.^34,35 The loss of ECM homeostasis and integrity is one of the hallmarks of cancer, resulting in cancer progression³⁶; thus, potential targeting and modifications of the ECM in PDAC are being intensely studied both in the preclinical and clinical settings.^37,38

Our results corroborate the importance of ECM organization in disease outcomes. Most importantly, we pinpoint critical interactions as potential targets for combination therapies. Fibronectin (FN1) is a principal constituent of the extracellular matrix that provides the scaffold to build up the entire ECM structure, affecting its architecture and composition. It is overexpressed in PDAC as part of the desmoplastic reaction, is associated with aggressive tumor characteristics,³⁹ and functions as a primary driver of tumor progression by enhancing growth factor signaling.⁴⁰ FN1 secreted by pancreatic stromal cells promoted gemcitabine resistance by downregulation of ERK1/2 phosphorylation.⁴¹ A synthetic FN1 blocking-peptide (Arg-Gly-Asp-Ser) abrogated the gemcitabine resistance in PDAC cell lines.⁴¹

Our datasets identified both FN1 and the matricellular protein SPARC remarkably overrepresented in preserved correlations across normal and tumor samples. High stromal expression of SPARC was linked to poor patient prognosis in PDAC,⁴² and pancreatic stellate–cell-derived SPARC promoted proliferation of PDAC cells in vitro.⁴³ Recent research indicated that the presence of plasma fibronectin determines the specific effect of SPARC on pancreatic cancer cells: depletion of fibronectin switched SPARC from promoting cancer cell proliferation to growth inhibition and induction of apoptosis.⁴³ In our data both FN1 and SPARC expression were correlated with MMP2 in normal and cancer tissue samples. MMP2 is a member of the matrix metalloproteinase (MMP) gene family, primarily located in the tumor stroma, capable of cleaving collagens and other components of the extracellular matrix involved in signal transduction.⁴⁴ In PDAC, MMP2 is among the most studied gelatinases, with a potential clinical relevance⁴⁵: MMP2 expression is linked to the activity and the invasive potential of pancreatic cancer cell lines, also associated with the development of the characteristic desmoplastic reaction.⁴⁶ The potential importance of MMP2 in PDAC progression is supported by studies using selective MMP2, MMP9, and MMP14 inhibitors in mice and Syrian hamsters, such as MMI-166, RO28-2653, OPB-3206, and SB-3CT.⁴⁵

The most abundant component of the mammalian ECM is collagen. Collagen Type IV Alpha 2 Chain (COL4A2) encodes one of the six subunits of collagen type IV, a major structural component of basal membranes.⁴⁷ In our dataset, COL4A2 expression was highly correlated with COL4A1, another type IV collagen alpha protein found in most connective tissues. In pancreatic cell lines and xenograft models of pancreatic cancer, type IV collagen production exceeded that of other ECM proteins examined, including collagen type I.⁴⁸ The abundantly produced collagen type IV is deposited near the pancreatic tumor cells, creating a basal membrane-like structure.⁴⁹ The structure colocalizes with integrin receptors on the cancer cell surface, stimulating the proliferation and migration of pancreatic cancer cells. Thus, type IV collagens provide essential survival signals to the pancreatic cancer cells through an autocrine loop.⁴⁹ High circulating level of postoperative collagen IV has been linked to dramatically reduced survival in PDAC patients undergoing curative surgery.⁵⁰

Substantiating previous findings, we found an association between high COL4A2 expression and poor prognosis, supporting the role of COL4A2 as a potential prognostic biomarker of PDAC. The drug candidate T12, a rationally designed inhibitor of collagen IV network formation, reduced both COL4A1 and COL4A2 expression and disrupted EMT-based chemo-resistance in mouse models of breast and lung cancers without apparent toxicity.⁵¹ Thus, targeting specific EMC proteins, including COL4A2 and COL4A1, may represent a novel therapeutic opportunity for treating PDAC.

Another protease, the ADAM Metallopeptidase Domain 10 (ADAM10), has also been heavily investigated by pharmaceutical companies. ADAM10 is ubiquitously expressed and cleaves many substrates involved in cancer signaling, such as PD-L1, NOTCH, EGFR/HER-ligands, among others.⁵² The evidence for ADAM10 as a suitable cancer biomarker increases: ADAM10 overexpression correlates with poor patient prognosis in cancers of the colon, lung, ovary, uterus, and stomach.⁵³ Moreover, it has a substantial value as a prognostic biomarker in glioblastomas, sacral chordomas, and triple-negative breast cancer and became a potential target.^52,54,55 In PDAC, ADAM10 was suggested to be a driver of tumor progression and fibrosis after radiation therapy. Pharmacological inhibition of ADAM10 delayed fibrosis and improved survival in radiation therapy-treated orthotopic and metastatic mouse models.⁵⁶ Targeting the active ADAM10 with a monoclonal antibody, such as mAb 8C7, that masks the substrate-binding pocket preferentially located in tumors compared to normal tissue may offer targeted inhibition to prevent tumor growth.⁵⁷

Our results suggest additional insight into the evolving tumor microenvironment during disease progression. The expression of MYL12A was associated with poor prognosis and correlated with MYL12B upregulation within both normal and tumor samples. The expression of both genes increased during tumor progression, from normal tissues to primary tumors to metastatic samples. Myosin Light Chain 12A (MYL12A) and 12B (MYL12B) regulate smooth muscle and non-muscle cell contractile activity. Both genes are implicated in cytokinesis and cell locomotion,⁵⁸ and intensifying expression levels indicate the acquisition of motility by cancer cells. The biological functions and rising expressions during cancer progression warrant the role of MYL12A both as a potential prognostic biomarker and a therapeutic target of PDAC.

We observed an increase in the expression of ENO1 and CLIC1 genes between normal samples and tumor tissue. Moreover, in our data, overexpression of both ENO1 and CLIC1 was associated with poor overall survival. Alpha enolase (ENO1) is a glycolytic enzyme localized on the cell surface with a complex role in cancer progression: ENO1 promotes glycolysis and ECM degradation, affects actin remodeling, facilitates tumor invasion, and migration, and activates cancer signaling pathways.⁵⁹ During inflammation, enolase migrates to the cell surface that promotes the production of pro-inflammatory cytokines, chemokines, reactive oxygen species, and nitric oxide.⁶⁰

Higher ENO1 expression has been associated with larger tumor size, increased invasion, shorter survival, and poor prognosis in several solid tumors, including hepatocellular carcinomas, breast- and lung cancers,⁵⁹ and PDAC.⁶¹ ENO1 has clinicopathological and diagnostic significance as a tumor-associated antigen as well: combining ENO1 autoantibodies with other biomarkers improved diagnostic sensitivity and accuracy of lung cancer detection.^62,63 In PDAC, several fundamental characteristics designate ENO1 to be an ideal biomarker candidate: it is localized on a cell surface, making it targetable; it is upregulated in tumors compared to normal tissues, is associated with clinical outcome, and has vital metabolic functions. Small molecule inhibitors (e.g., ENOblock, ENO1i) are already being investigated in various preclinical settings and disease models, including spinal cord injuries and multiple myelomas.^64,65 Additionally, DNA vaccinations targeting ENO1 in combination with chemotherapy may offer promising alternatives in PDAC treatment. In pancreatic cancer-prone mice, treatment with gemcitabine combined with anti-ENO1 vaccine elicited a robust antitumor response and reduced the size of PDAC lesions.⁶⁶ Administration of an anti-ENO1 mAb inhibited in vitro invasion of human PDAC cells and their metastatic spread in immunocompromised mice; thus, targeting ENO1 may be exploited as a novel therapy to increase the survival of metastatic PDAC patients.⁶⁷

PDAC presents unique challenges for chemotherapy and immunotherapy with an immunosuppressive microenvironment and a dense stroma, creating a physical barrier to drug delivery; thus, alternative therapeutic targets, such as ion channels, provide an increasingly evaluated option. Ion channels are involved in every aspect of oncogenesis, influencing critical processes within each hallmark of cancer. Chloride Intracellular Channel 1 (CLIC1) plays a vital role in many physiological functions, including ion homeostasis, acidity, and cellular volume, considered to be both the sensor and executor of oxidative stress. It can exist in a soluble form, but in pathological states, CLIC1 is expressed explicitly as a transmembrane chloride channel.⁶⁸ CLIC1 overexpression has been established in pancreatic cancer tissue and is indicated to be a putative oncogene.^69–71

In summary, our findings provide novel insight into intimate gene-networks corroborating previous knowledge about interactions within the tumor microenvironment. Although our results are based on in silico approaches, maintained correlations across RNA-seq and gene chip platforms substantiate the significance of our findings. Our results may be utilized to develop combination therapies against PDAC and offer novel ideas to be explored in future studies.