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.


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.
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 cancerassociated 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.

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.

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.

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 (

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).  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.

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).

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].

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 3 , 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.

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.

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 solidtumor 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. 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.

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].

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 oncogenicassociated genes are differentially impacted by serial PDX passaging.

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 COSMICassociated 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 COSMICassociated 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 COSMICassociated 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].    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 COSMICassociated gene are presented in Table S2A-H. At this level of analysis, it becomes evident  All data were filtered against normal control (NA12878). Copy # = copy numbers; Neutral = 2 copies, Deletion < 2 copies, Gains = 3-5 copies, Amplification > 6 copies. Note: HT72 and HT77 PDX were developed from the same patient at different time points. Due to low quantities of P0 HT77 tissue, P0 HT72 was used for both HT72 and HT77 analyses. All data were filtered against normal control (NA12878). Copy # = copy numbers; Neutral = 2 copies, Deletion < 2 copies, Gains = 3-5 copies, Amplification > 6 copies. Note: HT72 and HT77 PDX were developed from the same patient at different time points. Due to low quantities of P0 HT77 tissue, P0 HT72 was used for both HT72 and HT77 analyses.

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][17][18][19][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].

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).
found to either not occur in cancer-associated genes or to occur in genes w or not fully characterized functions.

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 tworing 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).  The data is normalized to total mappable variant counts specific to that particular sample. These data were normalized against sequencing depth.

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.

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.

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 (Figures 11-13, p-values of ≤2 × 10 −5 ) [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.  . Genes that were statistically significant (p < 0.05) and >1.5-fold change were selected for pathway enrichment analysis using GO, KEGG, and Reactome.
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 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.

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) . 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 (Figures 3 and 4, Table S2A-H) that were then validated for corresponding proteins by Western blot (Figure 4) and/or RPPA (Figure 14). 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) . 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 COS-MIC-associated CNV data indicated copy number gains/amplifications in MYC, RAD21, Cyclin D3, Cyclin E1, as well as loss of CDKN2A and TP53 genes (Figures 3-4, Table S2A- Figure 14. Integration of Western blot and RPPA data for OS. Expression of total and phosphoproteins, as evident by either RPPA or Western blot analysis (refer to figure key), were used to map out a network of biological pathways that may contribute to disease pathogenesis. Some of the selected proteins, such as MYC, PTEN, Cyclin D, and H2AX, were cross-validated by Western blot. For example, COSMIC-associated CNV data indicated copy number gains/amplifications in MYC, RAD21, Cyclin D3, Cyclin E1, as well as loss of CDKN2A and P53 genes that were then validated for corresponding proteins by Western blot and/or RPPA.
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. Figure 15. Activated pathways are preserved among P0 and their respective PDX passages. Heatmap of total and phosphoproteins via RPPA. Sarcomas (HT72, HT77, HT87, and HT74) and Wilms tumor (HT98, HT120, and HT139). Each data point represents the average of 3 replicate printed spots. All final data values were normalized to total protein following the subtraction of negative control. Coefficient of variation (CV) among the replicate spots of <10% [23] was used as the cut-off for positive signals. Arrows indicate targets that are downstream in the pathway.

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,.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 Figure 15. Activated pathways are preserved among P0 and their respective PDX passages. Heatmap of total and phosphoproteins via RPPA. Sarcomas (HT72, HT77, HT87, and HT74) and Wilms tumor (HT98, HT120, and HT139). Each data point represents the average of 3 replicate printed spots. All final data values were normalized to total protein following the subtraction of negative control. Coefficient of variation (CV) among the replicate spots of <10% [23] was used as the cut-off for positive signals. Arrows indicate targets that are downstream in the pathway.
Based on the RRPA data, there is evidence that a number of actionable cancerassociated 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.

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,.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 (Figures 11-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.

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 deidentified manner to validate targeted therapy such that earlier intervention with efficacious and safe molecularly guided therapy can be realized.

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 nonsomatic 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 (Figures 16 and 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.

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.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/cancers15010259/s1. Table S1: Comprehensive clinical information on pediatric and AYA patients donating solid tumors for PDXs, Figure S1: Histological characteristics are preserved in original tumors and their respective PDX passages. Tumor tissues from OS PDXs (HT72, HT77, HT87, and HT96), RMS PDXs (HT74), and Wilms tumor PDXs (HT98, HT120, and HT139) were harvested, formalin fixed, pro-cessed and stained with H&E. The H&E stains are representative of 3 mice per group. All images are at 20× magnification, Figure S2: Percentage of mouse reads filtered out from the established PDX passages. Highly conserved regions between mouse and human genome may result in false positive variants due to presence of mouse stromal tissue. Thus, to overcome this challenge any se-quencing reads that have been mapped to the mouse genome have been removed during WGS of each PDX passage for accurate depiction of the molecular profile. Percentage of mouse reads that were present and removed from each PDX passage are illustrated here. P0 represents the original human patient sample and will have 0% of mouse reads since it has not been introduced in the murine host, Figure S3: HT72 chromplots compare CNV distribution across genome of A P0 tumor and B-C corresponding PDXs. Orange = statistically significant amplification if >0 or deletion if <0; gray = background noise, Figure S4: HT77 chromplots compare CNV distribution across genome of A P0 tumor and B-D corresponding PDXs. HT72 and HT77 are from the same patient, thus P0 tumor of HT72 was used for comparison. Orange = statistically significant amplification if >0 or deletion if <0; gray = background noise, Figure S5: HT87 chromplots compare CNV distribution across genome of A P0 tumor and B-D corresponding PDXs. Orange = statistically significant amplification if >0 or deletion if <0; gray = background noise, Figure S6: HT96 chromplots compare CNV dis-tribution across genome of A P0 tumor and B-D corresponding PDXs. Orange = statistically sig-nificant amplification if >0 or deletion if <0; gray = background noise, Figure S7: HT74 chromplots compare CNV distribution across genome of A parental tumor and B-C corresponding PDXs. Orange = statistically significant amplification if >0 or deletion if <0; gray = background noise, Figure S8: HT98 chromplots compare CNV distribution across genome of A P0 tumor and B-D corresponding PDXs. Orange = statistically significant amplification if >0 or deletion if <0; gray = background noise, Figure S9: HT120 chromplots compare CNV distribution across genome of A P0 tumor and B-D corresponding PDXs Orange = statistically significant amplification if >0 or deletion if <0; gray = background noise, Figure S10: HT139 chromplots compare CNV distribution across genome of A P0 tumor and B-D corresponding PDXs. Orange = statistically significant amplification if >0 or deletion if <0; gray = background noise, Table S2A: Absolute copy numbers identified in COSMIC genes in HT72 cohort, Table S2B: Absolute copy numbers identified in COSMIC genes in HT77 cohort, Table S2C: Absolute copy numbers identified in COSMIC genes in HT87 cohort, Table S2D: Absolute copy numbers identified in COSMIC genes in HT96 cohort, Table S2E: Absolute copy numbers identified in COSMIC genes in HT74 cohort, Table S2F: Absolute copy numbers identified in COSMIC genes in HT98 cohort, Table S2G: Absolute copy numbers identified in COSMIC genes in HT120 cohort, Table S2H: Absolute copy numbers identified in COSMIC genes in HT139 cohort, Table S3A-B: Descriptive analysis of transversions and transitions in each P0-PDX cohort, Table S4: Descriptive analysis of total variants in P0-PDX subset, Figure Table S5A: HT72 oncoplot annotation part 1, Table S5B: HT72 oncoplot annotation part 2, Table S5C: HT77 oncoplot annotation part 1, Table S5D: HT77 oncoplot annotation part 2, Table S5E: HT87 oncoplot annotation part 1, Table S5F: HT87 oncoplot annotation part 2, Table S5G: HT96 oncoplot annotation part 1, Table S5H: HT96 oncoplot annotation part 2, Table S5I: HT74 oncoplot annotation part 1, Table S5J: HT74 oncoplot annotation part 2, Table S5K: HT98 oncoplot annotation part 1, Table S5L: HT98 oncoplot annotation part 2, Table S5M: HT120 oncoplot annotation part 1, Table S5N: HT120 oncoplot annotation part 2, Table S5O: HT139 oncoplot an-notation part 1, Table S5P: HT139 oncoplot annotation part 2, Figure S13: Gene expression profiles are maintained between P0 and its respective PDXs. This is evident by RNA-seq correlation plots for A HT72 (OS), B HT77 (OS), C HT87 (OS), D HT96 (OS), and E HT74 (RMS). The closer the correlation coefficient (R) is to 1 the more positive linear relationship exists between groups being compared [165], Table S6: Integration of transcript levels to copy numbers for prioritized targets, Figure S14: Activated pathways in Wilms tumors are preserved among P0 and their respective PDX passages. Heatmap of total and phosphoproteins via RPPA in Wilms tumor (HT98, HT120, and HT139) is presented. Each data point represents the average of 3 replicate printed spots. All final data values were normalized to total protein following the subtraction of negative control. Coefficient of variation (CV) among the replicate spots of <10% [23] was used as the cut-off for positive signals, Figure S15: Body weights are comparable between A vehicle vs palbociclib in HT77 and B between vehicle vs OTX015 treated groups and for HT96 PDXs. All body weights are represented as mean +/− SEM. Table S7: Intensity values of protein expression obtained from RPPA data following pathway groupings of proteins based on their function, Table S8A: Prioritized therapies for HT72 based on therapeutic response biomarkers and actionable targets, Table S8B: Prioritized therapies for HT77 based on therapeutic response biomarkers and actionable targets, Table S8C: Prioritized therapies for HT87 based on therapeutic response biomarkers and actionable targets, Table S8D: Prioritized therapies for HT96 based on therapeutic response biomarkers and actionable targets, Table S8E: Prioritized therapies for HT74 based on therapeutic response biomarkers and actionable targets, Table S8F: Prioritized therapies for HT120 based on therapeutic response biomarkers and actionable targets. to the Precision Genomics team at the Riley Hospital for Children for their support in obtaining pediatric solid-tumor patient samples as well as the patients and their families for their generous donation of tumor specimens for research. We send a special thanks to Brian Ashmore in the Indiana Pediatric Biobank for tissue processing and archiving. We thank Mary Murray for her discussions and input into clinical history annotations. We thank Bryan Helm for initial discussions and input on -OMICs analyses. We also thank Leila Fehme, Rebekah Addison, and Maggie Granatir for expert technical assistance with histology. WGS for Wilms tumor PDX HT98, HT120 and HT139 was conducted by the Center for Medical Genomics at IUSM. PDX generation and archiving was conducted in by the IUSCCC Preclinical Modeling and Therapeutics Core (PMTC), which are partially supported by the Indiana University Grand Challenges Precision Health Initiative. We also thank Matthew Essex and Art Baluyut for their support and insightful discussions.

Conflicts of Interest:
H.E.S. is a retiree and shareholder of Eli Lilly and Co. J.D.W. is a consultant to and shareholder of Theralink Technologies, Inc. E.F.P. is a consultant to and shareholder of Theralink Technologies, Inc. and Perthera, Inc. All other authors declare that they have no conflict of interest.