BLIMP1 Induces Transient Metastatic Heterogeneity in Pancreatic Cancer

Pancreatic ductal adenocarcinoma (PDAC) is an almost uniformly lethal cancer that is projected to become the second leading cause of cancer-related deaths in the United States by 2030 (1). Most patients with PDAC die from metastatic disease, underscoring the need to better understand the molecular mechanisms that drive disease progression and metastasis (2). Genomic analyses of PDAC have uncovered oncogenic KRAS and loss-of-function mutations in the CDKN2A, SMAD4, and TP53 tumor suppressors as key recurrent drivers of pancreatic cancer development (3–6). Although these studies have offered clues about metastatic progression, they have not uncovered consistent genetic alterations that explain the progression to a highly metastatic state (7–10).

Although genomic alterations create stable changes that increase cancer growth, transient alterations in the metastatic state of cancer cells can be induced by interactions with stromal cells, diverse physical cues, as well as by changes in the local tumor microenvironment. For example, the epithelial-to-mesenchymal transition (EMT) is a well-characterized transcriptional program that endows cancer cells with transient high metastatic ability (11). However, EMT might not be critical for PDAC dissemination or metastasis (12, 13). Subpopulations of PDAC cells with cancer stem cell (CSC)–like properties have also been described, but it is unclear whether these cells are the major source of metastases (14, 15).

In many cancer types, metastasis is thought to be driven by diverse extracellular cues that increase stem-like behavior as well as invasion and metastasis (16). PDAC in particular has an extensive desmoplastic stromal response that generates unique physical properties, including increased extracellular matrix stiffness and areas with limited oxygen and nutrient availability (17). However, whether PDAC metastasis is driven by features of the tumor microenvironment is unclear. Identification of key environmental factors could provide insights into the process of metastasis as well as aid in the development of novel therapeutic strategies.

Genetically engineered mouse models of PDAC recapitulate key genetic events of the human disease. Cre-mediated expression of oncogenic KRAS^G12D in the pancreatic cells of loxP-Stop-loxP Kras^G12D knock-in mice (Kras^LSL-G12D/+) leads to the development of early-stage pancreatic intraepithelial neoplasms (PanIN; ref. 18). Concomitant expression of point mutant Trp53 or deletion of Trp53, Cdkn2a, and/or Smad4 allows for the development of PDAC that can progress to gain multiorgan metastatic ability (19–23). Importantly, tumors arise in vivo from genetically defined lesions and evolve in their native context, providing the opportunity to identify the cancer cell–intrinsic and –extrinsic processes that contribute to tumor progression.

Here, we developed a novel mouse model of human PDAC, which enabled the isolation and molecular characterization of a highly metastatic subpopulation of pancreatic cancer cells. We demonstrate that these highly metastatic cancer cells exist within hypoxic tumor areas and that the transcription factor BLIMP1 drives their high metastatic potential. Gene expression signatures of the metastatic state, as well as of hypoxia-induced BLIMP1-dependent genes, predict PDAC patient outcome. These findings highlight microenvironment-induced heterogeneity as a driver of pancreatic cancer progression toward its deadly metastatic phase.

The chromatin-associated protein HMGA2 is a marker of increased malignancy in many tumor types, and high HMGA2 expression predicts poor prognosis in several major human cancer types, including PDAC (24–30). To determine whether neoplastic cells in genetically engineered mouse models of human PDAC also express HMGA2, we performed immunohistochemistry (IHC) on tumors at different stages of development. HMGA2 was not expressed in cells in the normal adult pancreas or PanINs in Kras^LSL-G12D;Trp53^LSL-R172H/+;Pdx1-Cre (KPC) mice, but was expressed in a subset of PDAC cells (Supplementary Fig. S1A and data not shown). In human PDAC, HMGA2 expression correlates with metastasis to lymph nodes and poor prognosis, and we confirmed that high HMGA2 expression in patients with PDAC predicts shorter survival (Supplementary Fig. S1B–S1D; refs. 28, 31). Together, these results document the expression of HMGA2 in a subset of cancer cells in mouse models of PDAC and confirm the correlation of the presence of cancer cells in the HMGA2⁺ state with poor outcome in patients with PDAC.

To uncover the cellular and molecular features of HMGA2⁻ and HMGA2⁺ cancer cells, we generated a mouse model that would allow the isolation of these PDAC cell subpopulations. We incorporated two additional alleles into the KPC mouse model: a Cre-reporter allele (R26^LSL-Tomato) to fluorescently mark all neoplastic cells, and an Hmga2 knock-in allele, which is converted by Cre from its wild-type conformation (Hmga2^CK) into a GFP reporter (Hmga2^eGFP; Fig. 1A; refs. 20, 32). In the heterozygous state (Hmga2^CK/+), the potential for GFP expression is restricted to cells in which Cre has inverted a loxP-flanked region and GFP expression remains under control of all endogenous Hmga2 regulatory elements (20). In KPC;R26^LSL-Tomato/+;Hmga2^CK/+ mice (referred to as KPC^colors mice), all cancer cells were Tomato positive, and HMGA2-expressing cancer cells were both Tomato and GFP positive (Supplementary Fig. S1A and S1E).

The dual fluorescent marking of cancer cells in KPC^colors mice provided us with the ability to isolate Tom⁺GFP⁻ and Tom⁺GFP⁺ cancer cells by fluorescence-activated cell sorting (FACS; Supplementary Fig. S1F and S1G). Consistent with HMGA2 expression observed by IHC, variable percentages of cancer cells in individual tumors were GFP⁺ (Fig. 1B and C). In the KPC model, progression from PanINs to adenocarcinoma is driven by loss of the wild-type (WT) Trp53 allele (19). Tom⁺GFP⁻ and Tom⁺GFP⁺ samples contained less than 10% remaining Trp53^WT allele, and loss of the Trp53^WT allele led to the stabilization of mutant p53 protein in both GFP⁻ and GFP⁺ cells (Supplementary Fig. S1H–S1J). Thus, Tom⁺GFP⁻ and Tom⁺GFP⁺ cells represent two distinct subpopulations of pancreatic cancer cells.

We next performed cell culture and transplantation-based in vivo metastasis assays on GFP⁻ and GFP⁺ PDAC cells. GFP⁺ cells consistently formed more spheres when plated into ultra-low-attachment plates and formed more colonies when plated at low density under standard tissue culture conditions (Supplementary Fig. S1K and data not shown). Most importantly, for 8 out of 8 tumors from KPC^colors mice, the GFP⁺ PDAC cells formed more metastases than their GFP⁻ counterparts when transplanted intravenously into recipient mice (Fig. 1D–F). On average, GFP⁺ cells were more than 10 times more metastatic than GFP⁻ cells (P < 0.008; Fig. 1F). Interestingly, the tumors that arose from GFP⁺ cells almost always had heterogeneous GFP expression, suggesting that GFP⁺ cells may be in a transient state with the potential to give rise to both GFP⁻ and GFP⁺ cells (Supplementary Fig. S1L and S1M).

To uncover prometastatic programs within the highly metastatic GFP⁺ PDAC cell state, we performed RNA sequencing (RNA-seq)–based gene expression profiling on six pairs of GFP⁻ and GFP⁺ cells (Fig. 2A and Supplementary Fig. S2A and S2B). Global clustering of all samples did not clearly separate GFP⁻ from GFP⁺ samples (Fig. 2B). However, direct pairwise comparison of GFP⁻ and GFP⁺ cells uncovered hundreds of genes with consistent and significant differences (Fig. 2C and D). Neither canonical epithelial markers nor genes related to EMT were consistently different between GFP⁻ and GFP⁺ cells (data not shown). We also did not observe enrichment for previously described gene expression signatures of PDAC metastasis or putative CSCs in GFP⁺ cells (refs. 13, 14; data not shown). Using flow cytometry, we confirmed that both GFP⁻ and GFP⁺ cancer cells had heterogeneous expression of the ductal/CSC marker CD133 and the epithelial marker EPCAM (Supplementary Fig. S2C and S2D; refs. 33–35). Histologic features and IHC for differentiation markers confirmed that Hmga2 expression is largely independent from differentiated state (Fig. S2E).

In addition to the paired GFP⁻ and GFP⁺ PDAC samples, we performed RNA-seq analyses on FACS-sorted, bulk Tom⁺ cancer cells from primary tumors and metastases (Fig. 2A and Supplementary Fig. S2F). If metastases had stable gene expression differences from primary tumors, this approach could identify gene expression alterations that contribute to metastatic ability or growth at secondary sites. Interestingly, comparison of primary tumors to all metastases identified very few significant differentially expressed genes (Supplementary Fig. S2G). Comparisons of primary tumors to liver metastases, but not to lymph node metastases, uncover several genes that were significantly differentially expressed in the liver metastases (Supplementary Fig. S2H and S2I). Consistent with a recent report on human PDAC metastasis (36), gene set analysis uncovered a trend toward enrichment for programs related to glucose metabolism in liver metastases (Supplementary Fig. S2J). Importantly, genes that were differentially expressed between GFP⁻ and GFP⁺ PDAC cells were not consistently different between primary tumors and metastases, consistent with the transient nature of the GFP⁺ cell state (Fig. 2E). Finally, high expression of a gene signature composed of genes that were more highly expressed in metastatic GFP⁺ cancer cells predicted worse outcome in patients with PDAC (Fig. 2F and G).

To gain further insight into the metastatic process and identify potentially prometastatic factors, we focused on several genes that were among the most significantly and dramatically upregulated in GFP⁺ cells (fold change > 2; P < 10⁻⁶; Supplementary Fig. S3A). We stably knocked down five top candidate genes (Ero1l, Slc16a3, Glut1, Hilpda, and Blimp1/Prdm1) in a PDAC cell line (688M) derived from liver metastasis from a KPC;R26^LSL-Tomato/+ (KPCT) mouse (Supplementary Fig. S3B–S3F). We assessed the importance of these genes in metastasis by quantifying the number of metastases that formed from subcutaneously and orthotopically transplanted tumors. These experiments suggested that the transcription factor BLIMP1/PRDM1 could have prometastatic functions in PDAC (Supplementary Fig. S3G–S3L). BLIMP1 is a transcription factor that was of particular interest due to its well-established role as a master regulator of cell fate determination during plasma B-cell differentiation and primordial germ cell development (37, 38). Blimp1 was one of the most highly upregulated genes in GFP⁺ cells, being 4- to 27-fold higher in GFP⁺ cells (P < 0.05; Fig. 3A). We confirmed increased BLIMP1 protein expression in sorted GFP⁺ cells relative to GFP⁻ cells (Fig. 3B). Blimp1 expression was not consistently different between bulk cancer cells from primary tumors and metastases from the KPCT mice, consistent with the unstable nature of the metastatic state (Supplementary Fig. S3M).

To further assess whether BLIMP1 contributes to metastatic ability, we knocked down Blimp1 using two independent shRNAs in 688M cells (Supplementary Fig. S3N). Blimp1 knockdown reduced the number of metastases seeded from subcutaneous tumors by >50-fold (P < 0.005; Fig. 3C–E). Blimp1 knockdown in a second metastasis-derived PDAC cell line (1004M) also significantly reduced metastatic ability (Fig. 3F and G and Supplementary Fig. S3O). Interestingly, although Blimp1 appeared to be required for metastatic ability, overexpression of BLIMP1 in multiple PDAC cell lines did not consistently enhance metastatic ability, suggesting that it is not sufficient to drive PDAC metastasis (Supplementary Fig. S3P–S3S).

We next used a Blimp1 conditional knockout allele to investigate BLIMP1 function in autochthonous PDAC (37). Blimp1^flox/flox;Pdx1-Cre mice were viable and their pancreata did not show obvious histologic changes, suggesting that Blimp1 is not required for pancreas development or homeostasis (data not shown). KPCT;Blimp1^flox/flox mice had similar overall pancreatic tumor burden but shorter survival compared with control KPCT;Blimp1^+/+ mice (Supplementary Fig. S3T and S3U). Pancreata from KPCT;Blimp1^flox/flox mice contained PanINs as well as adenocarcinomas that were similar to PDAC in control KPCT;Blimp1^+/+ mice (Supplementary Fig. S3V and S3W and data not shown). To assess the effect of Blimp1 deficiency on metastatic progression in vivo, we carefully quantified the number of Tom⁺ disseminated tumor cells (DTC) in the peritoneal cavity as well as metastases in KPCT;Blimp1^flox/flox and control mice. Fourteen out of 15 control mice (KPCT;Blimp1^+/+ and KPCT;Blimp1^flox/+) developed metastases, which were often numerous and widespread in many different sites, including the lymph nodes, diaphragm, lungs, and liver (Fig. 3H–K). Conversely, only 3 out of 14 KPCT;Blimp1^flox/flox mice developed metastases (Fig. 3K). Additionally, peritoneal DTCs could be detected in only half of KPCT;Blimp1^flox/flox mice, but were present in all control mice (Fig. 3I–K). Together with our observations from cell lines, these data suggest that Blimp1 promotes metastatic proclivity of PDAC.

To place BLIMP1 in a pathway involved in metastasis, we next used gene set enrichment analysis (GSEA) and gene ontology (GO) enrichment analysis to identify pathways altered in the more metastatic GFP⁺ cells. These analyses uncovered an overwhelming enrichment for hypoxia-induced genes in GFP⁺ cells (Fig. 4A–D and Supplementary Table S1). Genes expressed more highly in GFP⁺ cells were also enriched for HIF1-binding motifs near their transcription start sites, and our analyses identified significant enrichment of both HIF1 and HIF2 regulated genes in GFP⁺ cells (Fig. 4B; Supplementary Fig. S4A and S4B; Supplementary Table S1). Conversely, genes downregulated in GFP⁺ cells were enriched for cell-cycle processes, consistent with hypoxia-induced cell-cycle arrest (Supplementary Table S1; ref. 39). We confirmed the upregulation of the canonical HIF1 target gene ERO1L at the protein level in sorted GFP⁺ PDAC cells (Supplementary Fig. S4C).

Pimonidazole-defined hypoxic areas were significantly enriched for HMGA2⁺ cells (Supplementary Fig. S4D and S4E). We also employed multicolor sequential immunofluorescence staining to show that Hmga2⁺ areas were enriched for the expression of the canonical hypoxic target protein GLUT1 (Supplementary Fig. S4F and S4G; ref. 40).

Based on the striking enrichment of HIF targets in GFP⁺ PDAC cells from KPC^colors mice, we determined whether Hmga2 expression is regulated by hypoxia. Under hypoxia, we noted only a slight increase in HMGA2 protein levels in PDAC cell lines (Supplementary Fig. S4H). Although HIF target genes were enriched in HMGA2-expressing PDAC cells, Hmga2 knockdown had no effect on the hypoxia-induced expression of canonical HIF1 target genes (Supplementary Fig. S4I). Thus, it remains unclear why HMGA2 marks this highly metastatic PDAC subpopulation, but these data suggest that other aspects of the in vivo microenvironment either in conjunction with, or independent from, hypoxia induce HMGA2 expression in these cells.

To determine whether BLIMP1 expression is regulated by hypoxia in human and murine PDAC, we assessed BLIMP1 mRNA and protein expression in PDAC cell lines exposed to hypoxia (0.5% O₂ for 24 hours). Hypoxia led to the induction of multiple canonical HIF target genes, HIF1α stabilization, and the prominent and consistent induction of BLIMP1 in two mouse and four human PDAC cell lines (Supplementary Fig. S4J and S4K and Fig. 4E and D). Hypoxia-mediated induction of BLIMP1 in mouse and human PDAC cells was attenuated by HIF1α knockdown, suggesting that HIF1α is at least partially required for BLIMP1 induction under these conditions (Fig. 4G and H). BLIMP1 induction in human PDAC cell lines was also partially HIF2 dependent (Supplementary Fig. S4B). Expression of stable HIF1α was sufficient to increase BLIMP1 expression in PDAC cells (Fig. 4I). Finally, human PDACs with the highest BLIMP1 expression are enriched for hypoxia signatures relative to those with the lowest BLIMP expression (Supplementary Fig. S4L and data not shown).

We next investigated how hypoxia and HIF induce Blimp1 expression. To determine whether Blimp1 can be induced indirectly by secreted factors, we measured Blimp1 levels in PDAC cells cultured with conditioned media from hypoxia-treated cells or recombinant VEGFA, which has been shown to induce BLIMP1 in endothelial cells (41). In both cases, we did not observe robust Blimp1 induction (Supplementary Fig. S5A–S5C). Blimp1 was induced rapidly after exposure to hypoxia, paralleling the kinetics of canonical HIF target genes, suggesting that Blimp1 might be induced directly by HIF (Supplementary Fig. S5D). Analysis of chromatin accessibility around the Blimp1 locus (see below) enabled the prioritization of multiple putative distal regulatory regions that contained hypoxia-response elements (HRE; Supplementary Fig. S5E). HIF1a ChIP qPCR identified a cluster of 3 adjacent HREs upstream of Blimp1 that were bound by endogenous HIF1α in PDAC cells under hypoxia (Fig. 4J). This HRE-containing putative distal regulatory region conferred hypoxia responsiveness in a heterologous reporter system, which was abolished by mutation of its HRE motifs (Fig. 4K and Supplementary Fig. S5F and S5G). Furthermore, Blimp1 knockdown significantly reduced the ability of PDAC cells cultured under hypoxia to form spheres and had a variable effect of migratory ability in cell culture (Supplementary Fig. S6A–S6I). These results suggest a role for BLIMP1 in cellular behaviors related to metastatic ability.

To characterize Blimp1's function in hypoxic cells, we profiled the gene expression and genome-wide chromatin accessibility of shControl and shBlimp1 PDAC cells cultured under normoxic and hypoxic conditions (Fig. 5A). Hypoxia can induce changes in chromatin state, and BLIMP1 has been implicated in both plasma cell precursors and primordial germ cells as a regulator of chromatin structure (42–44). We uncovered widespread hypoxia-induced changes in chromatin accessibility, with differentially accessible regions being enriched for HIF-binding elements (Fig. 5B and Supplementary Fig. S6J and S6K). In addition, hypoxia induced genes associated with newly open chromatin regions more than those with constitutively open or closed regions, suggesting that hypoxia likely regulates target gene induction in part through chromatin accessibility changes (Supplementary Fig. S6L). Interestingly, Blimp1 knockdown had minimal impact on hypoxia-induced changes in chromatin accessibility, indicating that the function of BLIMP1 is largely independent of its ability to recruit factors that lead to changes in chromatin state (Supplementary Fig. S6M and S6N).

Our parallel RNA-seq analysis identified many genes that were dramatically and significantly altered by hypoxia (Fig. 5C and D). As expected, canonical genes related to hypoxia were induced, whereas cell cycle–related programs were suppressed (Fig. 5E). Consistent with the induction of Blimp1 by hypoxia, Blimp1 knockdown affected the expression of more genes when the cells were cultured under hypoxic conditions (Fig. 5C; Supplementary Fig. S6O, and comparison between Fig. 5F and G). BLIMP1 was required for both the induction and repression of subsets of hypoxia-regulated genes (Supplementary Table S2). Under hypoxia, cell cycle–related programs were enriched in shBlimp1 cells compared with shControl cells, suggesting that BLIMP1 might play a role in hypoxia-induced cell-cycle arrest (Fig. 5H). Approximately 12% of hypoxia-repressed genes required BLIMP1 for their full suppression (N = 95 of 825 hypoxia-repressed genes; Fig. 5I and Supplementary Fig. S6P and S6Q). The majority of these hypoxia-repressed, BLIMP1-dependent genes were related to cell-cycle processes, consistent with the role of BLIMP1 in suppressing proliferation during plasma B-cell differentiation (Supplementary Fig. S6R; refs. 45, 46).

Additionally, approximately 35% of hypoxia-induced genes required BLIMP1 for their full induction and were less induced under hypoxia in shBlimp1 cells (N = 833 of 2,342 hypoxia-induced genes; Fig. 5J). Genes encoding proteins involved in hypoxic responses and cell mobility were reduced in shBlimp1 cells compared with shControl cells (Fig. 5K; Supplementary Fig. S7A and S7B and Supplementary Table S3). We found that accessible distal regulatory regions within 500 kb of the transcription start sites of BLIMP1-dependent, hypoxia-induced genes were enriched for transcription factor–binding motifs that closely resemble the BLIMP1 motif (ref. 47; IRF1/IRF2; Supplementary Fig. S7C and S7D). Although the regulation of these BLIMP1-dependent genes is likely to be multifaceted, the enrichment of these motifs suggests that at least a subset of these genes may be regulated directly by BLIMP1. Finally, high expression of a gene expression signature composed of hypoxia-induced, BLIMP1-dependent genes predicted worse outcome for patients with PDAC (Fig. 5L). These results suggest that Blimp1 is a hypoxia-regulated gene that regulates a defined subset of hypoxia-controlled genes in PDAC cells.

To gain additional insight into the function of BLIMP1 in PDAC, we integrated our ex vivo RNA-seq data from GFP⁻ and GFP⁺ PDAC cells with our in vitro RNA-seq data from shControl and shBlimp1 cells cultured under normoxia and hypoxia. As anticipated, a vast majority of genes that are more highly expressed in GFP⁺ cells were also upregulated by hypoxia in PDAC cells in cell culture (Fig. 6A and Supplementary Fig. S8A and S8B). Furthermore, many hypoxia-induced genes that were more highly expressed in GFP⁺ cells in vivo required Blimp1 for optimal induction under hypoxic conditions in vitro (Supplementary Fig. S8A and S8C). These results underscore the strong hypoxia signature in the GFP⁺ cells and highlight the contribution of BLIMP1 to the expression of these genes.

To further relate these gene expression programs with BLIMP1 expression in human PDAC, we defined a 36-gene signature of BLIMP1-dependent, hypoxia-induced genes that are also higher in the GFP⁺ state. Across multiple human PDAC gene expression datasets, this BLIMP1 signature correlated with BLIMP1 expression, suggesting conserved mechanism of BLIMP1 function in human PDAC in vivo (Fig. 6B and Supplementary Fig. S8D and S8E).

Our gene expression profiling suggested that BLIMP1 might be required for hypoxia-induced cell-cycle arrest. To directly test this, we cultured shControl and shBlimp1 cells at 0.5% and 20% O₂ and assessed proliferation by short-term BrdUrd labeling. Although shControl cells almost completely arrested under hypoxia, shBlimp1 cells continued to proliferate (Supplementary Fig. S8F and S8G). To determine whether BLIMP1 reduces the proliferation of PDAC cells in tumors in vivo, we assessed the proliferation of cancer cells in pancreatic tumors in KPCT;Blimp1^flox/flox and control KPCT;Blimp1^+/+ mice. Cancer cells in autochthonous Blimp1-deficient tumors had a higher mitotic index (Fig. 6C and Supplementary Fig. S8H). The higher proliferation of cancer cells in tumors from KPCT;Blimp1^flox/flox mice is also consistent with the shorter survival of KPCT;Blimp1^flox/flox mice (Supplementary Fig. S3U).

Many of the genes that were hypoxia-induced, BLIMP1-dependent, and expressed at higher levels in the more metastatic GFP⁺ PDAC cells have been previously implicated as prometastatic factors in other cancer types. These genes included Pgf, Dusp1, Hmox1, Car9, Glut1, and Hilpda (48–53). Consistent with our RNA-seq data, we observed reduced GLUT1 and CAR9 protein expression in PDACs in KPCT;Blimp1^flox/flox mice compared with KPCT mice (Fig. 6D and E and Supplementary Fig. S8I–S8M). High expression of the lipid droplet–associated protein Hilpda in other cancer types correlates with disease progression and metastasis (53, 54). Hildpa expression was higher in GFP⁺ PDAC cells, induced by hypoxia in murine and human PDAC cells, and its induction was partially Blimp1-dependent (Fig. 6F and G and Supplementary Fig. S4J). Hilpda knockdown reduced metastasis in our initial analysis, and we further confirmed that Hilpda knockdown in PDAC cells significantly reduced their metastatic ability (Fig. 6H–J; Supplementary Fig. S8N and Supplementary Fig. S5G–S5L). These data suggest that Hilpda is a Blimp1-regulated prometastatic factor in PDAC.

To uncover molecular mechanisms that contribute to the metastatic ability of PDAC, we initially took two unbiased gene expression-profiling approaches: analysis of HMGA2-GFP⁻ and HMGA2-GFP⁺ PDAC subpopulations as well as analysis of bulk cancer cells from large primary tumors and macrometastases. In both cases, we specifically isolated cancer cells at high purity by FACS to avoid confounding our analyses with contaminating stromal cell populations. Analysis of bulk cancer cells from primary tumors and metastases uncovered few significant gene expression changes, implying that cancer cells in the largest primary tumors possess most of the molecular features required for metastatic spread.

These findings are in stark contrast to the extensive gene expression differences between large primary tumors and metastases that we uncovered in a parallel study on a Kras^G12D-driven, Trp53-deficient mouse model of lung adenocarcinoma (55). In the lung cancer model, large primary tumors often existed in an earlier nonmetastatic state that had profound gene expression differences from metastases. In the lung, oncogenic Kras^G12D alone can drive extensive tumor growth, and even tumors in Kras^LSL-G12D;Trp53^flox/flox mice do not immediately receive benefit from being Trp53 deficient (55–58). Thus, pancreatic tumors may be forced into a potentially metastatic state by the selective pressures of primary tumor growth, thereby explaining the high likelihood of metastatic spread even in patients with relatively small tumors (59).

Despite these observations, multiple lines of evidence suggest that the metastatic ability of PDAC is still an acquired phenotype. We previously noted mice with widespread PanIN lesions that lacked any DTCs in their peritoneal cavities (60). Additionally, we and others have generated mice with clonally marked pancreatic tumors and documented that not all tumors give rise to metastases (60, 61). Although we did not observe gene expression differences between large primary tumors and metastases, we have documented microenvironment-driven metastatic heterogeneity. Our results support a model in which the development of hypoxic regions generates cells with increased metastatic ability (62). Consistent with results from autochthonous mouse models, human PDAC is a highly hypoxic cancer type (63) and the metastatic ability of orthotopically grown, patient-derived PDAC xenografts is predicted by their level of hypoxia (64).

Hypoxia has been shown to induce metastasis in multiple cancer types through various mechanisms (reviewed in refs. 64, 67, 68). Hypoxia has been linked to alterations in EMT/MET, angiogenesis, local invasion and intravasation, and extravasation, as well as the formation of the premetastatic niche (65). Although some consequences of hypoxia may be relatively generalizable across cancer types, some outputs of hypoxia may also be cancer type-specific; thus, the importance of BLIMP1 in these different steps of the metastatic cascades as well as in different cancer types remains to be determined.

Hypoxia also has a tremendous impact on the self-renewal and differentiation of progenitor/stem cell lineages. For example, hypoxia potentiates the engraftment of human hematopoietic stem cells in recipient mice (66, 67) and also helps maintain the stemness of embryonic stem cells and iPS cells in culture (68, 69). In several cancer types, hypoxia has also been shown to play important roles in maintaining CSCs. In brain tumors, hypoxia promotes and/or maintains cancer-cell stemness similar to the effect of hypoxia on bona fide stem cells (70, 71). Several studies have identified subpopulations of murine and human PDAC cells with CSC characteristics based on their ability to generate new tumors upon transplantation (14, 72). Interestingly, the highly metastatic PDAC state that we identified is not directly related to previously reported CSC populations, the differentiation state of the cancer cells, or EMT. Thus, whether the highly metastatic PDAC cell state and these CSC states represent parallel or partially overlapping programs will be an important area for future study.

We initially anticipated that the high metastatic ability of HMGA2-expressing cells would be driven by HMGA2 itself. Surprisingly, this is not the case, as Hmga2 deficiency has no impact on the metastatic ability of tumors in the KPC PDAC mouse model, nor does it influence the induction of canonical hypoxia target genes (BMG, S-HC, MMW; manuscript in preparation and Supplementary Fig. S4I). HMGA2 could play a subtle role in the later stages of metastatic outgrowth or may simply be a marker of the metastatic state.

Mechanistically, our results uncover hypoxia/HIF-mediated induction of the transcription factor BLIMP1 as one molecular link between the tumor microenvironment and transient induction of prometastatic gene expression programs in PDAC. Although our data show that BLIMP1 can be induced through hypoxia-mediated stabilization of HIF, other factors within the tumor microenvironment may also affect HIF activity (Fig. 6K). In PDAC, BLIMP1 functions as a molecular switch that promotes metastatic ability while suppressing cell division under hypoxia (Fig. 6K). Our results are consistent with the link between BLIMP1 and migratory ability of human lung and breast cancer cell lines in vitro (73, 74). Blimp1 has not been described as a hypoxia/HIF target gene in normal cell types, but hypoxia may also influence BLIMP1 expression in those settings. In early embryos, where oxygen levels are low prior to the formation of major blood vessels (75), BLIMP1 is expressed in primordial germ cells, where it represses somatic programs and helps maintain pluripotency (38, 76). BLIMP1 is also critical for the differentiation of plasma cells that are generated in secondary lymphoid organs and maintained in bone marrow, both of which have hypoxic regions (37, 77, 78).

In summary, our findings support the concept of microenvironmental, rather than mutational, evolution being a critical factor that fosters PDAC metastatic ability. We found that intratumoral hypoxia, which is an inevitable feature of advanced human PDAC, induces the expression of the prometastatic transcription factor BLIMP1. The co-option of this master regulatory transcription factor promotes metastatic ability, and the molecular output of BLIMP1 expression is the modulation of discrete hypoxia-induced gene expression programs. A greater understanding of the origins and molecular features of cancer cells with transient high metastatic ability could provide therapeutic opportunities to reduce metastatic spread and further our appreciation of the obligate plasticity of these cells during the metastatic process.

Kras^LSL-G12D, Trp53^LSL-R172H, Blimp1^flox, Pdx1-Cre, Rosa26^LSL-tdTomato, and Hmga2^CK mice have been described (18, 20, 32, 37, 79, 80). Mice with the Kras^LSL-G12D and the R26^LSL-tdTomato alleles in cis on chromosome 6 were used to maximize retention of the R26^LSL-tdTomato allele even in genomically unstable tumors. Six- to 10-week-old NOD/SCID/γc (NSG) mice (The Jackson Laboratory; stock number: 005557) were used for transplantation experiments. The Stanford Institutional Animal Care and Use Committees approved all animal studies and procedures.

All histologic staining was performed on paraffin-embedded, formalin-fixed sections as described previously (60). Briefly, 4-μm sections were rehydrated and subjected to antigen retrieval before IHC using Vector Lab ABC Vectastain Kit. We used custom FIJI macro scripts for the quantification of IHC. See Supplementary Experimental Procedures for the detail of staining procedures and IHC quantification.

RNA and genomic DNA samples were extracted from 10⁴ to 5 × 10⁴ sorted cancer cells using the Qiagen AllPrep DNA/RNA Micro Kit. RNA from ex vivo FACS-purified cells (15 ng/sample) was converted to cDNA and amplified with the NuGEN Ovation RNA-seq system. Subsequently, amplified cDNA was sonicated and subjected to library preparation using the Illumina TruSeq DNA sample preparation kit. Total RNA from shControl or shBlimp1 688M cells cultured in 0.5% or 20% O₂ was used for the preparation of RNA-seq libraries with Illumina's TruSeq RNA Library Prep Kit v2 according to the manufacturer's protocol. Sequencing was performed on Illumina HiSeq 2000 for 100-bp paired-end (ex vivo samples) and single-end (in vitro samples) reads. See Supplementary Experimental Procedures for details of RNA-seq analysis.

Murine PDAC 688M cells cultured in 0.5% or 20% O₂ were also used for ATAC-seq library preparation. Briefly, nuclei were extracted before incubation with TDE1 Tn5 transposase (Illumina). The fragmented genomic DNA was PCR amplified and ATAC-seq libraries were sequenced on an Illumina NextSeq with paired-end 76 bp reads using an Illumina High Output Kit. ATAC-seq data were processed as previously described with some modifications (81). See Supplementary Experimental Procedures for details of ATAC-seq analysis.

Cell lysates were prepared with RIPA buffer plus protease inhibitors. Proteins were separated by PAGE before being transferred onto a Bio-Rad PVDF membrane. Primary antibodies were incubated in the presence of 5% skim milk at 4°C overnight, followed by staining with horseradish peroxidase–conjugated secondary antibodies. Enhanced chemiluminescence was performed to visualize the proteins of interest. See Supplementary Experimental Procedures for more details of western blot analyses.

To induce hypoxia in vitro, cancer cells were seeded at subconfluency and cultured in a hypoxia chamber (Invivo2-400, Ruskin Technologies) with 0.5% O₂ for 24 hours. Cells were subsequently lysed with TRIzol (Thermo Fisher Scientific, 15596-018) directly on tissue culture dishes for RNA extraction. RNA concentration was quantified on a NanoDrop spectrophotometer (Thermo Fisher Scientific, NanoDrop 2000 UV–Vis Spectrophotometer) and converted to cDNA according to the manufacturer's protocol (Thermo Fisher Scientific, 4368814). For the quantification of transcripts, SYBR green (Sigma-Aldrich, S9194) was used with specific primer pairs. β-actin was used as internal control. See Supplementary Data for more detailed information.

None of the cell lines used in this study were authenticated. The years when the PDAC cell lines were obtained are as follows: murine PDAC cell 688M, 2014; 1004M, 2014; 887M, 2017; 1814, 2015; 1810, 2015; human PDAC cell Panc1, Colo357, BxPC1, AsPC1, and Capan1 were all obtained in 2014. All PDAC cell lines used in experiments in this study were early passage, and aliquots were stored in liquid nitrogen. Thawed cells were used within 1 to 2 months of thawing.

The 688M and 1004M PDAC cells were cultured at subconfluency shortly before harvest for transplantation. All cells used in the transplantation experiments were validated for knockdown efficacies of targeted genes. Briefly, cells were trypsinized and washed 3× in cold PBS before subcutaneous injection. Cells (2.5 × 10⁵ per injection) were injected into the dorsal flank. The numbers of Tom⁺ metastases in the lung were quantified by direct counting using a fluorescence dissecting scope. Alternatively, hematoxylin and eosin sections were used to quantify lung metastases seeded by Tomato-negative 1004M cell line. Metastases in the lung were validated by histology.

The 688M PDAC cell derivatives validated for efficient knockdown were washed 3× in cold PBS before resuspension in 100% Matrigel (Corning, 356231). A surgical procedure was performed with direct injection of the cells/Matrigel mixture into the pancreas of NSG mice. See Supplementary Experimental Procedures for more detail of the orthotopic transplantation.

For comparison between two quantitative variables, we used the Student t test when samples were not paired and the paired t test for paired samples. When more than two variables were compared, either one-way ANOVA or Kruskal–Wallis test were used. For comparison of survival in Kaplan–Meier analyses, we used the log-rank test for univariate survival analyses. The Fisher exact test was used in the analysis of contingency tables. Analyses were performed using Prism 6.0 (Graphpad Software Inc.).

The accession number for all the next-generation sequencing data is included in the following superseries: GSE90825.

M.M. Winslow has received honoraria from the speakers bureaus of Genentech and Merck. No potential conflicts of interest were disclosed by the other authors.

Conception and design: S.-H. Chiou, M.M. Winslow

Development of methodology: S.-H. Chiou, A.S. Kathiria, P. Chu, L. Castellini, A.C. Koong, M.M. Winslow

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): S.-H. Chiou, G.X. Wang, D. Yang, B.M. Grüner, A.S. Kathiria, R.K. Ma, M. Kozak, L. Castellini, A.C. Koong

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): S.-H. Chiou, V.I. Risca, G.X. Wang, B.M. Grüner, D. Vaka, A.C. Koong, M.M. Winslow

Writing, review, and/or revision of the manuscript: S.-H. Chiou, V.I. Risca, B.M. Grüner, R.K. Ma, M. Kozak, E.E. Graves, P. Mourrain, A.C. Koong, A.J. Giaccia, M.M. Winslow

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): S.-H. Chiou, R.K. Ma

Study supervision: M.M. Winslow

Other (reviewed the pathology specimens): G.E. Kim

We thank the Stanford Shared FACS Facility and the Protein and Nucleic Acid Facility for expert assistance; Carolyn Sinow, Santiago Naranjo, and Shashank Cheemalavugu for technical assistance; Chen-Hua Chuang and Nicholas Denko for experimental advice; Teri Longacre for the analysis of human PDAC IHC samples; Justin Kenkel for help with the pancreatic orthotopic transplantation procedure; Stephano Mello, Edward LaGory, and Julia Arand for reagents; Louis Leung, Chris Probert, Peyton Greenside, Andrew Seung-Hyun Koh, and Xun Lan for bioinformatics advice; Laura Attardi, Julien Sage, Jennifer Brady, Kenneth Olive, David Feldser, the Winslow laboratory, the Greenleaf laboratory, and the Stanford pancreatic cancer research community for helpful comments; and Sean Dolan and Alexandra Orantes for administrative support.

This work was supported by a Pancreatic Cancer Action Network–AACR Award in memory of Skip Vinagh (13-20-25-WINS to M.M. Winslow), NIH grant R00CA151968 (to M.M. Winslow), and in part by the Stanford Cancer Institute support grant (P30-CA124435) from the National Cancer Institute. S. Chiou and B.M. Grüner were supported by Stanford Dean's Fellowships. B.M. Gruner was additionally supported by a Pancreatic Cancer Action Network-AACR fellowship in memory of Samuel Stroum (14-40-25-GRUE). V.I. Risca was supported by a Walter V. and Idun Berry Fellowship. E.E. Graves was supported by NIH grant R01CA197136. G.X. Wang and P. Mourrain were supported by the John Merck Fund and NIH grant R01MH099647. P. Mourrain was additionally supported by the Brightfocus Foundation. A.J. Giaccia was supported by NIH grants P01CA067166 and R35CA198291.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.