The inflammatory response of human pancreatic cancer samples compared to normal controls

Pancreatic ductal adenocarcinoma (PDAC) is a poor prognosis cancer with an aggressive growth profile that is often diagnosed at late stage and that has few curative or therapeutic options. PDAC growth has been linked to alterations in the pancreas microbiome, which could include the presence of the fungus Malassezia. We used RNA-sequencing to compare 14 matched tumor and normal (tumor adjacent) pancreatic cancer samples and found Malassezia RNA in both the PDAC and normal tissues. Although the presence of Malassezia was not correlated with tumor growth, a set of immune- and inflammatory-related genes were up-regulated in the PDAC compared to the normal samples, suggesting that they are involved in tumor progression. Gene set enrichment analysis suggests that activation of the complement cascade pathway and inflammation could be involved in pro PDAC growth.


Introduction
Pancreatic cancer (PDAC) is the 9 th most common cancer in the US but is the 4 th most common cause of cancer related death (~54,000/year and ~44,000/year respectively).The median 5 year survival for stage 4 disease is 9% [1].The high death rate with respect to the prevalence rate is due to poor early detection and a lack of meaningful advancement in systemic therapeutics.The most common somatic mutations, Kirsten rat sarcoma viral oncogene (KRAS), tumor protein p53 (TP53), cyclin dependent kinase inhibitor 2 A (CDKN2A), and SMAD family member 4 (SMAD4) have been identified by whole-exome and -genome sequencing of large PDAC cohorts and form the majority of unique mutations in patients with PDAC [2].
The tumor microenvironment (TME), which represents a complex ecosystem involving interactions between immune cells, cancer cells, stromal cells, and the extracellular matrix, can support tumor proliferation, survival, and metastasis and can be highly immunosuppressive [3][4][5].
In one paper, they found in patients with a resected PDAC and a higher tumor microbial diversity did better than resected PDAC with a lower microbial diversity.They also showed that patients who carried the long term survival (LTS) microbiome signature, namely (Pseudoxanthomonas, Streptomyces, and Saccharopolyspora) lived a long time.There was no microbiome signature for the patients who had a short term survival (STS).Patients with high diversity and low diversity had overall survival of 9.66 vs. 1.66 years respectively using univariate Cox proportional hazard models.
They presented data in a mice model showing the benefit of a fecal microbiota transplantation of the LTS compared to the STS gut microbiome.Thus, tumor micro-environment is an important part of the prognosis of patients with resected PDAC [6].
On the other hand, bacterial ablation in the pancreas was shown to alter the immune system, allowing the immune system response to decrease growth of the tumor by a reduction in myeloid-derived suppressor cells and an increase in M1 macrophage differentiation, promoting Th1 differentiation of CD4+ T cells and CD8+ T cell activation.Bacterial ablation also enabled efficacy for checkpoint-targeted immunotherapy by up-regulating programmed death (PD-1) expression [7].More recently the same group looked at the role of the fungal microbiome in PDAC [8].In a mouse model they showed a significant migration of fungi from the lumen of the gut into the lumen of the pancreas in both mice and humans with PDAC compared to non-malignant controls.Both in mice and humans, the Malassezia fungus was highly enriched in the malignant pancreatic tumor.They demonstrated that the binding of mannosebinding lectin (MBL) to glycans of the fungal wall, and lectin pathway activation was required for oncogenic progression in pancreatic cancer.
Furthermore, deletion of MBL or C3 in the extra-tumoral compartment or knockdown of C3ar in tumor cells were both protective against tumor growth in mouse models.There was thus a suggestion that using anti-fungal agents at some point of the treatment of pancreatic cancer may result in a shrinking or slowing down the growth of the tumor.In a paper reviewing expression of complement in various cancers most cancers over-expressed C3 but neutrally-expressed C5, and the most over expressed gene was CD59, suggesting efficient protection of malignant cells from complement-mediated killing [9].
To document the expression of C3 and Malassezia and associated genes in our patient population, we decided to perform studies on paraffin embedded pancreatic cancer biopsies using pancreatic cancer and their normal tissue counterpart samples from the University of New Mexico (UNM) Tissue Repository.We obtained exemption to perform the study from the local Human Research Protections Program.

Microbiome results
Using optimized methods [10], we used RNA-seq analysis on matched tumor-normal PDAC tissue samples derived from formalin-fixed paraffin embedded (FFPE) slices for 15 patients obtained from UNM Tissue Repository, generating high quality data for 14 patients, with an average of 32 x 10^-6 reads per sample (Tables 1 and 2).After sequencing, reads were taxonomically classified using kraken2 [11][12][13] and braken [14] against a library containing human, viral, bacterial, and fungal genomes, including the genomes of several Malassezia species.The fungus Malassezia was present in all tumor samples at varying concentrations (Tables 1 and 2).We also detected Malassezia in all of the normal tissue samples, which is in contrast to a previous study (Tables 1 and 2) [8].Despite reports that PDAC tumors have a lower microbial diversity than normal pancreatic tissues [6], we found no significant difference in the number of microbial genera present in normal samples versus tumor samples (Table 1).

Transcriptome analysis
Unsupervised multidimensional scaling (MDS) was used to compare the normal and tumor samples.The MDS plot (Fig 1A) revealed two tumor sample outliers, which were more similar to the normal samples than the other tumor samples.In addition, most of the normals clustered tightly together, while the tumor samples were more dispersed, indicating that the tumor gene expression pattern were more heterogeneous while the normals were more homogenous.Differential gene expression (DGE) analysis of the 14 matched tumor-normal patient samples  Malassezia is proposed to promote tumor growth via the complement pathway, disruption of which activates inflammation [8,15].Gene set enrichment analysis (GSEA) can reveal pathway and disease phenotypes within sets of differentially expressed genes.GSEA of the set of genes differentially expressed between tumor and normal samples against the curated (C2), ontology (C5), and hallmark (H) gene sets at the molecular signatures database [16,17], revealed enrichment for genes involved in the complement cascade, complement activation, and inflammatory response.Additionally, our differentially expressed genes are enriched for expected gene sets such as pancreatic cancer, epithelial mesenchymal transition, and KRAS signaling (Table 3).
There is an extensive literature of the linkage between the complement/inflammation genes and cancer, including PDAC, focusing primarily on the tumor microenvironment and how immune cells can affect tumor growth through the inflammatory process [18,19].Since each of these pathways were significantly enriched in our dataset (Table 3), we specifically examined which genes were differentially expressed, and their directionality.As expected, only a subset of the genes in each gene set were significantly differentially regulated: We were particularly interested in the expression pattern of genes known to play a role in pancreatic cancer, complement, inflammation and EMT.Using the genesets described above, as well as a pancreatic cancer genesets (Gruetzmann pancreatic cancer up and Gruetzmann pancreatic cancer down [20], we looked for genes involve in multiple processes that indicate  networks of influential genes using a CNET plot (Fig 4).This method showed that there were very few genes that overlapped the down regulated pancreatic cancer genes list and the complement pathway (1 gene), inflammatory response (1 gene), and EMT (3 genes) lists (Fig 4).However, we observed a number of genes that were up-regulated in pancreatic cancer and involved in both EMT (229 genes) and the complement pathways (14 genes), but a relative small number of inflammatory response genes (6 genes), all of which were up-regulated in our tumor samples.Interesting, PLAUR (urokinase plasminogen activator receptor) overlapped with all of these pathways, while FN1 (fibronectin 1) overlapped with both complement and EMT processes.
The upregulation and secretion of interleukin-33 (IL-33) in KRAS G12D PDAC subtype cells, and the interaction of these cells with intra-tumor fungi was reported.This resulted in the recruitment of T H 2 lymphoid cells and innate lymphoid cells (ILC2 cells) in the PDAC tumor micro-environment (TME).The intra tumor fungi identified were Malassezia globoza and A Altermata.Dectin-1 was activated by fungal components, which resulted in a dectin-1 mediated Src-Syk-CARD9 pathway activation.This in turn resulted in the secretion of IL-33, resulting in a pro-type 2 immune response and tumor progression in this PDAC subset.Reversal of IL-33 or anti-fungal treatment caused PDAC tumor regression [21].The role of intra-tumor fungus is distinct from prior observation where a fungus was shown to activate the complement system [8].
Since our differentially expressed genes were also enriched for the Hallmark KRAS signaling pathway (Table 3), we also looked for overlaps between the KRAS signaling pathway geneset, and the other three sets used above.complement receptor 2 (CR2/CD21) is a part of the complement pathway which activates C5 and the membrane attack complex (MAC).C3dg, a cleavage product of C3, can increase the complement receptor 2(CR2/CD21) binding to the B Cell Receptor (BCR) through co-ligation of C3dg/Ag complexes.CR2-mediated stimulation of peripheral B cell subpopulations demonstrated that CR2 can amplifying BCR signal transduction in the B-2 cells, but not B-1 cells.Furthermore, anti-IgM/C3dg conjugates did not stimulate B cells from CR2-deficient or CD19-deficient mice [22].
Human uPAR is a 55-60 kDa glycosylated protein that is bound to the outer lipid bilayer of the cell membrane by a glycosyl-phosphatidylinositol (GPI) anchor.Both uPA and/or uPAR allow malignant cells to dissolve the extracellular matrix (ECM) barriers and metastasize to distant sites due to the high proteolytic and migratory potential.UPAR is also believed to interact with cancer-associated intracellular signal-transduction pathways regulating other tumor related pathways, including, proliferation, survival, migration, invasion, metastasis, angiogenesis, and epithelial-mesenchymal transition (EMT).In the paper they showed a graph of uPAR gene expression profile in human cancer.UPAR expression was highly overexpressed in PDAC [23].
We also found over-expression of Galectin-3 (Gal-3) and Dectin-1 in the tumor samples compared to the normal samples (Fig 6).Galectins are a group of conserved proteins with the ability to bind β-galactosides through certain carbohydrate-recognition domains (CRD).Gal-3 is structurally unique as it contains a C-terminal CRD linked to an N-terminal protein-binding domain, being the only chimeric galectin.Gal-3 also contributes in several ways to inflammation and to the innate immune response.Gal-3 is a multifunctional protein with roles in tumor cell adhesion, proliferation, differentiation, angiogenesis, metastasis, and apoptosis [24].One group demonstrated that Gal-3 was highly expressed in human PDAC samples: 83 of 125 (66.4%), while only 32 of 108 (29%) of non PDAC pancreatic samples were positive for Gal-3 high expression.By using K-RAS mutated pancreatic cell lines genetically engineered to express high or low levels of Gal-3, they demonstrated that down regulation of Gal-3 decreased PDAC cell proliferation and invasion in vitro and reduced tumor volume in a mouse model.Conversely, transfection of Gal-3 cDNA into low-level Gal-3 PDAC cells, increased K-Ras activity.Gal-3 was found to bind to K-Ras and maintained K-Ras activity [25].In another paper, translocation of Gal-3 from the nucleus to the cytoplasm or plasma membrane led to K-Ras stabilization, and a decrease in the K-Ras suppressor miRNA let-7 (let-7).Using a Gal-3-knockout mouse embryonic fibroblasts (MEFs) model, they demonstrated that Gal-3 down regulated let-7 which resulted in increased expression of K-Ras [26].
Dectin-1 recognizes β-glucans with its carbohydrate recognition domains (CRD) and transduces signals through its immune-receptor tyrosine-based activation motif (ITAM)-like motif in the cytoplasmic domain [27].A recent study found that Dectin-1 is highly expressed on macrophages in pancreatic ductal adenocarcinoma (PDAC).Dectin-1 signaling resulted in PDAC progression, while Dectin-1 deletion or blockade of its downstream signaling decreased tumor growth.They showed that Syk phosphorylation was a downstream signal for Dectin-1 [28].Furthermore the recent paper describing the effect of the fungus on PDAC subset (Kras G12D) via an alternative pathway (the IL-33-ILC2/T H 2 axis) is activated by a dectin-1 mediated Src-Syk-CARD9 pathway [21].

Discussion
We explored the TME in the samples and found overexpression of C2, Dectin-1, and Galectin-3 represented in box plot analysis.There was no literature about the role of C2 in pancreatic cancer, but there was data showing C2 has a role in amplifying BCR signal transduction in each subpopulation of B-2 cells, but not B-1 cells, as noted above.Dectin-1 and Galectin-3 can augment PDAC progression as noted above.These two molecules can propagate the growth of the PDAC in vivo and mice models, but this must be looked at in the context of a multi molecular array, some inhibiting and some promoting PDAC growth.
We also compared the gene expression profiles of PDAC patient tumor samples to the normal sample and found several differential expressed genes.These genes were enriched for genes that are known to be involved in the complement cascade and inflammation, as well as genes involved in epithelial-mesenchyme transition, KRAS signaling, and known PDAC genes.From this we were able to identify one gene (PLAUR) that sits at the intersection of these pathways and may be a potential target for PDAC treatment.In addition, we found a handful of genes that represent additional attractive drug targets because they overlap subsets of complement, inflammation, KRAS and EMT pathways.Additionally, these results suggest that the TME contributes to PDAC progression by activation of these critical pathways.Examples of clinical correlation between PDCA prognosis and these genes include: 1) A correlation was noted between the overall survival of the pancreatic cancer patients with plasma Galectin-3 levels (n = 21).The high-galectin-3 group (n = 10) showed a poorer prognosis than the lowgalectin-3 group (n = 11) in pancreatic cancer patients (p = 0.006) [29], and 2), PLAU upregulation was associated with the basal type of PDAC, and an increased PLAU protein expression was associated with poor prognosis (P = 0044), and within the basal subtype, the clinical outcome in the high PLAU expression group is significantly worse than the low PLAU expression group (P = 0.018) [30].
We did not have the data on the clinical outcome of the patients.Thus no biological correlative analysis was performed.Unlike the Aykut B et al [8] and Alam A et al [21] groups we noted the fungus Malassezia in the pancreatic cancer and normal control in relatively the same proportions.However, despite the fact that there was no difference in Malassezia in tumors versus normal tissue, there was a significant difference in expression of a number of genes related to inflammation, infection and EMT, suggesting other causes of the differences such as a the tumor itself as an example.This was a single institution study, with only 14 samples from the Hospital Tissue Bank.It was felt that tools like DeSeq2 and EdgeR are specifically designed to analyze RNA-seq counts, even with low sample sizes.They have consistently identified differential expression and pathways in previous cases with smaller sizes with unpaired samples.Our application of these tools successfully detected differentially expressed genes even after multiple hypothesis correction.However, we acknowledge that this may not have enough power to generalize the results.Based on our data we have no way of determining if this is a reaction to a fungus or the tumor itself.In summary, we have demonstrated the interaction between proteins involved in inflammation can also have a role in propagation of growth in pancreatic cancer growth.

RNA isolation and sequencing
Total RNA was extracted from FFPE slices using the RNeasy FFPE kit (Qiagen) and the manufacture's protocol.Synthesis of cDNA and library preparation were performed using the SMARTer Universal Low Input RNA Kit for Sequencing (Clontech) and the Ion Plus Fragment Library Kit (ThermoFisher) as previously described [10,31,32].Sequencing was performed using the Ion Proton S5/XL systems (ThermoFisher) in the Analytical and Translational Genomics Shared Resource at the UNM Comprehensive Cancer Center.RNA sequencing data is available for download from the NCBI BioProject database using study accession number PRJNA940178.

Data analysis
Prior to alignment, non-human RNA-seq reads were identified and removed from analysis using the kraken2 taxonomic sequence classification system against a library containing human, fungal, bacterial and viral genomes [11][12][13] and final genera level abundances were calculated using Bracken (v2.5.0, [14].High-quality, trimmed reads classified as human, were aligned to GRCh38 (hg38) using tmap (v5.10.11).Exon counts were calculated using HTseq (v0.11.1, [33] against a BED file containing non-overlapping exons from UCSC genome hg38 and gene counts were generated by summing counts across exons.Samples were normalized for library size using DEseq2 (v1.34.0, [34]) and low expressing genes were excluded from the final analysis using a filtering threshold of 20 reads in 9 samples.Multi-dimensional scaling (MSC) was performed using plotMDS from the limma package (v3.50.0) and the base stats package was used for unsupervised heirachrical clustering [35].EdgeR (v 3.36.0 was used for the differential expression (crosswise comparison of the groups) using the glm method with an adjusted p-value (method = BH) of cutoff of 0.05 and requiring a minimum fold change of 1.5.ClusterProfiler (v 4.2.1), was used to analyze the differentially expressed genes for gene set enrichment (GSEA) against the Molecular Signatures Database and visualize differentially expressed genes in gene sets of interest [16,17,36].Default values were used at all points [37].For additional details see [10,33].All analysis was done using R version 4.1.2.

Ethics statement
We performed studies on paraffin embedded pancreatic cancer biopsies using pancreatic cancer and their normal tissue counterpart samples from the UNM Tissue Repository.These samples were anonymized and thus we obtained no clinical data from these samples.The UNM Human Research Protection Program (HRCC) reviewed our submission on 12/6/2020.They determined that informed consent and HIPAA authorization addendums were not applicable to this study.We thus obtained no verbal or written consent.

Fig 1 .
Fig 1. (A) Multidimensional scale plot of PDAC tumor (red bar) and normal (blue bar) RNA-seq data.(B) Volcano plot summarizing the differential gene expression analysis, showing log2 fold change vs. log10 of the p-value (BH adjusted).Red indicates genes that were > = 1.5 log2 fold change and had an adjusted p-value < = 0.05.(C) The heatmap summarizes the unsupervised clustering and differential gene expression analysis comparing the normal samples to the tumor samples.Dendrograms show the unsupervised hierarchical clustering of samples (top) and genes (left side).A larger version of this heatmap with samples and genes labeled is provided in the supplementary results (S1 Fig in S1 Appendix).https://doi.org/10.1371/journal.pone.0284232.g001 64 of 200 for complement (Fig 2A), 52 of 200 for inflammatory (Fig 2B), and 104 of 200 for EMT (Fig 2C).Interestingly, overall expression of C3 was higher than C2 expression on both tumors and normals, but only C2 was significantly up-regulated in tumors (Fig 3).

Fig 4 .
Fig 4. CNET plot showing gene overlaps in enriched genesets.Geneset enrichment analysis of genes differentially expressed between tumor and normal samples indicated enrichment for sets of genes up (olive) or down (orange) regulated in pancreatic cancer, the complement cascade (green), epithelial mesenchymal transition (blue), and inflammatory response (purple).Genes at the end spokes indicate differential expressed in our dataset that were enriched in each geneset, color of dots indicated whether the gene was up (red tones) or down (blue tones) regulated in tumors relative to normal samples [17, 21, 35].https://doi.org/10.1371/journal.pone.0284232.g004

Fig 5 .
Fig 5. CNET plot showing gene overlaps in enriched genesets.Geneset enrichment analysis of genes differentially expressed between tumor and normal samples indicated enrichment for sets of genes involved in KRAS signaling (purple), inflammatory response (turquoise), complement cascade (red), and epithelial mesenchymal transition (green).Black dashed circle indicates PLAUR, which connects to all four genesets.Genes at the end spokes indicate differential expressed in our dataset that were enriched in each geneset, color of dots indicated whether the gene was up (red tones) or down (blue tones) regulated in tumors relative to normal samples [17, 21, 35].https://doi.org/10.1371/journal.pone.0284232.g005

S1
Appendix.S1 Fig.The heatmap summarizes the supervised clustering and differential gene expression analysis comparing the normal samples to the tumor samples.Dendrograms show the unsupervised hierarchical clustering of samples(top) and genes (left side).S2 Fig.CNET plot showing gene overlaps in enriched genesets.Geneset enrichment analysis of genes differentially expressed between tumor andnormal samples indicated enrichment for sets of genes involved in KRAS signaling (purple), inflammatory response (blue), complement cascade (green),and epithelial mesenchymal transition (turquoise).
Figure also includes genes known to up (orange) or down (olive) regulated in PDAC.Black dashed circle indicates

Table 2 . RNA sequencing statistics by sample. Sample Name Tissue Malassezia Total Reads Total Reads Aligning to hg38 Percent of Total Reads Aligning to Malassezia
https://doi.org/10.1371/journal.pone.0284232.t002