This single-center prospective study was conducted at Zhejiang Provincial People’s Hospital. From May 2016 to November 2016, a total of 50 patients primarily diagnosed with PDAC and received surgery were enrolled in this study. The database was locked for follow-up and analyses on November 2018. Patients with a concurrent malignant neoplasm were excluded. The histopathological status was evaluated by at least two experienced pathologists. TNM staging was defined according to American Joint Committee on Cancer (AJCC) TNM staging system for pancreatic cancer (10). Radiographic assessment using the Response Evaluation Criteria in Solid Tumors (RECIST) version 1.1 was performed and based on standard of care clinical guidelines. We followed the Strengthening the Reporting of Observational Studies in Epidemiology reporting guideline statement to ensure the quality of data reported in this study (11). This study was approved by the ethical committee at Zhejiang Provincial People’s Hospital (No. 2016KY129). All patients provided informed written consent before undergoing any study-related procedures. This study was performed in accordance with the Declaration of Helsinki.

Surgical tumor tissue samples and blood lymphocytes were collected from each patient. Tissue samples were fixed by formalin-fixed paraffin-embedded (FFPE). H&E staining was performed for each section according to the manufacturer’s instructions (Beyotime, Suzhou, China) before dissection. As a guide for stromal lesion, we performed immunohistochemistry staining for vimentin (Cat #5741 purchased from Cell Signaling Technology, Danvers, MA, USA) in the adjacent section. Ten to 15 sections of FFPE (thickness: 5–10 μm) were cut using a microtome (RM2265, Leica, Germany). LCM was performed to separate the neoplastic components, fibrotic stroma, and normal pancreatic tissue on a Leica LMD7000 microscope as previously described (12). The microdissection was performed by at least two senior pathologists, and any disagreement between these pathologists was resolved by discussion. All specimens were then stored in phosphate-buffered saline (PBS) before subsequent processing.

Genomic DNA was extracted from neoplastic, stromal, and normal components using a QIAamp DNA Mini Kit (Qiagen, Hilden, Germany). Germline DNA was extracted from blood lymphocytes using the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany). DNA quality was estimated using a Qubit fluorometer and a Qubit dsDNA (BR) Assay Kit (Invitrogen, Carlsbad, CA, USA). Genomic DNA from normal tissues and germline DNA from blood lymphocytes were used as negative control to eliminate the interference of germline variants and contaminating cancer cells.

Extracted DNA was sheared into 200- to 250-bp fragments via a Covaris S2 instrument (Woburn, MA, USA). KAPA LTP Library Preparation Kit for Illumina (KAPA Biosystems, Boston, MA, USA) was used to prepare indexed NGS libraries. NEBNext FFPE DNA Repair Mix (Ipswich, UK) was used for FFPE DNA repair during library construction, and the detailed protocol can be obtained from https://international.neb.com/protocols/2015/01/16/protocol-for-use-with-nebnext-ffpe-dna-repair-mix-m6630-and-other-user-supplied-library-construction-reagents. Additional information regarding library preparation was described by Lv et al. (13).

Libraries were hybridized to custom-designed biotinylated oligonucleotide probes (Integrated DNA Technologies, IA, USA). The captured genomic regions included the most common driver genes of solid tumors (14). We chose their entire exome regions to construct the basic panel. Next, genomic regions related relevant to the effects of chemotherapy, targeted drugs, and immunotherapy per available clinical and preclinical research were added to the panel. Finally, high-frequently mutant regions recorded in the Catalogue of Somatic Mutations in Cancer (COSMIC, http://cancer.sanger.ac.uk/cosmic) and The Cancer Genome Atlas (TCGA; https://cancergenome.nih.gov/) were involved. Overall, 1,021 genes were involved in this panel. Sequencing was carried out using Illumina 2 × 100-bp paired-end reads on an Illumina HiSeq 3000 instrument according to the manufacturer’s recommendations using a TruSeq PE Cluster Generation Kit v3 and a TruSeq SBS Kit v3 (Illumina, San Diego, CA, USA). Additional detailed information regarding library preparation was described by Lv et al. (13). The median sequencing depth of stromal and neoplastic components was 941× (360× to 1,626×) and 1,045× (345× to 1728×), respectively.

Terminal adaptor sequences and reads with more than 50% low-quality base reads, or those with more than 50% N bases, together with their mate pair were removed from raw reads. Subsequently, Burrows–Wheeler Aligner (BWA; version 0.7.12-r1039, http://bio-bwa.sourceforge.net/) tool was used to align clean reads to the reference human genome (hg19) with default parameters. Duplicate reads were identified and marked with Picard’s Mark Duplicates tool (https://software.broadinstitute.org/gatk/documentation/tooldocs/4.0.3.0/picard_sam_markduplicates_MarkDuplicates.php). The Gene Analysis Toolkit (GATK, https://www.broadinstitute.org/gatk/) was used to perform local realignment and base quality recalibration.

Somatic single-nucleotide variations (SNVs) and insertions or deletions of small fragments (indels) were called using the MuTect2 algorithm (https://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/org_broadinstitute_gatk_tools_walkers_cancer_m2_MuTect2.php). The filter criteria included 1) variants supported by fewer than five high-quality reads (base quality ≥30, mapping quality ≥30) were filtered; 2) variants were filtered as cross-contamination if present in >1% samples in custom single-nucleotide polymorphism (SNP) databases (dbsnp, https://www.ncbi.nlm.nih.gov/projects/SNP/; 1000G, https://www.1000genomes.org/; ESP6500, https://evs.gs.washington.edu/; ExAC, http://exac.broadinstitute.org/) and self-built SNP database; 3) synonymous mutations (also listed in Table S1 ) were removed; 4) variants with allele frequency less than 1% were removed; and 5) variants detected in matched blood lymphocytes and normal tissue were removed. The final candidate variants were all manually verified in the Integrative Genomics Viewer (IGV; https://igv.org/). Remaining mutations were considered validated somatic variants.

Two steps were performed to determine driver or passenger mutations. First, evidential driver genes of PDAC were determined according to Bailey et al. (15). Second, Polymorphism Phenotyping v2 (PolyPhen-2, http://genetics.bwh.harvard.edu/pph2/) and Sortig Intolerant From Tolerant (SIFT, http://sift.bii.a-star.edu.sg/) were used to predict whether the protein structural change derived by one mutation was harmful or not. Those with PolyPhen-2 Score >0.85 or SIFT Score <0.05 were defined as harmful mutations. Generally, harmful mutations in driver genes were defined as driver events, and the others, including harmless mutations in driver genes and all mutations in passenger genes, were defined as passenger events.

PyClone algorithm was used to determine the clonal clusters (16). The essential parameters included the variant allele frequencies (VAFs) and copy numbers of non-synonymous mutations in both tumor and stromal components. Copy number was estimated by Contra algorithm (http://contra-cnv.sourceforge.net).

Human pancreatic cancer cell line Panc1 was purchased from American Type Culture Collection (Manassas, VA, USA), and it was grown in DMEM culture medium (Gibco, Carlsbad, CA, USA) supplemented with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin, in a humidified atmosphere at 37°C and 5% CO₂.

Total RNA was extracted from cells using TRIzol reagent (Invitrogen, California, USA) according to the manufacturer’s instruction. cDNA was subsequently synthesized using PrimeScript™ RT Master Mix kit (Takara). RT-PCR was performed using SYBR Premix Ex Taq™ kit (Takara, Dalian, China). The primers are listed in Table S2 .

Panc1 cells were seeded on coverslips and cultured under different glycemic conditions for 3 days. Then medium was removed, and vimentin was stained using the rabbit anti-human vimentin antibody (Cat #5741) purchased from Cell Signaling Technology (Danvers, MA, USA). Evaluation was performed using confocal laser scanning microscope (Leica, Wetzlar, Germany).

Descriptive statistics was performed using SPSS 22.0 (IBM, Armonk, NY, USA). Spearman’s correlation analysis was performed using GraphPad Prism 7 (GraphPad Software, La Jolla, CA, USA) to assess the relevance of VAFs between mutants KRAS and TP53, as well as the correlation between tumor cell fraction (TCF) and maximal VAF in stromal or neoplastic component. The parameter comparison for different patient groups was performed using the Mann–Whitney U-test (two groups) or one-way ANOVA test (≥3 groups) (SPSS 22.0). The Kaplan–Meier survival analysis was used to compare disease-free survival (DFS) between different subgroups, and Cox regression was performed to determine the influence of multiple factors on DFS. Both analyses were conducted by SPSS 22.0. A two-tailed p-value <0.05 was considered statistically significant.

Clinical characteristics of patients are summarized in Table 1 . All of 50 patients were diagnosed with primary PDAC. The median age at diagnosis was 65 years (ranged from 37 to 84 years). The number of male and female patients was 30 (60.0%) and 20 (40.0%), respectively. The majority of enrolled patients were stage I/II (n = 44, 88.0%). Histologically, 25 patient specimens (50.0%) were poorly differentiated in terms of cellular morphology, and the others were moderately differentiated. The maximal diameter of tumor in situ was >4 cm in seven patients (14.0%). Regional lymph nodes were involved in 19 patients (38.0%). The adjacent nerve and vasculature were invaded in 41 (82.0%) and 21 (42.0%) of patients, respectively.

Surgical resection was performed in all patients as the sole treatment for 20 patients (40.0%). Thirty (60.0%) received adjuvant chemotherapy after surgery. Eight (16.0%) patients had ≥2 lines of chemotherapy. At the time of last follow-up, 19 patients (38.0%) experienced local (two, 10.5%) or distant recurrences (17, 89.5%) postoperatively.

Cellular morphology was determined via H&E staining, and the incisal margin of LCM was kept away from nests as much as possible to avoid the contamination of neoplastic cells. As a specific biomarker of stroma-derived cells, vimentin was enriched in stromal components isolated from the adjacent section ( Figure 1A ). We evaluated the mutant prevalence in the stromal components. In total, 127 somatic mutations (median = 3, ranged from 1 to 7) were detected in stromal components from 39 patients ( Table S1 ). KRAS (n = 28, 71.8%), TP53 (n = 24, 61.5%), and CDKN2A (n = 9, 23.1%) were the most recurrent mutant genes ( Figure 1B ). All of KRAS mutations were located in the hot spot codons 12 and 61, including G12C (n = 1), G12D (n = 11), G12R (n = 2), G12V (n = 12), Q61H (n = 2), and Q61R (n = 1) ( Figure S1 ). Two KRAS mutations (G12R and G12V) co-existed in P05. Twenty-one TP53 mutations (87.5%) occurred in DNA-binding domain ( Figure S1 ). As the top two prevalent mutant genes, KRAS and TP53 were co-altered in 21 patients ( Figure 1B ).

Droplet digital PCR (ddPCR) was used to verify the credibility of stromal mutations, and the most common mutation site (KRAS G12D/V) from seven stromal specimens with different VAF range was selected. As the results, all of seven mutations were indeed repeated via ddPCR, and the VAFs estimated by ddPCR and NGS demonstrated strong consistency (R² = 0.9067, p = 0.0009, Figure S2 ).

Correlation between clinical characteristics and genomic alterations of the stroma was then explored. The stroma was more likely to be genomically altered in patients without vascular invasion than in those with vascular invasion (p = 0.0464, Table 1 ). Mutants KRAS (p = 0.0036) and TP53 (p = 0.0101), as well as co-mutants KRAS and TP53 (p = 0.0029), more likely occurred in older patients than in younger ones, although the distribution of age at diagnosis suggests no discrepancy between patients with and without stromal mutations ( Table 1 ). Nevertheless, there was no significant correlation between other baseline characteristics and mutations in stroma ( Table 1 ).

Next, a total of 248 mutations were detected from 50 neoplastic specimens (median = 5, ranged from 1 to 11, Table S1 ). Commonly altered genes in neoplastic components included KRAS (n = 47, 94.0%), TP53 (n = 43, 86.0%), CDKN2A (n = 12, 24.0%), and SMAD4 (n = 9, 18.0%) ( Figure S3 ). We validated the mutation prevalence of both stromal and neoplastic components using TCGA data. Altered KRAS and TP53 were the most commonly altered genes and demonstrated concordance in three cohorts, with less prevalence in the stroma than in the other two cohorts (71.2% and 61.5% in stroma, 94.0% and 86.0% in tumor, and 90.7% and 69.3% in TCGA) ( Figure 1C ).

Stromal mutations were further verified in matched neoplastic components. Overall, 114 mutations co-existed in both components. Thirteen stromal mutations were absent in matched neoplastic components, while 134 mutations were exclusive to the neoplastic components ( Figures 1D , S3 ). Critical driver gene analysis and function prediction showed 71 candidate driver events (62.3%, 71/114) that occurred in nine genes, including KRAS, TP53, CDKN2A, ARID1A, GNAS, KDM6A, RNF43, SMAD4, and TGFBR2 ( Figure S4 ). All mutations of KRAS and TP53 in stroma were also identified in matched neoplasm, except one KRAS mutation in P05, which harbored two different KRAS mutations in the stroma ( Figure S3 ). Furthermore, KRAS mutations identified in the neoplastic components were absent in matched stroma for 19 patients (38.0%). The same TP53 mutations were identified in matched neoplastic components of only 19 patients (38.0%) ( Figure S3 ). For 13 stroma-specific mutations, there were only two driver events (15.4%, KRAS G12V for P05, SMAD4 W302* for P43, Table S1 ), while 57 driver events (42.5%) were underlined in 134 neoplasm-specific mutations. We further performed pathway enrichment analysis for stroma-specific and neoplasm-specific mutations using Kyoto Encyclopedia of Genes and Genomes (KEGG) resource. As a result, a more intimate connection with tumorigenesis and development was highlighted for neoplasm-specific mutations compared with stroma-specific mutations ( Figure S5 ). Those co-existed mutations in both components indicated that cancer cells affected the genomic features of stroma via EMT. On the other hand, those non-tumor-related mutations in stromal components demonstrated that genomic variants indeed occurred in TME, which might be attributed to clonal evolution.

The TCF could impact the analysis of clonality. The correlation between TCF estimated by microimaging and the mutational status of the neoplasm and stroma were evaluated ( Figure 2A ). The TCF in patients with mutant stroma was significantly higher than in those with normal stroma (p = 0.0019, Figure 2B ). However, the VAFs in neither stroma ( Figure S6A ) nor neoplasm ( Figure S6B ) were correlated with the TCF. Patients with KRAS/TP53-mut stroma also demonstrated a higher TCF than did those with wild-type stroma (p = 0.0371, p = 0.0014, Figures 2C, D ). Nevertheless, no discrepancy was seen for patients with or without KRAS/TP53-mut neoplasm ( Figure S7 ), possibly due to the extremely high detection rate of KRAS and TP53 mutation in neoplastic components.

To eliminate the frequency bias between neoplastic and stromal mutations, we normalized the VAFs of each mutation (absolute VAF/maximum VAF in the same specimen) in paired stromal and neoplastic components. Mutations with ≥50.0% normalized VAFs were more likely to be clonal events and occurred early during development of stroma or neoplasm than those with <50.0% normalized VAFs. For 114 mutations common in both components, 99 mutations (86.8%) were clonal events in stroma, among which 96 (84.2%) presented ≥50.0% normalized VAFs in both neoplasm and stroma ( Figure 3A ). However, for 134 neoplasm-specific mutations, only 72 (53.7%) were clonal in neoplastic components ( Figure 3A ), significantly lower than the proportion of common mutations (χ² p < 0.0001).

In neoplastic and stromal components, 95.8% (46/48) and 89.7% (26/29) of KRAS mutations were clonal, respectively ( Figure 3A ). For mutant TP53, 97.7% (42/43) in neoplasm and 87.5% (21/24) in stroma were clonal variants ( Figure 3A ). The VAFs of KRAS and TP53 mutations in the same stroma presented statistically significant consistency (Spearman’s r² = 0.7481, p < 0.0001, Figure 3B ), indicating that mutants KRAS and TP53 might have co-occurred around the same early period. Moreover, stroma with co-occurred KRAS and TP53 tended to harbor more mutations in stroma than the others (p < 0.0001) ( Figure 3C ). Similarly in neoplastic components, KRAS and TP53 often co-altered and showed significant consistency (Spearman’s r² = 0.6841, p < 0.0001, Figure 3D ). More mutations were detected in neoplasm with co-mutations than those without co-mutations (p = 0.0272, Figure 3E ). Altogether, these results indicated that the neoplasm and adjacent stroma were likely to share the similar clonal trunk events, especially for mutants KRAS and TP53, during the early tumorigenesis.

Subsequently, PyClone strategy was utilized to reconstruct the evolutionary trajectory for neoplastic components. Ultimately, four subtypes were clarified for patients according to the genetic and driver imprinting derived from neoplasm upon stroma. The first type (A) included 12 patients characterized by the overall genomic and evolutional concordance, indicating that all neoplastic clones and driver events could be identified in stroma ( Figures 4A , S8A ). The second type (B) included 20 patients, and all stroma harbored the neoplastic initial clones and driver events but lacked some of subsequent clones ( Figures 4B , S8B ). The third type (C) involving five patients demonstrated total absence of neoplastic mutations in stroma ( Figures 4C , S8C ). The fourth type (D) included only two patients. For the first patient, the initial clone was absent, but a driver mutation (RNF43 K568Sfs*132) involved in a latter clone was identified in stroma ( Figure 4D ). For the second patient, the initial neoplastic clone indeed expressed in stroma, but no driver mutation was identified ( Figure S8D ). Of note, patients with mutants KRAS and TP53 in stroma were either type A or B. There was no significant discrepancy of TCF among the four different subtypes ( Figure S9 ). Overall, type A/B was defined as neoplasm-like stroma based on the similar genetic performance and evolutional entanglements between neoplastic and stromal components.

Epithelial cancer cell were observed to migrate into the stroma and transformed to mesenchymal phenotype via EMT in a PanIN mouse model carried KRAS mutation (9). Based on this study and aforementioned findings, multi-lineage differentiation potential and mesenchymal transformation of KRAS-mutant cell line Panc1 were evaluated to validate this result. After 3 days of cell culture in DMEM with 5 and 25 µM of dextrose, five regulator markers of cell pluripotency involving OCT4, BMI1, NANOG, SOX2, and CD24 were detected by RT-PCR. These five markers exhibited significant increase in the 25 µM group compared with the 5 µM group ( Figure 5A ). Furthermore, the expression of E-cadherin was significantly decreased while β-catenin was significantly increased in a high-nutrient environment, indicating the progress of EMT and cell migration ( Figure 5A ). Subsequently, immunofluorescent staining revealed an enhanced expression of vimentin in the 25 µM dextrose group but absent in the 5 µM dextrose group ( Figure 5B ). Those findings indicated that a fraction of tumor cells harbor multi-lineage differentiation potential and that KRAS mutation might facilitate this process.

The association between DFS and genomic status of stroma, as well as multiple clinicopathologic risk factors of PDAC, was analyzed for 32 patients who were followed up over 1 year after surgical operation. The Kaplan–Meier analysis showed that 32 patients with neoplasm-like stroma had a markedly reduced DFS time (median = 3.9 months) than had the other patients (median DFS unreached, hazard ratio = 3.079, 95% CI 1.126 to 7.215, p = 0.0193, Figure 6A ). We further evaluated whether stromal KRAS mutations were associated with postoperative survival. As the results, patients with KRAS-mutant stroma had a significantly poorer DFS time (median = 3.9 months) than those with KRAS-wild-type stroma (median DFS unreached, hazard ratio = 3.304, 95% CI 1.247 to 8.751, p = 0.0162, Figure 6B ), which was also confirmed by multivariate analysis (hazard ratio = 2.962, 95% CI 1.174 to 7.471, p = 0.021, Table 2 ). Patients with TP53-mutant stroma also showed poorer DFS time (median = 3.8 months) than those with TP53-wild-type stroma (median DFS unreached, hazard ratio = 3.143, 95% CI 1.112 to 8.880, p = 0.0307). However, although exhibiting a certain trend in univariate analysis, all clinicopathologic risk factors of PDAC, including clinical stage, tumor size, histological differentiation, lymph node involvement, nerve, and vasculature invasion, showed no significant association with DFS in univariate and multivariate analyses ( Figure S10 , Table 2 ).