Pancreatic ductal adenocarcinoma (PDAC) is among the deadliest malignancies and potentially curable only with radical surgical resection at early stages. The tumor microenvironment has been shown to be central to the development and progression of PDAC. A better understanding of how early human PDAC metabolically communicates with its environment and differs from healthy pancreas could help improve PDAC diagnosis and treatment. Here we performed deep proteomic analyses from diagnostic specimens of operable, treatment-naïve PDAC patients (n = 14), isolating four tissue compartments by laser-capture microdissection: PDAC lesions, tumor-adjacent but morphologically benign exocrine glands, and connective tissues neighboring each of these compartments. Protein and pathway levels were compared between compartments and with control pancreatic proteomes. Selected targets were studied immunohistochemically in the 14 patients and in additional tumor microarrays, and lipid deposition was assessed by nonlinear label-free imaging (n = 16). Widespread downregulation of pancreatic secretory functions was observed, which was paralleled by high cholesterol biosynthetic activity without prominent lipid storage in the neoplastic cells. Stromal compartments harbored ample blood apolipoproteins, indicating abundant microvasculature at the time of tumor removal. The features best differentiating the tumor-adjacent exocrine tissue from healthy control pancreas were defined by upregulation of proteins related to lipid transport. Importantly, histologically benign exocrine regions harbored the most significant prognostic pathways, with proteins involved in lipid transport and metabolism, such as neutral cholesteryl ester hydrolase 1, associating with shorter survival. In conclusion, this study reveals prognostic molecular changes in the exocrine tissue neighboring pancreatic cancer and identifies enhanced lipid transport and metabolism as its defining features.

Significance:

In clinically operable pancreatic cancer, regions distant from malignant cells already display proteomic changes related to lipid transport and metabolism that affect prognosis and may be pharmacologically targeted.

Pancreatic ductal adenocarcinoma (PDAC, here pancreatic cancer) is projected to become the second leading cause of cancer deaths by 2030 (1, 2). Radical surgery combined with chemotherapy has modestly improved the prognosis, but recurrent tumors and distant metastases are common, and the 5-year survival rate is only about 10% (3). The identified driver mutations in KRAS, p53, SMAD, and CDKN2A are present in other cancers with better prognoses (4, 5), signifying the need to identify additional factors that contribute to this lethal malignancy.

Among the modifiable risk factors associated with the development of PDAC, obesity with dyslipidemia and associated cardiometabolic complications are globally on a rapid rise. Obesity and fatty infiltration of the pancreas have been recognized as risk factors for pancreatic precancerous lesions (6), but how lipid metabolism contributes to the development of human PDAC is not well understood. Aberrant accumulation of cholesteryl esters, i.e., storage forms of cholesterol generated by the enzyme sterol O-acyltransferase (SOAT1, a.k.a. ACAT1) and deposited in intracellular lipid droplets, was reported in human pancreatic cancer specimen (7). Moreover, SOAT1 inhibition suppressed pancreatic cancer cell growth and was proposed as a potential treatment strategy in PDAC (7, 8). Additional studies in mice suggested that pancreatic cancers exhibit increased low-density lipoprotein (LDL) receptor–dependent cholesterol uptake (9) and that stromal cells active in cholesterol synthesis supply cholesterol to tumors that in turn express cholesterol import proteins (10).

The fibrotic microenvironment, known as the desmoplastic reaction, represents one of the distinctive features of PDAC from early on (1, 11). It is related to cancer development in multiple ways, for instance, by promoting metabolic rewiring of tumor cells and generating a barrier to perfusion and treatment, but also by limiting tumor progression and dissemination (12). Understanding the coevolution of the tumor–stroma interplay may represent one of the keys to early detection and treatment of PDAC (13). Yet, systematic information on the compartmentalization of lipid metabolic activities between the neoplastic cells and the surrounding stroma in human pancreatic cancer is currently lacking.

The present study was undertaken to investigate the distribution of metabolic activity across the exocrine and stromal areas in human PDAC, with a focus on lipid metabolism. We probed the expression of candidate cholesterol uptake and synthesis proteins from 24 tumor microarrays and analyzed lipid deposition from 16 fresh cryosections. We also conducted an unbiased, laser-capture microdissection-based deep proteomic exploration from 14 diagnostic PDAC specimens. In addition to the neoplastic cell clusters and their intervening stromal regions, we isolated morphologically benign exocrine glands and their neighboring connective tissue from the tumor margin, i.e., the area that histologically resembles healthy exocrine tissue and surrounds the cancerous lesion. We considered these samples not only to serve as intraspecimen controls for the neoplastic areas but possibly also to reveal potential incipient changes related to the malignant transformation process. With this in mind, we compared our findings to reported 12 healthy pancreatic proteomes. Our results reveal widespread lipid metabolic alterations in the neoplastic areas, abundant vasculature-related pathways in the stromal regions and more subtle, but prognostically relevant changes in the clinical tumor margin.

Patient materials

All patients were operated on in the Helsinki University Hospital for pancreatic cancer between 2000 and 2020, and the tissue materials were collected from the surgeries either as snap-frozen fresh biopsies or as diagnostic FFPEs. Of the 14 samples used for proteomic analysis, 12 originated from pancreaticoduodenectomies and 2 from distal pancreatectomies with random anatomic orientation of the diagnostic specimen. Healthy control pancreatic samples originated from the Helsinki Biobank. The study was conducted in accordance with the Declaration of Helsinki, patients provided written informed consent, and the institutional ethical committee approved the research with licenses 226/E6/2006 (continued 2013) and HUS/2075/2016.

Proteomic sample preparation

After excluding outliers on bilirubin, amylase, age, and survival time, the patients’ diagnostic slides were prescreened for suitable histology and sufficient tissue material for harvesting. Original FFPE blocks for the patients were cut to 10 μm or 20 μm thickness (batches 1 and 2) using a standard floating microtome and sections collected on MembraneSlide NF 1.0 PEN -slides. UV-pretreating the MembraneSlides (14) followed by prolonged drying (16 hours at 56°C) ensured that the samples adhered properly before deparaffination rehydration in xylene-alcohol series. The slides were then stained manually in excess Lillie-Mayer hematoxylin for 2 minutes without differentiation, followed by eosin (0.1% in EtOH). The stained samples were dehydrated and mounted on coverslips with glycerol, as this allowed gentle removal of the coverslips after digitalizing the slides. Immersing the slides in H₂O for up to 16 hours removed the coverslips. Finally, the slides were air-dried and stored at +4°C for up to 2 weeks before laser microdissection. A pathologist carefully outlined the segments to be dissected from the samples by drawing ROIs on the digitalized slides using ZeissBlue software (RRID:SCR_013672). The created data set including the segmented tissues of interest that were used for the proteomic analysis is openly available. Note that as the segmented regions originate from a single section, they do not fully capture the 3D distribution of tissue types in the tumor environment. An in-house software, CL2M, now renamed and published as BIAS, performed the alignment and conversion of the ROIs to LMD platform–compatible versions (15).

Brightfield imaging

Brightfield images of the slides were collected using a Zeiss AxioImager 2 microscope equipped with a color camera (Zeiss AxioCam 105), a motorized stage, and a 20× air-immersion objective (Plan Apochromat, 20×/0.8). For image acquisition and stitching, we defined a tile-scanning experiment with manual focus mapping for each slide individually using Zeiss Zen 2 Blue software. The digitalized slides have roughly 0.22 μm/px lateral resolution and negligible stitching artifacts.

Laser-capture microdissection

To precisely excise the ROIs selected by a pathologist, we registered the data between the scanner (Zeiss AxioImager 2) and the laser-microdissection microscope (Zeiss PALM), according to (15, 16). Three landmark positions were defined and used as position and orientation references to transform data from one microscope to another. Using these reference points, coordinates defining the cutting region for an ROI can be transformed from the source image coordinates to the target microscope coordinates. After generating ROIs, the selected regions were isolated using the Zeiss PALM system with a 63× LCM-compatible magnification objective (LD Plan-Neofluar, 63×). The cutting was carried out using the ultraviolet N2 laser microbeam system of Zeiss PALM, emitting 3-ns pulses. Finally, the isolated regions were catapulted into a cap of PCR tubes containing 5 μL PBS.

Mass spectrometry

Deparaffinization of FFPE tissue, antigen retrieval, and protein extraction

FFPE-derived protein extracts were digested using the S-Trap method following the manufacturer's instructions with a few modifications (17). FFPE samples were incubated with 10 μL of 1:50 diluted EnVision FLEX Target Retrieval Solution High pH (Agilent Dako) at 97°C for 10 minutes with shaking (500 RPM). Then, the samples were subjected to a brief centrifugation at 14,000 × g at 4°C for 3 minutes, followed by the removal of the EnVision solution and the paraffin. These steps were repeated up to three times, ensuring that the samples were paraffin-free.

The protein extraction was performed by adding 20 μL of extraction buffer per section, 100 mmol/L TEAB containing 25 mmol/L DTT, and 10 w/v% SDS (pH 8). The samples were then incubated at 99°C for 1 hour with shaking (500 RPM) to induce antigen retrieval and to break the amino acid cross-links. The samples were then centrifuged to remove the Tris–HCl buffer. Then, 50 μL of extraction buffer 100 mmol/L AmBic was added per sample. The samples were then sonicated to reduce viscosity and shear the DNA in the Bioruptor Plus UCD-300 (Diagenode) for 20 cycles (15 seconds on and 15 seconds off) at 4°C followed by centrifugation at 20,000 × g at 18°C for 10 minutes for clarification of the lysates.

After extraction and reduction, protein determination was performed using the Pierce 660 nm Protein Assay (Thermo Scientific) according to the manufacturer's instructions. The protein determination for the samples containing SDS included the addition of an ionic detergent compatibility reagent (Thermo Scientific). For the different experiments, 50, 100, or 150 μg of proteins were digested, and all digestions were performed in triplicate. One neoplastic stromal sample was lost during the mass spectrometry analysis.

S-Trap digestion

The lysates were alkylated by adding 50 mmol/L iodoacetamide followed by a 30-minute incubation at room temperature in the dark. This step was followed by the addition of 1.2% phosphoric acid (final concentration) to acidify the samples. Then, S-Trap binding buffer (90% methanol, 100 mmol/L TEAB) was added to each sample to a final volume of eight times the volume of the lysate. Digestion was performed on the S-Trap Micro Spin Columns (ProtiFi).

Briefly, for the tryptic digestion, the samples were loaded onto the S-Trap filter followed by a short centrifugation. The captured proteins were washed with S-Trap binding buffer and centrifuged between each wash. After that, the S-Trap column was moved into a new receiver tube/plate. Digestion buffer (50 mmol/L TEAB) containing trypsin in a 1:50 w/w ratio (enzyme/protein) was added on top, and the samples were incubated overnight at 37°C. The following day, the peptides were eluted first with the digestion buffer. The second elution was performed with 0.2% aqueous formic acid, and the final elution consisted of 50% acetonitrile containing 0.2% formic acid. The eluted peptides were then acidified with 50% formic acid to a final pH of ∼3 and dried down in a centrifugal evaporator.

nLC-MS/MS analysis and database search

The nLC-MS/MS analysis was performed with an Ultimate 3000 RSLC nano pump (Thermo Scientific) coupled to a Q-Exactive HF-X (Thermo Scientific) mass spectrometer equipped with an EASY-Spray ion source. The peptides were reconstituted in 2% ACN with 0.1% TFA, and their concentration was measured with Pierce Quantitative Colorimetric Peptide Assay according to the manufacturer's instructions (Thermo Scientific).

Approximately 1 μg of peptides was injected and loaded onto an Acclaim PepMap 100 C18 (75 μm × 2 cm, 3 μm, 100 Å, nanoViper) and then separated on an Acclaim PepMap RSLC C18 column (75 μm × 25 cm, 2 μm, 100 Å; Thermo Scientific). A flow rate of 300 nL/minute and a column temperature of 45°C were used. The solvents used for the nonlinear gradient were A (0.1% formic acid) and B (0.08% formic acid in 80% ACN).

For the detail of LC gradients and data acquisition methods, a 145-minute method was used, with a gradient developing from 2% to 25% of solvent B in 115 minutes, then increasing to 32% in the next 10 minutes, and to 45% in 7 minutes, and then further increasing to 90% in 8 minutes, staying at 90% for the last 5 minutes. The samples were analyzed with a top 20 DDA (data-dependent acquisition) method. Full MS scans were acquired in the Orbitrap mass analyzer, with a range of 375–1,500 m/z, a resolution of 120,000 (at 200 m/z), a target AGC value of 3e6, and a maximum injection time of 100 ms. The 20 most intense peaks (charge state ≥ 2) were fragmented in the HCD cell (NCE of 28), and MS2 scans were acquired with a resolution of 15,000 (at 200 m/z), a target AGC value of 1e5, and a maximum injection time of 50 ms. The precursors were isolated with a 1.2 m/z window, the ion selection threshold was set to 8e3, and the dynamic exclusion was 40 seconds.

All samples were analyzed in random order. HeLa protein digest standard (200 fmol, Pierce) was run regularly as a QC standard to check the system performance. All raw files were analyzed with Proteome Discoverer 2.3 (Thermo Scientific, RRID:SCR_014477) using SEQUEST HT against the UniProtKB human database (downloaded 2020-03-06, canonical and isoforms RRID:SCR_004426; ref. 18). The mass spectrometry and proteomic data that support the findings of this study are openly available at ProteomeXChange with reference number PXD026371. Reprocessed data from Kushner and colleagues are available in Supplementary File S1.

Proteomic data processing

Normalization within data sets

All proteomic data analysis was performed in R (v4.0.3, RRID:SCR_001905), mostly using the tidyverse package collection (v1.3.1, RRID:SCR_019186). Raw protein intensities were first log2 transformed and then median centered sample-wise to generate normally distributed log-median normalized protein abundances. Batches were identified manually based on the date of the mass spectrometry, and the degree of the batch effect was verified using principal component analysis (PCA).

Batch effects were corrected using ComBat (RRID:SCR_010974) by using tissue type as a covariate, in addition to the batch itself, to retain tissue type-specific features (19, 20). Abundances were shifted to the positive direction with 1 prior projection to log-space to retain 0-intensity observations. Proteins were required to be present in >50% of samples in all batches and in >50% of samples of at least one tissue type for batch-effect correction. This approach limits the coverage of proteins studied, ensuring that the pathway-level comparisons are sufficiently independent of the number of identified proteins while still permitting tissue-type enrichment of proteins. These log-median normalized and batch-effect corrected protein abundances are referred to as “normalized values” throughout the article.

Pathway expression quantification, circle plots

Reactome pathways and their constituent protein identifiers were used as the basis for pathway-level analysis (21). For each pathway, its expression value was quantified by averaging the normalized abundances of proteins comprising the pathway. Pathways with <3 detected proteins, or with less detected proteins than the ceiling of the square root, or with <15% of the number of constituting proteins were omitted from the later analysis. The square-root filter is effective only for pathways containing 4 to 53 proteins and limits the uncertainty of pathway expression values related to pathways with only a few constituent proteins, e.g., for pathways with less than 7 constituent proteins, a coverage of 15% is achieved with only 1 detected protein, which we do not consider meaningful. The Reactome database (RRID:SCR_003485) numbers 15,580 distinct UniProt ids against 3,127 in our normalized data, meaning that we have covered 21.2% of the Reactome proteome, and conversely, that the 15% filter necessitated that the number of detected proteins for a pathway was at least 70% of the expected number.

Metabolic circle plots are based on Reactome hierarchical structure, individual pathways represented by nodes and generated using ggraph package. The mean difference in pathway expression values between tissue types is visualized with a red–blue color axis.

Volcano plots

The volcano plots are based on normalized protein abundances or pathway expressions. The mean difference between tissue types is calculated and plotted on the x-axis, whereas the two-tailed false discovery rate corrected, negated, and log-transformed P value is plotted on the y-axis. On the right-hand y-axis, we also display the density of the observations in the background.

Sparse PCA analysis and associated heatmaps

The input data, i.e., the pathway expression values or normalized protein abundances, were subject to robust sparse PCA (22). The analysis was run with default parameters except for the dimensions used and the alpha-sparsity controller, which was set to 2.25 × 10−4. For pathway-level analysis, the first two PCA dimensions were used, whereas the first three were used for protein-level analyses. This helped to resolve interpretable tissue type–specific features and vastly improved computation time.

The heatmaps were generated using the ComplexHeatmap-package (RRID:SCR_017270). The observations, i.e., proteins or pathways, with nonzero weight for the analyzed dimensions were a subset from the original input data. The subset data were then input to the heatmap function (RRID:SCR_017270). Running 10,000 repetitions ensured a robust fit for columns clustered using k-means with 3 or 4 centers and Euclidean distance. The rows were hierarchically clustered using the average method with default parameters.

Survival analysis on proteomic data

Patients were split into long (n = 7) and short (n = 7) survival groups by a two-year survival threshold. The difference in normalized protein abundances between these groups was then assessed with absolute difference in abundances (Over2 – Under2), and their significance was quantified with heteroscedastic two-tailed Student t test without correction for multiple testing. Most prognostic features were defined as having the highest Euclidean distance from origo in volcano plot with scaled axes, i.e., the proteins with the highest combined effect and significance were selected. Housekeeping-normalized values were used for plotting when the control pancreas was present in the plot. Survival data for this study are available in Supplementary File S2.

Normalization of proteomic data across data sets

Protein abundances between different studies or data sets were normalized using what we call housekeeping normalization. First, we defined a set of proteins that have a consistent expression across samples and then normalized all values within an experiment to the median and standard deviation of this reference set. The housekeeping proteins were defined as the proteins with the least coefficient of variation between two multitissue systematic proteomic screening studies (23, 24). Up to 500 housekeeping proteins were used for normalization in both cases.

Additionally, we utilized an alternative approach that does not require data to be batch-effect corrected or complete. Group factor analysis (GFA) is a Bayesian matrix factorization method that aims to disentangle the observed data into components (factors) that are present only in a specifically defined subset of samples (groups; ref. 25). Formally, GFA aims at reconstructing the observed data matrix Y from the product of two matrices Y(m)  =  XW(m)T, where X is called the latent variable matrix and W(m) is called the projection or weight matrix for sample group m. By defining the sample groups to be the 5 tissue types (adjacent and neoplastic parenchyma, control parenchyma, adjacent, and neoplastic stroma), we exploited GFA to extract factors from the protein abundance matrix Y (proteins × samples) that show similar characteristics (weights) for samples within a tissue type but are different across them. The hierarchical model of GFA includes a unique variance parameter for the distribution of the factor weights within a group, hence the model prefers decompositions that are active across whole tissue types.

We used the GFA package in R with K = 50 as the initial number of factors and all other parameters set to default (26). The package handles missing values intrinsically, and hence no data filtering was needed for the analysis. During training, excessive components were eliminated, retaining K = 28 components from the initial 50. The components were ordered according to decreasing effect size measured as | $\overline {| {{X}_f} |} \times \overline {| {{W}_f} |} $ |⁠, where Xf and Wf denote the columns for factor f in the latent variable and projection matrix, respectively, and | $\overline {\ | . |\ } $ | denotes averaging over absolute values. To get a better understanding of the identified factors, we extracted the proteins and pathways that characterized the latent variable space of the factors. We used the genome-wide pathway enrichment analysis tool of String v11.0b (available at string-db.org; ref. 27) on the latent variable values with default settings (FDR threshold 1%) to determine the pathways related most closely to the factor. The Kaplan–Meier survival analysis was conducted on the weights of the factor using the median to divide patients into two groups.

Proteomic data processing of verification studies

Processed proteomic abundances were used directly from the supplementary data of Le Large and colleagues (28) and Cao and colleagues (29), and subjected to the same processing pipeline that was used for our data. In the case of the Le Large and colleagues study, only patients with available paired neoplastic parenchymal and stromal samples were kept, resulting in 14 patients analyzed. In the case of the Cao and colleagues study (29), 37 patients were analyzed upon filtering with the following criteria: sufficient tumor cellularity (105 samples in the original study), available paired adjacent sample, and at least a 30-day long follow-up time. For estimating survival significance, we split the patients into long (n = 16) and short (n = 15) survival groups with a threshold of 400-day postoperative survival. Individuals censored from the study earlier than this threshold were excluded from the analysis (n = 6).

IHC of NCEH1 and AP2M1

Four-micrometer slices were cut with microtomy from FFPE diagnostic blocks (the same blocks as used for the proteomic analysis extended by FFPE blocks for the control samples) and dried for 1 to 2 hours at most 60°C. We proceeded directly with antigen retrieval by incubating the slides for 15 minutes at 98°C in high pH (pH = 9) EnVision Flex target retrieval solution (DM828) using an Agilent Dako Pretreatment Module. Endogenous peroxidase was blocked with EnVision Flex peroxidase-blocking reagent (SM801) for 15 minutes. Primary antibodies (NCEH1: Sigma-Aldrich, HPA026888, RRID:AB_10601557 and AP2M1: Proteintech, 27355-1-AP, RRID:AB_2880853) were incubated overnight at 5°C in Dako REAL Antibody Diluent (S2022) at 1:200 dilution. Samples were incubated for 20 minutes in EnVision Flex/HRP (SM802) secondary antibody. Chromogen development was performed by incubating the samples for 10 minutes in EnVision Flex DAB (DM827), followed by 1-minute incubation in Dako Mayer's hematoxylin (S3309) to stain nuclei. The entire staining process was done with Autostainer 480S, LabVision Corp. Thereafter, the specimens were dehydrated in aqua-70%-96%-AA-ZYL and mounted in Pertex Histolab media. The slides were screened with Pannoramic 250 Flash III (3DHISTECH Ltd.) using a single focal plane at 20× resolution. QuPath (v0.3.0, RRID:SCR_018257; ref. 30) was used to manually define 3 classes of regions of interest (neoplastic, adjacent, and control) from the slides (2–34 ROI/class/slide, median 6) and exported with automatic scripts for further processing. Color deconvolution was performed with manually identified parameters in Matlab (R2020a v9.8 and RRID:SCR_001622). Figures were created in R (v4.0.3, RRID:SCR_001905) and ggplot2 (v3.3.5, RRID: RRID:SCR_014601). All codes including the groovy scripts for exporting from QuPath and a user-friendly GUI to define color-deconvolution parameters in Matlab are publicly available at https://bitbucket.org/szkabel/pirhonen-et-al-2021.

Sample preparation and imaging by coherent anti-Stokes Raman scattering

Special attention was paid to the preservation of lipids, as they are easily disturbed by sample processing, such as elution with solvents or mechanical stress. First, fresh tissue biopsies were snap-frozen immediately and kept at −80°C before utilization. Samples were equilibrated to −20°C before cryosectioning to 20 μmol/L thickness with a Leica CM3050 S cryostat. Sections were gently gathered on SuperFrost (Thermo Scientific J3800AMNZ) slides and air-dried for 30 minutes at room temperature. The prepared slides were stored at −20°C before utilization. Upon imaging, slides were defrosted and immediately fixed using 4% paraformaldehyde solution at 4°C for 30 minutes letting the sample equilibrate to room temperature. Slides were then rinsed with milli-Q water, autofluorescence was quenched using 50 mmol/L NH₄Cl, and slides were rinsed again. Two drops of pure glycerol were used for mounting the coverslips, sealed in place using nail polish.

Samples were imaged within a couple of hours after thawing using a Leica SP8 coherent anti-Stokes Raman scattering (CARS) confocal microscope controlled with Leica LAS X Life Science software. The instrument is equipped with an ultra-short pulsed light source (picoEmerald, APE), producing synchronous 5–7 ps pulses required for CARS microscopy. The Stokes beam was set at 1,064 nm, and the pump laser was tuned to 816.5 nm using an optical parametric oscillator to correspond to CH₂ bonds' vibrational frequency at 2,845/cm. Pulses were aligned spatially and temporally to the confocal plane for the selected objective (Leica 25×/0.95 HCX IR APO L, water immersion). The forward propagated signal was then diverted and passed through bandpass filters to external detectors (Leica HyD). 10 × 10 tile scans were captured using focus mapping with 6 manually selected focus points.

Nail polish, still malleable, was removed by incubating the slides in pure EtOH overnight. Samples were then developed using 2 minutes Lillie-Meyer progressive hematoxylin staining, followed by 0.1% Eosin-Y in EtOH with 0.01% acetic acid. The samples were finally rinsed, mounted on coverslips with Pertex (Histolab) and imaged as described in the brightfield imaging section.

The captured images were aligned with Zen Zeiss 2 Blue, and the signal was quantified from the regions of interest, i.e., neoplastic parenchyma, neoplastic stroma, and the adjacent fatty tissue. Background CARS intensity, as detected from outside the sample area, was subtracted from the signal before reporting the measured values.

IHC of tissue microarrays

Tissue microarray (TMA) blocks were assembled with 1-mm tumor punctures collected from representative areas based on reevaluated hematoxylin and eosin–slides. Four punctures were collected per patient and inserted into the donor block with a semiautomatic microarrayer (MTABooster Version 1.01 for Beecher Manual Arrayer; Alphelys Beecher Instruments).

Subsequently, 4-μm thick parallel sections were cut from these FFPE TMA slides, using a standard floating microtome. Quickly dried (1 hour at 56°C) slides underwent deparaffination rehydration in a xylene-alcohol series. We used Triton 0.1% for 15 minutes to permeabilize the cell membranes. Antigens were retrieved by incubating the samples for 20 minutes at 90°C in citrate buffer (10 mmol/L citric acid, 0.05% Tween 20, pH 6.0). Endogenous peroxidase was blocked with 0.3% H₂O₂ in PBS for 15 minutes, and nonspecific antibody binding with 0.2% BSA (10% FBS in TBS for DHCR7 antibody) for 30 minutes.

Primary antibodies were incubated at +4°C for 16 hours in their respective blocking solutions with the following concentrations: LDL-receptor (LDLR) antibody 1:100 (Santa Cruz sc-18823, RRID:AB_627881), smooth muscle actin (SMA) antibody 1:100 (Abcam ab5694, RRID:AB_2223021), and DHCR7 1:50 (in-house, details below). For LDLR and SMA, MaxPO-multi (Histofine, 414151F) was used as the secondary antibody with a 30-minute incubation. For DHCR7 antibody, samples were incubated at +37°C for 1.5 hours with 1:250 goat anti-rabbit IgG HRP (Bio-Rad, 1706516, RRID:AB_11125547). AEC (Histotec, 415184F) was used for chromogen development: 30 minutes for αLDLR and αSMA, for αDHRC7 only 1 minute was sufficient. Finally, after a 2-minute Lillie-Meyer progressive hematoxylin staining, the samples were rinsed, and coverslips mounted using 45 μL Mowiol-DABCO.

DHCR7 antibody generation and validation

The coding region of the N-terminal region of human DHCR7 (amino acids 1–36 in tandem, of NM_001360) was synthesized using oligonucleotides and cloned into the pTAT-5 expression vector (Addgene, plasmid: #112593). The His-TRX–DHCR7 fusion protein was expressed in BL21(DE3) cells and purified by Ni-NTA Sepharose (according to the manufacturer, GE Healthcare). The protein was used for the immunization of three rabbits (Pineda, Antikörper-Service). The antisera were affinity-purified by using His-TRX-DHCR7 linked to Pierce NHS-Activated Agarose Slurry, according to the manufacturer (Thermo Scientific).

For epitope blocking, the DHCR7 antibody was diluted 1:10 in 0.2% BSA and incubated for 2 hours at room temperature with or without the antigenic peptide (2 mg/mL). The solution was centrifuged 13,000 RPM for 5 minutes and the supernatant was diluted 1:25 prior to the usual staining regime.

TMA data processing

The individual microscope images for the CK7 staining were corrected against uneven illumination using CIDRE (31). The corrected tiles were stitched using Fiji (RRID:SCR_002285) and Grid/Stitching Plugin (32, 33) to create a single image per spot. The registration of the corresponding tissue sections started with preprocessing. The sample areas were detected using Ilastik v1.3.2 (RRID:SCR_015246; ref. 34) by training models for each stain separately on randomly selected five images/stain. The resulting grayscale probability images for sample areas were used to determine pairwise initial registrations by utilizing an affine framework relying on the Alpha AMD similarity measure (35). The images were downscaled to 10% of their original size for these preprocessing steps.

Subsequently, the original images were transformed with the initial registration and registered by using the same Alpha AMD framework for a second time. RGB images were internally converted to grayscale using the scikit-image module v0.17.2 (RRID:SCR_021142) of Python v3.6.6 (RRID:SCR_008394). In both registration rounds, the performance of the automatic alignment is assessed by computing the 2D correlation coefficient between the saturation channels of the images in HSV color space, and the transformation is kept only if the value has increased compared with the previous setup. Using the same correlation values as weights for each pairwise registration, the maximum spanning tree was identified over the complete graph of the sections.

Finally, the images were transformed into a joint space by combining the registrations over the paths to the section with the highest degree in the spanning tree. The accuracy of the automatic registration was checked manually, and 11 low-quality spots were discarded. The representative colors of the nuclear and antibody stains were manually identified separately for each section and used for performing color deconvolution on the RGB images (36). The colocalization of the stains was quantified with Manders' split coefficients, and their significance was assessed with Costes' P value on the registered and deconvolved images (37, 38). The full TMA image analysis workflow except figure creation was performed in Matlab (R2020a v9.8, RRID:SCR_001622) and the source code of our image analysis workflow is freely available at https://bitbucket.org/szkabel/pirhonen-et-al-2021.

Data availability

The proteomic data set generated in this study is available at ProteomeXChange at: http://proteomecentral.proteomexchange.org/cgi/GetDataset?ID=PXD026371.

The diagnostic slides including the digital segmentation of the four tissue types is available at: https://doi.org/10.23729/2a5e6d68-8f2e-4b93-b0e0-cedf52c3291c.

Localization of cholesterol biosynthesis in the neoplastic parenchyma

We initially performed IHC stainings of tissue microarray slides for major cholesterol biosynthesis and uptake-related proteins in a cohort of 24 PDAC patients (Table 1). Parallel slides were stained with antibodies against 7-dehydrocholesterol reductase (DHCR7, the last protein in the Kandutsch-Russell pathway for cholesterol biosynthesis; Supplementary Fig. S1A), LDL receptor as well as cytokeratin 7 (CK7, tumor marker) and SMA (stromal marker), and the images were overlaid using a robust image registration framework (Supplementary Fig. S1B). Anti-LDLR positive areas were sparse, with the strongest signals observed in areas morphologically detected as Langerhans islets (Supplementary Fig. S2A). However, the results revealed a prominent DHCR7 immunoreactivity in the CK7-positive neoplastic areas (Supplementary Fig. S2B and S2C). As a specificity control, the overlap of CK7 with SMA was marginal and nonsignificant based on Costes' randomization of colocalization (Supplementary Fig. S2C).

Table 1.

Clinical data of patients involved in proteomics and histologic analyses.

TMAsCryosectionsProteomics
n 24 16 14 
Sex (male), n (%) 13 (54.2) 9 (60.0) 7 (50.0) 
Age (years), mean ± SD 63.64 ± 8.52 68.74 ± 8.10 65.07 ± 8.58 
Survival time (days), mean ± SD 954.12 ± 897.17 1,044.07 ± 1,016.86 830.36 ± 649.14 
Hospital stay (days), mean ± SD 16.08 ± 20.87 10.50 ± 4.95 11.00 ± 3.72 
Postoperative stay (days), mean ± SD) 14.71 ± 20.83 12.38 ± 5.77 10.57 ± 4.62 
T, mean ± SD 2.62 ± 0.71 2.93 ± 0.47 3.14 ± 0.36 
N, mean ± SD 0.75 ± 0.44 0.69 ± 0.48 1.00 ± 0.00 
M, mean ± SD 0.12 ± 0.34 0.00 ± 0.00 0.00 ± 0.00 
Grade, n (%) 
 1 2 (8.7) 2 (18.2) 0 (0.0) 
 2 14 (60.9) 8 (72.7) 14 (100.0) 
 3 7 (30.4) 1 (9.1) 0 (0.0) 
Radicality (R1), n (%) 6 (25.0) 4 (36.4) 3 (23.1) 
Size of primary (mm), mean ± SD 37.17 ± 19.12 37.09 ± 18.51 38.57 ± 13.63 
Metastasized nodules, mean ± SD 2.04 ± 2.12 2.85 ± 3.11 3.50 ± 1.99 
Nodules collected, mean ± SD 23.29 ± 10.11 20.77 ± 6.39 20.14 ± 7.69 
Preoperative chemotherapy (true), n (%) 2 (8.3) 1 (7.7) 0 (0.0) 
Bilirubin (μmol/L), mean ± SD 35.12 ± 55.04 26.00 ± 14.82 20.79 ± 13.27 
CA19-9 (IU/L), mean ± SD 2,331.62 ± 7,816.17 292.64 ± 309.79 459.14 ± 412.74 
CEA (μg/L), mean ± SD 4.67 ± 5.60 2.53 ± 2.03 2.79 ± 1.48 
Albumin (g/L), mean ± SD 36.43 ± 5.22 NaN ± NA 38.16 ± 2.44 
WBC (E9/L), mean ± SD 7.57 ± 3.42 NaN ± NA 7.77 ± 2.54 
Thrombocytes (E9/L), mean ± SD 292.87 ± 108.09 NaN ± NA 268.21 ± 57.90 
CRP (mg/L), mean ± SD 6.12 ± 5.56 0.61 ± 0.09 4.00 ± 2.68 
Statin use, n (%) 6 (25.0) 6 (37.5) 1 (7.1) 
TMAsCryosectionsProteomics
n 24 16 14 
Sex (male), n (%) 13 (54.2) 9 (60.0) 7 (50.0) 
Age (years), mean ± SD 63.64 ± 8.52 68.74 ± 8.10 65.07 ± 8.58 
Survival time (days), mean ± SD 954.12 ± 897.17 1,044.07 ± 1,016.86 830.36 ± 649.14 
Hospital stay (days), mean ± SD 16.08 ± 20.87 10.50 ± 4.95 11.00 ± 3.72 
Postoperative stay (days), mean ± SD) 14.71 ± 20.83 12.38 ± 5.77 10.57 ± 4.62 
T, mean ± SD 2.62 ± 0.71 2.93 ± 0.47 3.14 ± 0.36 
N, mean ± SD 0.75 ± 0.44 0.69 ± 0.48 1.00 ± 0.00 
M, mean ± SD 0.12 ± 0.34 0.00 ± 0.00 0.00 ± 0.00 
Grade, n (%) 
 1 2 (8.7) 2 (18.2) 0 (0.0) 
 2 14 (60.9) 8 (72.7) 14 (100.0) 
 3 7 (30.4) 1 (9.1) 0 (0.0) 
Radicality (R1), n (%) 6 (25.0) 4 (36.4) 3 (23.1) 
Size of primary (mm), mean ± SD 37.17 ± 19.12 37.09 ± 18.51 38.57 ± 13.63 
Metastasized nodules, mean ± SD 2.04 ± 2.12 2.85 ± 3.11 3.50 ± 1.99 
Nodules collected, mean ± SD 23.29 ± 10.11 20.77 ± 6.39 20.14 ± 7.69 
Preoperative chemotherapy (true), n (%) 2 (8.3) 1 (7.7) 0 (0.0) 
Bilirubin (μmol/L), mean ± SD 35.12 ± 55.04 26.00 ± 14.82 20.79 ± 13.27 
CA19-9 (IU/L), mean ± SD 2,331.62 ± 7,816.17 292.64 ± 309.79 459.14 ± 412.74 
CEA (μg/L), mean ± SD 4.67 ± 5.60 2.53 ± 2.03 2.79 ± 1.48 
Albumin (g/L), mean ± SD 36.43 ± 5.22 NaN ± NA 38.16 ± 2.44 
WBC (E9/L), mean ± SD 7.57 ± 3.42 NaN ± NA 7.77 ± 2.54 
Thrombocytes (E9/L), mean ± SD 292.87 ± 108.09 NaN ± NA 268.21 ± 57.90 
CRP (mg/L), mean ± SD 6.12 ± 5.56 0.61 ± 0.09 4.00 ± 2.68 
Statin use, n (%) 6 (25.0) 6 (37.5) 1 (7.1) 

Lipid deposition in pancreatic cancer samples

Considering the reported lipid storage in pancreatic cancers observed by stimulated Raman scattering microscopy (7), we investigated potential lipid accumulation in fresh cryosections of 16 operable PDAC patients (Table 1) using nonlinear label-free imaging. The sections were first imaged with CARS microscopy for the detection of lipids and second harmonic generation microscopy for a structural overview provided by fibrillar collagens (Supplementary Fig. S3A), followed by conventional hematoxylin and eosin stainings (Supplementary Fig. S3B). We have previously used this strategy to detect lipid accumulation and fibrosis in the human fatty liver (39). Prominent CARS signals were obtained, as expected, from the visceral adipose tissue surrounding the exocrine glandular areas (Supplementary Fig. S3A and S3C–S3E). However, both the neoplastic and adjacent parenchymal areas yielded very weak CARS signals (Supplementary Fig. S3A and S3C–S3E). Thus, despite the high cholesterogenic enzyme expression in the neoplastic exocrine pancreas, no major acinar cell lipid accumulation could be detected in the PDAC samples. Of note, due to the limited sensitivity of CARS imaging, this does not exclude the presence of minor subcellular lipid deposition in these areas.

Laser-capture microdissection and high-resolution spatial proteomics

To obtain a more holistic view of the metabolic landscape in the malignant and neighboring tissue regions, we carried out untargeted, high-resolution proteomics of spatially microdissected PDAC specimens. We screened 40 PDAC patients' routine diagnostic histologic slides to identify cases suitable for the study. A set of 14 patients with a confirmed diagnosis, available clinical data (Table 1), and sufficient tissue, especially morphologically benign parenchyma around the cancerous areas, were selected. None of the patients had received preoperative chemotherapy that may alter tissue architecture and protein composition. Standard diagnostic formalin-fixed paraffin-embedded blocks of the tumors were first sectioned to 10 μm (batch 1) or 20 μm (batch 2) thickness and stained with hematoxylin and eosin. The sections were imaged, and the areas of interest were digitally segmented by a pathologist for microdissection. The corresponding segmentations were then superimposed on the laser-microdissection platform (Fig. 1A).

Figure 1.

Overview of the proteomic analysis. A, Overview of sample processing steps. B, An exemplary sample for microdissection. White boxes present insets displayed in the right bottom corner of the image. Line colors represent different tissue types: blue, adjacent parenchyma; cyan, adjacent stroma; red, neoplastic parenchyma; orange, neoplastic stroma. C, Overview of the data processing steps. Blue, protein-level data; red, pathway-level data. N indicates the number of proteins/pathways present in the step. The numbers and letters in parentheses in the Outcome column indicate the figure and the panel where the plot is shown. D, PCA plot of log-median normalized and batch effect–corrected protein abundances. Tissue types: AP, adjacent parenchyma; AS, adjacent stroma; NS, neoplastic stroma; NP, neoplastic parenchyma. Individual patients, P1–P14. H&E, hematoxylin and eosin.

Figure 1.

Overview of the proteomic analysis. A, Overview of sample processing steps. B, An exemplary sample for microdissection. White boxes present insets displayed in the right bottom corner of the image. Line colors represent different tissue types: blue, adjacent parenchyma; cyan, adjacent stroma; red, neoplastic parenchyma; orange, neoplastic stroma. C, Overview of the data processing steps. Blue, protein-level data; red, pathway-level data. N indicates the number of proteins/pathways present in the step. The numbers and letters in parentheses in the Outcome column indicate the figure and the panel where the plot is shown. D, PCA plot of log-median normalized and batch effect–corrected protein abundances. Tissue types: AP, adjacent parenchyma; AS, adjacent stroma; NS, neoplastic stroma; NP, neoplastic parenchyma. Individual patients, P1–P14. H&E, hematoxylin and eosin.

Close modal

The neoplastic cells typically formed irregular ductal structures or hyperchromatic cell clusters embedded in the connective tissue. Such cell clusters, as well as regions of the intervening stroma, were separately collected to represent the neoplastic parenchymal and neoplastic stromal samples (NP, NS; Fig. 1B). Histologically normal-appearing exocrine pancreatic tissue and its neighboring stromal areas were additionally collected from each specimen to identify molecular changes reaching the visually unaffected regions. We term these samples as (tumor) adjacent parenchyma and adjacent stroma (AP, AS; Fig. 1B). To obtain good proteomic coverage and to limit tissue damage in laser microdissection, a few hundred cells per cluster were collected. The total tissue volume collected per sample averaged 30 × 106 μm³ (corresponding to a cube with an edge of 311 μm).

Deep mass spectrometric (MS) analysis of the samples identified a total of 6,979 unique proteins, representing a 72% coverage of pancreatic reference proteomes (Supplementary Fig. S4A–S4C). The raw protein abundances were first log-median normalized and corrected for batch effect after excluding low coverage observations (Fig. 1C). The code used and a more detailed description of the analysis are openly accessible at: https://bitbucket.org/szkabel/pirhonen-et-al-2021.

Distinct proteomic signatures in the microdissected samples

The 3,127 proteins available from all the samples formed three distinct clusters in PCA: neoplastic parenchyma, adjacent parenchyma, and a combination of the stromas (Fig. 1D). The first PCA dimension differentiated the parenchymal samples from the stromal ones and the second dimension separated the neoplastic parenchyma from the other tissue compartments. Together, these two dimensions accounted for about one third of all the variance in the proteomic data.

A pairwise comparison of the compartmental proteomes revealed that typical exocrine pancreas-specific proteins, such as alpha amylase (AMY2A) and chymotrypsin-like elastase 3A (CELA3A), were more abundant in the adjacent than in the neoplastic regions, in line with dedifferentiation of the histologically malignant areas. On the other hand, keratins, such as keratin 18 and 19, were enriched in the neoplastic parenchyma (Fig. 2A and B). Connective tissue-related proteins, such as type I collagen (COL1A2), fibril-associated collagen type XIV (COLA14A1), and fibrillin-1 (FBN1), were expectedly more prevalent in the stromal than parenchymal regions (Fig. 2B; Supplementary Fig. S4D and S4E), demonstrating the desmoplastic reaction. Interactive graphs of the data to search, annotate, and compare the distributions of individual proteins between tissue compartments (NP vs. NS 6,300 proteins; AP vs. AS 6,277 proteins; NP vs. AP 6,613 proteins; NS vs. AS 5,997 proteins) are available as Supplementary File S3.

Figure 2.

Exploration of tissue compartmental proteomes. A and B, Illustration of the interactive data available as Supplementary File S3. Differential protein expression between adjacent and neoplastic parenchyma (A), and between neoplastic stroma and parenchyma (B). C, Heatmap of the pathway enrichment scores after batch correction and pruning. Samples are clustered with k-means = 3, as suggested by the PCA plot, whereas pathways are organized with hierarchical clustering. Tissue type and patient identification as in Fig. 1.

Figure 2.

Exploration of tissue compartmental proteomes. A and B, Illustration of the interactive data available as Supplementary File S3. Differential protein expression between adjacent and neoplastic parenchyma (A), and between neoplastic stroma and parenchyma (B). C, Heatmap of the pathway enrichment scores after batch correction and pruning. Samples are clustered with k-means = 3, as suggested by the PCA plot, whereas pathways are organized with hierarchical clustering. Tissue type and patient identification as in Fig. 1.

Close modal

Tissue compartmental proteomes reveal differential pathway expressions

To better understand pathway-level differences between the tissue types, we focused on Reactome pathways based on their protein coverage: in the pathways included, at least 50% of the expected proteins were present in our data set. This resulted in the estimation of pathway enrichment for 1,095 Reactome pathways. Pairwise comparisons of the tissue compartments at the pathway level are exemplified in Supplementary Fig. S4F and S4G. All pairwise comparisons are shared as interactive graphs in Supplementary File S4 that can be queried similarly as Supplementary File S3.

We hypothesized that despite patient-based variation, pruning the number of pathways might reveal general biological features of the tissue compartments shared between patients. Indeed, the initial PCA constructed using 1,095 pathways was essentially reconstituted using a robust sparse PCA algorithm with only 55 pathways (Fig. 1C; Supplementary Fig. S5A and S5B; Supplementary Table S1). A heatmap demonstrating the expressions of these pathways in the four tissue compartments was constructed (Fig. 2C). This revealed that the pathway “Digestion of dietary lipid” as well as “Cobalamin transport and metabolism” (that contains shared proteins with “Digestion of dietary lipid”) were heavily downregulated in all neoplastic parenchymal regions compared with the adjacent parenchymal ones. Moreover, pathways related to rRNA processing and protein translation were also reduced. Considering the high digestive enzyme synthesis and secretory activity typical for the exocrine pancreas, these results agree with the loss of normal pancreatic functions in the cancerous areas. This may be attributed to the dedifferentiation of acinar cells and/or proliferation of neoplastic ductal cells producing fewer digestive enzymes than the adjacent acinar regions. The two stromal compartments were more alike and distinguished from the parenchymas by vasculature-related functions, such as platelet activation (GP1b-IX-V activation signaling) and coagulation, scavenger receptor activity, as well as connective tissue remodeling–related pathways (Fig. 2C). Additionally, although the “Cellular response to hypoxia” Reactome pathway was not included in the pruned set identified by sparse PCA, it was more enriched in the parenchymal than in stromal samples (Supplementary Fig. S4F). This speaks for relative hypoxia in the parenchymal regions compared with the stromal areas.

Lipid metabolism and transport proteins

To focus on lipid-related metabolism, we visualized the differences in protein expressions between parenchymal and stromal compartments for Reactome lipid-related pathways, such as “Fatty acid metabolism” or “Triglyceride metabolism.” Those “Lipid metabolism” pathways that displayed differential expression between the two compartments all showed an enrichment in the parenchymal tissues, both neoplastic and adjacent ones, relative to the stromas (Fig. 3A). In contrast, for “Plasma lipoprotein assembly, remodeling, and clearance” pathways, neoplastic and adjacent stromal samples showed higher enrichments than the corresponding parenchymal ones (Fig. 3B).

Figure 3.

Lipid-related Reactome pathway enrichments and individual protein abundances. A and B, Reactome pathway enrichment differences between stroma and parenchyma. The center node represents a Reactome pathway (A, Lipid metabolism; B, Plasma lipoprotein assembly, remodeling and clearance) and nodes connected to it represent the two underlying hierarchical levels. On each level, dots closer to or further from the center of the circle compare neoplastic or adjacent tissues, respectively. The color indicates the difference in pathway expression between parenchymal and stromal tissues. For example, fatty acid metabolism and triglyceride metabolism were more active in both of the parenchymal regions compared with their neighboring stromas. Gray dot and text indicate that there was no significant representation for a pathway. Pooled data from 14 patients. C, Log median-normalized abundances of proteins from the cholesterol biosynthesis Reactome pathway (top) and manually selected proteins related to lipoprotein uptake (bottom) in patients from batch 2. Lipoprotein receptors checked for but not detected in our data are ABCA1, GPIHBP1, LOX1, LRP1B, LRP2, LRP4, LRP5, LRP6, LRP8, SCARB1, SCARF1, STAB2, and VLDLR. Tissue type and patient identification as in Fig. 1.

Figure 3.

Lipid-related Reactome pathway enrichments and individual protein abundances. A and B, Reactome pathway enrichment differences between stroma and parenchyma. The center node represents a Reactome pathway (A, Lipid metabolism; B, Plasma lipoprotein assembly, remodeling and clearance) and nodes connected to it represent the two underlying hierarchical levels. On each level, dots closer to or further from the center of the circle compare neoplastic or adjacent tissues, respectively. The color indicates the difference in pathway expression between parenchymal and stromal tissues. For example, fatty acid metabolism and triglyceride metabolism were more active in both of the parenchymal regions compared with their neighboring stromas. Gray dot and text indicate that there was no significant representation for a pathway. Pooled data from 14 patients. C, Log median-normalized abundances of proteins from the cholesterol biosynthesis Reactome pathway (top) and manually selected proteins related to lipoprotein uptake (bottom) in patients from batch 2. Lipoprotein receptors checked for but not detected in our data are ABCA1, GPIHBP1, LOX1, LRP1B, LRP2, LRP4, LRP5, LRP6, LRP8, SCARB1, SCARF1, STAB2, and VLDLR. Tissue type and patient identification as in Fig. 1.

Close modal

To further scrutinize lipid metabolic parenchyma–stroma cross-talk, we generated a heatmap of protein abundances related to cholesterol biosynthesis and uptake. The protein abundances in the “Cholesterol biosynthesis” Reactome pathway confirmed that the neoplastic parenchyma indeed exhibited the highest cholesterol biosynthetic activity across the examined tissue types (Fig. 3C), confirming and extending our DHCR7 IHC results (Supplementary Fig. S2). Regarding cholesterol uptake, the most prominent protein in the LDLR family was not LDLR itself, but LDLR-related protein 1 (LRP1) that was abundant throughout the four tissue compartments (Fig. 3C). This matches with the relatively weak LDLR immunoreactivity in the histochemical stainings. Of the high-density lipoprotein (HDL) receptor family, only one protein was identified: the multifunctional protein known as high-density lipoprotein binding protein (HDLBP)/vigilin that showed the highest abundance in the adjacent parenchymal samples (Fig. 3C).

To assess the parenchyma–stroma compartmentalization of lipid metabolism and transport in an independent data set, we performed a matching pathway analysis using microdissected pancreatic cancer and stromal proteomes recently reported by Le Large and coworkers (28). The proteomic coverage of this data set was somewhat smaller than ours (Supplementary Fig. S4B and S4C) and separate samples representing exocrine and stromal components adjacent to the neoplastic regions were not collected. However, the malignant parenchymal and their neighboring stromal proteomes were comparable with ours (Supplementary Fig. S6A–S6C). Importantly, we found that the pathway analysis of the Le Large and colleagues data reproduced the pattern observed in our samples: an overwhelming dominance of the neoplastic parenchyma over stroma in lipid metabolic pathway protein expressions (Supplementary Fig. S6D) and higher plasma lipoprotein transport activity in the neoplastic stroma compared with the neoplastic parenchyma (Supplementary Fig. S6E).

Prognostic features in pancreatic cancer

We next analyzed possible prognostic factors in our proteomic data set, by dividing the 14 patients into two groups based on postoperative survival over or under two years (7 and 7 patients) and plotting the pathways most prognostic for survival. In the neoplastic areas, both the parenchyma and stroma, pathways related to O2 and CO2 exchange, and thereby most likely reflecting the presence of capillaries, contributed toward a worse prognosis. The most significantly prognostic pathway was, however, that of “Plasma lipoprotein clearance” and, somewhat unexpectedly, this effect was linked to the histologically benign adjacent parenchyma (Fig. 4A). In fact, this tissue compartment harbored overall the highest number of pathways significant for survival (Fig. 4B).

Figure 4.

Protein abundances and pathway expressions are most significant for survival in adjacent parenchyma. A, Prognostic effects of the most significant survival predicting pathways. Each bar shows the difference between the medians of pathway expression in long (over 2 years) and short (under 2 years) surviving patients. Color represents the significance of the separation between the groups. B, Number of prognostic (P < 5%) pathways in each tissue type. Colors show the distribution of these pathways across root nodes in Reactome. C, Protein-level volcano plot of prognostic effects in the adjacent parenchyma. The x-axis indicates the difference in median protein abundance between patients surviving over and under 2 years. Proteins enriched in patients with short survival have negative values. D, Representative images showing IHC stainings of NCEH1 (top) and AP2M1 (bottom). Top left shows an overview of a larger area, and white boxes indicate the regions shown as insets in the bottom row. Colored regions indicate examples of ROIs used for the analysis: red, neoplastic parenchyma; blue, adjacent parenchyma. Top right shows a representative region from a separate control sample. E, Quantification of the IHC stainings in D, measured as the average color-deconvolved brown channel intensity over the defined ROIs. n = 14 patients for neoplastic and adjacent parenchyma and n = 7 for control parenchyma (6–50 ROI/sample). Wilcoxon signed-rank test was performed and P values are indicated at the top. F, Each point represents an adjacent parenchymal ROI (n = 70, 90 for NCEH1 and AP2M1, respectively). The plot shows the average intensity of the color-deconvolved brown channel over the ROI plotted against its distance to the closest neoplastic parenchymal ROI (distance measured between the centroids of the ROIs).

Figure 4.

Protein abundances and pathway expressions are most significant for survival in adjacent parenchyma. A, Prognostic effects of the most significant survival predicting pathways. Each bar shows the difference between the medians of pathway expression in long (over 2 years) and short (under 2 years) surviving patients. Color represents the significance of the separation between the groups. B, Number of prognostic (P < 5%) pathways in each tissue type. Colors show the distribution of these pathways across root nodes in Reactome. C, Protein-level volcano plot of prognostic effects in the adjacent parenchyma. The x-axis indicates the difference in median protein abundance between patients surviving over and under 2 years. Proteins enriched in patients with short survival have negative values. D, Representative images showing IHC stainings of NCEH1 (top) and AP2M1 (bottom). Top left shows an overview of a larger area, and white boxes indicate the regions shown as insets in the bottom row. Colored regions indicate examples of ROIs used for the analysis: red, neoplastic parenchyma; blue, adjacent parenchyma. Top right shows a representative region from a separate control sample. E, Quantification of the IHC stainings in D, measured as the average color-deconvolved brown channel intensity over the defined ROIs. n = 14 patients for neoplastic and adjacent parenchyma and n = 7 for control parenchyma (6–50 ROI/sample). Wilcoxon signed-rank test was performed and P values are indicated at the top. F, Each point represents an adjacent parenchymal ROI (n = 70, 90 for NCEH1 and AP2M1, respectively). The plot shows the average intensity of the color-deconvolved brown channel over the ROI plotted against its distance to the closest neoplastic parenchymal ROI (distance measured between the centroids of the ROIs).

Close modal

To scrutinize the prognostic features in the adjacent parenchyma, we first analyzed the association of individual proteins with survival in this tissue type (Fig. 4C). The protein showing the strongest predisposing individual effect was serine and arginine-rich splicing factor 4 (SRSF4; Fig. 4C). This suggests that alterations in the splicing machinery may contribute to poor prognosis in PDAC, in line with recent studies reporting the association of alternative splicing with bad prognosis in PDAC (40, 41). The top-scoring hits from the “Plasma lipoprotein clearance” pathway included adaptor-related protein complex 2 subunit mu (AP2M1) and neutral cholesterol ester hydrolase 1 (NCEH1). The former plays a role in the endocytic uptake of cargo, such as lipoprotein receptors, whereas the latter hydrolyzes cellular cholesterol esters into free, i.e., unesterified cholesterol, thus leading to cholesterol relocation from lipid droplets to membranes. Together, these data suggest that the tumor-adjacent parenchymal tissue contains molecular changes important for disease prognosis and based on the pathway-level evidence, lipid metabolism plays a role in this effect.

Distribution of select prognosis-related proteins by IHC

We next assessed the distribution of two top-scoring proteins from the “Plasma lipoprotein clearance” pathway, i.e., AP2M1 and NCEH1, by IHC staining of the 14 PDAC patient samples used for the proteomic analysis. For comparison, we also included seven pancreatic samples that histologically exhibited no signs of malignancy. The signals of both NCEH1 and AP2M1 immunoreactivity were most pronounced in the neoplastic regions (Fig. 4D and E). There was also considerable reactivity in the adjacent parenchymal areas, especially in the vicinity of the tumor. However, the overall signal in the regions furthest away from the tumor was comparable to those of the control samples, although in the case of NCEH1, small patches of cells in the tumor-adjacent parenchyma were often positive (Fig. 4D and E). Correspondingly, for both NCEH1 and AP2M1, immunoreactivity in the tumor-adjacent regions decreased with increasing the distance from the overtly neoplastic areas (Fig. 4F), suggesting that these changes may somehow be related to the malignant transformation process.

Verification of findings in external data sets

As our small patient cohort does not allow firm conclusions regarding the prognostic significance of markers, we also surveyed other PDAC proteomic studies in the literature. A schematic overview of all the proteomic data sets used to verify our findings is provided in Fig. 5A. The recently published resource article on pancreatic cancer by Cao and colleagues (29) contained 37 additional PDAC patients with sufficient clinical information, neoplastic tissue cellularity, and paired normal adjacent tissue for comparison with our findings. We observed a reasonable correlation of protein abundances between Cao and colleagues (29) and our samples (Supplementary Fig. S7A) and similar differential expressions in neoplastic versus adjacent tissues, supporting the loss of normal digestive functions in the neoplastic areas (Supplementary Fig. S7B). Most importantly, also in the cohort of Cao and colleagues (29), the adjacent tissue harbored the highest number of pathways significant for survival (Fig. 5B) and increased expression of AP2M1, ESYT1, and NCEH1 in this tissue (rather than in the neoplastic area) associated with decreased survival (Fig. 5C). Altogether, these data argue that changes in the tumor-adjacent exocrine tissue are critical determinants of survival in pancreatic cancer.

Figure 5.

Overview of related studies used to verify findings of this study. A, Schematic figure depicting the origin of samples in the studies included in this paper. Le Large et al. (28) and this study applied laser microdissection while the others used bulk samples. On the right, the figure numbers in parentheses indicate where in this paper data supporting the finding can be found. We integrate our samples with control pancreata from Kushner et al. (42) in the subsequent parts of the paper. B, Number of prognostic pathways (P < 5%) in each tissue type in the Cao et al. study (29). Colors show the distribution of these pathways across root nodes in Reactome. C, Kaplan–Meier survival plots on select lipid-related protein abundances in the Cao et al. study (29).

Figure 5.

Overview of related studies used to verify findings of this study. A, Schematic figure depicting the origin of samples in the studies included in this paper. Le Large et al. (28) and this study applied laser microdissection while the others used bulk samples. On the right, the figure numbers in parentheses indicate where in this paper data supporting the finding can be found. We integrate our samples with control pancreata from Kushner et al. (42) in the subsequent parts of the paper. B, Number of prognostic pathways (P < 5%) in each tissue type in the Cao et al. study (29). Colors show the distribution of these pathways across root nodes in Reactome. C, Kaplan–Meier survival plots on select lipid-related protein abundances in the Cao et al. study (29).

Close modal

Tumor-adjacent exocrine pancreas differs from healthy pancreas

To investigate how the tumor-adjacent parenchyma may differ from a nondiseased pancreas, we reprocessed proteomics raw data from 12 healthy control human pancreata (42) and compared these to our data set (Fig. 5A). In general, these healthy control samples had fewer detected proteins of lower total abundance than our samples. To enable comparisons between them, we used two alternative strategies: normalization against housekeeping, i.e., most constantly expressed, proteins in an independent proteomics data set of multiple tissues (24) and GFA that relies on matrix factorization of minimally processed data (Supplementary Fig. S8A).

Following the conventional housekeeping proteome normalization, PCA of the 12 control pancreatic proteomes and our samples clearly separated the control pancreata from the other tissue types (Supplementary Fig. S8A and S8B). On PCA dimension 1, the control pancreata were closest to the adjacent parenchyma of our data set, whereas on dimension 2, they were closer to the stromal samples, in keeping with the fact that stromal areas were not separately dissected in the control samples. To further study how the control pancreata differed from our samples, we constructed a sparse PCA heatmap. This showed that the adjacent parenchyma had retained a similar abundance of pancreatic enzymes (such as CELA3A and AMY2A) as the healthy pancreas (Supplementary Fig. S8C). However, regarding some proteins, the adjacent parenchyma differed from the healthy control samples: for example, select cytoskeletal proteins, such as keratins 9 and 10, were more abundant, whereas others, such as lamin A/C (LMNA), were less abundant in the adjacent parenchyma (Supplementary Fig. S8C).

GFA highlights factors differentiating most between tumor-adjacent and healthy pancreas

Similar to PCA, GFA decomposes the data into components (factors) whose sum reconstructs the data, but enables a more in-depth exploration, as missing values are not excluded, and no batch-effect correction is conducted. Hence, the first two factors in GFA (according to decreasing effect size) captured batch-based variation (Fig. 6A and B), whereas later factors revealed biologically meaningful properties of the data. Factor 3 captured similarities between control and adjacent parenchymal tissues and explained this similarity with retained healthy pancreatic functions as represented by high abundance of pancreatic proteins and the “Pancreatic secretion” pathway (Fig. 6C). Factor 4 separated the neoplastic parenchyma from other tissues and explained the difference with the spliceosome Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, emphasizing the importance of splicing differences also at the pathway level (Fig. 6D).

Figure 6.

GFA of the proteomics data identifies key differences between tissue types. A, The projection matrix WT visualizes similarities and differences between the tissue types. Columns are samples, rows are factors identified by GFA, and the colors indicate the value of the weight. The first two factors show bimodal distributions according to batches in our data, followed by factors that capture tissue-specific characteristics of the data. B, Factor weights for factors 1 to 5, grouped by tissue type and colored by the sample batch. Weights for factors 1 and 2 show bimodal distribution according to the batches, whereas later factors are aligned within a tissue type. C–E, Interpretation of the factors 3 to 5. A factor's overall contribution to a sample is the product of the latent variable and the factor's weight. Consequently, the most extreme latent variable values determine the proteins that are the most significant for the definition of the factor. Left, latent variable values with the most extreme proteins from the top-scoring pathways (shown on the right) are highlighted. Note that in E, LMNA and KRT9 are shown as a comparison with housekeeping normalization; they are not part of the pathways from the right. Middle, factor weights showing the effect direction of the factor for tissue types. For instance, in factor 3 (C), both the adjacent parenchyma and the control parenchyma have positive weight, meaning that the proteins with positive latent variable (e.g., pancreatic lipase: PNLIP) were highly abundant in those samples. Conversely, negative factor weights indicate a high abundance of proteins with a negative latent variable value, for instance, in factor 5, apolipoproteins were more abundant in the tumor samples compared with control. Right (text next to arrows), pathway-level analysis of the factor latent variables. The list includes in both directions the top-scoring three KEGG pathways according to the STRING enrichment score for factors 3 and 4, and the single significantly enriched annotated UniProt keyword for factor 5. Rightmost part of E, Kaplan–Meier survival plot based on the weights of factor 5 in adjacent parenchyma. Tissue type and patient identification as in Fig. 1 and additionally control pancreas (CP).

Figure 6.

GFA of the proteomics data identifies key differences between tissue types. A, The projection matrix WT visualizes similarities and differences between the tissue types. Columns are samples, rows are factors identified by GFA, and the colors indicate the value of the weight. The first two factors show bimodal distributions according to batches in our data, followed by factors that capture tissue-specific characteristics of the data. B, Factor weights for factors 1 to 5, grouped by tissue type and colored by the sample batch. Weights for factors 1 and 2 show bimodal distribution according to the batches, whereas later factors are aligned within a tissue type. C–E, Interpretation of the factors 3 to 5. A factor's overall contribution to a sample is the product of the latent variable and the factor's weight. Consequently, the most extreme latent variable values determine the proteins that are the most significant for the definition of the factor. Left, latent variable values with the most extreme proteins from the top-scoring pathways (shown on the right) are highlighted. Note that in E, LMNA and KRT9 are shown as a comparison with housekeeping normalization; they are not part of the pathways from the right. Middle, factor weights showing the effect direction of the factor for tissue types. For instance, in factor 3 (C), both the adjacent parenchyma and the control parenchyma have positive weight, meaning that the proteins with positive latent variable (e.g., pancreatic lipase: PNLIP) were highly abundant in those samples. Conversely, negative factor weights indicate a high abundance of proteins with a negative latent variable value, for instance, in factor 5, apolipoproteins were more abundant in the tumor samples compared with control. Right (text next to arrows), pathway-level analysis of the factor latent variables. The list includes in both directions the top-scoring three KEGG pathways according to the STRING enrichment score for factors 3 and 4, and the single significantly enriched annotated UniProt keyword for factor 5. Rightmost part of E, Kaplan–Meier survival plot based on the weights of factor 5 in adjacent parenchyma. Tissue type and patient identification as in Fig. 1 and additionally control pancreas (CP).

Close modal

Importantly, factor 5 revealed differences between the tumor-adjacent and healthy control samples. First, proteins related to neuroendocrine secretion [such as chromogranin A (CHGA), neurosecretory protein VGF (VGF) and glucagon (GCG)] were retained in the healthy controls but missing from our samples (Fig. 6E). This agrees with the fact that Langerhans islets were not included in our microdissected parenchymal samples. Moreover, proteins such as Sorbin and SH3 domain-containing protein 2 (SORBS2, a.k.a. ArgBP2) and LMNA (also observed by PCA) were markedly downregulated in the adjacent parenchymas compared with controls (Fig. 6E). Interestingly, ArgBP2 was reported to be repressed during the oncogenic transformation of the pancreas and linked to the control of cell adhesion and migration (43). In contrast, lipid transport–related proteins, especially extended synaptotagmin 1 (ESYT1) and several apolipoproteins, were more abundant in the tumor-adjacent than in healthy control parenchyma (Figs. 6E and 7A). Of note, ESYT1, implicated in intracellular lipid transport via membrane contact sites (44), also showed a tendency for enrichment in patients with a short survival in our patients (Fig. 4C).

Figure 7.

Lipid metabolic landscape of pancreatic cancer. A, Heatmap of batch effect–corrected and housekeeping-normalized protein abundances for select plasma lipoproteins and cellular lipid-related proteins differentiating the tumor and control samples. Tissue type and patient identification as in Fig. 1. Note that the stromal abundance of apolipoproteins is also recapitulated by canonical plasma proteins (as exemplified by fibrinogens FGA, FGB, and FGG), suggesting that the apolipoproteins originate from blood supply rather than stromal synthesis. B, Cartoon summarizing the main findings of this study. Dark purple, exocrine cells; light purple, connective tissue cells. Top left, healthy control pancreas, with intact pancreatic secretion and cellular architecture. Middle, tumor-adjacent parenchyma and stroma in pancreatic cancer with still intact morphology but increased apolipoprotein supply from blood, presence of LRP1, and increased intracellular lipid transport proteins. Right, dedifferentiated neoplastic parenchyma and fibrotic stroma in pancreatic cancer. Besides high lipid uptake, lipid (cholesterol) biosynthesis is also upregulated.

Figure 7.

Lipid metabolic landscape of pancreatic cancer. A, Heatmap of batch effect–corrected and housekeeping-normalized protein abundances for select plasma lipoproteins and cellular lipid-related proteins differentiating the tumor and control samples. Tissue type and patient identification as in Fig. 1. Note that the stromal abundance of apolipoproteins is also recapitulated by canonical plasma proteins (as exemplified by fibrinogens FGA, FGB, and FGG), suggesting that the apolipoproteins originate from blood supply rather than stromal synthesis. B, Cartoon summarizing the main findings of this study. Dark purple, exocrine cells; light purple, connective tissue cells. Top left, healthy control pancreas, with intact pancreatic secretion and cellular architecture. Middle, tumor-adjacent parenchyma and stroma in pancreatic cancer with still intact morphology but increased apolipoprotein supply from blood, presence of LRP1, and increased intracellular lipid transport proteins. Right, dedifferentiated neoplastic parenchyma and fibrotic stroma in pancreatic cancer. Besides high lipid uptake, lipid (cholesterol) biosynthesis is also upregulated.

Close modal

We, therefore, assessed if the factors identified by GFA were associated with survival by conducting Kaplan–Meier analyses on the weights of the factors. Of the first 6 factors (contributing to >70% of the cumulative effect size), factor 5 in the adjacent parenchyma was associated most significantly with survival (Fig. 6E). Thus, the proteomic signature that classifies factor 5, particularly upregulation of proteins related to lipid transport, not only differentiates between healthy and adjacent parenchyma but also appears as prognostically relevant. In summary, these data show that although the tumor-adjacent tissue largely retains normal exocrine functionalities, its proteomic signature differs from that of a healthy pancreas, and this difference is most systematic for lipid metabolism-related proteins.

Here, we generated the spatially best resolved proteomic data set of early (i.e., deemed operable) human PDAC to date. This enabled systematic characterization of pathways, including those represented by relatively scarce proteins, such as metabolic enzymes, in different tissue compartments. An important aspect was to include histologically benign regions from the tumor margin and to compare these to healthy pancreatic proteomes. Although PDAC is traditionally considered to arise from ductal cells, the mechanisms of PDAC initiation and the cell lineage of origin remain unsettled, with increasing evidence pointing to remarkable plasticity between ductal and acinar cells (45, 46). Interestingly, combined analysis of all the data sets in the present study suggests that the proteomes of the tumor-adjacent exocrine and healthy control ductal regions have rather similar characteristics distinct from those of the neoplastic areas (Supplementary Fig. S9). Hence, the changes in the proteome of the overtly transformed regions dominated over acinar–ductal lineage differences.

The malignant exocrine regions showed a broad downregulation of pathways related to normal pancreatic secretory functions (Fig. 7B). In parallel, these regions exhibited high cholesterol biosynthetic activity, as demonstrated by the upregulation of the mevalonate pathway enzymes. Recently, SOAT1 was identified as a key factor maintaining mevalonate pathway hyperactivity in PDAC organoid and mouse models (47). However, in our clinical PDAC patients, SOAT1 levels were lower than those of NCEH1, the enzyme catalyzing the reverse reaction. Accordingly, no prominent deposition of cholesteryl esters in lipid droplets was observed by imaging of the neoplastic areas. Although enzyme levels alone do not dictate pathway activities, it appears that in early human PDAC, the neoplastic cells aim to maximize membrane cholesterol availability. This is likely needed for increased membrane biogenesis of dividing cells and for promoting lipid raft-dependent functions, such as angiogenesis, adhesion, and migration (48).

When comparing lipid metabolic activities between the exocrine and stromal regions based on protein representations at the pathway level, it was evident that both the neoplastic and adjacent parenchymal areas dominated in lipid metabolic processes over the neighboring stromas. On the other hand, plasma lipoprotein–related pathways were more abundant in the stromal areas. We confirmed these trends in the currently only available other human PDAC tumor–stroma proteomic data set of Le Large and colleagues (28). The prominent plasma lipoprotein-related pathways and apolipoproteins, such as ApoAI, ApoB, and ApoE, in the stroma most likely reflect the presence of microvasculature (Fig. 7A). This suggests that the neoplastic cells in operable PDAC are relatively well perfused, offering nutrients and extravasation routes, but also access for therapies. The expression of canonical lipoprotein receptors, such as LRP1 (also known as ApoE receptor), suggests efficient lipoprotein uptake throughout tissue compartments. Instead, what the functional role of the abundant RNA-binding protein vigilin, a.k.a. HDL-binding protein in the pancreas, is remains to be addressed. In the liver, it was reported to modulate ApoB translation (49).

The adjacent parenchymal regions showed distinct differences to healthy control pancreata, indicating more widespread pathologic changes in the tumor margin than histologically evident. Remarkably, the adjacent parenchyma contained most of the changes relevant for disease prognosis. Based on the pathway-level evidence, altered lipid metabolism contributed most prominently to this effect, with lipid-binding proteins such as ESYT1 and NCEH1 being among the upregulated ones and also prognostic in the independent data set of Cao and colleagues (29). This implies that besides enhanced synthesis and uptake of lipids, intracellular lipid handling is also modified in the affected tissues. In line with our findings, it is of interest that a high NCEH1 mRNA level is unfavorable for survival in 176 pancreatic cancer patients from TCGA (P < 0.001; www.proteinatlas.org) and identified as one of the prognostic biomarkers in pancreatic cancer based on bioinformatics analysis (50). Together, these findings suggest that by exploiting liabilities created by altered lipid metabolism, combinatorial strategies targeting not only lipid synthesis and uptake but also intracellular transport and storage may, in the future, offer additional options for intervention (51).

Finally, our findings raise the interesting question of the functional significance of the proteomic changes found in the tumor margin, i.e., whether they represent adaptations of the exocrine tissue to cues deriving from the neoplastic areas and/or possibly early signs of a progressively malignant phenotype in a pancreatic system level disorder. In the latter scenario, the changes found in the adjacent parenchyma might precede those observed in the neoplastic regions. In any case, it is critical to note that in operable PDAC, not only the surrounding stroma but also exocrine pancreatic tissue at more distant sites is affected and exhibits prognostic molecular changes. This urges us to more thoroughly investigate the close-to-normal tissue areas removed in radical tumor resections, taking advantage of the rapid development of microdissection omics approaches, to learn about therapeutically amenable targets for adjuvant therapy at the time of operation.

Á Szkalisity reports grants from the Magnus Ehrnrooth Foundation during the conduct of the study. E. Ikonen reports grants from the Academy of Finland, the Sigrid Juselius Foundation, Fondation Leducq, and Jane and Aatos Erkko Foundation during the conduct of the study. No disclosures were reported by the other authors.

J. Pirhonen: Data curation, software, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. Á. Szkalisity: Data curation, software, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. J.K. Hagström: Resources, formal analysis. Y. Kim: Data curation, investigation, methodology, writing–original draft. E. Migh: Investigation, methodology, writing–original draft. M. Kovács: Investigation, methodology. M. Hölttä: Conceptualization. J. Peränen: Methodology. H. Seppänen: Resources. C. Haglund: Resources. J. Gil: Investigation, methodology. M. Rezeli: Investigation, methodology. J. Malm: Resources. P. Horvath: Resources. Gy. Markó-Varga: Resources. P. Puolakkainen: Conceptualization, resources. E. Ikonen: Conceptualization, resources, supervision, writing–original draft, project administration, writing–review and editing.

The authors thank Anna Uro and Pia Saarinen for excellent technical assistance, Digital Pathology services of the Helsinki Biobank for slide scanning and Biomedicum Imaging Unit, supported by Biocenter Finland and HiLIFE, for infrastructure support. The authors wish to thank the Finnish Grid and Cloud Infrastructure (FGCI) for supporting this project with computational resources. They thank Thermo Fisher Scientific for their generous support. This work was done under the auspices of a Memorandum of Understanding between the European Cancer Moonshot Center in Lund and the U.S. National Cancer Institute's International Cancer Proteogenome Consortium (ICPC). ICPC encourages international cooperation among institutions and nations in proteogenomic cancer research in which proteogenomic datasets are made available to the public. This work was also done in collaboration with the U.S. National Cancer Institute's Clinical Proteomic Tumor Analysis Consortium (CPTAC). This study was supported by the Academy of Finland (307415 and 324929 to E. Ikonen), University of Helsinki (HiLIFE Fellowship and Centre of Excellence matching funds, E. Ikonen), Sigrid Juselius Foundation (E. Ikonen and H. Seppänen), Fondation Leducq (E. Ikonen), Jane and Aatos Erkko Foundation (E. Ikonen), Emil Aaltonen Foundation (J. Pirhonen), Orion Research Foundation (J. Pirhonen), Ida Montin's Foundation (J. Pirhonen), Magnus Ehrnrooth Foundation (Á. Szkalisity), Cancer Foundation Finland (H. Seppänen and J. Pirhonen), Helsinki University Hospital Research Funds (P. Puolakkainen and H. Seppänen), Mary and Georg Ehrnrooth Foundation (H. Seppänen), LENDULET-BIOMAG Grant (2018-342 to P. Horvath), European Regional Development Funds (GINOP-2.3.2-15201600006, GINOP-2.3.215201600026, GINOP-2.3.2-152016-00037 to P. Horvath), H2020 (ERAPERMED-COMPASS, DiscovAIR; to P. Horvath), Chan Zuckerberg Initiative (Deep Visual Proteomics; to P. Horvath), and Berta Kamprad Foundation, Lund, Sweden (FBKS-2020-22 to M. Rezeli; FBKS-2020-18 to Y. Kim and J. Gil).

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.

Note: Supplementary data for this article are available at Cancer Research Online (http://cancerres.aacrjournals.org/).

1.
Mizrahi
JD
,
Surana
R
,
Valle
JW
,
Shroff
RT
.
Pancreatic cancer
.
Lancet
2020
;
395
:
2008
20
.
2.
Hidalgo
M
.
Pancreatic cancer
.
N Engl J Med
2010
;
362
:
1605
17
.
3.
Neoptolemos
JP
,
Kleeff
J
,
Michl
P
,
Costello
E
,
Greenhalf
W
,
Palmer
DH
.
Therapeutic developments in pancreatic cancer: current and future perspectives
.
Nat Rev Gastroenterol Hepatol
2018
;
15
:
333
48
.
4.
Zhang
J
,
Bajari
R
,
Andric
D
,
Gerthoffert
F
,
Lepsa
A
,
Nahal-Bose
H
, et al
.
The International Cancer Genome Consortium data portal
.
Nat Biotechnol
2019
;
37
:
367
9
.
5.
Heath
AP
,
Ferretti
V
,
Agrawal
S
,
An
M
,
Angelakos
JC
,
Arya
R
, et al
.
The NCI genomic data commons
.
Nat Genet
2021
;
53
:
257
62
.
6.
Rebours
V
,
Gaujoux
S
,
D'Assignies
G
,
Sauvanet
A
,
Ruszniewski
P
,
Lévy
P
, et al
.
Obesity and fatty pancreatic infiltration are risk factors for pancreatic precancerous lesions (PanIN)
.
Clin Cancer Res
2015
;
21
:
3522
8
.
7.
Li
J
,
Gu
D
,
Lee
SSY
,
Song
B
,
Bandyopadhyay
S
,
Chen
S
, et al
.
Abrogating cholesterol esterification suppresses growth and metastasis of pancreatic cancer
.
Oncogene
2016
;
35
:
6378
88
.
8.
Li
J
,
Qu
X
,
Tian
J
,
Zhang
JT
,
Cheng
JX
.
Cholesterol esterification inhibition and gemcitabine synergistically suppress pancreatic ductal adenocarcinoma proliferation
.
PLoS One
2018
;
13
:
1
11
.
9.
Guillaumond
F
,
Bidaut
G
,
Ouaissi
M
,
Servais
S
,
Gouirand
V
,
Olivares
O
, et al
.
Cholesterol uptake disruption, in association with chemotherapy, is a promising combined metabolic therapy for pancreatic adenocarcinoma
.
Proc Natl Acad Sci U S A
2015
;
112
:
2473
8
.
10.
Nicolle
R
,
Blum
Y
,
Marisa
L
,
Loncle
C
,
Gayet
O
,
Moutardier
V
, et al
.
Pancreatic adenocarcinoma therapeutic targets revealed by tumor–stroma cross-talk analyses in patient-derived xenografts
.
Cell Rep
2017
;
21
:
2458
70
.
11.
Erkan
M
,
Hausmann
S
,
Michalski
CW
,
Fingerle
AA
,
Dobritz
M
,
Kleeff
J
, et al
.
The role of stroma in pancreatic cancer: diagnostic and therapeutic implications
.
Nat Rev Gastroenterol Hepatol
2012
;
9
:
454
67
.
12.
Özdemir
BC
,
Pentcheva-Hoang
T
,
Carstens
JL
,
Zheng
X
,
Wu
CC
,
Simpson
TR
, et al
.
Depletion of carcinoma-associated fibroblasts and fibrosis induces immunosuppression and accelerates pancreas cancer with reduced survival
.
Cancer Cell
2014
;
25
:
719
34
.
13.
Storz
P
,
Crawford
HC
.
Carcinogenesis of pancreatic ductal adenocarcinoma
.
Gastroenterology
2020
;
158
:
2072
81
.
14.
Maurer
HC
,
Olive
KP
.
Laser capture microdissection on frozen sections for extraction of high-quality nucleic acids
.
Methods Mol Biol
2019
;
1882
:
253
9
.
15.
Mund
A
,
Coscia
F
,
Kriston
A
,
Hollandi
R
,
Kovács
F
,
Brunner
A
, et al
.
Deep visual proteomics defines single-cell identity and heterogeneity
.
Nat Biotechnol
2022
;
40
:
1231
40
.
16.
Brasko
C
,
Smith
K
,
Molnar
C
,
Farago
N
,
Hegedus
L
,
Balind
A
, et al
.
Intelligent image-based in situ single-cell isolation
.
Nat Commun
2018
;
9
:
1
7
.
17.
Kuras
M
,
Woldmar
N
,
Kim
Y
,
Hefner
M
,
Malm
J
,
Moldvay
J
, et al
.
Proteomic workflows for high-quality quantitative proteome and post-translational modification analysis of clinically relevant samples from formalin-fixed paraffin-embedded archives
.
J Proteome Res
2021
;
20
:
1027
39
.
18.
Bateman
A
.
UniProt: a worldwide hub of protein knowledge
.
Nucleic Acids Res
2019
;
47
:
D506
15
.
19.
Leek
JT
,
Johnson
WE
,
Parker
HS
,
Jaffe
AE
,
Storey
JD
.
The SVA package for removing batch effects and other unwanted variation in high-throughput experiments
.
Bioinformatics
2012
;
28
:
882
3
.
20.
J
WE
,
Li
C
,
Rabinovic
A
.
Adjusting batch effects in microarray expression data using empirical Bayes methods
.
Biostatistics
2007
;
8
:
118
27
.
21.
Jassal
B
,
Matthews
L
,
Viteri
G
,
Gong
C
,
Lorente
P
,
Fabregat
A
, et al
.
The reactome pathway knowledgebase
.
Nucleic Acids Res
2020
;
48
:
D498
503
.
22.
Erichson
NB
,
Zheng
P
,
Manohar
K
,
Brunton
SL
,
Kutz
JN
,
Aravkin
AY
.
Sparse principal component analysis via variable projection
.
SIAM J Appl Math
2020
;
80
:
977
1002
.
23.
Kim
MS
,
Pinto
SM
,
Getnet
D
,
Nirujogi
RS
,
Manda
SS
,
Chaerkady
R
, et al
.
A draft map of the human proteome
.
Nature
2014
;
509
:
575
81
.
24.
Wang
D
,
Eraslan
B
,
Wieland
T
,
Hallström
B
,
Hopf
T
,
Zolg
DP
, et al
.
A deep proteome and transcriptome abundance atlas of 29 healthy human tissues
.
Mol Syst Biol
2019
;
15
:
1
16
.
25.
Klami
A
,
Virtanen
S
,
Leppaaho
E
,
Kaski
S
.
Group factor analysis
.
IEEE Trans Neural Networks Learn Syst
2015
;
26
:
2136
47
.
26.
Leppäaho
E
,
Ammad-Ud-Din
M
,
Kaski
S
.
GFA: exploratory analysis of multiple data sources with group factor analysis
.
J Mach Learn Res
2017
;
18
:
1
5
.
27.
Szklarczyk
D
,
Gable
AL
,
Lyon
D
,
Junge
A
,
Wyder
S
,
Huerta-Cepas
J
, et al
.
STRING v11: protein–protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets
.
Nucleic Acids Res
2019
;
47
:
D607
13
.
28.
Le Large
TYS
,
Mantini
G
,
Meijer
LL
,
Pham
TV
,
Funel
N
,
van Grieken
NCT
, et al
.
Microdissected pancreatic cancer proteomes reveal tumor heterogeneity and therapeutic targets
.
JCI Insight
2020
;
5
;
e138290
.
29.
Cao
L
,
Huang
C
,
Zhou
DC
,
Hu
Y
,
Lih
TM
,
Savage
SR
, et al
.
Proteogenomic characterization of pancreatic ductal adenocarcinoma
.
Cell
2021
;
184
:
5031
52
.
30.
Bankhead
P
,
Loughrey
MB
,
Fernández
JA
,
Dombrowski
Y
,
McArt
DG
,
Dunne
PD
, et al
.
QuPath: open source software for digital pathology image analysis
.
Sci Rep
2017
;
7
:
1
7
.
31.
Smith
K
,
Li
Y
,
Piccinini
F
,
Csucs
G
,
Balazs
C
,
Bevilacqua
A
, et al
.
CIDRE: an illumination-correction method for optical microscopy
.
Nat Methods
2015
;
12
:
404
6
.
32.
Schindelin
J
,
Arganda-Carreras
I
,
Frise
E
,
Kaynig
V
,
Longair
M
,
Pietzsch
T
, et al
.
Fiji: an open-source platform for biological-image analysis
.
Nat Methods
2012
;
9
:
676
82
.
33.
Preibisch
S
,
Saalfeld
S
,
Tomancak
P
.
Globally optimal stitching of tiled 3D microscopic image acquisitions
.
Bioinformatics
2009
;
25
:
1463
5
.
34.
Berg
S
,
Kutra
D
,
Kroeger
T
,
Straehle
CN
,
Kausler
BX
,
Haubold
C
, et al
.
Ilastik: interactive machine learning for (bio)image analysis
.
Nat Methods
2019
;
16
:
1226
32
.
35.
Ofverstedt
J
,
Lindblad
J
,
Sladoje
N
.
Fast and robust symmetric image registration based on distances combining intensity and spatial information
.
IEEE Trans Image Process
2019
;
28
:
3584
97
.
36.
Ruifrok
AC
,
Johnston
DA
.
Quantification of histochemical staining by color deconvolution
.
Anal Quant Cytol Histol
2001
;
23
:
291
9
.
37.
Manders
EMM
,
Verbeek
FJ
,
Aten
JA
.
Measurement of co-localization of objects in dual-colour confocal images
.
J Microsc
1993
;
169
:
375
82
.
38.
Costes
SV
,
Daelemans
D
,
Cho
EH
,
Dobbin
Z
,
Pavlakis
G
,
Lockett
S
.
Automatic and quantitative measurement of protein-protein colocalization in live cells
.
Biophys J
2004
;
86
:
3993
4003
.
39.
Pirhonen
J
,
Arola
J
,
Sädevirta
S
,
Luukkonen
P
,
Karppinen
S-M
,
Pihlajaniemi
T
, et al
.
Continuous grading of early fibrosis in NAFLD using label-free imaging: a proof-of-concept study
.
PLoS One
2016
;
11
:
e0147804
.
40.
Yang
C
,
Wu
Q
,
Huang
K
,
Wang
X
,
Yu
T
,
Liao
X
, et al
.
Genome-wide profiling reveals the landscape of prognostic alternative splicing signatures in pancreatic ductal adenocarcinoma
.
Front Oncol
2019
;
9
:
1
15
.
41.
Kawalerski
RR
,
Leach
SD
.,
Escobar-Hoyos
LF
. Pancreatic cancer driver mutations are targetable through distant alternative RNA splicing dependencies
.
Oncotarget
2021
;
12
:
525
33
.
42.
Kushner
IK
,
Clair
G
,
Purvine
SO
,
Lee
JY
,
Adkins
JN
,
Payne
SH
.
Individual variability of protein expression in human tissues
.
J Proteome Res
2018
;
17
:
3914
22
.
43.
Taieb
D
,
Roignot
J
,
André
F
,
Garcia
S
,
Masson
B
,
Pierres
A
, et al
.
ArgBP2-dependent signaling regulates pancreatic cell migration, adhesion, and tumorigenicity
.
Cancer Res
2008
;
68
:
4588
96
.
44.
Giordano
F
,
Saheki
Y
,
Idevall-Hagren
O
,
Colombo
SF
,
Pirruccello
M
,
Milosevic
I
, et al
.
XPI(4,5)P2-dependent and Ca2+-regulated ER-PM interactions mediated by the extended synaptotagmins
.
Cell
2013
;
153
:
1494
.
45.
Neuhöfer
P
,
Roake
CM
,
Kim
SJ
,
Lu
RJ
,
West
RB
,
Charville
GW
, et al
.
Acinar cell clonal expansion in pancreas homeostasis and carcinogenesis
.
Nature
2021
;
597
:
715
9
.
46.
Backx
E
,
Coolens
K
,
Van den Bossche
JL
,
Houbracken
I
,
Espinet
E
,
Rooman
I
.
On the origin of pancreatic cancer: molecular tumor subtypes in perspective of exocrine cell plasticity
.
Cell Mol Gastroenterol Hepatol
2022
;
13
:
1243
53
.
47.
Oni
TE
,
Biffi
G
,
Baker
LA
,
Hao
Y
,
Tonelli
C
,
Somerville
TDD
, et al
.
SOAT1 promotes mevalonate pathway dependency in pancreatic cancer
.
J Exp Med
2020
;
217
:
e20192389
.
48.
Greenlee
JD
,
Subramanian
T
,
Liu
K
,
King
MR
.
Rafting down the metastatic cascade: the role of lipid rafts in cancer metastasis, cell death, and clinical outcomes
.
Cancer Res
2021
;
81
:
815
7
.
49.
Mobin
MB
,
Gerstberger
S
,
Teupser
D
,
Campana
B
,
Charisse
K
,
Heim
MH
, et al
.
The RNA-binding protein vigilin regulates VLDL secretion through modulation of Apob mRNA translation
.
Nat Commun
2016
;
7
;
12848
.
50.
Bai
R
,
Rebelo
A
,
Kleeff
J
,
Sunami
Y
.
Identification of prognostic lipid droplet-associated genes in pancreatic cancer patients via bioinformatics analysis.
Lipids Health Dis
2021
;
20
:
1
12
.
51.
Snaebjornsson
MT
,
Janaki-Raman
S
,
Schulze
A
.
Greasing the wheels of the cancer machine: the role of lipid metabolism in cancer
.
Cell Metab
2020
;
31
:
62
76
.

Supplementary data