Chromatin run-on sequencing analysis finds that ECM remodeling plays an important role in canine hemangiosarcoma pathogenesis

Background Canine visceral hemangiosarcoma (HSA) is a highly aggressive cancer of endothelial origin that closely resembles visceral angiosarcoma in humans, both clinically and histopathologically. Currently there is an unmet need for new diagnostics and therapies for both forms of this disease. The goal of this study was to utilize Chromatin run-on sequencing (ChRO-seq) and immunohistochemistry (IHC) to identify gene and protein expression signatures that may be important drivers of HSA progression. Results ChRO-seq was performed on tissue isolated from 17 HSA samples and 4 normal splenic samples. Computational analysis was then used to identify differentially expressed genes and these factors were subjected to gene ontology analysis. ChRO-seq analysis revealed over a thousand differentially expressed genes in HSA tissue compared with normal splenic tissue (FDR < 0.005). Interestingly, the majority of genes overexpressed in HSA tumor tissue were associated with extracellular matrix (ECM) remodeling. This observation correlated well with our histological analysis, which found that HSA tumors contain a rich and complex collagen network. Additionally, we characterized the protein expression patterns of two highly overexpressed molecules identified in ChRO-seq analysis, podoplanin (PDPN) and laminin alpha 4 (LAMA4). We found that the expression of these two ECM-associated factors appeared to be largely limited to transformed endothelial cells within the HSA lesions. Conclusion Outcomes from this study suggest that ECM remodeling plays an important role in HSA progression. Additionally, our study identified two potential novel biomarkers of HSA, PDPN and LAMA4. Interestingly, given that function-blocking anti-PDPN antibodies have shown anti-tumor effects in mouse models of canine melanoma, our studies raise the possibility that these types of therapeutic strategies could potentially be developed for treating canine HSA.


Background
Angiosarcomas (AS) are highly aggressive malignant tumors originating from endothelial cells. They account for approximately 2% of all soft tissue sarcomas in humans with the number of cases increasing significantly over the past 30 years [1][2][3][4][5]. For patients presenting with non-metastatic AS, the reported five-year survival rate is~35%. These patients also have a 75% chance of local recurrence within 24 months and a 50% likelihood of metastases developing despite local treatment [4,6]. When metastases are already present at the initial presentation, the five-year survival rate is poor, with a median survival time of just 3 months [4]. Although published research on this rare tumor is increasing, we still know very little about the pathogenesis of this disease.
Canine hemangiosarcoma (HSA) is histopathologically similar to angiosarcoma, with both forms of this disease following a similar clinical course [7]. Canine HSA is most commonly observed in the spleen, but also occurs in other organs such as the heart, liver, and dermis with the later form often being associated with ultraviolet radiation-associated oncogenesis. Prognosis for the visceral forms of HSA is poor, with most dogs dying from this disease within months of their diagnosis [8,9]. For dogs diagnosed with splenic HSA, surgical removal of the spleen can increase the life expectancy up to 6 months, and when combined with chemotherapy, can prolong their life for up to 12 months [8,9].
While all breeds are susceptible to hemangiosarcoma, German Shepherds, Golden Retrievers, and Labrador Retrievers, are particularly predisposed to developing this disease [10], with the estimated lifetime risk of Golden Retrievers developing HSA being 20%. This breed predisposition suggests a genetic component for HSA [11]. Genome-wide association studies have been carried out to identify risk loci [12,13]. SNP array analysis of genomic DNA from Golden Retrievers with HSA identified a risk locus on chromosome 5 that was shared by~20% of the cases [12]. In another study, microarray analysis of genomic DNA from dogs with HSA revealed gained copy number aberrations on several genes such as PDGFRA, KDR and VEGFA [13]. A recent whole exome sequencing study of dogs with HSA across a range of breeds revealed somatic mutations in tumor suppressor genes, including TP53, and two genes (PIK3CA and PIK3R1) in the Phosphoinositide 3-kinase (PI3K) pathway [14]. Comparison of somatic copy number aberration profiles in human angiosarcoma and canine hemangiosarcoma identified recurrent copy number gains in KDR (31% in human, 22% in canine) and KIT (17% in both) [15]. These genome-wide studies reveal that, while specific genetic aberrations are associated with HSA in some populations, these alterations are not sufficient to explain the majority of HSA cases. These observations support the hypothesis that HSA pathogenesis is heterogeneous in nature.
Transcriptome analysis has been previously performed on cell lines and tumor tissues to identify molecular features that define canine HSA. Gene expression profiling of HSA cell lines and non-malignant proliferating endothelial cells revealed that HSA cell lines appeared to overexpress genes associated with inflammation, angiogenesis, adhesion, invasion, metabolism, cell cycle progression, and patterning [16]. Microarray and RNA-seq analysis of visceral HSA tumors identified three distinct molecular subtypes; angiogenesis, inflammation and lipogenesis. These molecular subtypes did not appear to be associated with a specific breed or tumor morphologic subtype [17].
A variety of analytic tools exist for assaying the transcriptome and molecular alterations in tissue. One of these tools, chromatin run-on and sequencing (ChROseq), uses RNA polymerase activity to measure transcription and, as such, provides for a highly-sensitive, base-pair level readout of gene expression [18,19]. Given the lack of clarity regarding the molecular underpinnings of HSA, the major goal of this study was to utilize ChRO-seq to document changes in gene expression between normal splenic tissue and HSA tumor tissue. Our results show that the majority of genes that are upregulated in HSA appear to be associated with extracellular matrix (ECM) remodeling. Additionally, we further characterize two ECM-associated molecules that were highly overexpressed in HSA tumor tissue, podoplanin (PDPN) and laminin alpha 4 (LAMA4). We show by immunohistochemistry (IHC) that the expression of these two cancer-associated factors appears to be primarily limited to HSA tumor cells.

ChRO-Seq analysis of transcription in HSA and normal splenic tissue
For our genome-wide study, we performed ChRO-seq analysis on tumor tissue from dogs histopathologically diagnosed with HSA (20) and on normal splenic tissue (4) ( Table 1). We quantified the similarity of transcription between these tissue samples using Pol II abundance in annotated gene bodies. Three HSA samples (B297, B675, B788) were excluded due to a low number of mappable reads (< 2 million reads). Sequencing data from the remaining samples was then analyzed to create a Spearman's rank correlation matrix (Fig. 1a). Seven thousand seven genes (> 20 rpkm) were used to calculate the correlation coefficients by GENE-E. Four normal tissues and 3 HSA tissues (B307, B829 and C349) were found to form one cluster, with 2 HSA samples (C001 and C034) also being similar to the Normal/HSA cluster.
DESeq2 analysis was then performed on the samples to identify differentially expressed genes. When compared to normal splenic tissue, a total of 906 genes were upregulated in the HSA samples (FDR < 0.005) while 358 genes were downregulated (FDR < 0.005). An MA plot showing the differentially expressed gene set is shown in Fig. 1b. In this plot, the X axis represents the average expression over all samples and the Y axis represents the log2 fold change between the tumor samples and normal spleen. Highly upregulated extracellular matrix (ECM)-associated genes and downregulated heme synthesis-associated genes are labelled with blue dots (Fig. 1b). Differentially expressed genes (FDR < 0.001, 659 genes) were then assessed to evaluate expression levels in individual samples (Fig. 1c). Two hundred and ninety-one genes (rpkm > 20) are shown in the resulting heatmap which also contains an associated hierarchical clustering dendrogram (one minus spearman rank correlation with complete linkage). Results show that the expression profiles of HSA samples appear to be distinct from normal splenic tissue. Extracellular matrix remodeling appears to be a major feature of HSA When evaluating the top 50 most highly upregulated genes in the HSA group, we observed that many of these molecules are associated with the ECM, examples of which include podoplanin (PDPN), laminin alpha 4 (LAMA4), ADAM12, and fibronectin 1 (Fig. 1b). To test the prevalence of ECM-associated molecules in our dataset, we carried out gene ontology (GO) analysis using PANTHER (protein analysis through evolutionary relationships) (Fig. 2). We utilized a total of 369 of the most highly upregulated genes (FDR < 0.0005) as an input list, with 282 of these genes then being identified as uniquely mapped genes in a Canis lupus familiaris reference database. A total of 20 annotation categories were detected (FDR < 0.05) from the reactome pathway database, with the top 10 annotation categories from this dataset being found to be related to ECM regulation. A more detailed list of the ECM-associated factors that were upregulated in HSA tumor tissue is shown in Table 2, along with their fold change (log 2 ) and FDR. A full list of Reactome pathways terms and the list of genes is shown in the supplemental table. Interestingly, when we queried genes in the downregulated dataset (124 genes with an FDR cut-off of 0.0005), we did not observe statistically significant functional associations using GO analysis. Therefore, we adjusted the FDR cut-off to 0.005 and this adjustment identified 358 genes as being downregulated in the HSA samples with 301 of these genes being uniquely mapped in the reference database. As shown in Fig. 2, five annotation categories were detected in this dataset, with the most enriched reactome pathway term being heme biosynthesis. Examples of specific genes that are associated with heme biosynthesis include, UROS, FECH, UROD, JCHAIN, and ALAD (Fig. 1b). This result suggests that the expression of genes involved in normal splenic function is suppressed in HSA tumor tissue.

Characterization of ECM-associated proteins in HSA tumor tissue
To more directly test the hypothesis that ECMassociated factors are overexpressed in HSA tumor tissue, we first stained paraffin sections of tumor and normal tissue with Masson's trichrome to label collagen fibers (Fig. 3). Staining of the more vascular regions of the HSA tumor tissue show that collagen fibrils are abundant throughout these regions. This observation fits well with our finding that multiple collagen isoforms were highly overrepresented in our HSA-overexpression ChRO-Seq dataset. Interestingly, in these vascular regions, the malignant tumor cells appear to surround the collagen bundles, forming what resembles inverted blood vessels. To further support these findings, images and We next further characterized two ECM-associated factors that were highly overexpressed in HSA tumor tissue. For this analysis we focused on PDPN (FDR 6.43E-05, fold change 3.62) and LAMA4 (FDR 1.58E-09, fold change 3.77) because these molecules have been previously associated with cancer progression in other tumor types [20][21][22][23][24][25]. For further characterization, we first investigated PDPN and LAMA4 expression levels in our ChRO-seq dataset. A genome browser view of transcription levels from one normal and three HSA samples is shown in Fig. 4a. Expression levels from all ChRO-seq samples are shown in supplemental figure S2. Results show that, at the PDPN locus, two of the HSA cases (C085 and C356) showed increased RNA polymerase activity compared to normal tissue, while expression levels of the third HSA sample (C340) was comparable with normal tissue. LAMA4 was transcribed from the minus strand and all HSA samples showed higher transcriptional activity than normal tissue. We next investigated the expression of mature PDPN and LAMA4 transcripts in tumor tissue and normal spleen using RT-PCR. Results show that, while PDPN or LAMA4 transcripts were either absent or weakly expressed in normal spleen, LAMA4 was significantly expressed in all three tested HSA tumor samples while PDPN was robustly expressed in one (C356) of the three HSA samples (Fig. 4b and supplemental Figure S7). Next, we performed IHC analysis of PDPN and LAMA4 on HSA tumor sections to confirm that these molecules are expressed in tumor tissue at the protein level and to investigate their localization. HSA samples B297, B783 and B848 are shown in Figs. 5 and 6 with images and quantitative data from additional samples being shown in supplemental Figure S3 and S4. IHC analysis found that anti-PDPN staining in normal splenic tissue was low to absent. In HSA tumor tissue, however, anti-PDPN staining was robust in the cytoplasm of tumor cells encircling the collagen fibrils ( Fig. 5c and d,  arrow). In addition to tumor cells, anti-PDPN staining was also evident in the cytoplasm of endothelial cells in tumor-adjacent blood vessels (Fig. 5b, arrowhead). The localization pattern of anti-LAMA4 staining was similar to that of PDPN, with little staining being seen in normal splenic tissue and strong anti-LAMA4 staining being observed in the cytoplasm of HSA tumor cells encircling collagen fibrils. Again, as with PDPN, anti-LAMA4 staining was also observed in the cytoplasm in the endothelial cells from tumor-adjacent blood vessels ( Fig. 6c and d).
Serial sections of additional HSA samples were then stained with Masson's trichrome and immunohistochemistry was performed on these additional tissues using anti-PDPN and LAMA4 antibodies (Supplemental figure S5 and S6). Sections from one sample (B176, Figure S5) showed solid/cavernous features with blood vessels (Figure S5 a). Masson's trichrome staining highlighted the presence of collagen fibrils filling spaces between regions of neoplastic cells ( Figure S5 a and b). Anti-PDPN staining appeared to be primarily cytosolic in these neoplastic cells and in endothelial cells lining adjacent blood vessels (Figure S5 c). Anti-LAMA4 staining was strong in extracellular regions of the tumor and Sections from sample B554 primarily showed capillary/ cavernous features ( Figure S6). Masson's trichrome staining of this sample revealed a core collagen region surrounded by malignant endothelial cells with spindle shaped nuclei ( Figure S6 a and b). Anti-PDPN and LAMA4 staining appeared to be primarily cytosolic in tumor cells from this sample.

Discussion
In this study, we first utilized ChRO-seq to identify RNA polymerase activity in hemangiosarcoma tumor tissue and normal splenic tissue. Analysis of the correlation matrix for our ChRO-seq dataset finds that 14 of the HSA tumor samples clustered together, while 3 of the HSA samples appeared to be more similar to normal splenic tissue (Fig. 1). While speculative, it is possible . GAPDH primers were used as a positive control. Negative control (NC) lacks cDNA template. Results show that PDPN and LAMA4 mRNA do not appear to be expressed in normal spleen. LAMA4 transcripts were observed in all three HSA samples while PDPN showed strong expression in one of the three HSA samples that these 3 samples represent a different subtype of HSA and, as such, displayed a gene expression profile that is more similar to normal tissue. Alternatively, it is also possible that these 3 presumptive HSA tissues were actually tumor-adjacent normal tissue or a mix of normal and tumor tissue. While grading of canine hemangiosarcoma is not often utilized due to the aggressive nature of the neoplasm, a grading scheme does exist and differences in gene expression between the samples may correlate with progression and eventual outcome [26].
Previous HSA transcriptomic profiling studies have implicated a number of signaling pathways in HSA pathogenesis. Tamburini et al., for example, found that, when compared to cell lines derived from splenic hematomas, HSA cell lines exhibited several different distinct gene expression profiles, including signatures for angiogenesis, inflammation, adhesion, invasion metabolism, cell cycle and signaling. In another study, Gorden et al. performed microarray and RNA-seq analyses on visceral HSA tumor tissues (12 spleen, 7 heart, 4 liver, 1 lung) and identified three distinct tumor subtypes that were associated with either angiogenesis, inflammation, and adipogenesis. Our finding that ECM remodeling appears to be a major gene expression signature in HSA differs from these previous studies. One possible reason for the difference in outcomes between our study and the Tamburini study is that our report compared gene expression patterns in HSA splenic tissue with normal splenic tissue while the Tamburini study was comparing the gene expression profiles of HSA cell lines with splenic hematoma cell lines [16]. Thus, the differences may be due to the fact that the expression profile of normal splenic samples differs significantly from hematoma samples. Additionally, these differences may also be because the expression profile of cell lines may differ significantly from primary tissue samples due to decreased cellular complexity and prolonged cell culture. One potential reason why outcomes from our study differed from the Gorden study is that we evaluated differences in gene expression patterns between HSA and normal splenic tissue while the other study analyzed gene expression patterns within HSA tissues. Lastly, another reason why outcomes from our study may have differed from both of these previous studies is that ChRO-seq analysis detects nascently synthesized RNA as opposed to mature transcripts, whose levels can be affected by a variety of factors including transcript stability.
The extent to which ECM genes were upregulated in HSA tumors in our study is highlighted by our GO analysis which found that 9 of the top 10 biological process categories were ECM-related (Fig. 2). Interestingly, 5 of these 9 categories related to collagen function. We more directly tested for the abundance of collagen fibers in the tumor samples using Masson's trichrome stain and found extensive collagen deposition throughout the tumor tissue (Fig. 3). Collagen deposition was observed both in the more solid areas of the tumor and, in particular, in the tumor regions filled with malformed vascular channels. In these vascular regions, neoplastic endothelial cells are often found to encircle the collagen bundles forming "gumball"-shaped structures that look like inverted blood vessels. While stromal fibroblasts are presumably primarily responsible for the synthesis of these collagen fibers, it is also possible that the tumor cells may also be partly responsible for synthesis of these collagen fibers.
In addition to collagen, we also found that PDPN and LAMA4 were highly overexpressed in HSA tumor tissue when compared to normal tissue. These molecules were of particular interest to us given their close association with cancer progression in other types of cancer. Podoplanin is a mucin-type protein that contains an extracellular region, transmembrane domain, and intracellular tail. It is widely known as a marker for lymphatic endothelial cells and also to play a critical role in heart and lung development and in development of the lymphatic endothelial system [43][44][45]. PDPN appears to play several roles in cancer progression. A number of studies have shown that PDPN expression in cancer cells promotes tumor cell proliferation and invasion [20][21][22]. In addition to cellintrinsic effects, PDPN is also believed to promote tumor metastasis by interacting with its receptor, CLEC2, on the Fig. 6 Immunohistochemical localization of LAMA4 protein in HSA and normal splenic samples. a LAMA4 detection was low to absent in normal splenic tissue (B297, normal white pulp). b anti-LAMA4 staining appears weak in endothelial cells of blood vessels that are adjacent to tumor masses (B783). c and d anti-LAMA4 staining is strong in tumor cells within the tumor's vascular regions (B848). In these regions the LAMA4positive tumor cells are seen to surround what appears to be collagen bundles. LAMA4 was also detected in tumor cells from the more avascular regions of the tumor (Supplemental materials). c = low magnification, d = high magnification. Scale bar represents 50 μm platelets. This interaction then promotes the coating of tumor cells by platelets, thereby protecting tumor cells from the immune system [46]. In addition to its role in human cancers, PDPN is also overexpressed in canine squamous cell carcinomas and melanomas [47] and PDPN mAbs were recently found to have potent anti-tumor activity in mouse xenograft models of canine melanoma [48]. Interestingly, overexpression of PDPN in mice leads to disseminated intravascular coagulation [49], a condition that is strikingly similar to that seen in many dogs with HSA [50]. Our ChRO-seq data demonstrated that PDPN expression varies significantly between HSA samples ( Figs. 1 and 4). IHC analysis found that PDPN was robustly expressed in transformed endothelial cells in certain HSA tumor samples. However, similar to our ChRO-seq finding, we did find that PDPN appears to only be expressed in a subset of HSA tumors. In future studies, we plan to test whether PDPN expression in HSA tumors correlates with disease severity.
Laminins form large heterotrimeric αβγ protein complexes and are a prominent component of basement membranes. LAMA4 is distinct from other laminin isoforms in that it lacks a polymerization domain with the loss of this domain potentially facilitating tumor cell migration [23]. Interestingly, LAMA4 was recently described as an "oncolaminin" due to its strong association with cancer cell migration and tumor progression in a range of cancers [24]. These links to cancer include recent studies which found that LAMA4 and MCAM (melanoma cell adhesion molecule) are highly enriched in tumor blood vessels in renal cell carcinoma and colorectal carcinoma. Additionally another study found that expression of these molecules is enhanced in locally invasive and metastatic clear cell renal cell carcinoma [25]. Further, antisense oligonucleotides against laminin-8 (LAMA4 and LAMB1) were found to block the invasion of glioma cells and neovascularization in vitro [51]. Taken together, these published studies indicate that LAMA4 plays a key role for vascular development, tumor progression and metastasis. In our study, deseq2 analysis found that LAMA4 is highly transcribed in HSA tumor tissue. This observation was supported by PCR analysis of mRNA isolated from HSA tumor tissue ( Fig. 4 and supplemental Figure S7). Additionally, our IHC analysis found that LAMA4 protein expression appears to be primarily limited to malignant endothelial cells (Fig. 5). Given the previously defined roles for LAMA4 in cancer progression and tumor metastasis, our findings raise the possibility that LAMA4 is an important mediator of canine HSA pathobiology.

Conclusions
Outcomes from our studies found that the majority of upregulated genes in HSA tumor tissue appear to be associated with ECM remodeling. This finding was supported by IHC analysis which found a robust collagen network throughout HSA tumor tissue. Additionally, we further characterized two ECM-associated factors, PDPN and LAMA4, and found that their expression was largely limited to tumor cells in HSA tissue. Both of these molecules have been previously identified as potential biomarkers in other types of cancer. Given the recent previous preclinical mouse studies showing that anticanine PDPN antibodies can block the growth of melanomas, our findings raise the possibility that similar types of therapies may have utility for treating canine HSA.

Samples
Canine tissue samples were obtained from the Canine Comparative Oncology & Genomics Consortium and Cornell University Hospital for Animals. Pathology was independently confirmed, and patient demographics are described in the supplemental material (Table 1). ChROseq libraries were made from splenic hemangiosarcoma (n = 20) and normal splenic tissues from HSA dogs (n = 4). Normal splenic samples from three beagles were either snap-frozen or fixed in PFA and Paraffin-embedded (FFPE) for sectioning. Additional FFPE sections of HSA cases were obtained from Cornell University Hospital for Animals.

Chromatin-run-on and sequencing
Chromatin was extracted from each tissue sample and chromatin run-on was performed as described previously [18,19,52]. ChRO-seq library preparations were executed according to the Illumina protocol and were sequenced using Illumina NextSeq500 sequencing. Raw sequence FASTQ files were assessed by FastQC for quality control [53]. Single-read sequencing data were preprocessed and mapped to the canine genome assembly using the proseq 2.0 pipeline (https://github.com/Danko-Lab/proseq2.0) [54]. Briefly, Single-read sequencing data were preprocessed to remove the adapter sequences and trimmed based on base quality, and deduplicating the reads based on unique molecular identifiers in RNA adapters. Sequencing reads were mapped into a canine genome assembly (CanFam 3.1) using the Burrows-Wheeler Aligner software package [55]. The number of mappable reads are listed in Table 1. Aligned BAM files were converted into bigWig format, which were used to create the matrix of read counts mapping to each annotated gene (CanFam3 ensGene). Seven thousand seven genes (> 20 rpkm) were used to calculate the Spearman's rank correlation coefficients to create the matrix and dendrogram by GENE-E (https://software.broadinstitute.org/GENE-E/). The list of differentially expressed genes and an MA-plot were produced by running DESeq2 (false discovery rate < 0.01) [56]. Upregulated genes (369 genes, FDR < 0.0005) and down regulated genes (358 genes, FDR < 0.005) in HSA were subjected to PANTHER (Protein Analysis Through Evolutionary Relationships) overrepresentation gene ontogeny analysis (GO) on Reactome pathway dataset (version 65, released 2019-12-22) [57]. Differentially expressed genes (FDR < 0.001, 659 genes) were assessed to show the expression level in individual samples. Heatmaps and dendrograms were created by GENE-E.

RT-PCR
Snap frozen tissues (3 normal spleens and 3 HSA tissues, see Table 1 for demographics) were pulverized in a tissue crusher and total RNA was extracted with Trizol reagent. 1 μg of RNA was reverse-transcribed to cDNA using High-capacity cDNA reverse transcription kit (Applied Biosystems) following the manufacture's protocol (25°C for 10 min, 37°C for 120 min, 85°C for 5 min). All cDNA reactions were diluted 5-fold with water prior to use in PCR. 1 μl of diluted cDNA was used for PCR using Go-Taq master mix (Promega Corporation) with canine specific primers (0.4 μM each, sequences are shown in supplemental material), products were amplified using the following parameters: 95°C for 5 min, 30 cycles of 3 steps (95°C for 30 s, 58°C for 30 s, and 72°C for 20 s) and 72°C for 5 min. The amplified products were analyzed using agarose gel electrophoresis (1.2%).