Abstract
Immunotherapy has limited efficacy in glioblastoma (GBM) due to the blood–brain barrier and the immunosuppressed or “cold” tumor microenvironment (TME) of GBM, which is dominated by immune-inhibitory cells and depleted of CTL and dendritic cells (DC). Here, we report the development and application of a machine learning precision method to identify cell fate determinants (CFD) that specifically reprogram GBM cells into induced antigen-presenting cells with DC-like functions (iDC-APC). In murine GBM models, iDC-APCs acquired DC-like morphology, regulatory gene expression profile, and functions comparable to natural DCs. Among these acquired functions were phagocytosis, direct presentation of endogenous antigens, and cross-presentation of exogenous antigens. The latter endowed the iDC-APCs with the ability to prime naïve CD8+ CTLs, a hallmark DC function critical for antitumor immunity. Intratumor iDC-APCs reduced tumor growth and improved survival only in immunocompetent animals, which coincided with extensive infiltration of CD4+ T cells and activated CD8+ CTLs in the TME. The reactivated TME synergized with an intratumor soluble PD1 decoy immunotherapy and a DC-based GBM vaccine, resulting in robust killing of highly resistant GBM cells by tumor-specific CD8+ CTLs and significantly extended survival. Lastly, we defined a unique CFD combination specifically for the human GBM to iDC-APC conversion of both glioma stem-like cells and non–stem-like cell GBM cells, confirming the clinical utility of a computationally directed, tumor-specific conversion immunotherapy for GBM and potentially other solid tumors.
Introduction
Immunotherapeutic approaches such as immune checkpoint inhibitors have shown high benefits for several types of solid tumors (1, 2). However, their efficacy is limited in glioblastoma (GBM; ref. 3). Development of new treatments is complicated by the immunosuppressed or “cold” tumor microenvironment (TME) of the GBM tumor (3), which is (i) depleted of dendritic cells (DC), the most potent antigen-presenting cells (APC), and tumor-infiltrating CD8+ CTL although dominated by immunosuppressive cells (e.g., regulatory T cells and myeloid-derived suppressor cells; refs. 4, 5) and (ii) characterized by a high degree of intra- and inter-patient molecular heterogeneity (6) and loss of heterozygosity at human leukocyte antigen (HLA) loci and antigenicity, leading to immune escape (7). To “heat up” the cold TME, recent efforts have focused on tumor cell-extrinsic approaches, albeit with mixed results, such as immune checkpoint blockade, re-tuning the TME cytokine milieu, disrupting the blood–brain barrier to recruit immune cells to the TME, or tumor vaccines and cellular therapy that rely heavily on peripheral immune activation (8). Although CTL priming traditionally is thought to occur primarily in tumor-draining lymph nodes, intratumor CTL priming by TME-resident DCs is increasingly recognized as being critical for reversing TME immunosuppression (9). However, circulating DCs have difficulty accessing the GBM TME. Thus, methods that can raise the number of functional intratumor DCs in GBM may be effective at reheating its cold TME.
Recently, Linde and colleagues (10) reported a cancer vaccination approach using a combination of the granulocyte and macrophage master regulator CCAAT/enhancer-binding protein alpha (C/EBPα) with the myeloid differentiation factor PU.1 (SPI1) to induce general myeloid lineage reprogramming in a variety of tumor cell types. However, whereas these myeloid APCs could present endogenous antigens to T cells, they lacked the ability to cross-present exogenous antigens to prime naïve CD8+ CTLs, a functional hallmark associated with DCs and critically required for antitumor immunity (11). In addition, the general myeloid lineage reprogramming may not be an effective approach in GBM and other non-immunogenic tumors, in which myeloid-associated epigenetic programs are frequently adopted by these tumors to escape immune detection (12). Yet, these observations lend credence to the hypothesis that additional cues could be identified to propel the myeloid-associated GBM cells further toward a terminal DC-like state capable of antigen cross-presentation. More recently, Zimmermannova and colleagues (13) applied a three-factor combination of PU.1, IRF8, and BATF3 and successfully converted different types of cancer cells, including GBM cells, to cells with DC-like functions, thereby restoring tumor immunogenicity. Of note, the same combination of transcription factors has previously been identified to reprogram several different normal somatic cells to DC-like cells (14). Therefore, it remains unclear whether optimal and safer conditions for the translational development of this approach will call for a unique set of factors for each tumor type rather than a tissue-agnostic general combination of factors.
A major barrier to fate conversion has been the lack of knowledge about cell fate determinants (CFD), traditionally determined empirically, that is, specific CFDs expressed temporally and/or deterministically to force a desired fate trajectory. Many researchers have instead resorted to the time-consuming process of first converting cells to an induced pluripotent stem cell (iPSC) intermediary, followed by differentiation to the desired cell type using a tightly orchestrated cocktail of tissue-specific growth factors or empirically determined CFDs (15). Fate conversion in cancer is particularly challenging due to the genetic instability and fluctuating fates of the cells within a tumor. Nevertheless, the high degree of fate plasticity of cancer cells may provide a rare opportunity that with the right CFD combination cancer cells could be intentionally reprogrammed to acquire new and potentially therapeutic functions. Here, we developed and then applied a computationally validated machine learning (ML) method to rapidly identify unique CFDs to specifically convert GBM cells into APC-like cells with DC functions (iDC-APC) capable of antigen cross-presentation. We demonstrate that intratumor iDC-APCs effectively reheated the cold TME, reduced tumor growth, prolonged survival, and synergized with existing immunotherapy such as immune checkpoint blockade and DC-based tumor vaccination.
Materials and Methods
Study approval
All animal experiments were approved by the IACUC at the University of Southern California and performed in compliance with all ethical regulations about animal research. For glioma stem-like cells (GSC) and GBM tumor tissues, the study was approved by the IRB of the University of Florida and conducted in accordance with the Declaration of Helsinki ethical standards. Written informed consent was obtained from five patients at first diagnosis for collection and immediate processing of fresh brain tumor samples at the time of surgical resection.
Network-based ML algorithm for CFD prediction
Gene network Reconstruction pipeline (GeneRep)
Inferring regulatory networks from gene expression data is instrumental in discovering new master regulatory genes, especially in the context of cancer, as it can also reveal large-scale disruptions of transcriptional regulation. One of the best characterized and experimentally validated algorithms to reconstruct regulatory gene networks is the mutual information (MI)-based ARACNe (Algorithm for the Reconstruction of Accurate Cellular Networks; ref. 16). In ARACNe, the MI of all edges is computed, and an edge is eliminated if (i) its MI is below a manually preset P value threshold or has the lowest value in any triplet gene combination using the data processing inequality theorem or (ii) its support from a bootstrap procedure is below a threshold. The challenge is how to set an optimal P value threshold for MI filtration—too low and the number of remaining edges is too small, whereas too high and false positive edges dramatically increase. ARACNe’s authors suggest two strategies to determine the P value threshold, based on a null distribution or extrapolation, respectively. The default is extrapolation as it requires less computational power but often leads to high numbers of false positives, whereas null distribution produces the opposite effect with too few edges. To this end, we developed GeneRep, a fully automated pipeline for gene network reconstruction, which uses a three-step filtration process (support, MI, and sum of MI) to produce low rates of false positive edges while maximally recovering true positives by using (i) FDR-based filtration (17) instead of P value; (ii) an adaptive instead of fixed threshold by automatically searching for the best threshold between a real dataset and randomly generated datasets; and (iii) a new sum-of-MI filter. GeneRep is implemented using the NextFlow dataflow management framework, making it suitable for large-scale work in high-performance computing settings.
Performance of GeneRep’s filtration method on six RNA-seq datasets from TCGA
To evaluate the new multi-level filtration approach, we compared its performance to the original ARACNe method using six RNA sequencing (RNA-seq) datasets in TCGA from five cancer types with the largest number of RNA-seq samples—breast carcinoma, renal clear cell carcinoma, lower-grade glioma, head and neck squamous cell carcinoma, and lung adenocarcinoma—as well as GBM (Supplementary Table S1). For each cancer type, we directly compared GeneRep with ARACNe with default (extrapolation) or generated null distribution settings, in which MI thresholds were calculated by extrapolation (ARACNe-extrapolation) or by random reshuffling of input datasets for all genes (ARACNe-null distribution), respectively. For each real dataset, three shuffled datasets were generated with n × (n − 1)/2 null MI pairs, where n represents the number of genes in each dataset (about 20,000), for a total of approximately 1.2 × 109 null MI values, enough to accurately map P value to MI with low P value (10−8).
The MI threshold for ARACNe-extrapolation in each of the six cancers was low (mean of 0.06), and that for ARACNe-null distribution was high (mean of 0.47) even though they are based on the same P value threshold (Supplementary Fig. S1A and S1B). In contrast, GeneRep’s MI thresholds were between extrapolation and null distribution (mean of 0.17). Differences in MI thresholds of the three methods led to broad differences in the numbers of total edges recovered, estimated numbers of false discovered edges, and estimated FDR (Supplementary Tables S1 and S2).
The estimated number of false positive edges is the number of edges in the final network of randomly shuffled datasets that have undergone the same filtration pipeline as the real original dataset. Networks generated by ARACNe-extrapolation with the lowest MI threshold have the highest number of false positive edges as well as FDR (mean of 90%; Supplementary Fig. S1C and S1D; Supplementary Table S1). Due to the inclusion of false edges, the total numbers of ARACNe-=extrapolation networks were also the highest (mean of 757,479; Supplementary Fig. S1E and S1F). In contrast, the ARACNe-null distribution method with the highest MI thresholds produced an FDR close to 0 but at the expense of a significantly reduced number of recovered edges (mean of 90,766). GeneRep, with optimal MI thresholds and filtration steps, generated networks with acceptable FDR (mean of 1%) and significantly higher numbers of recovered edges (mean of 474,299) in comparison with ARACNe-null distribution (Supplementary Tables S1 and S2; Supplementary Fig. S1E and S1F).
Generation of integrated networks for 28 human cancer types
Of the 33 cancer types present in TCGA, we selected 28 cancer types with more than 75 independent RNA-seq samples to generate a reference network for each cancer using GeneRep, integrating major features including RNA-seq, proteomics, microRNA, and clinical data. Contained in each reported reference network are the number of real and shuffled edges, FDR, and the top hubs by number of edges in each cancer-type network (Supplementary Table S3). Several top hub genes include well-known drivers in their respective cancers. For example, RUNX2, SMARCA1, and MYT1L are well described in GBM biology and featured prominently in the top 100 hubs by the number of edges in the GeneRep-generated GBM reference network, and FOXC1, GATA3, and ESR1 in the breast cancer reference network (Supplementary Fig. S2).
network Systems Calculation of Optimal Ranking Engine (nSCORE)
From about 20,000 genes, there are only limited sets of targets that are important for disease development. There are several algorithms to identify those targets based on network theory. Genes in a gene regulatory network can be considered as nodes, and the edge between nodes represents the relationship between genes. Centrality measures used in graph theory and network analysis such as degree, betweenness, Google PageRank, and eigenvalue can be applied to quantify the influence of a node in a network. The simplest degree-centrality measures the number of direct neighbors. In cell networks, -degree centrality identifies hub genes that interact with many other genes and are important in cell biology. Betweenness measures the number of shortest paths from all nodes to all others that pass through the node of interest. Google PageRank and eigenvector are based on the concept that the most important node is the one that connects to the highest number of important nodes, which can be calculated iteratively or algebraically. In iterative calculation, the output score is used as the input for successive calculations until convergence is assumed. A similar approach is adopted in nSCORE. Various centralities can be used as inputs in nSCORE.
Of note, centralities identify hubs in a static, undisturbed network. One approach is to apply centralities to differentially expressed subnetworks extracted from the global network, in which nodes in subnetworks are selected from the top of differentially expressed genes (e.g., by FDR) or if their FDR is less than a threshold of 0.05. Another approach is to measure the level of change in neighboring nodes surrounding the source node, as successfully used in the ranking algorithms CellNet (18, 19) and Mogrify (20) to identify a set of transcription factors that enhances cell fate conversion, and in the VIPER algorithm (21). Some of the limitations of current methods are (i) the exclusive use of networks of known relationship, (ii) only direct targets allowed, (iii) all available node information not fully leveraged, and (iv) iterative scoring not implemented. Iterative scoring allows better capturing of network-wide information. To this end, we developed nSCORE, a generalized automated node importance scoring framework that incorporates limitless scoring schemes using a set of parameters (Supplementary Fig. S3). nSCORE combines many existing parameters known individually to influence network properties, and thus can apply to any type of network and node statistics input. The node importance score (niscore) is the aggregation of the source node and neighborhood scores. The score is calculated iteratively with the output of the previous calculation serving as the input for the next and so on. Inputs include network (e.g., GeneRep and STRING) and node statistics (e.g., logFC, FDR, P value, LR, or centralities; Supplementary Table S4).
Step 1: Extract subnetworks of differentially expressed genes from the whole network using -top_genes_proportion and -g.
Step 2: Calculate individual node scores (indscore): If-k = TRUE, then the input node statistics are converted to rank value. Indscore of gene i: indscorei, whereas stat is k-th statistics in -l of gene i.
Step 3: Calculate the neighborhood score of gene i: nbhscorei = , where wi,j is weight of the closeness between nodes i and j (MI or rho), and ns is the number of neighbors of node i that requires s steps to reach the source node i.
If -a is TRUE, then nbhscorei = .
Step 4: Combine nbhscore and indscore to obtain niscorei.If -d = “n,” then niscorei = nbhscorei. If -d = “m,” then niscorei = indscorei. If -d = “s,” then niscorei = nbhscorei + ngene × indscorei, where ngene is the number of genes in the neighborhood of gene i in case of -a = FALSE to balance indscore with the nbhscore. If -a = TRUE, then ngene = 1. If -d = “p,” then niscorei = nbhscorei × indscorei.
Step 5: Iterative refinement of niScores. The niscores from step 4 are used as input to step 2 for iterative calculation, repeated -r times or until convergence is reached. The sum of gene-level differences in ranking between consecutive iterations will be used as the objective function for monitoring convergence.
nSCORE model parameter selection
To select the parameters for the nSCORE model, we used the cross-training and validation approach based on 48 phenotype pair comparisons curated from public datasets (Supplementary Table S5). The phenotypic comparisons were composed of three main types of gene perturbation: 11 direct cell-type-to-cell-type conversions (trans-differentiation), 24 wild-type gene versus gene depletion or knockout phenotypes, and 13 wild-type gene versus gene overexpression (OE). Ground truths for CFDs in the 11 cell–cell conversions are based on public literature of successful experimental trans-differentiation. For gene knockdown (KD), knockout (KO), or OE comparisons, ground truths are the genes that were deliberately perturbed experimentally. In total 2,882 configurations of parameter sets were generated by combining different parameter options in nSCORE’s input models (See Supplementary Table S6–Excel Format for the entire set of configurations and all the datasets used). To cross-validate, we divided the datasets into training and validation cohorts at a ratio of 5 training to 1 validation and optimized for the highest ranking scores of the ground truth genes. Supplementary Table S6 contains the list of molecular perturbations curated for the nSCORE validation and the median rankings of the ground truth genes of the 2,882 configurations. Supplementary Table S7 shows the robustness of nSCORE’s predictions in 11 gene perturbation datasets as compared with the ground truths. Supplementary Figure S4 showed the mean ranking of ground truth genes in six cross-validation experiments. Configuration 1,796 was demonstrated to be the most robust when all datasets were combined and considered and thus was chosen as the default configuration for all the subsequent applications including those in the GBM to iDC-APC conversions. Configuration 1,796 uses log fold change, P value, and PageRank together to derive node scores. It should be noted that the nSCORE algorithm achieved high accuracy for KD/KO/OE experiments because in these perturbations, the ground truths were firmly established. The trans-differentiation/fate reprogramming’s ground truths have lower predicted rankings than those in the KD/KO/OE comparisons due to the uncertain nature of the ground truths. The empirically determined genes used to convert one cell type to another cell type may be in fact not the most efficient for the conversion and additional functional knowledge of the predicted genes will need to be leveraged to further select optimal combinations for experimental validation.
Generation of lentiviral CFD combination constructs for viral production
ORF expression clones for mouse and human CFDs were purchased from Genecopoeia, and sequence verified. To create the fusion construct, cDNAs encoding either mouse or human IRF8, BATF3, and ID2 were synthesized; linked by P2A and T2A cleavage sequences; and cloned into the lentiviral vector pSF-simple (pSF-simple-mIRF8-P2A-mBATF3-T2A-mID2 for mF3 and pSF-simple-hIRF8-P2A-hBATF3-T2A-hID2 for hF3). For PU.1 constructs, cDNAs encoding mouse and human full-length PU.1 (Spi1) and PU.1 without the proline-glutamate-serine-threonine (PEST) domain [delta-PEST PU.1 (dP), a kind gift from Thomas Graf; ref. 22] were cloned into the lentiviral response plasmid pSF-Lenti, which contains PGK-puro, to generate pSF-Lenti-PU.1 and pSF-lenti-detaPU.1. Human IKZF1, CEBPD, IRF4, CTSZ, and MITF were cloned into the pSF-simple vector individually. For lentiviral production, 7.5 × 106 HEK 293T cells were plated overnight in intact DMEM in a 10-cm Poly-D lysine hydrobromide (SIGMA) coated dish, followed by transfection with a 2:1 ratio of total DNA in PEI (SIGMA) together with the viral packaging and envelop plasmids PSPAX2 and PMD2.G, respectively in advanced DMEM supplemented with 1.25% FBS, 1× pyruvate, 10 mmol/L HEPES and 10 mmol/L sodium butyrate (all from Gibco). Media were replaced at 16 hours after transfection. Viral supernatants were collected every 24 hours and centrifuged for 20 hours at 25,000 g to sediment viral particles, which were resuspended in opti-MEM (Gibco), and viral titers were measured.
Cell culture
Mouse GBM cell lines KR158, KR158-luc, and GL261 (generous gifts from Duane Mitchell in 2018); human GBM cell lines LN428, LN308, and LN827 (generous gifts from Joshua Rubin in 2015; ref. 23); U87 (ATCC, 2020); 4T1 (ATCC, 2016; ref. 24); and HEK 293T cells (ATCC, 2018) were cultured in DMEM media with 10% FBS and 1% pen/strep. GSCs (gift from Brent Reynolds, 2016–2018, previously characterized in refs. 23, 25) derived from GBM samples were grown in stem cell media (STEMCELL, Cat#05750 with basic fibroblast growth factor, EGF, and heparin). DC2.4 (SIGMA, 2020) were cultured in RPMI 1640 supplemented with 10% FBS and 1% pen/strep. Cells were used for experiments within 7 days of thawing or no more than two to three passages dependent on the cell line. Cells were authenticated by the short tandem repeats method on a representative aliquot with the most recent performed in 2022. Mycoplasma screen was performed yearly. For transduction of CFD-containing lentiviruses, GBM cells in baseline DMEM media were incubated with indicated lentiviral constructs for 24 hours, followed by media replacement to RPMI 1640 supplemented with 10% FBS, 1% pen/strep, and 3 µg/μL puromycin for selection of dP- or hP-transduced cells. On day 5 after transduction, cells were passaged into RPMI 1640 media supplemented with 10% FBS and 1% pen/strep without puromycin until analysis.
Quantitative RT-PCR
Total RNA was extracted from KR158, GL261, LN428, LN308, LN827, U87, GSC1, GSC2, and GSC3 cells using the QIAGEN RNeasy Mini Kit (Cat#74106) in accordance with the manufacturer’s protocol. Reverse transcription of 1 µg total RNA was performed using the iScript cDNA Synthesis kit (Bio-Rad, Cat#1708891). qPCR was performed in triplicate using PowerUp SYBR Green Master Mix (Applied Biosystems, Cat#A25741) and Applied Biosystems’ QuantStudio 3 and 5 for indicated genes in relevant figure panels (See Supplementary Table S8 for primer sequences). GAPDH was used for normalization. Gene expression was determined using the ΔΔCt method.
Western blotting
Cells were lysed for 20 minutes on ice with RIPA buffer (150 mmol/L NaCl, 1% NP-40, 0.5% sodium deoxycholate, 0.1% SDS, 25 mmol/L, PH 7.4 Tris) containing protease inhibitors (Roche) and centrifuged at 13,000 g at 4°C for 20 minutes. Supernatants were collected, and protein concentration was determined using a protein assay dye reagent (Bio-Rad). Equal amounts of proteins were resolved by SDS-PAGE and transferred onto polyvinylidene difluoride membranes. Membranes were blocked with 5% nonfat milk in Tris-buffer saline with 0.1% Tween-20 (TBST), incubated with indicated primary antibodies (1:500) at 4°C overnight, washed with TBST, probed with horseradish peroxidase-conjugated secondary antibodies (1:500) at room temperature for 1 hour, and signals detected using a ChemiDoc MP Imaging system (Bio-Rad). For all antibodies used, refer to Supplementary Table S9.
Flow cytometry
Single-cell suspensions were generated by treating adherent cultures of mouse and human adherent GBM cells (see “Cell culture” section above) with 0.25% trypsin (Gibco) for 5 minutes at 37°C and then Fc-blocked for 20 minutes before incubation with designated fluorescent-labeled antibodies (Supplementary Table S9) at 4°C in the dark for 2 hours. Cells were washed once with FACS buffer (2% FBS, DPBS). FACS was performed on BD FACS Canto II and analyzed by FlowJo_V10. Live cells were separated from debris using an SSC-A (y) versus FSC-A (x) dot plot. FSC-A (x) lattice plots, FSC-H (y), and FSC-A (x)/SSC-H (y) versus SSCA (x) lattice plots were used to exclude doublets, as previously described (23). Singlets were then gated and analyzed.
Splenic DC isolation
Splenic DCs were isolated from spleens from healthy C57BL/6 mice (Jackson Lab) using a pan DC isolation kit (Miltenyi Biotec, Cat#130-100-875) following the manufacturer’s instructions. For the single-cell analysis, isolated DCs were subjected to maturation by culturing in KR158-derived tumor lysate for 3 hours and further treated with TNFα (100 ng/mL, Peprotech Cat#315-01A) overnight; 600 immature and 600 mature splenic DCs were mixed together with other FACS-sorted reprogrammed populations for single-cell RNA-seq library construction as described below in the single-cell method. For the cross-presentation experiment, isolated splenic DCs were subjected to maturation by culturing with media containing ovalbumin for 3 hours and further treated with TNFα (100 ng/mL) overnight and used for cross-presentation to OT1 CD8+ T cells as described below.
Antigen presentation assay
Murine GBM cells were transduced with lentiviral CFD constructs on day 1. On day 4, cells were reseeded. For direct presentation, reseeded cells were transduced with lentivirus expressing full-length ovalbumin (OVA; Addgene Plasmid #64599) on day 5. On day 9, 10,000 OVA-expressing reprogrammed cells were cocultured in a 96-well U bottom plate with 100,000 CD8+ T cells isolated from the spleen of OT1 mice (Jackson Lab) using the CD8+ T isolation kit (Cat#130-104-075, 130-104-454, Miltenyl) in T-cell medium (Gibco) supplemented with 50-U/mL mIL2 (Gibco) for 24 hours. Supernatants were then collected for IFNγ ELISA and cells for blocking with anti-mCD16/32 (Cat#553142, Biolegend) for 10 minutes at 4°C before cytometry for CD8, CD3, and CD69 as detailed above, and 15 minutes incubation with a viability dye at 4°C (Cat#65-0865-14, Invitrogen). For cross-presentation, 10,000 day 9 reprogrammed GBM cells were cocultured with 100,000 CD8+ T cells from OT1 mice supplemented with ovalbumin (10 μg/mL, SIGMA) for 24 hours. Supernatants were then collected for IFNγ ELISA and cells for flow cytometry as above.
ELISA for INFγ production by T cells
Supernatants from cocultured media of (i) the direct presentation and cross-presentation of OVA peptide experiments by iDC-APCs (from KR158 and GL261 cells) and OT1 CD8+ T cells and (ii) KR158-specific CTLs and KR158-derived iDC-APCs were analyzed using the DuoSet ELISA Development System per the manufacturer’s instructions (R&D, Cat#DY008). Briefly, plates were coated with a primary antibody (R&D, Cat#DY485-05) at 4°C overnight. At room temperature, samples and standards were incubated in triplicate in the primary antibody-coated plates for 2 hours, followed by 20 minutes of incubation with HRP-conjugated streptavidin and biotinylated antibody for 2 hours, and an HRP chromogenic reagent for 20 minutes, and with three washes in between each step. After the addition of the termination solution, optical density was measured at 450 to 570 nm using a PICO microplate reader (Thermo). Quantification of INFγ concentration in each sample was measured and calculated based on a standard curve.
Phagocytosis assay
In 96-well flat well plates, 1 × 105 reprogrammed cells were plated at 37°C overnight. Media was then carefully removed; cells were washed once with PBS and incubated with 50 μL of prepared pHrodo-red Zymosan bioparticle solution for 1.5 hours at 37°C, followed by fluorescence microscopy using a Zeiss Axios Imager fluorescence microscope (20–40× objectives) or flow cytometry. The pHrodo-red Zymosan bioparticle solution was prepared according to the manufacturer’s instructions (Thermo Fisher, Cat#P35364). Briefly, one vial of Zymosan bioparticles was resuspended in 2 mL live cell imaging solution (Thermo Fisher, Cat#A14291DJ), transferred into a clean 5 mL FACS tube, and sonicated for 5 minutes (20 amplitude).
Immunofluorescence and confocal microscopy
Frozen samples from KR158 brain tumors in C57BL/6J mice transduced with ev-ev, dP-ev, or dP-mF3 set in OCT were sectioned into 12- to 14-μm-thick slides and mounted onto positively charged glass slides, dried in a desiccator for at least 1 hour, permeabilized twice with 0.05% Triton-X in PBS, blocked with a solution of 10% normal goat serum (SIGMA) for 30 minutes with gentle shaking, incubated with primary antibody solution (Supplementary Table S9) in 2% normal goat serum in PBS 0.01% Triton-X in a covered humid chamber in the dark at room temperature overnight. The next day, stained sections were washed three times for 5 minutes with PBS 0.01% Triton-X, incubated with fluorescently labeled secondary antibody (Supplementary Table S9) diluted in 2% normal goat serum PBS 0.01% Triton-X100 (SIGMA) for 1 hour at room temperature, washed three times for 5 minutes with PBS 0.01% Triton X-100, nuclear counterstained with 4'6'-diamidino-2-phenylindole (DAPI, SIGMA) for 10 minutes, mounted an anti-fade medium and allowed to cure for at least 24 hours in the dark until ready for imaging. The stained and fixed samples were imaged using a Zeiss 800 inverted confocal microscope. Images were captured at 63× oil immersion objective either at 0.5× or 1× zoom, keeping all the conditions of the microscope, exposure, and software settings identical for all samples. For all other analyses, Zen software was used.
Intracranial iDC-APC conversion protocol
KR158-luc cells were infected with the indicated CFD lentiviruses for 24 hours prior to surgery and resuspended at a density of 3 × 105 cells per 3-µL PBS. 3 × 105 cells were slowly (1 µL/minute) implanted into the posterior frontal lobe of the brain of 6-week-old C57BL/6J mice or NOD.Cg-PrkdcscidIl2rgtm1Wjl/SzJ (NSG; Jackson Lab) using an automated mouse stereotaxic localizer (Stoelting’s) at 2 mm lateral and 3.5 mm deep on the right side, with the fontanelle as the reference point. Orthotopic tumor growth was monitored by bioluminescence imaging (BLI) and overall survival was recorded [see In vivo Imaging System (IVIS) Spectrum]. For the tumor RNA-pulsed DC vaccine experiment, following the intracranial implantation of 3 × 105 reprogrammed KR158-luc cells into the posterior frontal lobe of the brain of 6-week-old C57BL/6J mice as above, KR158 total RNA pulsed DCs (see DC vaccine for how these were generated) were injected subcutaneously into mice behind the ear at days 1, 8 and 15 at 0.5 × 106/100 μL/mouse. Orthotopic tumor growth was monitored by BLI, and overall survival endpoints were recorded.
Although no differences in tumor penetrance were observed between the sexes in the KR158 model, we have noticed a more delayed latency (by up to 7–10 days) in female mice. Although this difference did not impact overall survival analysis, it may introduce variables that are difficult to control for short-interval readouts such as BLI and TME immunophenotyping at 3 weeks post-intratumor iDC-APCs. For this reason, we used only male mice for all in vivo studies of KR158 that required short-term readouts. For all other experiments, both sexes were used when appropriate.
IVIS pectrum
The IVIS system (Xenogen) was used to monitor brain tumor growth in animals by BLI. Mice were anesthetized with isoflurane (5% induced and 2% maintained). RediJect D-fluorescin bioluminescent substrate (PerkinElmer, Cat#UL08RV01) was injected subcutaneously into mice, and images were repeatedly taken until the bioluminescence signal reached its peak. Data were analyzed using in vivo imaging software (Caliper Life Sciences).
Soluble PD-1 decoy construction and intracranial assay
The soluble PD1 decoy (sPD1) and negative control (NC) constructs were synthesized using the Genescript service. sPD1 was constructed by ligating the synthesized cDNA encoding the codon-optimized extracellular domain protein of mouse wild-type PD1 [amino acids (aa) 1–169] and 6 His tag into the pGenLenti vector (pGenLenti-mouse PD1, Lot#U089YHG070-2/J201884, Genescript). The NC was constructed by ligating the synthesized cDNAs encoding the codon-optimized extracellular domain protein of sPD1 without the programmed death ligand 1 (PDL1) binding motif (aa 1–69 and 78–169) and 6 His tag into the same pGenLenti vector. Culture supernatants of KR158 mouse GBM cells transduced with sPD1 or NC were collected and subjected to Western blotting with mouse anti-His and anti-mouse PD1. For the intracranial assay with iDC-APCs, 3 × 105 KR158-luc cells cotransduced with the reprogramming combination dP-mF3 (dP plus mF3) or the empty vector (ev-ev) control and the sPD1 or NC construct were implanted into the posterior frontal lobe of the brains of 6-week-old C57BL/6J mice on day 0. On day 10, an equal viral titer of either sPD1 or NC lentivirus was injected into the developing tumor through the same needle track using the same coordinates. Orthotopic tumor growth was monitored by BLI, and overall survival endpoints were recorded.
DC vaccine
Tumor RNA-based DC vaccines were electroporated as previously described (26). Briefly, C57BL/6 bone marrow DCs (Cell Biologics, Cat#57-6200) were resuscitated in DC medium [RPMI1640, 5% FBS, 1-mol/L 4-(2-hydroxyethyl)-1-piperazineethanesulfonic acid, 100 mmol/L sodium pyruvate, 55-µmol/L β-mercaptoethanol (Life Technologies), 200 mmol/L glutamine, and 10 mmol/L nonessential amino acids, supplemented with 1% penicillin/streptomycin (Life Technologies), and 10-μg IL4 and 10-μg GM-CSF (R&D Systems) per 500-mL media]. After 24 hours, 25-μg total KR158 tumor RNA was isolated (RNeasy, Qiagen) and then electroporated into five million DCs in 200 μL Opti-MEM by BTX electroporator (ECM 830, BTX), with settings as follows: LV, 300 V, 500 μs, 1 pulse. C57BL/6J mice were vaccinated with a total of 250,000 electroporated DCs at a concentration of 0.5 × 106/100 μL/mouse.
Generation of tumor-specific T cells and their co-culture with reprogrammed GBM cells
Tumor-specific T cells were generated as previously described (26). Briefly, C57BL/6 mice of equal number of males and females received intradermal DC vaccination with 250,000 total KR158 RNA pulsed DCs. Seven days later, the spleens were dissociated into a single-cell suspension. These cells were then re-stimulated with total KR158 RNA electroporated DCs for an additional 5 days to expand tumor-specific T cells. T cells were further isolated using a pan T-cell isolation kit (Miltenyi Biotec, Cat#130095130). The enriched tumor-specific T cells (effectors) were cocultured for 24 hours with day 9 reprogrammed KR158 GBM cells (targets) at E:T ratios of 1:1 and 1:10 in 96-well U bottom plates in T-cell medium supplemented with 50-U/mL mIL2. Supernatants were then collected for IFNγ ELISA and cells for flow cytometry as above. Because parent KR158 GBM cells have severely reduced MHCI and MHCII expression, to confirm that the generated T cells were specific to parent KR158 GBM cells, we pretreated parent KR158 GBM cells with 100-ng/mL IFNγ (R&D Systems, Cat#485-MI-100/CF) for 3 hours to upregulate MHCI and MHCII expression and thoroughly washed off the excess IFNγ prior to coculturing them with the KR158-specific T cells at the E:T ratio of 10:1. Supernatants were then collected for IFNγ ELISA as above.
Cytotoxicity assay
The cytotoxicity assay was performed according to the manufacturer’s protocol (APC Annexin V Apoptosis Detection Kit with 7-AAD, BioLegend, Cat#640930). KR158 cells transduced with ev-ev, dP-ev, or dP-mF3 and cocultured with KR158-specific CTLs were washed twice with cold Cell Staining Buffer (BioLegend), then resuspended in Annexin V Binding Buffer at 0.25 to 1.0 × 107 cells/mL. A total of 100 µL of cell suspension was transferred into a 5-mL test tube containing 5 µL of APC Annexin V and 5 µL of 7-AAD viability dye solution, gently vortexed, incubated for 15 minutes at room temperature in the dark, and 400 µL of Annexin V Binding Buffer added. Fractions of cells with Annexin V binding were identified using BD FACS Canto II flow cytometry with debris and doublets excluded from the analysis.
Single Cell RNA-seq
Single-cell suspensions of KR158 expressing ev-ev, dP-ev, dP-mF3, and spleen-derived mouse DC were prepared and applied to the Chromium Single Cell Chip (10× Genomics) as per the instructions provided by the manufacturer. Only CD45+/MHCII+ sorted cells from dP-ev and dP-mF3 were used. For the ev-ev culture and mouse DCs, whole cell populations were used. Subsequently, single-cell RNA-seq (scRNA-seq) libraries were generated using the Chromium Next GEM Single Cell 3′ (single index). Sequencing output and Q/C are listed in Supplementary Table S10. The main operations were performed using the Seurat R package (3.2.2; refs. 27, 28) unless otherwise stated. When option parameters for the function deviated from the default values, we provided the details of the changes accordingly. Most of the changes to the default options were made to accommodate and leverage the large size of the dataset.
Cell ranger aggregation
Conversion of the raw sequencing data from the bcl to fastq format and the subsequent alignment to the reference genome GRCh38 (GENCODE v.24) and gene count were performed using the cellranger software (10× Genomics, version 4.0.0) with the command cellranger mkfastq, the STAR aligner, and the command cellranger count, respectively. Results from all libraries and batches were pooled together using the command cellranger aggr without normalization for dead cells as it will be handled downstream. The filtered background feature barcode matrix obtained from this step was used as input for sequential analysis.
Normalization of unique molecular identifier
Using the global scaling normalization method, the feature expression for each cell was divided by the total expression, multiplied by the scale factor (10,000), and log-transformed using the Seurat R function NormalizeData with the method “Log Normalize.”
UMAP dimension reduction
The integrated multiple batch dataset was used as input for Uniform Manifold Approximation and Projection for Dimension Reduction (UMAP) dimension reduction (29). The feature expression was scaled using the Seurat function ScaleData, followed by a PCA run using the function RunPCA (Seurat) with the total number of principal components (PC) to compute and store option of 100. The UMAP coordinates for single cells were obtained using the RunUMAP function (Seurat) with the top 75 PCs as input features (dims = 1:75) with minutes; dist = 0.75 and the number of training epochs n.epochs = 2000.
Clustering of cells
We relied on a graph-based clustering approach implemented in the Seurat package, which embeds cells in a K-nearest neighbor graph with edges drawn between similar cells and partition nodes in the network into communities.
Bulk RNA-seq
Total bulk RNAs were extracted from 5 × 106 each of the three biological replicates of LN428, LN308, LN827, and three technical replicates of human monocyte-derived DCs, utilizing QIAGEN RNeasy Midi Kit (Cat#75144) according to the manufacturer’s instructions. Bulk RNA-seq libraries were constructed using an NEBNext Ultra II kit (New England Bio), pooled, and sequenced on a NovaSeq 6000 Illumina instrument. For sequencing analysis, paired-end 150-bp reads were trimmed with trimmomatic v/0.36. Alignment and gene counts were generated against the GRCh38.p12 genome assembly using the annotation GeneCode release 28 by STAR v2.6.0b with default options and quantmode = GeneCounts. Bulk RNA-seq outputs are listed in Supplementary Table S11. Analysis was performed using the ML algorithm GeneRep and nSCORE (see above).
NetMHC 4.0 analysis for epitope identification
The identification of potential neoantigens introduced by four vectors (pSF-Lenti, pSF-mdP, pSF-sim, and pSF-mF3) under the C57BL/6J background was performed using NETMHC 4.0 as previously described (30). Briefly, amino acid sequences from the vectors were formatted into FASTA files. The relevant MHCI alleles for C57BL/6J mice, H-2Kb and H-2Db, were selected. The sequences were uploaded, and the peptide lengths were set to 8 to 11 amino acids, with default binding affinity thresholds.
Statistical analyses
Statistical analysis was performed using GraphPad Prism 8 software. All Student t tests were two-sided, and P values ≤ 0.05 (with 95% confidence intervals) were considered statistically significant for each specific statistical comparison (∗, P < 0.05; ∗∗, P < 0.01; ∗∗∗, P < 0.001; ∗∗∗∗, P < 0.0001). One-way ANOVA was used to adjust for multiple comparisons when appropriate, except for the comparison of TME genome expression in which two-way ANOVA was applied. Data with continuous results are expressed as mean ± SEM.
Data availability
All public datasets used in this study can be found in Supplementary Table S12 for the mouse GBM conversion and Supplementary Table S13 for the human GBM conversion. The raw expression data are available in Gene Expression Omnibus (GEO) with accession numbers GSE270855 and GSE270857. The GeneRep-nSCORE algorithm code developed and used in this study to generate global regulatory gene networks and predict CFDs for the GBM to iDC-APC conversions can be accessed on GitHub at https://github.com/TranLabUSC/NETZEN-classic. All other data are available in the article and its supplementary files or from the corresponding author upon reasonable request.
Results
Identification of CFDs for converting murine GBM cells to induced APCs with a DC-like profile (iDC-APC)
To identify CFDs for the conversion from GBM to iDC-APCs, we applied a network-based ML algorithm to 30 independent RNA-seq profiles of mouse DCs of all subtypes [conventional (cDC1/2) and plasmacytoid (pDC) DCs] from GEO (31) and the three murine high-grade glioma models KR158, CT2A, and GL261 representing several clinicopathologic and genetic features of human GBM with poor (KR158 and CT2A) and moderate (GL261) immunogenicity and responsiveness to immunotherapy (32). To account for the impact of treatment on gene expression and ranking, we also included RNA-seq datasets from the same three GBM lines that had been treated with the standard chemotherapy temozolomide for a total of six independent GBM expression datasets. For comparison, we curated 3 to 25 independent RNA-seq datasets each from other murine immune cell types, hematopoietic stem cells (HSC), and embryonic stem cell/iPSC (Fig. 1A and B; Supplementary Table S12).
Conversion of murine GBM cells to iDC-APCs. A, tSNE 2D expression map of three murine GBM cell lines KR158, GL261, and CT2A, murine DCs, and other immune and non-immune cells, showing a direct versus indirect—via iPSCs and HSC–conversion of GBM to DCs. B, 2D image of a Gene-Rep-nSCORE-generated 1,000-gene regulatory network for the conversion from murine GBM cells to iDC-APCs showing the top 10 ranked CFDs and the four large pathways regulating various aspects of the myeloid and APC states. C, Diagrams of the dP-ev and dP-mF3 lentiviral constructs (right) and radiographs of immunoblotting for expression of indicated CFD combinations in total lysates of KR158 cells. D, Surface expression of CD45, MHCII, CD80, and CD86 in KR158 and GL261 cells at days 7 and 9, respectively, after being transduced with indicated CFD combinations. The gate shown for CD45+MHCII+ is as a percentage of single, live cells. E, Phase contrast photographs of KR158 cells expressing the indicated CFD expression and the positive control mouse DC2.4 cells. Blue arrows denote dendrite-like membrane protrusions. Scale bar, 10 μm. F, Representative histograms of MHCI expression (top) and bar graphs (bottom) of mean fluorescence intensity of MHCI in KR158 and GL261 cells expressing the indicated CFD combinations. G, Bar graphs showing expression of the indicated key components of the antigen processing and presenting machinery and DC-associated cytokines in KR158 (top) and GL261 (bottom) cells expressing the indicated CFD combinations. H, Representative histograms (left) of the cytosolic CFSE dye intensity over 5 days and bar graphs (right) of the percentage of low proliferative cells (maximal CSFE intensity) in subpopulations of dP-mF3-expressing KR158 and GL261 cells. All experiments in triplicate were repeated at least 3 times. Data are represented as mean ± SEM. Analyses were performed using two-way ANOVA. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 00001.
Conversion of murine GBM cells to iDC-APCs. A, tSNE 2D expression map of three murine GBM cell lines KR158, GL261, and CT2A, murine DCs, and other immune and non-immune cells, showing a direct versus indirect—via iPSCs and HSC–conversion of GBM to DCs. B, 2D image of a Gene-Rep-nSCORE-generated 1,000-gene regulatory network for the conversion from murine GBM cells to iDC-APCs showing the top 10 ranked CFDs and the four large pathways regulating various aspects of the myeloid and APC states. C, Diagrams of the dP-ev and dP-mF3 lentiviral constructs (right) and radiographs of immunoblotting for expression of indicated CFD combinations in total lysates of KR158 cells. D, Surface expression of CD45, MHCII, CD80, and CD86 in KR158 and GL261 cells at days 7 and 9, respectively, after being transduced with indicated CFD combinations. The gate shown for CD45+MHCII+ is as a percentage of single, live cells. E, Phase contrast photographs of KR158 cells expressing the indicated CFD expression and the positive control mouse DC2.4 cells. Blue arrows denote dendrite-like membrane protrusions. Scale bar, 10 μm. F, Representative histograms of MHCI expression (top) and bar graphs (bottom) of mean fluorescence intensity of MHCI in KR158 and GL261 cells expressing the indicated CFD combinations. G, Bar graphs showing expression of the indicated key components of the antigen processing and presenting machinery and DC-associated cytokines in KR158 (top) and GL261 (bottom) cells expressing the indicated CFD combinations. H, Representative histograms (left) of the cytosolic CFSE dye intensity over 5 days and bar graphs (right) of the percentage of low proliferative cells (maximal CSFE intensity) in subpopulations of dP-mF3-expressing KR158 and GL261 cells. All experiments in triplicate were repeated at least 3 times. Data are represented as mean ± SEM. Analyses were performed using two-way ANOVA. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 00001.
The top 20 ranked CFDs required to convert each of the three murine GBM cell lines to DCs and other myeloid cell types (i.e., macrophages, microglia, and neutrophils) are listed in Supplementary Table S14. The top 10 ranked CFDs exhibited significant overlaps, especially PU.1, IRF5, IRF8, and ID2—all known to play important roles in myeloid lineage development and functions. ATOX1, a copper chaperone implicated in wound healing (33), was also consistently identified in both GBM-DC and GBM-eosinophil conversions (Supplementary Tables S14 and S15), indicating that it is unlikely an essential driver for the final DC-like state. However, there were also unique differences. For example, BATF3, a critical master regulator of cDC1 (34) and required for TME-resident DCs to regulate CTLs (35), was consistently ranked in the top 12 CFDs across all six GBM-DC conversions, but not present in the top 20 CFDs in the conversions of GBM to other myeloid-cell types. Consistent with the general myeloid lineage reprograming approach (10), conversions from GBM to non-DC myeloid cells featured high ranks for PU.1 and C/EBPΑ or C/EBPΒ, which was absent in the GBM-DC conversions (Supplementary Tables S14 and S15). Among the six independent GBM-DC conversions, the top 10 CFDs were virtually identical, attesting to the accuracy and robustness of the ML-directed method, and can be segregated into three functional groups based on their known functions, including general immune regulation [IRF5 (36)], myeloid [IRF8 (37) and PU.1 (38)], DC [ID2 (39), MYCL (40), and BATF3 (34)], and HSC [CBFA2T3 (41), BCL11A (42), and HHEX (43)] development.
To identify a CFD combination that could efficiently convert murine GBM cells to DCs, we examined the conversion network in detail. Both PU.1 and IRF8 projected direct and parallel connections exclusively to the myeloid differentiation and general immune regulation pathways (Fig. 1B), consistent with PU.1’s known master myeloid regulatory role (38) and cooperativity with IRF8 in regulating MHCII expression in APCs (44), indicating that these two CFDs may be necessary for the GBM-DC conversion. Within the BATF3, ID2, IRF5, and MYCL cluster, we observed significant overlaps in connections to the four regulatory pathways and selected BATF3 and ID2 for further testing because (i) BATF3, a critical DC regulator, was uniquely associated with all six GBM-DC conversions, (ii) ID2 has been reported to cooperate with BATF3 and IRF8 to regulate cDC development (45), and (iii) both IRF5 and MYCL have more general immunoregulatory roles rather than in specific myeloid lineages (46, 47). Lastly, the three CFDs CBFA2T3, BCL11A, and HHEX (41–43) straddled the HSC and myeloid differentiation pathways (Fig. 1B) and thus were predicted to be less important for a direct conversion to DCs without an HSC intermediary. Consequently, we selected the PU.1/IRF8/ID2/BATF3 combination as a candidate for the GBM-DC conversion and compared its ability to induce CD45 and MHCII expression in KR158 and GL261 cells with those that lacked one or more of these four CFDs or included other CFDs from the top 10 list. Combinations lacking PU.1 produced lower frequencies of CD45+ and no double CD45+MHCII+ cells, confirming its critical requirement for the conversion (select combinations shown in Supplementary Fig. S5A–S5D). Given that wild-type murine PU.1 is unstable due to the PEST-rich degradation domain, we instead used a PEST-negative PU.1 (PU.1 delta PEST or dP; ref. 22) and confirmed that dP was much more efficient at inducing CD45+ cell and upregulating MHCI, a known target of PU.1 (Supplementary Fig. S5B). As we predicted, any combination containing dP but lacking BATF3, ID2, or IRF8—including the dP/IRF8/BATF3 (PIB) combination reported to reprogram many somatic and cancer cells to DC-like cells (13, 14)—reduced fractions of CD45+ and CD45+MHCII+ cells relative to PU.1/IRF8/ID2/BATF3 (Supplementary Fig. S5), reaffirming the cooperativity of all four CFDs and the ML-directed, tumor-specific CFD identification algorithm. Accordingly, the PU.1/IRF8/ID2/BATF3 combination was advanced to functional validation. In subsequent experiments, we linked dP to a puromycin resistance cassette via an internal ribosomal entry site constructed into the pSF-Lenti vector for the enrichment of dP-expressing cells with puromycin. The other three CFDs (IRF8/ID2/BATF3) were synthesized as a fusion construct using viral peptide linkers (designated mF3) into the pSF-core vector without the puromycin cassette. This approach ensured all four CFDs were efficiently expressed (Fig. 1C), as combining them and the puromycin cassette into the same vector resulted in an excessive payload and suboptimal transgene expression. In addition, to rule out a potential new epitope imbalance in the four vectors (i.e., the empty pSF-Lenti and pSF-core, dP, and mF3) that could confound our subsequent analyses, we performed an epitope prediction inquiry of these vectors using the NetMHC 4.0 algorithm (30). No high binders that could be presented by MHCI H2kb and H2db of the KR158’s and GL261’s C57BL/6J background were identified in any of the vectors, whereas only one weak binder in the pSF backbone was shared in all four vectors (Supplementary Table S16).
In both KR158 and GL261 cells, dP plus mF3 (dP-mF3) expression produced similar to higher fractions of CD45+ cells, compared with dP plus the empty virus pSF-core (dP-ev). Of the CD45+ cells, dP-mF3 generated up to 2–3-fold higher fractions of CD45+MHCII+ cells compared with dP-ev (Fig. 1D). Expression of the co-stimulatory receptor ligands CD80 and CD86 on CD45+MHCII+ cells was also up to 4-fold higher in cells expressing dP-mF3 compared with dP-ev. In addition, expression of the cDC1-associated markers CD11c, XCR1, and CLEC9A (48) was 2–8-fold higher in dP-mF3 compared with dP-ev. The empty virus pSF-Lenti plus empty pSF-core combination (ev-ev) produced negligible to no CD45+ and CD45+MHCII+ cells with cDC1 markers (Fig. 1D; Supplementary Fig. S6). CD45+MHCII+ cells generated using dP-mF3, when compared with those generated using dP-ev, displayed morphologic changes from adherent, elongated cells to non-adherent, round cells with dendrite-like protrusions, reminiscent of natural DCs (Fig. 1E; Supplementary Fig. S7). Compared with the ev-ev control, robust upregulation of MHCI and the immune checkpoint receptor PDL1 was also observed with both dP-mF3 and dP-ev, reflecting the increased inflammation associated with reprogramming, with dP-mF3 again exhibiting a significantly higher degree (Fig. 1F; Supplementary Fig. S8). To determine whether dP-mF3-induced CD45+MHCII+ cells were propelled further toward an APC state compared with those expressing dP-ev, we measured mRNA expression of the canonical antigen processing and presenting machinery for both MHCI and MHCII [e.g., transporter associated with antigen presentation (TAP)1, TAP2, immunoproteasomes low-molecular-mass protein-2 and LMP7, endoplasmic reticulum aminopeptidase 1, and TAP-associated glycoprotein; ref. 49] and representative cytokines [e.g., interleukin-12beta (IL12β), IL6, IL15, and TNFα; ref. 50] frequently associated with APCs, especially DCs. In both conversions in KR158 and GL261 cells, dP-mF3-induced CD45+MHCII+ cells consistently unregulated several antigen processing components and DC-associated cytokines to 2−10-fold higher levels compared with those with dP-ev, suggesting that they may be closer in phenotype to APCs, especially DCs (Fig. 1G). The ev-ev control expressed negligible levels of these markers. In addition, dP-mF3-induced CD45+MHCII+ cells had the highest fractions (40% in KR158 and 80% in GL261) of low-proliferative cells as measured by retention of the intracellular fluorescent cell tracker CFSE over a course of 5 days, compared with 20% to 33% and 20% to 46% in partially (CD45+MHCII− and CD45−MHCII+), 4% and 8% in non-reprogrammed (CD45−MHCII−), and 1% and 10% in parent cells, respectively (Fig. 1H), consistent with the dP-mF3-induced CD45+MHCII+ cells approaching closer to the terminally differentiated APC state.
Next, to scrutinize the reprogramming dynamics as GBM cells transition toward the APC/DC state, we performed scRNA-seq analysis of sorted CD45+MHCII+ cells from dP-mF3 and dP-ev cultures and unsorted ev-ev controls in KR158 cells as compared with murine splenic cDCs/APCs from three independent healthy C57BL/6J mice. For batch effect control and normalization across different sample types, we added 5% of a standardized murine breast cancer cell line 4T1 (24) from the BALB/C background. Sequencing output is shown in Supplementary Table S5. In total, approximately 28,537 single cells were resolved using a graph-based clustering technique in the Seurat R package (28) and UMAP (29) for dimension reduction (Fig. 2A). CD45+MHCII+ cells from both dP-mF3 and dP-ev formed a distinct cluster from ev-ev control cells and lay between the ev-ev control and splenic cDCs. The two clusters of CD45+MHCII+ cells from dP-mF3 and dP-ev had significant overlaps in UMAP positions and markers found on general APCs (e.g., Fcer1g, Cd52, Cd74, Cxcl16, Ccl5, Asprv1, HLA class II—H2-aa, H2ab1, and H2-eb1, HLA class I—H2-k1 and H2-q7), predominantly on DCs (Cd11c), or on macrophages (Cybb, Il4i1, Cd11b, Clec4d, Ccl22, and Ccr7), expression of markers more specific to DCs, specifically cDC1, like Clec9a and Xcr1 (48) was observed only in dP-mF3-induced CD45+MHCII+ (Fig. 2A; Supplementary Figs. S9 and S10). In addition, KR158 cells reprogrammed by either dP-mF3 or dP-ev, lacked B- and T-cell markers regardless of the extent of reprogramming (Supplementary Fig. S11), ruling out lymphocytes, especially lymphocytes in a B-cell state as an alternate fate. When compared with cDCs, dP-mF3-induced CD45+MHCII+ cells shared broad similarity in global expression landscape that was lacking in dP-ev-induced CD45+MHCII+ cells and ev-ev controls (Fig. 2B), especially in several key DC functional pathways defined in Gene Ontology (51)—for example, cDC differentiation (GO:0043011), antigen processing and presenting via MHCI (GO:0002474), and DC-associated cytokine regulation (GO:0002732; Fig. 2C).
iDC-APCs share expression profile overlaps with murine cDCs. A, scRNA-seq UMAP cluster maps at resolution 1 of sorted CD45+MHCII+ cells from KR158 GBM cells expressing dP-mF3 or dP-ev as compared with unsorted single-cell ev-ev KR158 cells and enriched splenic cDCs in two clusters with overlapping markers for DCs and macrophages. 4T1 cells were used for batch effect normalization. Expression of indicated mRNAs of DC- and macrophage-specific genes is shown in both composite and split views of individual clusters. B, An expression heatmap of a 200 gene cluster in CD45+MHCII+ cells from dP-mF3 and dP-ev as compared with ev-ev controls and enriched cDCs. C, Expression heatmaps of three functional pathways in the DC state in CD45+MHCII+ cells from dP-mF3 and dP-ev as compared with ev-ev controls and enriched DCs.
iDC-APCs share expression profile overlaps with murine cDCs. A, scRNA-seq UMAP cluster maps at resolution 1 of sorted CD45+MHCII+ cells from KR158 GBM cells expressing dP-mF3 or dP-ev as compared with unsorted single-cell ev-ev KR158 cells and enriched splenic cDCs in two clusters with overlapping markers for DCs and macrophages. 4T1 cells were used for batch effect normalization. Expression of indicated mRNAs of DC- and macrophage-specific genes is shown in both composite and split views of individual clusters. B, An expression heatmap of a 200 gene cluster in CD45+MHCII+ cells from dP-mF3 and dP-ev as compared with ev-ev controls and enriched cDCs. C, Expression heatmaps of three functional pathways in the DC state in CD45+MHCII+ cells from dP-mF3 and dP-ev as compared with ev-ev controls and enriched DCs.
Overall, these results indicate that the PU.1/IRF8/ID2/BATF3 combination is sufficient and more efficient than PU.1 alone at reprograming murine GBM cells toward iDC-APCs, with the resultant cells showing DC-associated global phenotypes and growth arrest commensurate with the degree of reprogramming characteristics of professional APCs.
iDC-APCs induced with dP-mF3 exhibit DC-like functions in vitro
Next, we sought to validate the antigen processing and presenting functions of iDC-APCs by measuring three key properties of professional APCs, including phagocytosis, direct presentation of endogenous antigens, and cross-presentation of exogenous antigens to naïve CD8+ T cells, the latter being a hallmark function of DCs critical for antitumor immunity (11). To evaluate phagocytic activity, we measured the ability of iDC-APCs to engulf fungal particles loaded with the pHrodo red dye that turns red in the lysosomal acidic environment once phagocytosed (52). CD45+MHCII+ cells produced with either dP-mF3 or dP-ev exhibited comparably high phagocytotic capability (Fig. 3A), much higher than in partially reprogrammed (CD45+MHCII− or CD45−MHCII+) and non-reprogrammed (CD45−MHCII−) cells, consistent with PU.1 being the master myeloid regulator with its associated phagocytotic upregulation (Supplementary Fig. S12). The ev-ev controls possessed negligible phagocytosis capacity. To assess iDC-APCs’ antigen processing and presenting capacity, we turned to the well-established chicken OVA antigen assay, in which the SIINFEKL epitope is presented in the context of the MHCI H2kb on APCs and recognized by the transgenic OT1 T-cell receptor CD8+ T cells (53). To facilitate direct presentation, we transfected a full-length OVA construct into reprogrammed GBM cells to bypass antigen internalization, followed 3 days later by measurement of the surface level of H2kb-SIINFEKL using an antibody (Fig. 3B). Consistent with the relative upregulation of MHCI under these conditions (Fig. 1F), the fraction of H2kb-SIINFEKL+ cells in the dP-mF3 culture was 3-fold higher than in dP-ev, which in turn was more than 15-fold higher compared with the ev-ev control (Fig. 3C). When only the CD45+ fraction was examined, 25% of dP-mF3-induced cells were positive for H2kb-SIINFEKL. In contrast, only a negligible fraction of the dP-ev-induced cells were (Fig. 3C), even though nearly all KR158 cells greatly upregulated MHCI with either dP-ev or dP-mF3 (Fig. 1F), likely reflecting in part the much lower expression of the antigen processing and presenting machinery in dP-ev cells compared with dP-mF3 (Fig. 1G). Moreover, the higher direct presentation of intracellular SIINFEKL in dP-mF3 cells translated to more efficient priming and activating of naïve OT1 CD8+ T cells compared with dP-ev and ev-ev control cells, as measured by the early T-cell activation marker CD69 (54) and IFNγ secretion (Fig. 3D).
Validation of DC-like properties of murine iDC-APCs. A, Representative fluorescence images and bar graphs of pHrodo-red particle phagocytosis in CD45+MHCII+ KR158 and GL261 cells expressing the indicated CFD combinations. DC2.4 cells serve as positive controls. Scale bar, 10 μm. B, An assay diagram of direct presentation of OVA to OT1 CD8+ CTLs. C, Bar graphs of H2kb-SIINFEKL expression in KR158 cells expressing the indicated CFD combinations as a percentage of live (left) or CD45+ cells (right). D, Representative dot plots (left) and bar graphs of OT1 CD8+ T-cell activation (right) as measured by CD69 and IFNγ expression through direct presentation of endogenous OVA by KR158 cells expressing the indicated CFD combinations. E, An assay diagram of cross-presentation of exogenous OVA to OT1 CD8+ CTLs. F, Representative dot plots (left) and bar graphs of OT1 CD8+ T-cell activation (right) as measured by CD69 and IFNγ expression through cross-presentation of exogenous OVA by KR158 and GL261 cells expressing the indicated CFD combinations and by spDCs with or without premixing with ev-ev controls at the same ratio of CD45+MHCII+ as in dP-mF3. All experiments in triplicate were repeated at least three times. Data are represented as mean ± SEM. Analyses were performed using two-way ANOVA. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 00001.
Validation of DC-like properties of murine iDC-APCs. A, Representative fluorescence images and bar graphs of pHrodo-red particle phagocytosis in CD45+MHCII+ KR158 and GL261 cells expressing the indicated CFD combinations. DC2.4 cells serve as positive controls. Scale bar, 10 μm. B, An assay diagram of direct presentation of OVA to OT1 CD8+ CTLs. C, Bar graphs of H2kb-SIINFEKL expression in KR158 cells expressing the indicated CFD combinations as a percentage of live (left) or CD45+ cells (right). D, Representative dot plots (left) and bar graphs of OT1 CD8+ T-cell activation (right) as measured by CD69 and IFNγ expression through direct presentation of endogenous OVA by KR158 cells expressing the indicated CFD combinations. E, An assay diagram of cross-presentation of exogenous OVA to OT1 CD8+ CTLs. F, Representative dot plots (left) and bar graphs of OT1 CD8+ T-cell activation (right) as measured by CD69 and IFNγ expression through cross-presentation of exogenous OVA by KR158 and GL261 cells expressing the indicated CFD combinations and by spDCs with or without premixing with ev-ev controls at the same ratio of CD45+MHCII+ as in dP-mF3. All experiments in triplicate were repeated at least three times. Data are represented as mean ± SEM. Analyses were performed using two-way ANOVA. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 00001.
To assess iDC-APCs’ cross-presentation capability, we added wild-type, full-length OVA protein to a coculture of iDC-APCs and purified naïve OT1 CD8+ T cells and measured H2kb-SIINFEK-specific activation of OT1 CD8+ T cells (Fig. 3E). In both KR158 and GL261 cells, dP-mF3-induced iDC-APCs were more efficient at OVA cross-presentation and activation of OT1 CD8+ T cells than those with dP-ev and the ev-ev controls (Fig. 3F), indicating that dP-mF3-induced iDC-APCs are the most DC-like. To further appraise the cross-presenting efficiency of dP-mF3-induced iDC-APCs, we obtained isolated syngeneic splenic DCs (spDC), confirmed their surface marker identity and CD80/CD86 co-receptor upregulation following electroporation with the eGFP antigen (Supplementary Fig. S13), and compared dP-mF3-induced iDC-APCs to these natural counterparts. For proper comparison, we premixed the same number of splenic DCs with ev-ev control cells to match the fraction of CD45+MHCII+ in the dP-mF3 condition and then added the mixture to the coculture with OT1 CD8+ T cells. The OVA cross-presentation capacity of dP-mF3-induced iDC-APCs was similar to superior compared with ev-ev+spDCs and reached 60% to 70% efficiency of spDCs added alone without the ev-ev GBM reconstitution (Fig. 3F).
Taken together, the PU.1/IRF8/ID2/BATF3 combination reprograms murine GBM cells to iDC-APCs with canonical DC-like phenotypes and functions that approach those of natural DCs.
iDC-APCs reheat the cold TME of GBM
To determine the immune impact of iDC-APCs and rule out other non-immunologic effects of the CFDs on tumor growth, we orthotopically implanted an equal number of luciferase-tagged KR158 (KR158-luc) cells expressing dP-mF3, dP-ev, or ev-ev into the right posterior front lobe of immunodeficient NSG and syngeneic immunocompetent C57BL/6J mice, and measured tumor growth rates by luciferase BLI, survival, and immunologic changes in the TME (Fig. 4A). The KR158 model was selected for the extreme coldness of its TME and lack of responsiveness to immunotherapy (23). In the absence of an adaptive immune system in the NSG hosts, tumor growth rates were comparable in all three populations at day 7 post-implantation. By day 14, whereas no difference was observed between dP-mF3 and dP-ev, there was a statistically significant and comparable decrease in tumor size in both dP-ev and dP-mF3 animals compared with the ev-ev control (Fig. 4B), which may reflect the impact of the two CFD iterations on tumor cells’ proliferation. However, this effect was negligible because there was no difference in median overall survival in the three cohorts, with no survivor past 30 days (Fig. 4D). In contrast, in the immunocompetent C57BL/6J hosts, dP-mF3-expressing tumors grew at a consistently slower rate compared with the dP-ev and ev-ev control tumors (Fig. 4C). This translated into a statistically significant survival advantage favoring dP-mF3 over both dP-ev and ev-ev animals, whereas no difference was observed between dP-ev and ev-ev animals (Fig. 4E). These results indicate that the major in vivo effects of dP-mF3-induced iDC-APCs were immune-mediated.
Intratumor iDC-APCs reactivate the cold TME of murine GBM tumors. A, A schema detailing the intratumor iDC-APC reprograming in immunosuppressed NSG and immunocompetent C57BL/6J hosts. B–E, Antitumor effects of intratumor iDC-APCs in KR158-luc cells is dependent on adaptive immunity: Representative photographs and dot plots showing orthotopic growth by BLI of GBM expressing the indicated CFD combinations in NSG (B) and C57BL/6J (C) hosts; n = 10 per group; Kaplan–Meier estimates showing survival rates after implantation in NSG (D) and C57BL/6J (E) hosts. Two-way ANOVA was used to compare tumor size differences and a log-rank test to compare survival rates. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ns, not significant. F–L, iDC-APCs reheat the cold TME of GBM: The immune TME of the indicated reprogrammed tumors are summarized in an expression heatmap of a 24-immune-gene profile by qRT-PCR (n = 5 per cohort; F) and representative images of immunofluorescence and bar graphs of the mean numbers (yellow arrows) per 20× field of CD11c+CLEC9A+ cDC1 (G), F4/80+ macrophages (H), CD3+CD4+ T cells (I), CD3+CD8+ CTLs (J), GZMB+CD8+ CTLs (K), and CD25+CD8+ CTLs (L) with 4,6-diamidino-2-phenylindole (DAPI) counterstain. Three to five independent slides were used for enumeration. Scale bar, 50 μm. Data are represented as mean ± SEM. Two-way ANOVA was used to compare the immune TME profile and IF image difference. ****, P < 00001.
Intratumor iDC-APCs reactivate the cold TME of murine GBM tumors. A, A schema detailing the intratumor iDC-APC reprograming in immunosuppressed NSG and immunocompetent C57BL/6J hosts. B–E, Antitumor effects of intratumor iDC-APCs in KR158-luc cells is dependent on adaptive immunity: Representative photographs and dot plots showing orthotopic growth by BLI of GBM expressing the indicated CFD combinations in NSG (B) and C57BL/6J (C) hosts; n = 10 per group; Kaplan–Meier estimates showing survival rates after implantation in NSG (D) and C57BL/6J (E) hosts. Two-way ANOVA was used to compare tumor size differences and a log-rank test to compare survival rates. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ns, not significant. F–L, iDC-APCs reheat the cold TME of GBM: The immune TME of the indicated reprogrammed tumors are summarized in an expression heatmap of a 24-immune-gene profile by qRT-PCR (n = 5 per cohort; F) and representative images of immunofluorescence and bar graphs of the mean numbers (yellow arrows) per 20× field of CD11c+CLEC9A+ cDC1 (G), F4/80+ macrophages (H), CD3+CD4+ T cells (I), CD3+CD8+ CTLs (J), GZMB+CD8+ CTLs (K), and CD25+CD8+ CTLs (L) with 4,6-diamidino-2-phenylindole (DAPI) counterstain. Three to five independent slides were used for enumeration. Scale bar, 50 μm. Data are represented as mean ± SEM. Two-way ANOVA was used to compare the immune TME profile and IF image difference. ****, P < 00001.
To determine the molecular basis of the positive clinical impact of dP-mF3-induced iDC-APCs in the C57BL/6J hosts, we profiled the immune TME in these animals at 3 weeks after implantation for transcripts of 24 key markers of DCs, TILs, and immune checkpoint receptors. Markers for cDCs and pDCs (23) were greatly upregulated in dP-mF3 tumors compared with the cold TME of the ev-ev tumors (Fig. 4F). The relative increases in cDCs and other APCs such as macrophages in dP-mF3 tumors were confirmed by immunostaining for CD11c+CLEC9A+ cells (Fig. 4G) and F4/80+ cells (Fig. 4H), respectively. Moreover, there was a significant surge in markers for TILs and specifically CD8+ CTL activity and chemokines [INFγ, granzyme B (GZMB), perforin 1 (PRF1), CX3CR1 and CCL4; Fig. 4F; ref. 23], which was recapitulated by immunostaining, which showed a large influx into the TME specifically in dP-mF3 tumors of CD3+CD4+ T cells (Fig. 4I) and CD3+CD8+ T cells (Fig. 4J), with pervasive expression of early (CD25) and mature (GZMB) effector activation markers (Fig. 4K and L)–both have been particularly associated with antigen-specific CD8+ T-cell stimulation (55, 56). Given that the conversion program did not produce lymphocytes (Supplementary Fig. S11), the increase in T-cell recruitment to the TME was likely a direct result of the iDC-APC-induced reheating of the TME. Several immune checkpoint receptors, which are early phase markers of immune activation, were also upregulated in dP-mF3 compared with ev-ev control tumors. For all these markers and immune cell TME infiltrates, tumors expressing dP-ev exhibited modest to negligible increases compared with ev-ev and dP-mF3 tumors with no detectable recruitment of activated CD25+ and GZMB+CD8+ CTLs (Fig. 4F and G), which may explain the lack of survival benefit in these animals.
Taken together, these observations confirm and paint a coherent in vivo picture, in which dP-mF3-induced iDC-APCs with DC-like properties created directly in the TME were sufficient to reheat the cold TME of GBM, reduce tumor growth, and improve survival. However, we cannot rule out contributions of iDC-APC–independent effects of the dP-mF3 to the observed phenotypes.
Synergy between dP-mF3-reprogrammed GBM and an intratumor sPD-1
The significant upregulation of several immune checkpoints during the dP-mF3-induced reprogramming process and in dP-mF3 tumors, especially the PD1/PDL1 axis (Supplementary Fig. S8; Fig. 4F), raises the possibility that they may restrict the antitumor functions of tumor-infiltrating activated CD8+ CTLs (Fig. 4J–L). Although immune checkpoint blockade has shown high benefits for other solid tumors (1, 2), its efficacy in GBM is limited, however, due in part to the impeded access of these agents to the cold TME of GBM (3, 57). As a result, we asked whether an sPD1 expressed directly in the TME could synergize with dP-mF3-induced iDC-APCs and TME reprogramming to improve tumor control and survival. To generate the sPD1, we removed the transmembrane domain (TM) within exon 3 and the cytoplasmic domain from wild-type murine PD1, as previously described (58). For a soluble NC, we also deleted the 7-aminoacid PDL1-binding motif within the immunoglobulin V sequence in exon 2 of sPD1 (Fig. 5A). Both the sPD1 and NC constructs contain a C-terminal His tag. Expression of sPD1 and NC in supernatants of KR158 GBM cells cotransduced with a lentivirus encoding dP-mF3 and either sPD1 or NC was confirmed with anti-His, whereas only sPD1 was detected using the neutralizing anti-murine PD1 RMP14 (59), which is specific for the PDL1-binding motif (Fig. 5B). In the KR158 GBM model, which is known to be highly resistant to anti-PD1 immunotherapy (60), the addition of sPD1, but not NC, significantly extended both orthotopic tumor control (Fig. 5C–E) and overall survival (Fig. 5F) in the dP-mF3 animals (dP-mF3 + NC vs. dP-mF3 + sPD1), whereas neither sPD1 nor NC alone conferred any clinical benefit to the ev-ev controls. These findings support a strong therapeutic synergy between dP-mF3-reprogrammed GBM TME and anti-PD1 immunotherapy.
GBM reprogramming synergizes with an intratumor sPD1. A, A schematic of the sPD1 and the NC constructs. The NC construct lacks the PDL1-binding motif in the immunoglobulin V (IgV) sequence. B, Representative immunoblots of his-tagged sPD1 and NC proteins in supernatants of KR158 GBM cells co-transfected with dP-mF3 and sPD1 or NC. Only sPD1 is detectable by the anti-PD1 mAb RMPI14 specific for the PDL1-binding motif absent in NC. C–F, A schema detailing the combination of dP-mF3-reprogrammed KR158 GBM and intratumor sPD1 (C), specifically creating a therapeutic synergy that delayed tumor growth as measured by BLI (D) and quantified in a line graph (E), and significantly extended survival as measured by Kaplan–Meier estimates (F). N = 10 per group. Log-rank test was used to compare survival rates. *, P < 0.05; **, P < 0.01; ****, P < 0.0001; ns, not significant. wtPD1, wild-type PD1.
GBM reprogramming synergizes with an intratumor sPD1. A, A schematic of the sPD1 and the NC constructs. The NC construct lacks the PDL1-binding motif in the immunoglobulin V (IgV) sequence. B, Representative immunoblots of his-tagged sPD1 and NC proteins in supernatants of KR158 GBM cells co-transfected with dP-mF3 and sPD1 or NC. Only sPD1 is detectable by the anti-PD1 mAb RMPI14 specific for the PDL1-binding motif absent in NC. C–F, A schema detailing the combination of dP-mF3-reprogrammed KR158 GBM and intratumor sPD1 (C), specifically creating a therapeutic synergy that delayed tumor growth as measured by BLI (D) and quantified in a line graph (E), and significantly extended survival as measured by Kaplan–Meier estimates (F). N = 10 per group. Log-rank test was used to compare survival rates. *, P < 0.05; **, P < 0.01; ****, P < 0.0001; ns, not significant. wtPD1, wild-type PD1.
Synergy between dP-mF3-reprogrammed GBM and a tumor RNA-based DC vaccine
A leading cause of failure of current immunotherapy in GBM and other non-immunogenic tumors is the frequent intratumor loss of heterozygosity at HLA loci, especially MHCI, leading to immune escape (7). Therefore, we asked whether the considerable MHCI upregulation induced by dP-mF3 or dP-ev alone in GBM cells (Fig. 1F) may enhance tumor-cell detection and killing by tumor-specific CD8+ CTLs generated by a tumor DC-based vaccine, in addition to the contribution of the intratumor iDC-APCs in turning the GBM TME immunoreactive. We again focused on the KR158 model because of its resistance to DC-based vaccines and killing by tumor-specific CD8+ CTLs due in large part to its immunosuppressed TME and profound MHCI downregulation (23, 26), similar to human GBM. To produce KR158-specific CD8+ CTLs, splenocytes from healthy syngeneic C57BL/6J mice that had been vaccinated with total KR158 RNA-pulsed syngeneic splenic DCs were incubated ex vivo with the same KR158 RNA-pulsed DCs in the presence of exogenous IL2 to specifically activate and expand KR158-specific effector CTLs, followed by purification using magnetic beads for CD3 (26). Expanded KR158-specific effector T cells, of which greater than 96% were CD8+ (Supplementary Fig. S14A) were then cocultured with reprogrammed KR158 target cells for 24 hours at increasing E:T ratios of 1:10 and 1:1, and KR158-specific CTL activation and cytotoxicity were measured by CD69 and INFγ production and rates of early apoptotic (Annexin V binding) and residual live KR158 cells, respectively (Fig. 6A). The activated fraction of KR158-specific CD8+ CTLs and the intensity of KR158-specific cytotoxicity strongly correlated with the E:T ratio and the relative expression levels of MHCI (absent/low, intermediate, and high) in the three target KR158 populations (ev-ev, dP-ev, and dP-mF3), respectively (Fig. 6B and C and Fig. 1F). The low recognition of the parent KR158 cells by KR158-specific effector CTLs was most likely due to the profound downregulation of MHCI in KR158 GBM cells, as these enriched effector CTLs re-engaged their parent KR158 targets much more efficiently after these KR158 had been treated with 100-ng/mL INFγ for 3 hours to induce upregulation of MHCI (Supplementary Fig. S14B; ref 61). In addition, a small population of KR158-specific CD4+ T cells was also found to be activated, presumably due to the iDC-APCs present in the target dP-ev- and dP-mF3-expressing KR158. Although we cannot rule out a contribution of iDC-APCs and these few activated CD4+ T cells to the observed CD8+ CTL activation and KR158 target cytotoxicity, these effects were deemed to be minimal because dP-mF3-expressing target KR158 cells with a much higher percentage of functional iDC-APCs with cross-presenting capability led to only a small increase in CD8+ CTL activation compared with dP-ev alone in the high E:T ratio of 1:1 and no increase in CD4+ T-cell activation at E:T of 1:1 compared with 1:10 was detected despite a significant increase in cytotoxicity (Fig. 6B and C).
GBM reprogramming synergizes with a DC-based tumor vaccine (DC Vax). A–C, A schema details the ex vivo production of purified KR158-specific CTLs following a KR158 RNA-pulsed DC vaccine to assess cytotoxic synergy with reprogrammed KR158 target cells (A). Bar graphs showing CD69 and IFNγ expression of KR158-specific CD8+ and CD4+ T cells and KR158-specific cytotoxicity measured by apoptotic (surface Annexin V binding) and residual live tumor cells as percentages of live CD45− cells in cocultures of increasing ratios (E:T) of CTL effectors (E) to targets (T) expressing the indicated CFD combinations at 1:10 (B) and 1:1 (C). All experiments in triplicate were repeated at least three times. Data are represented as mean ± SEM. Analyses were performed using two-way ANOVA. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 00001. ns, not significant. D–F, A schema details the combination of reprogrammed KR158 GBM and a KR158 RNA-based DC vaccine (D), demonstrating a specific synergy between the DC Vax and dP-mF3-reprogrammed TME in delaying tumor growth as measured by BLI (E) and extending survival as measured by Kaplan–Meier estimates (F). Twenty-five percent of the maximal dP-mF3-reprogrammed TME condition was sufficient to create a synergy with the DC Vax in delaying tumor growth (G) and extending survival (H). N = 10 per group. Log-rank test was used to compare survival rates. *, P < 0.05; **, P < 0.01; ****, P < 0.0001; ns, not significant.
GBM reprogramming synergizes with a DC-based tumor vaccine (DC Vax). A–C, A schema details the ex vivo production of purified KR158-specific CTLs following a KR158 RNA-pulsed DC vaccine to assess cytotoxic synergy with reprogrammed KR158 target cells (A). Bar graphs showing CD69 and IFNγ expression of KR158-specific CD8+ and CD4+ T cells and KR158-specific cytotoxicity measured by apoptotic (surface Annexin V binding) and residual live tumor cells as percentages of live CD45− cells in cocultures of increasing ratios (E:T) of CTL effectors (E) to targets (T) expressing the indicated CFD combinations at 1:10 (B) and 1:1 (C). All experiments in triplicate were repeated at least three times. Data are represented as mean ± SEM. Analyses were performed using two-way ANOVA. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 00001. ns, not significant. D–F, A schema details the combination of reprogrammed KR158 GBM and a KR158 RNA-based DC vaccine (D), demonstrating a specific synergy between the DC Vax and dP-mF3-reprogrammed TME in delaying tumor growth as measured by BLI (E) and extending survival as measured by Kaplan–Meier estimates (F). Twenty-five percent of the maximal dP-mF3-reprogrammed TME condition was sufficient to create a synergy with the DC Vax in delaying tumor growth (G) and extending survival (H). N = 10 per group. Log-rank test was used to compare survival rates. *, P < 0.05; **, P < 0.01; ****, P < 0.0001; ns, not significant.
Next, we sought to confirm the potential synergy between the dP-mF3-reprogrammed GBM TME, including intratumor iDC-APCs and MHCI upregulation, and a systemic DC-based vaccine. We first established orthotopic KR158-luc GBM tumors with (dP-mF3) or without (ev-ev) iDC-APC reprogramming in syngeneic C57BL/6J mice, then administered three weekly intradermal KR158 RNA-pulsed DC vaccines (DC Vax) or saline control (PBS), and measured tumor growth by BLI and overall survival (Fig. 6D). No clinical benefit was observed in the DC Vax alone cohort when compared with the no treatment (ev-ev + PBS) control, consistent with the resistant phenotype of this non-immunogenic, MHCI-low/negative GBM model (26). dP-mF3 alone resulted in a modest but statistically significant improvement in survival as seen earlier (Fig. 4E). When dP-mF3 was combined with DC Vax, synergistic effects in both tumor growth control and overall survival were observed (Fig. 6E and F). To determine whether the synergy between intratumor dP-mF3 and the DC Vax was dependent on a minimum threshold of dP-mF3-induced TME reprogramming, we implanted mixtures of dP-mF3- and ev-ev-expressing KR158 cells at different ratios, followed by the DC Vax. A TME ratio of dP-mF3:ev-ev as low as 25%:75% was sufficient to synergize with the DC Vax to delay tumor growth (Fig. 6G) and prolong survival (Fig. 6H), with both efficacy parameters further improving by increasing the dP-mF3 fraction to 50% and 100%, thus supporting a reprogramming dose-dependent synergistic response.
In summary, the PU.1/IRF8/ID2/BATF3 (dP-mF3) combination is sufficient to reprogram murine GBM cells into iDC-APCs with canonical DC-like functions, thereby reheating the cold TME of GBM while also inducing widespread PD1/PDL1 and MHCI upregulation. Working together, these positive changes produced by dP-mF3 create a therapeutic synergy with existing immunotherapy like immune checkpoint blockade and DC-based tumor vaccination.
PU.1 plus IKZF1 are sufficient to convert human GBM cells to iDC-APCs with a DC-like profile
The observations of murine GBM-DC conversion led us to ask whether the conversion of human GBM cells to iDC-APCs could similarly be achieved using the same ML-directed approach. Given the notable fate plasticity of human GBM, our initial investigation was focused on ascertaining whether the transition from GBM cells to cells with DC-like markers could occur naturally in vivo to a sufficient degree and with enough distinctiveness to enable the identification of potential CFDs for this transformation. We collected five primary treatment-naïve GBM tumors and performed scRNA-seq using a graph-based clustering technique in the Seurat R package (28) and UMAP (29) for dimension reduction as above. Non-immune tumor cells and immune cells were well separated by CD45 expression. Tumor-specific markers had low to negligible expression in the immune compartments and vice versa. Although a small population of CD45+MHCII+ cells with weak tumor marker expression—specifically GFAP—existed, intermixing with the microglia and MDSC clusters, few, if any, of these cells carried other GBM markers (e.g., SOX11, MAP1B, EGFR, and PCSK1N) or DC-specific markers (e.g., FCER1A, CLEC9A, CLEC10A, CLEC4C, and LILRA4) that could allow for defining tumor-derived DC-like APCs (Supplementary Fig. S15), indicating that such spontaneous transition is exceedingly uncommon in vivo.
Due to the substantial differences between mice and humans, we suspected that a unique set of human CFDs would be required for the human GBM cell to iDC-APC conversion. In fact, whereas the human equivalents of dP-mF3—hP-hF3 and hdP-hF3—successfully induced a small increase in CD45+ cells, similar to PU.1 alone in mice, they failed to consistently generate CD45+MHCII+ cells in several human GBM cell lines tested. In addition, hdP was found to be toxic in prolonged culturing in these cells (Supplementary Fig. S16). To identify CFD combinations unique for the human GBM-DC conversion, we obtained bulk RNA-seq expression profiles of six independent human monocyte-derived DCs and three independent human GBM cell lines including LN428, LN827, and LN308 (23), and applied the same network-based ML algorithm. Sequencing output is shown in Supplementary Table S11.
Of the top 10 ranked CFDs predicted to convert the three human GBM cell lines to iDC-APCs (Supplementary Table S17), the myeloid lineage master regulator PU.1 was ranked first. This was followed by IRF4, a paralog of IRF8 and a master regulator for general immune functions thought to cooperate with PU.1 in promoting CD8+ DC development (62); the IKAROS family zinc finger 1 (IKZF1) and cathepsin Z (CTSZ), which are implicated in HSC differentiation, especially lymphocyte (63) and myeloid (64) development, respectively; three regulators of HSC maintenance, including the two members of the microphthalmia transcription factor (MITF) family TFEC (transcription factor EC) and MITF, and cellular repressor of E1A stimulated genes 1 (CREG1; refs. 65–67); two downstream factors Zinc finger protein 366 (ZNF366, also known as DC-SCRIPT or DC-specific transcript) and growth arrest-specific protein 7 (GAS7), which are thought to regulate conventional DC specification via IRF8/IRF4 and growth arrest and general phagocytosis machinery, respectively (62, 68); and CCAAT enhancer binding protein delta (CEBPD), a paralog of CEBPA, which forms a heterodimer with CEBPB upstream of PU.1 (69). Like the murine GBM-DC conversion, hPU.1 regulated two large subnetworks of the myeloid and DC fates (Fig. 7A). Of the remaining CFDs, IKZF1 and IRF4 had major links mostly to the myeloid and DC regulatory clusters, thus overlapping with PU.1-controlled subnetworks and suggesting possible cooperation between these CFDs in the human GBM-DC conversion (Fig. 7A). This prediction was consistent with recent reports demonstrating that IKZF1 uniquely cooperates with PU.1 depending on the timing of its expression during hematopoiesis to guide differentiation toward the myeloid lineage, specifically DCs, and away from the default lymphoid trajectory in the absence of PU.1 (70). As we predicted, IKZF1 showed strong cooperation with PU.1 in generating up to 3-fold and 30-fold higher CD45+ and double CD45+MHCII+ cells, respectively, compared with PU.1 alone in the human GBM cell lines LN428 and LN827 (Supplementary Fig. S17), whereas not inducing a lymphocytic lineage (Supplementary Fig. S18). Next, to test the contribution of the other CFDs to the PU.1/IKZF1 (PI) combination, we added one or more of the remaining CFDs and assessed the production of iDC-APCs (CD45+MHCII+). Expression of CFDs was confirmed by immunoblotting or qRT-PCR (if antibodies were unavailable; Supplementary Fig. S19). In all the combinations tested, PI remained the most efficient and consistent across the three GBM cell lines plus U87 GBM cells, which were not included in the ML prediction, in producing the highest fractions of CD45+ and CD45+MHCII+ cells with the highest expression of the co-receptors CD80/CD86, CD11c, and the general inflammatory state markers MHCI and PDL1 (Fig. 7B–H; Supplementary Fig. S20A). In fact, in U87 GBM cells, when compared with hP–hF3—the human equivalent of the murine dP-mF3—and to hPIB (PU.1/IRF8/BAT3)—a tumor-agnostic and somatic conversion combination (13, 14)—the PI combination was 50% to 100% more efficient at inducing CD45+ and CD45+MHCII+ cells with higher expression of MHCI and CD11c (Supplementary Fig. S21), thus attesting to the accuracy and efficiency of the ML approach.
Conversion of human GBM cells to iDC-APCs by PU.1 plus IKZF1 (PI). A, A 2D image of a 300 gene regulatory network for the conversion from human GBM cells to iDC-APCs showing the top 10 ranked CFDs and the three large pathways regulating various aspects of the myeloid and DC states. B–H, Surface expression of CD45, MHCII, CD80, and CD86 (B) and bar graphs showing frequency of CD45+ (C), CD45+MHCII+ (D), CD11c+ (E), CD86+ (F), CD80+ (G), and MHCI+ (H) cells in the four indicated human GBM cell lines at days 7–9 after being transduced with the indicated CFD combinations. PI = PU.1/IKZF1; PIIr = PU.1/IKZF1/IRF4; PIC = PU.1/IKZF1/CEBPD; hF6 = PU.1/IKZF1/IRF4/CEBPD/CTSZ/MITF. I, Bar graphs showing the frequency of CD45+, CD45+MHCII+, CD80+, CD86+, and MHCI+ cells in the three independent human GSC lines at days 7 after being transduced with the indicated CFD combinations.
Conversion of human GBM cells to iDC-APCs by PU.1 plus IKZF1 (PI). A, A 2D image of a 300 gene regulatory network for the conversion from human GBM cells to iDC-APCs showing the top 10 ranked CFDs and the three large pathways regulating various aspects of the myeloid and DC states. B–H, Surface expression of CD45, MHCII, CD80, and CD86 (B) and bar graphs showing frequency of CD45+ (C), CD45+MHCII+ (D), CD11c+ (E), CD86+ (F), CD80+ (G), and MHCI+ (H) cells in the four indicated human GBM cell lines at days 7–9 after being transduced with the indicated CFD combinations. PI = PU.1/IKZF1; PIIr = PU.1/IKZF1/IRF4; PIC = PU.1/IKZF1/CEBPD; hF6 = PU.1/IKZF1/IRF4/CEBPD/CTSZ/MITF. I, Bar graphs showing the frequency of CD45+, CD45+MHCII+, CD80+, CD86+, and MHCI+ cells in the three independent human GSC lines at days 7 after being transduced with the indicated CFD combinations.
To confirm that the predicted role of PI in the GBM-DC conversion is applicable to both GBM cells and GSCs and DC subtypes other than cDCs, we curated independent RNA-seq profiles of 10 human DCs of various subtypes, 55 human GBM cell lines, and 13 GSC lines representing different GBM subtypes (Supplementary Table S13), processed and analyzed them using the same ML workflow. Human GBM-DC conversion networks generated using these broader datasets showed large CFD similarity as the above and again confirmed the critical role for the top-ranked PU.1 and IKZF1 with overlapping regulatory roles in myeloid and DC-associated hubs. The other top-ranked CFDs shared analogous connections with PU.1 and IKZF1 (Supplementary Fig. S22), suggesting that they are likely dispensable for the conversion as predicted above. Indeed, human PI was confirmed to efficiently upregulate CD45+ and CD45+MHCII+, CD80/86+, and CD11c+ and the general inflammatory state MHCI and PDL1 in three independent human GSC lines (25). Again, adding other supporting CFDs to the PI combination did not change the production of iDC-APCs or the overall inflammatory condition (Fig. 7I; Supplementary Fig. S20B), further confirming the utility of the ML-identified PI combination.
Consistent with its being an optimal combination for the human GBM-DC conversion, PI expression produced the most consistently across the four human GBM lines the highest fraction of low proliferative cells (Supplementary Fig. S23A), the highest phagocytic capacity (Supplementary Fig. S23B), and the highest expression of key antigen processing and presenting machinery components (Supplementary Fig. S23C) and select DC-associated cytokines (Supplementary Fig. S23D), compared with all other combinations tested.
Taken all together, the ML-identified PU.1 plus IKZF1 (PI) combination is sufficient to produce iDC-APCs from human GBM cells and GSCs with markers and global expression patterns closely similar to natural human DCs.
Discussion
Our findings contribute to an emerging paradigm in cancer immunotherapy by demonstrating how GBM cells can be reprogrammed into iDC-APCs. This reprogramming empowers tumor cells to assume the critical functions of professional APCs, leveraging their position within the TME to engage directly with a wide spectrum of tumor neoantigens and, from there, recruit, prime, and activate tumor-specific CTLs. Such a transformation of the TME enables the body’s antitumor immune response to progress unimpeded inside the tumor, potentially harmonizing with and enhancing the effects of existing cancer immunotherapies.
Importantly, unlike the general myeloid lineage reprogramming process (10) with the C/EBPα and PU.1 combination and the tissue type-agnostic conversion method employing the PIB combination (13, 14), our approach relies on a robust network-based ML algorithm to provide higher precision and speed of CFD identification for conversions of specific cancer cells to iDC-APCs with cross-presentation capability, and potentially a higher safety margin. This is evidenced by the fact that CFDs for the murine and human GBM to iDC-APC conversions shared only partial overlaps, as would be expected from the differences between human de novo GBM and chemically induced murine GBM models. Differences in in vivo conversion efficacy and safety between these methods will need to be compared empirically in the future, however.
Given that fully reprogrammed (CD45+MHCII+) iDC-APCs have the lowest proliferative rates in the background of fast-dividing non-reprogrammed GBM cells, the measured in vitro conversion rates were likely an underestimate. In addition, in vivo conversion rates for the same CFD combinations may depend on the efficiency of CFD delivery and the tissue milieu in which iDC-APCs are generated and remain to be determined. However, professional APCs, especially DCs, are apex immune inducers and as such having a high frequency of these cells in the TME may not be necessary for effective antitumor immune activation. Indeed, even in highly immunogenic tumors like melanoma and renal cancer, DC ratios in the TME are often less than 1% to 2% of all tumor-infiltrating immune cells and well below 0.1% of total tumor cells (71). In keeping with this notion, our results showed that despite the modest conversion rates, fully reprogrammed iDC-APCs successfully reheated the TME and enhanced tumor control and survival, but only when they possessed antigen cross-presentation capability even if they were not of a full DC fate. To that end, our ML computational suite can be used to identify additional CFDs that can bridge the iDC-APC state with natural DCs, including CFDs that may need to be repressed rather than expressed to further enhance conversion depth and efficiency. Despite the high fate plasticity observed in GBM, the wide chasm separating the GBM and DC states seems too high a huddle for the GBM-DC transition to occur naturally, at least in the five primary GBM tumors examined. However, our results and those from others indicate that with the right CFDs, this transition can be forced and potentially harnessed for therapeutic objectives. A fundamental question that remains unclear is whether iDC-APCs may eventually revert to their prior malignant state because the underlying genetic defects are not deleted during conversion. To determine the long-term fate of iDC-APCs in the context of the immune reactions they induce will require a detailed cell fate tracking system. However, leukemic cells reprogrammed through the general myeloid lineage reprogramming method have been shown to retain their new myeloid state (72). Although these observations are reassuring, correcting genetic mutations is neither the focus for these technologies nor practical due to the heterogeneity of GBM. Indeed, the novelty and attractiveness of the iDC-APC reprogramming approach are its mutation-agnostic nature.
The synergy between the widespread MHCI upregulation produced by the iDC-APC reprogramming and a tumor DC-based vaccine embodies a potential practical application to reverse MHCI downregulation, a common mechanism of immune escape by tumors (7). Other cancer immunizing platforms like peptide or nanoparticle-based vaccines, or adoptive transfer of tumor-specific T cells can also be combined with the iDC-APC reprogramming to enhance efficacy. In both murine and human GBM conversions, there was also significant upregulation of immune checkpoints, especially the PD1/PDL1 axis—changes that reflect the general inflammatory activation of the iDC-APC reprogramming. Persistent expression of these receptors in the TME often leads to T-cell exhaustion (72). To that end, our data demonstrating a therapeutic synergy between iDC-APC reprogramming and an intratumor sPD1 further supports the integration of the reprogramming approach with existing tumor immunotherapy. Although the effect of standard radiation, chemotherapy, and tumor-treating fields on the iDC-APC TME reprogramming was not explored in this work, the timing and sequence of iDC-APCs with one or more of these treatment modalities may have substantial impacts on outcomes. For instance, the presence of intratumor iDC-APCs may help prepare the TME to capitalize on new neoantigens induced by these therapies (23, 73), especially in cancers with low mutational burdens like GBM. Developing a reliable mode of delivery will also be needed to move the iDC-APC approach into the clinic. Several recent vectors such as replication-competent retrovirus with high selectivity for cancer cells, non-integrating viruses like adenovirus and adeno-associated virus (74, 75) with favorable safety profiles in patients with brain tumors, or nanoparticles are excellent options. Lastly, the use of ML-based CFD identification for specific cell type conversion may help further reduce off-target conversions when compared with the less specific methods like cytokines and tissue type-agnostic conversion methods. Thus, this ML-based cell conversion platform can be applied to other tumors to create a new class of fate-based gene immunotherapy.
Authors’ Disclosures
S.B. Le and D.D. Tran report a patent for 62/952,725 pending and a patent for 62/586,655 pending. D.D. Tran reports grants and personal fees from Novocure and Monteris and grants from Sarepta, Lacerta, Merck, Novartis, Northwest Biotech, Stemline, and Tocagen outside the submitted work. No disclosures were reported by the other authors.
Authors’ Contributions
T. Liu: Conceptualization, data curation, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. D. Jin: Conceptualization, data curation, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. S.B. Le: Conceptualization, data curation, software, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. D. Chen: Resources, data curation, investigation, visualization, methodology, writing–review and editing. M. Sebastian: Data curation, investigation, visualization, writing–review and editing. A. Riva: Software, methodology. R. Liu: Data curation, investigation. D.D. Tran: Conceptualization, resources, data curation, software, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing.
Acknowledgments
We would like to thank the D.D. Tran laboratory for their insightful feedback and assistance. We would like to acknowledge the support of the Norris Comprehensive Cancer Center of USC and its core facilities and of the University of Florida’s core facilities. This work was supported in part by award number R42CA228875 to D.D. Tran, award number P30CA014089 to the USC Norris Comprehensive Cancer Center, and award number F30CA232641 to M. Sebastian from the National Cancer Institute. This work was also supported in part by award number 6BC04 from the Bankhead Coley Research Program of the Florida Department of Health to D.D. Tran.
Note: Supplementary data for this article are available at Cancer Immunology Research Online (http://cancerimmunolres.aacrjournals.org/).
References
Supplementary data
GeneRep’s performance (Supporting the Materials and Methods Section)
Conversion of murine GBM cells to iDC-APCs (Supporting Main Figure 1)
iDC-APCs upregulate the immune checkpoint PD-L1 (Supporting Main Figure 1F).
Murine iDCs-APCs share overlaps with murine DCs. (Supporting Main Figure 2A)
Murine iDCs-APCs lack B and T cell markers. (Supporting Main Figure 2A)
All antibodies used