Quantitative proteomic analysis of pancreatic cyst fluid proteins associated with malignancy in intraductal papillary mucinous neoplasms.

BACKGROUND: The application of advanced imaging technologies for identifying pancreatic cysts has become widespread. However, accurately differentiating between low-grade dysplasia (LGD), high-grade dysplasia (HGD), and invasive intraductal papillary mucinous neoplasms (IPMNs) remains a diagnostic challenge with current biomarkers, necessitating the development of novel biomarkers that can distinguish IPMN malignancy.
METHODS: Cyst fluid samples were collected from nine IPMN patients (3 LGD, 3 HGD, and 3 invasive IPMN) during their pancreatectomies. An integrated proteomics approach that combines filter-aided sample preparation, stage tip-based high-pH fractionation, and high-resolution MS was applied to acquire in-depth proteomic data of pancreatic cyst fluid and discover marker candidates for IPMN malignancy. Biological processes of differentially expressed proteins that are related to pancreatic cysts and aggressive malignancy were analyzed using bioinformatics tools such as gene ontology analysis and Ingenuity pathway analysis. In order to confirm the validity of the marker candidates, 19 cyst fluid samples were analyzed by western blot.
RESULTS: A dataset of 2992 proteins was constructed from pancreatic cyst fluid samples. A subsequent analysis found 2963 identified proteins in individual samples, 2837 of which were quantifiable. Differentially expressed proteins between histological grades of IPMN were associated with pancreatic diseases and malignancy according to ingenuity pathway analysis. Eighteen biomarker candidates that were differentially expressed across IPMN histological grades were discovered-7 DEPs that were upregulated and 11 that were downregulated in more malignant grades. HOOK1 and PTPN6 were validated by western blot in an independent cohort, the results of which were consistent with our proteomic data.
CONCLUSIONS: This study demonstrates that novel biomarker candidates for IPMN malignancy can be discovered through proteomic analysis of pancreatic cyst fluid.

Background

Intraductal papillary mucinous neoplasms (IPMNs) are precancerous lesions that grow in the pancreatic ducts and are characterized by papillary growth of the ductal epithelium. The production of thick mucinous fluid, another hallmark of IPMNs, causes cystic dilation and can progress into pancreatic ductal adenocarcinoma [1-4]. Depending on the malignancy, IPMN is classified as low-grade dysplasia (LGD), intermediate-grade dysplasia (IGD), high-grade dysplasia (HGD), and invasive IPMN. According to the official guidelines for managing pancreatic IPMN, only patients with HGD or invasive IPMN require surgery, because they are at higher risk of their disease developing into cancer [5]. Milder forms of IPMN can be managed with active surveillance and do not warrant surgical intervention. However, current methods for assessing the histological grades of IPMNs are unreliable, and as a result, patients with milder IPMN are often subjected to unnecessary operations [6-10].

In clinical practice, MRI and CT scans, cytological examination of cyst fluid, measurement of tumor markers such as carcinoembryonic antigen (CEA) and carbohydrate antigen 19-9 (CA 19-9), and analysis of GTPase Kras (KRAS) and guanine nucleotide-binding protein alpha subunit (GNAS) mutations are used to categorize patients with pancreatic cysts [6, 7, 10–16]. Features of pancreatic images in MRI or CT scans are generally used to assess the potential malignancy of cysts but have low diagnostic accuracy—up to 40% of neoplastic cysts are misdiagnosed as pseudocysts, and the overall accuracy ranges from 20 to 80% [17-19]. Cytological examination of pancreatic cyst fluid is an alternative approach, but it has difficulties in identifying the existence of malignancy when sufficient sample volumes are unavailable [16, 20–23]. Differentiating mucinous cysts from other cystic lesions by measuring carcinoembryonic antigen levels in cyst fluid has relatively low accuracy (79% sensitivity, 73% specificity) [17, 24]. Similarly, as shown by Frossard et al. [25], CA 19-9, a pancreatic cancer marker, also performs poorly in distinguishing mucinous cysts and other lesions, with 15% sensitivity and 81% specificity [16]. Analyzing GNAS mutations are only applicable for samples that are acquired during the early stages of IPMN [20, 23, 26, 27]. The general consensus is that existing methods for diagnosing IPMN histological grades are imprecise and unreliable, even when used in tandem [6, 7, 10, 16, 17, 20].

Because pancreatic cyst fluid contains secreted proteins from tumor cells at higher proportions, several groups, such as Poersch et al. [28], have concluded that it is a better experimental model of IPMN histological grades than serum and plasma [16, 29–32]. Consequently, pancreatic cyst fluid has been widely favored in recent research on IPMN, because it is obtainable by endoscopic ultrasound-guided fine needle aspiration biopsy, which is minimally invasive [6, 25, 33]. Many studies have focused on discovering protein markers that differentiate mucinous from nonmucinous cyst fluid and cyst fluid that is related to IPMN dysplasia, based on DNA methylation and telomerase activity, as demonstrated by Hata et al. [20]. Diagnosing histological grades of IPMN using pancreatic cyst fluid by proteomic analysis is a relatively unexplored area [6, 7, 20, 34, 35]. Thus, the IPMN dysplasia proteome has not been characterized extensively.

Cuoghi et al. [36] performed a cursory profiling study of the proteomic patterns of pancreatic cyst fluids from various cystic lesions, including IPMN, MCN, serous cystadenomas, pancreatic neuroendocrine tumors, and pseudocysts, identifying 220–727 proteins in these fluids. Specifically, 243 proteins were identified in the IPMN groups. Gbormittah et al. [37] characterized glycoproteins and nonglycoproteins in mucinous and nonmucinous pancreatic cyst fluid to identify DEPs as potential biomarker targets. They found 230 proteins in mucinous subtypes and 290 proteins in nonmucinous subtypes; the DEPs between mucinous and nonmucinous cyst fluid were associated with lipid metabolism, energy metabolism, and stress responses. These studies were unable to determine the IPMN histological grades, merely differentiating between mucinous and nonmucinous cyst fluid. These recent studies demonstrate that the current cyst fluid proteome lacks the coverage to extrapolate meaningful conclusions on the molecular and biological activities of the identified proteins, which ultimately impedes our understanding of IPMN histology in terms of proteomic differences and biological functions.

In this report, we aimed to comprehensively identify pancreatic cyst fluid proteins and discover differentially expressed proteins in accordance with histological grades of IPMN. Recently, we reported a platform for in-depth profiling of pancreatic cyst fluid [38]. Using this platform, the protein expression patterns of pancreatic cyst fluid were analyzed on a high-resolution mass spectrometer to discover potential biomarkers of IPMN histological grades. Subsequently, we validated some of the 18 candidate markers by western blot. We report here that pancreatic cyst fluid is a valuable source for biomarker studies as it contains putative markers related to IPMNs and that bioinformatics analyses using identified proteins of cyst fluid enhance our understanding of IPMNs at the molecular level. We ultimately intend to discover marker candidates that can help patients avoid unnecessary operations.

Methods

Clinical samples

Cyst fluid samples were collected from 9 IPMN patients during their pancreatectomies at Seoul National University Hospital (Seoul, South Korea) from April 2013 to December 2015. At least 200 μL of cyst fluid was aspirated from each patient. The samples were then snap-frozen in liquid nitrogen and stored at − 80 °C. All patients consented to participation in the study in accordance with Institutional Review Board guidelines (IRB No. 1301-095-458). IPMN samples were divided into low-grade dysplasia (LGD, n = 3), high-grade dysplasia (HGD, n = 3), and invasive IPMN (n = 3).

Pancreatic cyst fluid protein sample preparation

Each pancreatic cyst fluid sample was transferred to an Eppendorf tube. Viscous samples that could not be pipetted were sonicated briefly (Sonics and Materials Inc., USA) to remove the mucus. All samples were centrifuged at 15,000 rpm for 20 min at 4 °C, and the supernatant was placed into a new tube. The protein concentration was estimated using a BCA reducing agent compatibility assay kit (Thermo Scientific, Rockford, IL, USA). Equal portions of each sample were pooled to create a peptide library from 600 µg of proteins. One hundred micrograms of individual protein samples were used for label-free quantification. Cold acetone (Sigma-Aldrich, USA) was added to the supernatant to the ratio of 5:1 (v/v) to precipitate the proteins. The mixture was vortexed gently and incubated overnight at − 20 °C. The precipitate was centrifuged for 10 min (15,000 rpm at 4 °C), and the supernatant was carefully decanted, after which 500 µL cold acetone was added to the pellet. After this wash step, the pellet was centrifuged for 10 min (15,000 rpm at 4 °C). The remaining acetone was poured off, and the pellet was air-dried for 2 h.

Protein digestion and desalting

The pellet was dissolved in 30 μL of lysis buffer (4% SDS, 0.1 M DTT, 0.1 M Tris–Cl, pH 7.4). The mixture was gently vortexed and boiled for 30 min at 95 °C. The boiled mixture was then transferred through a 30-kDa cutoff filter (Amicon® Ultra, Millipore, USA) with 300 μL 8 M urea (8 M Urea, 0.1 M Tris–Cl, pH 8.5) and centrifuged (14,000g, 15 min, 20 °C). This filtration step was repeated twice to dilute and lower the SDS concentration. Next, 200 μL 50 mM IAA (50 mM IAA, 8 M urea, 0.1 M Tris–Cl, pH 8.5) was added to each sample and incubated for 1 h at 25 °C. Each sample was then centrifuged and washed twice with 300 μL 8 M urea and then three times with 300 μL 40 mM ammonium bicarbonate (ABC).

After the samples were centrifuged, 100 μL 40 mM ABC and 0.1 μg/μL trypsin (at a trypsin:sample ratio of 1:80, wt/wt) were added to each sample and incubated for 18 h at 37 °C. Next, the filters (nine individual samples, one pooled sample) were transferred to new collection tubes, which were centrifuged after 100 μL 40 mM ABC was added. Fifty microliters NaCl was added to each individual sample, and 50 μL water was added to the pooled sample. The pooled sample underwent an additional digestion step [39, 40]. Again, the filter unit was transferred to a new tube and centrifuged after 200 μL 8 M urea was added. Then, the unit was centrifuged twice with 300 μL 40 mM ABC. One-tenth of the concentration of trypsin that was used in the first digestion step was added with 100 μL 40 mM ABC, and the unit was incubated for 18 h at 37 °C. Next, the filter was transferred to another tube, and the peptides were collected by sequential centrifugation with 100 μL 40 mM ABC and 50 μL 0.5 M NaCl.

Prior to acidification and desalting, all tryptic peptides were measured by tryptophan fluorescence assay to determine the volume that was required to extract the same amount of peptides from each sample [41]. The equalized amounts of peptides were then set aside for label-free quantification. The measured peptides were acidified with 10 μL 10% TFA and desalted with homemade C18-StageTip columns as described [42]. The desalted peptides were then lyophilized on a speed-vacuum centrifuge and stored at − 80 °C.

Peptide fractionation by high-pH reverse phase fractionation

To increase the number of identified proteins, the pooled cyst fluid sample was fractionated using two methods: modified stage-tip-based high-pH peptide fractionation [43, 44] and offline HPLC high-pH fractionation on an Agilent 1260 Bio-inert. For stage-tip fractionation, half of the lyophilized peptides were dissolved in 200 μL of loading buffer (15 mM ammonium hydroxide solution, pH 10, and 2% acetonitrile) and separated on a pipette-based C18 RP microcolumn. The column was constructed by plugging the bottom of a 200 μL transparent pipette tip with C18 Empore disk membrane (3 M, Bracknell, UK) and packing the tip with POROS 20 R2 resin. The plugged tip was rinsed three times with 100 μL 100% methanol and then three times with 100 μL 100% acetonitrile (ACN). The column was then conditioned with 100 μL of loading buffer using a syringe. The peptides were loaded onto the column at pH 10. An ACN gradient of 2, 5, 7.5, 10, 12.5, 15, 17.5, 20, 22.5, 25, 27.5, 30, 32.5, 35, 40, 50, 60, 70, 80, and 100% was used to elute 20 fractions, which were collected into six tubes discontinuously to distribute eluents of varying hydrophobicity. These six fractions were lyophilized in a speed-vacuum centrifuge and stored at − 80 °C.

The remaining half of the lyophilized peptides was dissolved in 80 μL of loading buffer (15 mM ammonium hydroxide in water, pH 10). The peptides were loaded onto the column, and 96 (2 mL Square Collection Plate, Waters, UK) fractions were eluted by applying an ACN gradient (pH 10, 5–35%) for 40 min at a flow rate of 0.2 mL/min and washing the column with 90% ACN for 10 min at 0.2 mL/min. The ACN gradient was established by mixing varying proportions of solution A (0.1% formic acid in HPLC-grade distilled water) and solution B (0.1% formic acid in ACN). The 96 fractions were concatenated according to the column number of the plate to produce 12 pooled fractions. The resulting 12 tubes were lyophilized in a speed-vacuum centrifuge and stored at − 80 °C.

LC–MS/MS analysis

The peptide samples were analyzed using an LC–MS/MS configuration, comprising an Easy-nLC 1000 (Thermo Fisher Scientific, Waltham, MA, USA) that was coupled to a Q Exactive mass spectrometer with a nanoelectrospray ion source (Thermo Fisher Scientific, Waltham, MA, USA), per our established protocol [38, 44, 45]. Peptides were separated on a 2-column system that was composed of a trap column (75 μm I.D. × 2 cm, C18 3.0 μm, 100 Å) and an analytical column (50 μm I.D. × 15 cm, C18 3.0 μm, 100 Å).

Fractionated peptides were subjected to an ACN gradient (6–60%) for 235 min. The gradient was created by mixing solvent A (2% ACN and 0.1% v/v formic acid) and solvent B (100% acetonitrile and 0.1% v/v formic acid) at various proportions. The spray voltage was set to 2.0 kV in positive ion mode, and the temperature of the heated capillary was set to 320 °C. Mass spectra were acquired in data-dependent mode by top 20 method on an Orbitrap analyzer with a mass range of 350–1700 m/z and a resolution of 70,000 at m/z 200. HCD scans were acquired at a resolution of 17,500. HCD peptide fragments were acquired at a normalized collision energy (NCE) of 27. The maximum ion injection time for the survey scan and MS/MS scan was 20 and 80 ms, respectively. All samples were analyzed in three technical replicates.

Raw data search

The MS data from the Q Exactive were processed in MaxQuant (version 1.5.5.1 with built-in Andromeda search engine) [46]. Precursor MS signal intensities were determined, and HCD MS/MS spectra were de-isotoped and filtered, such that only the 20 most abundant fragments per 100 m/z range were retained. Protein groups were identified by searching the MS and MS/MS data of the peptides against the Uniprot human database (2014 December, 88,717 entries). Both the forward and reverse amino acid sequences were taken into account when calculating the false discovery rate (FDR). Following established target-decoy search procedures [47], search results were filtered at FDR < 1% for identifying peptides, modification sites, and proteins. The search was conducted in digestion mode trypsin/P, which assumes cleavage at carboxyl sides of lysine and arginine, including cases where the subsequent residue is a proline.

The following parameters were used in the database search: precursor and HCD fragment mass tolerances of 6 and 20 ppm, respectively; tolerance of up to two missed cleavages; carbamidomethylation of cysteine as a fixed modification; and oxidation of Met and acetylation of protein N-term as variable modifications. The minimum peptide length was set to six residues. Peptides were assigned to protein groups by the principle of parsimony [48-50]. The principle is applied to derive the smallest list of probable protein groups that adequately represent the identified peptides, which reduces sequence redundancy issues. All proteomics data in this report have been deposited in the ProteomeXchange Consortium (http://proteomecentral.proteomexchange.org/) through the PRIDE partner repository: dataset identifier PXD008302 [51, 52].

Label-free quantification and statistical analysis

Label-free quantification (LFQ) and statistical analysis were performed in MaxQuant (version 1.5.5.1) and Perseus (version 1.5.8.5), respectively, according to our previous studies [43, 45]. Protein abundance was obtained from LFQ intensity values. LFQ intensity was calculated as described by the equation by Cox et al. [53]. Each of the three histological groups in this study had three biological replicates, which in turn had three technical replicates each. Thus, a total of 9 LFQ intensity values exist per histological group (three biological replicates × three technical replicates). LFQ intensity values greater than zero were deemed valid. Proteins with at least six valid values within a histological group were used in statistical analysis for label-free quantification. This criterion was used to reduce the possibility of analyzing proteins that are nonspecific to histological grades. After log2-transformation of protein intensities, the missing values were replaced with expected intensities based on the normal distribution (imputation width = 0.3, shift = 1.8) of log2-transformed LFQ intensities [43]. Student’s t test was applied to the preprocessed dataset of matched proteins to detect DEPs across grades of IPMN dysplasia. The comparative pairs for the statistical analysis were LGD versus HGD (comparison 1), HGD versus invasive IPMN (comparison 2), and LGD versus invasive IPMN (comparison 3). A Benjamini–Hochberg FDR threshold of 0.05 was applied to each pair to find significantly changed proteins. Subsequently, the expression patterns of overlapping DEPs across two or more pairs were analyzed to screen for biomarker candidates. DEPs that had expression patterns that varied based on the malignancy of IPMN were selected as final biomarker candidates. The resulting DEPs were subjected to hierarchical clustering in Perseus (version 1.5.8.5) with the following parameters: Euclidean distance, average linkage, the number of clusters of 100, maximal number of iterations of 10, the number of restarts of 1, and k-means preprocessing prior to clustering.

Bioinformatics analysis

The gene ontologies (GOs) of all DEPs were annotated using the DAVID bioinformatics resource tool (https://david.ncifcrf.gov/) and the UniprotKB database (http://www.uniprot.org/). The GO analysis included information on biological process (BP), cellular component (CC), and molecular function (MF). Pathway analysis was performed using the KEGG database (http://www.genome.jp/kegg/). Secretory protein prediction and functional annotation were performed using SignalP 4.1 (http://www.cbs.dtu.dk/services/SignalP/), SecretomeP 2.0 (http://www.cbs.dtu.dk/services/SecretomeP/), and TMHMM, server 2.0 (http://www.cbs.dtu.dk/services/TMHMM/). Ingenuity pathway analysis (IPA) was used to conduct functional analysis (Ingenuity Systems, http://www.ingenuity.com/). The plasma proteome database (PPD) was used to confirm the association between the proteins that were identified in human plasma and the proteins that were identified in this study [54, 55]. The proteins that were identified in our dataset were crossreferenced with mRNA and protein expression in pancreatic sections in the Human Protein Atlas (http://www.proteinatlas.org/).

Western blot analysis

A total of 19 pancreatic cyst fluid samples—10 LGD, 4 HGD, and five invasive IPMN—were used to validate the candidate markers. Equal volumes of a pooled cyst fluid sample were loaded onto each gel to correct for the intensity of the blots. Pancreatic cyst fluid samples were mixed with 5× SDS loading dye (250 mM Tris–Cl, pH 6.8, 10% SDS, 50% glycerol, 0.5 M DTT, 0.1% bromophenol blue). Proteins (40 μg, as measured by BCA assay) were separated on 10% SDS-PAGE gels and transferred to polyvinylidene fluoride (PVDF) membranes (Hybond-P, GE Healthcare, Pittsburgh, PA, USA). The membranes were stained with Ponceau S dye (P7170, Sigma-Aldrich, USA), blocked with 5% BSA for 2 h at RT, and incubated overnight at 4 °C with the following primary antibodies: rabbit monoclonal anti-HOOK1 (ab150397, Abcam, Cambridge, UK) at 1:250, mouse monoclonal anti-PTPN6 (sc-7289, Santa Cruz Biotechnology, USA) at 1:1000, and mouse polyclonal anti-SERPINA5 (ab67368, Abcam, Cambridge, UK) at 1:100. The membranes were then washed five times with Tris-buffered saline and Tween-20 (TBS-T) before being incubated with the following HRP-conjugated secondary antibodies: anti-rabbit (ab6721, Abcam, Cambridge, UK) at 1:1000 and anti-mouse (ab6789, Abcam, Cambridge, UK) at 1:2500 for 2 h at RT. The membranes were developed with ECL solution (West-Q chemiluminescent substrate Kit-plus, GenDEPOT, Barker, TX, USA), and the signals were visualized on an LAS-4000 (Fujifilm, Japan).

Results

Clinical sample characterization

Pancreatic cyst fluid samples from nine patients were classified into three groups: LGD (n = 3), HGD (n = 3), and invasive IPMN (n = 3). The samples did not differ significantly in composition, with the exception of serum CEA level and CA 19-9 concentration measured by chemiluminescent microparticle immunoassay and cyst size (Table 1). The invasive IPMN patient group had the highest average CEA and CA19-9 concentrations at 7.67 ± 7.06 and 117.17 ± 142.78 mg/L, respectively. CEA and CA19-9 levels were generally higher in the more severe forms of IPMN. The average CEA concentration was approximately 3-fold higher for HGD than LGD subjects and 7-fold higher in invasive IPMN versus LGD. In addition, the average CA19-9 level was approximately 2-fold and 30-fold greater for these comparisons. Our samples were consistent with several publications that have reported that malignant cysts tend to be larger, as evidenced by our invasive IPMN samples (6.63 ± 3.74 cm) being twice as large as LGD (2.93 ± 0.54 cm) and HGD (2.50 ± 0.41 cm) samples on average [5, 56–59].

Table 1

Demographic and clinical characteristics of the study population used in the label-free quantification

Group	Pancreatic cyst fluids
	LGD(n = 3)	HGD(n = 3)	Invasive IPMN(n = 3)
Age (years)	69.00 ± 1.41	66.33 ± 8.58	58.33 ± 11.09
Gender
Male	1	2	1
Female	2	1	2
Gland type
Gastric	2	2	1
Intestinal	0	1	1
Oncocytic	0	0	1
Unknown	1	0	0
Duct type
Main	0	0	1
Branch	2	1	0
Mixed	0	2	2
Unknown	1	0	0
Cyst focality
Single	2	3	3
Multiple	1	0	0
Unknown	0	0	0
Mural nodule
Y	0	3	3
N	3	0	0
Unknown	0	0	0
Cyst location
Head	1	0	2
Body/tail	1	3	1
Mixed	1	0	0
CEA concentration (mg/L)	1.13 ± 0.53	3.07 ± 1.27	7.67 ± 7.06
CA 19-9 concentration (mg/L)	4.00 ± 1.48	6.87 ± 7.17	117.17 ± 142.78
Cyst size
Cyst size (cm)	2.93 ± 0.54	2.50 ± 0.41	6.63 ± 3.74
< 3.0	1	2	1
≥ 3.0	2	1	2

In-depth analysis of pancreatic cyst fluid

The overall scheme of the study was based on a proteomic platform of cyst fluids that we established earlier [38]. In this study, nine individual pancreatic cyst fluid samples of various types [LGD (n = 3), HGD (n = 3), and invasive IPMN (n = 3)] were used for label-free quantification. All samples were centrifuged, and only the supernatant was used. The same portions of individual samples were pooled and fractionated to generate a peptide library, which increased the depth of the identified proteins. In contrast, the nine individual samples were not fractionated. After a series of sample preparation steps, LC–MS/MS analysis was performed on a Q Exactive mass spectrometer. The MS data were processed in MaxQuant (version 1.5.5.1), and the statistical analysis was performed in Perseus (version 1.5.8.5) (Fig. 1; Additional file 1: Fig. S1a).

Fig. 1

Detailed experimental workflow. Pancreatic cyst fluid samples from nine individuals (3 LGD, 3 HGD, and 3 invasive IPMN) were included in this study. After centrifugation, pellets and debris were discarded, and the supernatant was collected for proteomic analysis. Proteins in each individual sample and pooled cyst fluid (comprised of equal portions of individual samples) were precipitated with cold acetone. Following FASP digestion, the pooled cyst fluids were fractionated using two types of high-pH fractionation methods. Every prepared sample was analyzed on a Q Exactive mass spectrometer

In total, 2992 proteins were identified and 2938 proteins were quantified (Additional file 2: Table S1). A total of 28,728 peptides were identified, and 553,199 peptide spectra matches were found. In the peptide library and the nine individual samples, 2778 and 2963 proteins were identified, respectively. Comparing the peptide library with the individual samples, 2749 proteins (91.9% of all identified proteins) were shared (Fig. 2a). In the nine individual samples, most of the identified proteins (95.7%) were usable for quantitative analysis, as evidenced by the 2963 and 2837 proteins that were identified and quantified (Fig. 2b). Approximately 2200–2500 proteins were quantified in each sample group. The three IPMN groups were similar with regard to the number of quantified proteins (Additional file 1: Fig. S1b). In contrast, there was substantial individual variation in the number of identified and quantified proteins within the same histological subgroups. This pattern was observed across all nine samples (Additional file 1: Fig. S1c).

Fig. 2

Overlap of identified proteins in the individual samples and peptide library and comparison of identified and quantified proteins. a All identified proteins in the nine individual samples and peptide library; 91.9% of proteins were identified both in the individual samples and peptide library. b All identified proteins and quantified proteins in the nine individual samples; 95.7% of quantifiable proteins overlapped with the identified proteins

On average, the number of identified and quantified proteins increased by 129 and 83, respectively, in individual samples when matched with the peptide library. In addition, the number of identified peptides rose by 752 on average in individual samples with HGD 1 displaying the greatest improvement of 2109. (Additional file 3: Table S2). As shown by the Venn diagram, approximately 77% of identified and 63% of quantified proteins overlapped in all histological groups and 337 additional proteins were identified exclusively when the search was performed with the generated peptide library. Whereas the number of proteins that overlapped in the three histological groups decreased by approximately 6% when searched without the peptide library (Additional file 1: Fig. S1d, e, Additional file 2: Table S1). This result implies that the number of proteins that were common between individual samples rose due to the contribution of the peptide library. As shown in Fig. 3, the percentage of overlapping proteins from the biological replicates in each histological grade ranged from 27 to 46%. The dynamic range of the proteome spanned over seven orders of magnitude overall, but most proteins (95%) were expressed within four orders (Fig. 4). Overall, the proteins with lower orders of magnitude were analyzed, and tumor marker proteins, such as MUC5AC, MUC1, and CEA, were quantified with high intensity in the dataset.

Fig. 3

Identified and quantified proteins in three individual samples for each histological grade. Identified proteins of three biological replicates in LGD (a), HGD (b), and invasive IPMN (c). Quantified proteins of three biological replicates in LGD (d), HGD (e), and invasive IPMN (f) (INV: invasive IPMN)

Fig. 4

Dynamic range of quantified proteins. Distribution of expression intensities of quantified proteins show a large dynamic range of abundance, but 95% of the proteins were expressed within four orders of magnitude. Several tumor marker proteins, such as MUC2, CEA, and KRAS, were quantified

Comparative analysis using other proteome databases

Our bioinformatics analysis showed that secreted proteins accounted for 60.5% (1810 proteins) of the 2992 identified proteins (Additional file 1: Fig. S2a, Additional file 2: Table S1). Across SecretomeP, SignalP, and TMHMM, 1527, 682, and 381 proteins were identified, respectively (Additional file 1: Fig. S2b, Additional file 2: Table S1). Protein accession numbers were converted into gene names, and the resulting redundancy was discarded prior to comparative analysis. We compared our dataset with the Human Plasma Proteome Database to assess the likelihood that the discovered proteins are potential blood markers [54, 55]. As a result, 2299 (79.7%) of the identified proteins were plasma proteins (Additional file 1: Fig. S2c, Additional file 2: Table S1). To determine whether the discovered proteins are expressed in the pancreas, the dataset was crossreferenced with The Human Protein Atlas (http://www.proteinatlas.org, May 31, 2017)—2613 genes had corresponding mRNA entries and 2021 genes had corresponding protein entries in the pancreas (Additional file 1: Fig. S2e, Additional file 2: Table S1).

Variation in individual cyst samples

Coefficient of variation (CV%) values were calculated to evaluate the reproducibility of the technical and biological replicates. The CV% values of log2-transformed LFQ intensity sums of technical triplicates of individual samples ranged from 0.32 to 6.45% (Additional file 4: Table S3). All CV% values of log2-transformed LFQ intensities of each quantified protein in technical triplicates of individual samples were less than 20%. Moreover, the average CV% value of individual samples ranged from 1.085 to 1.524% (Additional file 1: Fig. S3a). Pearson correlation coefficients of the LFQ intensities of technical triplicates and their averages were greater than 0.9 (Fig. 5a–c, Additional file 1: Fig. S3b–d). These data suggest that the variation between technical replicates was low. In contrast, the variation between biological triplicates was generally high, based on the Pearson correlation coefficients, which ranged from 0.370 (between LGD1 and LGD2) to 0.789 (between HGD1 and HGD2) (Fig. 5d–f, Additional file 1: Fig. S3b–d).

Fig. 5

Pearson correlation coefficients of technical replicates (TRs) and biological replicates in each histological group. a–c Pearson correlation coefficients of technical replicates in each histological group. Three types of marks indicate the Pearson correlation coefficients of each comparison between technical replicates (◆TR1 vs. TR2, ■TR1 vs. TR3, ▲TR2 vs. TR3). Red dots represent the average Pearson correlation coefficient in each of the 3 comparisons. a–c indicate LGD, HGD, and invasive IPMN, respectively. Box plot representation of Pearson correlation coefficients between each biological replicate in LGD (d), HGD (e), and invasive IPMN (f)

Differentially expressed proteins in IPMN dysplasia

The 1751 proteins that had at least six valid values within a histological group were used for the statistical analysis (Additional file 5: Table S4). By Student’s t test (Benjamini–Hochberg FDR = 0.05), 149, 48, and 98 proteins were differentially expressed between comparisons 1 (LGD versus HGD), 2 (HGD versus invasive IPMN), and 3 (LGD versus invasive IPMN), respectively (Additional file 5: Table S4, Fig. 6), 75, 32, and 64 of which were upregulated. By unsupervised hierarchical clustering, the DEPs clustered by IPMN histology (Additional file 1: Fig. S4).

Fig. 6

Volcano plots based on p values in all comparison groups. To determine the differentially expressed proteins, Student’s t test was performed with a Benjamini–Hochberg FDR value of 0.05. The colored dots indicate the proteins that passed the t test for significance between LGD versus HGD (a), HGD versus invasive IPMN (b), and LGD versus invasive IPMN (c). The blue dots represent downregulated proteins, and the red dots denote upregulated proteins

There were 243 DEPs across comparisons 1, 2, and 3. Among the 243 DEPs (Fig. 7), 142 were upregulated and 91 were downregulated in at least 1 comparison group (Additional file 6: Table S5). Enriched DEPs were used to conduct GO and KEGG pathway analyses to identify their overrepresented terms in biological process, molecular function, and cellular component. The DEPs from comparisons 1 and 3 were analyzed by ingenuity pathway analysis (IPA) bioinformatics tool to track biological processes that became activated or more pronounced in aggressive malignancy.

Fig. 7

Venn diagram of differentially expressed proteins in three comparison groups. By Student’s t test (Benjamini–Hochberg FDR = 0.05), 149, 48, and 98 proteins were differentially expressed between LGD and HGD, between HGD and invasive IPMN, and between LGD and invasive IPMN, respectively. A total 243 proteins were found to be DEPs when overlapping components of the Venn diagram were excluded; 49 proteins were shared in at least two comparison groups and are highlighted in white

By GO enrichment analysis, 243 DEPs were involved primarily in vesicle-mediated transport and cell surface receptor signaling with regard to biological process. The analysis also found that 76.6% of DEPs were extracellular region proteins. The molecular functions of the DEPs were primarily associated with peptidase activity and regulation (Additional file 1: Fig. S5a–c, Additional file 7: Table S6). By KEGG pathway analysis, the 142 upregulated proteins were associated with the ribosome, oxidative phosphorylation, and endocytosis, whereas the 91 downregulated proteins were linked to leukocyte transendothelial migration, focal adhesion, and ECM-receptor interaction (Additional file 1: Fig. S5d).

The significantly changed proteins from comparison 1 and 3 were examined by IPA with regard to biological processes that are related to pancreatic cysts and aggressive malignancy. Core analysis in IPA was used to evaluate the biological functions that were most likely to be affected by changes in expression of proteins in our dataset. As a result, 149 DEPs in comparison 1 and 98 DEPs in comparison 3 were associated with such terms as cellular movement and angiogenesis in Diseases and Functions, which are indicative of malignancy; the biological terms that correlated with aggressive malignancy are highlighted in yellow (Additional file 1: Fig. S6a, b). Cell growth and vasculogenesis were upregulated among the DEPs in comparison 1. A total of 98 DEPs in comparison 3 were upregulated in most Diseases and Functions terms, except for apoptosis of tumor cell lines—particularly metastasis-related terms, such as cell spreading and angiogenesis.

Comparison analysis is usually performed to visualize trends in protein expression across several analyses. Consistent with our expectations, the Diseases and Bio functions heat map of the comparison analysis demonstrated that the DEPs that were associated with cell movement of endothelial cells and angiogenesis were more highly expressed in comparison 3 versus 1. The term “apoptosis of tumor cell lines” was downregulated in comparison 3 but unchanged in comparison 1 (Additional file 1: Fig. S6c). A higher percentage of DEPs in comparison 3 was associated with pancreas-specific diseases, such as chronic pancreatitis and associated diseases than DEPs in comparison 1 (Additional file 1: Fig. S6d).

Biomarker candidates for IPMN malignancy

Proteins that were shared by at least two comparison groups were chosen as the initial marker candidates. Our rationale was that significantly changed proteins that are common to several comparison groups are more likely to be associated with the malignancy of IPMNs [60]. A total of 49 candidates expressed in at least two comparison groups were selected from 243 DEPs (Fig. 7). Then, the DEPs that had expression patterns that were consistent with the degree of IPMN malignancy were selected as the final candidates. Table 2 details the results of the statistical analysis of the 49 DEPs, including the p value, fold-change, and t test significance for each comparison group. Of the 49 DEPs, 38 were secreted proteins and 33 were confirmed to be expressed in the pancreas as mRNA or proteins in The Human Protein Atlas. In addition, 35 proteins were confirmed to be expressed in plasma, according to the Human Plasma Proteome Database (Table 2).

Table 2

Detailed statistical analysis of 49 DEPs that were shared in at least two comparison groups

Protein name	LGD versus HGD			HGD versus invasive IPMN			LGD versus invasive IPMN
	t test significance	p value	Fold-change	t test significance	p value	Fold-change	t test significance	p value	Fold-change
EH domain-containing protein 1	+	0.00594	− 15.43	+	0.01955	7.35		0.54376	− 2.10
Immunoglobulin lambda variable 3–21	+	0.00586	− 18.35	+	0.00001	41.24		0.54057	2.25
Collagen alpha-1(VI) chain	+	0.00224	− 22.78	+	0.01403	5.61		0.19350	− 4.06
Vacuolar protein sorting-associated protein 37B	+	0.01350	− 6.78	+	0.04669	4.56		0.76562	− 1.49
Ig lambda-6 chain C region	+	0.03265	− 4.96	+	0.02993	20.47		0.31759	4.13
Ras-related protein Ral-B	+	0.00571	− 5.53	+	0.00111	12.71		0.33130	2.30
Lysozyme C	+	0.03497	3.95	+	0.03528	− 5.79		0.66323	− 1.47
Immunoglobulin heavy variable 1/OR15-1		0.00000	1.00	+	0.02490	5.06	+	0.00046	7.47
Immunoglobulin lambda variable 5–45		0.93977	− 1.15	+	0.03513	17.43	+	0.02740	15.11
bAldo–keto reductase family 1 member B10		0.97902	− 1.08	+	0.00808	− 131.01	+	0.02709	− 140.91
Ig heavy chain V–III region TIL		0.22132	− 1.77	+	0.00288	6.15	+	0.00022	3.48
Ig heavy chain V–III region BUT		0.61619	1.58	+	0.01276	4.09	+	0.03193	6.47
Ig heavy chain V–III region BUR		0.61583	1.72	+	0.03530	5.43	+	0.03177	9.36
Ig heavy chain V–III region ZAP		0.81620	1.29	+	0.01443	6.03	+	0.00587	7.80
Ig alpha-2 chain C region		0.14392	− 2.13	+	0.02303	3.70	+	0.02778	1.74
aThymidine phosphorylase		0.89045	− 1.25	+	0.00774	48.91	+	0.01333	39.00
aTestis-expressed sequence 12 protein		0.75598	− 1.40	+	0.01204	16.69	+	0.00379	11.94
aProtein Hook homolog 1		0.33603	2.29	+	0.00053	50.05	+	0.00000	114.50
Collagen alpha-1(XV) chain	+	0.00001	− 11.37		0.54289	1.75	+	0.01252	− 6.49
Prosaposin	+	0.01353	− 2.27		0.69798	1.16	+	0.00131	− 1.95
bPhosphoinositide-3-kinase-interacting protein 1	+	0.00707	− 7.41		0.00000	1.00	+	0.01999	− 7.35
bThy-1 membrane glycoprotein	+	0.00966	− 8.39		0.00000	1.00	+	0.01176	− 10.48
Latent-transforming growth factor beta-binding protein 2	+	0.01295	− 30.47		0.83621	1.13	+	0.02886	− 26.88
Protein disulfide-isomerase A3	+	0.00534	27.04		0.64689	− 1.78	+	0.00176	15.23
Agrin	+	0.01010	− 95.04		0.64723	1.58	+	0.03007	− 59.97
Alpha-1-antichymotrypsin	+	0.00041	− 8.93		0.52129	1.65	+	0.01226	− 5.40
Collagen alpha-1(IV) chain	+	0.00000	− 23.08		0.36550	2.54	+	0.01235	− 9.09
bTrefoil factor 1	+	0.03148	− 4.27		0.32370	− 4.75	+	0.03214	− 20.28
bPlasma serine protease inhibitor	+	0.00000	− 3.00		0.07044	− 2.33	+	0.00013	− 7.00
Guanine nucleotide-binding protein G(k) subunit alpha	+	0.00446	− 7.09		0.57151	1.44	+	0.01032	− 4.92
Guanine nucleotide-binding protein G(k) subunit alpha	+	0.00042	− 12.58		0.94988	− 1.07	+	0.00077	− 13.45
Collagen alpha-3(VI) chain	+	0.00016	− 41.81		0.61089	2.11	+	0.00316	− 19.85
aTyrosine-protein phosphatase non-receptor type 6	+	0.02363	4.39		0.50095	1.58	+	0.00391	6.95
bKallistatin	+	0.01495	− 6.30		0.68656	1.42	+	0.01178	− 4.45
bFibrillin-1	+	0.00000	− 57.02		0.43700	1.93	+	0.00012	− 29.49
bClaudin-18	+	0.00001	− 14.94		0.67515	− 1.49	+	0.00101	− 22.20
Coronin-7	+	0.00858	6.09		0.85877	− 1.18	+	0.01207	5.15
aNuclease-sensitive element-binding protein 1	+	0.04238	6.17		0.62145	1.67	+	0.00138	10.30
Histone H3.1	+	0.00023	10.62		0.25248	− 2.30	+	0.01205	4.62
bMucin-5AC	+	0.02741	− 15.38		0.69722	− 1.94	+	0.02010	− 29.87
aMucin-2	+	0.00006	7068.07		0.76813	− 3.54	+	0.03333	1996.84
bWAP four-disulfide core domain protein 2	+	0.00554	− 13.23		0.83040	− 1.38	+	0.00690	− 18.23
bCystatin-M	+	0.01152	− 19.49		0.00000	1.00	+	0.03235	− 13.72
BPI fold-containing family B member 1	+	0.00026	− 27.59		0.49944	1.78	+	0.00180	− 15.53
Intelectin-1	+	0.01505	844.39		0.76240	− 3.48	+	0.01605	242.93
aTalin-1	+	0.03653	8.08		0.60098	1.46	+	0.02980	11.78
Basement membrane-specific heparan sulfate proteoglycan core protein	+	0.00003	− 114.37	+	0.00582	13.69	+	0.03552	− 8.35
Histone H3	+	0.00007	22.86	+	0.01650	− 3.15	+	0.00843	7.26
Cadherin-related family member 2	+	0.00024	− 52.14	+	0.02477	13.15	+	0.02747	− 3.96

aDominantly expressed in invasive IPMN

bDominantly expressed in LGD

Of the 49 shared DEPs between groups, 18 had expression patterns that were consistent with the degree of malignancy. PTPN6, MUC2, TLN1, and YBX1 were expressed in lower amounts in LGD but gradually elevated in HGD and invasive IPMN. Conversely, SERPINA5, AKR1B10, and TFF1 expression decreased as IPMN histological grade progressed. Other proteins, such as HOOK1, TYMP, TEX12, FBN1, CLDN18, THY1, MUC5AC, CST6, WFDC2, PIK3IP1, and SERPINA4, were predominantly expressed in LGD or invasive IPMN but not in other groups (Fig. 8, Additional file 1: Fig. S7). Based on these results, these 18 proteins were selected as potential biomarkers of IPMN dysplasia.

Fig. 8

Six biomarker candidates among 18 proteins that had expression patterns that were consistent with the degree of IPMN malignancy. HOOK1 (a), PTPN6 (b), and MUC2 (c) were predominantly expressed in invasive IPMN. FBN1 (d), CLDN18 (e), and SERPINA5 (f) were primarily expressed in LGD (*

Six biomarker candidates among 18 proteins that had expression patterns that were consistent with the degree of IPMN malignancy. HOOK1 (a), PTPN6 (b), and MUC2 (c) were predominantly expressed in invasive IPMN. FBN1 (d), CLDN18 (e), and SERPINA5 (f) were primarily expressed in LGD (*