Large-Scale Comparative Analysis of Canine and Human Osteosarcomas Uncovers Conserved Clinically Relevant Tumor Microenvironment Subtypes Open Access

Osteosarcoma (OS) is a rare malignancy of the bone that mostly affects children and young adults. Despite aggressive treatment regimens that typically include additional surgery, radiotherapy, and chemotherapy, 5-year survival rates of patients with OS have plateaued around 60% to 70% (1). Hence, there is a clinically unmet need for biomarkers that can stratify patients into meaningful subgroups and open avenues for precision oncology. With the absence of a clear druggable, driving molecular target in OS (2), several groups have turned their attention towards the tumor microenvironment (TME) aiming to identify immune-related features associated with favorable prognosis (3–11). However, because of shortage of data, the prognostic value of most omic signatures has not been comprehensively evaluated in multiple cohorts, especially in the context of clinically established prognostic features such as the stage of disease at diagnosis and percent tumor necrosis following neoadjuvant chemotherapy (12). Furthermore, despite demonstrating immunogenic characteristics, recent immunotherapy-based clinical trials in OS have yielded disappointing results (13–17), suggesting that we still lack an in-depth understanding of the primary and metastatic OS TME.

To better understand disease progression in OS several preclinical murine models have been developed (18). However, even with increasing sophistication, they still fail to capture the sheer complexity of the tumor itself and a co-evolving microenvironment, inclusive of different tumor-associated stroma and immune components. Naturally occurring canine OS, on the other hand, have strong biological and molecular similarities with human OS (19–24) and are estimated to be at least 10 times more common in pet canines compared with humans (25). Pet canines develop tumors spontaneously in the face of an intact, educated immune system and comparable environmental exposures to humans. Hence, they are a potentially useful resource to study the OS TME and test the robustness and generalizability of biomarkers (11, 26).

In this work, we perform bulk transcriptomic deconvolution of 213 primary and 100 metastatic OS samples derived from 210 pet canines enrolled in COTC021/022 clinical trials conducted by the Comparative Oncology Trials Consortium (COTC; ref. 27). This analysis uncovers distinct clinically relevant TME subtypes, which we validate using a targeted canine-specific gene expression profiling tool (NanoString canine IO) and IHC staining of canine OS specimens. In contrast to prior work, which has extensively analyzed the immune infiltrate of OS, our analysis sheds insights into the role of additional overlooked components of the OS TME such as the extracellular matrix and tumor vasculature on overall disease progression. Most importantly, with the help of machine learning, we build a TME subtype-based risk stratification model and show that it generalizes from canine to human OS. Unlike previous studies, we comprehensively benchmark our proposed TME subtype-based risk stratification model against varying clinical features of both canine and human patients and demonstrate its prognostic value independent of previously established prognostic factors. Overall, this work highlights the tremendous potential of utilizing canines to derive clinically translatable biomarkers for human OS and inform personalized immunotherapeutic treatment strategies.

Canine OS tumor samples were curated from a large multisite clinical trial of 324 canines (27), which we refer to as the DOG² cohort. Prior to enrollment, patient tumor biopsies were examined and diagnosed as OS by veterinary anatomic pathologists at COTC institutions (https://ccr.cancer.gov/comparative-oncology-program/consortium, last accessed May 13, 2022). Each participating COTC member institution obtained and maintained approval from their respective Institutional Animal Care and Use Committee prior to enrolling patients in this trial. Dogs were actively recruited into the clinical trial over a span of 31 consecutive months, with patient-specific finalized outcomes reported up to 3 years postenrollment. Following diagnosis, additional treatment-naïve tumor tissue was collected at the time of surgical limb amputation by COTC investigators as a part of the trial schema. Canines were randomized to receive either standard of care or standard of care + adjuvant sirolimus (rapamycin) therapy. After completing sample QA/QC analysis and medical record review, a total of 190 frozen primary tumor samples and 8 frozen metastatic tumor samples from 195 out of 324 canines were deemed to be of sufficiently good quality (RIN > 8 and a total RNA quantity > 100 ng) for bulk RNA sequencing analysis. RNA was isolated from canine frozen tumor tissue in RNAlater using Qiagen Allprep DNA/RNA Mini Kit (Cat. #80204). The total RNA quality and quantity were assessed using Nanodrop 8000 (Thermofisher) and Agilent 4200 Tapestation with RNA Screen Tape (Cat. # 5067-5576) and RNA Screen Tape sample Buffer (Cat. #5067-5577). The mRNA library preparation and RNA Sequencing analysis were performed as previously described (28). The raw read counts for each gene were used for downstream differential expression analysis, using edgeR (v3.40.2; RRID:SCR_012802) and gene set enrichment analysis using the clusterProfiler package (v4.6.2; RRID:SCR_016884). The normalized read counts in log Transcripts Per Million [log₂(TPM + 1)] scale were utilized for downstream deconvolution analysis using MCP counter (29) and xCell (30).

Approximately 25% (83/324) of all canines enrolled in COTC021/022 trials (27) underwent a postmortem examination (31). During necropsy, additional frozen and formalin-fixed samples were collected from multiple tissues including any metastases. Formalin-fixed tissue samples were collected in 10% neutral buffered formalin and embedded in paraffin. Diagnosis of OS metastases was based on the finalized necropsy report, which included both gross and histologic examination of canine tissues by veterinary anatomic pathologists at the COTC institutions. Histologic analysis of each sample was performed to identify regions of high tumor content. For the majority of formalin-fixed, paraffin-embedded (FFPE) blocks, RNA was extracted from scrolls. Samples with <50% viable tumor tissue were collected via hand microdissection. RNA isolation from FFPE tissue was performed in the Molecular Histopathology Laboratory. Samples were placed into 1.5-mL microcentrifuge tubes. One drop of mineral oil was added to each tube and heated to 65°C for 2 minutes 160 μL of 0.2-mg/mL Proteinase K (Qiagen) and 320 U of RNaseOUT (ThermoFisher) in Buffer PKD (Qiagen) were added to each sample. Samples were placed into a Thermomixer R (Eppendorf) for 5 minutes at 65°C and then at 56°C for 68 hours with 5 minutes of shaking at 1,400 rpm every 30 minutes, adding 10 μL of 20 mg/mL Proteinase K (Ambion) to each sample every 24 hours. Lysates were incubated at 80°C for 15 minutes, centrifuged for 3 minutes at 16,000 g at room temperature, and transferred to fresh tubes leaving the paraffin behind. Lysates were DNA digested by adding 16 μL of DNase Booster Buffer (Qiagen) and 10-μL DNase I (Qiagen) with a 15-minute incubation at room temperature; 394 μL of Buffer RBC (Qiagen) and 12.5 μg of GenElute-LPA (Sigma cat#56575; RRID:SCR_000488) were added to the lysate. RNA was extracted with a QIAcube Connect MDx using the RNeasy Mini reagents and protocol. The addition of an equal amount of 70% ethanol excludes the small fragments of RNA. RNA was eluted in 30 μL of RNase-free water. Yield and purity were determined by a NanoDrop ONE spectrophotometer (ThermoFisher). RNA quality was determined by Agilent 2100 Bioanalyzer’s DV200 analysis. RNA samples with DV200 > 50% moved forward for NanoString analysis. This resulted in 91 FFPE-derived samples (39 dogs; 14 primary tumors, 77 metastases) and 24 frozen samples (11 dogs; nine primary, 15 metastases) from a total of 46 dogs. RNA was profiled using the Canine IO Panel and analyzed with the nCounter System (NanoString Technologies), according to the manufacturer’s protocol. RNA was loaded at 150 ng per sample. All hybridizations were 17 to 22 hours long, and all counts were gathered by scanning on HIGH mode for 280 fields of view per sample.

In total, 119 COTC021/022 canine cases with archived FFPE tissue and matched bulk transcriptomics were sent for hematoxylin and eosin (H&E), IHC, and Masson’s Trichrome staining at the College of Veterinary Medicine, Cornell University (Supplementary Fig. S1). The smallest section included in this study had an overall tissue area of 221,778.2 μm². Sequential FFPE sections from each case were immunolabeled for CD3 (T cells) and CD204 (Monocyte/Macrophage lineage; See Supplementary Table S1). Each FFPE section was manually reviewed by a board-certified veterinary pathologist (CJR) to assess whole slide image quality and identify regions of interest (ROI) containing tumors while excluding areas containing necrosis and artifacts such as tissue folds, background staining, dust, and bubbles. Out of 119 cases evaluated, 108 cases passed the manual review for CD204 IHC, 107 cases passed the manual review for CD3 IHC, 113 cases passed the manual review for Masson’s Trichrome, and 106 cases passed the review for H&E. For osteoid quantification, whole slide images were evaluated using HALO AI v3.6 (Indica Labs) with algorithms to detect and segment tumor osteoid within canine OS. Annotations were created to represent the spectrum of histologic features within tumors and to differentiate osteoid within tumor samples and training was performed on a DenseNet v2 classifier. Subsequent classified regions were validated and iteratively improved, as needed, by board-certified veterinary pathologists at the Molecular Histopathology Laboratory. Immune cells (CD3, CD204) and collagen (Masson’s trichrome) were quantified in an automated fashion within tumor ROIs using QuPath v0.5.0 (RRID:SCR_018257; See Table 1; Supplementary Tables S2–S5; Supplementary Fig. S2; ref. 32).

A total of eight publicly available human OS bulk gene expression and single-cell RNA sequencing (RNA-seq) datasets were analyzed in this study: TARGET (N = 87; ref. 33), GSE39055 (N = 37; ref. 34), GSE33383 (N = 84; ref. 35), GSE32981 (N = 23; ref. 36), GSE30699 (N = 76; ref. 37), GSE21257 (N = 53; ref. 38), GSE16091 (N = 34; ref. 11), and GSE152048 (N = 11; ref. 39). Preprocessed bulk RNA sequencing data of Human OS from the TARGET study was obtained from the UCSC Xena browser (www.xenabrowser.net/datasets; RRID:SCR_018938). Specifically, gene expression was quantified from RNA sequencing data using STAR (2.7.10b; RRID: SCR_004463) and gencode v22 human genome annotation and normalized to log₂(TPM + 1) scale. For more details on the expression quantification pipeline, see https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/. Normalized microarray gene expression datasets were curated based on an initial search on the Gene Expression Omnibus (GEO) database using the following keywords: “Osteosarcoma,” “biopsy,” and “homo sapiens.” Eventually, all datasets with sample size >20 or matched clinical metadata documenting progression-free survival outcomes were selected for further analysis. Besides bulk gene expression, preprocessed single-cell RNA-seq UMI count data of 11 human OS was obtained from the GEO database using the following study accession ID: GSE152048 (39). Single-cell RNA-seq UMI count data was further processed using Seurat v4.3.0.1 (RRID:SCR_007322; ref. 40) following the procedure described in the original study to remove cells that were potentially doublets or contained <300 expressed genes. Filtered single-cell gene expression profiles from each sample were then averaged and normalized to log₂(TPM + 1) scale to generate one bulk gene expression profile per sample.

To facilitate meaningful comparative analysis across multiple datasets while controlling for differences associated with different gene expression measurement platforms and sample preservation techniques, all samples were grouped together either by dataset of origin or sample preservation/platform and scaled by the mean and standard deviation of gene expression values in each group (41). The standardized z-scores were then provided as input to the MCP counter, a gene set enrichment-based deconvolution algorithm that estimates the relative abundance of each cell type or lineage based on the average normalized expression of a predefined set of marker genes unique to each cell type or lineage (29). MCP counter by default estimates the relative cell type abundances of the following cell types commonly found in the TME: CD8⁺T cells (CD8B), T cells (CD28, CD3D, CD3G, CD5, CD6, CHRM3-AS2, CTLA4, FLT3LG, ICOS, MAL, MGC40069, PBX4, SIRPG, THEMIS, TNFRSF25, and TRAT1), cytotoxic lymphocytes (CD8A, EOMES, FGFBP2, GNLY, KLRC3, KLRC4, and KLRD1), B cells (BANK1, CD19, CD22, CD79A, CR2, FCRL2, IGKC, MS4A1, and PAX5), NK cells (CD160, KIR2DL1, KIR2DL3, KIR2DL4, KIR3DL1, KIR3DS1, NCR1, PTGDR, and SH2D1B), monocytic lineage (ADAP2, CSF1R, FPR3, KYNU, PLA2G7, RASSF4, and TFEC), myeloid dendritic cells (CD1A, CD1B, CD1E, CLEC10A, CLIC2, and WFDC21P), neutrophils (CA4, CEACAM3, CXCR1, CXCR2, CYP4F3, FCGR3B, HAL, KCNJ15, MEGF9, SLC25A37, STEAP4, TECPR2, TLE3, TNFRSF10C, and VNN3), endothelial cells (ACVRL1, APLN, BCL6B, BMP6, BMX, CDH5, CLEC14A, DIPK2B, EDN1, ADGRL4, EMCN, ESAM, ESM1, FAM124B, HECW2, HHIP, KDR, MMRN1, MMRN2, MYCT1, PALMD, PEAR1, PGF, PLXNA2, PTPRB, ROBO4, C, SHANK3, SHE, TEK, TIE1, VEPH1, and VWF) and fibroblasts (COL1A1, COL3A1, COL6A1, COL6A2, DCN, GREM1, PAMR1, and TAGLN), which we refer to as ECM/stroma in this study given close relationship between expression of MCP counter-defined fibroblast markers and ECM production. The resulting relative abundances of each cell type correspond to the TME profile of each sample. For additional validation using a more extended reference cell atlas, we performed deconvolution using xCell (30), another highly cited state-of-the-art gene set enrichment-based deconvolution algorithm that defines marker genes for 64 distinct cell types or lineages expected to be present in normal or cancerous tissue samples. Besides MCP counter and xCell, there are several other more sophisticated deconvolution algorithms published in the literature (42–47). However, we report our results based on the MCP counter for the following reasons: (i) is simple and standardized implementation across various gene expression platforms compared with other deconvolution approaches, (ii) effectively mitigates issues associated with spill-over of markers from closely related cell types or lineages (48), and (iii) can be applied to bulk expression datasets of cancer types lacking a complete reference cell type atlas (48).

Following deconvolution analysis with an MCP counter, the standardized TME profiles of each sample from the discovery cohort were analyzed using the ConsensusClusterPlus package (RRID:SCR_016954; v1.62.0; ref. 49). Consensus clustering is a popular quantitative evidence-based approach frequently used to determine the number of clusters and robustness of cluster assignments in gene expression datasets while accounting for the impact of sampling bias (50). The method works by repeatedly subsampling a given dataset at random and applying a clustering algorithm (e.g., K means) to the subsampled version of the dataset to generate multiple clustering results. The final cluster assignments are determined based on a consensus, thereby enhancing their robustness compared with cluster assignments obtained from a single clustering run on the entire dataset. To quantify the stability of cluster assignments, pairwise consensus values are computed. These values represent the proportion of times two data points are grouped into the same cluster across different clustering runs. By analyzing the cumulative distribution function (CDF) of these consensus values for varying numbers of clusters (K), the optimal number of clusters can be determined, as shown by Monti and colleagues, through extensive validations on simulated and real-world data. Specifically, the optimal number of clusters (K) corresponds to the point where the area under the CDF curve begins to plateau, indicating limited evidence for the existence of additional stable clusters. The following parameters were specified as input for the consensus clustering analysis: distance = “euclidean,” clusterAlg = “kmeans,” pItem = 0.8, pFeature = 1, seed = 2,425, maxK = 10, and reps = 100. Following consensus clustering analysis, further inspection of the relative change in area under the CDFs for each choice of K (i.e., the number of clusters) suggested the existence of three stable clusters (Supplementary Fig. S3) labeled as follows: IE (immune enriched), IE-ECM (immune enriched dense extracellular matrix-like), and ID (immune desert) subtypes. For more details on the process of unsupervised class discovery using consensus clustering analysis, please see the ConsensusClusterPlus R package documentation. To enable robust classification of new samples into one of the three identified clusters in independent datasets, a naïve Bayes classifier (51) was trained on canine tumor data to predict its assigned cluster from its TME profile. We chose the naïve Bayes classifier because of its simplicity in training and implementation. For training, we only considered cell types whose abundances could be estimated across all datasets analyzed in this study. For added robustness against outliers, we incorporated a standard data augmentation strategy of randomly injecting gaussian noise (mean: 0, standard deviation: 0.01) into the training data prior to feeding it as input to the naïve Bayes classifier (52).

Single-cell clustering and classification data were utilized to estimate cell type proportions of each human OS sample from the following single-cell RNA-seq study: GSE152048 (39). Following cell type proportion calculation, all cell type proportions were standardized to z-scores for visualization in a heatmap (Supplementary Fig. S4A and S4B).

Raw read count data and the inferred TME subtype of each sample from the discovery cohort were provided as input to edgeR (v3.40.2; RRID:SCR_012802; ref. 53) using default parameter settings for subtype-specific differential expression analysis. For each TME subtype, all genes were ranked by their log fold change in expression estimated by edgeR. For a given subtype, genes ranked at the top of the list correspond to genes with high positive log fold change in expression associated with that subtype (i.e., have relatively higher expression in that subtype vs. other subtypes), whereas genes ranked at the bottom of the list correspond to genes with high negative log fold change in expression associated with that subtype (i.e., have relatively lower expression in that subtype vs. other subtypes). For each subtype, the ranked list of all genes was then provided as input to the standard gene set enrichment analysis pipeline implemented in clusterProfiler (v4.6.2; ref. 54) with the following parameters specified: (nPermSimple: 100,000, minGSSize = 10, and maxGSSize = 500) to estimate the relative enrichment of cancer hallmark pathways in each subtype. For each subtype, the top five enriched pathways ranked in decreasing order by the normalized enrichment score were selected for downstream biological interpretation.

To perform Kaplan–Meier and Cox proportional hazards survival analysis of canine clinical trial and human clinical datasets, we utilized the survival (v3.5.5) and survminer (v0.4.9) packages implemented in R. In the case of canines, adjuvant sirolimus treatment did not affect any of the trial endpoints, inclusive of disease-free interval (DFI) or overall survival. A post hoc analysis of the pharmacokinetic parameters from canines receiving standard of care + adjuvant sirolimus indicates a high degree of variability in exposures and a lower-than-expected trough blood level, below the predicted target of 10 μmol/L (27). Thus, these canine patients were considered together in the survival analyses presented herein. Four cases were excluded from the trial: 1702, 1706, 902, and 1701, hence had missing clinical data. Cases 902 and 1701 had metastases found at the time of surgery, Case 1702 had an erroneous diagnosis of osteosarcoma instead of fibrosarcoma, and Case 1706 was determined to be ineligible because of a preexisting heart condition. In the case of human OS, all 405 patients analyzed in this study received standard-of-care neoadjuvant chemotherapy, with the majority of these patients (N = 383) receiving a diagnostic biopsy of their primary tumor prior to neoadjuvant chemotherapy treatment. Hence these patients were considered together in downstream survival analysis. 185 out of these 383 patients had matched progression-free survival outcome information measured, which was standardized to the time interval (in days) from the time of diagnosis to the time of relapse or metastasis. This information was utilized to perform Kaplan–Meier survival analysis and plot progression-free survival trajectories of patients stratified by their TME subtype. Of these 185 patients, 83 patients had information on percent tumor necrosis following histopathological evaluation after neoadjuvant chemotherapy treatment (>90%, ≤ 90%) along with age, gender, and stage of disease at initial diagnosis (localized vs. metastatic). This information was utilized to fit a Cox proportional hazards regression model utilizing the following variables as features: age, gender, stage, percent necrosis, and TME subtype. For an independent cohort of patients [N = 53 out of 84 patients from GSE33383 (35)], a binary outcome indicating whether the patient progressed to metastatic disease within 5 years from initial diagnosis was measured. This information was utilized to estimate the proportion of patients developing metastasis within 5 years for each TME subtype.

The canine clinical trial datasets generated during and/or analyzed during the current study are available on the ICDC web portal (https://caninecommons.cancer.gov/#/study/COTC022) and GEO database (RRID:SCR_005012; accession ID: GSE238110). The human osteosarcoma datasets analyzed in this study are available from Genomic Data Commons (https://gdc-hub.s3.us-east-1.amazonaws.com/download/TARGET-OS.star_counts.tsv.gz) and GEO database (RRID:SCR_005012; accession ID: GSE16091, GSE21257, GSE30699, GSE32981, GSE33383, GSE39055, and GSE152048). Supplementary Table S6 details how each dataset was utilized for generating each figure in this study. All the codes required to reproduce the results of this study are deposited at https://github.com/spatkar94/doghuman_TME_deconvolution.git. Additional data queries may be directed to the corresponding author.

To comprehensively characterize different cellular components of the OS TME, frozen tumor samples were collected from dogs enrolled in the COTC021/22 clinical trials (27). A total of 190 primary tumor samples and eight metastatic tumor samples acquired from 195 dogs yielded sufficient RNA quality for RNAseq (See “Materials and Methods”), which we refer to as the DOG² RNA-seq cohort (Fig. 1A, See additional details in “Materials and Methods”; ref. 28). We performed bulk transcriptomic deconvolution of this cohort followed by consensus clustering analysis (49) to uncover three distinct TME subtypes of canine OS (Fig. 1B; Supplementary Fig. S3A–S3D): (i) IE, which consists of tumors that are enriched for cytotoxic T and NK cell populations; (ii) IE-ECM, which consists of tumors with a diverse milieu of infiltrating immune cells and in addition are enriched for endothelial and ECM/Stromal subpopulations; and (iii) ID, which are predominantly cold tumors with low levels of immune infiltrate. We observe similar clustering of patients using another state-of-the-art deconvolution algorithm, xCell (30) which considers a wider range of cell types and cell states (Supplementary Fig. S5A–S5D). Moreover, the observed TME subtypes are independently uncovered with NanoString data of 23 primary and 92 metastatic tumor samples from a subset of 46 canines from the DOG² cohort that passed quality control for NanoString analysis (Fig. 1C; Supplementary Fig. S6; “Materials and Methods”). Notably, we do not observe statistically significant associations between TME subtype and sample preparation (FFPE vs. Frozen: χ² test P value = 0.14), as well as tissue site (χ² test P value = 0.42). On performing subtype-specific differential expression and functional enrichment analysis (see “Materials and Methods”), we observe that the top five hallmark pathways enriched in the IE subtype are as follows: ALLOGRAFT_REJECTION, INTERFERON_ALPHA_RESPONSE, INTERFERON_GAMMA_RESPONSE, INFLAMMATORY_RESPONSE, and TNFA_SIGNALING_VIA_NFKB. The top five hallmark pathways enriched in the IE-ECM subtype are as follows: IL6-JAK-STAT3_SIGNALING, COAGULATION, INFLAMMATORY_RESPONSE, KRAS_SIGNALING_UP, and TNFA_SIGNALING_VIA_NFKB. Conversely, the ID subtype shows a complete lack of an immune response with the following top five enriched hallmark pathways: E2F_TARGETS, MYC_TARGETS_V1, G2M_CHECKPOINT, MYC_TARGETS_V2, and OXIDATIVE_PHOSPHORYLATION (Supplementary Tables S7–S9). Given that chromosomal instability (CIN) is also a characteristic feature of OS (55), we additionally estimated the level of CIN in each sample using a previously defined gene expression signature of CIN (56). Interestingly, we observe significant differences in the estimated expression of CIN across subtypes with the ID subtype having the highest estimated expression of CIN, closely followed by the IE-ECM and IE subtypes (Fig. 1E).

The IE subtype, which is characterized by significantly high levels of cytotoxic T and NK cells constitutes a minority of all immune infiltrated tumors (N = 11 from the DOG² RNA-seq cohort, N = 17 from the DOG² NanoString cohort) and has significantly higher expression levels of immune checkpoint genes: PD-1, PD-L1 CTLA4, LAG3, TIM3, and TIGIT in comparison to other subtypes (See Supplementary Fig. S7A and S7B). The IE-ECM subtype (N = 82 DOG² RNA-seq cohort, N = 54 DOG² NanoString cohort) constitutes a majority of the immune infiltrated tumors and is uniquely characterized by relatively high expression of ECM, endothelial cell markers, and in addition, genes involved in oncogenic KRAS signaling and epithelial–mesenchymal transition. We additionally performed quantitative imaging analysis of a subset of 119 cases from the DOG² cohort with bulk transcriptomics and matched whole slide H&E, Masson’s trichrome staining, and immunolabeling for CD3 and CD204 of their primary tumors (see “Materials and Methods”). Despite challenges associated with intratumor heterogeneity, this analysis reveals strong differences in overall levels of intratumor T cells and macrophages across the three subtypes. However, the differences in overall levels of intratumor collagen and osteoid were modest (Fig. 2A–E). Nevertheless, the observed trends in overall levels of each quantitative imaging feature were consistent with the characteristic features of transcriptomic-derived TME subtypes (ANOVA test P value < 0.05).

We next investigated if the TME subtype of the treatment-naïve primary tumor is predictive of clinical outcomes of canines from the DOG² cohort. All canines were metastasis-free at the time of diagnosis with most canines eventually developing metastatic disease during the clinical follow-up period. See “Materials and Methods” for additional details on sample curation and Supplementary Table S10 for all pertaining clinical metadata. Remarkably, the inferred TME subtypes of the primary tumor strongly stratify the overall survival and DFIs of canines with the IE subtype having the most favorable prognosis, closely followed by the IE-ECM subtype, and finally the ID subtype having the worst prognosis (Fig. 3A and B; overall survival log-rank test P value: 0.0012, DFI log-rank test P value: 0.0017). We additionally fit a multivariate Cox proportional hazards regression model to assess the predictive value of the primary TME subtype in the context of other clinically measured variables such as age, weight, primary tumor location (proximal vs. nonproximal to humerus), gender, alkaline phosphatase (ALP) levels and clinical trial treatment arm (standard of care chemotherapy with or without adjuvant sirolimus). Overall, the primary TME subtype has a significant effect on DFIs even after adjusting for other clinical factors (Fig. 3C and D). This suggests that the primary TME composition could serve as an important prognostic biomarker for metastatic disease progression and risk stratification in OS.

Despite the fact that OS frequently metastasizes to distinct organs, very little is known about its metastatic TME landscape and how it relates to primary tumors of patients receiving standard-of-care treatment. This is mainly because of challenges associated with curating tumor samples from human patients (1, 2). However, fewer barriers exist for collecting and analyzing tumor samples from canines. This makes canine OS a natural model for studying the evolution of the TME in these relatively rare tumors from primary to metastatic disease. Hence, we performed a comparative analysis of a small subset of 28 canine patients from the DOG² cohort, where both the primary and metastatic tumor samples were collected and passed RNA quality control (See “Materials and Methods”). To assess whether the primary tumor TME subtype influences the metastatic TME subtype, we constructed a contingency table of TME subtype assignments for paired primary and metastatic tumors from all 28 dogs. Overall, when comparing the TME composition of primary tumors to their matched metastases from a variety of geographic locations, no distinct trends emerged, suggesting that the observed TME subtypes could be both conserved and switched as tumors evolve from primary to metastatic disease (Fig. 4A; χ² test P value: 0.1059). However, it is important to note the presence of intratumor heterogeneity within primary as well as metastatic tumor samples in this cohort and differences in platforms, which could potentially mask clinically significant trends in TME evolution (See Fig. 4B). Nevertheless, we observed that the dogs with more favorable clinical outcomes (overall survival times > median of cohort) had significantly higher levels of cytotoxic T-cell infiltration in various metastatic sites compared with canines with less favorable outcomes (overall survival times < median), irrespective of metastatic site (Fig. 4B). These results suggest that in addition to the primary tumor, an increased presence of cytotoxic T cells at different metastatic sites leads to more favorable overall survival outcomes.

Although canine OS has been shown to have strong similarities with human OS, it is yet unknown whether biomarkers derived from the analysis of canine data can reliably stratify clinical outcomes of human patients. To test this hypothesis, we trained a simple naïve Bayes classifier (51) to predict the TME subtype of each tumor sample based on its MCP counterestimated TME composition (Fig. 5A; see “Materials and Methods” for more details). The TME subtype classifier, on average, achieved a 10-fold cross-validation accuracy of 91.5%, CI, 90%–93%, and kappa score of 0.872, CI, 0.85–0.89, when trained and validated on canine tumor samples from the DOG² cohort. We then applied the trained naïve Bayes classifier on seven independent publicly available human OS bulk gene expression datasets (N = 394) and in addition one single-cell dataset (N = 11; refs. 11, 34–39, 57). The majority of primary tumor samples curated from public datasets represented prechemotherapy diagnostic samples (N = 383 of 405). Detailed clinical metadata for each dataset evaluated is provided in Supplementary Table S11. Figure 5B shows a heatmap of MCP counter-estimated relative cell type abundances in human OS samples along with the classifier-predicted TME subtype. Overall, we observe conserved trends in estimated cell type composition, CIN, and immune checkpoint gene expression levels across classifier-predicted TME subtypes (Fig. 5B and C; Supplementary Figs. S8A, S8B, and S9). We additionally looked at the TME composition of 11 human OS tumor samples that were previously profiled at a single-cell level (39). Despite the low sample size, we observe that tumor samples classified as the IE subtype (N = 1) are uniquely characterized by elevated levels of cytotoxic CD4⁺ T cell, CD8⁺ T cell, and NK cell populations at the single-cell resolution (Supplementary Fig. S4A and S4B) and, in addition, antigen-presenting dendritic cell subpopulations (CD14⁺/CD163⁺ and CD1C⁺ DCs), which have been previously linked with favorable responses to immunotherapies (58). Samples classified as the IE-ECM subtype (N = 4) are uniquely characterized by elevated levels of mesenchymal stem cell and malignant cell populations of osteoblastic and chondroblastic lineage (Osteoblastic_3, Osteoblastic_4, Osteoblastic_5, Osteoblastic_6, Chondro_proli, Chondro_hyper1, Chondro_hyper2, and Chondro_trans). These populations were previously shown to express diverse transcriptional programs involved in an inflammatory response and, in addition, angiogenesis, myogenesis, oxidative phosphorylation, hypoxia, MYC, TGF-beta, TP53, and KRAS signaling (39). Lastly, the ID subtype (N = 6) is uniquely characterized by elevated levels of a highly proliferative subpopulation of osteoblastic cancer cells (Osteoblastic_1 and Osteoblastic_2).

Next, when comparing clinical outcomes, we observe that the classifier predicted TME subtypes stratify progression-free survival of human patients (N = 185) similar to what was seen in canines (Fig. 5D). Importantly, we find that the associations between primary tumor TME subtypes and progression-free survival outcomes remain significant (IE: HR: 0.13; 95% CI, 0.019–0.91; P value: 0.04, IE-ECM: HR: 0.47; 95% CI, 0.23–0.97; P value: 0.042) even after adjusting for survival differences associated with conventional prognostic variables such as disease stage at diagnosis (HR: 2.49; 95% CI, 1.25–4.97; P value: 0.009) and percentage necrosis postchemotherapy (HR: 0.25; 95% CI, 0.113–0.56; P value <0.001; Fig. 5E; ref. 12). Furthermore, in an independent dataset of human OS patients (GSE33383, N = 53), the inferred TME subtypes are strongly predictive of developing metastasis within 5 years (Fig. 5F; χ² test P value: 0.021). In summary, these results suggest that canines could serve as a useful real-world model for deriving clinically translatable biomarkers for human OS.

In this work, we performed an unsupervised analysis of the TME of canine OS patients to uncover three distinct TME subtypes that are conserved in human OS. To the best of our knowledge, this study represents the largest and most comprehensive comparative analysis of the OS TME in both canines and humans spanning multiple datasets. Importantly, in contrast to prior works, we show that the TME subtypes are predictive of progression-free survival outcomes of both canine and human patients, independent of previously studied prognostic factors such as tumor location and serum alkaline phosphatase levels in the case of canines, and tumor stage and percentage necrosis in the case of humans. Furthermore, we observe a significant correlation between the TME subtype and CIN in both canine and human OS. Chromosomal instability is a hallmark feature of multiple cancers, including OS, and is one of the central players in accelerating tumor evolution, drug resistance, and metastasis (59, 60). In OS, CIN has been shown to arise from global hypomethylation of DNA (61). However, the effect of CIN on the surrounding TME is complex and an ongoing area of investigation. Although CIN can trigger innate immune-mediated clearance of cancer cells via the activation of the cGAS-STING pathway (60), several immune escape mechanisms have been reported (62). One of them is the silencing of genes encoding type 1 IFNs (63–65) downstream of the cGAS-STING pathway. In both canine and human OS samples analyzed in this study, expression of type 1 IFN response genes is lowest in the ID subtype, the one with the highest CIN (see Supplementary Fig. S10A and S10B). This is also recapitulated at the single-cell level, where we see tumor samples belonging to the ID subtype have elevated levels of proliferative cancer cell subpopulations (Osteoblastic_1 and Osteoblastic_2; Supplementary Fig. S4B) shown to lack expression of type 1 IFN response genes (39).

When comparing against recent literature, we find that the TME subtypes uncovered in this work bear the closest resemblance to three of four recently identified pan-cancer TME subtypes shown to be conserved across several different adult cancer types (41). The IE subtype uncovered in this work recapitulates the IE/nonfibrotic subtype identified across multiple adult cancer types and was shown to have the most favorable responses to immune checkpoint blockade therapies irrespective of tissue of tumor origin (41). In both canine and human OS, only a small fraction of tumors (∼5% in canines, ∼10% in humans) belong to this subtype. The majority of immune infiltrated tumors in OS in fact belong to the IE-ECM subtype, which corresponds to the immunosuppressive IE/Fibrotic subtype from pan-cancer analysis of Bagaev and colleagues (41). These findings potentially explain the lack of efficacy of immune checkpoint inhibitors in recent clinical trials of OS despite the presence of an immune infiltrate (13–15), implying that targeting ECM production in the OS TME in combination with immunotherapy may be an effective strategy to improve patient outcomes (66, 67). The ID subtype in OS recapitulates the Depleted subtype observed in other adult cancer types and is linked with intrinsic resistance to immunotherapy (68). The high levels of CIN observed in this subtype coupled with low type 1 interferon gene expression suggest that these tumors could potentially be susceptible to combination immunotherapy treatments that reverse the silencing of type 1 interferon genes (69). However, sustained type 1 IFN signaling has also been shown to promote immune suppression in the TME and resistance to immunotherapy (70). Hence, additional research is required to precisely determine in which context type 1 IFN treatment can enhance the immune responses of patients with immunologically cold tumors. When comparing the TME subtypes with pan-cancer intratumor immune cytotoxicity scores defined previously by Rooney and colleagues (71), we find that they are strongly correlated (Supplementary Fig. S11A and S11B). However, we find that immune cytotoxicity scores are prognostic in human OS, not canine OS (Supplementary Fig. S11C and S11D). These findings suggest that although immune cytotoxicity scores are closely linked with the TME subtypes, careful consideration of species-specific differences is crucial when using these scores as surrogate markers for TME subtypes.

Despite the merits of this study, there are some limitations that should be noted and improved upon in future work. First, the issue of pet owners’ decisions for euthanasia can confound an accurate estimation of the effect of the TME subtype on overall survival outcomes. This complicates the comparative assessment of prognostic biomarkers in canines and humans. However, in contrast to previous comparative studies of OS, which focus on overall survival (11, 72), we focused on DFIs and progression-free survival outcomes of canines and humans, respectively. This alleviates biases associated with the wishes of individual pet owners and enhances the translational value of the DOG² cohort. However, additional prospective validation is still required to fully realize the translational potential of the canine-derived biomarkers uncovered in this work. Second, the majority of data analyzed in this study was derived from bulk transcriptomic profiling of tumor biopsies masking the underlying intratumor heterogeneity present in each patient (Fig. 4B). Further comparative analysis of canine and human OS based on spatial transcriptomics sequencing of tumors, for example, using the NanoString GeoMx platform, is likely to provide more detailed insights into the TME landscape and its impact on patient outcomes. Third, the number of matched primary and metastatic tumor pairs analyzed in this study was limited. Future studies using additional paired primary and metastatic tumor samples from canines and more sophisticated mathematical modeling approaches such as those inspired by optimal transport theory are likely to provide deeper insights into TME evolutionary patterns.

In summary, this work demonstrates the tremendous potential of harnessing data from naturally occurring cancers in canines to model the complexities of the TME in humans and uncover novel biomarkers of treatment resistance and metastasis. Although this work primarily focuses on addressing an unmet clinical need in OS, the comparative analysis approaches outlined herein hold broad relevance across various other cancer types. Consequently, the insights gleaned from the analysis of canine clinical trial data hold the potential to identify new transformative biomarkers and treatment strategies for human patients, ultimately leading to improved clinical outcomes.

No disclosures were reported.

The content of this publication does not necessarily reflect the views or policies of the Department of Health and Human Services, nor does mention of trade names, commercial products, or organizations imply endorsement by the U.S. Government. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

S. Patkar: Conceptualization, data curation, software, formal analysis, supervision, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing. J. Mannheimer: Resources, writing–review and editing. S.A. Harmon: Writing–review and editing. C.J. Ramirez: Resources, writing–review and editing. C.N. Mazcko: Resources, writing–review and editing. P.L. Choyke: Supervision, writing–review and editing. G.T. Brown: Resources, data curation, writing–review and editing. B. Turkbey: Supervision, writing–review and editing. A.K. LeBlanc: Conceptualization, resources, data curation, supervision, writing–original draft, writing–review and editing. J.A. Beck: Conceptualization, resources, data curation, supervision, validation, investigation, visualization, writing–original draft, project administration, writing–review and editing.

Nucleic acid isolation and QA/QC were performed in the NCI CCR Genomics Technology Laboratory, with mRNA sequencing and initial data analysis conducted at the NCI CCR Sequencing Facility, through the Frederick National Laboratory for Cancer Research (FNLCR), Frederick, MD 21701. We would like to acknowledge Elijah Edmondson for assistance with the analysis of digital pathology slides. Finally, we would like to thank the NCI-COTC member institutions, canine clinical trials teams, and the families who enrolled their dogs into COTC-021/022. This work was supported by the Intramural Program of the National Cancer Institute, NIH (Z01-BC006161) and the Intramural Research Programs of the National Center for Advancing Translational Sciences, NIH (Z01-TR000249).