Chromatin structure regulates cancer-specific alternative splicing events in primary HPV-related oropharyngeal squamous cell carcinoma

ABSTRACT

Human papillomavirus-related oropharyngeal squamous cell carcinoma (HPV+ OPSCC) represents a unique disease entity within head and neck cancer with rising incidence. Previous work has shown that alternative splicing events (ASEs) are prevalent in HPV+ OPSCC, but further validation is needed to understand the regulation of this process and its role in these tumours. In this study, eleven ASEs (GIT2, CTNNB1, MKNK2, MRPL33, SIPA1L3, SNHG6, SYCP2, TPRG1, ZHX2, ZNF331, and ELOVL1) were selected for validation from 109 previously published candidate ASEs to elucidate the post-transcriptional mechanisms of oncogenesis in HPV+ disease. In vitro qRT-PCR confirmed differential expression of 9 of 11 ASE candidates, and in silico analysis within the TCGA cohort confirmed 8 of 11 candidates. Six ASEs (MRPL33, SIPA1L3, SNHG6, TPRG1, ZHX2, and ELOVL1) showed significant differential expression across both methods. Further evaluation of chromatin modification revealed that ASEs strongly correlated with cancer-specific distribution of acetylated lysine 27 of histone 3 (H3K27ac). Subsequent epigenetic treatment of HPV+ HNSCC cell lines (UM-SCC-047 and UPCI-SCC-090) with JQ1 not only induced downregulation of cancer-specific ASE isoforms, but also growth inhibition in both cell lines. The UPCI-SCC-090 cell line, with greater ASE expression, also showed more significant growth inhibition after JQ1 treatment. This study confirms several novel cancer-specific ASEs in HPV+OPSCC and provides evidence for the role of chromatin modifications in regulation of alternative splicing in HPV+OPSCC. This highlights the role of epigenetic changes in the oncogenesis of HPV+OPSCC, which represents a unique, unexplored target for therapeutics that can alter the global post-transcriptional landscape.

KEYWORDS:

• Human papillomavirus
• oropharyngeal cancer
• alternative splicing
• super-enhancers

Introduction

Head and neck squamous cell carcinoma (HNSCC) is the sixth leading cause of cancer deaths worldwide. While these tumours have been traditionally associated with heavy tobacco and alcohol use, the risk factors are changing, and there has been a significant increase in oropharynx tumours related to human papillomavirus (HPV). The incidence of HPV-positive (HPV+) oropharyngeal squamous cell carcinoma (OPSCC) has more than doubled over the past decades making it one of the most rapidly growing cancer populations in the United States, while the incidence of HPV-negative (HPV−) OPSCC has decreased [1]. These trends are likely attributable to decreasing tobacco use as well as rising prevalence of HPV within the population [2]. Notably, within TCGA, HPV+ tumours show relatively few genetic alterations in cancer drivers at the exome level compared to HPV− tumours [3]. This may be related to fewer carcinogenic exposures to cause DNA alterations, but it may also suggest that changes in post-transcriptional regulation or epigenetic changes may take precedence in HPV-related head and neck carcinogenesis. Indeed, post-transcriptional changes have been shown to have an important role in head and neck cancer in mediating oncogenesis, tumour progression and drug resistance [4–10]. Furthermore, the epigenetic landscape (on the level of both DNA methylation and histone modification) of HPV+ HNSCC samples correlates with the expression of cancer-related genes and cancer prognosis [10–13]

Alternative splicing is a post-transcriptional process that allows for diversity in protein production in the absence of genetic alterations. Previous literature has shown that variations in the splicing of a wild-type (WT) isoform of a gene can contribute to carcinogenesis. For example, alternatively spliced forms of genes BCL2L1 and TNR6 have been shown to be either pro or anti-apoptotic in human cancers, based on exon inclusion and expression [14]. Within head and neck cancer, multiple genes, such as LAMA3, DST, and TP63 are alternatively spliced [8]. Recent work has also revealed functionally active alternatively spliced variants of AKT3 and GSN in HPV+ OPSCC [9,15]. In addition, over 100 additional splicing events were identified in a previous publication, but their function and regulation remain unstudied [9].

Epigenetic DNA methylation has been well established as a method of gene expression regulation in tumours including in the head and neck. In HNSCC, epigenetic silencing of tumour suppressors have been shown to play a role in tumorigenesis [5–7], and histone modifications may mediate drug resistance [4]. Moreover, a growing body of literature also has implicated epigenetic changes in splicing regulation. In other tumour types, it has been observed that DNA methylation is higher in exons and intergenic regions compared to introns, suggesting that DNA methylation may regulate alternative splicing [16,17]. It has also been proposed that the generation of aberrant isoforms with alternative starts is a potential mechanism by which cancer cells circumvent suppression of oncogene transcription by hypermethylation [9,18,19].

Prior literature has established the interplay between chromatin state and gene splicing. Histone methylation marks are significantly enriched at exons, and studies have shown H3K4me3 to recruit spliceosomal factors through adaptor proteins [20]. In addition, histone acetyltransferases such as the SAGA protein complex interact with and recruit splicing machinery, and facilitate recognition of intron branchpoints [21]. In addition, these histone modifications associated with open chromatin, including H3K27-acetylation, are associated with enhancers which regulate the transcriptional activity of multiple neighbouring gene targets, up to 1–1.5 megabases away, through the formation of chromatin loops [22,23]. These super-enhancers (SEs) are defined by recruitment of tissue-specific transcription factors, H3K27 acetylation, BRD4 and MED1 proteins [24].

Herein, we propose that epigenetic modification in HPV-related OPSCC, specifically regions of chromatin, enriched by H3K27ac, regulate expression of alternative splicing variants which contribute to oncogenesis in these tumours. To evaluate this hypothesis, we first performed validation of 11 alternatively spliced genes from a list of previously identified candidate 109 ASEs [9] to confirm both their presence in OPSCC primary tumours and the validity of ASEs identified in silico. We then utilize these ASEs to evaluate the correlation between the epigenetic landscape and ASE expression. HPV+ head and neck cell lines were subsequently treated with JQ1, the inhibitor of BRD4, which recognizes the H3K27ac histone mark. JQ1 treatment resulted in alteration of splice variant balance as well as growth inhibition in multiple head and neck cell lines. This is the first study to demonstrate the role of H3K27ac landscape and SE positioning in the epigenetic regulation of ASEs in HPV+ HNSCC.

Results

Validation of 11 candidate ASEs by qRT-PCR and with TCGA

A total of 11 ASEs were selected from 109 previously published ASEs [9] for qRT-PCR validation (Table 1). The ASEs selected for validation were among the most prevalent ASEs identified in the tumour cohort, and which also showed correlation with DNA methylation and chromatin state, as previously published [9]. The 11 ASEs included 8 alternative start events, 2 deletion events, and one insertion event. All ASEs except ELOVL1 exhibited outlier over-expression in the tumour group relative to normal tissues. Of the 11 ASEs, only 4 were canonical (CTNNB1, MKNK2, SYCP2, and ZNF331). Figure 1 shows the sashimi plots of 2 representative HPV+ OPSCC and 2 non-cancer samples for the ZHX2 (Figure 1(a)) and MKNK2 (Figure 1(b)) ASEs.

Table 1. Eleven alternative splicing events selected for validation .

Gene	Full name		Coordinates	Type of splicing event	Canonical
GIT2	GIT ArfGAP 2	chr12	110,418,984–110,421,415	Insertion	No
CTNNB1	Catenin beta 1	chr3	41,280,845–41,281,310	Deletion	Yes
MKNK2	MAP kinase interacting serine/threonine kinase 2	chr19	2,037,828–2,040,133	Deletion	Yes
MRPL33	Mitochondrial ribosomal protein L33	chr2	27,975,814–28,002,300	Alternative start	No
SIPA1L3	Signal-induced proliferation associated 1-like 3	chr19	38,423,332–38,519,729	Alternative start	No
SNHG6	Small nucleolar RNA host gene 6	chr8	67,834,960–67,838,433	Alternative start	No
SYCP2	Synaptonemal complex protein 2	chr20	58,496,505–58,497,445	Alternative start	Yes
TPRG1	Tumour protein P63 regulated 1	chr3	189,027,966–189,028,175	Alternative start	No
ZHX2	Zinc fingers and homeoboxes 2	chr8	123,875,561–123,875,716	Alternative start	No
ZNF331	Zinc finger protein 331	chr19	54,058,142–54,059,099	Alternative start	Yes
ELOVL1^a	ELOVL fatty acid elongase 1	chr1	43,831,294–43,831,939	Alternative start	No

^aASE under-expressed in tumours

Genomic coordinates, as aligned to GRCh37/hg19, are shown of the junction correlating with the alternative splicing event identified.

Figure 1. Two validated splice variants ZHX2 and MKNK2 ASEs.

Alternative start site in ZHX2 (A.) and the deletion event in MKNK2 (B.) Sashimi plots of representative controls (N = Normal, light) and representative tumours (T = Tumour, dark) are shown along with the designed WT and ASE primer-probe assays. Sashimi plots display downsampled RPM counts using IGV viewer of given junction transcripts. Boxplots show the results of ASE qRT-PCR validation in n = 22 UPPP non-cancer controls and n = 31 HPV+ OPSCC samples. Mann–Whitney U test revealed significant overexpression of ASE (normalized to WT expression) within tumours for both splice variants (p-values <0.001).

In vitro validation by qRT-PCR revealed 9 out of the 11 genes to have significantly higher rates of outlier ASE expression (FDR < 0.05) in tumour samples compared to normal samples (Table 2). TPRG1 ASE had the highest prevalence of outlier overexpression in tumour samples at 94% (n = 29 out of 31), followed by MRPL33 at 90% (n = 28 out of 31), and MKNK2 at 77% (n = 24 out of 31, Figure 1(d), FDR<0.00001). Notably, for most ASEs investigated, qRT-PCR showed higher sensitivity for detection of ASEs with a higher prevalence of outlier expression compared to RNA-Seq within this cohort (Table 2).

Table 2. Validation of 11 ASEs in vitro (qRT-PCR) and in silico (TCGA).

		Discovery			Validation
		RNA-seq JHU cohort (nT = 47, nN = 25)			PCR (in vitro) JHU cohort (nT = 31, nN = 22)			RNA-Seq (in silico) TCGA cohort (nT = 44, nN = 16)
		nT (%)	nN (%)	FDR	nT (%)	nN (%)	FDR	nT (%)	nN (%)	FDR
Associated gene	GIT2	19 (41.4)	1 [4]	6.92E-04	14 (45.2)	1 (4.6)	1.70E-04	5 (11.4)	0 (0)	3.11E-01
	CTNNB1	18 (39.2)	1 [4]	1.50E-03	6 (19.4)	1 (4.6)	1.04E-01	25 (56.9)	0 (0)	4.56E-05
	MKNK2	29 (63.1)	1 [4]	6.19E-07	24 (77.5)	1 (4.6)	1.12E-09	7 (16)	0 (0)	1.73E-01
	MRPL33	26 (56.6)	1 [4]	6.59E-06	28 (90.4)	1 (4.6)	2.29E-13	31 (70.5)	0 (0)	4.54E-07
	SIPA1L3	15 (31.9)	0 (0)	6.52E-04	16 (51.7)	0 (0)	2.46E-06	22 (29)	0 (0)	1.78E-04
	SNHG6	24 (52.2)	0 (0)	2.39E-06	19 (61.3)	0 (0)	4.26E-09	28 (63.7)	0 (0)	4.22E-06
	SYCP2	18 (39.2)	0 (0)	1.16E-04	7 (22.6)	1 (4.6)	5.29E-02	35 (79.6)	0 (0)	1.37E-08
	TPRG1	29 (63.1)	0 (0)	2.74E-08	29 (93.6)	1 (4.6)	6.41E-14	12 (27.3)	0 (0)	2.53E-02
	ZHX2	20 (43.5)	1 [4]	3.30E-04	22 (71)	1 (4.6)	2.12E-08	31 (70.5)	0 (0)	4.54E-07
	ZNF331	21 (45.7)	0 (0)	1.81E-05	17 (54.9)	0 (0)	6.19E-07	0 (0)	0 (0)	1.00E+00
	ELOVL1^a	33 (71.8)	1 [4]	1.28E-08	8 (25.9)	1 (4.6)	2.26E-02	40 (91)	0 (0)	3.24E-11

n_N is the number of normal samples.

n_T is the number of tumour samples.

^aASE underexpressed in tumours.

Bold indicates significant (FDR < 0.05).

RNA-Seq analysis of ASE prevalence was originally performed on the discovery ‘JHU cohort’ (left-most column) with 25 non-cancer controls (n_N) and n = 47 HPV+ HNSCC (n_T). In vitro validation by qRT-PCR using custom assays (Table S1) was performed on a partial JHU cohort with available tissue for analysis (n_T = 31 and n_N = 22), and in silico RNA-seq validation was performed on an independent TCGA cohort (n_T = 44 and n_N = 16). All p-values are adjusted comparing ASE outlier prevalence between the normal and tumour groups. All ASEs except ELOVL1 exhibited outlier over-expression in the tumour group relative to normal tissues.

Cross-validation was then performed in silico on an independent cohort (TCGA) of 16 normal samples and 44 HPV+ OPSCC. Within the TCGA data set, 8 out of the 11 ASEs had significant differential outlier ASE expression (FDR < 0.05) between the tumour and normal samples (Table 2). The ASEs with the highest prevalence in the tumour group in the TCGA cohort were ELOVL1 at 91% (n = 40 out of 44), followed by SYCP2 at 80% (n = 35 out of 44). ZNF331 ASE was not identified as an outlier in any of the tumour samples. In addition, when comparing ASE expression among all three platforms, outlier detection in normal samples remained consistently at a frequency below 5% (1 sample or less, Table 2). All 11 candidate ASEs were validated either in vitro or in silico. Consistent differential expression was found in six genes (MRPL33, SIPA1L3, SNHG6, TPRG1, ZHX2, and ELOVL1).

Spatial correlation between chromatin modifications and ASE

Next, GenometriCorr (projection statistics) analysis was performed to evaluate the potential relationship between ASE expression and chromatin modification. This analysis showed that the spatial distribution of the 109 ASEs significantly correlated with tumour-specific H3K27ac super-enhancer (SE) regions as defined by ROSE analysis (Figure 2(a), p < 0.001) but not with normal-specific H3K27ac SEs. Interestingly, GenometriCorr analysis also -showed that tumour-specific SEs were significantly correlated with the location of alternative start sites (p = 0.021), as well as both canonical (p = 0.011) and non-canonical alternatively spliced tumour-related gene isoforms (p = 0.005), while no correlation was seen with normal tissue super-enhancers (p > 0.2) (Supplemental Table 2). One example of H3K27ac enrichment in relationship with a splicing event is shown in the TPRG1 gene (Figure 2). ROSE analysis identified nearby SE regions, enriched with H3K27ac, within the 1 million bp region surrounding TPRG1 in PDX and cell lines, but not normal tissue (highlighted in Figure 2(b)). Within this region, the presence of the TPRG1 ASE expression is confirmed by RNA-Seq analysis only in samples (SCC090 & X2) that demonstrated H3K27ac enrichment with absence in normal tissue. Notably, both cell lines had relatively low expression of canonical WT TPRG1.

Figure 2. Spatial correlation of chromatin reorganization and alternative splicing events.

a. GenometriCorr adjusted p-values showing the significant overlap between tumour-specific H3K27ac chromatin modification super-enhancer regions compared to splicing events in tumours which is not seen in normal samples. b. Super-enhancer enrichment in the 1 Mb region surrounding TPRG1 ASE, enriched with H3K27ac (highlighted area) and correlates with the presence of the TPRG1 ASE within RNA-Seq data (c) shown in a zoomed-in view of the specific splicing junction with the TPRG1 gene.

In addition to H3K27ac, three other chromatin marks (H3K4me3, H3K9ac and H3K9me3) were studied in relationship to alternative splicing expression. These data show that active marks H3K27ac and H3K4me3 both show consistent spatial enrichment in regions of tumour-related ASEs, only in tumour-specific marks. H3K9ac showed spatial correlation that was not specific to normal or tumour tissue, an effect that was observed also when relating to gene expression in previously published work [10]. H3K9me3 is a repressive mark, and as expected did not show any significant relationship to ASE expression in either normal or tumour-specific marks.

Next, we evaluated the spatial relationship of enrichment of these histone marks to specific ASE types (Supplemental Table 2). Interestingly, in all active histone marks, there was a significant spatial correlation between histone marks and non-canonical ASEs. Next, alternative start site ASE isoforms were evaluated and all three histone marks showed tumour-specific spatial correlation with alternative start site splicing events. As previously published work has shown, the active histone marks (H3K27ac, H3k4me3, and H3K9 c) demonstrate enrichment in transcription start site regions which may relate to these alternative start site isoforms. Analysis of relationship of these chromatin marks for insertion and deletion events was deferred, given the lower prevalence of these splicing events.

Epigenetic-driven changes in ASE landscape decrease cell proliferation

Two HPV+ cell lines (UM-SCC-047 and UPCI-SCC-090) were then treated with the bromodomain inhibitor, JQ1, which is expected to downregulate the activity of SEs [25]. Out of 109 ASE previously detected for HPV+ OPSCC, 26 ASEs were detected in UM-SCC-047 and 60 ASEs in UPCI-SCC-090 using RNA-Seq analysis. When compared to the number of ASEs in tumours, UM-SCC-047 was within the lower quartile and UPCI-SCC-090 was in the upper quartile (Figure 3(a)). JQ1 suppressed proliferation in both UPCI-SCC-090 and UM-SCC-047 (p < 0.0001 and p < 0.005, respectively, Figure 3(b,c). However, the growth suppression effect of JQ1 was more pronounced in ASE-rich cell line, UPCI-SCC-090.

Figure 3. JQ1 treatment of HPV+ HNSCC cell lines.

(a) Alternative splicing event (ASE) count in each of the 47 samples of discovery cohort, along with the ASE count in UM-SCC-047 and UPCI-SCC-090 cell lines as shaded in black. Samples where ASE count ≥ median of 46 ASEs per tumour are shown in dark grey, with UM-SCC-047 in the lower half (ASE count of 26) and UPCI-SCC-090 in the upper half (ASE count of 60) (b and c) Treatment with JQ1 resulted in a decrease in cell proliferation in 090 (p < 0.0001) and 047 (p < 0.005). (d and e) Fold change of ASE and WT isoform expression in response to JQ1: A fold change greater than 1 indicates an increase, while a fold change less than 1 indicates a decrease. Error bars denote standard deviation.

JQ1 treatment resulted in changes in the splicing balance of several of the 11 genes studied, producing varying effects on the expression of the ASE and WT isoforms. In the UM-SCC-047 cell line, nearly all genes showed a decrease in ASE expression after JQ1 treatment, while wild-type gene expression was not significantly changed in five genes, and decreased in six of those studied. In UPCI-SCC-090, again a majority of genes (8 of 11) showed decrease in splice variant expression after JQ1 treatment. In this cell line, JQ1 had more variable impact on wild-type gene expression, with four genes showing decreased expression, and the remainder with stable or increased expression. In addition, the splice variant ELOVL1, which was suppressed in tumours, showed enhanced expression after JQ1 treatment in this cell line. It is also noted that the splice variant of MRPL33 showed enhanced expression after JQ1 treatment in both cell lines.

RNA-Seq analysis of both UM-SCC-047 and UPCI-SCC-090 cell lines was conducted comparing JQ1 treatment to control (DMSO). Among 109 ASE candidates [1], 75 splice junctions were captured by RNA-seq with adequate coverage involving 68 unique genes for which overall gene expression was available. The UM-SCC-047 cell line had relatively low overall expression of ASEs, with 24 (32%) ASEs showing no expression in both DMSO or JQ1 treated cells. There were 12 junctions with significant differential expression after JQ1 treatment, and 83% (n = 10) showed decreased expression after JQ1 treatment (Supplemental Table 3). However, all but one of these (ASE in BAI2) were accompanied by significant decreases on overall gene expression after JQ1 treatment. In contrast, after JQ1 treatment of UPCI-SCC-090, 91% (n = 10 of 11) of junctions with significant differential expression showed decreased expression. Half of these (n = 5) were accompanied by either no change or increase in overall gene expression, showing splicing-specific changes induced by JQ1 treatment (Supplemental Table 4).

To better understand the specificity of JQ1 treatment, we performed JQ1 treatment in an additional HPV+ HNSCC cell line (93-VU-147 T), as well as immortalized oral keratinocyte cell lines (OKF6 & NOKSI). Significant growth inhibition was seen in all cell lines treated with JQ1. As in UM-SCC-047 and UPCI-SCC-090, JQ1 also significantly decreased a majority of ASEs (9 of 11) in 93-VU-147 T and wild-type expression remained primarily unchanged or increased (7 of 11 genes, Supplemental Figure 1). In immortalized normal cell lines (Supplemental Figure 2), greater off-target effects were seen in response to JQ1 treatment. There was greater inhibition of wild-type gene expression in immortalized oral keratinocytes, with 72–81% of genes of interest showing decrease in wild-type expression with JQ1 treatment in the non-cancer cell lines. Notably, JQ1 did decrease expression of ASEs as well; however, only about half of the studied ASEs were able to be detected in the normal cell lines, even in untreated controls. Therefore, growth inhibition in the normal cell lines may be more related to off-target impact on gene expression.

Evaluation of specificity to HPV+ disease

To understand whether validated ASEs were specific to the HPV+ population of the disease, additional analysis was performed on an HPV− HNSCC cohort from TCGA RNA-sequencing analysis. Normalized junction expression for the 11 identified ASEs was compared between 411 tumours and 34 normal samples. This analysis revealed that about half of the ASEs (n = 6) were also significantly differentially expressed in HPV− tumours (CTNNB1, MKNK2, MRPL33, SNHG6, SYCP2, and ZHX2, FDR<0.05, Supplemental Table 5), suggesting that not all ASEs are specific to the HPV+ sub-type. However, five ASEs did not have a significant differential expression in HPV− tumours and may be more specific to HPV-related tumours (GIT2, SIPA1L3, TPRG1, ZNF331, and ELOVL1).

Furthermore, we sought to understand splicing events in the context of HPV viral integration. Of tumours with detectable HPV16 genome, 16 tumours showed reliably detectable viral integration, 19 did not show evidence of integration and 10 were equivocal with limited HPV−human fusion detected. We then compared the rate of alternative splice variant expression, utilizing the previously published 109 candidates [9], in tumours with integrated HPV viral genome (n = 16) and non-integrated tumours (n = 19). The number of splice variants expressed in integrated tumours was significantly higher than non-integrated tumours (mean 55.1 vs. 41.3 splicing events, Student’s T-test, p = 0.001, Supplemental Figure 3).

Discussion

The scarcity of mutational alterations in HPV+ OPSCC has limited our understanding of the underlying biology of this disease and the detection of oncogenic targets. This study therefore seeks to elucidate the post-transcriptional process of alternative splicing in HPV+ OPSCC through the characterization of cancer-specific ASEs and investigation of their epigenetic control. In other disease processes, alternative splicing of the WT isoform of a gene has been shown to contribute to carcinogenesis [18–20]. For example, alternatively spliced forms of genes BCL2L1 and TNR6 are known to be either pro or anti-apoptotic in human cancers, based on exon inclusion and expression [12]. Similarly, we have recently demonstrated that cancer-specific splicing in AKT3 and GSN leads to oncogenic splice variants in head and neck SCC which promote growth and invasion [9,15].

Alternative splicing variants allow for variant protein forms with oncogenic potential to be expressed in tumours, even in the absence of DNA mutations. However, therapeutic targeting of individual splice variants may be challenging unless the mechanisms allowing for their expression in tumours is better understood. What is the process by which the balance between splice variants is determined; and what alterations in tumours allow for expression of tumour-specific splicing events? The data from this study suggest that chromatin modifications may contribute to the regulation of aberrant cancer-specific splice variants.

In order to study the impact of chromatin modifications on ASE expression, we selected 11 ASEs from a list of 109 previously published ASE candidates for additional validation. These selected candidates offer a representative sample of 10% of the ASEs identified by previously published algorithm, with high prevalence and correlation with DNA methylation [9]. Within validation, the use of outlier statistic methods with multiple testing correction for detection of aberrant splicing events is extremely stringent, only counting those samples with outlier expression at least two standard deviations above average of any normal samples. Despite this stringency, we identified six ASEs (MRPL33, SIPA1L3, SNHG6, TPRG1, ZHX2, and ELOVL1, as having significant tumour-specific expression in both biological evaluation and within TCGA.

These validated splice variants occur in genes of potential biological interest. Upregulation of SNHG6 expression predicts poor prognosis in colorectal cancer [26] and poor clinical features in gastric cancer [27]. Similarly, ZHX2 expression correlates with worse prognosis in renal cell carcinoma [28]. These publications suggest the potential biological role of ASEs in HPV+ OPSCC requiring further in-depth functional analysis that is outside of the scope of the current manuscript. However, multiple recent studies have shown the functional role of tumour-specific splicing variants that contribute to oncogenesis within head and neck cancer. In prior work, a functionally active splice variant of AKT3 was identified, in which its impact on growth in HNSCC was isoform specific [9]. Other splicing variants have recently been reported showing isoform-specific functional activity in head and neck cancer, including an isoform of GSN which promotes migration [15] and an LOXL2 isoform which promotes growth and progression [29].

Interestingly, some of the genes have also been studied in the context of HPV integration. A study by Parfenov et al. found HPV integration to occur on chromosome 3, and that it is associated with chromosomal rearrangements leading to the amplification of TPRG1 along with TP63 [30]. TP63 has been a frequent site of HPV integration in both oropharynx and cervical cancers and play a potential role in SE regulation in those tissues, and the ASE in TPRG1 was found to be specific to HPV-positive tumours  [10,31,32]. The MRPL33 ASE is particularly of interest because HPV integration has been noted to occur in the vicinity of MRPL33 – within a 3 million base pair region surrounding the gene; however, the MRPL33 ASE was also seen in HPV− tumours [33], and so may not be exclusively driven by HPV integration. While not all ASEs were specific to HPV+ tumours, we found that tumours with integration of the viral HPV genome exhibited significantly more ASEs.

While epigenetics, such as chromatin modifications and DNA methylation, are widely known to contribute to oncogenesis, mounting evidence over recent years has pointed to an interrelation between chromatin and epigenetic modifications with regulation of alternative splicing in disease. The histone acetyltransferase Gcn5 was demonstrated to interact with two U2 snRNP proteins and its HAT activity is required for their recruitment to pre-mRNA branch points [21]. Chromatin compaction is known to have a kinetic effect on differential splicing in which weaker splice sites can be detected [34]. Likewise, the subunit Brm of the SWI/SNF chromatin remodeller regulates splicing by pausing transcription, thereby promoting inclusion of exons with weaker splice sites [35]. The epigenetic modification of H3K4me3 regulates splicing through binding protein CHD1 recruiting of spliceosome component U2 snRNP. CHD1 knockdown reduced splicing efficiency and combined knockdown with H3K4me4 further reduced U2 snRNP recruitment [20].

In this study, H3K27ac modifications were found to significantly spatially correlate with ASEs, suggesting a stronger link between alternative splicing and histone acetylation than previously understood. H3K27-acetylation is associated with super-enhancers, that can regulate and activate gene expression in large groups of genes through changes in chromatin loop structures and recruitment of tissue-specific transcription factors to significantly change intracellular biology [22,23]. In prior work, tumour-specific gene expression changes were correlated with both active histone marks H3K27ac as well as H3K4me3 [10]. Similarly in this study, these chromatin marks showed tumour-specific enrichment related to alternative splicing events not seen in normal tissue. These relationships were strongest when relating to ASEs which were non-canonical, and ASEs relating to alternative start sites, which may be due to enrichment in transcriptional start sites which has been previously demonstrated [10].

A variety of evidence implicates bromodomain containing protein 4 (BRD4), an acetylated histone-binding protein, in splicing regulation [36]. BRD4 has been associated with U2AF65 activity and the production of mature spliced transcripts, and its mutant form leads to a reduction in U1 snRNP recruitment at intron-containing genes. In the context of this study, assays of the ASEs and their associated WT junction revealed a clear change in splice variant balance after treatment with bromodomain inhibitor JQ1. Both cell lines showed that a majority of the tumour-specific ASEs were downregulated after JQ1 treatment to a greater degree than wild-type isoforms, including SNHG6, ZHX2, TPRG1 and SIPA1L3 which were validated across all methods.

Treatment with JQ1 significant decreased cell growth of both HPV+ head and neck cancer cell lines UM-SCC- 047 and UPCI-SCC-090. There are known off-target growth effects of JQ1 which may be related to JQ1’s suppression of MYC [37]. It is notable that growth suppression was twice as great in UPCI-SCC-090 compared to UM-SCC-047, and UPCI-SCC-090 had a greater degree of aberrant ASE expression. In addition, JQ1 had a more specific suppression of aberrant ASE expression compared to wild-type gene in UPCI-SCC-090 in both PCR and RNA-Seq.

Previous studies [9,15] have shown that ASEs within HPV+ OPSCC (AKT3 and GSN) are functionally relevant and confer oncogenic properties. Given the changes in ASE expression seen after JQ1 treatment, some of the greater growth inhibition seen in the ASE-rich UPCI-SCC-090 cell line may be due to the alteration of its oncogenic ASE expression. Super-enhancers represent exciting therapeutic targets in cancer, with small molecule inhibitors such as JQ1 with the potential to affect global tumour oncogene expression [38]. Chromatin-regulating small molecules such as histone deacetylase inhibitors (vorinostat) and hypomethylating agents (azacitidine and decitabine) have already been FDA approved for treatment of haematologic malignancies [39,40], and recent evidence suggests a role for super-enhancer in the oncogenesis of head and neck cancer [10,41].

One limitation with JQ1 treatment is that it demonstrated greater off-target effects in immortalized oral keratinocytes, compared to treated tumour cell lines, with alteration in wild-type genes. In comparison, in tumour cell lines, expression alteration was more specific to alternative splice site. Indeed, prior work has shown that these histone modifications may also modulate disease-specific gene expression in additional to splicing variants [10]. Therefore, the growth suppression seen from JQ1 treatment may be related to both ASE changes and off-target effects.

We acknowledge the additional limitation of this study in methodology and scope. In the PCR validation of splice variants, most PCR primers were designed de novo, so non-specific annealing of primers may impact reliability and sensitivity of detection of specific splicing events. In addition, there was observed to be variation between the different validation approaches in the prevalence of tumour outliers. These variations may be explained by differences in cohorts and/or techniques. Within PCR validation, several ASEs showed higher prevalence than in the RNA-Seq data, which may relate to its ability to detect specific splice junctions in an unbiased manner dictated by the depth of sequencing. With regard to JQ1 treatment experiments, as discussed previously, off-target effects with impact on gene expression may also impact growth effects observed, particularly in non-tumour models. The changes in ASE expression and balance were observed in only a subset of ASEs consisting of the selected 11 ASEs and changes to the full ASE landscape were not interrogated in all cell lines.

In conclusion, eleven splicing variants were validated in HPV+ OPSCC in vitro, and in silico on an independent TCGA cohort. Six ASEs were confirmed as differentially expressed between normal and tumour samples across methodologies, involving genes which have evidence of potential oncologic impact based on prior literature. In addition, HPV+ head and neck cancer cell lines treated with JQ1 demonstrate how epigenetic-driven changes can both alter the splice variant landscape of these tumours, and potentially inhibit tumour growth, which was more prominent in an ASE-rich cell line. Overall, this study provides novel evidence that ASEs are one mechanism by which the activity of relevant proteins is altered in HPV+ OPSCC and are under regulation by epigenetic changes like histone acetylation and SE machineries influenced by HPV activity. Understanding the ASE landscape of HPV+ OPSCC provides an opportunity for the discovery of novel molecular targets which alter the epigenetic landscape which may result in changes in ASE expression and subsequent tumour inhibition.

Patients and methods/materials

Study cohorts and patient-derived xenografts

The initial discovery cohort was comprised of data from a recently published study of 47 primary tumour tissue samples from patients with HPV+ OPSCC and 25 normal mucosal tissue from uvulopalatopharyngoplasty (UPPP) surgical samples in non-cancer-affected patients [9]. All tissue samples were collected under an institutional review board-approved protocol (#NA_00036235). Due to limited sample RNA quantity, a partial subset of this discovery cohort (n = 22 UPPP non-cancer controls, and n = 31 HPV+ OPSCC samples) was used for the biological validation. Two patient-derived xenografts (PDX) models were generated from two primary tumour samples from this discovery cohort [10]. From The Cancer Genome Atlas (TCGA) [42], additional data for a separate cohort of 44 HPV+ samples (with 16 corresponding normal oropharyngeal tissues) and 411 HPV− samples (with 34 corresponding normal tissues) was obtained for in silico validation. Junction quantification and RSEM gene expression from RNA-sequencing data were obtained from Broad GDAC firebrowse (gdac.broadinstitute.org) utilizing mRNASeq from Illumina HiSeq platform for head and neck squamous cell carcinoma.

Cell culture and drug treatment

Three HPV+ head and neck cell lines were utilized: UM-SCC-047, UPCI-SCC-090 [43], and 93-VU-147 T provided by Dr Thomas Carey (University of Michigan) and Dr Susanne Gollin (University of Pittsburgh). Cell lines were authenticated using the Short Tandem Repeat (STR) Profiling Service by the Johns Hopkins School of Medicine Genetic Resources Core Facility DNA Services. Cells were grown on high-glucose DMEM media (Clontech, Mountain View, CA), supplemented by 10% fetal bovine serum (FBS) and 1% Penicillin-Streptomycin at 37°C in 5% CO2 [10]. Additionally, immortalized oral keratinocyte cell lines (OKF6 and NOKSI) were cultured and drug treated in a similar manner.

HPV+HNSCC and immortalized oral keratinocyte cell lines were grown in 6 well tissue culture dishes and were treated with 500 nM JQ1 (bromodomain inhibitor, Selleck Chemicals, Houston, TX, USA) or with 0.1% of dimethyl sulphoxide (DMSO), as a control, added to high-glucose DMEM media for 72 h at 37°C. The cellular growth-monitoring experiments were performed in pentaplicates using 4-h incubation with 1:10 diluted Alamar Blue (Bio-Rad) in high-glucose DMEM Media measured every 24 h.

RNA isolation, sequencing and qRT-PCR

RNA from treated and non-treated cell lines, as well as primary samples, was extracted using mirVana miRNA Isolation Kit (Ambion, Forster City, CA). RNA-Seq data for this cohort were previously published and analysed (GSE112027) [9,10,44]. MapSplice2 version 2.0.1.9 was used to quantify read counts during the alignment [44,45]. Gene expression values were quantified from RNA-sequencing data using RNA-Seq by Expectation-Maximization (RSEM) version 1.2.9 [46].

Due to limited sample RNA quantity, a partial subset of this discovery cohort (n = 22 UPPP non-cancer controls, and n = 31 HPV+ OPSCC samples) was used for qRT-PCR validation. All RNA samples were reverse transcribed into cDNA using the High Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Forster City, CA). QRT-PCR primers and probes of the ASEs and their WT controls were designed using Integrated DNA Technologies PrimerQuest tool, with the probe spanning across the ASE junction or its associated WT junction (Supplemental Table 1). These were validated by touchdown PCR (data not shown). All qRT-PCR experiments were normalized to GAPDH expression and performed on an Applied Biosystems 7900 Real-Time PCR system in triplicates.

H3K27ac chromatin immunoprecipitation (ChIP), ChIP-Seq analysis, and the detection of H3K27ac super-enhancers

H3K27ac-specific ChIP-Seq data for 2 PDX samples (PDX1 and PDX2) and 2 UPPP samples (UPPP1 and UPPP2) from the cohort described above, and the UM-SCC-047 and UPCI-SCC-090 cell lines, were previously published and analysed (GSE112027) [10]. Chromatin was extracted from these six samples, as previously described [10]. The H3K27ac (8173, CST) antibody was used to isolate DNA segments bound by this histone modification, using SimpleChIP Plus Enzymatic Chromatin IP Kit (9005, CST). MACS (Model-based Analysis of ChIP-Seq algorithm, version 1.4.2) called ChIP-Seq peaks for each sample using the input DNA in that sample as a control [47]. MACS peaks for H3K27ac were further ranked with the ChIP-Seq signal intensity and stitched if necessary by the Ranking Of Super Enhancers (ROSE) analysis software [23] for each chromatin mark and for each sample. The resulting set of merged peaks with heavy ChIP-Seq abundance is referred to as enhancers. The top ranking of these were classified as super-enahcers (SE).

Detection of heavy ChIP-Seq abundance for H3K27ac, H3K4me2, H3K9ac, and H3K9me3

In order to evaluate whether other histone marks have spatial correlation with ASE expression, we completed similar wet-lab and dry lab analysis for additional active histone marks (H3K4me3, and H3K9ac), and repressive mark (H3K9me3). Just as in the case of H3K27ac, we performed chromatin immunoprecipitation of H3K4me3 (9751, CST), H3K9ac (9649, CST), and H9K9me3 (13,969, CST) using the same kit (9005, CST). Similarly, high-throughput sequencing and MACS analysis were completed for these additional marks. To recapitulate H3K27ac analysis, we employed the same ROSE algorithm to call regions of heavy ChIP-Seq abundance peaks, which are normally called ‘super-enhancers’ and ‘typical enhancers’ by ROSE in the case of H3K27ac. However, we intentionally added ‘typical enhancer’ calls for non-enhancer marks H3K4me3, H3K9ac, and H3K9me3, as these would not strictly be referred to as ‘super-enhancers.’ This analysis was repeated for all four marks including H3K27ac in parallel. The ROSE analysis was completed for each mark and each sample.

GenometriCorr spatial correlation analysis

The R package GenometriCorr [48] was used to compare the positional distribution of the 109 ASEs with signal-enriched peaks for 4 chromatin marks (H3K27ac, H3K4me3, H3K9me3, and H3K9ac). The peaks were detected by ROSE software, and the term ‘super-enhancers’ was defined based on prior publication methods for all four chromatin marks.(REF) The package is intended to test the positional correlations (co-localization) of pairs of genome-wide annotations, which are represented as the sets of genomic intervals. For each of the four marks, super-enhancer datasets were categorized into normal-specific and tumour-specific events, via removal of ROSE-called peaks detected in both tumour and normal samples, and the analysis was performed separately for normal-specific and tumour-specific SEs. P-values were adjusted using the Bonferroni correction based on the number of tests performed.

Outlier and median statistical analyses

Differential ASE expression between tumour and normal groups was evaluated with outlier analysis using the OGSA Bioconductor package [49], on ASE expression normalized to gene expression. Samples were categorized as outliers based on expression distribution in normal samples, and presence of outliers was compared between groups using Fisher’s exact test, and FDR adjusted p-values utilizing the Benjamini–Hochberg method [9]. qRT-PCR data were performed in triplicate with two-tailed Mann–Whitney test [50] at a significance level of 0.05

HPV integration

HPV integration status of tumours was determined based on RNA-sequencing analysis using MapSplice, as previously published [10]. In short, to detect integration, a reference genome was created as a chimera including the human genome and HPV16 genome and tumours were mapped to this chimera genome. In this way, viral integration was detected as a fusion site between the human and HPV16 genome. We considered the viral genome integrated if there were at least three discordant pairs (one paired-end read mapped to the viral genome, and its mate pair mapped to the human genome) and one split read (read spanned the human-viral junction). These seven total reads support integration at the same locus [10,30].

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplemental Materials

Supplemental data for this article can be accessed here

Additional information

Funding

This work was supported by the Johns Hopkins University [Catalyst Award]; National Institutes of Health (NIH) [T32 DC000027]; NIH [P30CA006973]; NIH [5P50DE019032]; Russian Foundation for Basic Research (RFBR) [17-00-00208]; Russian Academic of Sciences [0112-2019-0001]; American Academy of Otolaryngology-Head and Neck Surgery (AAO-HNSF) [Resident Research Grant]; NIH [R21DE025398]; NIH [4R50CA243627]; NIH [R01DE027809].