Elucidating the transcriptional program of feline injection-site sarcoma using a cross-species mRNA-sequencing approach
BMC Cancer volume 19, Article number: 311 (2019)
Feline injection-site sarcoma (FISS), an aggressive iatrogenic subcutaneous malignancy, is challenging to manage clinically and little is known about the molecular basis of its pathogenesis. Tumor transcriptome profiling has proved valuable for gaining insights into the molecular basis of cancers and for identifying new therapeutic targets. Here, we report the first study of the FISS transcriptome and the first cross-species comparison of the FISS transcriptome with those of anatomically similar soft-tissue sarcomas in dogs and humans.
Using high-throughput short-read paired-end sequencing, we comparatively profiled FISS tumors vs. normal tissue samples as well as cultured FISS-derived cell lines vs. skin-derived fibroblasts. We analyzed the mRNA-seq data to compare cancer/normal gene expression level, identify biological processes and molecular pathways that are associated with the pathogenesis of FISS, and identify multimegabase genomic regions of potential somatic copy number alteration (SCNA) in FISS. We additionally conducted cross-species analyses to compare the transcriptome of FISS to those of soft-tissue sarcomas in dogs and humans, at the level of cancer/normal gene expression ratios.
We found: (1) substantial differential expression biases in feline orthologs of human oncogenes and tumor suppressor genes suggesting conserved functions in FISS; (2) a genomic region with recurrent SCNA in human sarcomas that is syntenic to a feline genomic region of probable SCNA in FISS; and (3) significant overlap of the pattern of transcriptional alterations in FISS with the patterns of transcriptional alterations in soft-tissue sarcomas in humans and in dogs. We demonstrated that a protein, BarH-like homeobox 1 (BARX1), has increased expression in FISS cells at the protein level. We identified 11 drugs and four target proteins as potential new therapies for FISS, and validated that one of them (GSK-1059615) inhibits growth of FISS-derived cells in vitro.
(1) Window-based analysis of mRNA-seq data can uncover SCNAs. (2) The transcriptome of FISS-derived cells is highly consistent with that of FISS tumors. (3) FISS is highly similar to soft-tissue sarcomas in dogs and humans, at the level of gene expression. This work underscores the potential utility of comparative oncology in improving understanding and treatment of FISS.
Feline injection-site sarcoma (FISS) is an aggressive subcutaneous soft-tissue cancer that occurs in approximately 1 in 10,000 domestic cats, frequently at vaccination sites . The etiology of FISS appears to involve non-resolving local inflammation leading to neoplastic transformation of fibroblasts or myofibroblasts and the development of a tumor mass . The triggers of this inflammation probably include vaccine adjuvants and mechanical insult to the subcutis (with latency estimates ranging from months to years ), although predisposing genetic factors for FISS have also been noted [2, 4, 5]. Due to its invasiveness, FISS is difficult to treat; with the current standard of care, radical excision with pre- or post-surgical radiation, recurrence occurs in 15–70% of cases and metastasis in 20% of cats [1, 6, 7]. Average survival may be as short as 6 months  but can be as long as 5 years with surgery and radiotherapy . Results from clinical trials of adjuvant doxorubicin suggest some benefit but are neither curative nor tailored to the patient-specific molecular signature of this histologically diverse neoplasm , and thus, new therapeutic approaches are needed. In the context of human soft tissue sarcoma, new molecular mechanisms and new therapeutic targets have been discovered through the use of comprehensive gene expression profiling of tumors and normal tissues in case-control studies [9,10,11]. Additionally, new insights into pathogenic mechanisms and driver mutations have been uncovered through comparative oncology approaches with transcriptome or exome profiling [12,13,14]. In the context of feline sarcoma, somatic copy number alterations (SCNAs) have been studied genome-wide, revealing recurrent copy number gains and losses , but it is unknown to what extent these alterations correlate with elevated or reduced tumor transcript abundances, though such changes would be expected based on current understanding of cancer biology . However, despite the promise of such systems-biology approaches, the tumor transcriptome of FISS has been neither systematically studied nor compared to the transcriptomes of soft tissue sarcomas in other species such as domestic dogs or humans.
In this study, we used high-throughput short-read mRNA sequencing (mRNA-seq) [15,16,17] to profile the transcriptomes of three FISS tumors and three normal patient-matched tissue samples as well as cultured FISS-derived cell lines and skin fibroblasts. With this approach we identified 3049 transcripts with altered abundance between FISS tumors and normal skin (of which 335 were also differentially expressed in the FISS cell lines vs. cultured fibroblasts), 17 biological processes or molecular pathways whose associated genes are enriched in the set of FISS differentially expressed genes, and nine genomic regions with coherent dysregulation that are probable FISS SCNAs (of which eight are novel per the current state of knowledge ). We then analyzed the degree of consistency between the FISS transcriptome and human cancer gene lists and patterns of differential expression in human and dog sarcomas. Among human orthologs of cat genes with altered expression in FISS vs. normal skin, we found significant enrichment of tumor suppressor genes (TSGs)  and oncogenes . Additionally, we found a 10 Mbp region of potential FISS SCNA whose human syntenic region has a recurrent SCNA in sarcoma. Using ortholog mapping, we compared the pattern of transcriptional dysregulation in FISS with the transcriptional patterns in soft tissue sarcomas in humans and in dogs. We found significant three-way overlap: 53 genes that are upregulated in sarcoma vs. normal tissue in all three species and 38 genes that are downregulated in sarcoma in all three species. Finally, we analyzed the sets of 53 and 38 genes that are consistently upregulated or downregulated using data from a cell line-based drug-to-transcriptome screen, identifying 11 drugs and four drug targets whose reported effects on cell lines suggest possible therapeutic benefit in treating FISS. This work represents the first mRNA-seq study of a feline neoplasm of which we are aware, and represents a template for how cross-species approaches could be used to study other types of feline cancers. The feline mRNA-seq data from this study are available in a public repository as a resource for the scientific community (see Methods).
Through the biobank at the Carlson College of Veterinary Medicine at Oregon State University, we obtained samples from three histologically confirmed FISS tissues and patient-matched normal tissues (two skin samples and one skeletal muscle sample) from the three cats aged 10–13 years. The FISS samples were high grade soft tissue sarcomas from neutered male domestic shorthair cats (thigh or between shoulder blades); none of the cats had metastasis at time of surgery. Tumor volumes ranged up to 500 cm3. The surgical treatment of FISS patients and the diagnostic work-up of tumor samples followed best clinical practices. Two cats were treated by amputation of the affected leg. In the third cat surgical excision included the dorsal processes of thoracic vertebrae. Clinical follow-up data were not available. Through the biobank at the Flint Animal Cancer Center at Colorado State University, we obtained samples from nine histologically confirmed canine soft tissue sarcomas (four grade 3 soft tissue sarcomas, two grade 3 peripheral nerve sheath tumors, two grade 1 myxomatous soft tissue sarcomas, and one grade 2 soft tissue sarcoma) and matched normal skeletal muscle samples. One of the grade 3 soft tissue sarcomas was annotated in the clinical record as possibly being histiocytic sarcoma given its anatomic location (joint capsule) and lymph node metastasis.
Feline FISS cell lines (derived from tumors from cat04, cat05, and cat07) were obtained from the College of Veterinary Medicine at the University of Wisconsin (courtesy of David Vail, Michael Huelsmeyer, and Ilene Kurzman). Primary cultures of feline skin fibroblasts were obtained from the biobank at the Carlson College of Veterinary Medicine at Oregon State University. Cells were grown in high glucose Dulbecco’s Modified Eagle’s Medium (DMEM) containing L-glutamine and sodium pyruvate (Corning #10-013CV, Axygen Incorporated, Union City, CA) with 10% (v/v) fetal bovine serum (FBS) (HyClone #SH3039603, HyClone Laboratories Inc., Logan, UT), 100 U/mL penicillin, and 100 μg/mL streptomycin (“pen/strep”) (Gibco #15140122, Invitrogen, Carlsbad, CA) at 37 °C and 5% CO2. Primary cell cultures were harvested for RNA at third passage and cells were harvested at 70% confluence.
Annotated genome information
For all computational analyses that used the complete sequence of the cat genome, we used the sequence for genome assembly Felis catus 6.2 (Broad Institute, Cambridge, MA, USA; released Sep. 2011) , which we downloaded from the Ensembl database (release 87, Dec. 2016). We obtained information about cat gene locations and exon structures in a Gene Transfer Format file from Ensembl (release 87). We obtained gene annotation information via the BioMart tool from Ensembl (release 87). For genome visualization we used the Integrated Genomics Viewer version 2.3.90 .
We isolated RNA from tissue samples and cell cultures using the Norgen Total RNA Purification Kit #17200 (Norgen Biotek, Thorold, ON), with elution using nuclease-free water.
FISS mRNA-seq profiling
RNA sample library preparation and high-throughput sequencing were performed by the Genomics Core at the Center for Genome Research and Biocomputing at Oregon State University. RNA samples were rRNA-depleted using Ribo-Zero Gold (Illumina, San Diego, CA, USA); strand-specific mRNA-seq libraries were prepared using the PrepX RNA-seq for Illumina Library kit on the Apollo 324 (Wafergen, Fremont, CA, USA); and barcoded libraries were sequenced on a HiSeq 3000 (Illumina) at 2 × 100 bp (paired-end sequencing) on one lane for the first batch of samples (see Additional file 1: Table S1). We generated sequence quality reports using FASTQC  and then aligned the reads to the annotated cat genome using the software tool STAR  (in the alignment, only uniquely aligned reads were retained, and we used basic two-pass mapping, with all first-pass junctions inserted into the genome indices). The alignment yielded an average of 1.0 × 108 mapped reads per sample. Next, we obtained counts of aligned reads per gene with featureCounts (version 1.5.1) using the Subread software program  with the minimum mapping quality score parameter set to the value 3.0 and genome-wide cat gene and exon annotations from Ensembl Release 87 . Given the fibrosarcoma histotypes of the FISS tumors in this study, for the supervised analysis of differential expression in primary tissue, we compared FISS to normal skin tissue only (not muscle). For testing individual genes for differential expression between the sample groups, we used DESeq2  with the Wald test and with p-value adjustment using the Benjamini-Hochberg method . For each gene, we also computed log-scale normalized counts as log2(1 + C), where C is the normalized expression level from DESeq2. We also re-analyzed the mRNA-seq data using the Felis catus 9.0 genome assembly and the Ensembl 95 gene annotations; we compared the gene-level FISS/skin log2 ratios that we obtained using FelCat9 with the gene-level ratios that we obtained using FelCat6.2; they were correlated at R = 0.995.
Principal component analysis (PCA)
For PCA, we used the function prcomp from the R software package stats ver. 3.4.2. For plotting loading scores we used the R package ggbiplot . We selected high-variance genes with minimum threshold of 2.0 for the variance of the normalized log2 expression levels across the samples for the PCA analysis.
Gene set enrichment analysis
For Gene Set Enrichment Analysis (GSEA) , we used the Java JAR file-deployed GSEA software (ver. 2.1.0) in GseaPreranked mode, with genes labeled by official gene symbol and ranked by log2(sarcoma/skin) or log2(sarcoma cells / fibroblast) as appropriate for the sample type (tissue or cultured cells). For the gene sets, we used the C5 gene sets from the MSigDB database (release 5.2) . For the in silico compound screen, we screened gene sets from our study against ten drug-to-cell-line datasets using the Enrichr web tool (ver. 06/28/2017)  with a cutoff of FDR < 0.05. Redundant functional annotations were reduced to representative annotations by manual curation.
Oncogenes and tumor suppressor genes (TSGs) datasets
We obtained lists of human oncogenes from the ONGene database . Additionally, we obtained lists of tumor suppressor genes from the TSGene database . For the oncogene/TSG analysis, we separated all identified feline genes differentially expressed in tissues into orthologs of human TSGs, oncogenes, or genes that are neither TSGs nor oncogenes. For obtaining a “background” gene list, filtered for all genes whose normalized log2 expression count exceeded 2.0 in either the sarcoma or skin sample for cat01.
Reverse Transcription was performed in a 7.5 ng RNA/μL reaction volume using the High Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Foster City, California). For qPCR, 3 μL cDNA was loaded in duplicate with TaqMan primer probes and TaqMan Master Mix (Applied Biosystems, Foster City, California) and reactions were carried out on a StepOnePlus (Applied Biosystems) for 40 cycles. In each sample, the Ct value of each of eight genes (BARX1, BARX2, FBN2, GAPDH, ITM2B, LATS1, LEF1, MMP13) was normalized to the average of the Ct measurements of two endogenous normalizer genes (ATE1 and SHOC2; selected based on their observed stable expression across mRNA-seq samples) from the same sample . The TaqMan assay product numbers (Applied Biosystems) used were as shown in Table 1. For LATS1, a custom TaqMan assay was used based on a 166-mer submitted sequence (provided in Additional file 2: Supplementary Note S1).
Total proteins were isolated from cultured cells in radioimmunoprecipitation assay (RIPA) buffer containing protease and phosphatase inhibitors (EZBlock, BioVision, Milpitas, CA). A 27 gauge needle was used to homogenize cells and cellular debris was removed with two 12,000 g centrifugations for 10 min each at 4 °C. Ten micrograms total protein, diluted in Laemmli sample buffer containing β-mercaptoethanol, were loaded onto 7.5% (for FN1) or 12% (for BARX1) Mini-PROTEAN® TGX™ acrylamide gels (BioRad Laboratories, USA) and separated by electrophoresis in a tris-glycine buffer system. Proteins were transferred to 0.2 μm nitrocellulose (BioRad), stained with Pierce™ Reversible Protein Stain Kit (Thermo Scientific, Rockford, IL), and imaged. The membranes were de-stained, rinsed in TBST, and blocked with 5% bovine serum albumin in tris-buffered saline+Tween 20 (TBST) for 1 hour at room temperature. Anti-Fibronectin (Abcam, ab6328) and Anti-Barx1 (Abcam, ab26156) were diluted at 1:2000 and 1:1000 respectively in blocking buffer, and membranes were incubated overnight at 4 °C. Following TBST rinses, secondary staining with 1:10,000 Goat Anti-Mouse-HRP (Abcam, ab6789 for FN1) and 1:40,000 Goat Anti-rabbit-HRP (Abcam, ab205718, for BARX1) occurred at room temperature for 45 min followed by rinsing in TBST, and detection with ECL Select™ Western Blotting Detection Reagent (Amersham™ #RPN2235, Italy). Image Quant LAS 4000 software was used to acquire the chemiluminescent images and digitization was performed with Image Quant TL.
Genomic block-based test for potential somatic copy number alterations (SCNAs)
We tested for concordant up- or down-regulation (in FISS vs. skin) of mRNA levels for genes within 10 Mbp regions of the genome as a proxy measure for detecting genomic DNA copy number alterations. We aggregated gene-level RNA-seq count data sets into adjacent, non-overlapping 10 Mbp windows, genome-wide. We averaged the gene-level log2(sarcoma/normal tissue) ratios of the genes within each window and then used a genome-wide permutation test (1000 independent shufflings of the genes’ window assignments), rejecting the null hypothesis if the empirical p value (computed by comparing the window-average based on the unshuffled assignments to the sorted vector of window-averages based on the shuffled assignments) satisfied p < 0.001. We further filtered windows to call a “possible SCNA” only if the window-level, absolute average log2(sarcoma/normal tissue) value exceeded 1.0. We mapped between these window regions and previously reported recurrent FISS SCNAs  using chromosome coordinate overlap.
We carried out all statistical tests using the R statistical computing environment (version 3.3.3) and using R software packages from the Bioconductor system  release 3.4 (with the exception of gene functional enrichment analysis; see above sections “Gene Set Enrichment Analysis” and “Enrichr Analysis”). Statistical p-values were adjusted for multiple hypothesis tests using the false discovery rate (FDR) method .
Human soft tissue sarcoma datasets
A mixture of microarray and mRNA-seq transcriptome datasets derived from human soft-tissue sarcoma and normal skeletal muscle were obtained from NCBI GEO  via accession numbers GSE54734, GSE6481, GSE2719, GSE2361, GSE2719, GSE7905, and GSE2193. Gene-level expression data were normalized on each sample (by mean expression level) and then merged. Each gene’s log expression levels across the samples from the merged datasets were tested for a difference of means assuming homoscedasticity (Student’s t-test) or heteroscedasticity (Welch’s t-test) as indicated by an F-test. For each gene for which the equal-means null was rejected at FDR ≤ 0.05), the effect size was computed as average(tumor)/average(normal).
Canine soft tissue sarcoma mRNA-seq profiling
100 ng – 50 μg total RNA from tissue-derived RNA samples were rRNA-depleted using Dynabeads mRNA DIRECT Micro Kit; mRNA-seq libraries were prepared using Ambion RNA-Seq library construction technologies; and barcoded libraries were sequenced on an Ion Proton System at 2 × 100 bp and then aligned to the CanFam3 genome assembly using StrandNGS (Strand Life Sciences Pvt. Ltd., Bengaluru, India) and quantified and normalized using the DEseq algorithm.
Three-species (feline, canine, and human) permutation analysis of soft-tissue sarcoma transcriptomes
For the three-species analysis of the sarcoma transcriptome, we obtained species-specific, primary tissue differential gene expression values (log2(sarcoma/normal)) in three studies: feline (FISS and normal skin), human (soft tissue sarcoma and normal skeletal muscle), and canine (soft tissue sarcoma and normal skeletal muscle). For the 138 genes that were differentially expressed in sarcoma vs. normal tissue in all three species [consisting of 53 genes that are upregulated in sarcoma vs. normal in all three species; 38 genes that are downregulated in sarcoma vs. normal in all three species; and 47 genes that have mixed differential expression], we constructed a 138 × 3 binary matrix (138 genes by three species) based on the differential expression values. The binary matrix contained 1, if the row’s gene is upregulated in sarcoma vs. normal tissue in FISS in the particular column’s species, and 0 otherwise. We then tested the number of rows containing (1,1,1) in each of 100,000 permuted matrices using the R function sample. We computed an empirical p-value as the number of times the count of (1,1,1) rows in a permuted matrix exceeded 53. We performed an identical permutation analysis for downregulated genes, assigning 1 to an entry in the binary matrix if the gene corresponding to the row is downregulated in the species corresponding to the column, and using 100,000 permutations and using 38 as the reference value.
Other genomic datasets
We identified cat-to-human synteny blocks using the VISTA Comparative Genomics web tool  via the Ensembl portal (synteny analysis from Ensembl release 87). For the circos plot, we obtained cytogenetic band data for the cat genome from the UCSC genome browser portal , and used the cytogenetic band data for the human genome that was included with the RCircos software package . We mapped cat genes to human orthologs (and vice-versa), obtained cat gene annotations, and obtained mapped Ensembl gene identifiers to official gene symbols using Ensembl Release 87. We obtained human sarcoma SCNA genes information from the cBioPortal database . For the cross-species analysis of sarcoma SCNAs, we downloaded from cBioPortal a table of cytogenetic band-level amplification/deletion significance test results that were derived from a SCNA study of adult soft tissue sarcomas that was carried out by the Cancer Genome Atlas (TCGA) consortium . The cBioPortal results table was based on analysis of TCGA-published gene-level SCNA data from 265 tumor samples, comprising 58 dedifferentiated liposarcomas, two desmoid/aggressive fibromatosis, 105 leiomyosarcomas, 11 malignant peripheral nerve sheath tumors, 25 myxofibrosarcomas, two pleomorphic liposarcomas, one not otherwise specified, four soft tissue sarcomas, ten synovial sarcomas, and 49 undifferentiated pleomorphic sarcoma/malignant fibrous histiocytoma/high-grade spindle cell sarcomas.
Cells were plated at 10,000/well in a 96-well plate containing DMEM supplemented with 10% FBS and pen/strep and incubated for 24 h at 37 °C and 5% CO2. Media with GSK-1059615 (#11569, Cayman Laboratories, Ann Arbor MI) at 60, 30, 15, 7.5, 3.75, 1.88 μM was added to triplicate wells. Control cells were exposed to the vehicle reagent (dimethyl sulfoxide, DMSO) at concentrations equivalent to the DMSO concentrations in the corresponding GSK-1059615 samples’ media. After 24 h, media was replaced with new GSK-1059615 or control media. At 48 h, media was removed and 0.004% resazurin (Acros Organics #189900010, NJ) in phosphate-buffered saline (PBS) was added. The plate was incubated for 1 h at 37 °C and 5% CO2. Fluorescence was measured at 530 nm excitation, 590 nm emission with a Tecan Infinite F200. Relative Fluorescence Units (RFU) were used to calculate percent of baseline metabolic activity. At each concentration, measurements were acquired in three biological replicates.
All cat mRNA-seq data files for this study have been deposited in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA)  (see Declarations; Availability of Data and Material).
Differential expression in feline injection site-sarcoma (FISS) vs. normal tissue
In order to study transcriptome differences between FISS and normal tissue, we first compared the transcriptomes of FISS tumor tissue, normal skin, and normal skeletal muscle (with muscle included as an outgroup for the normal tissue samples; N = 6 samples total). In an unsupervised analysis (principal components analysis, PCA) of the 6252 genes with significant variance (see Methods), the muscle and skin samples were easily separable from the FISS tumor samples by principal components 1 and 2 (PC1 and PC2) (Fig. 1a, b), which together explained the majority (72%) of the total variance of gene expression. The skin samples and the FISS samples were each tightly grouped in the PCA analysis, and given the fibrosarcoma histotypes of the FISS tumors in this study, for the subsequent supervised analysis of differential expression, we compared FISS to normal skin tissue only, identifying 3049 differentially expressed genes (q < 0.05). In order to identify molecular functions and biological processes that are associated with the sets of genes that were up- and down-regulated in FISS vs. skin, we ranked genes by their log2(sarcoma/skin) values, mapped genes to human orthologs, and then analyzed sets of human genes that are annotated with specific functional ontology terms (Gene Ontology terms)  for enrichment among up- or down-regulated genes using Gene Set Enrichment Analysis (GSEA) . For genes that are upregulated in FISS vs. skin, we found enrichment (q < 0.05 and absolute GSEA enrichment score > 2.0) of gene sets for five categories of gene annotation terms (Fig. 1c): Toll-like receptor signaling, regulation of interleukin 12 (IL-12) production, DNA helicase activity, replication fork, and regulation of interleukin 2 (IL-2) synthesis. For genes that are downregulated in FISS vs. skin, we found enrichment of downregulated genes in four categories of gene sets (Fig. 1d): intermediate filament, long-chain fatty acid metabolism, cell-cell adherens junction, and thioester metabolism.
Differential expression in FISS-derived cell lines vs. skin-derived fibroblasts
As a complement to direct mRNA-seq of FISS tumor and normal tissue samples, we compared the mRNA transcriptomes of cultured FISS tumor-derived cells and skin-derived fibroblasts (two FISS-derived biological replicates and two fibroblast biological replicates each from different cats; N = 8 samples total; see Methods). We identified 3553 genes that were differentially expressed between the FISS-derived cell line and fibroblast sample groups (q < 0.05). By PCA analysis (Fig. 2a, b), the fibroblast samples were more tightly grouped than the FISS-derived samples, with PC1 and PC2 together explaining 73% of the variance in gene expression. Using GSEA, we analyzed genes ranked by their log2(sarcoma cell / fibroblast) expression values in cultured cells, detecting enrichment (q < 0.05, absolute GSEA enrichment score > 2.0) of upregulated genes in three categories of annotation gene sets (Fig. 2c): sister chromatid segregation, DNA replication (consistent with findings from the primary tissue analysis), and double-strand break repair. Additionally, we found enrichment of downregulated genes in five annotation gene sets (Fig. 2d): lysosomal lumen, heparin binding, smooth endoplasmic reticulum, platelet degranulation, and regulation of vascular endothelial growth factor (VEGF) production. For all 21,890 cat genes profiled, the mRNA-seq count data, gene annotation information, and differential expression test results for both the primary tissue analysis and the cultured cells analysis are provided online (Additional file 3: Table S2).
To assess the extent to which the sarcoma cell lines preserve significant molecular characteristics of FISS tumors, we investigated the degree of concordance between the set of genes that are differentially expressed between FISS tumors and skin tissue and the set of genes that are differentially expressed between the FISS-derived cell lines and fibroblasts (Fig. 3). We found 335 genes that were differentially expressed in both the tissue mRNA-seq analysis and the cell culture mRNA-seq analysis, representing a significant overlap (p < 10− 16; Fisher’s exact test; odds ratio 0.60, 95% c.i. 0.53–0.68). Moreover, among the common set of 335 “FISS marker” genes, the degree of concordance in the direction of the differential expression (up in both, or down in both, or up in one analysis and down in the other) was high (Fig. 3a), with an odds ratio of 6.3 (95% c.i. 3.8–10.6), and significantly differs from 1.0 at p < 10− 14 (Fisher’s Exact Test on 2 × 2 contingency table).
Based on magnitude of effect size, consistency of up- or down-regulation in sarcomas in multiple species (see Sec. Gene-level cross-species transcriptome analysis), and/or annotated gene functions suggesting regulatory importance, we selected—from the 241 genes that are consistently (with the same direction of cancer/normal fold-change in both primary tissue and in cultured cells) differentially expressed in FISS vs. normal cells—four FISS-upregulated genes (BARX1, FBN2, MMP13, and LATS1) and four FISS-downregulated genes (BARX2, LEF1, GAPDH, and ITM2B) for targeted validation by quantitative polymerase chain reaction (qPCR). The qPCR measurements were consistent with the mRNA-seq measurements, in terms of the direction and magnitude of the fold-change (sarcoma vs. normal) for seven out of the eight genes tested (Table 1 and Additional file 4: Figure S1; r 2 = 0.879).
Cross-species TSG/oncogene analysis
Insofar as little is known about the genes whose dysregulation in FISS presumably drives disease progression, we tested whether cat orthologs of (i) known human oncogenes, (ii) known human tumor suppressor genes (TSGs), or (iii) human genes that are neither oncogenes nor TSGs have different distributions of log2(sarcoma/normal) expression values. Based on a similar analysis of cancer transcriptome data in our previous study of canine bladder cancer , we hypothesized that cat orthologs of human TSGs would be more likely to be downregulated in FISS, and similarly, that cat orthologs of human oncogenes would be more likely to be upregulated in FISS. Indeed, kernel density analysis of the log2(sarcoma/skin) values for the 1542 cat genes that we had found to be differentially expressed in FISS tumor vs. skin, grouped by the annotation status of the human ortholog of the cat gene (TSG, oncogene, or neitherFootnote 1) were consistent with our hypothesis (Fig. 4a; Table 2). The odds ratios for a gene to be upregulated in FISS (conditioned on the annotation of its human ortholog) differ from 1.0 at p < 0.001 (Fisher’s Exact Test (FET) on 2 × 3 contingency table). We also tested the hypothesis for the cultured cells mRNA-seq data by plotting the kernel density-smoothed distributions of log2(sarcoma cell / fibroblast) values for the 3049 genes that were differentially expressed in sarcoma-derived cells vs. fibroblasts, conditioned on the human ortholog’s annotation (TSG, neither, or oncogene). In the resulting density plots (Fig. 4b) the effect of the human ortholog’s annotation on the density distribution of log ratios was not as strong but it was still evident (Table 2). In this case, the FET was not significant (p = 0.22) but distributions of ranks of log2(sarcoma cell / fibroblast) values in cultured cells differ significantly between the “TSG” and “oncogene” sets of feline genes (p < 0.001, Kolmogorov-Smirnov test). For both the tissue (Fig. 4c) and cultured cells (Fig. 4d) mRNA-seq datasets, the effect of the cat gene’s human ortholog’s annotation status on the proportions of genes that are up- or down-regulated in FISS is clearly evident.
Somatic copy number alteration analysis
Recurrent somatic copy number alterations (SCNAs) have been reported in FISS and for feline soft-tissue sarcomas that are not injection-site-associated  and previous studies of other cancers have successfully detected known SCNAs through region-based analysis of tumor mRNA-seq data . Thus, we analyzed the gene-level differential expression data for FISS vs. skin samples to identify large chromosomal regions of coherent up- or down-regulation at the transcript level. We divided the cat genome into 10 Mbp windows and within each window, we compared the average of the log2(sarcoma/skin) values [denoted by angle brackets, ⟨log2(sarcoma/skin)⟩] for all the genes in the window to a background distribution (see Methods). This permutation analysis yielded a p-value for each window, and we found nine putative SCNA regions for which p < 0.001 and | ⟨log2(sarcoma/skin)⟩| > 1.0. Of these regions, four show coherent upregulation (in sarcoma vs. normal skin) of the genes within them and five show coherent downregulation (Fig. 5a). Of these regions, the region (50–60 Mbp on chromosome Fc-D3) that was the most downregulated transcriptionally in our study was previously reported to have a recurrent somatic deletion in FISS .
Cross-species SCNA analysis
Hypothesizing that SCNAs that drive cancer progression in one species may correspond to SCNAs in the syntenic genomic region for a related cancer in another species, we mapped putative SCNA regions of the cat genome to their syntenic regions in the human genome and then tested whether the human regions have recurrent SCNAs in human soft tissue sarcomas. We found a feline genomic region harboring an SCNA in FISS whose human syntenic region was previously reported to have a recurrent SCNA in soft-tissue sarcoma  (Fig. 5b). The 50–60 Mbp region of Fc-D3 coincides with a previously reported recurrent chromosomal deletion in FISS and in non-injection-site-associated feline sarcoma ; its human syntenic region spans 20.9–80.3 Mbp on Hs-18 (Hs-18q23), which is recurrently deleted in human soft-tissue sarcoma (q < 0.02 based on analysis of SCNA data from ref.  for the 13 genes in Hs-18q23). Detailed information on the nine putative FISS SCNAs is provided online (Additional file 5: Table S3).
Gene-level cross-species transcriptome analysis
Next, we investigated the extent of overlap of the FISS transcriptome with the transcriptomes of human and dog soft-tissue sarcomas. To do so, we obtained differential gene expression (log2(sarcoma/normal)) values for a human study (a published transcriptome profiling study of soft sarcoma and normal skeletal muscle tissues; see Methods) and a canine study (mRNA-seq profiling of banked soft tissue sarcoma and normal skeletal muscle samples; see Methods). Via ortholog mapping, we obtained log2(sarcoma/normal) values for each orthologous gene triple (cat/dog/human) across all three species. We found 138 genes that were differentially expressed in sarcoma vs. normal tissue in all three species: 53 genes that are upregulated in sarcoma vs. normal in all three species (Fig. 6a); 38 genes that are downregulated in sarcoma vs. normal in all three species (Fig. 6c); and 47 genes that have mixed differential expression in the three species (Fig. 6b). Using a permutation approach, we determined that the observed degree of concordance (91 genes out of 138 having consistent fold-change signs in all three species) was extremely unlikely to occur by chance (by permutation test, p < 10− 5 for the upregulated genes and p < 10− 5 for the downregulated genes; see Methods).For the 47 genes with inconsistent directions of differential expression across the three species, the FISS and canine soft-tissue sarcoma had more ortholog pairs with consistent fold-change directions (19 ortholog pairs) than FISS and human soft-tissue sarcoma (13 ortholog pairs) and canine and human soft-tissue sarcoma (15 ortholog pairs).
In silico screen for potential new treatment approaches FISS
Results from previous cancer drug repositioning studies have supported the idea that drugs whose transcriptional effects in human cell lines are substantially opposite (in terms of treatment-vs.-vehicle fold changes, across genes that are transcriptionally altered by the drug treatment) the transcriptional alterations that are measured in tumor vs. normal tissue are more likely to have cancer-inhibiting effects [41, 42]. Thus, we screened the sets of 53 and 38 genes that are consistently dysregulated (upregulated and downregulated, respectively) in soft-tissue sarcoma in cats, dogs, and humans (Fig. 6a, c) against ten previously published drug-to-cell-line response datasets  using the web tool Enrichr  in order to identify drugs and drug targets that may be therapeutically beneficial for sarcoma treatment. For the sarcoma-upregulated genes, the screen identified (FDR < 0.05) five drugs and two drug targets (Fig. 6d; Table 3), and for the sarcoma-downregulated genes, the screen identified six drugs and two drug target protein families (Fig. 6e; Table 3).
Immunoblot analysis of BARX1 and FN1 expression in FISS
Among the 53 genes that are transcriptionally upregulated in FISS as well as in soft-tissue sarcomas (vs. normal tissue) in dogs and humans, we validated at the protein level the upregulation of two genes (BARX1, BarH-like homeobox 1; and FN1, fibronectin 1) in FISS-derived cells vs. skin-derived fibroblasts. By immunoblot with quantitative image analysis we measured that BARX1 protein is three-fold more abundant in FISS-derived cells than in skin-derived fibroblasts (Fig. 6f, Additional file 6: Table S4, Additional file 7: Figure S2), and that FN1 protein is 2.1-fold more abundant in a FISS-derived cell line (cat05) than in skin-derived fibroblasts (Fig. 6g, Additional file 8: Table S5, Additional file 7: Figure S2), though FN1 was not consistently upregulated at the protein level in all FISS-derived cell lines (Additional file 7: Figure S2).
IC50 assay of GSK-1059615
>Because the compound GSK-1059615 (an inhibitor of phosphoinositide 3-kinases (PI3K) and mammalian target of rapamycin (mTOR)) had the strongest enrichment in the in silico screen (Fig. 6d) and based on its molecular targets as described in Discussion, we selected it for an in vitro growth inhibition assay. We found that 48 h incubation with GSK-1059615 potently inhibited growth of FISS-derived cells in vitro, with an IC50 of 4.6 μM (Fig. 7).
This work was principally motivated by three questions: (1) What genes are differentially expressed in FISS, and what are their biological pathways and functions? (2) In which genomic regions are genes coherently differentially expressed in FISS vs. normal tissue? and (3) To what extent does the FISS transcriptome overlap with those of dog and human sarcomas?
The 3049 genes that we detected with altered transcript abundances in FISS tumors vs. normal skin comprise genes in multiple pathways that could plausibly play a role in the pathogenesis or progression of FISS.
Sarcoma-upregulated genes that are common to sarcomas in the three species
Among the genes that are upregulated in sarcoma in all three species, several genes stand out: FAP (fibroblast activation protein α) is known to be expressed in soft tissue sarcomas and is involved in control of fibroblast growth , WT1 (Wilms tumor suppressor) is expressed in over half of soft tissue sarcomas and has prognostic significance , and PRAME (preferentially expressed antigen in melanoma) promotes cell growth, is expressed in tumors including soft tissue sarcomas, and has been proposed as a potential immunotherapy target for sarcoma . Additionally, several enriched gene annotations are notable. Many of the upregulated genes are associated with the immune system or complement system (TLR2, CD14, C3AR1, ITGA4, CXCR4, IFI44, CSF3R, CSF1R, CD53, CD48, IL6, CCL13, CD86, LY86, LCP1, OAS, SRGN, CTSS, TREM1) or the extracellular matrix (ECM) (TNC, SPP1, MMP13, LRRC15, LAMA4, FN1, FBN2, FAP, COL6A3, COCH, ADAM28). Other cellular functions or biological processes with which multiple consistently upregulated genes are annotated include cell proliferation (PTH1R, PRAME, MCM5, IGFBP4, GPC3), cytoskeleton (MAP1B, SEPT11, KRT7, ASAP1), gene regulation (WT1, BICC1, BARX1), lysosome (LAPTM5, CTSS, CTSK), signaling (GPC3, RGS1, RAB31), trafficking (TGONL2, SNX, COPB1), and metabolism (FFAR2, ALDH18A1).
The finding of enrichment of ECM-related genes among the 53 genes that are upregulated in sarcoma is consistent with the abundant ECM typically seen in histopathological analysis of FISS tumors  and the upregulation of both ECM structural constituents (COL6A3, FBN2, LAMA4), as well as degradative enzymes (MMP13, ADAM28) is suggestive of ECM remodeling. Of note, increased expression of orthologous genes regulating ECM remodeling was also seen in human and canine sarcomas (Fig. 6a), which may suggest that some molecular mechanisms of FISS invasiveness are shared with soft-tissue sarcomas in dogs and humans.
At the level of specific cytokine pathways, we found enrichments for genes associated with two specific cytokines (interleukin-2 and interleukin-12) among the 53 genes that are upregulated in sarcoma. Genes annotated as having functions related to regulation of production of interleukin-2 (IL-2) (e.g., CD86, CD28, PRKCQ, CD80, IL1A, LAG3, IL1B, IL17F, CD276, GLMN) were enriched at FDR < 0.05 with a GSEA enrichment score of 2.0. Interleukin-2 is primarily produced by activated CD4+ T lymphocytes, and it promotes the differentiation of immature T cells into regulatory T cells. Lymphocytic aggregates have been reported in approximately 59% of FISS cases, and these are typically enriched for T cells . IL-2 has antitumor activity and, in humans, recombinant IL-2 immunotherapy has been used clinically for treating melanoma. In the context of feline fibrosarcoma, adjuvant immunotherapy with local injection of recombinant human IL-2 extended survival times  and local canarypox-vectored recombinant IL-2 immunotherapy significantly reduced recurrence rates in cat fibrosarcomas [64, 65]. Genes annotated as having functions related to the regulation of production of the cytokine interleukin-12 (IL-12) (e.g., ACP5, CCR7, MAST2, IDO1, RIPK2, TLR4, IL10, MEFV, TLR2, IRF8, NFKB1, IFNG, JAK3, TIRAP, IL12B) were enriched among genes upregulated in FISS at FDR < 0.02 with a GSEA enrichment score of 2.2. Interleukin-12 is produced by activated phagocytic cells such as macrophages, dendritic cells, and neutrophils; it influences both innate and adaptive arms of the immune system, including augmenting cytotoxic CD8+ T-cell activation , reprogramming T helper 17 (TH17) cells into a TH1 phenotype, and reprogramming tumor-associated antigen-presenting cells to prime cytotoxic T-cells for antitumor activity [66, 67]. Inflammation is a hallmark of FISS, with lymphoplasmacytic infiltrates always present in the tumor as well as macrophages and neutrophils; furthermore, significant aggregates of macrophages are present in or at the periphery of approximately 25% of tumors . In a mouse xenograft model of human rhabdomyosarcoma, administration of anti-histone antibody-conjugated IL-12 caused long-term remission and increased survival . In cats, systemic IL-12 administration for treatment of spontaneous soft-tissue sarcoma has been evaluated in a phase-one clinical trial and achieved dose-dependent delivery of IL-12 to the tumor . Overall and together with our transcriptional findings for genes related to IL-12 regulation, these results are suggestive that IL-12 immunotherapy may prove beneficial in the treatment of FISS.
Sarcoma-downregulated genes that are common to sarcomas in the three species
Among the genes that are downregulated in sarcoma in all three species, there are three known tumor suppressors: TP63 , EPHA1 , and DUSP26 (which may alternately function as an oncogene depending on the cancer context) . Many of the downregulated genes are associated with membrane transport (UNC93A, SLC9A2, SLCL5A1, SLC39A2, SLC38A4, PRSS8, CLDN8, CLDN4, AQP3, ABCD2, SMPD3) or metabolism (MAOB, HSD17B2, HMGCS2, FA2H, DCT, ALOX12B). Other cellular functions or biological processes for which multiple consistently downregulated genes are annotated include gene regulation (TFCP2L1, RORC, CPEB3, BARX2), signaling (PTPN3, LY6G6C), immune function (MBP, AIM1L), cytoskeleton (SNTB1, PACSIN3, MAPT, KRT16, GDPD2, COBL; see also Fig. 1d, “intermediate filament”), and cell adhesion (ITGB6, CDH15, CARD14; see also Fig. 1d, “cell-cell adherens junction”). The latter gene category is perhaps unsurprising given that there is altered and often reduced cell-cell adhesion and cell-extracellular matrix (cell-ECM) interaction in many cancers  and given the invasive nature of FISS. Among the cross-species downregulated genes, there are two additional genes that particularly noteworthy: SERPINB13 is an angiogenesis inhibitor  that is downregulated in head and neck cancers , and the neutral sphingomyelinase SMPD3 is a cell cycle regulator .
The finding that many genes with functions in fatty acid (FA) metabolism, such as multiple FA elongases (Additional file 3: Table S2), fatty acid 2-hydroxylase (FA2H), and arachidonate 12-lipoxygenase (ALOX12B, which produces an eicosanoid signaling molecule that stimulates epidermal lipid envelope synthesis; see Fig. 6c) are downregulated in sarcoma vs. normal skin (see also Fig. 1d, “long-chain fatty acid metabolism”) is interesting. Many tumors are lipogenic in order to sustain cell proliferation and tumor growth  and skin has a modest overall rate of FA synthesis, on a per-gram basis, compared to other normal tissues . On the other hand, the finding that the angiogenesis-inhibiting gene SERPINB13 is expressed at lower levels in soft-tissue sarcoma vs. normal tissue is not surprising (Fig. 6c), given the strong evidence for the role of the neovasculature in supporting sarcoma tumor progression .
Parallel examination of tissue samples and cultured cells reveals that there is an overall high degree of concordance in differential expression of genes in FISS tissue versus normal skin compared to FISS-derived cell lines versus cultured fibroblasts. Overall, this finding strengthens the utility of this in vitro model with one caveat. Our cross-species analysis of tumor suppressor gene (TSG)/oncogene expression suggests that bias in differential expression of feline orthologs to known oncogenes and TSG is not as strong in cultured cells as it is in primary tissue (Table 2); this could be a reflection of shared selective pressure imposed by environmental parameters of the cell culture systems.
Over 158 recurring focal SCNAs have been identified in a recent survey of copy number alterations in human cancers and cancer cell lines . Although not all of these SCNAs contain known cancer-target genes, most occur in more than one human cancer type. Human soft tissue sarcoma is known to have many recurrent SCNAs, including deletions and amplifications, that range from focal to broad in scale . Our region-based analysis of tumor mRNA-seq data for coherently up- or down-regulated genes in FISS versus skin identified nine putative recurring SCNAs. By our mRNA-seq-based method, we found one putative SCNA (a probable deletion on Fc-D3) that overlapped with a recurrent SCNA detected in a previous DNA-based analysis of feline sarcoma . Overall, this modest amount of overlap is not unexpected given the limited size of the primary tissue dataset used in this study and the differing measurement modality. However, our finding does bolster the likely functional significance of the recurrent Fc-D3 deletion in FISS. One of the nine putative SCNA regions, Fc-C1 (70–80 Mbp) which had the highest-magnitude overall effect on transcription of any of the putative SCNA regions in FISS, is particularly noteworthy. Its human syntenic region, Hs-18q23, is recurrently deleted in sarcomas and contains two genes, GALR1 and CYB5A, that are annotated as tumor suppressors .
The eight genes whose differential expression in FISS vs. normal skin we chose to validate by qPCR include three transcription factors (BARX1 and BARX2 , which encode BarH-like homeobox transcription factors 1 and 2 which are normally expressed in gastrointestinal tissue; and LEF1, which encodes lymphoid enhancer binding factor 1, a transcription factor downstream of the Wnt / β-catenin pathway ); two genes encoding extracellular proteins (MMP13, which encodes matrix metallopeptidase 13, which is involved in degrading extracellular matrix proteins such as fibrillar collagen and fibronectin and which has a role in controlling angiogenesis in some cancers ; and FBN2, which encodes an extracellular protein, fibrillin-2, that is a structural component of microfibrils in the extracellular matrix ); and a tumor suppressor (LATS1, which encodes large tumor suppressor kinase 1, a serine/threonine kinase that complexes with cell division cycle protein 1 (CDC1) in early mitosis ). Overall, we validated seven out of the eight genes originally identified as differentially expressed through mRNA-seq that were assayed by qPCR, with the exception being LATS1 (notably, it was the only gene for which a custom hydrolysis probe-based qPCR assay had to be designed due to the absence of an off-the-shelf commercial assay). Among all eight genes, the expression ratios as measured by mRNA-seq and qPCR were highly concordant, with r ≠ 0 significant at P < 0.0001. Given that BARX1 transcript abundance is nearly forty-fold higher in FISS than in skin, and given our interest in the FISS-upregulated gene FN1 (a five-fold FISS-upregulated gene that encodes fibronectin-1, an ECM protein that has an important role in cell adhesion, growth, migration, and differentiation and whose expression is associated with more advanced disease in renal cancer ), we selected BARX1 and FN1 to assay at the protein level. The findings that the gastrointestinal transcription factor BARX1 is upregulated (three-fold) in FISS vs. normal skin and that the ECM protein FN1 is upregulated (two-fold) in some FISS-derived cell lines vs. normal fibroblasts are worthy of further focused investigation to determine their specific functions in the pathogenesis of FISS.
Of the fifteen drugs and drug targets that were identified through the screen of drug-to-cell-line databases using the sets of genes that are coherently upregulated or downregulated in sarcoma (across studies in three species), the compound GSK-1059615 is particularly interesting. GSK-1059615 inhibits PI3K and mTOR and it has been investigated as a potential antineoplastic in humans . The finding (Fig. 7) that GSK-1059615 potently inhibits growth of FISS-derived cells suggests that dual PI3K and mTOR targeting is worthy of further study as a potential therapeutic approach for FISS. Additionally, the platelet-derived growth factor BB monomer (PDGF-BB) is noteworthy as a potential target. Members of the platelet derived growth factor family, including PDGF-BB, act as mitogens on mesenchymal cells, often in an autocrine fashion . Depending on the cell type, PDGF-BB signaling can activate cell proliferation in the absence of inflammation. Tumor-derived PDGF-BB promotes pericyte recruitment and leads to increased stability of intratumoral vasculature, tumor cell proliferation, and survival. Notably, the anti-PDGF-A antibody olaratumab has shown promise in treatment of advanced human soft tissue sarcoma . Together, these previous and new findings suggest that PDGF-BB inhibition may be worthy of study in the treatment of FISS. Also intriguing is the identification of the drugs AS601245, a Janus N-terminal kinase inhibitor, and rosiglitazone, a PPARɣ selective agonist, as their co-administration to colon cancer cells has been reported to synergistically reduce cell migration .
In addition to providing a comprehensive molecular picture of FISS that will inform future targeted studies, this work has three principal conclusions. First, our results provide strong evidence that window-based analysis of mRNA-seq data can uncover somatic copy number alterations. Second, the results show that the transcriptome of FISS-derived cells is highly consistent with that of FISS tumors, bolstering the relevance of FISS-derived cultured cells as a model of FISS tumor cancer cells. Finally, the results of this study show that at the level of gene expression, FISS has a high degree of similarity to soft-tissue sarcomas in dogs and humans. While the use of systems biology approaches to map the molecular basis of feline cancers is still in its early stages, the present work illustrates the promise of integrating mRNA-seq with cross-species analytical approaches for uncovering candidate mechanisms underlying FISS pathogenesis and for uncovering new therapeutic approaches (such as dual targeting of the PI3K and mTOR pathways) for this challenging disease.
This included cat genes without human orthologs.
feline injection site sarcoma
gene set enrichment analysis
principal component analysis
somatic copy number alteration
Zabielska-Koczywąs K, Wojtalewicz A, Lechowski R. Current knowledge on feline injection-site sarcoma treatment. Acta Vet Scand. 2017;59:47.
Porcellato I, Menchetti L, Brachelente C, Sforna M, Reginato A, Lepri E, et al. Feline injection-site sarcoma. Vet Pathol. 2017;54:204–11.
Martano M, Morello E, Buracco P. Feline injection-site sarcoma: past, present and future perspectives. Vet J. 2011;188:136–41.
Thomas R, Valli VE, Ellis P, Bell J, Karlsson EK, Cullen J, et al. Microarray-based cytogenetic profiling reveals recurrent and subtype-associated genomic copy number aberrations in feline sarcomas. Chromosom Res. 2009;17:987–1000.
Kirpensteijn J. Feline injection site-associated sarcoma: is it a reason to critically evaluate our vaccination policies? Vet Microbiol. 2006;117:59–65.
Hendrick MJ, Brooks JJ. Postvaccinal sarcomas in the cat: histology and immunohistochemistry. Vet Pathol. 1994;31:126–9.
Phelps HA, Kuntz CA, Milner RJ, Powers BE, Bacon NJ. Radical excision with five-centimeter margins for treatment of feline injection-site sarcomas: 91 cases (1998-2002). J Am Vet Med Assoc. 2011;239:97–106.
Rossi F, Marconato L, Sabattini S, Cancedda S, Laganga P, Leone VF, et al. Comparison of definitive-intent finely fractionated and palliative-intent coarsely fractionated radiotherapy as adjuvant treatment of feline microscopic injection-site sarcoma. J Feline Med Surg. 2018:1098612X18758883.
Network CGAR. Comprehensive and integrated genomic characterization of adult soft tissue sarcomas. Cell. 2017;171:950–65.e28.
Takahashi A, Nakayama R, Ishibashi N, Doi A, Ichinohe R, Ikuyo Y, et al. Analysis of gene expression profiles of soft tissue sarcoma using a combination of knowledge-based filtering with integration of multiple statistics. PLoS One. 2014;9:e106801.
Baird K, Davis S, Antonescu CR, Harper UL, Walker RL, Chen Y, et al. Gene expression profiling of human sarcomas: insights into sarcoma biology. Cancer Res. 2005;65:9226–35.
Ramsey SA, Xu T, Goodall C, Rhodes AC, Kashyap A, He J, et al. Cross-species analysis of the canine and human bladder cancer transcriptome and exome. Genes Chromosomes Cancer. 2017;56:328–43.
Fowles JS, Brown KC, Hess AM, Duval DL, Gustafson DL. Intra- and interspecies gene expression models for predicting drug response in canine osteosarcoma. BMC Bioinformatics. 2016;17:93.
Dhawan D, Paoloni M, Shukradas S, Choudhury DR, Craig BA, Ramos-Vara JA, et al. Comparative gene expression analyses identify luminal and basal subtypes of canine invasive urothelial carcinoma that mimic patterns in human invasive bladder Cancer. PLoS One. 2015;10:e0136688.
Xiong Y, Wu S, Du Q, Wang A, Wang Z. Integrated analysis of gene expression and genomic aberration data in osteosarcoma (OS). Cancer Gene Ther. 2015;22:524–9.
Lister R, O’Malley RC, Tonti-Filippini J, Gregory BD, Berry CC, Millar AH, et al. Highly integrated single-base resolution maps of the epigenome in Arabidopsis. Cell. 2008;133:523–36.
Bainbridge MN, Warren RL, Hirst M, Romanuik T, Zeng T, Go A, et al. Analysis of the prostate cancer cell line LNCaP transcriptome using a sequencing-by-synthesis approach. BMC Genomics. 2006;7:246.
Zhao M, Sun J, Zhao Z. TSGene: a web resource for tumor suppressor genes. Nucleic Acids Res. 2013;41 Database issue:D970–6.
Liu Y, Sun J, Zhao M. ONGene: a literature-based database for human oncogenes. J Genet Genomics. 2017;44:119–21.
Pontius JU, Mullikin JC, Smith DR, Sequencing Team A, Lindblad-Toh K, Gnerre S, et al. Initial sequence and comparative analysis of the cat genome. Genome Res. 2007;17:1675–89.
Thorvaldsdottir H, Robinson JT, Mesirov JP. Integrative genomics viewer (IGV): high-performance genomics data visualization and exploration. Brief Bioinform. 2012;14:178–92.
Andrews S. FASTQC: a quality control tool for high throughput sequence data. 2015. https://www.bioinformatics.babraham.ac.uk/projects/fastqc/. Accessed 2018.
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21.
Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30:923–30.
Cunningham F, Amode MR, Barrell D, Beal K, Billis K, Brent S, et al. Ensembl 2015. Nucleic Acids Res. 2015;43(Database issue):D662–9.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing on JSTOR. J R Stat Soc Ser B Methodol. 1995;57:289–300.
Vu V. ggbiplot: a biplot based on ggplot2. 2015. https://github.com/vqv/ggbiplot. Accessed 2015.
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102:15545–50.
Kuleshov MV, Jones MR, Rouillard AD, Fernandez NF, Duan Q, Wang Z, et al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. 2016;44:W90–7.
Litvak KJ, Shmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods. 2001;25(4):402–8.
Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, et al. 10.1186/gb-2004-5-10-r80. Genome Biol 2004;5:R80. doi:https://doi.org/10.1186/gb-2004-5-10-r80.
Barrett T, Wilhite SE, Ledoux P, Evangelista C, Kim IF, Tomashevsky M, et al. NCBI GEO: archive for functional genomics data sets--update. Nucleic Acids Res. 2013;41(Database issue):D991–5.
Frazer KA, Pachter L, Poliakov A, Rubin EM, Dubchak I. VISTA: computational tools for comparative genomics. Nucleic Acids Res. 2004;32 Web Server issue:W273–9.
Kuhn RM, Haussler D, Kent WJ. The UCSC genome browser and associated tools. Brief Bioinform. 2013;14:144–61.
Zhang H, Meltzer P, Davis S. RCircos: an R package for Circos 2D track plots. BMC Bioinformatics. 2013;14:244.
Gao J, Aksoy BA, Dogrusoz U, Dresdner G, Gross B, Sumer SO, et al. Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci Signal. 2013;6:l1.
Leinonen R, Sugawara H, Shumway M, Collaboration INSD. The sequence read archive. Nucleic Acids Res. 2010;39:D19–21.
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. The gene ontology consortium. Nat Genet. 2000;25:25–9.
Patel AP, Tirosh I, Trombetta JJ, Shalek AK, Gillespie SM, Wakimoto H, et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science. 2014;344:1396–401.
Chen H-R, Sherr DH, Hu Z, DeLisi C. A network based approach to drug repositioning identifies plausible candidates for breast cancer and prostate cancer. BMC Med Genet. 2016;9:51.
Shigemizu D, Hu Z, Hung J-H, Huang C-L, Wang Y, DeLisi C. Using functional signatures to identify repositioned drugs for breast, myelogenous leukemia and prostate Cancer. PLoS Comput Biol. 2012;8:e1002347.
Lamb J, Crawford ED, Peck D, Modell JW, Blat IC, Wrobel MJ, Lerner J, Brunet JP, Subramanian A, Ross KN, Reich M, Hieronymus H, Wei G, Armstrong SA, Haggarty SJ, Clemons PA, Wei R, Carr SA, Lander ES, Golub TR. The Connectivity Map: using gene-expression signatures to connect small molecules, genes, and disease. Science. 2006;313(5795):1929–35.
Raj L, Ide T, Gurkar AU, Foley M, Schenone M, Li X, et al. Selective killing of cancer cells by a small molecule targeting the stress response to ROS. Nature. 2011;475:231–4.
Lim D, Kim KS, Kim H, Ko K-C, Song JJ, Choi JH, et al. Anti-tumor activity of an immunotoxin (TGFα-PE38) delivered by attenuated Salmonella typhimurium. Oncotarget. 2017;8:37550–60.
Finch PW, Rubin JS. Keratinocyte growth factor/fibroblast growth factor 7, a homeostatic factor with therapeutic potential for epithelial protection and repair. Adv Cancer Res. 2004;91:69–136.
Hennequin LF, Allen J, Breed J, Curwen J, Fennell M, Green TP, et al. N-(5-chloro-1,3-benzodioxol-4-yl)-7-[2-(4-methylpiperazin-1-yl)ethoxy]-5- (tetrahydro-2H-pyran-4-yloxy)quinazolin-4-amine, a novel, highly selective, orally available, dual-specific c-Src/Abl kinase inhibitor. J Med Chem. 2006;49:6465–88.
Hirai H, Sootome H, Nakatsuru Y, Miyama K, Taguchi S, Tsujioka K, et al. MK-2206, an allosteric Akt inhibitor, enhances antitumor efficacy by standard chemotherapeutic agents or molecular targeted drugs in vitro and in vivo. Mol Cancer Ther. 2010;9:1956–67.
Mulvihill MJ, Cooke A, Rosenfeld-Franklin M, Buck E, Foreman K, Landfair D, et al. Discovery of OSI-906: a selective and orally efficacious dual inhibitor of the IGF-1 receptor and insulin receptor. Future Med Chem. 2009;1:1153–71.
Xie J, Li Q, Ding X, Gao Y. GSK1059615 kills head and neck squamous cell carcinoma cells possibly via activating mitochondrial programmed necrosis pathway. Oncotarget. 2017;8:50814–23.
Heldin C-H. Targeting the PDGF signaling pathway in tumor treatment. Cell Commun Signal. 2013;11:97.
Pang X, Shu Y, Niu Z, Zheng W, Wu H, Lu Y, et al. PPARγ1 phosphorylation enhances proliferation and drug resistance in human fibrosarcoma cells. Exp Cell Res. 2014;322:30–8.
Imura Y, Nakai T, Yamada S, Outani H, Takenaka S, Hamada K, et al. Functional and therapeutic relevance of hepatocyte growth factor/c-MET signaling in synovial sarcoma. Cancer Sci. 2016;107:1867–76.
Cerbone A, Toaldo C, Minelli R, Ciamporcero E, Pizzimenti S, Pettazzoni P, et al. Rosiglitazone and AS601245 decrease cell adhesion and migration through modulation of specific gene expression in human colon cancer cells. PLoS One. 2012;7:e40149.
Barrett SD, Bridges AJ, Dudley DT, Saltiel AR, Fergus JH, Flamme CM, et al. The discovery of the benzhydroxamate MEK inhibitors CI-1040 and PD 0325901. Bioorg Med Chem Lett. 2008;18:6501–4.
Kruidenier L, Chung C-W, Cheng Z, Liddle J, Che K, Joberty G, et al. A selective jumonji H3K27 demethylase inhibitor modulates the proinflammatory macrophage response. Nature. 2012;488:404–8.
Spitzenberg V, König C, Ulm S, Marone R, Röpke L, Müller JP, et al. Targeting PI3K in neuroblastoma. J Cancer Res Clin Oncol. 2010;136:1881–90.
Davies BR, Logie A, McKay JS, Martin P, Steele S, Jenkins R, et al. AZD6244 (ARRY-142886), a potent inhibitor of mitogen-activated protein kinase/extracellular signal-regulated kinase kinase 1/2 kinases: mechanism of action in vivo, pharmacokinetic/pharmacodynamic relationship, and potential for combination in preclinical models. Mol Cancer Ther. 2007;6:2209–19.
Dohi O, Ohtani H, Hatori M, Sato E, Hosaka M, Nagura H, et al. Histogenesis-specific expression of fibroblast activation protein and dipeptidylpeptidase-IV in human bone and soft tissue tumours. Histopathology. 2009;55:432–40.
Kim A, Park EY, Kim K, Lee JH, Shin DH, Kim JY, et al. Prognostic significance of WT1 expression in soft tissue sarcoma. World J Surg Oncol. 2014;12:214.
Roszik J, Wang W-L, Livingston JA, Roland CL, Ravi V, Yee C, et al. Overexpressed PRAME is a potential immunotherapy target in sarcoma subtypes. Clin Sarcoma Res. 2017;7:11.
Couto SS, Griffey SM, Duarte PC, Madewell BR. Feline vaccine-associated fibrosarcoma: morphologic distinctions. Vet Pathol. 2002;39:33–41.
Quintin-Colonna F, Devauchelle P, Fradelizi D, Mourot B, Faure T, Kourilsky P, et al. Gene therapy of spontaneous canine melanoma and feline fibrosarcoma by intratumoral administration of histoincompatible cells expressing human interleukin-2. Gene Ther. 1996;3:1104–12.
Jasa D, Soyerc C, Fornel-Thibaudb PD, Oberlia F, Vernesa D, Guigala P-M, et al. Adjuvant immunotherapy of feline injection-site sarcomas with the recombinant canarypox virus expressing feline interleukine-2 evaluated in a controlled monocentric clinical trial when used in association with surgery and brachytherapy. Trials in Vaccinology. 2015;4:1–8.
Jourdier T-M, Moste C, Bonnet M-C, Delisle F, Tafani J-P, Devauchelle P, et al. Local immunotherapy of spontaneous feline fibrosarcomas using recombinant poxviruses expressing interleukin 2 (IL2). Gene Ther. 2003;10:2126–32.
Henry CJ, Ornelles DA, Mitchell LM, Brzoza-Lewis KL, Hiltbold EM. IL-12 produced by dendritic cells augments CD8+ T cell activation through the production of the chemokines CCL1 and CCL17. J Immunol. 2008;181:8576–84.
Kerkar SP, Goldszmid RS, Muranski P, Chinnasamy D, Yu Z, Reger RN, et al. IL-12 triggers a programmatic change in dysfunctional myeloid-derived cells within mouse tumors. J Clin Invest. 2011;121:4746–57.
Schilbach K, Alkhaled M, Welker C, Eckert F, Blank G, Ziegler H, et al. Cancer-targeted IL-12 controls human rhabdomyosarcoma by senescence induction and myogenic differentiation. Oncoimmunology. 2015;4:e1014760.
Siddiqui F, Li C-Y, Larue SM, Poulson JM, Avery PR, Pruitt AF, et al. A phase I trial of hyperthermia-induced interleukin-12 gene therapy in spontaneously arising feline soft tissue sarcomas. Mol Cancer Ther. 2007;6:380–9.
Guo X, Keyes WM, Papazoglu C, Zuber J, Li W, Lowe SW, et al. TAp63 induces senescence and suppresses tumorigenesis in vivo. Nat Cell Biol. 2009;11:1451–7.
Mäki-Nevala S, Sarhadi VK, Knuuttila A, Scheinin I, Ellonen P, Lagström S, et al. Driver gene and Novel mutations in Asbestos-exposed lung adenocarcinoma and malignant mesothelioma detected by exome Sequencing. Lung. 2016;194:125–35.
Patterson KI, Brummer T, Daly RJ, O’Brien PM. DUSP26 negatively affects the proliferation of epithelial cells, an effect not mediated by dephosphorylation of MAPKs. Biochimica et Biophysica Acta (BBA) - Molecular Cell Research. 2010;1803:1003–12.
Cavallaro U, Christofori G. Cell adhesion in tumor invasion and metastasis: loss of the glue is not enough. Biochim Biophys Acta. 2001;1552:39–45.
Shellenberger TD, Mazumdar A, Henderson Y, Briggs K, Wang M, Chattopadhyay C, et al. Headpin: a serpin with endogenous and exogenous suppression of angiogenesis. Cancer Res. 2005;65:11501–9.
de Koning PJA, Bovenschen N, Leusink FKJ, Broekhuizen R, Quadir R, van Gemert JTM, et al. Downregulation of SERPINB13 expression in head and neck squamous cell carcinomas associates with poor clinical outcome. Int J Cancer. 2009;125:1542–50.
Stoffel W, Hammels I, Jenke B, Binczek E, Schmidt-Soltau I, Brodesser S, et al. Neutral sphingomyelinase (SMPD3) deficiency disrupts the Golgi secretory pathway and causes growth inhibition. Cell Death Dis. 2016;7:e2488.
Menendez JA, Lupu R. Fatty acid synthase and the lipogenic phenotype in cancer pathogenesis. Nat Rev Cancer. 2007;7:763–77.
Gandemer G, Pascal G, Durand G. Lipogenic capacity and relative contribution of the different tissues and organs to lipid synthesis in male rat. Reprod Nutr Dev. 1983;23:575–86.
Rocchi L, Caraffi S, Perris R, Mangieri D. The angiogenic asset of soft tissue sarcomas: a new tool to discover new therapeutic targets. Biosci Rep. 2014;34:e00147.
Beroukhim R, Mermel CH, Porter D, Wei G, Raychaudhuri S, Donovan J, et al. The landscape of somatic copy-number alteration across human cancers. Nature. 2010;463:899–905.
Barretina J, Taylor BS, Banerji S, Ramos AH, Lagos-Quintana M, Decarolis PL, et al. Subtype-specific genomic alterations define new targets for soft-tissue sarcoma therapy. Nat Genet. 2010;42:715–21.
Misawa Y, Kanazawa T, Mochizuki D, Imai A, Endo S, Carey TE, Mineta H. Epigenetic inactivation of galanin and GALR1/2 is associated with early recurrence in head and neck cancer. Clin Exp Metastasis. 2016;33(2):187–95.
Jones FS, Kioussi C, Copertino DW, Kallunki P, Holst BD, Edelman GM. Barx2, a new homeobox gene of the Bar class, is expressed in neural and craniofacial structures during development. Proc Natl Acad Sci U S A. 1997;94:2632–7.
Przybyl J, Kidzinski L, Hastie T, Debiec-Rychter M, Nusse R, van de Rijn M. Gene expression profiling of low-grade endometrial stromal sarcoma indicates fusion protein-mediated activation of the Wnt signaling pathway. Gynecol Oncol. 2018;149:388–93.
Lederle W, Hartenstein B, Meides A, Kunzelmann H, Werb Z, Angel P, et al. MMP13 as a stromal mediator in controlling persistent angiogenesis in skin carcinoma. Carcinogenesis. 2010;31:1175–84.
Quondamatteo F, Reinhardt DP, Charbonneau NL, Pophal G, Sakai LY, Herken R. Fibrillin-1 and fibrillin-2 in human embryonic and early fetal development. Matrix Biol. 2002;21:637–46.
Xia H, Qi H, Li Y, Pei J, Barton J, Blackstad M, et al. LATS1 tumor suppressor regulates G2/M transition and apoptosis. Oncogene. 2002;21:1233–41.
Waalkes S, Atschekzei F, Kramer MW, Hennenlotter J, Vetter G, Becker JU, et al. Fibronectin 1 mRNA expression correlates with advanced disease in renal cancer. BMC Cancer. 2010;10:503.
Demicco EG, Maki RG, Lev DC, Lazar AJ. New therapeutic targets in soft tissue sarcoma. Adv Anat Pathol. 2012;19:170–80.
Andrae J, Gallini R, Betsholtz C. Role of platelet-derived growth factors in physiology and medicine. Genes Dev. 2008;22:1276–312.
Andrick BJ, Gandhi A. Olaratumab: a novel platelet-derived growth factor receptor α-inhibitor for advanced soft tissue sarcoma. Ann Pharmacother. 2017;51:1090–8.
The authors thank Mark Dasenko, Aaron Trippe, Matthew Peterson, and Chris Sullivan of the Oregon State University Center for Genome Research and Biocomputing and Cheri Goodall, Milan Milovancev, and Shay Bracha of the Oregon State University Carlson College of Veterinary Medicine for technical assistance. We thank David Vail, Michael Huelsmeyer, and Ilene Kurzman for kindly providing the FISS cell lines.
SAR acknowledges support from the Animal Cancer Foundation. CVL acknowledges support from the National Institutes of Health (P01CA090890-06A2), the Oregon State University College of Veterinary Medicine (internal research grant), and the Oregon State University Center for Genome Research and Biocomputing. CK acknowledges support from the Brooke’s Blossoming Hope for Childhood Cancer Foundation. The funding agencies played no role in the design of the study; the collection, analysis, or interpretation of the data; or in writing the manuscript.
Availability of data and materials
All cat mRNA-seq data files for this study have been deposited in the NCBI SRA database under accession number SRP134052. The dog mRNA-seq data files are available from the corresponding author upon request.
Ethics approval and consent to participate
All animal work on this project was conducted under animal use protocols approved by the institutional animal care and use committees at Oregon State University (FISS samples) and Colorado State University (canine soft tissue sarcoma samples). Tissue samples had been originally obtained as diagnostic specimens and were banked for research purposes with owner consent.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Table S1. Sample information for the FISS transcriptome profiling. Each row corresponds to an RNA sample. Columns are as follows: sample name, the coded sample name used within this manuscript; biosample type, indicates whether the sample came from primary tissue or cultured cells; normal or cancer, distinguishes whether the sample is normal or cancerous; tissue type, gives the tissue subtype, such as sarcoma, skin, or muscle; Cat, the coded cat patient ID for this study; batch, the mRNA-seq batch number (see Methods); “#reads”, the total number of reads derived from the sample; “#uniquely mapped reads”, the number of reads that mapped to a unique location in the cat genome; %alignment, the percentage of the total number of reads that aligned to a unique location in the cat genome; #genes with non-zero count, the number of genes for which the count of aligned reads is greater than zero. (XLSX 11 kb)
Supplementary Note 1: This file contains the cDNA sequence used for the qPCR assay design for the gene LATS1. (PDF 42 kb)
Table S2. Differentially expressed genes in FISS. mRNA-seq gene expression data for 21,890 cat genes in all feline tissue and cell culture samples in this study. Columns as follows: Geneid, Ensembl identifier for the gene; Chrom, identifier of the cat chromosome (MT, mitochondrial chromosome; chromosome identifiers starting with “JH” are unplaced scaffolds) to which the gene has been mapped; Gene Start, leftmost chromosomal coordinate of the gene in the genome assembly; Gene End, rightmost chromosomal coordinate of the gene; Strand, the direction of transcription of the gene (1 denotes positive strand, − 1 denotes negative strand); Symbol, official gene symbol; Description, gene description field; Columns 8–21, absolute RNA-seq tag counts within the indicated gene (row), in the indicated sample (column); log2_fold_change_tissue, log2(sarcoma/skin) value (or “NA” if not statistically significantly different from zero); p_adj_tissue, adjusted p-value for expression level comparison between the sarcoma and skin sample groups; log2_fold_change_cultured_cells, log2(sarcoma cell / fibroblast) value (or NA if not detected as nonzero with statistical significance); p_adj_cultured_cells, adjusted p-value for expression-level comparison between the sarcoma cell line and fibroblast sample groups. (XLSX 4618 kb)
Figure S1. Concordance of qPCR measurements with mRNA-seq measurements of relative gene expression in FISS vs. skin samples. Each mark represents a gene (genes are as described in Table 1) and marks are labeled by official HGNC gene symbol. (PDF 26 kb)
Table S3. Regions of coherent up- or down-regulation in sarcoma vs. normal tissue. Each row corresponds to a 10 Mbp region of the cat genome for which the expression levels of the genes within the region are (together) consistently up- or down-regulated in sarcoma tumor tissue vs. normal skin. Columns as follows: Fc chr., cat chromosome; Fc pos. (Mbp), the chromosomal coordinate of the start of the region, in Mbp; Genes within region, the Ensembl gene identifiers of all genes that are annotated within the region; log2(sarcoma/normal), the average of the log2(sarcoma/normal) expression values for all genes in the region; Hs chr., the chromosome of the human genome region that is syntenic to the indicated cat genome region; Hs pos, the coordinates (in the GRCh38 genome assembly) of the human genome region that is syntenic to the indicated cat genome region. (DOCX 128 kb)
Table S4. Quantification of protein expression of BARX1 in fibroblasts and FISS cells. Columns as follows: “27 kDA BARX1”, integrated intensity of the BARX1 band at 27 kDa; “40 kDa Memcode”, integrated intensity of the normalization control chain at 40 kDa; “BARX1/Mem”, ratio of intensity of the band at 27 kDa to the intensity of the band at 40 kDa, for the indicated row; “FISS/fibrobl”, ratio of “BARX1/Mem” for the indicated row, to “BARX1/Mem” for the first row (“cat01 fibrob.”). Rows as follows: “cat01 fibrobl.”, fibroblasts from skin sample from cat01; “cat04 FISS”, cells derived from FISS tumor sample from cat04; “cat05 FISS”, cells derived from FISS tumor sample from cat05. (DOCX 44 kb)
Figure S2. Full membrane images for BARX1 and FN1 immunoblots. (A) Anti-BARX1 stained Western with ladder. (B) Same membrane stained for total protein with Pierce Removable Total Protein. (C) Anti-FN1 stained Western with ladder. (B) Same membrane stained for total protein with Pierce Removable Total Protein. (PDF 132 kb)
Table S5. Quantification of protein expression of FN1 in fibroblasts and FISS cells. Columns as follows: “220 kDA FN1”, integrated intensity of the FN1 band at 220 kDa; “40 kDa Memcode”, integrated intensity of the normalization control chain at 40 kDa; “FN1/Mem”, ratio of intensity of the band at 220 kDa to the intensity of the band at 40 kDa, for the indicated row; “FISS/fibrobl”, ratio of “FN1/Mem” for the indicated row, to “FN1/Mem” for the first row (“cat01 fibrob.”). Rows as follows: “cat01 fibrobl.”, fibroblasts from skin sample from cat01; “cat04 FISS”, cells derived from FISS tumor sample from cat04; “cat05 FISS”, cells derived from FISS tumor sample from cat05. (DOCX 38 kb)
About this article
Cite this article
Wei, Q., Ramsey, S.A., Larson, M.K. et al. Elucidating the transcriptional program of feline injection-site sarcoma using a cross-species mRNA-sequencing approach. BMC Cancer 19, 311 (2019). https://doi.org/10.1186/s12885-019-5501-z