Integrative Multi-OMICs Identifies Therapeutic Response Biomarkers and Confirms Fidelity of Clinically Annotated, Serially Passaged Patient-Derived Xenografts Established from Primary and Metastatic Pediatric and AYA Solid Tumors

Simple Summary

Solid tumors account for ~60% of pediatric, as well as adolescent and young adult (AYA), cancers, and outcomes for patients with these progressive diseases remain poor. This highlights the critical need to develop tumor models from patients with aggressive cancers so that oncogenic signatures can be identified for therapeutic testing. Thus, patient-derived xenografts (PDXs) were established from sarcoma and Wilms tumor patients at diagnosis or following treatment. Overall, the molecular landscape of serially passaged PDXs recapitulated the original tumor based on an integrated multi-OMICS pipeline that cross-validated cancer-associated pathways. Actionable mechanisms of tumor progression were identified. CDK4/6 and BETs were prioritized as biomarkers of therapeutic response for in vivo validation. In osteosarcoma PDXs harboring pertinent molecular signatures, inhibition of CDK4/6 or BETs decreased growth. This systematic approach that links patient disease history to data generated from its corresponding PDX provides a foundation to discover improved therapies for patients with high-risk cancers.

Abstract

Establishment of clinically annotated, molecularly characterized, patient-derived xenografts (PDXs) from treatment-naïve and pretreated patients provides a platform to test precision genomics-guided therapies. An integrated multi-OMICS pipeline was developed to identify cancer-associated pathways and evaluate stability of molecular signatures in a panel of pediatric and AYA PDXs following serial passaging in mice. Original solid tumor samples and their corresponding PDXs were evaluated by whole-genome sequencing, RNA-seq, immunoblotting, pathway enrichment analyses, and the drug–gene interaction database to identify as well as cross-validate actionable targets in patients with sarcomas or Wilms tumors. While some divergence between original tumor and the respective PDX was evident, majority of alterations were not functionally impactful, and oncogenic pathway activation was maintained following serial passaging. CDK4/6 and BETs were prioritized as biomarkers of therapeutic response in osteosarcoma PDXs with pertinent molecular signatures. Inhibition of CDK4/6 or BETs decreased osteosarcoma PDX growth (two-way ANOVA, p < 0.05) confirming mechanistic involvement in growth. Linking patient treatment history with molecular and efficacy data in PDX will provide a strong rationale for targeted therapy and improve our understanding of which therapy is most beneficial in patients at diagnosis and in those already exposed to therapy.

1. Introduction

Pediatric cancers, which account for approximately 1% of new cancer diagnoses each year, remain the leading cause of disease-related mortality in children [1,2]. While solid tumors are quite rare in this patient population, solid tumors represent approximately 60% of all malignant neoplasms observed in children (0–14 years), as well as adolescents and young adults (AYA; 15–39 years) in the United States [3,4]. Pediatric and AYA solid tumors highlighted in our study include osteosarcoma (OS), the most common malignant bone cancer in the pediatric and AYA populations [5]; rhabdomyosarcoma (RMS), the most common soft-tissue sarcoma in children [6]; and Wilms tumor, the most common renal cancer in children [7].

There is a critical need to identify efficacious and safe therapies to benefit pediatric and AYA patients with aggressive solid tumors (metastatic, relapsed, or recurrent cancers), wherein the 5-year survival rate is still only 0–30% [8]. In addition, there is a great desire in the oncology community to ultimately replace chemotherapies that can be overly toxic, mutagenic, and lead to long-term side effects that diminish quality of life in cancer survivors [1,9]. Despite this urgency, development of anti-cancer therapy drugs specifically for children significantly lags behind progress obtained in adults [8]. This is evident from the fact that, in 1980–2017, the US Food and Drug Administration approved only 11 anti-cancer therapy drugs for pediatric purposes [8]. The slow progress towards achieving improved treatment options for high-risk pediatric solid tumors can be, in part, attributed to factors such as the rarity of these diseases, which makes advancing potential therapeutic interventions into clinical trials rather difficult due to low patient numbers for clinical trials [8]. Furthermore, the vast amount of mechanistic and molecular genetic heterogeneity present among pediatric and AYA cancers further complicates development of targeted therapies [8].

Another major hurdle in improving pediatric cancer patients’ care is attributed to the limited availability of reliable biologically relevant preclinical models [10]. While cell lines are accepted as essential preclinical tools in research, they do not always faithfully recapitulate the molecular landscape observed in the primary human tumors [11,12]. There are some cases wherein the PDX model differentiates from the original human tumor. This can be attributed to sampling bias of genetically heterogenous tumors that may inadvertently occur [13,14]. The dissected piece of the tumor for PDX development may not fully represent the molecular landscape of the entire tumor [15], which, ultimately, may cause further clonal evolution and genetic drift compared to the original tumor [15,16,17,18,19]. It has also been reported that PDXs developed from advanced tumors have higher take-rate than those developed from less aggressive tumors (at early stages of disease progression), which may not fully mimic the full range of molecular landscapes present in the original tumor [13,20,21,22,23]. Additionally, since PDXs lack the immune system as well as the vasculature and stromal components from the host, tumors that predominately rely on stromal components for growth may not provide good fidelity between the PDX and the original tumor [15].

Despite some PDX models not truly recapitulating the molecular characteristics of the original human tumor, overall concordance between key oncogenes in PDXs and their respective original human tumor are evident in the majority of PDXs shown here, as well as in other comprehensive, large-scale genomic characterization studies conducted in PDXs encompassing a variety of solid tumors, such as carcinomas, sarcomas, and melanomas [24]. In many instances, these models can be established from limited tumor tissue through serial passaging in immunodeficient mice, providing an opportunity to interrogate tumor-specific molecular signatures and test new therapies precisely designed for that tumor [11,12,25]. Therefore, molecular characterization of the original human tumor specimens (P0) with their derived PDXs is an essential step in establishing model validity [12]. In the context of pediatric and AYA cancer PDX models, Rokita et al. [26] identified mutational signatures within the PDXs that recapitulate key alterations found in pediatric CNS tumors, extracranial solid tumors, and hematologic malignancies [26]. Furthermore, Nanni and colleagues proved that serial passaging of pediatric bone sarcomas, such as OS, recapitulated the morphology and the RNA expression profile of the original patient tumor [27].

While advancements in next-generation sequencing (NGS) have shed light on identification of targeted therapies based on genomic alterations of the DNA and/or RNA, not all observed genomic alterations will result in consequential outcomes in terms of protein function [28] relevant to cancer pathogenesis. Indeed, in recent clinical studies with molecularly targeted therapies against important targets, such as AKT, pre-treatment levels of the activation/phosphorylation levels of AKT itself predicted the response; however, genomic markers that infer AKT pathway activation such as mutations in PI3K, PTEN, and AKT1,2 did not [29].

Targeted therapies identified by single-level -OMICS analyses at the DNA, RNA, or protein/phosphoprotein level have the potential to miss therapeutic options that could be effective in heterogenous cancers like pediatric and AYA solid tumors. Therefore, to fully understand the molecular stability of PDXs following serial passaging and to exploit the power of using novel clinically annotated PDX models to uncover therapeutic options, we developed a multi-OMICs pipeline encompassing DNA, RNA, protein, and pathway-enrichment analysis. Therapeutic response signatures/biomarkers CDK4/6 and bromodomain extraterminal domain proteins (BETs) were selected as proof-of-concept targets to explore mechanisms of tumor growth and validate anti-tumor response to small-molecule-inhibitor monotherapy in OS PDXs. This will set the foundation for obtaining further insights into the molecular mechanisms that contribute to the pathogenesis of these diseases.

2. Materials and Methods

2.1. NOD.Cg-Prkdc Scid Il2rgtm1Wjl/SzJ (NSG) Mice

NSG mice were obtained from the on-site breeding colony maintained by the Preclinical Modeling and Therapeutics Core (PMTC) at Indiana University Simon Comprehensive Cancer Center (IUSCCC). All procedures were approved by the Institutional Animal Care (IACUC) and Use Committee at the Indiana University School of Medicine (IUSM, IACUC studies #19052 and 22028). Animals were maintained under pathogen-free conditions and maintained on a Teklad Lab Animal Diet (TD 2014, Harlan Laboratories, Indianapolis, IN, USA) with ad libitum access to sterile tap water under a 12 h light-dark cycle at 22–24 °C.

2.2. Development of PDXs from OS, RMS, or Wilms Tumor Specimens

Tumor samples for development of the PDXs were obtained from patients following consent under IUSM IRB protocol #1501467439. The PDXs were developed as previously described by Mattar et al. [30] with some modifications. The original tumor specimens from a patient (P0) were cut into fragments that were either flash-frozen (5 × 5 mm²), cryopreserved (3 × 3 mm²), or implanted (2 × 2 mm²) into the flanks of NSG mice gender-matched to that of the patient to develop PDX passage 1 (P1). Following implantation, growth kinetics were tracked by electronic caliper connected to StudyLog software (South San Francisco, CA, USA). Briefly, the mice were monitored at least twice weekly for overall health and more frequently if clinical signs were noted. During expansion phase, tumor volumes were measured twice weekly by electronic caliper for up to 300 days or until tumor size was 1000–1500 mm³. Tumor volumes in this range were preferred to minimize necrotic tissue in the tumors. Immediately after P1 PDX extraction, the tumors were cut to 2 × 2 mm² fragments under aseptic conditions, and then 1–2 PDX tumor tissue fragments were implanted into the flank of 6–8-week-old gender-matched NSG mice to generate PDX passage 2 (P2) tumors. The same process was conducted to obtain PDX passage 3 (P3) tissue. All P0 and corresponding PDXs were authenticated for their identity/species by DNA fingerprinting analysis using a 16- marker short-tandem repeat (STR) analysis (IDEXX BioResearch, Columbia, MO, USA). The first priority for tumor samples at our institution is for diagnosis and specimens required for ongoing clinical trials. Generation of PDX is the next priority. Tissue quantity can vary greatly from patient to patient. The PDXs characterized in this study are a subset of our larger pediatric and AYA PDX panel of rare tumors. Over the past several years, our success rate was defined as expansion of PDX at least through passage 2. In OS, 51.4% cases resulted in successful PDX generation. For RMS cases, 53.8% led to successful PDX generation, and 80% of Wilms tumor cases gave rise to successful Wilms tumor PDXs. For the biopsy samples used for PDX generation, the surgeon’s expertise and tumor scans were used to identify viable tumor tissue. For resection samples used for PDX generation, the pathologists release viable tumor based on their expert analysis of the tissue sent by the surgical team. To transfer the tissue sample to the lab for further processing, MACS tissue storage medium was used (Miltenyi Biotec, cat # 130–1100-008).

2.3. DNA Extraction

Total genomic DNA was isolated from flash-frozen P0 and PDX samples (n = 30 samples) using the QIAamp DNA Mini Kit (Qiagen, Germantown, MD, USA) or DNeasy Blood and Tissue Kit (Qiagen, Germantown, MD, USA), as indicated by the manufacturer’s instructions. Subsequently, DNA quantity and was assessed by nanodrop, QUBIT Fluorometer, and Agilent 2100 Bioanalyzer.

2.4. Bioinformatic Data Analysis

2.4.1. Somatic Whole-Genome Sequencing (WGS) Analysis

WGS was performed at New York Genome Center (NYGC, New York, NY, USA) on the following PDXs: HT72, HT77, HT87, HT96, HT74, and HT98. For these samples, TruSeq PCR-free preparation was done. WGS on HT120 and HT139 was performed at IUSM Center for Medical Genomics core using the Illumina Nextera DNA Flex Library Prep Kit. For all P0 specimens and PDX samples, the sequencing depth/coverage was 90×. All SNV and CNV data will be uploaded to the dbGaP database (https://www.ncbi.nlm.nih.gov/gap/). OS, RMS, and Wilms tumor cohorts were sequenced using the NovaSeq 6000 Illumina platform to generate raw reads in FASTQ format. The raw sequencing data for the control sample, NA12878, were obtained from the International Genome Sample Resource and used as the normal reference human genome. Adapter sequences, primers, and other unwanted sequences were first identified and trimmed using Cutadapt [31] and later filtered against mouse genome using BBsplit [32] from BBtools suite to remove mouse-originated reads. The raw sequencing reads were aligned to the human reference genome hg38 with Burrows–Wheeler aligner (BWA) [33].

The objective of the bioinformatics analyses was to compare the copy number variations (CNVs) and single-nucleotide variations (SNVs) of each P0 versus respective PDX passages as well as to compare within each set of PDX passages.

2.4.2. Copy Number Variation (CNV) Analysis

CNV analysis was performed using CNVkit [34]. The copy number profile of NA12878 was used as the reference to correct the biases in regional coverage and GC content in P0 and its respective PDX passages of HT72, HT77, HT87, HT96, HT74, HT98, HT120, and HT139. Using the circular binary segmentation, genome-wide segmented copy number estimates for each sample were generated based on fixed thresholds. Genes with significant deletions and duplications in each call set were determined using the “genemetrics” command with default thresholds on log2 ratios and minimum number probes per segment. Within each P0–PDX cohort with NA12878 as the reference, the genes with differential copy number estimates in at least one of the PDX passages were called. To identify cancer-associated genes, 723 genes from the Catalogue of Somatic Mutations in Cancer (COSMIC) database [35] were used to further filter the copy number call set of each PDX cohort. R packages such as “ape” was used to plot the dendrogram to understand the genomes of PDXs, “pheatmap” to generate heatmaps of differential copy number profiles for each cohort and “ggplot” to visualize the overview of genome-wide copy number alteration distribution on the chromosomal level. Correlation bar plots were created to visualize the P0–PDX cohorts. A copy number of 3–5 was designated as copy number gain, whereas a copy number ≥ 6 was referred to as copy number amplification [36]. Anything fewer than two copies was designated as copy number loss or deletion.

2.4.3. Single-Nucleotide Variation (SNV) Analysis

Sentieon Suite was used for data analysis to generate the SNV/Indels calls in each P0–PDX cohort for HT72, HT77, HT87, HT96, HT74, HT98, HT120, and HT139 [37]. The SNVs identified by the variant calling process, and their variant consequences were annotated using Ensembl Variant Effect Predictor (VEP) with loaded genomic prediction algorithms and databases, such as Ensembl’s Variant Consequences, SIFT, PolyPhen, dbSNP, ClinVar, and COSMIC. For SIFT, if the probability score is <0.05, then the variant is predicted to alter protein function and score of >0.05 is predicted to be tolerated and not impact protein function [38,39]. Polyphen scores > 0.908 for variants are characterized as “probably damaging” [39,40]. Variants with scores between 0.556–0.908 are designated as “possibly damaging [40],” and scores ≤ 0.446 suggests that those variants are “benign”, indicating no impact on protein function [39]. Variant haplotypes were called using NA12878 as the reference, which included any variants in P0 or the corresponding PDX passages in each tumor. They were verified by their presence of nucleotide mutations and/or significant shifts in variant allele frequencies, such as losses of heterozygosity (LOH), when compared to NA12878. We further filtered the variants according to the calls’ quality and the impact on function. Therefore, only the variants of high quality (depth ≥ 10, genotype quality ≥ 20, allelic balance between 0.2–0.8) that are predicted to be of high impact and to have a deleterious effect on functions (polyphen score > 0.15 and SIFT ≤ 0.05) were retained [41]. The mutations in genes involved in important molecular processes, as well as oncogenic pathways, were determined. An in-house R package was developed to generate oncoprint plots to visualize a differential mutation profile in tumor samples compared to NA12878. MAFtools and ggplot2 were used to summarize the SNV statistics, i.e., count of variants per variant class, consequences, and transitions/transversions.

2.4.4. RNA-Seq Data Analysis

RNA extraction and paired-end RNA-seq sequencing were conducted at Quick Biology, Inc (Monrovia, CA USA; Illumina sequencing of 30M, PE150). Human Adipose Tissue Poly A+ RNA (Takara Cat# 636162) was used as the filtering control. The sequencing reads were first mapped to the latest UCSC transcript set using Bowtie2 version 2.1.0 [42], and the gene expression level was estimated using RSEM v1.2.15 [43]. The BBsplit function in the BBmap tool was utilized to remove the mouse reads. Differentially expressed genes were identified using the edgeR program [44]. Genes showing altered expression with p < 0.05 and more than 1.5-fold changes were considered differentially expressed. Pearson correlation analysis was performed to study the similarity between the somatic sample, P0, and the PDXs, for each dataset. Goseq [45] was used to perform the gene ontology (GO) enrichment analysis, and Kobas was used to perform the pathway analysis [46]. Additionally, two other pathway databases (KEGG and Reactome) were also utilized for enrichment analysis.

2.5. Western Blot Analysis

Protein extraction of PDX tumor samples was conducted using either UREA or radioimmunoprecipitation assay (RIPA), based on reports suggesting that protein extraction from solid tumor samples require different methods (i.e., UREA and/or RIPA) due to differences in the cytoplasmic location, solubility, and molecular weights of the protein of interest [47]. Due to the limited availability of P0 samples, they were not used for Western blot analysis. Flash-frozen PDX tumor samples were pulverized by smashing twice with a cryo-PREP^® Automated Dry Pulverizer (Covaris, Woburn, MA, USA). Approximately 25 mg PDX powder was lysed in 300 ul 8 M Urea Buffer (8 M urea, 100 mM Tris pH 8, 5 mM DTT) or 300 ul RIPA buffer containing protease and phosphatase inhibitors and processed with a Bioruptor ^® Plus sonication device on high setting for 5 cycles [30 s on; 30 s off] at 4 °C (Diagenode, Denville, NJ, USA). Debris was spun out, and the protein concentration of the cleared lysate was determined by RCDC Protein Assay (Bio-Rad, Hercules, CA, USA) and quantified using BioTek Synergy H4 (BioTek, Winooski, VT, USA). For Western blot analysis, proteins were separated on TGX Stain-FreeTM gels (Bio-Rad, Hercules, CA, USA) along with Precision Plus ProteinTM All-Blue Standards (Bio-Rad, Hercules, CA, USA) and transferred to LF PVDF membrane using the Trans-Blot Turbo Transfer System (Bio-Rad, Hercules, CA, USA). Membranes were blocked for 1 h at room temperature in 5% non-fat dry milk in TBS-T (137 mM NaCl, 20 mM Tris, and 0.1% Tween 20). Membranes were washed following antibody incubation with TBS-T for a total of three 10 min washes. The correct molecular weight for each protein was confirmed by the Precision Plus All Blue Standard (Bio-Rad, Hercules, CA, USA). The following antibodies were diluted in either 5% non-fat dry milk or 5% BSA in TBS-T per manufacturer’s instruction and used for detection: rabbit anti-RAD21 (130 kDa, cat# 4321, Cell Signaling Technology, Boston, MA, USA); mouse anti-c-MYC [9E10] (67 kDa, cat# sc-40, Santa Cruz Biotechnology, Inc., Dallas, TX, USA); mouse anti-p53 [D0-1] (53 kDa, cat# sc-126, Santa Cruz Biotechnology Inc.); mouse anti-Cyclin D3 (31 kDa, cat# 2936, Cell Signaling Technology, Boston, MA, USA); mouse anti-Cyclin E1 (48 kDa, cat# 4129, Cell Signaling Technology, Boston, MA, USA); rabbit anti-CDKN2A/p16^INK4A (17kDa, cat# ab108349, Abcam, Waltham, MA, USA); rabbit anti-PTEN (54 kDa, cat# 9559, Cell Signaling Technology, Boston, MA, USA); rabbit anti-RAC1 (21 kDa, cat# 4651, Cell Signaling Technology, Boston, MA, USA); rabbit anti-HSP90Aβ1 (96 kDa, cat# NBP2-68937, Novus Biologicals, Centennial, CO, USA); anti-phospho-H2AX serine 139 [γH2AX-Ser139] (15 kDa, cat# 2577, Cell Signaling Technology, Boston, MA, USA); rabbit anti-total H2A.X [D17A3] (15 kDa, cat# 7631, Cell Signaling Technology, Boston, MA, USA); rabbit anti-RB1 [D20] (110 kDa, cat# 9313, Cell Signaling Technology, Boston, MA, USA); rabbit anti-phospho-RB1 Ser795 (110 kDa, cat# 9301, Cell Signaling Technology, Boston, MA, USA); rabbit anti-BRD4 [E2A7X] (200 kDa, cat# 13440, Cell Signaling Technology, Boston, MA, USA); rabbit anti-CDK4 [D9G3E] (30 kDa, cat# 12790, Cell Signaling Technology, Boston, MA, USA); rabbit anti-CDK6 [D4S8S] (36 kDa, cat# 13331, Cell Signaling Technology, Boston, MA, USA); rabbit anti-vinculin [E1E9V] (124 kDa, cat# 13901, Cell Signaling Technology, Boston, MA, USA), and rabbit anti-GAPDH [14C10] (37 kDa, cat# 2118, Cell Signaling Technology, Boston, MA, USA). Blots were then incubated with the appropriate horseradish peroxidase-conjugated secondary antibody, diluted 1:5000 in 5% non-fat dry milk–TBST for 1 h at room temperature (Anti-mouse IgG HRP-conjugate, cat# W4021, Promega, Madison, WI, USA; Anti-rabbit IgG HRP-conjugate, cat# W4011, Promega). Membranes were again washed following secondary antibody incubation with TBS-T for a total of three 10 min washes. Proteins were detected using SuperSignal Western Chemiluminescent Substrate (Thermo Fisher Scientific; Waltham, MA, USA) and imaged using the Bio-Rad ChemiDoc Imaging System (Bio-Rad, Hercules, CA, USA). Human vertebral mesenchymal stem cells (HVMSCs; CP7550 ScienCell Research Laboratories, Carlsbad, CA, USA), as well as human osteoblasts-femoral (HOB-F; CP4610 ScienCell Research Laboratories, Carlsbad, CA, USA), were utilized as our controls.

2.6. Histologic Characterization of Original Tumor Specimen and PDX

For analysis of tissue integrity, tumor tissues from P0 and their corresponding passaged PDXs were fixed in 10% neutral buffered formalin at 4 °C for 72 h and then embedded in paraffin. As previously described by our group, 5 µm sections were cut and flattened on a heated water bath, floated onto microscope slides, and dried as previously described [5]. For routine assessment, slides were stained with hematoxylin and eosin [5]. Microscopic analysis of hematoxylin and eosin (H&E) stains was conducted in a blind fashion by a board-certified pathologist. Slides were scanned using Aperio whole-slide imaging software (Sausalito, CA, USA).

2.7. Protein Pathway Activation Mapping via Reverse Phase Protein Array (RPPA) Analysis of FFPE Xenograft Samples

2.7.1. Xenograft Tissue Processing and Generation of Whole-Tissue Lysates

FFPE-embedded xenograft tissue blocks were sectioned at 8 μm onto glass slides. RPPA analysis was conducted on original P0 tumors when available and their respective PDXs for HT72, HT77, HT87, HT96, HT74, HT98, HT120, and HT139. Between 2 and 6 sections per sample were de-paraffinized in xylene, rehydrated in 100% and 95% ethanol, briefly fixed in 70% ethanol containing protease inhibitors (Complete Mini EDTA-free, Roche, Indianapolis, IN, USA), then dehydrated in 95% ethanol, 100% ethanol and then lysed directly from slides in an appropriate volume of extraction buffer containing 50% Tissue Protein Extraction Reagent (T-PER, Thermo Fisher Scientific; Waltham, MA, USA), 47.5% 2× Tris-Glycine SDS sample buffer (Invitrogen, Waltham, MA, USA), and 2.5% β-mercaptoethanol (Thermo Fisher Scientific, Waltham, MA, USA). The resulting whole-tissue lysates were heated for 8 min at 100 °C. Samples were stored at −80 °C and reheated at 100 °C for 2 min, vortexed, and briefly centrifuged just prior to printing.

2.7.2. Array Printing and Analysis

The total protein concentration in each sample was estimated by printing onto a nitrocellulose slide (Avid, Grace Biolabs, Bend, OR, USA) along with a BSA standard concentration curve using a Quanterix 2470 Arrayer equipped with 180 μm pins at a deposition volume of ~11 nL per sample. Total protein levels were assessed by staining with Sypro Ruby Protein Blot Stain (Invitrogen) according to the manufacturer’s instructions. Tissue lysate concentrations were adjusted to a maximum of 250 μg/mL with extraction buffer. RPPA printing for antibody staining was performed as described above, and antibody staining and image analysis of xenograft samples was conducted as previously described [48,49,50,51]. Briefly, signaling pathway activation was evaluated by staining the arrays with 123 antibodies against signaling endpoints, mainly phosphorylated and cleaved protein products. Before use for RPPA analysis, antibody specificity was confirmed by Western blot and analysis as previously described [51]. Antibody staining intensities were generated using the MicroVigene v5 Software Package (Vigenetech; Carlisle, MA, USA).

2.7.3. RPPA Statistical Analysis

Each data value represents the average of 3 replicate printed spots. All final data values were normalized to total protein following the subtraction of negative control. Additionally, all data values within an endpoint, such as a specific kinase, can be compared to each other, but data values between different endpoints such as different kinases or signaling proteins cannot be directly compared. If a zero value is obtained, then those are true zero values where spot intensities were either not above local background or not above negative-control spot intensities. To ensure data reliability, we prioritized RPPA data wherein coefficients of variation (CV) among the replicate spots were ≤10% [52].

2.8. Compounds

Palbociclib (CDK4/6 inhibitor) and OTX015 (BET inhibitor) were both purchased from Chemietek (Indianapolis, IN, USA). For in vivo studies, palbociclib was formulated in 50 mmol/L sodium lactate, pH 4, and OTX015 was formulated in 10% DMSO, 20% PEG400, 5% Tween 80, and 65% water.

2.9. Screening of OS PDX with Small-Molecule Inhibitors Based on Therapeutic Response Biomarkers

Through the multi-OMICS analytical pipeline described in this study, high-risk molecular signatures in OS PDXs derived from patients with aggressive disease were prioritized based on our interrogation of the Drug Gene Interaction Database (DGIdb; https://www.dgidb.org/ (accessed on 7 July 2022)), as well as evidence from the literature for in vivo screening. For these in vivo efficacy studies, HT77and HT96 OS PDX tumor fragments were implanted into gender-matched NSG mice. Once tumor volumes reached 100–200 mm³, mice were randomized and treated with either 120 mg/kg palbociclib or 25 mg/kg OTX015. HT77 PDX mice were treated with vehicle (n = 8) or palbociclib (n = 7) at 120 mg/kg for 5-consecutive days PO, SID, daily for 3 weeks. HT96 PDX mice were treated for ~6 weeks with 25 mg/kg OTX015 (n = 5) or vehicle control (n = 6), BID, PO, daily for 5 consecutive days. For all in vivo studies, tumor volumes and body weights were evaluated 2–3 times per week.

2.10. Statistical Analysis

All data are reported as mean ± SEM. Data were considered significant at p < 0.05. For in vivo studies, tumor growth data was evaluated by two-way ANOVA with a Holm–Sidak multiple comparison test where * p < 0.05, # p <0.001, ^ p < 0.0001. GraphPad Prism Software (GraphPad Inc., San Diego, CA, USA) was used for data analysis and visualization.

3. Results

3.1. Comprehensive Clinical Annotations of PDX Models Established from OS, RMS, and Wilms Tumor Patients at Different Phases of Therapy

The workflow of the PDX pipeline consists of the direct implant of the original tumor samples into NSG mice, followed by serial passaging of the samples to expand PDX tissue for study (graphical abstract). This is followed by STR authentication and multi-OMICS analyses of P0 versus the corresponding PDX passages (P1, P2, and/or P3 when available). Cross-validation of -OMICS data such as CNVs at the protein level provides further confirmation of the molecular signatures. Importantly, this pipeline links the patient disease journey prior to and after sample collection to -OMICS data, integrated pathway analysis, and actionable signatures with target-specific drugs through the Drug Gene Interaction Database. This provides a solid platform to prioritize targeted therapies and evaluate in vivo efficacy and safety. In partnership with the Pediatric Precision Genomics Program at Riley Children’s Hospital, we established a panel of pediatric and AYA solid-tumor PDXs from tumor specimens of naïve and pre-treated patients diagnosed with OS, RMS, or Wilms tumor. Comprehensive clinical annotations for each patient including diagnoses, demographics, and disease status are shown in

Table 1. Overview of patient characteristics for the PDX panel.

PDXs	Diagnosis	Primary or Progressive Sample	Biopsy or Resection	Gender	Race	Age	Disease Status
HT72	OS	Progressive	Biopsy	Male	Caucasian	18	Deceased
HT77	OS	Progressive	Resection	Male	Caucasian	18	Deceased
HT87	OS	Primary	Biopsy	Female	Caucasian	17	High-grade
HT96	OS	Primary	Biopsy	Male	Caucasian	9	Deceased
HT74	RMS	Primary	Biopsy	Female	Caucasian	14	Deceased
HT98	Wilms tumor	Primary	Biopsy	Male	Caucasian	3	C.R.
HT120	Wilms tumor	Progressive	Resection	Female	Caucasian	9	Deceased
HT139	Wilms tumor	Primary	Resection	Male	Caucasian	8	C.R.

C.R. = complete response.

. Details of treatment plans and therapeutic responses prior to and following the acquisition of tumor samples for PDX development are provided in Table S1.

HT72 and HT77 PDXs were derived from metastatic sites in a heavily pre-treated male OS patient at different time points during salvage therapy. HT87 PDX was established from a treatment-naïve female OS patient that was presented with localized disease classified histologically as high-grade OS, received MAP therapy (methotrexate, doxorubicin, and cisplatin), and at the time of publication was off therapy. HT96 PDX was derived from a localized tumor site from a treatment-naïve male OS patient that presented with metastatic disease at the time of diagnosis. After the sample was acquired for PDX development, the patient underwent frontline and multiple salvage therapies but eventually succumbed to disease and is now deceased. HT74 PDX was generated prior to therapy from a female patient with an unclassified high-grade RMS with widely metastatic and advanced disease at diagnosis that did not respond to therapy and is deceased. Wilms tumor models were selected from patients with disease classifications ranging from Stages II–IV. Both Wilms tumors (HT98 and HT139) exhibited localized disease; however, complete response to chemotherapy and nephrectomy was achieved in both cases. Wilms tumor HT120 was derived from a metastatic lesion in the liver of a pre-treated female patient that has presented with multiple recurrences of metastatic disease since the PDX was established.

3.2. Maintenance of Histological Fidelity in PDX Models

Sections of the serially passaged PDXs were analyzed by H&E and fidelity evaluated by a blinded pathologist. Overall, the PDX-passaged tumors maintained the histological fidelity of their original tumors (Figure S1). Though, as has been reported by others, human stromal components are lost following engraftment and serial passaging [53]. Once solid tumors engraft and propagate in the murine host, the PDX tumors are comprised of human tumor cells and mouse stromal tissues [54,55,56]. WGS analysis indicated that the percentage of mouse reads was similar in the majority of the PDX passages (8–12%) with one exception. HT87 mouse reads were 21–27% mouse across the HT87 P1–P3 passages (Figure S2). It has been reported that murine stromal content can either be conserved or vary among tumors irrespective of whether they originate from similar or different cancers [57].

3.3. Genome-Wide Analysis of P0 Tumor Specimens versus Their Respective PDX Serial Passages

Genome-wide variant analysis was done to identify and compare the molecular signatures at the DNA level within each P0–PDX cohort. Multiple levels of comparative DNA analysis were conducted using WGS data of each P0–PDX cohort. In the first level of analysis, an overview of the genomic landscape where similarities and differences in DNA sequences of the P0 versus the respective passages and among the PDX passages for each individual model are described irrespective of potential impact on protein expression and function. This is complemented with a focused analysis on to what extent oncogenic-associated genes are differentially impacted by serial PDX passaging.

3.4. Comparative Anlysis of CNV Profiles across the Genome

Circular binary algorithm was used to segment the genomic regions with similar average log2 ratios across all genes for P0 tumor specimens and the corresponding PDXs for HT72, HT77, HT87, HT96, HT74, HT98, HT120, and HT139. NA12878 served as the normal human genome filtering control, as it is a well-recognized benchmark dataset [58]. To evaluate overall genomic stability of serially passaged PDXs, the chromosomal distribution of all CNVs across the genome were visualized by chromplots (Figures S3–S10) [59]. Compared to OS PDXs (HT72, HT77, HT87, HT96) and RMS PDX (HT74), the Wilms tumor PDXs (HT98, HT120, and HT139) had fewer CNV events (losses/deletions or gains/amplifications), consistent with the fact that Wilms tumors tend to be more genomically stable tumors, as illustrated by chromplots (Figures S3–S10) and in the detailed analysis below.

CNV similarities versus differences between the P0 sample and its subsequent passages (P1, P2, and/or P3), as well as among the passages for each model, were quantified by Pearson’s correlation coefficient (PCC) score and plotted on correlation heatmap plots (Figure 1) [60,61]. A Pearson’s positive correlation coefficient (PPCC) score of 1 signifies complete concordance between the datasets and score of 0 indicates complete discordance [17,61]. PPCC scores are classified as weak (0–0.3), moderately positive (0.3–0.5), strongly positive (0.5–0.7), or very strongly positive (0.7–1) [61]. OS PDXs HT72 and HT77 (from same patient at different phases of therapy), as well as HT96 cohorts, exhibited a very strong positive correlation between P0 and passaged PDXs for total CNVs (PPCC scores, 0.9–0.98 for HT72 and HT77; PPCC scores of 0.72–0.94 for HT96 Figure 1A), and COSMIC-associated CNVs (PPCC scores 0.90–0.98 for HT72 and HT77; PPCC sores of 0.75–0.95 for HT96 Figure 2A). While OS PDX HT87 exhibited relatively more discordance than the other OS PDX, a strong positive correlation between P0 and its respective PDX passages was still evident for total CNVs (PPCC scores, 0.57–0.97; Figure 1A) and HT87 COSMIC-associated CNVs (PPCC scores 0.56–0.97; Figure 2A). Whether this discordance is attributable to the lower number of human reads that could be called due to high murine stromal levels (21–27% mouse reads across the passages, Figure S2) and/or clonal divergence is not clear. The highest degree of divergence between the P0 and the PDX passages was found in the RMS HT74 PDX. PPCC scores for the RMS HT74 PDX were moderately positive for the total CNVs (P0 versus P1, PPCC = 0.32; P0 versus P2, PPCC = 0.34; Figure 1B). Moreover, the COSMIC-associated CNVs were even more discordant than the total CNVs presenting weak PPCC scores (P0 versus P1, PPCC = 0.25; P0 versus P2, PPCC = 0.21; Figure 2B), consistent with genetic drift of the serially passaged PDX at least at this level of DNA analysis. Of note, once the sarcoma PDXs were established, CNVs across the passages of all PDXs stabilized. For Wilms tumor PDXs (HT98, HT120, and HT139), very strongly positive correlation was observed between each P0 and respective PDX passages for total CNVs (PPCC scores = 0.94–0.96 for HT98, PPCC scores of 0.91–0.98 for HT120, and PPCC scores of 0.93–0.98 for HT139; Figure 1C). Similarly, COSMIC-associated CNVs also showed strong positive correlation in Wilms tumor samples (PPCC scores = 0.96–0.99 for HT98, PPCC scores of 0.89–0.97 for HT120, and PPCC scores of 0.92–1 HT139; Figure 2C), thus highlighting the relatively low levels of inherent genome instability in these tumors [52].

To provide additional insight into the degree of changes found in the P0–PDX sets, the absolute number of losses and gains/amplifications of the CNVs in each COSMIC-associated gene are presented in Table S2A–H. At this level of analysis, it becomes evident that most differences in copy number losses, as well as gains/amplifications in the COSMIC-associated genes, are modest. As an example, all amplifications (>6 copies) and copy number gains (<6 copies) exhibited only minor differences of 1–2 copies across the P0 and the respective PDX passages (Table S2A–H). Allelic losses were overall maintained with some exceptions highlighted below.

CNV profiles of 25 genes from the COSMIC database were prioritized based on their expression and the role they have in that particular cancer type. For example, across pediatric and AYA OS PDXs, ATRX gene deletions (0–1 copies present) were evident (Figure 3A–D; Table S2A–D). Previous studies have shown that deletions in ATRX resulting in loss of protein expression can contribute to pediatric and AYA OS tumorigenesis [62,63]. OS PDXs also had increased MYC and RAD21 copy number gains (3–7 copies; Figure 3A–D; Table S2A–D), both of which are known cancer markers with multi-functions leading to increased DNA replication and oncogenic-induced replication stress, a hallmark of aggressive OS [5]. Notably, HT72, HT77, and HT87 acquired a deletion of RB1 in one allele during the PDX serial passaging (Figure 3A–C; Table S2A–C), which was not seen in the respective P0 sample, whereas HT96 maintained the RB1 deletion found in the original tumor (Figure 3D; Table S2D). A bi-allelic PTEN deletion was evident in the OS HT72 PDX that was not found in the P0 sample or the OS PDX HT77 derived from the same patient at a later resection (Table S2A,B). It should be noted that not every gene will appear in each PDX cohort. As an example, ATM has a copy number gain present in P0, P1, P2, and P3 of HT96 that was not observed in NA12878 normal control sample that was used as the filter. However, ATM copy number gain, amplification, or deletion in Figure 3A–C (HT72, HT77, and HT87) is not listed, since the copy number change was not different between NA12878 vs. P0, P1, P2, and/or P3. In the unclassified RMS HT74, CCND1 gain (3–5 copies) and a CDKN2A deletion (0 copies) were found in both the P0 and the respective passages, which is a common alteration found in RMS [64,65] (Figure 3E; Table S2E). CCND1 had copy numbers of 3 and 5 in HT74-P1 and HT74-P2, respectively, and both alleles of CDKN2A were deleted in HT74-P1 and HT74-P2 (Figure 3E; Table S2E). In primary encapsulated Wilms tumors (HT98 and HT139), deletion of both alleles of AMER1 was detected in both the P0 and respective PDXs (Figure 3F,H; Table S2F,H). WTX has been reported to be mutated, truncated, or deleted in 30% of Wilms tumor cases [66,67]. However, clinical significance of WTX alterations are still unclear. When assessing the CNV profile of HT120, which was derived from a metastatic Wilms tumor of the liver (Figure 3G; Table S2G), no deletions for WTX were evident. Though CNVs in genes involved with cell cycle regulation were evident, such as CDK4 and CCND2 copy number gains (three copies of each gene) compared to primary Wilms tumor samples [HT98 (Table S2F) and HT139 (Table S2H)]. Of note, copy number gains of KRAS (3 copies) were found in the metastatic Wilms tumor HT120-P0 and corresponding passaged PDXs (Figure 3G, Table S2G). Other groups have reported activating mutations of KRAS in human Wilms tumor samples that correlate with hyperactivation of the WNT pathway [66,68].

3.5. Integration of Prioritized Cancer-Associated CNVs in OS with Corresponding Protein Expression

To validate the impact of copy number gains or amplifications and losses on gene expression, protein levels were analyzed for cancer-associated pathways in the OS PDXs (Figure 4). Due to limited availability of P0 samples, they were prioritized for WGS and RNAseq. Specific CNV-protein associations (Figure 4) were selected due to their importance in OS pathogenesis. RAD21, MYC, γH2AX, total H2AX, and P53 all have been reported to play roles in the DNA replication, replication stress, and genome instability in OS [5,69]. Additionally, amplification of cell cycle proteins cyclin D3 (CCND3) and cyclin E1 (CCNE1), as well as CDKN2A deletion, are common signatures found in high-risk OS [70,71,72], all of which result in uncontrolled cell cycle progression in OS [72]. Moreover, roles for RAC1 [73] and HSP90Aβ1 [74] have also been reported in OS pathogenesis.

Cross validation of CNVs (Table S2) and protein expression (Figure 4) demonstrated congruence of copy numbers and protein expression for CCND3 in the HT72 and HT77 OS PDXs (Figure 4, copy number of 16–20), as well as for CCNE1 in the HT87 PDXs (Figure 4, copy number of 5–6). In some cases, copy number gains/amplifications or losses/deletions (Table S3A,B) do not correlate with the corresponding protein expression (Figure 4). For example, as mentioned above, PTEN in HT72 is reported to have zero copies in both passages (Table S2A). However, PTEN protein was still present in HT72-P1 but absent in HT72-P2, which correlates with the lack of both PTEN alleles (Figure 4). While WGS data for TP53 indicated copy numbers of 3–6 (Table S3A,B), the P53 protein was not detectable (Figure 4). Detailed analysis of the TP53 gene locus revealed that multiple SNVs in the upstream and downstream regions of TP53 gene were present, which likely disrupt the expression of P53. Additionally, the RB1 gene was found to have two wild-type copies in OS HT87-P2 and HT87-P3, but protein expression levels were low relative to other PDX models, highlighting the fact that intact RB1 alleles do not necessarily equate to higher RB1 protein levels. In OS HT96-P1-P3, one allele of RB1 was deleted based on CNV data; however, protein levels were still maintained (Figure 4). BET BRD4 protein was relatively stable across the OS PDX passages. Moreover, while CDK4 protein levels were maintained across the OS PDX passages, CDK6 levels varied. In the paired HT72 and HT77 OS PDXs, CDK6 was quite high in HT72 but absent in the HT77 PDXs. Furthermore, CDK6 levels in OS PDX HT87 and OS PDX HT96 decreased during the serial passaging, suggesting that CDK6 is not necessary for growth and survival of the PDX. Differential expression of CDK4 and CDK6 was evident in the OS PDXs. Dai et al. previously demonstrated that CDK4 and CDK6 can have distinct biological functions that regulate cancer growth and metastasis [75]. Functional diversity of CDK4 and CDK6 have not been investigated in sarcomas and warrants further investigation. Notably, these disparities between linking CNVs and corresponding protein expression indicate the need to validate protein expression since the efficiency of sequencing reads, coupled with factors regulating transcription, translation, and post-translation, could lead to inconsistencies between CNV and protein levels [76,77].

3.6. Comparative Analysis of Single-Nucleotide Variations (SNVs), Somatic Variants, and Single-Base-Pair Substitution Variants in P0 Tumor Specimens versus Their Respective PDXs

Quantitative data with descriptive analysis [78] on the characterization of common genomic features (SNVs, somatic variants, and transitions/transversions) was conducted. In this analysis, genomic deviations are depicted in the P0 and the corresponding PDX passages. While some alterations were predicted to be pathogenic, many deviations were found to either not occur in cancer-associated genes or to occur in genes with unknown or not fully characterized functions.

Stability of small (SNVs, Indels) and large (insertions, deletions) structural variants, as well as somatic base substitutions for COSMIC genes, were overall preserved between the P0 tumor and the corresponding passaged PDX cohorts (Figure 5).

Scheme 10 base pairs of DNA [79]. At the total SNV level, all sarcoma PDXs exhibited a very strongly positive correlation between P0 and respective passaged PDXs and among the PDX passages for each individual OS model (PPCC scores of, 0.8–1 obtained HT72, HT77, and HT87 cohorts; PPCC sore of 0.8–0.9 for HT96; Figure 6A). PPCC scores for the RMS HT74 PDX were strongly positive for SNVs (P1 versus P2, PPCC = 0.94; Figure 6B). However, P0 versus P1 or P2 had weak PPCC scores of 0.10 and 0.17, respectively; Figure 6B). For Wilms tumor PDXs (HT98, HT120, and HT139), very strongly positive correlation was observed for SNVs between each P0 and respective PDX passages for total CNVs (PPCC scores = 0.96–0.97 for HT98, PPCC scores of 0.92–0.94 for HT120, and PPCC scores of 0.94–0.95 for HT139; Figure 6C).

In P0 tumors and their respective passaged PDXs, most of the SNVs in P0 (Figure 7) were preserved in successive PDX passages. UpSet plots demonstrate intersections of SNVs that are the same between the P0 and its subsequent PDX passages. For example, in HT72, 87.4% of SNVs are common between P0, P1, and P2 (Figure 7A), whereas HT77 illustrates 88.9% common SNVs among P0, P1, P2, and P3 (Figure 7B). In HT87 and HT96, 78.4% (Figure 7C) and 76.8% (Figure 7D) SNVs, respectively, are commonly found between P0 and their subsequent passages. Similarly, RMS PDX HT74 also has about 76.9% SNVs occurring between the P0, P1, and P2 (Figure 7E). For Wilms tumors, SNVs that are commonly occurring between P0 and their respective passages are all > 90% (HT98 = 94.7%, HT120 = 91.9%, HT139 = 96.3%; Figure 7F,H). In addition, there were small numbers of SNVs that were only found in subsets of a particular cohort (Figure 7).

3.7. Single-Base-Pair Substitutions Patterns across PDX Passages

This variant type includes nucleotide transitions and transversions [81]. Transitions involve the exchange of a two-ring purine with another two-ring purine [i.e., adenine (A) with guanine (G) and vice versa] or one-ring pyrimidine with another one-ring pyrimidine [thymine (T) with cytosine (C) and the reverse] [81]. Transversions occur when a two-ring purine is interchanged with a one-ring pyrimidine or vice-versa [81]. Generally, transitions are less likely to alter the amino acid sequences, but they occur at higher frequency then transversions [81,82]. Our data was consistent with previous studies in colorectal cancers [83].

Across the P0 tumors and their corresponding PDX cohorts, the number of transversion events (Figure 8A–H) was lower than the transition events (Figure 8I–L). Overall frequency of transversions and transitions in all genes remained stable in the serial PDX passages compared to its respective P0 tumor (Figure 8) [78]. In terms of maintaining the same transversions or transitions within the same genes, 52–88% of the transversion and 55–88% of the transitions were maintained in OS P0 tumor and its corresponding OS PDX passages (Table S3A,B). Up to 49% of transversions and 51% of transitions from P0 were preserved between RMS P0 and its respective passaged PDXs (Table S3A,B). The Wilms tumor had a maximum of 94% and 97% of transversions and transitions respectively that were conserved in P0 tumor and its subsequent PDX passages (Figure 8; Table S3A,B).

3.8. Potential Functional Impact of Somatic Variants in the P0 Tumor and Respective Passaged PDXs

Potential functional effects of the somatic variants were assessed between the P0 tumor and the corresponding PDX passages (Figure 9). Consequences of these variants were predicted based on the ensemble plus sequence ontology predictions and organized in order of severity in terms of predicted effect on transcription (Figure 9) [84]. Splice acceptor variants (Figure 9A), stop-gains (Figure 9B), start losses (Figure 9C), and inframe deletions (Figure 9D) rarely occurred in the P0–PDX cohorts. In contrast, more frequent occurrences of in-frame insertions (Figure 9E), frameshift variants (Figure 9F), and missense variants (Figure 9G) were found. Notably, upstream and downstream gene variants were by far the most common types of variant consequences observed in the PDX models (Figure 9H,I). In the sarcomas, 50–82% of the variants in the P0 were maintained in the corresponding PDXs (Table S4). In Wilms tumors, high fidelity of variants was observed in the P0 and corresponding PDX passages with 91–96% of the variants retained (Table S4). Moreover, approximately 50% of the frameshift variants in OS, 60% in RMS, and 53% in Wilms tumor were passed from P0 to their corresponding PDX passages (Figure S11). As for missense variants, a high degree of variability was observed among the P0–PDX cohorts with 21–84% retention frequency of these variants in the P0 and all the corresponding PDXs (Figure S12), with the lowest retention frequency in the OS HT96 PDX and the highest in the Wilms tumor HT98 PDX.

3.9. Functional Predictions of Somatic Variants from the P0 Tumor and the Respective Serially Passaged PDXs

Somatic variants from the P0 tumor and the respective passaged PDXs were filtered against the reference NA12878 dataset. Thus, only variants found within the original tumor and/or the respective passaged PDXs were retained and evaluated by sorting tolerant from intolerant (SIFT) algorithm and [polymorphism phenotyping (PolyPhen) algorithm to predict the functional outcome of these variants [38,85].

For the PDX datasets, genes with variants considered to be damaging or deleterious based on SIFT and/or Polyphen scores are presented (Figure 10; Table S5). Variant allele frequencies (VAFs) and oncoplots are presented for each variant when the frequency was 0.20 across in at least one of the samples (P0 tumor and/or the respective passages). While oncoplots of each P0 tumor and the respective passaged PDXs depict deleterious missense mutations in many genes, it is noteworthy that deleterious missense mutations in MUC16 are present among all of the P0 tumors and their respective passaged PDXs (Figure 10A–H). A repeating peptide epitope of MUC16, a mucosal protein, known as cancer antigen 125 (CA125) is overexpressed in adult cancers, such as ovarian, lung, and pancreatic cancer, where it contributes to immunoprotection, metastasis, oncogenic signaling, and disease [86,87]. Missense variants in MUC16 are one of the most common mutations observed in this gene [88,89]. However, two deleterious single nucleotide polymorphisms observed in all our PDXs cohorts include: rs76810971 and rs769311065 (Table S5), and further investigation needs to be conducted to elucidate the roles of these polymorphisms, especially in the context of pediatric and AYA OS, RMS, and Wilms tumors. Additionally, in the Wilms tumor PDX HT139 oncoplot (Figure 10H), there is a “start-loss” variant for MYCL where the normal “A” of start codon “AUG” is switched with a “G” so the proper translation is not initiated [90]. Thus, this predicts that no MYCL protein is generated in this Wilms tumor PDX due to translational cessation.

3.10. Transcriptome-Based Enrichment Pathway Analyses in Sarcoma PDXs

Comparison of the sarcoma P0 and the respective PDX-P2 RNA-seq datasets indicated positive correlation of global gene expression with HT72, HT77, and HT96 exhibiting high transcriptome correlation (Figure S13, r = 0.87–0.96). P0 and PDX-P2 RNA-seq datasets derived from HT87 and HT74 were positively correlated but to a lesser extent (Figure S13, r = 0.7–0.78).

Prioritized gene targets of interest from the RNA-seq in P0 and PDX-P2 sarcomas were interrogated using pathway databases (GO, KEGG, and Reactome). The top 15 pathways exhibiting high statistical significance are illustrated with Reactome enrichment bubble plots for HT72 (Figure 11A,B), HT77 (Figure 11C), HT87 (Figure 12), and HT96 (Figure 13). Herein, gene ratio (number of genes in that particular pathway vs. total number DEG genes) and p-values are presented for the P0 tumor, as well as PDX-P2. Notably, genes for all passaged PDXs and original tumor pairs demonstrated high degree of enrichment in the cell cycle and mitotic cell cycle pathways that are commonly dysregulated in these tumor types (Figure 11, Figure 12 and Figure 13, p-values of ≤2 × 10⁻⁵) [91,92,93,94]. Furthermore, irrespective of the type of pediatric and AYA solid tumor present, all original P0 tumors and passage 2 (P2) of the PDXs indicated that cell cycle pathways were highly statistically significant and also had the highest numbers of differentially expressed genes.

Whether the transcript levels correlated with copy number of genes and/or protein expression from Figure 4’s Western blot really depended on the gene interest and which PDX was being investigated (Table S6). For example, in HT72, HT77, and HT96, transcript levels and copy numbers for CCND3 were increased. Additionally, an increase in Cyclin D3 protein expression for HT72-P2, HT77-P2, and HT96-P2 was also evident (Figure 4). For CDK6, increased transcript levels correlated within increased protein expression in HT72-P2 and HT87-P2 (Figure 4, Table S6). However, for other genes, such as TP53, transcript levels did not necessarily correlate with copy numbers and/or protein expression. Such discrepancies may be attributed to post-transcriptional and/or post-translational modifications.

Reactome pathway enrichment data is presented as bubble plots for each primary tumor and subsequent PDX passage 2 of HT72-P0 (Figure 11A; OS), HT72-P2 (Figure 11B; OS), HT77-P2 (Figure 11C; OS), HT87-P0 (Figure 12A; OS), HT87-P2 (Figure 12B; OS), HT96-P0 (Figure 13A; OS), and HT96-P2 (Figure 13B; OS). RNA-seq was conducted on the original tumor specimen (n = 1) and PDX-P2 (n = 3 biological replicates) of each PDX. Human Adipose Tissue Poly A+ RNA was used as the filtering control for the RNA-seq data. Differentially expressed genes (DEGs) with a >1 counts per million, fold change of >1.5, and p value of <0.05 were selected for subsequent pathway enrichment analysis via (GO, KEGG, and Reactome). Since all three databases showed overlapping pathway enrichment results, Reactome was selected as the representative database. The bigger the bubble plot the greater the number of genes that are affected in that specific pathway.

3.11. Comparative Analysis of Functional Protein/Phosphoprotein-Based Cell Signaling Activation Architecture and Pathway Enrichment in P0 Tumor Specimens vs. Their Respective PDX Passages

Different-OMIC datasets can be linked and used to cross-validate oncogenic pathway activation in cancer. The original P0 (when sufficient high-quality tissue was available) and subsequent PDX passages (P1–P3) were used for ≥ RPPA-based quantitative protein and phosphoprotein analysis and pathway activation mapping. RPPA was performed on 123 key proteins and phosphoproteins known to be “actionable” targets and/or associated with oncogenic pathways, with CV of <10% as the cut-off for positivity. Under these parameters, n = 13 total proteins and n = 58 phosphoproteins yielded positive signals of significance (Figure 14 and Figure S14) [95]. It should be noted that, in Figure 14 and Figure S14, while relative levels of proteins can be compared for each individual protein across a row, total/phosphoprotein levels cannot be compared between different proteins down each column. Profiling of select kinase activities via RPPA indicated significant activation of major pathways known to contribute to cancer survival, growth, and metastasis in the P0 and their passaged PDX. Of interest, Notch 1, a transmembrane receptor involved in growth and invasion was detected in all the P0 and respective PDX cohorts (OS and RMS Figure 14 and Wilms tumor, Figure S14) [96,97,98]. In contrast to the P0 tumor, Wilms tumor HT139 PDX passages display relatively decreased levels of the downstream effector for mTORC1, p70S6K, which is phosphorylated at serine 371 (S371) [99], as well as relatively increased cofilin phosphorylated at serine 3, which was previously demonstrated to contribute to cancer metastasis [100,101,102,103] and be associated with human glomerular disease [103]. OS HT77 PDX exhibited increased protein kinase C alpha subunit (PKCα) phosphorylated at T638 and protein kinase C beta II subunit (PKCβII) phosphorylated at T648 compared to P0. This signature is consistent with the invasive nature of the HT77 passaged PDXs; for PKC, it has been reported to promote metastasis in human and canine OS [104,105]. Using a stringent cut-off of CV ≤ 10 [52] as a positive signal for the RPPA data, multiple pathways contributing to autophagy, apoptosis, cell cycle, DNA damage, proliferation, survival, angiogenesis, invasion, migration, and cytoskeletal remodeling are activated, such as CDK4/6, DNA damage/oncogenic stress, MEK-ERK, and PI3K/AKT/mTOR.

We next integrated these data to develop a diagram depicting how pathway activation from the RPPA-based signaling analysis operative in both sarcoma and Wilms tumors can be used to construct biological networks that likely contribute to sarcoma pathogenesis (Figure 15). Many of the proteins are involved in multiple pathways, such as autophagy, apoptosis, cell cycle control, DNA damage, proliferation, survival, angiogenesis, invasion, migration, and cytoskeletal remodeling (Figure 15) [106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132]. Furthermore, Figure 15 integrates oncogenic-associated pathways identified by COSMIC-associated CNVs (Figure 3 and Table S2A–H) and/or proteome analyses. For example, in OS tumors COSMIC-associated CNV data indicated copy number gains/amplifications in MYC, RAD21, Cyclin D3, Cyclin E1, as well as loss of CDKN2A and TP53 genes (Figure 3 and Figure 4, Table S2A–H) that were then validated for corresponding proteins by Western blot (Figure 4) and/or RPPA (Figure 14).

Based on the RRPA data, there is evidence that a number of actionable cancer-associated pathways are indeed activated in all of the PDX samples. However, there are relative differences in the activation level of specific pathways across the PDX panel (Table S7). Whether these subtle differences in pathway activation equate to differential therapeutic responses in each PDX warrants further investigation.

3.12. Proof-of-Concept In Vivo Studies in OS PDX: Monotherapy Screens to Explore Mechanisms of Tumor Growth Based on Therapeutic Response Biomarkers

By utilizing the multi-OMICS analytical pipeline described in this study, we compiled high-risk molecular signatures in the PDX derived from patients with aggressive disease based on our interrogation of the Drug Gene Interaction Database (DGIdb; https://www.dgidb.org/ (accessed on 7 July 2022)), as well as evidence from the literature. Attractive candidates for therapeutic targeting are summarized for HT72 (Table S8A), HT77 (Table S8B), HT87 (Table S8C), HT96 (Table S8D), HT74 (Table S8E), and HT120 (Table S8F) [5,65,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154].Wilms tumor patients from which primary localized tumors were used to develop HT98 and HT139 PDXs responded well to frontline treatment (Table S1). Therefore, linking molecular signatures to therapeutic response biomarkers and drug candidates in these tumors was not prioritized. Based on our multi-OMICS pipeline, CDK4/6 hyperactivation and BETs were prioritized as targets of interest for in vivo screening since both signatures are present in OS PDX derived from patients with aggressive and metastatic disease. CDK4/6 plays an integral role in OS growth, survival, and invasion [5,65]. Targeting CDK4/6 was one priority, since these proteins function as major regulators of the cell cycle pathway, which is activated based on CNVs analysis (Figure 2). This includes CCND3 amplification, lack of CDKN2A/p16^INK4A, and presence of RB1+ where 1–2 wild-type alleles are present and the expression of RB1 (Figure 4). This is further supported by the fact that activation of cell cycle pathways is evident in the transcriptome-based enrichment pathway analyses (Figure 11, Figure 12 and Figure 13), as well as the protein pathway enrichment analysis (Figure 14). For example, in PDX HT77 as well as other OS PDXs, the RPPA data show increased total plus phosphorylated CDK4 (T172), as well as increased expression of phosphorylated RB (S780), which serve as therapeutic response biomarkers and predict sensitivity to CDK4/6 inhibitors such as palbociclib [65]. Additionally, lack of CDKN2A/p16^INK4A and overexpression of other components in this pathway, such as CCND3 (Cyclin D3) and RB1 (Figure 4), also predict sensitivity to CDK4/6 inhibitors as a means to block cell cycle entry [65]. PDX HT77 was derived from a metastatic lesion in a heavily pretreated 18-year-old OS patient (Figure 16A, Table S1) that was established prior to treatment with palbociclib. While the patient was treated late in the disease course with palbociclib, he eventually succumbed to the disease. In PDX HT77, palbociclib was well tolerated based on body weights (Figure S15A) and significantly blocked tumor growth compared to vehicle as long as the mice remained in treatment (Figure 16B). In this study, we were not only able to validate growth dependency on CDK4/6 but also provided evidence that the PDX HT77 was indeed responsive to palbociclib at a point in the disease history when the PDX was established.

Targeting BET proteins was also prioritized, since BET-BRD4 is a therapeutic target of emerging interest in sarcomas. BETs contribute not only to transcriptional regulation by binding to acetylated lysines, but also, in the appropriate molecular contexts, can lead to an imbalance in transcription–replication kinetics and induce non-tolerated levels of replication stress in cancer [155]. Studies by us demonstrated that BET inhibition can lead to increased replication stress and increased death in models of OS treated with CHK1 inhibitors [5]. Moreover, MYC, a common downstream target of BETs, is sensitive to BET inhibitors in some cellular contexts [5,156,157]. As previously mentioned, OS PDX HT96 is derived from a treatment-naïve 9-year-old OS patient that ultimately succumbed to metastatic disease (Figure 17A; Table S1). This model exhibits replication stress (i.e., BRD4 and γH2AX) and expresses MYC (Figure 4). Treatment with BET inhibitor, OTX015, lead to a transient decrease in OS PDX96 growth compared to vehicle control (Figure 17) and was well tolerated, as indicated by stable body weights over time (Figure S15). As mentioned above, molecular signature tables of the PDXs derived from patients with aggressive disease are provided in Table S8A–F. Data obtained from our multi-OMICS pipeline and pathway integration provide numerous biomarkers of therapeutic response and drug candidates of actionable pathways that could be used for development and testing of other targeted therapies, which will be the focus of future investigations.

4. Discussion

Precision genomics has emerged as a critical component of clinical decision making in treatment of patients with aggressive disease. Identifying the optimal targeted therapy for each patient based on the molecular signature of the tumor has the potential to change long-term treatment paradigms in hard-to-treat cancers. The field now needs to use relevant preclinical oncology models that are linked to the patient’s disease journey in a de-identified manner to validate targeted therapy such that earlier intervention with efficacious and safe molecularly guided therapy can be realized.

Therefore, there is great impetus to design feasible, cost-effective, and reliable approaches to uncover and cross-validate therapeutic targets as well as develop preclinical models that recapitulate the original tumor. Our objective for this study was two-fold. First, we sought to develop an integrated multi-OMICS approach for identifying molecular signatures (DNA, RNA, protein, and/or phosphoprotein) in PDX models developed at our institution in which disease history prior to and after tumor collection is linked in a de-identified manner to PDX datasets. In this study, PDX models derived from treatment naïve and pre-treated solid-tumor patients were evaluated. It is important to note that one limitation in our study is that, by chance, all of our PDXs were developed from Caucasian patients (Table 1) during the initial stage of program development at our institution. It has been reported that racial and ethnic diversities can have an impact on tumor mutational burden, as well as on response to immune-based therapies [158,159]. Furthermore, tumor mutational burdens for races, such as African American patients, have been reported to be inflated when publicly available germline variant databases are used for filtering non-somatic polymorphisms, demonstrating the need to use patient-specific germlines [159]. Second, to elucidate if genome instability impacted the molecular stability of original tumors following engraftment and serial propagation, PDXs derived from tumors with inherent genome instability (OS, RMS) and relatively more stable genomes (Wilms tumors) were selected for deep molecular characterization.

In depth surveillance and cross-validation of -OMICS data provides insights in gaining a holistic view of the molecular landscape of the cancer. Furthermore, analysis of pathway activation at the protein expression and protein activation/phosphorylation levels provides insights into prioritization of therapeutic targets [160]. RPPA-based pathway activation mapping is one such approach by which an unbiased multiplexed and comprehensive array of targeted total and phosphoproteins can be evaluated in a limited amount of sample [28,161]. This multiplex analytical method can allow more samples to be tested at one time, which has proven to be more cost-effective compared to other molecular proteome platforms. RPPA-based analysis of individual tumor samples could be a major advantage in the clinical setting, for it has the potential to provide a quantitative read-out of the expression and functional activation state of hundreds of protein drug targets that verify target and pathway activation using slides obtained from paraffin-embedded tissue blocks [28,161]. In our study, total and phosphoprotein data from RPPA analysis was used to highlight activation of a number of key signaling pathways involved in these diseases (Figure 14).

Generation of novel clinically-annotated PDXs provided the opportunity to have sufficient tissue to validate CNV effects on protein expression (Figure 4). Notably, based on the RPPA data, expression of total PTEN and MYC (c-MYC) was the highest in HT77-P3 and HT87-P3 (Figure 14), although Figure 4 illustrates that PTEN and MYC levels were detected in most OS PDXs by Western blot. As previously mentioned, PTEN in HT72 is reported to have zero copies in both passages (Table S2A). However, PTEN protein was still present in HT72-P1 but absent in HT72-P2, which correlates with the lack of both PTEN alleles (Figure 4). It is important to note that, in our RPPA data (Figure 14), while the presence of total PTEN in its inactive form based on its phosphorylation at S380 [162] is evident for HT72 P0-P3, these levels are relatively low (Figure 14). These variations may lie in the differences between the Western blot and RPPA detection techniques, such as the sample processing procedures, amount of input sample, number of freeze/thaw cycles, or antibodies for RPPA vs. Western blots recognizing different epitopes [163,164,165], rather than inter-slide variability, as the criteria positivity was stringent with a cut-off of CV ≤ 10 [52].

Our study highlights the value of integrating proteomic/phosphoprotein analysis as a multi-OMICS approach in conjunction with assessment of DNA and RNA alterations through PDX passages. While DNA and RNA can provide insight into impactful mutations and/or gene products involved in oncogenic pathway dysregulation, they do not always correlate with protein expression or function due to post-transcriptional and/or post-translational modifications, such as phosphorylation [166]. Our results show that, in general, serial propagation of these pediatric and AYA solid tumor PDXs, a congruence with the original tumors in terms of their histological, molecular, and proteomic characteristics was preserved, particularly in the detailed analyses of functionally impactful DNA sequences.

Our data indicate that eight out ten of our PDX models maintain the biological fidelity similar to that of the original patient tumors. In the two models [1 OS (HT87) and 1 RMS (HT74)] wherein the PDX passages show more discordance when compared to their respective original tumor specimens, it is important to note that many molecular signatures are still maintained, as is evident from CNV, SNV, RPPA, and Western blot data. The discordance at the DNA level observed in some of the P0–PDX cohorts may or may not be of significance in terms of protein function and pathway activation. However, further investigations are needed to determine if genetic alterations in cancer-associated and non-cancer-associated genes have deleterious/damaging impact on their respective protein function and if that contributes to disease pathogenesis in ways not understood at this point in time. There are many caveats to be considered when characterizing a PDX and determining its suitability for efficacy and safety testing from the standpoint of similarities and dissimilarities to the original P0 tumor. As previously mentioned in the introduction, a number of factors can contribute to the differences observed between the pairs of P0 tumor and passaged PDXs, such as clonal evolution, genetic drift, intra- and inter-tumoral heterogeneity, presence of human stroma and immune system in original tumors versus increased murine stromal content in passaged PDXs, and insufficient amount of the original tumor sequenced, as well as lineage splitting due to tumors’ spatial heterogeneity during PDX propagation [16,17,18,19]. Additionally, it should also be noted that it is impossible to know if clonal drift or evolution that appears in some PDXs over time might also naturally occur in that same patient at some point in time over the course of their disease.

Other challenges also exist that still hinder the use of PDXs in precision medicine [167], such as their lack of functional immune system and absence of a human tumor microenvironment [168]. While serial passaging of PDXs results in acceptable preservation of the human patient’s tumor cells, the primary tumor-associated stroma is replaced by murine stroma [57]. The ratio of tumor cells to stroma cells varies between individual tumors regardless of the cancer type [57]; thus, if a particular tumor, such as HT87, has more murine stromal components, this could indicate that the stromal compartments are critical for growth of this tumor. This may also provide insight into possible therapeutic interventions, as it may suggest that targeting the tumor microenvironment may be more beneficial in terms of decreasing growth for this particular tumor. Despite these challenges, PDXs provide a highly valuable preclinical tool for drug screens in pediatric cancers. These models can also be linked to clinically relevant information regarding the original tumor, as we have done here, which provides additional insights into the disease course and clinical outcomes [18] and helps us prioritize targets for therapeutic interventions.

As such, we utilized the integrated multi-OMICS approach to select targets of interest, such as CDK4/6 or BETs, to validate in vivo mechanisms of disease progression in OS PDX models (Figure 16 and Figure 17). Notably, clinical history for patient HT77 indicated that the CDK4/6 inhibitor palbociclib was given to the patient, albeit at a much later phase in the disease’s journey of progression (Figure 16A), such that the therapy may not have been able to compete against the increased tumor burden. However, when palbociclib was given at an earlier stage of the disease, as indicated by the HT77 PDX in vivo study (Figure 16B), tumor growth decreased. Studies are in progress to identify additional therapeutic vulnerabilities and explore combination targeted therapies in these models.

5. Conclusions

In this study, we established an integrated clinical- and research-based platform where molecular data can be obtained on pediatric and AYA patients with pediatric solid tumors to facilitate identification, prioritization, and validation of genomically guided treatments that can benefit patients. Oncogenic heterogeneity, as well as non-standardized methodologies to conduct -OMICS analyses on pediatric and AYA solid tumor PDXs, contribute to the significant lag in improving therapeutic and prognostic outcomes, especially in the relapsed setting [26,169]. Thus, novel clinically annotated, serially passaged PDX models that recapitulate the molecular landscape of the original tumor sample and can be cross-validated will help improve decision making of molecularly guided therapy. Linking patient treatment history with PDX data offers an exciting and underexplored opportunity to not only test new combination therapies based on -OMICS data but to also gain a better understanding of what therapy may be most beneficial to a patient who has already had therapy. Comprehensive identification and prioritization of actionable targets and signaling pathways using integrated multi-OMICS helps identify and cross-validate relevant therapeutic targets in aggressive pediatric and AYA solid tumors that could have been missed from evaluation of only a single layer of the tumor molecular landscape. Efficacy data obtained using these types of well-curated PDX models can provide rationale for moving targeted therapies to future clinical trials earlier in the disease course and at a time before development of chronic adverse events and disease progression, a stage that is nearly impossible to treat [170]. A major benefit of such study is to expand pediatric and AYA PDX panels for future patients and not necessarily for the patients from which the PDX was derived. The field in pediatric and AYA solid tumors is at the point where scientists are generating efficacy and safety data to provide rationale for earlier intervention in future patients. Depending on the clinical situation, tumors are sequenced either upfront or following resection. In the future, we envision the molecular signature and the clinical history of a new patient will be compared and aligned to similar datasets of previous patient(s) from which PDX(s) were also generated. While efficacy data from a PDX would not be the defining component of a clinical decision, it could provide supportive data to help guide decision making for the newly diagnosed or relapsed patient. As such, studies like this continue to be critical to identifying high-yield targets and focusing clinical trial efforts on innovative therapeutics most likely to benefit these underserved pediatric populations. With that said, however, as more PDX models are developed and data is linked to the corresponding patient’s clinical history outcomes, these models may eventually become a more prominent and accepted component of clinical decision making. As such, studies like this continue to be critical to identifying high-yield targets and focusing clinical trial efforts on innovative therapeutics most likely to benefit these underserved pediatric populations.