Clonal ZEB1-Driven Mesenchymal Transition Promotes Targetable Oncologic Antiangiogenic Therapy Resistance Free

On the basis of encouraging clinical trial results (1), bevacizumab, a humanized antibody targeting VEGF, is approved as monotherapy for recurrent glioblastoma (GBM), an aggressive brain malignancy with 90% mortality three years after diagnosis (2). Unfortunately, bevacizumab responses are typically transient, with 50% of GBMs that initially respond progressing soon thereafter (1), with acquired bevacizumab resistance associated with poor outcomes (3). Indeed, phase III trials of bevacizumab in newly diagnosed (4, 5) and recurrent (6) GBM revealed increased progression-free survival (PFS) but unchanged overall survival (OS).

Although some studies have suggested that antiangiogenic therapy resistance involves upregulated compensatory VEGF-independent angiogenesis (7, 8), others have suggested that resistance involves changes enabling tumor cells to thrive in the hypoxic microenvironment of a successfully devascularized tumor (9, 10). Resolving these discordant findings requires comprehensive study of the microenvironment of these resistant tumors. To address this knowledge gap in a manner accounting for the impact of treatment duration that can be prolonged in patients compared with in vivo studies in mice, we analyzed paired specimens from patients before and after tumor progression on variable duration of bevacizumab treatment and in two GBM xenograft models of bevacizumab resistance created by our group (11). These tissues were analyzed for serial treatment-associated changes in hypoxia and vascularity. We also used bulk and single-cell transcriptomics to screen for biomarkers of therapeutic window closure and regulators particularly upregulated in late clones associated with therapeutic resistance (12) whose targeting could prolong the bevacizumab therapeutic window in GBM.

HUVEC cells (ATCC) verified using short tandem repeat (STR) profiling were passaged under 6 times and cultured in EGM-2 (Lonza). Bevacizumab-sensitive and resistant U87-Bev^S/U87-Bev^R and patient-derived (Supplementary Table S1) xenografts were generated and cells for culture extracted as described previously (13, 14), confirmed Mycoplasma-free, and cultured in DMEM/F-12 plus 10% FBS and 1% P/S at 37°C, with some cells treated with 20 μmol/L honokiol (Sigma). For survival studies, trypan blue exclusion viability assay was performed (15). U87-Bev^S, U87-Bev^R and primary GBM neurospheres were cultured in Neurocult/Neurosphere media (StemCell Technologies) with 10 ng/mL bFGF (Thermo Fisher), 20 ng/mL EGF (Thermo Fisher) and B27 and N2 supplement (Thermo Fisher) at 37°C. Accutase (StemCell Technologies) was used to dissociate neurospheres into single cell suspensions.

Animal experiments approved by UCSF Institutional Animal Care and Use Committee (approval #AN105170-02) are in Supplementary Methods.

Morphology and form factor analyses were performed as described previously (11) and in Supplementary Methods.

See Supplementary Methods.

See Supplementary Methods.

ZEB1 was knocked out in U87-Bev^S, U87-Bev^R generations 1,4, and 9 by co-transfecting ZEB1 CRISPR/Cas9 KO and ZEB1 HDR Plasmids (sc-400201 and sc-400201-HDR; Santa Cruz Biotechnology) using FuGENE 6 Transfection Reagent (Promega). Transfections were done per the manufacturer's instructions, and verified using qPCR and microscopy for Red Fluorescent Protein (RFP).

Previously flash-frozen generational xenograft tumor chunks were retrieved and dissociated using a QiaShredder (Qiagen) and passage through a 21-gauge sterile syringe. Dissociated tissue was processed to obtain RNA using the RNeasy kit (Qiagen), following the manufacturer's protocol. RNA quality was tested using an RNA 6000 chip with the Bioanalyzer hardware (Agilent). RIN scores>8 were required for RNA quality. RNA was converted to labeled cRNA using the TargetAmp-Nano Labeling Kit for Illumina Expression BeadChip (EpiCentre), following the manufacturer's protocol. Labeled cRNA was kept at −20°C and given to the UCSF Genome Core Facility (GCF) for chip hybridization.

Tissue was dissociated by incubation in papain with 10% DNAse for 30 minutes. A single-cell suspension was obtained by manual trituration using a glass pipette. Cells were filtered via an ovomucoid gradient to remove debris, pelleted, and resuspended in Neural Basal Media with serum at 1,700 cells/μL, with 10.2 μL cells loaded into each well of a ×10 Chromium Single-Cell capture chip and two lanes captured. Single-cell capture, reverse transcription, cell lysis, and library preparation were performed per the manufacturer's protocol. Sequencing was performed on a HiSeq 2500 (Illumina, 100-bp paired-end protocol). Raw data were processed with CellRangeR 1.3. The resulting count table derived from CellRangeR was processed in R 3.4.1. Data normalization [log (CPM/100+1)] and subsequent t-SNE clustering was performed with the Seurat R package (16). Copy number inference was carried out with the CONICSmat R package (https://github.com/diazlab/CONICS/; ref. 17). Only tumor cells harboring at least one clonal mutation (chr10 loss or chr7 gain, determined by thresholding the posterior probability of the mixture model, pp<0.05) were accounted for in t-SNE clustering. For each clone, we detected genes specific to early (2–3 mutations) or late (>3 mutations) clones (P < 0.05, t test with Benjamini–Hochberg correction and log FC>0.3). We used Opossum (18) to calculate transcription factor binding motif enrichment in the gene sets specific to each clone set (Fisher test). The most significant enrichments were visualized in a bar graph.

GCF analysis of xenograft arrays (.idat files) were processed through the UCSF Bioinformatics Core and deposited in GEO (Accession number = GSE81465). We accessed our archived (19) microarray analysis of BRGs and their paired pretreatment GBMs from ArrayExpress (accession no. = E-MEXP-3296). To identify significantly dysregulated genes across generations, a two-component normal mixture-model was fitted by expectation-maximization to the Z-transformed Log(variance) distribution of gene expression values over generations 4 and 9. A posterior probability cutoff (≥0.95; ref. 20) yielded 717 significantly dysregulated genes. Next, a nonparametric bootstrap procedure (21) was performed using MATLAB to determine whether our data over-represented the Philips mesenchymal, proneural, and proliferative gene sets among these dysregulated genes. Microarray probes were mapped to HGNC gene names, and uniform random sampling of 717 of the 48,324 gene probes followed by matching to the mesenchymal, proneural, or proliferative gene sets was repeated 50,000 times, resulting in a putative null distribution of the number of matches to each gene set. A Poisson probability density function was fit by maximum likelihood estimation to this distribution and a 1-tailed P value calculated for the observed number of matches to the gene set in question. Genes were further clustered by k-means clustering based on gene expression trajectory across generations, and these expression trajectory clusters were manually curated into upregulated genes (sustained increases in expression through generation 4-9), downregulated genes (sustained decreases in expression through generation 4-9), and other genes. Expression of particular genes identified from microarrays was validated using qPCR.

After obtaining RNA in triplicate from xenografts as described in the microarray analysis, cDNA was created using qScript XLT cDNA Supermix (Quanta Bio), following the manufacturer's protocol. cDNA was diluted to a constant concentration for all samples to ensure similar nucleic acid loading. Quantitative PCR was carried out using Power Syber Green Master Mix (Applied Biosystems), primer sequences in Supplementary Tables S2–S4, and an Applied Biosystems StepOne Real-Time PCR cycler following Applied Biosystems Syber guidelines: 95°C for 10 minutes, followed by 40 cycles of 95°C/15 sec and 60°C/1 min. C_t values were calculated using StepOne software (Applied Biosystems).

See Supplementary Methods.

Live-cell imaging was performed using a Nikon Ti-E2000-E2 microscope equipped with a motorized, programmable stage (Applied Scientific Instrumentation), an incubator with constant temperature, humidity, and CO₂ (In vivo Scientific), a digital camera (Photometrics CoolSNAP HG II), and NIS Elements (Nikon) software. Images were taken at 5 ms exposure, 2 × 2 pixel binning using a ×10-objective (Nikon CFI Plan Fluor DLL). Cell motility was measured using phase contrast time-lapse images acquired every 15 minutes over 3 hours. ImageJ software (NIH) was used to track cell centroids from one frame to another to yield migration speeds, which were averaged over the experiment to yield the migration speed of a cell. Cells sticking to each other were excluded.

See Supplementary Methods.

For area analysis in HA devices, cells in devices were imaged every 3 days using Eclipse TE2000 Nikon Microscope with a Plan Fluor Ph1 ×10 objective. Images were acquired and stitched using NIS-Elements Software. For each device, total cell area was outlined in ImageJ and normalized to total day 1 cell area. To analyze morphology, invading edge protrusions were counted and normalized to cell mass length.

See Supplementary Methods.

Pyruvate concentrations, glucose uptake, glycolysis rates, and ATP levels were measured as described (9) and in Supplementary Methods.

See Supplementary Methods.

See Supplementary Methods.

Student t test and Mann–Whitney tests were used to compare means of continuous parametric and nonparametric variables, respectively. Paired t tests (parametric) and Wilcoxon signed-rank tests (nonparametric) were used to compare pre- and post-treatment variables from the same patients. Interobserver variability for manual immunofluorescence counting was assessed using SPSS VARCOMP analysis. Multiple comparisons used ANOVA, followed by Tukey Kramer multiple comparisons test. P values are two-tailed and P < 0.05 was significant.

We investigated the impact of bevacizumab treatment duration on patient bevacizumab-resistant glioblastomas (BRG). Immunohistochemistry for hypoxia marker CA9 and endothelial antigen CD31 revealed that increasing bevacizumab treatment duration before progression increased BRG hypoxia (Fig. 1A; R² = 0.63; P = 0.004) and decreased vessel density (Fig. 1B; R² = 0.66; P = 0.002), without affecting Ki-67 labeling (Fig. 1C; R² = 0.01; P = 0.7). We found similar hypoxic devascularization in two models of bevacizumab resistance established by our group. The first model transfers the effects of prolonged antiangiogenic therapy from the patient to the mouse via patient-derived xenografts (PDX) that maintain the sensitivity or resistance to bevacizumab of the original patient tumor, whereas the second model uses a multi-generational approach to painstakingly recapitulate in mice the prolonged antiangiogenic therapy that preceded resistance in patients (13, 14).

For our xenograft model based on PDXs, intracranial bevacizumab-responsive PDXs (SF8244; Supplementary Table S1) exhibited decreased vessel density (P = 0.02; Fig. 1D) and increased hypoxia (P = 0.01; Fig. 1E) in response to bevacizumab but bevacizumab-resistant PDXs (SF7796 and SF8106; Supplementary Table S1) were devascularized and hypoxic at baseline and those parameters did not change in response to bevacizumab (P = 0.5–0.9; Fig. 1D–E). We also found corroborative xenograft evidence that BRGs became hypoxic and devascularized in our second model of bevacizumab resistance, U87-Bev^R (Supplementary Fig. S1). During its multi-generational ectopic subcutaneous creation, U87-Bev^R exhibited decreased vessel density (P < 0.01; Supplementary Fig. S2) and increased hypoxia (P < 0.001; Supplementary Fig. S3) without altering proliferation (P = 0.07; Supplementary Fig. S4) relative to responsive U87-Bev^S xenografts. The same finding occurred when we implanted the final generation of U87-Bev^R and U87-Bev^S in the orthotopic intracranial microenvironment, where decreased vessel density (P = 0.03; Fig. 1F) and increased hypoxia (P = 0.001; Fig. 1G) were noted with bevacizumab treatment of intracranial bevacizumab-sensitive U87-Bev^S xenografts, whereas intracranial bevacizumab-resistant U87-Bev^R xenografts were devascularized and hypoxic at baseline and those parameters did not change in response to bevacizumab (P = 0.1–0.6; Fig. 1F–G). Thus, patient BRGs and our xenograft models exhibited increased hypoxia with decreased vessel density with increasing bevacizumab treatment duration but maintain their proliferative indices despite this therapy-induced harsh microenvironment, suggesting that BRGs grow despite successful bevacizumab-induced devascularization.

We used our previously published microarray analysis of 15 BRGs relative to paired pretreatment GBMs from the same patients (14, 19) to analyze expression of the 35 signature genes used to define GBM subtypes (proneural, mesenchymal, and proliferative; ref. 22). This analysis revealed a tendency for tumors to become more mesenchymal with increasing bevacizumab treatment duration before progression when bevacizumab was used as monotherapy (n = 5; P = 0.02; Fig. 2A) but not when it was used in combination with traditional chemotherapy (n = 10; P = 0.1; Supplementary Fig. S5).

Similarly, our bevacizumab-resistant PDX SF7796 exhibited increased mesenchymal gene expression by qPCR relative to our bevacizumab-responsive PDX GBM43 (P < 0.05; Fig. 2B). To serially investigate transcriptional changes associated with bevacizumab resistance evolution, we performed microarray expression analysis of U87-Bev^R generations 1, 4, and 9, a model established with the serial use of bevacizumab monotherapy. Bootstrapping analysis of microarray data revealed that 19 of the 170 genes from the extended lists of mesenchymal subtype markers (11%; ref. 22) were significantly dysregulated across U87-Bev^R generations (P < 0.001) and 9 of the 338 genes from the extended list of proneural subtype markers (2%; ref. 22) were dysregulated (P = 0.03). Furthermore, 37% of dysregulated mesenchymal genes were upregulated, whereas 44% of dysregulated proneural genes were downregulated (Fig. 2C and D). These findings were illustrated by a heatmap (Fig. 2E) and a three-dimensional plot created by averaging and normalizing expression of the 14 mesenchymal, 15 proneural, and 5 proliferative genes (22), to obtain mesenchymal, proneural, and proliferative gene expression scores on a +1 to −1 scale (Fig. 2F). The heatmap (Fig. 2E) illustrated early mesenchymal gene upregulation versus early proneural gene downregulation followed by a more modest later rise in proneural gene upregulation, and the three-dimensional plot revealed mesenchymal gene upregulation approaching +1 by generation 4 and persisting thereafter, versus more heterogeneous later changes in proneural gene expression (Fig. 2F). These findings mirrored gene expression changes occurring with increasing bevacizumab treatment duration in patient BRGs. To validate these findings, we performed qPCR for the 14 mesenchymal genes (22) and confirmed increased expression of all but one gene by generation 9 (Fig. 2G).

Because our microarray analysis revealed upregulated genes promoting cytoskeletal and extracellular matrix (ECM) re-organization and tumor cell migration (Supplementary Fig. S6), we investigated BRG cell invasiveness. We found more stellate morphology (P < 0.01; Supplementary Fig. S7; Fig. 3A) and invasiveness in Matrigel assays (P < 0.01; Fig. 3B) in GBM cells from our bevacizumab-resistant PDX models compared with GBM cells from our bevacizumab-responsive PDX models. Similarly, we found more stellate morphology (P < 0.01; Supplementary Fig. S8; Fig. 3C) and invasiveness in Matrigel assays (P = 0.0003–0.04; Supplementary Fig. S9; Fig. 3D) in U87-Bev^R versus U87-Bev^S cells, consistent with prior reports (3, 13, 14, 19). We then investigated BRG cell invasion in two different three-dimensional bioengineered models replicating different GBM invasion modes. The first model was hydrogel platforms (23) that model the gray and white matter along which peritumoral invasion occurs and contain the hyaluronic acid (HA) upregulated in bevacizumab-resistant tumor ECM (24). In this hydrogel model, U87-Bev^R cells proved more invasive than U87-Bev^S cells (P < 0.001; Fig. 3E). The second model was microchannel platforms in which tumor cells start in a tumor-like cell reservoir, which they migrate out of into a 3D matrix of HA conjugated to RGD peptide in a manner that models peritumoral invasion, after which, tumor cells migrate through the HA-RGD toward a parallel open channel “vessel” embedded in the 3D HA-RGD matrix in a manner that models perivascular invasion (25). U87-Bev^R cells were also more invasive than U87-Bev^S cells in these microchannel platforms (P < 0.05; Fig. 3F; Supplementary Figs. S10–S11).

Because residing in the perivascular niche could enable tumor-initiating stem cells to survive therapy-induced hypoxia (26), we investigated whether patient BRGs exhibited enriched tumor-initiating stem cells. Plating cells from patient BRGs (n = 4) versus bevacizumab-naïve recurrent GBMs (n = 3) in neurosphere medium that isolates tumor-initiating stem cells revealed increased cell counts from dissociated neurospheres in BRGs versus bevacizumab-naïve GBMs (P < 0.001; Fig. 4A). Consistent with this finding, our bevacizumab-resistant PDX SF7796 exhibited greater expression of a seven gene stem cell panel (Supplementary Table S2) by qPCR than our bevacizumab-responsive PDX GBM43 (P < 0.05; Fig. 4B). Similarly, U87-Bev^R tumors exhibited increased expression of this seven gene stem cell panel by qPCR, peaking at generation 4 versus U87-Bev^S tumors (P < 0.01 for generation 1; P < 0.001 for generation 4; and P < 0.01 for generation 9; Fig. 4C). Culturing U87-Bev^R cells in neurosphere medium revealed that, although U87-Bev^R cells formed fewer neurospheres than U87-Bev^S cells, U87-Bev^R neurospheres were larger (P = 0.002) and more cellular (P = 0.0002) than U87-Bev^S cells (Fig. 4D–H), reflecting less differentiated, more proliferative stem cells more resistant to hypoxia (27).

To determine whether these changes were arising homogeneously in individual cells or in single cells giving rise to multiple clones, we performed single-cell RNA sequencing of 857 cells from a BRG (Fig. 5A). When focusing on the tumor cells, and moving into the gene space of the 500 most differentially expressed genes identified by our BRG microarray analysis (14), we identified clusters of single cells (Fig. 5A), revealing that BRG cells differ in the expression of these genes. When analyzing these cells for mesenchymal gene expression (Fig. 5B) or copy number variation (Fig. 5C), five clones were identified, a high number indicative of significant evolutionary pressure from bevacizumab. Comparing gene expression of later clones with more mutations to early clones with fewer mutations revealed similar tumor cell proliferative potential based on KI67 expression (P = 0.8 early vs. late clones; Fig. 5D), with a trend toward a proneural to mesenchymal gene expression shift from early to late clones approximating statistical significance (proneural DLL3 downregulated, P = 0.06; mesenchymal YKL-40 upregulated, P = 0.09; Fig. 5D). To identify upstream regulators of these changes, we performed an enrichment test for transcription factor binding sites in the single-cell sequencing data, revealing that ZEB1-binding sites were enriched in genes from early and late clones, with the enrichment P value more significant in late clones (P = 1.0 × 10⁻⁶ vs. 5.1 × 10⁻⁸ for ZEB1 binding site enrichment in early vs. late clones; Fig. 5E), suggesting that ZEB1 could drive BRG gene expression, particularly in late, presumably more resistant clones.

Our single-cell sequencing revealed early clonal evolution of YKL-40 upregulation in patient BRG cells. Consistent with this finding, tumor lysates from our intracranial bevacizumab-resistant SF7796 PDXs exhibited greater YKL-40 gene expression than lysates from our intracranial bevacizumab-responsive GBM43 PDXs (n = 3/group; P < 0.001; Supplementary Fig. S12). Similarly, YKL-40 gene expression arose by qPCR across tumors from multiple U87-Bev^R generations (P < 0.001; Supplementary Fig. S13). Because YKL-40 is a secreted glycoprotein detectable in the circulation, we investigated YKL-40 as a potential biomarker of bevacizumab-induced mesenchymal change. We measured YKL-40 in sera of BRG (n = 8) versus bevacizumab-naïve patients with (n = 11) GBM immediately pre-operatively and found elevated serum YKL-40 in BRG patient serum versus serum from recurrent bevacizumab-naïve GBMs (P < 0.05; Fig. 5F).

We investigated the role of ZEB1, an upregulated gene in our single-cell sequencing expressing a transcription factor associated with the stem cells and mesenchymal change we identified with bevacizumab resistance (28–30), in the resistant phenotype. We found a trend of increased ZEB1 gene expression by qPCR across generations of U87-Bev^R xenografts compared with U87-Bev^S xenografts (P = 1 × 10⁻⁵ −0.02; Fig. 6A). Using CRISPR to disrupt ZEB1 expression in generations 1, 4, and 9 of U87-Bev^R and U87-Bev^S (Supplementary Fig. S14; Fig. 6B) caused complete loss of the mesenchymal gene expression (P < 0.001; Fig. 6C), altered morphology (P < 0.001; Fig. 6D; Supplementary Fig. S15), and parenchymal (P < 0.001; Supplementary Fig. S16) and perivascular (P < 0.001 protrusion density and P < 0.05 invasion area; Fig. 6E) invasiveness seen in U87-Bev^R cells. The role of perivascular invasiveness in bevacizumab resistance has been hypothesized to be due to it serving as a form of VEGF-independent angiogenesis that provides tumor cells with nutrients through direct transfer from tumor cell contact with endothelial cells. Along those lines, CRISPR-mediated ZEB1 knockout in U87-Bev^R cells reduced calcein dye transfer from HUVEC endothelial cells to tumor cells (P = 0.01; Fig. 6F) and reduced expression of tumor cell connexins (P < 0.001; Fig. 6G) that form gap junctions facilitating nutrient transfer between endothelial and tumor cells (31). ZEB1 knockout in U87-Bev^R also reversed the stem cell changes we noted in U87-Bev^R cells, by decreasing their expression of the stem cell gene panel described earlier (P < 0.001; Fig. 6H), stem cell counts (P < 0.001; Fig. 6I), stem cell diameters (P < 0.001; Fig. 6J), and neurosphere yield (P < 0.001; Supplementary Fig. S17).

Another aspect of the mesenchymal changes we noted that is a crucial component of how bevacizumab-resistant tumor cells thrive in the hypoxic devascularized microenvironment we describe here (Fig. 1) is metabolic reprogramming enhanced the Warburg effect. Using qPCR, we found increased expression of the GLUT3 glucose transporter, which we reported to drive metabolic reprogramming during bevacizumab resistance (9), in U87-Bev^R tumors compared with U87-Bev^S tumors (P < 0.001 generation four and P < 0.01 generation nine; Supplementary Fig. S18). Consistent with a lung cancer study (29), targeting ZEB1 with CRISPR in U87-Bev^R shutdown GLUT3 expression assessed by qPCR (P < 0.001; Fig. 7A) and western blot (Fig. 7B). Similarly, pharmacologically targeting ZEB1 with honokiol, a natural phenolic compound from seed cones that inhibits STAT3-mediated ZEB1 transcription (32), in U87-Bev^R cells reversed their elevated ZEB1 (P < 0.05) and GLUT3 (P < 0.05) expression (Supplementary Fig. S19). These findings suggested that ZEB1 was driving the GLUT3 expression that we have shown promotes the metabolic reprogramming of bevacizumab resistance. Consistent with this regulatory role of ZEB1 in metabolic changes associated with bevacizumab resistance, principal component analysis (PCA) of metabolites generated when U87-Bev^R and U87-Bev^S cells expressed CRISPR targeting ZEB1 or control (CTL) sequences were incubated in a low concentration (0.1 g/L) of ¹³C₆-glucose revealed identical metabolite profile between U87-Bev^S expressing CTL versus ZEB1 CRISPR, whereas ZEB1 CRISPR reversed many metabolite changes seen in U87-Bev^R, driving them toward the metabolite profile of U87-Bev^S (Supplementary Fig. S20). ZEB1 CRISPR also decreased glucose uptake (P < 0.001; Fig. 7C; Supplementary Fig. S21), decreased glycolysis (P < 0.001; Fig. 7D; Supplementary Fig. S22), and decreased ATP production (P < 0.001; Fig. 7E; Supplementary Fig. S23) in U87-Bev^R cells. We then used the Seahorse extracellular flux analyzer to dynamically assess oxygen consumption as a measure of oxidative phosphorylation and found that ZEB1 CRISPR increased the oxidative phosphorylation of U87-Bev^R cells (P = 0.002–0.008; Fig. 7F) at baseline (pre-injection 1) and during maximal respiration (post-FCCP) phases of the assay. Similarly, metabolomic comparison of U87-Bev^R CTL CRISPR and U87-Bev^R ZEB1 CRISPR cells in 0.1 g/L ¹³C₆-glucose revealed increased tricarboxylic acid (TCA) cycle intermediate metabolites in ZEB1 CRISPR cells, resembling the metabolic profile of U87-Bev^S cells (Fig. 7G). The translational impact of targeting ZEB1 in BRG was defined when honokiol proved more cytotoxic when treating U87-Bev^R cells than U87-Bev^S cells (P < 0.05; Fig. 7H) and reversed the mesenchymal morphology of surviving U87-Bev^R cells (P < 0.001; Supplementary Fig. S24 and Fig. 7I).

The failure of antiangiogenic agents like bevacizumab to achieve durable response in GBM (4) and other cancers (33) has led to debate about whether acquired bevacizumab resistance involves lost antiangiogenic effect due to compensatory upregulation of VEGF-independent angiogenesis pathways (7) or tumor adaptation to a devascularized microenvironment. We found in patient specimens and two novel xenograft models that BRGs exhibited hypoxia and devascularization that worsened with increasing treatment duration, indicating that the resistance occurs despite continued successful bevacizumab-induced devascularization. In some ways, our findings contrast with the work of Jain and colleagues (34) that demonstrated vascular normalization after antiangiogenic therapy. It is possible that our findings of therapy-induced devascularization in resistant tumors arise because prolonged antiangiogenic therapy beyond the duration of treatment needed to achieve vascular normalization leads to closure of the vascular normalization window. Further study would be needed to confirm this hypothesis, but it would be consistent with reports of greater bevacizumab effectiveness in tumors with greater starting vascularization (35) or milder therapy-induced devascularization (36), suggesting that excessive vascular pruning or rarefaction after bevacizumab therapy may negatively impact patient outcomes.

Our work revealed three adaptive changes in GBM in response to bevacizumab-induced hypoxic devascularization: Perivascular invasion, enrichment of tumor-initiating stem cells, and metabolic reprogramming. Using 3D bioengineered systems customized to reflect the ECM remodeling shown to occur in BRG (24), we demonstrated increased invasiveness, including perivascular invasiveness, in bevacizumab-resistant cells. Invasiveness after bevacizumab resistance has been described as perivascular, which supports tumor cells because the perivascular space is the entry point for nutrients into the brain (37–39). By invading alongside and engulfing preexisting cerebral microvasculature, perivascular invasion co-opts existing vasculature in a VEGF-independent and neo-angiogenesis-independent manner (40), providing a mechanism for continued GBM growth despite the VEGF blockade mediated tumor de-vascularization. Concomitant with this increased invasiveness, we found that the transcriptional changes during bevacizumab resistance enriched tumor-initiating stem cells, progenitor cells resistant to hypoxia both at a cellular level and due to these cells residing in the GBM perivascular niche (26), and described as enriched in other studies of antiangiogenic therapy resistance (41). A metabolic switch to glycolysis and away from oxidative phosphorylation is another adaptation of GBM to the hypoxia exacerbated by bevacizumab, as glycolysis enables more efficient ATP production (9, 42).

We found that mesenchymal change was crucial to these phenotypic changes in xenografts and patient specimens. Our work thus expands upon xenograft studies in which tumor treatment with anti-VEGF antibody promotes increased transcription of mesenchymal factors that portend a worse prognosis in GBM and other cancers (43, 44). GBM molecular subtypes have been defined on the basis of similarity to defined expression signatures. To study subtype changes associated with bevacizumab resistance, we chose the Phillips classification (22) rather than the Verhaak and colleagues (45) classification because the former uses more focused gene sets that change during GBM evolution. Using this classification, we demonstrated mesenchymal change in our xenograft models, which mirrored our finding of increased mesenchymal gene expression with longer bevacizumab treatment duration in patient BRGs. Beyond antiangiogenic therapy, our work adds to studies implicating mesenchymal gene expression in resistance to other anti-cancer therapies such as radiation in GBM (46) and EGFR inhibitors in lung cancer (47). These mesenchymal changes confer context-dependent advantages particular to the therapy that resistance is evolving against. Specifically, radiation resistance arises due to enhanced DNA damage repair in mesenchymal cells, EGFR-targeted therapy resistance arises due to the intrinsic refractoriness of signaling pathways to downregulation in mesenchymal cells, and antiangiogenic therapy resistance arises due to the perivascular invasion, stem cell enrichment, and metabolic changes we found associated with mesenchymal change in GBM. Thus, mesenchymal gene expression not only promotes resistance to antiangiogenic therapy but promotes a phenotype resistant to multiple therapeutic modalities.

To identify a clinically translatable biomarker of when the response window to bevacizumab may close and a translatable target that could prevent resistance, we analyzed bevacizumab resistance genes at a single-cell level. We found that mesenchymal gene expression was more prevalent in late BRG clones, implicating these changes as therapy-related because late clones identified in single-cell sequencing have been associated with therapeutic resistance (12). Further analysis of these late clones revealed that the resistant state had an associated biomarker, elevated serum glycoprotein YKL-40 in patients with BRG, and revealed transcription factor ZEB1 to be a promising regulator of mesenchymal change during bevacizumab resistance evolution. Targeting ZEB1 with CRISPR and honokiol disrupted bevacizumab resistance features and caused preferential toxicity in resistant cells. Although our study is the first to implicate ZEB1 as a bevacizumab resistance driver, serum YKL-40 elevation during bevacizumab treatment of patients with ovarian cancer predicts shorter PFS (48) and low baseline YKL-40 was associated with improved outcomes in patients with ovarian cancer (48) and GBM (49) receiving bevacizumab. Because honokiol could exert ZEB1-independent effects via inhibition of STAT3 (32), our honokiol findings are not as directly implicative of ZEB1's role in bevacizumab resistance as our CRISPR findings. Regardless, because bevacizumab resistance has been challenging to pharmacologically target (50), these findings and the ability to screen inhibitors of bevacizumab resistance in our models or other confirmatory isogeneic bevacizumab resistance models, that could be similarly established, to improve upon our work, which was done in a single isogeneic resistance model before validating in PDXs and patient specimens, could provide meaningful benefit to patients with GBM.

S. Müller is a scientist at Genentech. No potential conflicts of interest were disclosed by the other authors.

Conception and design: A. Chandra, A. Jahangiri, A.T. Nguyen, M.P. Pereira, S. Jain, J.W. Rick, S. Kumar, M.K. Aghi

Development of methodology: A. Chandra, A. Jahangiri, A.T. Nguyen, M.P. Pereira, S. Jain, J.H. Garcia, K.J. Wolf, A.A. Diaz, S. Kumar, M.K. Aghi

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): A. Chandra, A. Jahangiri, A.T. Nguyen, G. Yagnik, M.P. Pereira, S. Jain, J.H. Garcia, S.S. Shah, H. Wadhwa, R.S. Joshi, J. Weiss, K.J. Wolf, J.-M.G. Lin, J.W. Rick, A.A. Diaz, M.K. Aghi

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): A. Chandra, A. Jahangiri, W. Chen, A.T. Nguyen, G. Yagnik, M.P. Pereira, S. Jain, S.S. Shah, H. Wadhwa, K.J. Wolf, J.-M.G. Lin, S. Müller, J.W. Rick, A.A. Diaz, M.K. Aghi

Writing, review, and/or revision of the manuscript: A. Chandra, A. Jahangiri, W. Chen, G. Yagnik, M.P. Pereira, S. Jain, J.H. Garcia, S.S. Shah, K.J. Wolf, J.-M.G. Lin, J.W. Rick, A.A. Diaz, L.A. Gilbert, S. Kumar, M.K. Aghi

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): A. Chandra, A. Jahangiri, A.T. Nguyen, M.K. Aghi

Study supervision: A. Chandra, A. Jahangiri, A.T. Nguyen, S. Kumar, M.K. Aghi

M.K. Aghi was supported by the American Brain Tumor Association, McDonnell Foundation, American Cancer Society, University of California Cancer Research Coordinating Committee, and NIH (1R01CA227136 and 2R01NS079697). A. Chandra, A. Jahangiri, J.W. Rick, and S.S. Shah received Howard Hughes Medical Institute fellowships. A. Chandra received an Alpha Omega Alpha Carolyn L. Kuckein Student Research Fellowship. A. Jahangiri received an NIH Predoctoral Fellowship 1F31CA203372. W. Chen, J.H. Garcia, and J.W. Rick received UCSF School of Medicine Pathways Explore Grants. K.J. Wolf received an NSF Predoctoral Fellowship, NIH funding (F31CA228317), and a Siebel scholarship. S. Kumar received NIH funding (1R21CA174573, 1R01GM122375, and 1R21EB025017). UCLA Metabolomics Center performed metabolomics studies. The Seahorse analyzer was gifted from the S.D. Bechtel Junior Foundation to the Gladstone Institutes.

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.