The degree of metastatic disease varies widely among patients with cancer and affects clinical outcomes. However, the biological and functional differences that drive the extent of metastasis are poorly understood. We analyzed primary tumors and paired metastases using a multifluorescent lineage-labeled mouse model of pancreatic ductal adenocarcinoma (PDAC)—a tumor type in which most patients present with metastases. Genomic and transcriptomic analysis revealed an association between metastatic burden and gene amplification or transcriptional upregulation of MYC and its downstream targets. Functional experiments showed that MYC promotes metastasis by recruiting tumor-associated macrophages, leading to greater bloodstream intravasation. Consistent with these findings, metastatic progression in human PDAC was associated with activation of MYC signaling pathways and enrichment for MYC amplifications specifically in metastatic patients. Collectively, these results implicate MYC activity as a major determinant of metastatic burden in advanced PDAC.
Here, we investigate metastatic variation seen clinically in patients with PDAC and murine PDAC tumors and identify MYC as a major driver of this heterogeneity.
Tumor heterogeneity, most commonly studied in a primary disease setting, is a critical driver of phenotypic diversity, culminating in metastatic, lethal cancers (1–5). In most cancers, prognosis and therapeutic decisions are defined by the presence or absence of metastasis. However, tumor heterogeneity is increasingly being questioned at the level of metastatic disease, with recent studies in several cancer types suggesting that metastasis is not a binary phenotype but rather a disease spectrum ranging from oligometastatic (limited) to polymetastatic (widespread) disesase (6–8). Heterogeneity in the manifestation of metastatic disease can guide decisions on the use of local–regional versus systemic therapies with emerging evidence of its importance in clinical outcome (9–11). Despite its clinical significance, the mechanisms that underlie this spectrum of metastatic states remain unclear and largely understudied.
Pancreatic ductal adenocarcinoma (PDAC) represents a disease entity well suited for the study of metastasis, as most PDACs present with metastatic disease that is associated with a dismal prognosis (12). Genomic studies have comprehensively cataloged core mutations responsible for primary tumor development in PDAC (e.g., KRAS, TRP53, CDKN2A, and SMAD4), paving the path for genomic investigations of metastatic disease and the identification of metastasis-promoting alterations. Indeed, recent sequencing studies and functional analysis in model systems have associated genomic amplification in mutant KRAS alleles with progression from the nonmetastatic (stage III) to metastatic (stage IV) disease state (13). However, genetic factors mediating metastasic heterogeneity in patients and, importantly, the downstream cellular mechanisms remain largely undefined (14–18). Furthermore, it is unclear whether metastasis-associated alterations perturb the tumor microenvironment, whose influence on metastatic behavior is well documented (19–31). Therefore, understanding the interplay between genetic alterations that influence metastatic behavior and the tumor biology that promotes it—via cell-autonomous and/or non–cell-autonomous mechanisms—is crucial for understanding metastasis as a distinct disease state and critical for the development of more effective treatments.
One barrier to understanding metastatic heterogeneity has been a paucity of model systems that capture this natural variation and allow for direct assessments of paired primary tumors and metastases in vivo. This has limited the ability to define factors intrinsic to primary tumors that influence the extent of metastatic spread. We previously developed an autochthonous model of PDAC—the KPCX model—that employs multiplexed fluorescence-based labeling to track the simultaneous development of multiple primary tumor cell lineages and follow them as they metastasize (32). Importantly, this technique facilitates confirmation of lineage relationships in vivo, such that primary tumor clones with substantial metastatic potential can be distinguished from those having poor metastatic potential. Here, we show that this system recapitulates the variation in metastatic burden found in human PDAC, and we use it to dissect molecular and cellular features contributing to metastatic heterogeneity.
Metastatic Burden Is Variable in Human and Murine PDAC
Although most patients with PDAC have metastases (principally liver and lung), the number of metastases is highly variable from patient to patient (33, 34). Importantly, data regarding metastases have largely been obtained at autopsy and thus are confounded by varying treatment histories and reseeding due to end-stage disease (4). Thus, we first sought to characterize the burden of metastases in treatment-naive patients. To this end, we performed a retrospective analysis of initial computed tomography (CT) scans from 55 patients newly diagnosed with metastatic (stage IV) PDAC (Fig. 1A). The total number of lesions in the lung and liver was counted by examining both coronal and sagittal planes for both organs and binned into groups of 10, revealing a wide distribution of metastatic burden (Fig. 1B). K-means clustering identified two metastatic subgroups: a MetLow subgroup (≤10 metastases, 25/55) and a MetHigh subgroup (>10 metastases, 30/55; Fig. 1B; Supplementary Fig. S1A). Primary tumor size, age, sex, and race were not correlated with differences in metastatic burden (Fig. 1C; Supplementary Fig. S1B). However, having a greater number of metastases was associated with worse overall survival (Fig. 1D). Thus, even among patients with stage IV PDAC, metastatic burden is variable and correlates with clinical outcome.
We hypothesized that the differences in metastatic burden seen in human PDAC may also be present in autochthonous murine models. To test this, we used the KPCXY model—in which Cre-mediated recombination triggers expression of mutant KrasG12D and deletion of one allele of Trp53 in the pancreatic epithelium along with YFP and confetti (X) lineage tracers (Fig. 1E; Methods)—to measure metastatic heterogeneity in a cohort of tumor-bearing mice. By exploiting the multicolor features of the KPCX model, we previously showed that these mice harbor (on average) two to five independent primary tumor clones; importantly, the clonal marking of tumors with different fluorophores makes it possible to infer the lineages of primary tumors with different metastatic potential (32). In our earlier work with this model, we noted that in most tumor-bearing animals, even those with multiple primary tumors, liver and lung metastases were driven by a single tumor clone (Fig. 1E; Supplementary Fig. S1C). This suggested that tumor cell–intrinsic factors strongly influence the metastatic behavior of a tumor, even within a single animal.
To quantify differences in metastatic burden, we examined a panel of mice with at least two uniquely labeled fluorescent tumors in which most metastases could be attributed to a specific tumor on the basis of color (Fig. 1E; Supplementary Fig. S1C). A total of 85 primary tumors from 30 mice were examined, and gross metastases to the liver and lung arising from each tumor were then quantified by stereomicroscopy (Methods). Murine PDACs exhibited a wide distribution of metastatic burden, with a pattern resembling that of the human disease (Fig. 1F). Similarly, K-means clustering grouped murine samples into a low-metastasis subgroup (≤10 metastases, 58/85) and a high-metastasis subgroup (>10 metastases, 27/85), which we similarly refer to as MetLow and MetHigh, respectively (Fig. 1F; Supplementary Fig. S1D). As with the human disease, neither primary tumor size nor tumor cell proliferation correlated with metastatic burden (Fig. 1G; Supplementary Fig. S1E). Thus, the KPCXY model recapitulates the intertumoral metastatic heterogeneity seen in human PDAC and provides a unique experimental model for comparing highly metastatic and poorly metastatic tumor clones.
Individual Tumor Lineages in KPCXY Mice Correspond to Clones with Distinct Somatic Copy-Number Profiles
Although primary KPCXY tumors were easily distinguishable based on the expression of a distinct fluorophore, each tumor could have arisen via the clonal expansion of a single cell or through fusion of multiple tumors that happened to share the same color. Somatic copy-number alterations (SCNA) have been shown to provide an unambiguous picture of genomic heterogeneity and lineage relationships between primary tumors and matching metastases in human disease (35). Consequently, we performed copy-number analysis via genome sequencing on a set of 20 primary tumors, including multiregional sampling on a subset of the tumors in which sufficient tissue was available (nine tumors with two to four regions sampled per tumor; Fig. 2A; Supplementary Table S1). Tumors bearing different colors exhibited unique DNA copy-number profiles, indicating that they arose independently (Fig. 2B; Supplementary Fig. S2A; ref. 36). By contrast, multiregional sampling of monochromatic tumors revealed shared copy-number alterations, indicating that all subregions within a given tumor (defined by color) shared a common ancestral lineage (Fig. 2C; Supplementary Fig. S2B). In addition, subregion-specific alterations were also observed, suggesting that subclonal heterogeneity is also present in each tumor (Fig. 2C; Supplementary Fig. S2B). These results suggest that the monochromatic tumors observed in KPCXY mice are clonal in origin and continue to undergo subclonal evolution during tumor progression.
To ascertain the lineage relationships between primary tumors and metastases, we compared DNA copy-number profiles between liver metastases and primary tumors within a given mouse. This revealed that primary tumors and metastases of the same color shared common DNA copy-number profiles across the dataset, confirming on a genetic basis the fluorescence-based lineage relationships (Fig. 2D–F; Supplementary Fig. S2C). As most lung metastases were microscopic and difficult to isolate by dissection, they were not included in the molecular analysis. Together, these results indicate that the lineage history of metastases can be inferred by color and genomic analysis, allowing primary tumors with high versus low metastatic potential to be unambiguously classified.
Genomic and Transcriptional Analyses Identify Myc as a Potential Driver of Metastatic Phenotypes
We next sought to examine the molecular differences that distinguish primary tumors with high versus low metastatic potential. We began by examining large-scale (mega-base level as well as chromosome-wide) SCNAs in 20 MetHigh and MetLow primary tumor samples. This analysis revealed largely similar genome-wide copy-number patterns between MetHigh and MetLow primary tumors, with key PDAC-associated genes, such as loss of heterozygosity at Cdkn2a/b and Trp53 as well as chromosomal gain of Kras occurring at similar frequencies (Supplementary Fig. S3A). Thus, KPCXY tumors exhibit frequent copy-number alterations in canonical PDAC genes, but these alterations do not account for the variation in metastatic behavior between MetHigh and MetLow tumors.
We next asked whether other factors (genomic and/or transcriptional) may be acting to enhance metastasis in the MetHigh group. Focal amplifications in driver oncogenes—Cdk6 and Yap in breast cancer and mutant Kras in PDAC—have been linked to the acquisition of metastatic competence (13, 14, 37, 38). Consistent with prior studies, we observed focal amplicons at genomic regions encoding Cdk6, Yap, and Kras in our tumors (Fig. 3A; Supplementary Fig. S3A; refs. 13, 14, 37–39). However, in contrast to these amplifications, which occurred at equal frequencies in MetHigh and MetLow tumors, focal high-amplitude amplifications in Myc were found in 42.8% (3/7) of MetHigh tumors compared with 7.6% (1/13) of MetLow tumors (Fig. 3B). Thus, Myc amplifications are enriched in MetHigh tumors. In all cases, these amplifications were maintained in paired metastases (Supplementary Fig. S3B). In addition, RNA sequencing (RNA-seq) analysis demonstrated significantly higher levels of Myc transcripts in MetHigh tumors and metastases compared with MetLow tumors (Fig. 3C); overall, Myc was the third-most significantly upregulated gene in MetHigh tumors compared with MetLow tumors (Fig. 3D). Gene set enrichment analysis (GSEA) of the differentially expressed genes between MetHigh and MetLow tumors identified MYC and E2F signatures as highly enriched, along with other signatures that have been implicated in PDAC metastasis, including unfolded protein response, oxidative phosphorylation, and hypoxia (Fig. 3E; Supplementary Table S2; refs. 31, 40, 41). Moreover, MetaCore transcription factor enrichment analysis identified MYC as the transcription factor most significantly associated with genes overexpressed in MetHigh tumors (Fig. 3F), and Ingenuity Pathway Analysis placed Myc at the center of the interactome generated by these differentially expressed genes (Supplementary Fig. S4). Collectively, these results demonstrate a strong association between a tumor's metastatic behavior and the abundance and/or activity of MYC at the genomic and transcriptional levels.
Human PDAC can be grouped into two main transcriptomic subtypes—a well-differentiated classical/exocrine-like/progenitor (classical) subtype and a poorly differentiated squamous/quasi-mesenchymal/basal (basal-like) subtype (15, 16, 18, 42). We found that MetHigh tumors were associated with basal-like PDACs, in line with their more aggressive behavior (Fig. 3G). Likewise, applying murine MetLow and MetHigh signatures (see Methods) to human data from The Cancer Genome Atlas (TCGA) predicted a worse survival—indicative of disease recurrence—for patients with a MetHigh signature (Fig. 3H). These data indicate that murine MetHigh tumors correspond to the more aggressive subtypes of human PDAC.
A Panel of Cell Lines that Preserve the MetLow and MetHigh Phenotypes
To understand the mechanisms underlying these different metastatic properties, we generated a panel of cell lines from six MetHigh tumors and five MetLow tumors. Consistent with the parental in vivo tumors, Myc gene expression and MYC protein levels were higher in the MetHigh lines compared with the MetLow lines (Fig. 4A and B). SCNA analysis in these cells lines found that they retained most of the genomic alterations found in the matched primary samples, including Myc amplifications (Supplementary Fig. S5A and S5B). Furthermore, Myc amplifications were not found in any of the cell lines whose tumors were originally characterized as non–Myc amplified, indicating that in vitro culture does not select for this specific copy-number alteration. Importantly, elevations in MYC mRNA and protein were observed in both the Myc-amplified and nonamplified MetHigh lines, suggesting that elevated MYC expression is a stable phenotype of these cells in culture.
To investigate the metastatic properties of the MetHigh and MetLow lines in vivo, we performed orthotopic implantation of five MetHigh and five MetLow lines into the pancreas of NOD.SCID mice and examined distant organs for evidence of metastasis. Although the weights of MetHigh and MetLow tumors were not significantly different (Supplementary Fig. S6A), MetHigh tumors gave rise to 28-fold more liver and lung metastases compared with MetLow tumors (Fig. 4C). Consistent with the cell line expression differences, the orthotopic MetHigh tumors expressed higher levels of Myc compared with MetLow tumors (Supplementary Fig. S6B). To further confirm that differences in Myc expression were sufficient to drive the metastatic phenotype, we introduced a Myc overexpression (Myc_OE) construct into four MetLow lines and generated orthotopic tumors (Supplementary Fig. S6C). Myc overexpression led to a dramatic (22-fold) increase in liver and lung metastases (Fig. 4D). Thus, cell lines derived from spontaneously generated MetHigh and MetLow tumors retain their metastatic phenotypes upon implantation.
MYC Promotes Tumor Cell Intravasation through the Recruitment of Tumor-Associated Macrophages
To form distant metastases, cancer cells must navigate a series of events collectively referred to as the “metastatic cascade.” These events include (i) intravasation into the bloodstream or lymphatics, (ii) survival in the circulation, (iii) extravasation from the vessel, and (iv) growth and survival at the distant site (43). To determine the step(s) at which MYC was exerting its prometastatic effects, we began by measuring the number of circulating tumor cells (CTC) in orthotopically implanted MetHigh and MetLow tumors and in MetLow tumors engineered to overexpress Myc. Remarkably, CTCs arising from MetHigh and Myc_OE tumors were 38-fold and 17-fold more abundant than those arising from MetLow tumors (Fig. 4E), far greater than the approximately twofold increase in tumor weight resulting from Myc overexpression (Supplementary Fig. S6D). Next, we performed a tail vein metastasis assay, which bypasses the invasion step by introducing tumor cells directly into the bloodstream, and measured lung metastases. Surprisingly, in contrast to the orthotopic tumor experiment, there was no difference in the number of metastases between MetHigh and MetLow lines (Fig. 4F). Moreover, Myc overexpression had no effect on tumor cell survival in the circulation (Supplementary Fig. S6E–S6G). Taken together, these data suggest that MetHigh tumors achieve a higher metastatic rate principally by promoting cancer cell invasion into the circulation, which can be driven by increased Myc expression.
Beyond activation of tumor cell intrinsic programs, MYC can also affect tumor phenotypes by altering the tumor immune microenvironment (TiME; refs. 44–46). Thus, we sought to determine if differences in MYC levels between MetHigh and MetLow tumors were associated with distinct TiMEs. To this end, we examined the immune composition of parental primary tumors by staining for markers of immune cells previously implicated in metastasis of PDAC and other cancers. Although MetHigh and MetLow tumors had a similar degree of neutrophil infiltration, MetHigh tumors had lower numbers of CD3+ T cells but were highly enriched for F4/80+ macrophages (Fig. 5A). Thus, compared with MetLow tumors, the TiME of MetHigh tumors contains an increased number of tumor-associated macrophages (TAM).
To examine the ability of MetHigh and MetLow tumor cells to recruit macrophages, we cocultured these cell lines with primary bone marrow–derived murine macrophages (BMDM) in a transwell migration assay (39). Compared with MetLow cocultures, MetHigh cocultures exhibited greater macrophage migration toward the tumor cells (Fig. 5B). Consistent with these in vitro results, orthotopic tumors generated from MetHigh cell lines exhibited greater TAM infiltration than those generated from MetLow cell lines (Fig. 5C). Furthermore, MetLow lines overexpressing Myc gave rise to tumors with greater TAM infiltration compared with controls (Fig. 5D). These results demonstrate that MetHigh tumors exhibit enhanced macrophage recruitment and implicate Myc expression as a driver of macrophage infiltration. Macrophages have been reported to facilitate metastasis in several cancers (47). This is achieved, in part, through the activation of an “M2-like” polarization state in TAMs characterized by increased expression of Arg1 and CD206 expression (48). To determine whether Myc expression in pancreatic tumor cells alters macrophage phenotypes, we stained for these markers and found that both MetHigh tumors and Myc_OE tumors were enriched for Arg1+ and CD206+ TAMs compared with MetLow and EV control tumors (Fig. 5E–H). Thus, Myc overexpression is associated with an enrichment for M2-like macrophages in the tumor microenvironment.
In breast cancer, macrophages promote tumor cell invasion and metastasis through the development of specialized structures called tumor microenvironment of metastasis doorways, in which macrophages facilitate the movement of cancer cells across an endothelial barrier (49). To investigate whether macrophages might promote metastasis in PDAC through a similar mechanism, we performed an in vitro transendothelial migration (iTEM) assay, in which the proinvasive activity of Myc overexpression and macrophages could be directly assessed (50–52). Tumor cell intravasation across an endothelial monolayer was enhanced by either the addition of macrophages or Myc overexpression, an effect that was greatest when both stimuli were present (Fig. 5I).
To directly test whether TAMs are required for Myc-driven metastasis in PDAC, we performed a macrophage depletion experiment. We generated orthotopic tumors from MetLow Myc_OE cells and 10 days later treated mice with a combination of the colony-stimulating factor receptor inhibitor GW2580 and liposomal clodronate (Fig. 5J). Consistent with prior studies (29, 53–55), this regimen was highly effective at depleting both circulating and tumor-resident macrophages (Supplementary Fig. S7A and S7B) and caused a modest increase in tumor weight (Supplementary Fig. S7C). Depeletion also reduced tissue-resident macrophages in liver and lung premetastatic niches with no significant effect on neutrophil abundance (Supplementary Fig. S7D–S7F). Macrophage depletion resulted in a four- to sixfold reduction in metastases (Fig. 5K), an effect that had no impact on CTC viability, seeding, or outgrowth in the lung (Supplementary Fig. S7G–S7J). Taken together, these data suggest that Myc enhances metastatic spread at least in part by creating a TAM-rich environment in the primary tumor that increases tumor cell invasion. Although previous studies have implicated macrophages in tumor cell infiltration, these results directly link this process to the genomic and transcriptional activation of Myc, which occurs naturally in our model and is subsequently selected for as a driver of metastasis (21, 26, 28, 47, 53, 56–61).
MYC Ehances TAM Recruitment and Metastasis through Increased Expression of CXCL3 and MIF
To identify factors that might be responsible for the increased abundance of macrophages in MetHigh tumors, we mined our RNA-seq data to identify secreted factors that are differentially expressed between MetHigh and MetLow tumors [P < 0.01 and log fold change (logFC) >1]. This resulted in identification of six cytokines/chemokines upregulated in MetHigh tumors (Supplementary Fig. S8A). Interestingly, each of these factors has been previously implicated in regulating macrophage recruitment in pancreatic and other cancer types (59, 62–65). To determine which of these factors may be regulated by MYC in a clinically relevant setting, we examined gene expression data from the COMPASS trial cohort of patients with PDAC, comparing tumors with either high or low levels of MYC expression. This revealed that three of the six factors (MIF, CXCL3, and CCL3) were also enriched in human PDACs exhibiting elevated MYC expression (Fig. 6A). To determine which of these factors depend on MYC for their expression, we used short hairpin RNA to knock down Myc levels in three MetHigh tumor lines (Supplementary Fig. S8B). MIF and CXCL3 expression was significantly reduced following Myc knockdown (Fig. 6B), whereas the expression of the other factors was either unchanged or elevated following Myc knockdown (Supplementary Fig. S8C). Consistent with these findings, Myc overexpression in a MetLow cell line resulted in the upregulation of these factors (Fig. 6C). Furthermore, analysis of the Gene Transcription Regulation Database and Eukaryotic Promoter Database revealed direct evidence for Myc binding to the CXCL3 and MIF promoters (Supplementary Fig. S8D and S8E; refs. 66–69). These results suggest that MYC regulates macrophage recruitment in part through the regulation of Mif and Cxcl3.
To further examine the role of these chemokines in macrophage recruitment, we overexpressed each factor in MetLow lines (Supplementary Fig. S8F) and examined their effects in vivo. Compared with controls, orthotopic tumors from Cxcl3_OE and Mif_OE cells resulted in a three- to fourfold increase in intratumoral macrophages (Fig. 6D) and a significant increase in lung and liver metastases (Fig. 6E). Thus, Cxcl3 and Mif overexpression is sufficient to increase TAM abundance and promote metastasis of MetLow tumors.
Tumor-secreted factors can regulate immune cell phenotypes through specific recptor–ligand interactions and enzymatic activities. The cognate receptor for CXCL3 is CXCR2, which has previously been reported to modulate myeloid cell recruitment in prostate cancer and PDAC (59, 70). Similary, Mif can regulate immune cell function through binding various receptors and its tautomerase activity (71, 72). To examine the role of these factors in macrophage recruitment in our system, we treated Myc_OE tumor cells with a CXCR2 inhibitor (AZD5069) or MIF inhibitor (ISO-1) in an in vitro transwell migration assay (52, 70, 72, 73). Compared with vehicle controls, inhibition of Cxcr2 or Mif led to a 2- to 10-fold decrease in macrophage migration (Fig. 6F). Next, we assessed the impact of these drugs on TAM recruitment and metastasis in vivo (Fig. 6G). Treatment of Myc_OE orthotopic tumors with AZD5069 or ISO-1 reduced TAM recruitment and metastatic burden (Fig. 6H–J) and caused a slight decrease in tumor weight (Supplementary Fig. S8G). Interestingly, compared with each inhibitor alone, the combination of both inhbitors resulted in significantly lower levels of TAMs (F4/80+ and CD206+) and metastasis (Fig. 6H–J). Taken together, these data suggest that multiple Myc-regulated factors contribute to macrophage recruitment and metastasis in PDAC.
Metastasis in Human PDAC Is Associated with MYC Gene Amplification and Elevated Expression
Given the finding that genomic and transcriptional variation in Myc was associated with metastatic heterogeneity in murine PDAC, we sought to determine whether MYC is associated with similar metastatic phenotypes in human PDAC. Because most PDAC samples in the International Cancer Genome Consortium (ICGC) and TCGA are derived from resected stage I/II tumors (18), these datasets provide limited insight into the determinants of metastatic burden. Consequently, we analyzed data from the COMPASS trial cohort (NCT02750657), which is focused on patients with metastatic PDAC and uses laser capture microdissection (LCM) to enrich for tumor cells prior to whole-genome sequencing (WGS) or RNA-seq (14, 74). By comparing primary tumors and metastases, we found that 11.3% (n = 17/133) of metastases were enriched for MYC amplifications compared with 1.61% (n = 4/244) of resectable tumors (Fig. 7A and B; P = 7.6e-5, Fisher test). Likewise, advanced tumors (defined as either locally advanced or metastatic) were significantly enriched for MYC amplifications (9.22%; n = 19/206) compared with resectable tumors having no evidence of metastasis at diagnosis (1.04%; n = 2/192; Supplementary Fig. S9A; P = 1.33e-4). As predicted, amplification was associated with higher levels of MYC mRNA (Supplementary Fig. S9B). MYC-amplified tumors did not exhibit greater genomic instability compared with nonamplified tumors (Supplementary Fig. S9C), indicating that MYC amplification is not a proxy for more generalized chromosome-level events. Moreover, metastases expressed higher levels of MYC mRNA than primary tumors (Fig. 7C; P = 0.00312). These results indicate that MYC amplification and transcriptional upregulation are strongly associated with PDAC metastases.
Next, we examined a separate patient cohort (n = 20) in whom matched primary PDAC tumors and metastases were available for comparison. MYC amplifications were common in the primary tumors of patients with metastatic disease (35.0%; n = 7/20), and these amplifications were retained in the matching metastases (Supplementary Fig. S9D), similar to our mouse model. The enrichment of MYC amplifications in metastatic samples and the observed retention of the amplification when analyzing matched primary/metastasis samples suggests that amplification and/or transcriptional upregulation of MYC in primary PDACs are selected for and retained during tumor metastatic progression. Consistent with this notion, we identified a patient with PDAC in whom single-cell analysis of a paired primary tumor and metastasis revealed enrichment of a MYC-amplified subclone in the metastatic lesion compared with the primary (Fig. 7D and E; Supplementary Fig. S9E). Collectively, these data suggest that enhanced expression and/or genomic amplification of MYC is associated with metastatic spread in human PDAC, complementing our findings from the mouse model.
Phenotypic variation, the result of inter- and intratumoral heterogeneity arising during tumor progression, has made it challenging to understand the molecular mechanisms underlying tumor spread (1, 75). Consequently, the demonstration that certain genes function as “metastasis drivers”—promoting metastasis through mechanisms distinct from their roles in primary tumor growth—has proven elusive (76). In this study, we exploited an autochthonous PDAC model with varying degrees of metastatic spread to explore the molecular basis of naturally arising variation in metastatic burden. This system revealed a strong association between the level of MYC—at either the genomic or transcriptional level—and tumor metastasis, a relationship that was also observed in human PDAC samples. MYC exerts its prometastatic effect at least in part by recruiting proinvasive TAMs, leading to greater tumor cell intravasation into the bloodstream. These activities are not directly related to MYC's well-described role in primary tumor growth (46, 77, 78), as tumors with different levels of Myc expression grow at comparable rates despite dramatic differences in metastatic ability.
Prior work by us and others has examined the genetic events associated with PDAC metastasis. In one study, a comparison of matched primary tumors and metastases from four patients failed to reveal nonsynonymous mutations in driver oncogenes that distinguished primary tumors from metastases (79). By contrast, examination of DNA copy-number changes in mouse and human PDAC revealed a significant association between increased mutant KRAS gene dosage and metastatic progression (13, 14). Although our mouse studies did not detect an association between Kras focal amplification and metastatic potential, Kras gains (via entire gain of chromosome 6) were detected with comparable frequency in both MetHigh and MetLow tumors, consistent with the ability of both tumor populations to metastasize. Interestingly, among patients with PDAC who had KRAS amplifications or major allelic imbalances, the MYC_TARGETS gene set represented the most highly enriched signature (Fig. 7F) in mirroring the enrichment of this signature in MYC-amplified PDAC tumors (Supplementary Fig. S9F). Signaling through KRAS has long been known to affect MYC expression (44, 80–82), and thus our results are consistent with a model in which elevated MYC activity—as a result of MYC and/or KRAS amplification or some other mechanism—enhances metastatic activity. This interpretation is consistent with recent studies of lung, breast, and prostate cancers identifying a link between MYC amplification and brain or bone metastasis (37, 38, 83).
Although tumor formation in the KPCXY model results from shared founder mutations (KrasG12D activation and Trp53 loss), our genomic analysis revealed ongoing somatic events during tumor progression, resulting in heterogeneous patterns of genomic alterations within a given tumor. Such alterations were largely present at the level of copy-number gains and losses rather than point mutations or small insertions/deletions. Although the complexity of genomic rearrangements varied between tumors, the degree of genome instability did not correlate with metastatic burden. Thus, the increase in metastasis observed in MetHigh tumors is not a function of overall SCNA burden but is instead specific to MYC.
Although subregions within a tumor shared many genomic alterations, consistent with a clonal origin, distinct copy-number alterations were also present, suggesting ongoing sub-clonal evolution. Clonally related metastases exhibited unique (“private”) alterations; however, most copy-number gains and losses were shared with the parental primary tumor clone, suggesting that they were present prior to dissemination. In this respect, it is noteworthy that MYC amplifications in human PDAC were far more common in metastases than in primary tumors, including one case in which we were able to trace a metastatic lesion directly to a MYC-amplified subclone in the primary tumor. Collectively, these results suggest that subclonal MYC amplifications, which have been observed in human primary tumors, provide a selective advantage during metastatic progression (84, 85). Given that our analysis identified a Myc signature in tumor MetHigh clones without Myc amplifications, and MYC amplifications are present in only 11% of human PDAC metastases, other mechanisms for the increased expression of MYC mRNA in PDAC metastases are likely to exist.
As one of the best-studied oncogenes, MYC has been associated with multiple tumor-promoting activities (80). Given MYC's role in tumor cell growth and proliferation, one possible explanation for our results is that MetHigh tumors had an earlier onset and/or grew more rapidly, leading to increased metastasis by mass effect. Against this possibility, we found that tumor size and proliferation rates showed no correlation with metastatic burden in either mouse models or human patients. Likewise, our MetHigh and MetLow cell lines exhibited dramatic differences in metastatic ability despite giving rise to primary tumors of comparable size. Prior studies have shown that Myc overexpression in the pancreas in the context of tumor initiation, without KRAS activation, does not result in PDAC but rather insulinomas (86). By contrast, our work provides evidence that MYC hyperactivation—particularly in the setting of focal, high-amplitude amplification—confers metastatic properties after a primary PDAC is established. These results suggest that genetic context (e.g., mutant Kras status) and timing (e.g., early versus late) may determine whether enhanced tumor growth or promotion of metastasis is the prevalent phenotypic consequence of MYC hyperactivation.
Our studies implicate non–cell-autonomous mechanisms involving the recruitment of TAMs as contributors to metastatic heterogeneity. The ability of TAMs to promote tumor cell invasion is well documented (21, 26, 28, 47, 53, 56–61), a property that is in agreement with our finding that MetHigh and Myc_OE tumors exhibit enhanced vascular intravasation. Furthermore, MYC expression in tumor cells is known to shape the makeup of the surrounding immune microenvironment, making it more immunosuppressive (39, 45). In line with these observations, we find that MetHigh tumors are enriched for alternatively activated TAMs and have decreased T-cell infiltration—features that favor metastasis (21, 26, 28, 30, 31, 47, 53, 56–61). Our data thus support a model wherein stochastically arising tumor subclones with elevated levels of MYC alter the TiME to facilitate intravasation and metastasis.
Although the specific molecular mediator(s) of TAM recruitment remain to be fully elucidated, we speculate that MYC acts indirectly by regulating the expression of factors that regulate TAM migration and/or function. Consistent with this, we identified several chemokines/cytokines that were upregulated in MetHigh tumors and could also be induced by MYC overexpression. Functional validation of these factors identified CXCL3 and MIF as potential mediators of TAM recruitment and associated metastasis. Furthermore, combined inhibition of MIF and CXCR2 significantly reduced both metastasis and TAM invasion. Thus, we hypothesize that multiple secreted factors, as opposed to a single factor, act in concert to drive the MYC-associated increase in proinvasive macrophages.
Most patients with PDAC develop metastases. Our data show that even within this population, the extent of metastatic disease varies widely between patients and affects survival. Although many steps are required for tumor cells to metastasize, our data indicate that bloodstream invasion may be a rate-limiting event for metastasis in PDAC. Our work further suggests that in addition to its well-documented cell-autonomous role in tumor growth, MYC acts non–cell autonomously to promote metastasis. Supporting this, modest overexpression of MYC was previously shown to suffice with KRAS activation to drive PDAC metatsasis in an autochthonous mouse model (44). Given that MYC family members are focally amplified in 28% of human cancers (87), these results have broad implications for metastasis in tumor types other than PDAC (88).
All experiments were performed in accordance with the NIH policies on the use of laboratory animals and approved by the Institutional Animal Care and Use Committee of the University of Pennsylvania. KPCX mice were generated through a series of backcrosses as previously described (32). The RosaConfetti (“X”) reporter allele was introduced into mutant strains bearing Pdx1CreER (“C”), KrasG12D (“K”), and Trp53fl/+ (“P”) alleles to obtain Pdx1CreER; KrasG12D; Trp53fl/+; RosaConfetti (“KPCX”) mice (32). For most experiments, animals were heterozygous for the confetti reporter and also contained a RosaYFP allele in lieu of the second confetti allele to generate “KPCXY” mice. The YFP reporter was introduced to enable fluorescent lineage labeling of tumors that undergo a “no-color” recombination event in the confetti reporter, as previously described (89). To induce recombination, a suspension of tamoxifen (MP Biomedicals) in corn oil (Sigma-Aldrich) was administered to pups via lactation following oral gavage of the mother with 6 mg of the drug on postnatal days 0, 1, and 2. On average, tumor-bearing KPCXY mice were 14 to 16 weeks of age at time of sacrifice.
Multicolor Image Analysis
Pancreatic tumors and organs from tumor-bearing KPCXY mice were isolated and analyzed by fluorescent stereomicroscopy using a Leica M216FA fluorescent microscope with CFP, YFP, and dsRED filters (Chroma). As previously described (32), distinct colorimetric tumor clones in the primary tumor mass are defined as an anatomically contiguous region of monochromatic cells that share a border with adjacent clones of a different color (Fig. 1E; Supplementary Fig. S1C). To accurately quantify the contribution of different colorimetric tumors to metastases, we used the following criteria to identify KPCXY mice suitable for analysis: (i) presence of at least one metastatic lesion to the liver and/or lung, (ii) two or more tumors present, (iii) each metastatic primary tumor carries a unique fluorescent color, and (iv) metastatic lesions can be linked to a specific tumor based on a shared unique fluorescent lineage label. Using these criteria, we identified a panel of 30 mice with a total of 85 tumors. Metastases were quantified by fluorescent stereomicroscopy. Tumor size was determined using ImageJ (NIH) to measure the largest circumference of each fluorescent tumor.
Murine Tumor and Metastasis Sample Acquisition
Pancreatic tumors and associated liver and lung tissues were isolated from tumor-bearing KPCXY mice. Under fluorescent stereomicroscopy, individual colored tumors were identified and biopsied using a 6-mm punch biopsy. Initial biopsy specimens were placed in 750 μL RNAlater (Sigma-Aldrich) for downstream nucleic acid isolation. Subsequent biopsy specimens were submitted for cell line generation and histology. In tumors in which sufficient tissue was available, additional biopsy specimens were taken from anatomically distinct regions of the tumor to obtain subclonal biopsy specimens for genomic analysis. From seven KPCXY mice, we obtained biopsies from 20 tumors, eight of which were amenable to additional subregional biopsies. Paired primary tumors and metastases from each mouse were identified by shared fluorescent lineage labels. Metastases were harvested by microdissection under fluorescent stereomicroscopy and split, and portions were placed in 500 μL RNAlater for nucleic acid isolation or used directly for cell line generation. The remainder of the tissue was embedded for histology. In total, 56 metastases were isolated for genomic analysis. Although individual liver metastases were of sufficient size for microdissection, lung lesions typically were microscopic and could not be readily isolated for molecular analysis.
Tumor Digestion and Cell Lines
Pancreatic tumors were dissociated into single-cell suspensions through mechanical separation and enzymatic digestion as previously described (90). Murine PDAC cell lines 471_MetHigh_1, 832_MetHigh_1, 836_MetHigh_1, 850_MetHigh_4, 852_MetHigh_1, 853_MetHigh_1, 471_MetLow_2, 832_MetLow_2, 842_MetLow_2, 850_MetLow_1, and 852_MetLow_2 were derived from KPCXY primary tumors that were also evaluated by SCNA and RNA-seq. Murine cell lines were cultured in DMEM/F12 medium supplemented with 5 mg/mL D-glucose (Invitrogen), 0.1 mg/mL soybean trypsin inhibitor type I (Invitrogen), 5 mL/L insulin-transferrin-selenium (ITS Premix; BD Biosciences), 25 μg/mL bovine pituitary extract (Gemini Bio-Products), 5 nmol/L 3,3′,5-triiodo-L-thyronine (Sigma-Aldrich), 1 μmol/L dexamethasone (Sigma-Aldrich), 100 ng/mL cholera toxin (Sigma-Aldrich), 10 mmol/L nicotinamide (Sigma-Aldrich), 5% Nu-serum IV culture supplement (Thermo Fisher Scientific), and antibiotics (gentamicin 150 μg/mL, Gibco; amphotericin B 0.25 μg/mL, Invitrogen) at 37°C, 5% CO2, 21% O2, and 100% humidity. Cell lines were maintained and passaged according to ATCC recommended procedures and regularly tested for Mycoplasma using the MycoAlert Mycoplasma Detection Kit (Lonza).
Immunofluorescence and Histologic Analysis
Tissue samples were fixed in 4% paraformaldehyde (EMS) at room temperature for 45 minutes, followed by an overnight incubation in 30% weight/volume sucrose solution (Sigma-Aldrich). Samples were then embedded in optimal cutting temperature (Tissue-Tek) and frozen on dry ice. Staining was performed on 10-μm sections by first blocking with 5% donkey serum and 0.1% Tween-20 for 1 hour, followed by overnight incubation with primary antibody diluted in blocking buffer in a humidified chamber. Sections were washed three times in PBS containing 0.1% Tween-20. For immunofluorescence staining, slides were then incubated with DAPI (Life Technologies, 1:1,000) and Alexa fluorophore–conjugated antibodies (Jackson ImmunoResearch). For IHC, slides were first incubated with biotinylated secondary antibodies (Jackson ImmunoResearch) and developed using the ABC HRP and DAB kits per manufacturer protocols (Vectorlabs). Primary antibodies used were as follows: rat anti–Ki-67 (eBioscience, 14–5698–82), rabbit anti–c-Myc (Y69; Abcam, Ab32072), rabbit anti-CD3 (Invitrogen, PA1–29547), rabbit anti-F4/80 (Novus, NBP2–12506), rat antineutrophil (Abcam, NIMP_R14), anti-CD206 (R&D Systems, AF2535), and rabbit anti-Arg1 (Cell Signaling Technology, 93668).
Bone marrow immune cells were isolated as described (91), and 2 × 106 isolated immune cells were plated in a 6-well dish in Iscove's Modified Dulbecco's Medium (IMDM; Gibco, 12440053) supplemented with 10% FBS (Gibco), 1% L-glutamine (Corning, MT25005CI), 1% nonessential amino acids (Corning, 11140076), 1% sodium pyruvate (Gibco, 11360–070), 0.001% 2-mercaptoethanol (Gibco, 21985023), 1% penicillin–streptomycin (Gibco, 15140163), and 20 ng/mL recombinant murine macrophage colony-stimulating factor (M-CSF; PeproTech, 315–02) at 37°C, 5% CO2, 21% O2, and 100% humidity. BMDMs were differentiated for 7 days and used by gentle scraping before 10 days.
Macrophage Transwell Migration Assay
Macrophage invasion was assessed using a 12-well transwell chamber with 5- to 8-mm filter inserts (Corning). Tumor cells were plated in the bottom chamber 24 hours before addition of BMDMs. For Figs. 5B and 6F, tumor cells were plated on growth factor–reduced Matrigel (Corning, 356231), which was diluted 1:1 in PBS and plated onto the transwell. In total, 100,000 BMDMs were plated per transwell. After 24 hours, the nonmigrated BMDMs and Matrigel were gently removed with a swab. Cells in the bottom surface (migrated BMDMs) of the membrane were fixed in 4% paraformaldehyde for 15 minutes. DAPI (Fig. 5B; Invitrogen, D21490) or crystal violet stain (Fig. 6F; Sigma-Aldrich, 65092a-95) was added in PBS to the transwells. The membranes were imaged and number of macrophages counted in four random fields. For MIF and CXCR2 inhibition, AZD5069 (MedKoo, 206473; 0.1 mmol/L in DMSO) and ISO-1 (Santa Cruz, sc-204807B; 200 μm in DMSO) were added to tumor cells in transwells prior to coincubation with the macrophages. The experiment was performed in triplicate and repeated twice.
Macrophage Depletion In Vivo
A total of 10,000 tumor cells were orthotopically injected into the pancreas of NOD.SCID mice. Treatments started 10 days after implantation. GW2580, a CSF1R inhibitor (AdooQ Bioscience, A11959), was dissolved in 0.5% hydroxypropyl methyl cellulose and 0.1% Tween (HPMT) and dosed three times a week at 160 mg/kg by oral gavage. Then, 200 μL clodronate or control liposomes (Liposoma, CP-025–025) was given once per week by intraperitoneal injection. Blood was sampled to confirm depletion. Experimental and control mice were euthanized 14 days after initiation of treatment and analyzed for metastasis and immune cells by flow cytometry. Control and experimental groups were run in at least triplicate.
MIF and CXCR2 Inhibition In Vivo
In total, 50,000 tumor cells were orthotopically injected into the pancreata of NOD.SCID mice. Treatments started 10 days after implantation. The Mif inhibitor ISO-1 (Santa Cruz, sc-204807B) was dissolved in DMSO to make a 0.05-mg/μL stock solution and then diluted to 10 mg/kg in HPMT solution for gavage five times a week. The CXCR2 inhibitor AZD5069 (MedKoo, 206473) was dissolved in HPMT and dosed five times a week at 0.1 mg/kg by oral gavage. Experimental and control mice were euthanized 14 days after treatment start and analyzed for metastasis and immune cells by immunoflouresence staining. Control and experimental groups were run in at least triplicate.
The iTEM assay was performed as previously described (50–52). Briefly, transwells from EMD Millipore (MCEP24H48) were coated with 2.5 μg/mL Matrigel (356230, BD Biosciences) in a total volume of 50 μL. Then approximately 1 × 104 human umbilical vein endothelial cells (Lonza) in 50 μL EGM2 medium were plated on the inverted transwells previously coated with Matrigel and allowed to adhere for 4 hours at 37°C. Transwells were then placed into a 24-well plate with 1 mL EGM2 with all supplemental factors (Lonza) in the bottom well and 200 μL inside the top chamber and allowed to grow for 48 hours in order to form a monolayer. Pancreatic tumor cells (EV or Myc_OE) were labeled with CellTracker green dye and macrophages (Bac1.2F5) with CellTracker red (green, C7025, red, C34552; Invitrogen), resuspended in DMEM media (SH30253.01, Hyclone) without serum, and plated at 15,000 pancreatic cancer cells without macrophages or with 60,000 macrophages per transwell and allowed to transmigrate toward EGM2 containing 36 μg/mL of CSF1 for 4 hours. Samples were then fixed in 4% paraformaldehyde for 15 minutes, permeabilized with 1% Triton-X 100 for 5 minutes, and stained with ZO-1 (Sigma-Aldrich) to determine and locate the endothelial mono-layer formation. Transwell membranes were cut from the transwell chambers and mounted on a slide using Prolong Diamond antifade (Thermo Fisher Scientific). The slides were imaged using a Leica SP5 confocal microscope using a 60× 1.4 NA objective and processed using ImageJ (NIH). Quantitation was performed by counting the number of tumor cells that had crossed the endothelium within the same field of view (60×, 10 random fields) and represented as normalized values from at least three independent experiments.
Analysis of RNA-seq, Differential Gene Expression, GSEA, and Molecular Subtype
RNA-seq was performed on bulk tumor and metastasis samples from seven KPCX mice, resulting in 66 samples for analysis (primary tumor with subregional biopsy specimens and metastasis). RNA purity and integrity were verified on the Agilent Tapestation prior to library construction, followed by paired-end 50- to 75-bp sequencing on an Illumina HiSeq 4000 high-throughput sequencer. Alignment of fastq files was performed with STAR aligner v2.5.2b using mm10 as the reference genome (92). Gene-level expression data in terms of expected counts and fragments per kilobase of exon per million (FPKM) were obtained using RSEM v1.2.28 (93). Low-expressing genes were removed using a cutoff of 100 for count and 10 for FPKM. For primary tumor clones in which subregional biopsy specimens were taken, count data for the tumor were obtained by merging expression data of the subclone count data using the mean value. This resulted in a total of 54 samples for downstream analysis (MetLow = 13, MetHigh = 7, and metastasis = 34). For differential expression analysis, count data were normalized using the voom function in the limma R package followed by batch correction using the ComBat R package (94, 95). Then limma was used to perform differential expression between MetHigh, MetLow, and metastasis. Box plots of log2 FPKM values for genes were generated using the ggplot2 R package. To generate volcano plots, differential expression data comparing MetLow and MetHigh clones were plotted using ggplot2 with log base twofold change from MetLow compared with MetHigh tumors of each gene plotted on the x-axis and the adjusted P values plotted on the y-axis. Genes with adjusted P values < 0.01 and absolute log2 fold change >1 were highlighted. Differentially expressed genes were used as input for MSigDB GSEA (96). Transcription factor enrichment on all differentially expressed genes was performed using the Metacore software package (Clarivate Analytics). Network analysis was performed on all differentially expressed genes using Ingenuity Pathway Analysis software (Ingenuity Systems, Inc.). Molecular subtype classification using the Moffitt, Bailey, and Collisson (15, 16, 42) signature was performed on each sample by subtracting the sum of normalized expression for genes corresponding to specific classes within a particular molecular signature. We then took the maximum score observed across each class and assigned it to the samples. Heat maps were generated using differentially expressed genes with adjusted P value <0.05 and absolute logFC >1.
Survival Analysis of TCGA Data
TCGA pancreatic adenocarcinoma (PAAD) expression and patient- and sample-level clinical data were downloaded from cBioPortal. Samples were filtered to those classified as PAAD and having available expression data (162 of 186 samples). To develop a signature gene list associated with the MetHigh phenotype, we filtered differentially expressed genes between MetHigh and MetLow tumors using an adjusted P value <0.05 and absolute logFC >0.58, resulting in a set of genes that were upregulated or downregulated in the MetHigh tumors (736 up and 1,036 down). To calculate a signature score for each TCGA PAAD sample, we first z score normalized the TCGA PAAD expression data and then subtracted the sum of all downregulated gene expression from the sum of all upregulated gene expression values. We divided the signature score into high and low strata using a cutoff score >0. Kaplan–Meier analysis was done to compare survival between the two groups.
Human Stage IV Pancreatic Tumor and Metastasis Imaging Analysis
CT scans were obtained from patients with metastatic PDAC undergoing treatment at the University of Pennsylvania under a protocol approved by the Institutional Review Board (#822028). Patients were filtered to include only those with CT scan imaging of the abdomen and chest with intravenous contrast at the time of diagnosis and prior to any treatment and found to have stage IV disease. In total, 55 patients were included. CT images for each patient were reviewed, and metastatic lesions in liver and lung were counted. All metastases were examined in multiple planes to ensure accurate assessment. Tumor area was pulled from the initial radiologist report and measured at the largest diameter. The cutoff for high- and low-metastasis groups was determined using k-means with n = 2 clusters.
Human Pancreatic Cancer Patient Sample Acquisition with Genomic and RNA-seq Analysis
Sample acquisition resulted from patients recruited as part of the ICGC Pancreatic Cancer Ductal Adenocarcinoma Canadian sequencing initiative or the COMPASS trial as previously described (97). Tissue samples were collected at the University Health Network (Toronto), Sunnybrook Health Sciences Centre (Toronto), Kingston General Hospital (Kingston), McGill University (Montreal), Mayo Clinic (Rochester), or Massachusetts General Hospital (Boston) with written informed consent and approval from institutional review or research ethics boards. WGS and RNA-seq were performed on fresh-frozen tumor tissue samples that were enriched for tumor content by LCM. WGS and RNA-seq were performed at the Ontario Institute of Cancer Research as described previously (97). DNA read alignment and MYC copy-number variations were performed on paired-end whole-exome sequencing reads aligned to human reference genome hg19 using Burrows-Wheeler alignment (BWA) 0.6.2 (98). PCR duplicates were marked with Picard 1.90. Tumor cellularity, ploidy, and copy-number segments were derived using an in-house algorithm, CELLULOID (99). RNA reads were aligned to human reference genome hg38 and to transcriptome Ensemble v84 using STAR v2.5.2a (92). Duplicate reads were marked with Picard 1.121. Raw counts were obtained using HTSeq 0.6.1 (100). Differential gene expression analysis was performed with DESeq2 v.1.14.1 (101) using default settings. Briefly, RNA HTSeq count data were imported to generate a dispersion estimate and a generalized linear model. Wald statistical test was used to compare gene expression of MYC amplified to nonamplified cases. GSEA was performed using genes ranked based on the P value and sign of the log2 fold change from differential gene expression analysis. GSEA was run using GSEA Preranked 4.0.2 with default settings against hallmark gene sets (96). Statistical analyses included pairwise comparisons of quantitative variables performed using the Wilcoxon rank sum test. All tests were two-sided. Analyses were carried out in R 3.3.0.
SCNA in Murine Tumors
DNA purified from dissected murine tumors was processed for Illumina library preparation and sequencing using standard protocols. In brief, isolated DNA (100–1000 ng in total) was sonicated on a Covaris instrument. Sonicated DNA was then end-repaired and ligated to TruSeq dual-index library adapters. Index libraries were subsequently enriched by 10 cycles of PCR amplification, followed by pooling and multiplex sequencing targeting a coverage of roughly 2 million reads per sample (102). For data processing and copy-number inference, sequencing data were processed as previously described with mouse genomic bins computed in a manner similar to human bins (36). In brief, sequencing reads were mapped to the mouse reference genome built mm9. Sequencing reads were indexed and sorted, with PCR duplicates removed. Uniquely mapped reads were counted in each bin and normalized for guanine-cytosine content using lowess smoothing. Normalized read count data were then segmented using circular binary segmentation (103) with the profiles centered around a mean of 1. Chromosomal segments with variance that was above or below the mean were called gains or deletions, respectively. A threshold of 0.2 was used. For hierarchical clustering and lineage reconstruction, an analysis based on copy-number values and alteration breakpoints, in a genome-wide manner, was employed.
Bulk and Single-Cell Analysis of Matched Primary and Metastasis from Human Tissue
For 20 patients, tissue sections from flash-frozen samples were processed for bulk DNA purification using the Qiagen DNeasy Blood and Tissue Kit. Purified DNA was processed as described above for multiplex sequencing. A coverage of 2 million sequencing reads was similarly targeted. For a single case, matched pancreatic primary and liver metastasis tissue was retrieved and processed for single-nuclei isolation as previously described (36). Single nuclei were sorted based on DNA content from both diploid and polyploid populations of each tissue. Approximately 100 nuclei per tissue sample were amplified using the WGA4 kit (Sigma-Aldrich), with the resulting whole-genome amplified DNA processed for TruSeq Indexed sequencing library preparation, as described above. Sequencing data were processed as described above with the exception that a least squares fitting algorithm was used to calculate absolute integer copy number (102).
Quantification and Statistical Analysis
Statistical analysis used is indicated in each figure where relevant. All statistical analyses were performed with GraphPad Prism 8 (GraphPad Software). K-means clustering and survival analysis were carried out in R 3.3.0. Error bars show standard deviation or standard error of the mean shown as indicated in the legend, and P < 0.05 was considered statistically significant (*, P < 0.05, **, P < 0.005, and ***, P < 0.0001) unless otherwise indicated. ns denotes not significant.
Data and Code Availability
The datasets generated during this study are available through the Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra) under the accession numbers PRJNA647834, PRJNA646123, and PRJNA646156. Additional primary data needed for review or to replicate the findings will be made available on request. Code generated during this study is available at https://github.com/rmaddipati79/Maddipati_PDAC_metastasis.git. Additional detailed methods can be found in the Supplementary Materials and Methods.
C.A. Adkisson reports a Ruth L. Kirschstein Institutional National Research Service Award during the conduct of the study. J.R. Pitarresi reports grants from the NIH/NCI K99CA252153 and F32CA221094 and the NIH LRP, an American Gastroenterological Association Bern Schwartz Research Scholar Award, a Hopper Belmont Foundation Inspiration Award, and a UPenn P30 Center for Molecular Studies in Digestive and Liver Diseases Pilot and Feasibility Grant outside the submitted work. E.L. Carpenter reports other support from Merck, personal fees from Gordon Conference, Imedex, AHPBA, Bristol Myers Squibb, and GuardantHealth, and grants from UHG outside the submitted work. J.C. McAuliffe reports personal fees from Boston Scientific and nonfinancial support from Loki Therapeutics outside the submitted work. S.W. Lowe reports personal fees and other support from ORIC Pharmaceuticals and Blueprint Medicines, other support from Faeth Therapeutics, and personal fees from PMV Pharmaceuticals outside the submitted work. C.A. Iacobuzio-Donahue reports other support from Bristol Myers Squibb outside the submitted work. B.Z. Stanger reports grants from the NIH during the conduct of the study, as well as personal fees from iTeos Therapeutics and grants from Boehringer Ingelheim outside the submitted work. No disclosures were reported by the other authors.
R. Maddipati: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, writing–original draft, project administration, writing–review and editing. R.J. Norgard: Conceptualization, data curation, investigation, writing–review and editing. T. Baslan: Conceptualization, resources, data curation, investigation, visualization, methodology, writing–original draft, writing–review and editing. K.S. Rathi: Data curation, formal analysis, visualization, writing–review and editing. A. Zhang: Data curation, formal analysis, visualization, writing–review and editing. A. Saeid: Data curation, writing–review and editing. T. Higashihara: Data curation, writing–review and editing. F. Wu: Data curation, writing–review and editing. A. Kumar: Data curation, writing–review and editing. V. Annamalai: Data curation, writing–review and editing. S. Bhattacharya: Data curation, writing–review and editing. P. Raman: Conceptualization, data curation, supervision, investigation, visualization, writing–review and editing. C.A. Adkisson: Data curation, writing–review and editing. J.R. Pitarresi: Data curation, formal analysis, writing–review and editing. M.D. Wengyn: Data curation, writing–review and editing. T. Yamazoe: Data curation, writing–review and editing. J. Li: Data curation, writing–review and editing. D. Balli: Data curation, formal analysis, writing–review and editing. M.J. LaRiviere: Resources, data curation, visualization, writing–review and editing. T.-V.C. Ngo: Data curation, writing–review and editing. I.W. Folkert: Data curation, writing–review and editing. I.D. Millstein: Data curation, writing–review and editing. J. Bermeo: Data curation, writing–review and editing. E.L. Carpenter: Resources, data curation, writing–review and editing. J.C. McAuliffe: Data curation, supervision, writing–review and editing. M.H. Oktay: Data curation, supervision, writing–review and editing. R.A. Brekken: Data curation, supervision, writing–review and editing. S.W. Lowe: Resources, supervision, writing–review and editing. C.A. Iacobuzio-Donahue: Resources, supervision, writing–review and editing. F. Notta: Resources, formal analysis, supervision, writing–review and editing. B.Z. Stanger: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, writing–original draft, project administration.
We are grateful for helpful advice from Gregory Beatty, Chi Dang, David DeNardo, Rosalie Sears, and all members of the Stanger Laboratory. This work was supported by grants from the NIH (CA229803 to B.Z. Stanger, DK109292 to R. Maddipati, and CA236269 to R.J. Norgard), Cancer Prevention & Research Institute of Texas (RR190029 to R. Maddipati), and American Gastroenterological Association (Caroline Craig and Damian Augustyn Award in Digestive Cancer to R. Maddipati). J.C. McAuliffe is supported by grants from the Albert Einstein Cancer Center (#305631) and a Ruth L. Kirschstein T32 (CA200561). T. Baslan is supported by the William C. and Joyce C. O'Neil Charitable Trust, Memorial Sloan Kettering Single Cell Sequencing Initiative. F. Notta is supported by the Ontario Institute for Cancer Research (OICR), the Canadian Institutes of Health Research (no. 388785), the Cancer Research Society (no. 23383), and the Gattuso-Slaight Personalized Cancer Medicine Fund from Princess Margaret Cancer Centre. Additional support was provided by the Abramson Family Cancer Research Institute, the Abramson Cancer Center, and the NIH/Penn P30 Center for Molecular Studies in Digestive and Liver Diseases (P30DK050306).
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 Discovery Online (http://cancerdiscovery.aacrjournals.org/).