An orthogonal cell line with H3.3-G34W provide means for deep hnRNPA1L2 analysis

To analyze the intersection between hnRNPA1L2 and H3.3-G34W, we used the previously generated clones of HeLa cells that were gene-targeted to the H3F3A locus with a GFP-tagged H3.3-G34W (denoted isoH3.3-G34W). Parental HeLa cells or wild type isoH3.3-WT served as controls. Additionally, we used our panel of giant cell tumor of bone (GCTB) primary samples with their endogenous WT or G34W of H3.3. RNA-sequencing data from individual clones distinctly separated WT and G34W samples (Fig. 1a, c.f. Supplementary Fig. S1a). To determine the effect of H3.3-G34W on gene expression, we conducted hierarchical clustering of differentially expressed genes (DEG) by comparing GCTB stromal cells of mesenchymal origin with HeLa cells of epithelial origin. The results showed that approximately half of the DEGs shared between HeLa and GCTB had a similar expression pattern (Fig. 1a). The observed dichotomy is likely due to differences in the cells of origin, and H3.3 appears to contribute to DEGs depending on the cell type. A total of 1614 DEGs were statistically significant for GCTB, and 1276 for HeLa, with 209 genes shared between the two (see Fig. 1b). As per our previous study, we found that the DEGs in GCTB were targeted by the E2F transcription factor family [4]. Based on this, we estimate that 150 out of the 209 shared genes (72%) were targeted by E2F. The list of E2F target genes is available at the Broad Institute: Human Gene Set: HALLMARK_E2F_TARGETS and Encode. The E2F family of transcription factors (TFs) has been implicated in cancer due to its association with cell cycle regulation as a member of the CDK-RB-E2F nexus that promotes cell cycle entry through an E2F transcriptional program [15]. Hierarchical clustering indicated that the 209 genes targeted by E2F included both up-and down-regulated genes with approximately half (46%) being regulated similarly between GCTB and HeLa, and the other half (54%) showing an opposite pattern between the two cell types (Fig. 1c). A volcano plot indicates that isoH3.3-G34W results in both upregulated and downregulated genes (Fig. 1d). This suggests that H3.3-G34W governs a gene expression program of over one thousand genes, and that the overlap between the two diverse cell types converges to a substantial degree over E2F transcription factor regulation.

Fig. 1: Gene expression analysis by RNA-seq of GCTB vs. HeLa cells harboring H3.3-G34W.

a The results of hierarchical K-means clustering indicate distinct variations in DEGs. The heatmap depicts 12 GCTB (4 H3.3-WT and 8 -G34W) and 6 isogenic HeLa RNA-seq clones (three independent isogenic H3.3-WT and -G34W clones). b A comparison of WT vs. G34W in GCTB and HeLa extracted 209 overlapping DEGs of which 72% are E2F target genes. The total number of DEGs (adj. p < 0.05, |log2FC| > 1) is tabulated. c Hierarchical K-means clustering separates DEGs based on H3.3-G34W in GCTB and HeLa cells. The up- and down-regulated genes indicate cellular specificity of frequent E2F targets. d The volcano plot shows that there are 22% more down-regulated genes than up-regulated genes. e Multidimensional scaling was used to display the degree of separation between HeLa clones with WT and G34W genetic backgrounds, using a variance stabilizing transformation and showing the first two components. The knockout clones of hnRNPA1L2 cluster separately from either WT or G34W.

Knockout of hnRNPA1L2 is distinct from its functional role in H3.3-G34W

To investigate the role of hnRNPA1L2 and its link to H3.3-G34W, hnRNPA1L2 was knocked out using CRISPR/Cas9-mediated gene targeting in isoH3.3-WT (Supplementary Fig. S1b). RNA-sequencing analysis was performed on the knockout (designated hnRNPA1L2-KO) in isoH3.3-WT and compared to non-knockouts in its parental isoH3.3-WT and -G34W. The three samples were distinctly separated, as judged from principal component analysis (Fig. 1e). The results suggest that the G34W follows a distinct trajectory from the loss of function of hnRNPA1L2. However, there may be a partial overlap, which we tested in subsequent experiments.

The hnRNPA1L2 protein is an RNA-binding auxiliary factor that works with other hnRNPs to perform various functions, mainly by binding exons as an exon choice suppressor [6]. To compare isogenic HeLa cells with GCTB harboring G34W [11, 12], we performed RNA-seq and splicing analysis using the rMATS and SpliceR computational tools. The analysis of RNA-seq data from H3.3-G34W-containing GCTB and HeLa cells using rMATS revealed that exon skipping and inclusion (SE) were the most common splicing events when compared to WT. These events accounted for 60–70% of the total events (Fig. 2a). Interestingly, knockout clones of hnRNPA1L2 in isoH3.3-WT cells resulted in fewer SE events and instead showed more pronounced mutually exclusive exon (MXE) events (as indicated by the arrows in Fig. 2a). This suggests that the SE events observed in H3.3-G34W align with hnRNPA1L2 without altering its function. As shown in Fig. 2b, isoH3.3-G34W exhibited a higher frequency of SE events compared to the control group. This may be due to a redistribution of hnRNPA1L2 in H3.3-G34W samples, resulting in an increase in the number of SE events. Given our previous findings of global epigenetic alterations in GCTB [5], we further investigated whether these increased SE events were affecting epigenetic modifiers. It seems that certain epigenetic modulators and DNA modifiers were significantly altered, as shown in Supplementary Fig. S2. Similar to our previous analysis of H3.3-G34W in GCTB, isoH3.3-G34W also displayed significant changes in splicing patterns. These changes may reflect the broad impact of hnRNPA1L2, which includes essential gene regulators and chromatin remodelers.

Fig. 2: Splicing events in GCTB recapitulated in HeLa with H3.3-G34W vs. WT.

a SpliceR and rMATS analyses reveal that exon skipping/inclusion (SE) is the most common significant event. Note that hnRNPA1L2 knockout clones show a drop in SE events (arrow) and an increase in mutually exclusive events (MXE). b The plot shows significant splicing events that satisfied both FDR < 0.05 and |ΔPSI| > 0.1, where ΔPSI (percent spliced-in) is calculated by subtracting the average inclusion level of WT from the average inclusion level of G34W based on triplicate samples. Gray dots represent non-significant events, while red and blue dots represent significant events. c The rMAPS plot, based on rMATS output, identifies hnRNPA1L2 motifs as a significant motif at upregulated SE events (red) in the intron/exon junction of the target exon and the 3’end of the downstream intron in both GCTB and isogenic HeLa cells. The smallest p value is indicated with a yellow box (*p < 0.05, **p < 0.01), and a dashed line at 1.3 indicates a p value of 0.05 with significant peaks marked yellow. d A principal component analysis plot of HeLa cells was created using MASER based on the rMATS output of significant SE events. The plot successfully separates WT from G34W, which also display clonal heterogeneity.

Evidence of a significant overlap between hnRNPA1L2 and H3.3-G34W at target sites

We investigated whether the SE events in H3.3-G34W cells had hnRNPA1L2 binding motifs. To do this, we performed rMAPS analysis of the rMATS output to generate RNA splicing maps of RNA-binding proteins with defined binding motifs [16]. The hnRNPA1L2 protein binding to the canonical motif [AGT]TAGGG[AT], exhibited significant binding activity at the junction of target exons and downstream introns in both the GCTB and the HeLa transcriptomes. This finding further supports the association between H3.3-G34W and hnRNPA1L2, suggesting a potentially overlapping functional trajectory (see Fig. 2c). Out of the 115 RNA-binding proteins in the rMAPS test, hnRNPA1L2 was the third most significant RBP with a binding motif at and around the target exon (smallest p value in target exon 0.00784). The rMATS output in a multidimensional scaling plot also showed that H3.3-G34W clones dispersed while the H3.3-WT remained clustered (see Fig. 2d), suggesting that isoH3.3-G34W influenced exon target choice.

The analysis of splicing in both GCTB and isogenic HeLa with H3.3-G34W focuses on SE events. In both cases, there is significant evidence of hnRNPA1L2 binding in and around the target exon. Knockout of hnRNPA1L2 reduces SE events and favors mutually exclusive exon events.

RNA immunoprecipitation indicates extensive overlap between hnRNPA1L2 and H3.3-G34W

To determine the extent of overlap between hnRNPA1L2 and H3.3-G34W during transcription and RNA processing, we performed native, non-crosslinked RIP-seq [7] to immunoprecipitate and sequence hnRNPA1L2-associated RNA from isoH3.3-G34W and control cells. Our main objective was to identify the target RNA and broadly infer a link to the RNAs associated with H3.3-G34W. In this context, the ratio of hnRNPA1L2-peaks to the total number of reads was approximately four times higher in isoH3.3-G34W than in the parental control. This suggest that hnRNPA1L2 has novel G34W-assigned targets (Fig. 3a). In contrast, the RIP-peaks from GFP-tagged H3.3-G34W are less pronounced (rightmost bar in Fig. 3a). In accordance with the observed increase in SE events and the presence of distinct hnRNPA1L2 motifs, the hnRNPA1L2-peaks from the parental control strongly overlapped with both isoH3.3-G34W and its hnRNPA1L2-precipitates (see Fig. 3c). It is worth noting that isoH3.3-G34W (anti-GFP) precipitates a significantly larger number of transcripts than hnRNPA1L2 alone, but when used together, this number increases considerably (see Fig. 3c). The overlap between isoH3.3-G34W RIP-peaks and parental HeLa amounted to 80%, with a noticeable increase in isoH3.3-G34W peaks (Fig. 3d). Of the 19,023 H3.3-G34W peaks that overlapped with hnRNPA1L2-peaks at annotated genomic areas, 84% were found at gene bodies and only 2.2% at distal intergenic areas. Almost all of these peaks (96%) were over exons (Fig. 3e). In contrast to hnRNPA1L2 motifs, which appear to occur in both introns and exons, the estimate for their occurrence has not been quantified and may primarily be at exons (see Figs. 2c and 3e). We investigated whether there is an overlap between genes with SE events and genes with RIP-seq peaks. We found that 952 genes in isoH3.3-G34W had splicing events, and of those, 432 genes had RIP-seq peaks, which accounts for approximately 45% of aberrantly spliced genes with SE events. The estimated equivalent for GCTB was 42%, suggesting functional overlap between hnRNPA1L2 and H3.3-G34W.

Fig. 3: Close correlation of RIP-seq peaks of isoH3.3-G34W and hnRNPA1L2 in HeLa.

a The tally of peak numbers over the number of reads after RIP-seq with anti-hnRNPA1L2 antibody shows a four-fold increase in hnRNPA1L2 peak numbers in the presence of H3.3-G34W. The GFP-tagged isoH3.3-G34W in the third bar serves as a control. b Multidimensional scaling of RIP-seq peaks separated by genetic background. c Venn diagrams indicate strong overlap of parental hnRNPA1L2 RIP-seq peaks with both GFP and hnRNPA1L2 IPs in isoH3.3-G34W clones pointing to their similar targeting regimens. d The Venn diagram displays the overlap between H3.3-G34W and hnRNPA1L2 on the genome, as indicated by RIP-seq peaks in isogenic H3.3-G34W HeLa with IP of hnRNPA1L2 and GFP-tagged H3.3-G34W. e The annotated regions from the 19,023 RIP-seq peaks from H3.3-G34W indicate that peak location is primarily at gene bodies and exons. f The GO-term analysis of RIP-peaks aligns with RNA-seq around E2F target genes, cell cycle control, and RNA metabolism. The numbers indicate fraction of genes identified in each term.

The analysis of GO-terms for the 19,023 RIP-peaks once again indicated an association with E2F target genes and cell cycle regulation. This was previously suggested among the 209 overlapping DEGs from the RNA-seq analysis comparing HeLa and GCTB (refer to Figs. 1c and 3f). The close association of genomic targets between hnRNPA1L2 and isoH3.3-G34W, both intragenic and centered around E2F targets, suggests that the two proteins share a common trajectory and may contribute to tumorigenesis by dysregulating E2F function.

hnRNPA1L2 functionally aligns with H3.3-G34W at known and novel genes

In a previous study of isogenic HEK293 cells, hnRNPA1L2 was identified as a novel interacting partner to G34W after immunoprecipitation followed by mass spectrometry [4]. The significant overlap between hnRNPA1L2 and H3.3-G34W (Fig. 3a, c), the increased SE events that significantly carry hnRNPA1L2 motifs (Fig. 2a, c), and the redistribution of events in hnRNPA1L2-KO (Fig. 2a), suggests that the two molecules intersect and may share functional avenues related to GCTB tumorigenesis. The study aimed to determine if RIP-seq peaks were directed toward novel or known G34W targets. In this study, it was discovered that transcripts in the isoH3.3-G34W samples were more frequently associated with hnRNPA1L2 compared to parental cells. This was exemplified by the E3 ubiquitin ligase UBR4, which increased 18-fold (Fig. 4a and Supplementary Fig. 3a, b). Additionally, novel transcripts not targeted in parental cells, such as the E3 ubiquitin ligase HECTD4 were identified. However, some transcripts lost hnRNPA1L2 association in isoH3.3-G34W. It was hypothesized that genes with a high exon count would be more frequently affected, but an analysis of the correlation coefficient R analysis did not support this assumption (see Supplementary Fig. 3). The number of RIP-peaks unique to H3.3-G34W as seen in Fig. 4a was three times higher than in the parental group (red markers: 4094 genes, blue markers: 1281 genes), while the proportion of SE events in those genes was similar (5.8% vs. 6.3%, respectively). This indicates that more peaks lead to more events, indicating that hnRNPA1L2 is directed to new targets via its association with H3.3-G34W.

Fig. 4: Gene-by-gene analysis of hnRNPA1L2 RIP-seq peaks.

a The scatter plot displays the delta hnRNPA1L2 peak numbers of H3.3-G34W minus parental. The plot indicates peaks only detected in parental (blue dots), peaks only in H3.3-G34W (red dots), and peaks in genes of both parental and H3.3-G34W (black dots). For example, NBPF20 is a novel target that appears only in H3.3-G34W. b A comparison between splicing events from RNA-seq and RIP-seq events in H3.3-G34W in HeLa (red squares) and GCTB (blue square) over E2F targets. The 52 common genes were identified by matching the extracted SE events to annotated genes. c Waterfall plots show the inclusion level difference of these genes in GCTB. Ten genes in GCTB overlapped between RIP-seq peaks and SE events (i.e., coordinates over exons) and are marked in red text, while the remaining genes did not overlap over the same exons. d Of the 52 genes that showed perfect overlap between hnRNPA1L2 and SE events in both GCTB and HeLa, they were enriched predominantly over gene bodies. e Candidate genes from the 52 gene list showing overlap between SE events and RIP peaks. IGV traces of ZZ3, ACTN1, and KCTD17 show delta SE events of 0.03 in GCTB and −0.1 in HeLa, both at the false discovery rate of 0.01. f The pie-charts illustrate changes in isoforms in GCTB. Note that the parental does not have hnRNPA1L2 in these candidates.

An E2F nexus of cell cycle control emanates among candidate genes

To identify candidate genes targeted by hnRNPA1L2 in the H3.3-G34W background, we extracted E2F transcription factor targeted genes with RIP peaks that also had skipped exon (SE) events overlapping between HeLa and GCTB. We started with the SE events and hnRNPA1L2 peaks in HeLa, resulting in 389 genes to compare with 267 genes in GCTB. Together, they fulfilled the criteria of having hnRNPA1L2 peaks, being targeted by E2F, and displaying SE events. This comparison resulted in 52 genes that overlapped between HeLa and GCTB (Fig. 4b). The waterfall plot in Fig. 4c displays the SE events of 52 genes as candidates at the nexus of hnRNPA1L2-H3.3-G34W-E2F. The list includes genes encoding actin-associated proteins that are essential in shaping and anchoring actin, building intracellular structures centered on CD44, and cell-cell interactions (e.g., ACTN1, SVIL, TPM1 and FLNB). We also found genes encoding proteins linked to cytoskeleton remodeling (e.g., ENAH, ABI2, MBLN1, and TIAL1). The study identified genes related to ubiquitination and protein homeostasis, such as PSMD6 and UBE3A, as well as genes involved in protein transport and components of the nuclear pore complex, such as SEC31A, GOSR2, and SEC13. Additionally, genes associated with chromatin, such as ZZZ3, CTBP1, and EEF1D, and transcription and translation regulators, such as E2F6, PHF6, and MEF2A were among the targets. These genes are diverse examples of genes with SE events where hnRNPA1L2 and H3.3-G34W coincide based on their interaction, strongly suggesting a functional association. SE events align mainly with exons over annotated gene features (Fig. 4d). We also found instances of direct overlap between RIP-seq peaks and SE events (Fig. 4e) that resulted in quantitative changes in splice variants (Fig. 4f).

To gain a better understanding of where they would exert their joint function on the gene body, the RIP-seq peaks were graphed. The results showed that hnRNPA1L2 in parental cells mainly associates with the 3’UTR, while H3.3-G34W were distributed over both 3’UTR and 5’UTR (Figs. 5a and 4e). Notably, hnRNPA1L2 relocated almost half of its peaks to the upper coding region in the H3.3-G34W genetic background, which was further evident when plotting overlapping peaks of H3.3-G34W and hnRNPA1L2. It is noteworthy that our analysis has assigned most hnRNPA1L2 peaks to the exons which include 5’-exons (Fig. 3e). These results suggest a significant overlap between hnRNPA1L2 and H3.3 when altered to accommodate G34W.

Fig. 5: Genic distribution of hnRNPA1L2 uncovers H3.3-G34W mimicry.

a The RIP-seq peak density of hnRNPA1L2 over annotated genic regions accumulates over the 3’UTR in WT (purple line), but partially redistributes over the 5’UTR in H3.3-G34W (green line). b hnRNPA1L2 targets exons for suppression and in the process becomes hijacked by H3.3-G34W to the histone H3.3 normal distribution locations. Importantly, the number of hnRNPA1L2 peaks increases 4-fold in H3.3-G34W. c This illustration shows the increased E2F targeting at the 5’UTR in H3.3-G34W by hnRNPA1L2 (represented by green bars), which contributes to reduced SE events at the 3’UTR and elevated promoter suppression at the 5’UTR through E2F suppression (such as E2F4 and 6).