Gemcitabine Recruits M2-Type Tumor-Associated Macrophages into the Stroma of Pancreatic Cancer.

BACKGROUND: Pancreatic ductal adenocarcinoma (PDAC) is a very lethal disease that can develop therapy resistance over time. The dense stroma in PDAC plays a critical role in tumor progression and resistance. How this stroma interacts with the tumor cells and how this is influenced by chemotherapy remain poorly understood.
METHODS: The backbone of this study is the parallel transcriptome analysis of human tumor and mouse stroma in two molecular and clinical representative patient-derived tumor xenografts models. Mice (8 animals per group) were treated for 4 weeks with gemcitabine or control. We studied tumor growth, RNA expression in the stroma, tumor-associated macrophages (TAMs) with immunofluorescence, and cytokines in the serum.
RESULTS: A method for parallel transcriptome analysis was optimized. We found that the tumor (differentiation, gene expression) determines the infiltration of macrophages into the stroma. In aggressive PDAC (epithelial-to-mesenchymal transition high), we find more M2 polarized TAMs and the activation of cytokines and growth factors (TNFα, TGFβ1, and IL6). There are increased stromal glycolysis, reduced fatty acid oxidation, and reduced mitochondrial oxidation (tricarboxylic acid cycle and oxidative phosphorylation). Treatment with gemcitabine results in a shift of innate immune cells, especially additional infiltration of protumoral M2 TAMs (P < .001) and metabolic reprogramming.
CONCLUSIONS: Gemcitabine treatment of PDAC xenografts stimulates a protumoral macrophage phenotype, and this, in combination with a shift of the tumor cells to a mesenchymal phenotype that we reported previously, contributes to tumor progression and therapeutic resistance. Targeting M2-polarized TAMs may benefit PDAC patients at risk to become refractory to current anticancer regimens.

Introduction

Pancreatic ductal adenocarcinoma (PDAC) has a unique tumor microenvironment; it represents the most stroma-rich type of cancer. This stroma can comprise up to 90% of the tumor’s mass [1], and its components are very heterogeneous with multiple types of cells [stellate cells, endothelial cells, cancer-associated fibroblasts, lymphocytes, tumor-associated macrophages (TAMs), dendritic cells, and myeloid-derived suppressor cells (MDSCs)] and acellular components (physical dynamics of blood vessels, extracellular matrix, and soluble proteins) [[2], [3], [4]]. Physically, the dense stroma can block drug delivery by collapsing blood vessels and/or shielding the tumor cells from therapeutics, nutrients, and oxygen [5]. Its cellular components secrete cytokines and growth factors and, together with a hypoxic environment, drive tumor progression characterized by early invasion and metastasis [6,7]. Tumor cells and their neighboring stroma are engage in a two-way molecular communication that determines the tumor’s behavior and therapeutic resistance [1,8,9]. The importance of the stromal cells has recently been recognized but has not been investigated extensively [10,11]. To study these interactions will require biologically complex models such as genetically engineered mouse models or patient-derived xenograft (PDTX) models and techniques to discriminate the origin of the cells and their signals.

In the current study, the interaction between the stroma and tumor cells was investigated to explore the underlying mechanism of tumor progression and survival in two different human PDAC PDTX models with different disease biology. Tumor gene expression profiles of these cancer models were reported in our recent article [12]. The poorly differentiated and metastatic patient we investigated (EMT-high, PAC010) had a pronounced epithelial-to-mesenchymal transition (EMT) signature. In contrast, the moderately differentiated tumor model was from a localized tumor with an EMT-low profile (PAC006) [7,13].

In PDTX, the human stroma is replaced by murine stroma [14]. We used this as an opportunity to dissect the molecular characteristics of the stroma (mouse) and tumor cells (human). To be noted, in PDTX models, we can only investigate the tumor-infiltrating cells from the innate immune system as the host has a nonfunctional thymus whereby a large fraction of the stroma infiltrating cells are TAMs [15]. TAMs are a class of immune cells that are present in high numbers in the microenvironment of solid tumors and originate from bone marrow–derived blood monocytes or yolk sac progenitors. There is growing evidence that links increased recruitment of TAMs into the tumor tissue with that of tumor invasiveness, metastasis, immune escape, matrix remodeling, and therapeutic resistance [[16], [17], [18], [19], [20], [21]]. Conventionally, TAMs are classified as classical activated (M1) and nonclassical activated (M2) that reflect their distinct cellular function and metabolism. M1 macrophages are seen to have a proinflammatory and cytotoxic (antitumoral) function; M2 macrophages are anti-inflammatory (protumoral) and promote wound healing. This classification system might not fully reflect the actual plasticity of TAMs in the tumor microenvironment [[16], [17], [18], [19], [20], [21]].

By comparing the RNA expression of the stroma of EMT-high with that of EMT-low tumor model, we looked for specific pathways, regulatory molecules, and functions induced in the stroma by aggressive tumor cells, in particular related to the innate immune response. Subsequently, we investigated these markers in the stroma following gemcitabine treatment (widely used in the clinic) in our models that represent two distinct molecular and clinical presentations of PDAC. The present study combines basic biological understanding, unbiased molecular analysis of tumor and stroma, and a translational approach. A better understanding of TAMs in PDAC and the networks through which the tumor cells and TAMs communicate to promote tumor aggressiveness may contribute to new therapeutic targets that can augment current treatment.

Materials and Methods

Tumor Implantation, Treatment Procedure, Growth, and Sample Collection

The development of the PDTX models was recently reported [14]. For the current experiments, tumor tissues were implanted subcutaneously, and after the tumor had reached a volume of 100-200 mm3, the mice were randomly divided into groups of 8 mice. Group A (control) was treated with vehicle (0.9% NaCl), and group B (experimental) was treated for up to 28 days with intraperitoneal gemcitabine twice a week [gemcitabine 50 mg/kg (Hospira 38 mg/ml solution, Hospira Benelux, Belgium)]. The animal weight and tumor size were measured thrice a week, and the tumor volume was calculated [12]. Tumor tissue and serum were harvested as soon as the volume reached 1000 - 1500  mm3. At the time of harvesting, tissue samples were weighed, photographed, and stored for histological analysis and molecular profiling.

Animal care and all research procedures were executed in accordance with the applicable legal guidelines and under approval of the medical ethical committee for laboratory animals of the KU Leuven (P147/2012).

Histology and Immunohistochemistry

Immunohistochemistry was done on formalin-fixed, paraffin-embedded tumor tissue section to estimate human (tumor) and mouse (stroma) ratio using an antibody against human-specific Cytokeratin-pan antigen (Abcam, Cambridge, UK). In addition, immunofluorescence double staining was performed on formalin-fixed, paraffin-embedded tissue sections after antigen retrieval. As primary antibodies we used goat anti-mouse MMR/CD206 (R&D systems, AF2535-SP) or rat anti-mouse F4/80 antibody (Cl:A3-1 Bio-Rad, MCA497GA) monoclonal MHC Class II (I-A/I-E) (M5/114.15.2, eBioscience, 46-5321-82). Biotin-SP donkey anti-rat IgG (Jackson ImmunoResearch, 712-065-153) and donkey anti-Goat IgG (H+L) Alexa Fluor 568 (Invitrogen, A-11057) were used as secondary antibodies in combination with DAPI staining. Images were acquired on the Zeiss Axio Scan.Z1 using a ×20 objective and ZEN 2 software; images were processed using the software package QuPath (Version: 0.1.2) [22].

RNA Isolation, Next-Generation Sequencing, and Correct Classification of Reads per Organism in PDTX Samples

RNA was isolated with the RNeasy Kit (Qiagen, Chatsworth, CA) according to the manufacturer's instructions. RNA sequencing and processing were performed by the VIB Nucleomics Core (www.nucleomics.be). In short, the poly-A–containing mRNA molecules were purified, converted into cDNA, and sequenced on an Illumina HiSeq 4000 full flow cell. Performance of Xenome [23] and in silico combined reference genome (ICRG) [24] methods for reads classification was compared using a simulated mix dataset (see supplementary file).

RNA sequencing data are available in Gene Expression Omnibus (NCBI) under number GSE118197.

Gene Expression Analysis

The data were analyzed using two complementary software suites: Ingenuity Pathway Analysis (IPA) (http://www.ingenuity.com) and Gene Set Enrichment Analysis (GSEA) (http://software.broadinstitute.org/gsea/index.jsp). The core analysis of IPA identifies the most significant biological functions and/or diseases and the potential upstream regulators (genes, RNA, and proteins). GSEA identified enrichment of hallmark gene sets as well as KEGG-pathways. Hierarchical clustering (based on KEGG defined pathways, https://www.genome.jp/kegg/pathway.html, or EMT signature [12]) of the mRNA expression of the individual samples was done using PermutMatrix program (Version 1.9.3 EN).

Enzyme-Linked Immunosorbent Assay

Serum levels of mouse interleulin-6 (IL6), interferon-γ (IFNγ), and transforming growth factor-1β (TGFβ1) were determined in duplicate by enzyme-linked immunosorbent assay (ThermoFischer scientific/Invitrogen, Waltham, MA).

Statistical Analysis

All statistics on quantitative reverse-transcription polymerase chain reaction (RT-qPCR) were performed using SPSS v23 (IBM). Statistical differences between groups were assessed with a Student's t test or the Mann-Whitney rank sum test when appropriate. For differences in gene expression assessed by RT-qPCR, ANOVA test with post hoc Tukey's procedure was used. A P value below .05 was considered statistically significant. For survival analysis, we performed a log-rank test using MedCalc software (version 19.0.6, MedCalc Software bvba, Ostend, Belgium; https://www.medcalc.org; 2019).

Results

Design Parallel Transcriptome Analysis of Mouse Stroma Besides Human Tumor in a Single Sample

To test the hypothesis that tumor phenotype is linked to innate immune response, we analyzed the tumor (human) and stroma (mouse) of two PDTX-PDAC models by parallel transcriptome analysis in combination with immunofluorescence and protein analysis. Two tumor types from a previously established panel of PDAC-PDTX were selected (poorly differentiated PAC010 and the well/moderately differentiated tumor model PAC006) [12,14]. After engraftment in a group of immunocompromised mice, the tumor growth was monitored during treatment (Figure 1). To note, difference in growth rate for treated versus untreated PAC006 tumors demanded that we terminated the untreated animals 2 weeks in advance of the scheduled time. We cannot rule out that the moment of termination had an influence on our observations, either compensatory or enhancement; the expression analysis argues against this as we point out in the discussion.

Figure 1

Schematic representation of the study design: generation of patient-derived PDAC xenografts and their molecular characterization. (A) Establishment and validation design of patient-derived PDAC xenograft. (B) Molecular and histochemical characterization of the tumor cells and the stroma in models and in response to treatment.

The backbone of this study is the parallel transcriptome analysis of the human tumor and the mouse stroma. Sequenced reads from PDTX tissue samples are a mix of human and mouse reads that needed to be assigned to correct species prior to expression analysis. We investigated two approaches: 1) the Xenome classification [23] to discriminate human to mouse reads prior to mapping and 2) the In silico Combined human-mouse Reference Genome (ICRG) [24], where the reads are mapped to a combined genome and assigned to their species after mapping. To compare these approaches, we simulated different mix of reads from mouse and human skin fibroblasts sequenced datasets and compared their classification performances (see supplementary data). We found that Xenome discarded 1.2%-1.9% of the reads per sample due to ambiguity, cross-mapped reads, or unclassifiable reads. In our analysis, the ICRG and Xenome approaches on the experimental PDTX data showed similar results. Important to note, the ICRG is quite simple to implement compared to Xenome. We also investigated the effect of the read length (75 bp or 150 bp) and the use of single or paired-end reads, and the results show that the improvement is rather limited compared to increase in sequencing costs. We thus concluded that we can reliably quantify the gene expression for mixed species from a single sample using 75-bp single reads. To be noted, to perform the sequencing with enough coverage, we estimated on forehand the ratio of stroma (mouse) versus tumor (human) tissue using human specific cytokeratin staining (Supplementary Figure 2).

The separation of species-specific RNA reads resulted in two distinct transcriptomic profiles for each sample: one tumor profile based on human RNA sequences from the grafted tumor cells and one stromal profile from the mouse RNA sequences (Table 1). The comparison showed a big difference in gene expression between the tumors in the models (>4000 genes upregulated and >4000 genes downregulated).

Table 1

Number of Differentially Expressed Genes Between Conditions

	Tumor		Stroma
	Uncorr P Value < .001		Uncorr P Value < .001
	log2 ratio < −1	log2 ratio > +1	log2 ratio < −1	log2 ratio > +1
PAC010 vs PAC006 control	4249	4033	587	848
PAC006-GEM vs PAC006	708	1399	686	474
PAC010-GEM vs PAC010	956	926	871	832

EMT Classification of PDAC Models Used in This Study

The growth of the tumors in the PDTX mice, their morphology, and their gene expression indicate two behaviorally different tumor models as were seen in the original patients. Analysis of the expression of a panel of EMT-related genes (defined by Ingenuity consortium, retrieved January 2019, www.ingenuity.com) [12] revealed higher expression of mesenchymal markers in aggressive cancer PAC010 (e.g., ZEB1, STAT3, VIM, and SNAI2) with a reduced expression of epithelial markers (e.g., CDH1 and CTNNB1) (EMT-high signature) compared to the moderately differentiated model PAC006 (EMT-low signature). In a clinical dataset of 118 PDAC patients [25], we performed hierarchical clustering using the reported PDAssigner 62 gene set for pancreatic cancer [26] and using the IPA-EMT gene set; both gene sets resulted in two clusters that showed good agreement (see Supplementary Figure 5). The patients clustered based on EMT signature into EMT-high versus EMT-low and showed significant differences in terms of survival (Figure 2).

Figure 2

Survival analysis of human PDAC classified using IPA-EMT gene set. Hierarchical clustering was performed on 118 PDAC patients data retrieved from NCBI (GSE62165) using 55 EMT-associated genes [12]. H -settings: RNA expression = 2log values; dissimilarity = Pearson’s distance, HC = complete linkage, normalized rows = Z-score, seriation = multifragment heuristics. This signature could separate the patients into two groups (EMT-high patients and EMT-low patients). Using corresponding survival data and log-rank test, disease-free survival (DFS) and overall survival (OS) of these 118 patients were analyzed [25].

Identification of Top Enriched Processes in the Stroma for Aggressive PDAC

We performed a comprehensive gene expression analysis of the stroma components to identify critical regulatory molecules, pathways, cellular functions, and metabolic changes characteristic for an aggressive EMT-high PDAC. We found that between the stroma of EMT-high and EMT-low models (untreated), more than 1400 genes were differentially expressed (Table 1). Using GSEA, we could identify the Hallmark and Pathway gene sets that were most significantly different between both models (Table 2).

Table 2

Top Hallmark Gene Sets Differentially Expressed in the Stroma Between EMT-High (PAC010) and EMT-Low (PAC006) Identified by GSEA. A: Top 10 positively enriched gene sets. B: Top 10 gene sets that were negatively enriched. NES, normalized enrichment score; FDR, false discovery rate. Size is the number of genes of the pathway.

A	Gene Set (Positively Enriched)	Size	NES	FDR q-Val
1	TNFA signaling via NFKB	180	2.14	0.000
2	E2F targets	195	2.14	0.000
3	G2M checkpoint	188	2.09	0.000
4	Inflammatory response	163	1.93	0.000
5	IL6 JAK STAT3 signaling	73	1.92	0.000
6	Interferon alpha response	86	1.91	0.000
7	Unfolded protein response	107	1.90	0.000
8	MTORC1 signaling	197	1.84	0.000
9	Hypoxia	159	1.81	0.001
10	Glycolysis	166	1.79	0.001
B	Gene Set (Negatively Enriched)	Size	NES	FDR q-Val
1	Adipogenesis	174	−2.22	0.000
2	Myogenesis	122	−2.20	0.000
3	Bile acid metabolism	74	−1.98	0.000
4	Oxidative phosphorylation	185	−1.92	0.000
5	Hedgehog signaling	28	−1.84	0.002
6	Estrogen response early	151	−1.61	0.021
7	KRAS signaling DN	80	−1.61	0.018
8	Fatty acid metabolism	128	−1.59	0.017
9	Notch signaling	27	−1.59	0.017
10	Estrogen response late	147	−1.53	0.025

In the aggressive tumor, there was a positive enrichment of the hallmark gene set for hypoxia (normalized enrichment score: NES +1.81) and several cell-cell signaling pathways (TNFα, IL6, etc.). Important differences were also found for gene sets representing metabolism (e.g., glycolysis NES +1.79, oxidative phosphorylation NES −1.92 and fatty acid metabolism NES −1.59).

Using IPA, we also identified the regulator molecules for pathways enriched in the EMT-high tumor model. The top regulators included intercellular signaling molecules (TGFβ1, IL1β, and TNFα) found in both the tumor cells [12] and the stromal compartment (Table 3, A).

Table 3

IPA Analysis of Upstream Regulatory Molecules and Cellular Functions from Differentially Expressed Genes EMT-High Versus EMT-Low in the Stroma. A: Upstream regulator molecules identified by IPA. B: Functional annotation of processes for the cells in the stroma based on their underlying human tumor EMT-high (poorly differentiated, PAC010) versus EMT-low (well differentiated, PAC006) computed from differential gene expression of the mouse stroma.

A	Upstream Regulator	Molecule Type	Corr. P Value	Activation Statea
	IL1β	Cytokine	4.56E-29	Activated
	TNFα	Cytokine	5.22E-29	Activated
	TGFβ1	Growth factor	6.50E-29
	IFNγ	Cytokine	7.78E-27	Activated
	IL6	Cytokine	1.33E-16
	SMARCA4	Transcrip. regulator	3.60E-16
	IL10RA	Receptor	4.03E-16	Inhibited
	NFkβ (complex)	Complex	6.55E-16	Activated
	CEBPA	Transcrip. regulator	2.01E-15
	IL13	Cytokine	2.65E-15

Activation state by IPA; blank = no prediction (see supplementary file).

These signaling molecules can support tumor progression and indicate a possible paracrine communication between the two entities. In addition, importantly, functional analysis of the stromal reaction revealed the increased recruitment of leukocyte (Table 3, B) in EMT-high PDAC, which was subsequently investigated by immunofluorescence staining (see below).

Stromal Reaction in Response to Gemcitabine Treatment in PDTX-PDAC

We investigate the effects of drug treatment on the different models. Gemcitabine treatment resulted in a high number of differentially expressed genes in the tumor and also in the stroma (Table 1). Using GSEA, we identified the hallmark gene sets significantly enriched after gemcitabine treatment (Supplementary Tables 1, A and B, and 2, A and B). We focused on the processes identified for aggressive PDAC (Table 2, Table 3): leukocytes/TAMs, related metabolic pathways, and cell-cell signaling.

A graphic summary of how tumors/human (h) and stroma/mouse (m) responded to gemcitabine is presented in Figure 3. Two groups of related processes (cell-cell signaling and metabolism) are among the highest enriched gene sets. The individual processes were further analyzed through hierarchical clustering, IPA analysis, immunofluorescence, and/or protein determination in the serum.

Figure 3

Graphic presentation of top enriched Hallmark gene sets in tumors and in the corresponding stroma as result of gemcitabine treatment. We analyzed the differential expressed genes for four conditions: (A) hG6NES: tumor PAC006 GEM-treated versus control, (B) hG10NES: tumor PAC010 GEM-treated versus control, (C) mG6NES: stroma PAC006 GEM-treated versus control, and (D) mG10NES: stroma PAC010 GEM-treated versus control. GSEA was used to identify Hallmark gene sets (for full list, see Supplementary Tables 1, A and B, and 2, A and B). The majority of the top gene sets could be functionally grouped into three main groups; they are presented using the normalized enrichment score (NES with P value < .05). The NES is indicated by a color; the intensity is scaled within each row so that the highest enrichment score corresponds to bright red and suppression to bright green. Gray: NES not significant (uncorrected P value > .05).

An IPA analysis was performed for the stroma when treated with gemcitabine (Table 4). This predicts that gemcitabine treatment inhibits TGFβ1, TNF, and IFNγ regulating signals in the stroma (Table 4, A).

Table 4

IPA analysis using information from differentially expressed genes in the stroma of EMT-low (PAC006) and EMT-high (PAC010) in response to GEM treatment. Top upstream regulators (A) and associated cellular function (B) identified with their predicted activation state and P-value.

A Upstream Regulators
	PAC010 vs PAC006		PAC6-GEM vs PAC006		PAC010-GEM vs PAC010
	Act Statea	P-valueb	Act Statea	P-valueb	Act Statea	P-valueb
TGFβ1		10−28	INH	10−34	INH	10−48
TNFα	Act	10−28		10−16	INH	10−28
IL1β	Act	10−28		10−10		10−21
IFNγ	Act	10−26		10−10	INH	10−30
IL-6		10−16		10−17		10−35

Activation state by IPA, INH = inhibited, ACT = activation and blank = no prediction

Corrected P-value. (see supplementary file).

We performed RT-qPCR and enzyme-linked immunosorbent assay to determine relative gene expression of selected cytokines and growth factor: mRNA (IL6, IL10, Hif1α, and TGFβ1) in the stroma and protein level (IL6, TGFβ1, and IFNγ) in circulation. The RT-qPCR confirmed the RNA sequencing. Gemcitabine treatment indicates a suppression of IL6 and TGFβ1 mRNA (Supplementary Table 3). At protein level, of the tested cytokines and growth factors, only IL6 yielded detectable serum levels in the untreated EMT-high model. This level was significantly suppressed following gemcitabine treatment [PAC010 = 234 (75-518) pg/ml vs PAC010-GEM treated = 22.5 (16-31) pg/ml; median (25%-75%); P = .029].

M2-like, CD206 posF4/80 posTAMs Are Increased in EMT-High Tumors

One of the major cell types recruited from the blood are macrophages that infiltrate the tumor, become TAMs, and can acquire in this microenvironment a protumoral or antitumoral phenotype under influence of cytokines, interleukins, and metabolic factors. Macrophages play a critical role in tumor progression. RNA-sequencing analysis revealed an association of infiltration of leukocytes and cytokine signaling (Table 3, Table 4). The effect of gemcitabine treatment on blood cell/leukocyte infiltration was more significant in the EMT-high model compared to EMT-low (Table 4, B). We performed histochemical analyses on formalin-fixed tumor slices using immunofluorescence dual staining for M2 and M1 polarized TAMs. We used mouse specific macrophage marker F4/80, along with CD206high, as indicator of M2 macrophages (Figure 4, C) and the combination of MHCIIhigh and CD206low as indicators of M1 macrophages (Figure 4, A). Gemcitabine treatment caused significantly (P < .001) increased expression of M2-like macrophages in both models (Figure 4, D). Between the models, our result also indicated statistically significantly higher expression of M2 macrophages in the EMT-high model (P = .003), suggesting a link with tumor aggressiveness (Figure 4, D). The expression of M1 markers between the models and for EMT-low in response to gemcitabine treatment was not statistically significant (Figure 4, B). However, in EMT-high, gemcitabine treatment resulted also in an increase of M1-like TAM (P = .004). These results were supported by heat map clustering of gene expression of TAMs markers [27]. The heat map revealed that, in EMT-low, the M1 phenotype was highly variable between the animals. Both M1 and M2 signatures could separate the gemcitabine-treated samples from untreated samples in the EMT-high model, supporting a treatment-induced immunological shift in TAMs (Supplementary Figure 6).

Figure 4

Immunofluorescence staining of PDTX-PDAC tumor sections. Representative images for (A) DPI/MHCII/CD206 and (C) DPI/F4/80/CD206 from placebo- and gemcitabine-treated PDTX mice. Original magnification of histological images  × 10; scale bar 200  μm. The graphs indicate percentage of cells stained for (B) M1-type MHCIIposCD206low macrophages or (D) double-positive M2-type CD206posF4/80pos macrophages. Sections were stained with DAPI, and percentage was calculated to the total number of cells using QuPath software. Each staining is representative for the analysis of tumor sections of four animals per group.

Hierarchical clustering of gene sets representing dendritic cells and MDSCs indicates more MDSCs in the untreated, EMT-high model (Supplementary Figure 6, D). Gemcitabine treatment results in increased antigen presentation (H2 genes and PD-L1 (CD274) are upregulated) (Supplementary Figure 6, C). We see an enrichment of MDSCs following treatment with gemcitabine, accompanied in the EMT-high model by increased expression of the immune checkpoint markers Lag3 (CD223) and Vista (Supplementary Figure 7, E-F).

The Metabolic Profile of EMT-High Tumors Favors M2-Polarized Macrophages

The functional phenotype of M1 and M2 macrophage subtypes is highly regulated at the transcriptional and metabolic level. Proinflammatory M1 macrophages consume glucose and rely heavily on glycolysis for ATP production. They also exhibit a defect in tricarboxylic acid (TCA) cycle. On the other hand, anti-inflammatory M2 macrophages have an intact TCA cycle and favor fatty acid oxidation as mechanism to produce ATP. We used RNA-seq expression data to explore change in glycolysis/gluconeogenesis (KEGG 00010), citrate cycle (TCA) (KEGG 00020), and oxidative phosphorylation (KEGG 00190). When we concentrate on the rate-limiting enzymes for glycolysis, we see by hierarchical clustering that gemcitabine treatment downregulated the majority of these enzymes (Figure 5, A and B). The same downregulation was found when the full gene set was examined (data not shown). This downregulation of glycolysis is supported by the GSEA (Figure 5, C).

Figure 5

Changes in metabolic gene expression in PAC006 and PAC010 induced by gemcitabine. Hierarchical clustering of rate-limiting enzymes of glycolysis in (A) PAC006 and (B) PAC010. (C) Downregulation of glycolysis is confirmed by GSEA (e.g., PAC010). Key enzymes of TCA cycle are shown in (D) for PAC006 and (E) for PAC010. GSEA demonstrates enrichment of KEGG-TCA pathway (PAC006). The ratio of GLUL/GLS as marker for the use of glutamine to fuel the TCA cycle; (G) PAC006 and (H) PAC010. GSEA indicates a strong enrichment of (I) oxidative phosphorylation. The heat map shows the relative mRNA expression in glycolysis and TCA cycle of the mouse genes in the stroma in response to treatment of the xenografts. The red (high), black (middle), and green (low) colors indicate the relative expression intensity of each gene within a sample.

Discussion

A diagnosis of pancreatic cancer comes almost always with limited therapeutic perspectives and a bad prognosis for the patients. Several factors are considered important in this respect, one being the unique, stroma-rich microenvironment of a PDAC. The communication between the tumor and the stromal cells determines the evolution of the disease, but the underlying mechanism is only partly understood. To study the tumor-stroma interactions and to gain insight into possible therapeutic targets to improve treatment, we used two human tumor xenograft models of different clinical and morphological presentation. Gene expression analysis confirms the aggressive/poorly differentiated EMT-high model (PAC010) and the more moderate/differentiated tumor EMT-low model (PAC006). By comparing the gene expression of the stroma between these models, we looked for specific pathways, regulatory molecules, and functions characteristic for aggressive PDAC. Subsequently, we used these observations when we investigated the effect of gemcitabine treatment on the stroma.

In the untreated PDAC models, gene expression and functional analysis of the stroma indicated increased infiltration of leukocyte and macrophages. TAMs represent the most abundant leukocyte subpopulation in the stroma; they belong to the innate immune system [15] and can release signaling molecules facilitating cancer cell invasion, migration, angiogenesis, tumor progression, or metastasis [[28], [29], [30], [31], [32]]. In particular, TAMs secrete cytokines and growth factors that are involved in the induction of EMT in solid tumors [3,[33], [34], [35], [36]]. In PDTX, we could detect increased levels of IL6 in the serum of animals transplanted with EMT-high PDAC model in agreement with the gene expression analysis. Clinically, high serum levels of IL6 correlate with lower rates of survival and a higher chemoresistance [37]. The similarity of pathways found in the tumor [12] and in this study in the stroma suggests a possible paracrine communication between the two entities.

Within the tumor microenvironment, the TAMs come into contact with tumor-secreted factors that can polarize them toward M2 type macrophage [28,29,[38], [39], [40], [41], [42], [43]]. Immunofluorescence dual staining confirms a significantly higher presence of M2 polarized TAMs (CD206high on F4/80 mouse specific macrophage marker) in the EMT-high model. M2 macrophage infiltration is reported to be associated with a malignant tumor phenotype that correlates with poor prognosis [44,45].

Does the tumor-stromal interaction change when the tumor is treated with chemotherapy, and is such a change related to the phenotype of the PDAC? Gemcitabine-based combination chemotherapy is the standard therapy for most PDAC patients [46]. Resistance to gemcitabine and other drugs remains a major cause of therapeutic failure. Recently, we showed that gemcitabine treatment significantly induced the expression of mesenchymal markers in the residual tumor cells, indicating enrichment of metastatic and therapeutic resistance phenotype [12]. Other studies have also associated drug resistance of tumor cells with the acquisition of EMT phenotype [47,48]. To what extent the stroma under chemotherapy influences the survival is unknown. Gemcitabine treatment does not shift the microenvironment to a stroma as found for an aggressive PDAC (see GSEA or IPA) within the time frame we studied the PDTX. In contrast, cytokines like TNFα that are activated in EMT-high PDAC are inhibited by gemcitabine. However, in the present study, we clearly demonstrate that gemcitabine treatment induces an increased infiltration of M2 polarized TAMs in stroma of PDAC xenografts. In the more aggressive/EMT-high PAC010, following gemcitabine treatment, both TAM subtypes (M1 and M2) have increased. Interestingly, TAMs can contribute to EMT by exhibiting both pro- and anti-inflammation characteristics [49,50]. We further explored the metabolic pathways predominantly used in the stroma following gemcitabine treatment. The results support an increased M2 population. We found reduced glycolysis and upregulation of mitochondrial oxidative metabolism (TCA cycle and oxidative phosphorylation). Previous reports define them as metabolically M2 polarized TAMs [51,52]. In addition, upregulation of fatty acid oxidation [53] and increased GLUL/GLS ratio (indicating enhanced fueling of the TCA cycle through glutamine [27]) were identified. Unfortunately, in our study, we did not find a clear association between M2 macrophage infiltration and cytokines signaling (upstream regulators) as previously was reported (M1: commonly expresses higher levels of IL12, IL23, TNFα, and IL6 and M2: commonly expresses higher levels of IL10 and TGFβ1) [54]. The plasticity of the TAMs in the stroma allows also for cells with shared features of both M1 and M2, and the phenotype of TAMs can switch during different stages of tumor progression; these observations could explain this discrepancy [[55], [56], [57]]. There is some controversy in literature on how MDSCs react to gemcitabine treatment. Galluzi et al. report a reduction following treatment [58], while in the study by Plate et al. in patients with pancreatic cancer [59], after an initial reduction of the major MDSC subset (BDCA/CD1b), they found a rebound to levels equal or above baseline. In the present study, gene expression indicates that the MDSCs have increased and exhibit a stronger immune suppressive phenotype following gemcitabine treatment.

In a recent paper, it was reported that increased M2 TAMs in the PDAC stroma are associated with decreased survival under gemcitabine treatment. As a mechanism, they proposed that macrophage-released pyrimidines inhibit efficacy of gemcitabine therapy in pancreatic cancer [60]. This might be another mechanism by which M2 TAMs and tumor-stroma interactions have an impact on the patient’s prognosis. These observations as well as ours make M2 an attractive therapeutic target, especially in metastatic and therapy resistant PDAC [[61], [62], [63]]. Currently, a number of preclinical and clinical trials have been completed or are ongoing, targeting TAMs to treat different tumor types including pancreatic cancer (e.g., NCT03662412, NCT03184870, NCT01921699). These strategies aim to block monocyte recruitment [64], switch TAMs phenotype from M2 to M1 [65,66], deplete resident TAMs [67,68], or neutralize TAMs products [62,69].

Conclusion

Our PDTX models, representing two different clinically relevant phenotypes of PDAC, will contribute to a better understanding of tumor-stroma biology. When we take the expression analysis together of both PAC006 and PAC010, we think the effects in the stroma are firstly the result of the tumor EMT phenotype and the tumor-stroma interactions. Gemcitabine treatment results in a shift in the stroma to a predominantly immunosuppressive environment with M2 TAMs and MDSCs. This suggests a dynamic in the stroma that supports a mesenchymal phenotype of the tumor cells [12], influx of macrophages, and suppression of inflammation in the context of gemcitabine treatment. This molecular analysis identifies relevant pathways and molecular targets for new therapy development. Targeting M2-polarized TAMs may benefit PDAC patients at risk to become refractory to current anticancer regimens.

Conflict of Interest

All authors declare no potential conflicts of interest.

Funding Support

C. V. holds a mandate as Senior Clinical Investigator of the Research Foundation-Flanders (Belgium) (FWO). This study was partly supported by a research grant from “Kom op tegen Kanker” Belgium and VUYLSTEKE-FLIPTS FONDS LEVERKANKER.