Ewing Sarcoma and Osteosarcoma Have Distinct Immune Signatures and Intercellular Communication Networks

Primary bone sarcomas are rare pediatric cancers, with nearly 200 cases of Ewing sarcoma and nearly 600 cases of osteosarcoma diagnosed in the pediatric population annually in the United States (1, 2). Ewing sarcoma and osteosarcoma share a common tumor microenvironment (TME), with most primary tumors arising in the bone and most metastases occurring in the lung. Unfortunately, another commonality between these cancers is the exceedingly poor survival outcomes for patients presenting with metastatic disease or experiencing relapse. Despite these commonalities, Ewing sarcoma and osteosarcoma are quite genetically distinct (1, 2). Ewing sarcoma is driven by the fusion oncoprotein EWSR1::FLI1(1, 3). The presence of EWSR1::FLI1 causes transcriptional dysregulation and drives the sensitivity of these tumors to DNA damaging agents by promoting DNA R-loop formation (4). Ewing tumors have a low tumor mutational burden (∼0.5 mutations/MB). In contrast to Ewing sarcoma, osteosarcoma is not a fusion oncoprotein-driven cancer. It is a genetically chaotic tumor, displaying chromothripsis and frequent copy-number changes (5). This highly heterogeneous tumor does demonstrate some recurrent genetic aberrations, such as loss of TP53 and amplification of MYC (6).

Inherent challenges in studying bone sarcomas include: (i) the rarity of these cancers at single institutions, (ii) the historic lack of biopsies from metastatic and relapsed lesions, and (iii) the paucity of age-matched comparator datasets. Despite these challenges, understanding the TMEs of difficult-to-treat Ewing sarcoma and osteosarcoma is a priority, given the need for new treatment modalities to improve outcomes. In the past decade, immunotherapies have emerged as new modality for cancer treatment. Although checkpoint inhibitor–based immunotherapy has been successful in the treatment of many adult cancers such as melanoma, lung cancer (7), head and neck cancer, and some soft tissue sarcomas (8–11), these agents have not demonstrated meaningful clinical responses in adult patients with Ewing sarcoma or osteosarcoma to date (12–14). Notably, bone sarcomas are distinct biologically from soft tissue sarcomas. The examination of non-checkpoint inhibitor–based immunotherapeutic strategies or novel immunotherapy/chemotherapy combinations is still largely lacking for Ewing sarcoma and osteosarcoma.

Our current understanding of the immune profiles of adolescent bone sarcomas is limited, especially for Ewing sarcoma (15). Most studies of the TME in Ewing sarcoma have sought to link immune cell frequencies by IHC to prognosis. For example, an analysis of 27 therapy-naïve Ewing tumors suggested those with increased CD8⁺ T-cell infiltration correlated with improved outcomes (16). A larger study demonstrated that 15.4% of Ewing tumors demonstrated significant CD8⁺ T-cell infiltration but no correlation with outcome was noted (17). In addition to CD8⁺ T-cell infiltration, CD68⁺ cell (i.e., macrophage) infiltration was associated with poor outcome in a study of 41 primary Ewing tumors (18). Although informative, these studies lack a granular dissection of the immune subsets in the TME. A detailed understanding of the Ewing TME, especially in the setting of metastatic and relapsed disease, is a major gap in knowledge in the field.

In osteosarcoma, studies have examined the role of macrophages in promoting progression and metastasis (19–26). More comprehensive studies of the osteosarcoma immune landscape include a recent study of 48 osteosarcoma samples where whole-genome sequencing, bulk RNA sequencing, and IHC were used to demonstrate groups of osteosarcoma tumors with low (C1), intermediate (C2), or high (C3) immune infiltration. A negative association of PARP2 expression and immune infiltration was noted (27). Single-cell RNA sequencing (scRNAseq) was recently utilized to construct an osteosarcoma cell atlas using 11 osteosarcoma tumor samples at the time of local control. This work highlighted the TIGIT pathway as a promising single therapeutic target for osteosarcoma (28). Although the immunobiology of osteosarcoma is relatively better understood as compared with that of Ewing sarcoma, large gaps in knowledge persist. For example, scRNAseq analyses to date have not focused on immune cell populations, thus making the immune infiltrate data captured in osteosarcoma scRNAseq studies quite limited. Studies often analyze samples at the time of local control when recent chemotherapy is actively modulating the immune response. Finally, age-matched comparators are often not utilized, an important limitation to note when examining immune findings in adolescent patients given the age-related changes that occur in the TME.

Here, we conduct an immune-cell focused scRNAseq analysis of blood and tumors from adolescent patients with Ewing sarcoma or osteosarcoma without recent chemotherapy administration to define immunologic profiles. Given the limited knowledge about the immune landscape of Ewing sarcoma, we sought to interrogate peripheral and tumor-infiltrating leukocytes during primary disease and relapse and to cross-compare with osteosarcoma. Osteosarcoma samples were included as a comparator (given the shared TMEs and patient ages between patients with Ewing sarcoma and osteosarcoma) and to pursue an immune-centric profile of osteosarcoma. We hypothesized that prior exposure to DNA damage and/or the lung microenvironment may result in enhanced immune infiltration in Ewing tumors upon relapse. Given promising early findings suggesting that Ewing tumor immune infiltrates evolve over time, we sought to expand upon these preliminary findings and utilize a multi-modality approach to conduct a thorough analysis of Ewing sarcoma immunobiology. Given the rarity of these tumors even at large pediatric hospitals, we use our scRNAseq data to deconvolute additional bone sarcoma bulk RNA-seq datasets to bioinformatically gain a deeper understanding of factors driving immune cell composition in Ewing sarcoma and osteosarcoma. We leverage immunofluorescence imaging of primary and recurrent disease to evaluate immune infiltration and identify potential drivers of immune infiltration in recurrent disease using cell–cell communication analysis of chemokines and receptors across immune populations from our scRNAseq data. This work is an imperative step toward identifying and prioritizing logical immunotherapeutic interventions to pair with treatments such as chemotherapy or radiation for the treatment of advanced pediatric bone sarcomas.

Adolescent and young adult patients diagnosed with primary or recurrent Ewing sarcoma or osteosarcoma were enrolled in the Musculoskeletal Oncology Research Registry and Tissue Bank (STUDY20010034), a single institution biobank. For the samples utilized in this study, patient ages at the time of sample acquisition ranged from 13 to 19 years. In addition to acquisition of tumor tissue, blood specimens were obtained by venipuncture and collected using EDTA as an anticoagulant. Freshly resected tumor specimens were collected in phenol-free RPMI and were processed within hours of resection to analyze immune cell infiltration (STUDY19030108). All samples analyzed were collected prior to initiation of chemotherapy or minimally 6 months off chemotherapy to reduce the influence of chemotherapy on immune signatures. Tissue samples preserved in formalin and paraffin embedded were also utilized. Blood from healthy anonymous adolescent blood bank donors was acquired (IRB exemption PRO17090430). In addition, a cohort of formalin-fixed, paraffin-embedded paired samples (primary disease and relapse) from patients with Ewing sarcoma were obtained from the Children's Oncology Group Biorepository for immunofluorescence analysis of immune cell infiltration (IRB exemption PRO18030271). Some of these specimens were derived from blocks that had previously been used to make tissue microarrays, thus limiting tumor material available for analyses in those cases. Formalin-fixed, paraffin-embedded lung tissue from adolescent nonsmokers was acquired from The University of Pittsburgh CORID #451 as a control.

Peripheral blood mononuclear cells (PBMC) were isolated from whole blood by Ficoll–Hypaque density gradient centrifugation as described previously (29). Briefly, whole blood was diluted with Hanks-Buffered saline solution and was layered on Ficoll–Hypaque. Samples were spun for 20 minutes at 400 × g with the brake off. PBMC were collected from the cell/Ficoll interface and were washed with RPMI supplemented with 2% FBS. Tissue samples were manually dissected into 1 mm² pieces and were incubated for 15 minutes at 37°C in 5% CO₂ with 50 mg/mL Liberase DL (Millipore Sigma, Catalog No. 5466202001). In instances where bone fragments were present in the sample, tissue was resected from bone and the bone fragments were incubated along with the 1 mm² dissected tissue fragments. Single-cell suspensions were then passed over a 40 μmol/L filter to remove debris. Viable cells were counted from PBMC and tumor-infiltrating CD45⁺ cells using a Nexcelom Cellometer with AO/PI staining.

After creating single-cell suspensions, cells were stained with 1:100 anti-CD45:PE (BioLegend, Catalog No. 368510, RRID:AB_2566370) at 4°C for 15 minutes in PBS supplemented with 10% FBS. Cells were then washed in PBS + 10% FBS and were incubated at 4°C in PBS with 1:4,000 viability dye eFluor780 (eBioscience, Catalog No. 65–0863–14). Samples were sorted on a MoFlo Astrios high-speed cell sorter in the Hillman Cancer Center Flow Cytometry Core. Cells were sorted on the basis of forward scatter and side scatter parameters associated with leukocytes and for CD45 positivity and exclusion of viability dye.

Following sorting of live CD45⁺ (i.e., all immune) cells, cells were once again counted on the Nexcelom Cellometer with AO/PI for viability. Cells were then loaded immediately into the 10× Controller as per the manufacturer's protocol. Both 3′ v2 and 5′ v1 kits were used for library generation. Following droplet generation, samples were reverse transcribed as per the manufacturer's protocol. Drops were then broken, samples were cleaned up by SPRI selection, and libraries were amplified by PCR. The manufacturer's protocol was followed for generation of sequencing libraries with i7 single indices. For some samples, peripheral blood and tumor samples were multiplexed by cell hashing using total-seq A or total-seq C antibodies and run on the same lane of the 10× to generate multiplexed libraries (30). These CITEseq libraries were separated from gene expression libraries during SPRI selection after the initial cDNA amplification and were applied with i7 indices to generate sequencing-ready libraries.

Sequencing-ready libraries were diluted to 2 pmol/L and pooled for sequencing on either a NextSeq 550 at the University of Pittsburgh Genomics Core or a NovaSeq 6000 at the UPMC Genome Core. The read pattern used for all samples was as follows: read 1, 26 cycles; i7 index, 8 cycles; read 2, 98 cycles.

Raw Illumina runs were downloaded from the sequencing cores to the University of Pittsburgh High Throughput Computer cluster. Bcl2fastq was used to demultiplex raw runs into FASTQ files based on i7 indices. Following demultiplexing, samples were aligned to GRCh38 and feature barcode expression matrices were generated for downstream analysis using cellranger v3.1.0. CITEseq samples were aligned using CITEseqCount with known total-seq antibody barcodes as input.

For samples in which cell hashing was performed, gene expression libraries were first demultiplexed into individual samples. This was accomplished by reading the aligned CITEseq files into the R package Seurat (v3). CLR normalization was then performed, and k means clustering was used to identify two separate peaks for each total-seq barcode as described previously (31). The mean of these two peaks was used as the cutpoint to separate cells as positive or negative for each barcode. Cells that were single-positive for the expected barcodes were assigned to each respective multiplexed sample. For analysis of gene expression barcodes, aligned filtered feature barcode expression matrices for each sample were aggregated into a single object using Seurat. Samples were log-normalized based and integrated using Seurat's integration workflow (32) based on 5′ versus 3′ 10× Genomics’ chemistries. Briefly, the integration workflow was used to identify anchors across the top 2,000 shared variable features between samples and reciprocal PCA was used to combine samples generated with these two different chemistries. Following integration, principal component analysis was performed, and significant PCs were selected heuristically. All consecutive PCs were included until the variance explained was negligible. UMAPs were then generated using the significant PCs. Louvain clustering as implemented in Seurat v3 was used for the initial clustering of all cells. The number of clusters was selected on the basis of the minimum number of clusters that defined biological differences in the data. Clustree (33) was used to aid in the selection of the number of clusters by tracking cluster genesis and cluster stability with increasing resolution. Cell types were identified on the basis of canonical expression profiles within clusters. For cell subset analyses (i.e., myeloid cells and CD8⁺ T cells), we bioinformatically isolated cell subsets and performed integration as described above.

Differential gene expression analysis was used to identify significantly differentially expressed genes between clusters and cell types. The nonparametric Wilcoxon rank sum test was used to identify statistically significantly differentially expressed genes.

Variant call files from whole exome-seq were obtained from PeCAN. Mutation annotation format (MAF) files were created from VCF files using MAFtools. MAF files were directly downloaded from TARGET-OS. Tumor mutation burden was calculated by summing the total number of nonsynonymous mutations and dividing by the size of the protein-coding portion of the human genome (i.e., 30 megabases).

The R package ESTIMATE (34) was used to infer ImmuneScore, StromalScore, and TumorPurity from bulk RNA-seq datasets from PECAN and TARGET-OS.

The R package MCPcounter (35) was used to infer the abundance of T cells, CD8⁺ T cells, cytotoxic lymphocytes, NK cells, and cells from the monocytic lineage from PECAN and TARGET-OS.

CIBERSORTx (36) was used to characterize immune cell fractions from bulk mRNAseq data obtained from PECAN and TARGET-OS. A signature matrix was defined using cell types derived from our overall scRNAseq dataset. S-batch mode was employed in CIBERSORTx to batch correct the signatures matrix (generated with scRNAseq data) to apply to the mixture matrix (derived from bulk mRNAseq data). CIBERSORTx was run locally on the University of Pittsburgh high-throughput computer cluster via a singularity instance.

The R package singleseqgset (29) was used to assess gene set enrichment across myeloid cell subpopulations in the TME. Hallmark gene sets from the Molecular Signatures Database were used as the query gene sets. Progenitor and exhaustion gene sets were obtained from Beltra and colleagues (37) and tissue resident CD8⁺ T-cell gene sets were obtained from Kumar and colleagues (38).

Cell–cell communication networks were identified using an updated version of the previously described R package celltalker (29). Briefly, ligands and receptors were identified across all immune lineages. Genes for known ligands and receptors were chosen as described previously (39). Mean expression levels of ligands and receptors were considered jointly to prioritize interactions that were likely to occur in the TME. The R package circlize (40) was used to construct circos plots to indicate ligand and receptor interactions.

To determine ligands produced by CD14⁺CD16⁺ macrophages that interact with other immune cell types, we leveraged the R package nichenet (41). The ligand target matrix, ligand/receptor networks, and ligand/receptor weighted networks were derived from https://zenodo.org/record/3260758#.Yt6oeezMIb0. Ligand/receptor interactions were modeled using scRNAseq samples derived from recurrent disease, with CD14⁺CD16⁺ macrophages as the “sender” cells and all other immune populations as “receiver” cells. Genes of interest were defined as those that were upregulated in the TME from recurrent disease versus the peripheral blood. After active ligands and their downstream signaling pathways were identified, we created a gene set for the pathway downstream of each ligand by taking all genes with an inferred activity greater than 0 for that given ligand. We then used linear modeling across the recurrent microarray samples to infer drivers of immune infiltration associated with immune infiltration and with each disease (e.g., recurrent osteosarcoma and recurrent Ewing sarcoma).

Multiplexed IHC analysis was performed as described previously (29). Briefly, slide mounted formalin-fixed, paraffin-embedded tumor tissue was deparaffinized. AR6 (Akoya Biosciences, Catalog No. AR600250ML) or AR9 (Akoya Biosciences AR900250ML) were utilized for heat induced antigen retrieval in serial cycles to stain tissues for CD4 (Biocare Medical LLC, Catalog No. API3209AA), CD8 (Biocare Medical LLC, Catalog No. ACI3160A), CD20 (Leica Biosystems, Catalog No. NCL-L-CD20-L26, RRID:AB_563521), CD68 (Cell Signaling Technology, Catalog No. 76437, RRID:AB_2799882), and FOXP3 (Cell Signaling Technology, Catalog No. 12653, RRID:AB_2797979) using Opal detection fluorophores (Akoya Biosciences, Catalog No. NEL811001KT). Cell nuclei were also stained with DAPI (4′,6-diamidino-2-phenylindole). Immunofluorescence images were captured on the Vectra (Perkin Elmer) and analyzed using Phenochart and InForm software.

Wilcoxon rank sum test was used to test if two sets of measures were different from one another. Spearman correlation was used for correlation analyses. Cox proportional hazard analysis with a log-rank test was performed for survival analyses. Differentially expressed genes were identified by Wilcoxon rank sum test. Linear models in R were used to identify NicheNet-inferred pathways associated with ImmuneScore independent of disease and those associated with a specific disease. All statistical tests with a two-sided α less than 0.05 was considered significant. Multiple comparisons were controlled using a FDR, where values less than 0.05 were considered significant.

The raw and processed scRNAseq data generated for this study is publicly available through the Gene Expression Omnibus database under accession number: GSE198896. Custom code used for analysis in the manuscript has been previously described and is available at github.com/arc85. Whole-genome sequencing data for pediatric relapse tumor samples used for analysis in this study were obtained from St. Jude Cloud (https://www.stjude.cloud; ref. 42). The TARGET-OS results presented here are based upon data generated by the Therapeutically Applicable Research to Generate Effective Treatments (https://ocg.cancer.gov/programs/target) initiative, phs000218. TARGET-OS data used for this analysis are available at https://portal.gdc.cancer.gov/projects.

Given our interest in understanding Ewing sarcoma immunobiology during disease progression, paraffin-embedded tumor samples from patients with single or multiple biopsies (primary tumor, metastasis, first relapse, or second relapse) were analyzed for immune cell infiltration using multiplexed immunofluorescence for CD4, CD8, CD68, CD20, FoxP3, CD45, and DAPI. In total, 23 samples were available, including 13 Ewing sarcomas and 10 osteosarcomas. Longitudinal immune infiltrate data were able to be analyzed from four patients: two with Ewing sarcoma (Fig. 1A and B) and two with osteosarcoma (Fig. 1C and D). Immune infiltrate data from samples without longitudinal pairs is shown in Supplementary Fig. S1. Analyses from our local cohort demonstrate: (i) fewer infiltrating immune cells in primary tumors versus relapsed tumors and (ii) an evolution of immune infiltration during relapse. Although CD68⁺ (i.e., macrophage) cell infiltration was noted to predominate in osteosarcoma immune infiltrates, we noted an increase in CD8⁺ cells in Ewing tumors upon relapse. The most common site of bone sarcoma metastasis is the lung. We note that both Ewing sarcoma and osteosarcoma lung relapses in our cohort had more immune cell infiltration compared with bone primary tumors. The immune cell composition of nonsmoker adolescent lungs from postmortem samples is included (Supplementary Fig. S2) as a reference control.

Ewing sarcoma specimens at any single institution are limited, as this is a rare pediatric cancer. In addition, as compared with osteosarcoma, much less is understood about changes in the immune composition of Ewing sarcoma when comparing primary disease versus relapse in individual patients. Thus, to better understand immune cell infiltration specifically in Ewing sarcoma beyond samples available from our local institution, paired patient primary and relapse Ewing tumor samples were acquired from the Children's Oncology Group (study AEWS20B1-Q). In total, 36 paired tumor samples from patients with Ewing sarcoma were obtained from the COG Biorepository and underwent multiplexed IF staining and analysis for CD4, CD8, CD68, CD20, FoxP3, CD45, and DAPI (Fig. 1E). Following multi-IHC staining, immune infiltrates were quantified. Figure 1F compares primary (P) bone tumor CD4, CD8, CD68, CD20, and FoxP3 cell frequencies as compared with paired lung relapses of Ewing sarcoma (R) in 11 patients. There is a statistically significant (P-value = 0.04) increase in CD8⁺ cells in relapsed Ewing tumors when analyzing P and R tumors in aggregate. Aggregate data are meaningful but do not highlight individual variance in immune infiltration between patients. To address this, these data were also compared across paired tumor samples from individual patients with a primary Ewing sarcoma of the bone and a relapse in the lung (Fig. 1G). In most cases, more immune cells were present in the TME upon relapse in these patients. Although the most common site of primary Ewing sarcoma is the bone and the most common site of relapse is the lung, additional scenarios do arise (such as soft tissue primary, bone relapses, rare metastasis to the brain). Supplementary Figure S3 includes the immune infiltration data of these “nonclassic” paired patient Ewing tumor samples. Taken together, these data suggest that lung relapses and some primary soft tissue presentations of Ewing sarcoma demonstrate higher levels of immune cell infiltration compared with primary bone disease. From a clinical perspective, these data suggest that site and/or phase of disease (primary vs. relapse) may impact which immunotherapeutic approaches are the most logical to consider for future testing in patients with Ewing sarcoma.

We found the evolution of Ewing tumor immune infiltration intriguing and next wanted to perform a more detailed, immune-focused analysis of bone sarcomas at diagnosis and relapse. Adolescent patients with a suspected new diagnosis (prechemotherapy, diagnostic biopsy) or relapse of a bone sarcoma were asked for consent to participate in the Musculoskeletal Oncology Tissue Bank and Registry (see Materials and Methods). Fresh tumor samples from the prechemotherapy (or minimally 6 months from most recent chemotherapy in the setting of relapse) were acquired (see Materials and Methods). Four Ewing sarcomas and five osteosarcomas met these criteria (Supplementary Table S1). All viable CD45⁺ cells (i.e., all leukocytes) were immediately isolated from fresh biopsy cores and paired patient PBMC samples. Single-cell transcriptomes were generated from sorted, live CD45⁺ cells via 10× Genomics scRNAseq (Materials and Methods). In addition, four age-matched healthy (16- to 18-year-old) PBMC controls were acquired (see Materials and Methods), processed in an identical manner, and included for comparison (Fig. 2A). In total 29,032 viable CD45⁺ cells were analyzed and canonical gene expression profiles were used to identify 20 distinct immune cell populations in blood and within tumors (Fig. 2B; Supplementary Table S2). After identifying the major immune lineages, we next visualized UMAPs of PBMC and tumor-infiltrating CD45⁺ cells across all sample groups (Fig. 2C). Comparison of UMAPs between PBMC and tumor-infiltrating CD45⁺ cells revealed distinct cell states present within tumors compared with peripheral blood. Specifically, some CD8⁺ T-cell states and macrophages were restricted to tumors based on visual inspection of cell distributions across the UMAPs. Next, we quantified the frequencies of immune cells in PBMC and tumor-infiltrating CD45⁺ cells across samples (Fig. 2D and E; Supplementary Table S3). We found that immune cell frequencies were generally similar in PBMC across sample groups, with the exception of a trend toward increased CD14⁺CD16⁻ and CD14⁺CD16⁺ monocyte populations in patients with bone sarcomas as compared with healthy adolescent donors. Immune cell frequencies were similar in tumor-infiltrating CD45⁺ cells from both Ewing sarcoma and osteosarcoma, but there was higher variance in cell frequencies in the TME versus peripheral blood. Overall, we identified and quantified the major immune lineages in peripheral blood and within the TME of patients with Ewing sarcoma and osteosarcoma based on canonical transcriptional profiles.

Following identification and characterization of the major immune lineages in Ewing sarcoma and osteosarcoma, we next utilized our scRNAseq dataset to perform an integrative immunogenomic analysis of bulk RNA-seq data from two publicly available cohorts (Fig. 3A). The PeCan dataset from St. Jude Children's Research Hospital (https://pecan.stjude.cloud; ref. 42) has bulk RNA-seq and whole-exome sequencing data from both patients with Ewing sarcoma and osteosarcoma, whereas TARGET-OS (https://ocg.cancer.gov/programs/target/projects/osteosarcoma) includes data from osteosarcoma only. These two genomics datasets were generated from primary, pretreatment tumor samples. First, we accessed whole-exome sequencing data available for primary, pretherapy Ewing sarcoma and osteosarcoma tumors through PeCan. Analysis of this data demonstrated that patients with osteosarcoma had a significantly higher tumor mutation burden compared with patients with Ewing sarcoma (Fig. 3B). We also leveraged ESTIMATE to quantify the ImmuneScore, StromalScore, and TumorPurity in both Ewing sarcoma and osteosarcoma (see Materials and Methods; ref. 34). Briefly, ESTIMATE uses bulk mRNAseq expression data and preselected gene sets to quantify the relative amount of immune infiltrate, stromal cells, and tumor cells present within tumors (Supplementary Table S4). This analysis revelated higher ImmuneScore and StromalScore in patients with osteosarcoma, and higher TumorPurity in patients with Ewing sarcoma (Supplementary Fig. S4A). We found tumor mutation burden and ImmuneScore to be significantly correlated, suggesting that tumor mutations at least partially drive immune infiltration in pediatric sarcomas (Fig. 3C). Overall, we confirmed that osteosarcomas are more immunogenic and have higher levels of overall immune infiltration as compared with Ewing sarcoma.

Next, we evaluated whether any measures of immune infiltration were associated with overall survival in TARGET-OS. This analysis was limited to osteosarcoma and data available through TARGET-OS, as patient outcome data are not included in PeCan. We first evaluated whether ESTIMATE scores of tumor-infiltrating immune cells were related to overall survival. Using Cox proportional hazard analysis (Materials and Methods), we found that a higher ImmuneScore was significantly associated with longer overall survival (Fig. 3D). Because ImmuneScore is a crude measure of immune infiltration, we next assessed whether more granular signatures of immune infiltration were associated with outcome. We used MCPcounter (35) to infer more detailed immune cell abundances from TARGET-OS (Supplementary Table S5). We found that higher levels of monocytic lineage abundance was associated with better overall survival (Fig. 3E).

We next sought to deconvolute the proportions of individual immune cell subsets in the TME of osteosarcoma and Ewing sarcoma. To do this, we generated a reference matrix from our scRNAseq data to use as input into CIBERSORTx (Materials and Methods). In this manner, we could deconvolute immune cell frequencies from bulk RNA-seq data using immune signatures derived from the same tissue environments. We first compared immune frequencies between osteosarcoma and Ewing sarcoma from PeCan (Fig. 3F; Supplementary Figs. S4B and S4C; Supplementary Table S6). Interestingly, despite osteosarcoma demonstrating overall higher levels of immune infiltration, we found that Ewing tumors have a higher proportion of CD8⁺ T-cell infiltration (Fig. 3D; Supplementary Fig. S4D). Conversely, osteosarcomas have higher levels of macrophage subsets (especially CD14⁺CD16⁺ macrophages) as compared with Ewing sarcoma (Supplementary Fig. S4D). In addition to more CD8⁺ T cells, Ewing tumors also had significantly higher frequencies of B cells (Supplementary Fig. S4D). These data suggest that immune infiltration scores are largely driven by a myeloid-dense immune infiltrate in osteosarcoma, whereas CD8⁺ T cells are relatively more predominant in the Ewing sarcoma TME.

We next assessed the relationships between proportions of immune cell subsets, ImmuneScore, myeloid lineage abundance, and patient outcome. CIBERSORTx was utilized to quantify immune cell frequencies in osteosarcomas (Supplementary Fig. S4E), with ImmuneScore and quantified using ESTIMATE and monocytic lineage quantified by MCPcounter. We found significant correlations between ImmuneScore, monocytic lineage, and CD14⁺CD16⁺ macrophages (Fig. 3G). Because CD14⁺CD16⁺ appear to be driving the underlying signatures of immune infiltration and monocytic lineage, we assessed whether CD14⁺CD16⁺ macrophages were associated with overall survival. This analysis demonstrated that higher proportions of this specific macrophage subset in the TME is associated with better overall survival (Fig. 3H). Taken together, we have refined immune cell proportions using several different approaches with increasing granularity to arrive at the conclusion that CD14⁺CD16⁺ macrophages drive immune infiltration.

Given the therapeutic targetability of CD8⁺ T-cell and macrophage populations, we next chose to focus on these immune cell populations to determine functional differences in these populations in Ewing sarcoma and osteosarcoma. Currently, the transcriptional signatures of CD8⁺ T cells in the TME of Ewing sarcoma and osteosarcoma remain incompletely characterized. To assess the transcriptional states of CD8⁺ T cells in these two tumor bone sarcoma subtypes, we bioinformatically isolated CD8⁺ T cells and performed analysis only on this population. UMAPs of peripheral blood and tumor-infiltrating CD8⁺ T cells revealed nine clusters (Fig. 4A; Materials and Methods). The top differentially expressed genes across clusters are shown in Fig. 4B (all differentially expressed genes are in Supplementary Table S7). When assessing the relationship between clusters and sample types, clusters 6, 7, 0, 2, and 4 were associated with PBMC and clusters 3, 5, and 1 were associated with the TME (Fig. 4C; Supplementary Table S8). Clusters 2 and 7 expressed genes associated with naïve or memory CD8⁺ T cells (SELL, CCR7), whereas cluster 4 expressed an effector signature including GZMB. Cluster 0 was associated with expression of KLR molecules including KLRG1 and KLFB1 as well as the effector molecule PRF1. Cluster 1 expressed CD69, consistent with tissue residence or early activation. Cluster 3 co-expressed co-inhibitory receptors LAG3 and PDCD1 (gene for PD1) in addition to IFNG, whereas cluster 5 expressed high levels of interferon response genes (e.g., IFITM1, IFITM2, and IFITM3). Finally, cluster 8 had a signature that partially overlapped with cluster 6 and also expressed genes associated with the cell cycle (i.e., MKI67), the costimulatory molecule ICOS, and the co-inhibitory receptor CTLA4. Cluster 8 also expressed CXCL13, which is associated with exhausted tumor reactive CD8⁺ T cells and can contribute to tertiary lymphoid structure formation (43). Overall, unsupervised clustering revealed a spectrum of cell states associated with a variety of CD8⁺ T-cell functionality in PBMC and in the TME.

To better contextualize the CD8⁺ T-cell states across clusters, we utilized gene sets that describe progressive states of progenitor exhausted and exhausted CD8⁺ T cells and gene sets that describe tissue residence of CD8⁺ T cells (Fig. 4D; Supplementary Table S9). This analysis revealed strong peripheral tissue signatures in cluster 1 and secondary lymphoid signatures in clusters 6 and 7. Clusters 0 and 4 were associated with intermediate exhaustion states, and cluster 3 was most strongly associated with terminal exhaustion. Leading edge analysis is commonly used to identify genes that are most strongly upregulated among genes in each gene set, and we performed leading edge analysis across the clusters most strongly enriched in the TME and cluster 4 (Fig. 4E–H; Supplementary Table S10). This analysis revealed co-expression of the co-inhibitory receptors TIGIT, LAG3, PDCD1 (gene for PD1), and HAVCR2 (gene for TIM3) in the terminal exhaustion gene set representing cluster 3 (Fig. 4E). Cluster 4 expressed genes associated with an effector response including GZMB, PRF1, and GZMA in addition to the transcription factor PRDM1 (gene for BLIMP1) and the co-inhibitory receptor HAVCR2 (TIM3) but no other co-inhibitory receptors (Fig. 4F). Cluster 1 expressed the integrins ITM2C and ITGA1 and the chemokine receptor CXCR3 consistent with recruitment to peripheral tissues and also expressed LGALS3 (gene for Galectin-3), which has been shown to inhibit CD8⁺ T-cell function. Finally, cluster 5 expressed SOCS3 and a series of interferon response genes, suggestive of a downregulation of response to cytokine signaling. Dissection of the genes associated with CD8⁺ T-cell tissue residence and exhaustion gene sets revealed distinct functions and gene expression profiles associated with CD8⁺ T cells in these pediatric bone sarcomas.

We next leveraged our scRNAseq data to further characterize myeloid subpopulations in patients and healthy controls. To achieve this, we bioinformatically isolated and reclustered monocytes and macrophages from all samples, which permits a granular analysis of heterogeneity in this immune subpopulation. After bioinformatically isolating the myeloid cells, UMAPs revelated distinct peripheral blood versus tumor-infiltrating myeloid cell populations (Fig. 5A) consistent with monocytes (in PBMC) and macrophages (within tumors). We performed gene set enrichment analysis using Hallmark gene sets from the Molecular Signatures Database (see Materials and Methods) across the major myeloid populations (Fig. 5B; Supplementary Table S11). This analysis revealed enrichment of pathways associated with immune activation, inflammation, and hypoxia in CD14⁺CD16⁺ macrophages. These data suggest that the CD14⁺CD16⁺ macrophages likely play an antitumor role, consistent with their association with better overall survival in patients with osteosarcoma from TARGET-OS. Conversely, CD14⁺CD16⁻ macrophages had enrichment for MTOR signaling and oxidative phosphorylation and were more closely related to the other macrophage populations by hierarchical clustering (Fig. 5B). Interestingly, the CD14⁺CD16⁺ monocytes had the highest enrichment for IFN signaling. Overall, the different myeloid populations have distinct biological functions by gene set enrichment analysis.

We next extended our analyses beyond canonical cell types and performed Louvian clustering (see Materials and Methods) to identify related groups of cells in an unsupervised manner. This approach identified a total of 12 distinct clusters (Fig. 5C). A total of five clusters had substantial enrichment in cells derived from PBMC, whereas the remaining seven had substantial enrichment in the TME. Interestingly, the CD14⁺CD16⁺ macrophages constituted 3 clusters (clusters 2, 4, and 8), suggesting that there are distinct subsets of these cells. We next assessed the correspondence between clusters and sample groups (Fig. 5D; Supplementary Table S12). We found that clusters 0, 1, 3, and 7 were associated with patient PBMC and cluster 5 was associated with PBMC from healthy donors. Clusters 2, 9, 6, and 11 were associated with macrophages from Ewing sarcoma, and cluster 8, 4, and 10 were shared between Ewing sarcoma and osteosarcoma. We found that the clusters representing CD14⁺CD16⁺ macrophages were differentially enriched in the tumors of patients with Ewing sarcoma versus osteosarcoma. Specifically, cluster 2 CD14⁺CD16⁺ macrophages were more frequent in Ewing tumors and cluster 8 CD14⁺CD16⁺ macrophages were more frequent in osteosarcoma. These results suggest that there are different subpopulations of this macrophage subset in Ewing sarcoma and osteosarcoma.

Finally, we evaluated the top differentially expressed genes across myeloid clusters (Fig. 5E; Supplementary Table S13). Interestingly, there were three distinct sets of gene expression profiles that defined the CD14⁺CD16⁺ macrophages. Cluster 2 cells expressed the compliment receptors C1QA, C1QB, and C1QC along with RNASE1 and FOLR2. Cluster 4 expressed genes related to metallothioneins including MT1G and MT1X, suggestive of a role in responding to either DNA damage or oxidative stress (44). Finally, Cluster 8 expressed higher levels of HLA molecules, suggestive of a role in antigen presentation. Overall, our analysis of myeloid cells characterized biological roles of known subpopulations and identified three unique and distinct states of CD14⁺CD16⁺ macrophages.

After identifying distinct cell states in myeloid cells and CD8⁺ T cells, we next identified putative cell–cell communication networks in the TME of these two pediatric bone sarcomas using CellTalker (see Materials and Methods). We first assessed the number of interactions between pairwise cell types in primary tumors (Supplementary Figs. S5A and S5B). When comparing Ewing sarcoma and osteosarcoma, we found that CD14⁺CD16⁺ macrophages appeared to interact broadly with many cell types in osteosarcoma, while CD14⁺CD16⁻ and CD14⁻CD16⁻ macrophages appeared to communicate with a diverse array of cell types in Ewing sarcoma. We next leveraged circos plots to visualize the top 20 ligand/receptor interactions by a joint average measurement combining the expression level of the ligand and receptor (see Materials and Methods). Interactions between and within immune lineages for Ewing sarcoma and osteosarcoma were analyzed (Supplementary Figs. S5C and S5D; Supplementary Table S14). Consistent with the number of interactions visualized in Supplementary Fig. S5A, we found that CD14⁺CD16⁺ macrophages produced several chemokines including CCL8, CCL3, and CCL2 that could impact a broad array of immune cells including CD8⁺ T cells, CD4⁺T_conv (conventional T cells), B cells, and NK cells. Interestingly, when performing the same analysis with Ewing sarcoma, we found that many of the top 20 ligand/receptor interactions were derived from CD14⁺CD16⁻ monocytes in the TME. These included the chemokines CXCL5, CXCL3, CXCL1, CCL7, and CCL2, which interacted with CD8⁺ T cells, CD1C dendritic cells (DCs), and plasmacytoid dendritic cells (pDCs). Our intercellular communication analysis of primary osteosarcoma and Ewing sarcoma yielded a distinct array of myeloid cell–secreted chemokines and broadly expressed cognate receptors as among the top 20 interactions.

After interrogating the differences in intercellular communication in primary disease, we next sought to evaluate immunologic drivers of immune infiltration in recurrent disease. We first evaluated whether CD14⁺CD16⁺ macrophages were associated with immune infiltration in recurrent disease. To achieve this, we performed CIBERSORTx-based deconvolution of immune cell proportions and ESTIMATE to infer immune infiltration across three datasets of recurrent osteosarcoma and Ewing sarcoma (45–47). This analysis showed that, similarly to primary disease, CD14⁺CD16⁺ macrophages are significantly correlated with ImmuneScore (Fig. 6A). We next sought to interrogate the ligands produced by CD14⁺CD16⁺ macrophages and the downstream signaling pathways they activate from our scRNAseq dataset (Fig. 6B and C; Supplementary Table S15). Of the putative cytokines produced by CD14⁺CD16⁺ macrophages, we found that 107 ligands had significant downstream activities. We visualized the downstream signaling pathways activated by the top 50 ligands (Fig. 6C).

After identifying the top 50 ligands produced by CD14⁺CD16⁺ macrophages and their downstream signaling pathways, we evaluated the relationship between ligands and immune infiltration in recurrent osteosarcoma and Ewing sarcoma. To achieve this, we created gene sets that encapsulated the downstream signaling pathways associated with each ligand (Materials and Methods; Supplementary Table S16) and then built linear models to identify CD14⁺CD16⁺ ligand activities that were associated with immune infiltration (Materials and Methods). We found a total of 16 ligand activity pathways (Fig. 6D) that were significantly associated with immune infiltration independent of disease (i.e., Ewing sarcoma or osteosarcoma) across the three datasets utilized in Fig. 6A. Interestingly, we found that CXCL10 and CXCL12 ligand activities were associated with immune infiltration specifically in osteosarcoma (Fig. 6E). No ligands were significantly associated with immune infiltration in Ewing sarcoma. These data indicate that there are common ligands produced by CD14⁺CD16⁺ macrophages that contribute to immune infiltration, but there are also chemokines (i.e., CXCL10 and CXCL12) that are produced by these macrophages to a significantly higher degree in osteosarcoma.

The gene sets activated downstream of these significant ligands were partially overlapping, so we evaluated correlations between downstream activities and identified three distinct clusters of downstream gene expression (Fig. 6F). This analysis showed that cluster 1 gene signatures were largely associated with ligands that drove immune infiltration in both Ewing sarcoma and osteosarcoma, whereas cluster 3 was associated with those that drove immune infiltration in osteosarcoma (and CXCL2). We found FN1 was the sole member of cluster 3. Finally, we found that the aggregate of the Cluster 1 and Cluster 3 gene signatures were both associated with ImmuneScore (Fig. 6G). Overall, these results indicate that Cluster 1 ligands are reflective of those produced by CD14⁺CD16⁺ macrophages in recurrent Ewing and osteosarcoma and that Cluster 3 ligands are produced by CD14⁺CD16⁺ macrophages specifically in osteosarcoma. The Cluster 1 and Cluster 3 ligand activities are reflective of two distinct pathways of immune infiltration and boosting CXCL10 and CXCL12 levels in the TME of recurrent Ewing sarcoma could lead to higher levels of immune infiltration.

In this study, we examined the evolution of tumor immune cell infiltrates upon bone sarcoma progression. We address an important knowledge gap in the bone sarcoma field by conducting a multimodality, immune-focused analysis of Ewing tumors. Ewing sarcoma shares important microenvironmental commonalities with osteosarcoma (immune system/patient age, site of primary disease, sites of relapse) and by directly comparing the immunobiology of these tumors, we were able to determine logical preclinical interventions to enhance the anti-tumor immune response specifically for each tumor type and derive questions to pose as biological correlates on future clinical trials.

First, we will consider logical CD8⁺ T-cell-based therapeutic interventions derived from our findings. In depth analysis of CD8⁺ T cells in this work revealed a spectrum of progenitor and exhaustion cell states in Ewing sarcoma and osteosarcoma. In Ewing sarcoma, there was a strong effector CD8⁺ T-cell signature in peripheral CD8⁺ T cells and CD8⁺ T cells in the TME expressed LGALS3 [gene for Galectin-3 (Gal3)] and a signature of peripheral tissue residence. Gal3 has been shown to impair IFNγ secretion from CD8⁺ T cells, thus limiting antitumor immunity (48). Targeting Gal3 to help improve disease control in patients with lung metastatic Ewing sarcoma is of interest preclinically. Gal3 inhibitors are currently in clinical trials for some advanced adult cancers (e.g., NCT02117362, NCT05240131). An inhaled Gal3 inhibitor, TD139, has also been used in patients with idiopathic pulmonary fibrosis (49), a potentially interesting formulation when considering the treatment of lung metastases. We also observed a peripheral activated CD8⁺ T-cell population in Ewing sarcoma. Conversely, in osteosarcoma, a subset of CD8⁺ T cells in the TME of osteosarcoma expressed TIGIT (consistent with recent findings; ref. 28) and were terminally exhausted. However, in our analysis, TIGIT was only one of several co-inhibitory receptors expressed on this exhausted population, which included LAG3, PD-1, and TIM3. These results show that multiple inhibitory mechanisms likely contribute to CD8⁺ T-cell exhaustion in osteosarcoma and suggest that “rescue” of CD8 cells from this exhausted state in osteosarcoma would likely require targeting multiple inhibitory receptors simultaneously. Of note, despite the usual robust expression of CTLA4 in exhausted T cell states, CTLA4 was not part of this group of co-inhibitory receptors in osteosarcoma, thus suggesting one plausible reason for the historical lack of response of osteosarcoma to CTLA4 inhibition. In addition to co-inhibitory receptor expression, we noted a population of CD8⁺ T cells in the osteosarcoma TME that expressed SOCS3 and IFN-stimulated genes (ISG), suggesting that protracted inflammatory cytokine signaling may contribute to reduced CD8⁺ T-cell responses (50) in addition to exhaustion. In addition to blocking coinhibitory receptors, another approach is to consider agents that enhance costimulatory receptors, such as ICOS. A number of agents (vopratelimab, utomilumab, etc.) to enhance costimulatory receptor signaling are currently being tested in phase I and II trials (51).

In addition to CD8⁺ T cells, immunotherapies are increasingly being designed to modulate the function of tumor-associated macrophages. In both Ewing sarcoma and osteosarcoma, CD14⁺CD16⁺ macrophages were noted in recurrent disease (osteosarcoma relatively more than Ewing sarcoma); however cytokine signaling of these macrophages demonstrated a key difference. Ewing tumor-associated macrophages during relapse lack the strong CXCL12 and CXCL10 signaling that is noted in macrophages from osteosarcomas. The presence of CXCR3 ligands like CXCL10 in the TME is known to be a critical component of the antitumor immune response during treatment with immune check point inhibitors (52). The striking absence of CXCL10 (and CXCL12) signatures in Ewing tumors prompts the question as to whether “therapeutic” delivery of CXCL10 to the Ewing TME could improve immune infiltration and the antitumor immune response at key points in therapy? As noted in the discussion on T-cell findings above, we observed a peripheral activated CD8⁺ T-cell population in Ewing sarcoma. This demonstrates that cellular antitumor immunity is present in Ewing sarcoma, but that this population may not effectively migrate to the TME. Given the differences in recruitment mechanisms observed to drive immune infiltration in osteosarcoma versus Ewing sarcoma (Fig. 6), boosting intratumoral levels of CXCL10/CXCL12 may promote increasing infiltration of this CD8⁺ T-cell subset into the TME. In addition to macrophages, both monocytes and dendritic cells appear to play a role in cell–cell communication in the Ewing TME. Circulating immune cells in the monocyte lineage are important in cancer (53) and can infiltrate tumors and differentiate into macrophage populations. In lung cancer, an increase in circulating classical monocytes has been demonstrated (54), and in breast cancer, increases in classical monocytes has been associated with decreased rates of pathologic complete response following neoadjuvant chemotherapy (55). Classical monocytes have the potential to contribute to extra-cellular matrix remodeling, as occurs during cancer metastasis. Upon inflammation, classical monocytes have the potential to exit the circulation and differentiate into dendritic cells or macrophages (53, 56, 57). This relationship is of interest given the chemokines secreted by CD14⁺CD16⁻ monocytes in primary and recurrent Ewing sarcoma and the cytokine signature of CD1c⁺ DCs in recurrent Ewing sarcoma.

When considering macrophages and osteosarcoma, CD14⁺CD16⁺ macrophages appear to play an important role in driving immune infiltration in primary osteosarcoma, likely through secretion of CCL2, CCL3, and CCL8 in primary disease based on our analyses. Overall, promoting activation of, and ongoing chemokine secretion from, CD14⁺CD16⁺ macrophages in osteosarcoma could lead to enhanced immune cell infiltration. Therapeutically, this could perhaps be achieved through macrophage-specific stimulation with anti-CD40 (58) or anti-CD47 ligands (59, 60). Other creative manipulation of macrophages in osteosarcoma could center around engineered macrophages. Engineered macrophages can readily infiltrate solid tumors although they do not proliferate and persist like engineered T cells (61). Like breast cancer, some osteosarcomas express HER2 (ERBB2) and trials including HER2-based agents for the treatment of osteosarcoma are currently ongoing (e.g., NCT04616560). HER2-targeting macrophages have been generated to help target and destroy HER2-expressing cancer cells and would be of great interest to preclinically test for the treatment of osteosarcoma (62).

Preclinical testing of logical immunotherapies or chemo-immunotherapy combinations for Ewing sarcoma presents a challenge given the historic lack of syngeneic or transgenic (immunocompetent) animal models in the field (63). However, new immunocompetent models are emerging, with a promising report of a zebrafish model of Ewing sarcoma (64). Our ongoing work is prioritizing the development of humanized mouse models of Ewing sarcoma to better capture the contributions of the immune microenvironment to disease progression and response to therapy. Detailed human Ewing tumor immune signatures, such as those generated in this study, could be leveraged to validate immune composition profiles of humanized mouse models of Ewing sarcoma in the future. Preclinical testing of immunotherapies in osteosarcoma presents relatively fewer challenges than Ewing sarcoma given the availability of immunocompetent models.

Biologically, our findings stress the need for trials incorporating immunotherapy for the treatment of bone sarcomas to consider the specific stage, location of disease, and timing around neoadjuvant therapy (65), given the evolution of immune infiltration noted in this study. Patient age is also likely an important factor when considering immunotherapeutic agents for the treatment of bone sarcomas. Future immunobiology questions of interest could include defining immune infiltration pattern differences between subgroups of Ewing sarcomas or osteosarcomas and prospectively validating the association of ImmuneScore and survival for osteosarcoma. Given the rarity of these cancers, answering such questions would require large, national tissue collection, and sequencing efforts in future clinical trials.

This work presents a multimodality, in depth analysis of Ewing immunobiology and demonstrates that Ewing sarcoma immune infiltration is more complex than often thought (1), given its potential to evolve upon disease progression. As new immunotherapies continue to emerge, these data will serve as a platform and resource to continue to explore new therapeutic options worthy of preclinical testing for patients with aggressive bone sarcomas.

A.R. Cillo reports personal fees from Abound Bio outside the submitted work. J. Daley reports grants from NIH T32HD071834, NIH K12HD052892, and Burroughs Wellcome Fund during the conduct of the study. R. Watters reports other support from Natera, Inc. outside the submitted work. T.C. Bruno reports other support from Walking Fish Therapeutics outside the submitted work. D.A.A. Vignali reports grants, personal fees, and other support from Novasenta; personal fees from Tizona, Incyte, and Bicara; other support from Trishula and Oncorus; and personal fees and other support from Werewolf and BMS outside the submitted work. K.M. Bailey reports grants from NIH (NCI), Hyundai Hope on Wheels, Alex's Lemonade Stand Foundation, and Morden Foundation during the conduct of the study. No disclosures were reported by the other authors.

A.R. Cillo: Conceptualization, data curation, software, formal analysis, validation, investigation, visualization, writing–original draft, writing–review and editing. E. Mukherjee: Data curation, formal analysis, writing–review and editing. N.G. Bailey: Data curation, software, formal analysis, writing–review and editing. S. Onkar: Methodology, writing–review and editing. J. Daley: Formal analysis, methodology, writing–review and editing. C. Salgado: Resources, validation, writing–review and editing. X. Li: Software, methodology, writing–review and editing. D. Liu: Software, methodology, writing–review and editing. S. Ranganathan: Resources, validation, writing–review and editing. M. Burgess: Resources, writing–review and editing. J. Sembrat: Resources, writing–review and editing. K. Weiss: Resources, methodology, writing–review and editing. R. Watters: Resources, writing–review and editing. T.C. Bruno: Conceptualization, data curation, software, formal analysis, supervision, methodology, writing–review and editing. D.A.A. Vignali: Conceptualization, resources, supervision, funding acquisition, visualization, methodology, writing–review and editing. K.M. Bailey: Conceptualization, resources, data curation, software, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing.

First and foremost, we would like to thank all of the patients diagnosed with Ewing sarcoma and osteosarcoma and their families for their willingness to participate in clinical trials, registries, and research studies. We would like to acknowledge The Center for Organ Recovery and Education and donor families for their gifts. The authors would like to acknowledge the collaboration with the Children's Oncology Group Ewing Sarcoma Biology Subcommittee in order to conduct the analysis of paired patient tumor specimens (AEWS20B1-Q). We thank Drs. Katherine Janeway and Patrick Grohar for their review of this work. We thank the University of Pittsburgh Center for Research Computing for the computational resources to complete this study. This project used the UPMC Hillman Cancer Center and Tissue and Research Pathology/Pitt Biospecimen Core shared resource (which is supported in part by award P30CA047904). We also wish to thank all the current and former members of the Vignali Lab (Vignali-lab.com; @Vignali_Lab) and the Bruno lab (@BcellBruno) for their constructive comments and advice during this project.

This work was supported by the University of Pittsburgh Cancer Immunology Training Program (T32 CA082084 to A.R. Cillo) and a Hillman Postdoctoral Fellowship for Innovative Cancer Research (to A.R. Cillo). K.M. Bailey was supported by the NCI (1K08CA252178), a Hyundai Hope on Wheels Young Investigator award, and an Alex's Lemonade Stand Million Mile Award (2020, 2021). K.M. Bailey would like to thank the Morden Foundation and the UPMC Children's Hospital Foundation for their continued support. The University of Pittsburgh holds a Physician-Scientist Institutional Award from the Burroughs Welcome Fund (JDD). Research reported here has been supported by the National Institutes of Health under their award number T32HD071834 “Research Training Program for Pediatric Subspecialty Fellows” PI: Terence S. Dermody, M.D. This work was also supported by the NIH R35 CA263850 and P01 AI108545 (to D.A.A. Vignali) and SRAs from UPMC and Novasenta (to D.A.A. Vignali).

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.