Pan-cancer chromatin analysis of the human vtRNA genes uncovers their association with cancer biology

ATAC-seq data

The ATAC-seq data obtained by The Cancer Genome Atlas (TCGA) consortium were retrieved from UCSC Xena Browser (Goldman et al., 2020) on 03/10/2018. It comprises the genomic matrix TCGA_ATAC_peak_Log2Counts_dedup_promoter with normalized count values for 404 samples (385 primary solid tumor samples across 23 tissues). As is described in UCSC Xena Browser to calculate the average ATAC-seq values a prior count of 5 was added to the raw counts, then put into a "counts per million", then log2 transformed, then quantile normalized; the result is the average value in the file (log2(count+5)-qn values) across all technical replicates and all biospecimens belonging to the same TCGA sample group. The gene promoter for ATAC-seq is defined as a region within -1000 to +100bp from the TSS site. Assignment of promoter peak to gene mapping information derives from the peak summit within the promoter region. Peak location information of the ATAC-seq values for 500 bp gene promoter regions was retrieved from the file TCGA_ATAC_Log2Counts_Matrix.180608.txt.gz. All the analyses were performed for tissues with at least five primary tumor samples (21 tissues), which resulted in the exclusion of CESC (2 primary tumor samples) and SKCM (4 primary tumor samples). Correlations between ATAC-seq and DNA promoter methylation were performed for the 329 primary tumors that are studied by both strategies (Extended Data: Tables S1 and S2).

DNA methylation data

The DNA methylation data obtained from The Cancer Genome Atlas (TCGA) consortium was retrieved from UCSC Xena Browser (Goldman et al., 2020) on date 03/10/2018. It comprises the normalized beta-value of DNA methylation obtained using Illumina Infinium Human Methylation 450 BeadChip arrays of the total (9149 samples) normal adjacent tissues (746 samples across 23 tissues) and primary solid tumors (8403 samples across 32 tissues) of the Pan-Cancer TCGA cohort. The current study is limited to solid tumors since only AML DNA methylation is available at the TCGA. All the analyses were performed for tissues with at least five samples. Normal adjacent tissues comprise 16 tissues, excluding CESC (two samples), GBM (two samples), PCPG (three samples), SARC (four samples), SKCM (two samples), STAD (two samples) and THYM (two samples). Primary tumors comprise 32 tissues. The normalized promoter average beta-values (500 bp bin) comprise the following CpG sites for vtRNA1-1: cg12532653, cg05913451, cg14633504, cg16615348, cg25602765, cg13323902, cg18296956; for vtRNA1-2: cg21161173, cg13303313, cg00500100, cg05174942, cg25984996, cg11807153, cg15697852; for vtRNA1-3: cg23910413, cg19065177, cg02053188, cg01063759, cg07379832, cg07741016; and for vtRNA2-1: cg16615357, cg08745965, cg00124993, cg26896946, cg25340688, cg26328633, cg06536614, cg18678645, cg04481923 (Extended Data: Tables S2–S5).

Analysis of transcription factors associated to vtRNA genes

The transcription factors (TFs) occupancy analysis of vtRNAs was performed at a ± 3000 bp region centered at each vtRNA gene. The CHIP-seq data was retrieved from the UCSC genome browser (GRCh37/hg19) (Kent et al., 2002) and the K562 chronic myelogenous leukemia cell line was chosen since it has the richest TFs data of ENCODE dataset (ENCODE Transcription Factor Binding tracks - ENCODE 3 TFBS) (Euskirchen et al., 2007). The KEGG pathways enrichment analysis of the core 23 TFs common to all vtRNAs was done using STRING software (Szklarczyk et al., 2017) (Extended Data: Table S6).

Overall survival analysis

The Kaplan Meier curve analysis (Bland & Altman, 1998) was performed with software GraphPad Prism 6 using Log-ranked Mantel-Cox test. We analyzed the ATAC-seq normalized values and DNA methylation average normalized beta-values for vtRNA promoters and survival data until 4000 days (99% of data). The log-rank statistics was calculated for quartiles of vtRNAs promoter chromatin accessibility comparison groups (25^th and 75^th) (Extended Data: Tables S7 and S8).

Pathway enrichment and cluster chromosome localization analyses

The genome wide correlation analysis of promoter chromatin accessibility (ATAC-seq normalized values) for each vtRNA in comparison with all the annotated genes was performed using the genomic matrix TCGA_ATAC_peak_Log2Counts_dedup_promoter retrieved from UCSC Xena Browser (385 primary tumor samples across 23 tissues) (Goldman et al., 2020) on date 03/10/2018. The protein coding genes showing a Spearman correlation r ≥ 0.4 were selected for pathway and cluster location analysis. The pathway enrichment analysis was done using STRING software for the Gene Ontology Biological Process and KEGG pathways categories (Szklarczyk et al., 2017). The cluster localization analysis was performed with two approaches, the Cluster Locator algorithm (Pazos Obregón et al., 2018) and the Enrichr software (Kuleshov et al., 2016) (Extended Data: Tables S9 and S10).

Immune subtypes data

The data of Immune Subtypes defined by Thorsson et al. (Thorsson et al., 2018) is available at the Pan-Cancer TCGA and was retrieved from UCSC Xena Browser (Goldman et al., 2020) on date 03/10/2018 and from the supplementary material of Thorsson et al. article (Thorsson et al., 2018). ATAC-seq data of vtRNA promoters for the six Immune Subtype comprised a total of 364 samples: wound healing (105 samples), IFN-g dominant (93 samples), inflammatory (110 samples), lymphocyte depleted (43 samples), immunologically quiet (nine samples), TGF-b dominant (four samples). DNA methylation average normalized beta-values of vtRNA promoters for the six Immune Subtype comprised 7693 samples including wound healing (1963 samples), IFN-g dominant (2137 samples), inflammatory (2126 samples), lymphocyte depleted (933 samples), immunologically quiet (383 samples), TGF-b dominant (151 samples) (Extended Data: Tables S11–S14).

Statistical analysis

Unless specified all the variables are expressed as average value ± standard deviation (SD). Statistical analyses were performed with one-way ANOVA for multiple comparison tests, including Sidak’s Honest Significant Difference test as a post-hoc test, as referred in the legend of the Figures. Spearman equations were applied to determine the correlations. All the analyses were done in R software (version 3.6) using libraries ggcorrplot, corrplot, xlsx, heatmap.2 (clustering distance measured by Euclidean and Ward clustering algorithms) and software GraphPad Prism 6. The statistical significance of the observed differences was described using the p-value (* p < 0.05, ** p < 0.01, *** p < 0.001, **** p-value < 0.0001). Differences with a p-value of < 0.05 were considered significant.

VtRNA genes have different chromatin accessibility at the promoter region

The Cancer Genome Atlas consortium incorporated ATAC-seq data for tumor samples of the GDC Pan-Cancer dataset (385 primary tumors samples across 23 cancer types (Corces et al., 2018)). The ATAC-seq strategy is a genome-wide sequencing profiling approximation to chromatin accessibility using the hyperactive Tn5 transposase that simultaneously cut and ligate adapters preferentially to the accessible DNA regions. The counts of sequencing reads in a particular genomic region provides a direct measurement of chromatin accessibility (Buenrostro et al., 2015; Tsompana & Buck, 2014). For those gene promoters regulated by DNA methylation, ATAC-seq is a more direct approximation to chromatin accessibility than promoter CpG DNA methylation measurements because the latter is well correlated but not an exact measurement of chromatin accessibility. Indeed, the effect of DNA methylation on chromatin structure depends on the gene and the distribution of the methylation sites along it. It is also important to mention that gene expression interpretations based solely on chromatin accessibility are blind to post-transcriptional regulatory events such as processing, localization, and stability of the RNA transcript. Nevertheless, since the chromatin status of the different vtRNA promoters has been positively correlated with the abundance of their transcripts (Ahn et al., 2018; Fort et al., 2018; Helbo et al., 2015; Lee et al., 2014a; Lee et al., 2014b; Park et al., 2017a; Park et al., 2017b; Sallustio et al., 2016; Treppendahl et al., 2012), ATAC-seq is likely a good surrogate of vtRNA expression. DNA CpG methylation data is available for 746 normal adjacent (across 23 tissues) and 8403 primary tumor (across 32 tissues) tissue samples of the Pan-Cancer TCGA cohort (Extended Data: Tables S1–S5 and Figure S1). The full data is listed in Tables S1–S5 and the tissue composition of the three datasets of the cohort, comprising tumor ATAC-seq, normal and tumor DNA methylation, is shown in Figure S1. The DNA methylation data has nine more tissue types and a higher number of samples than the others, while the abundance of some tissues is skewed to twice enrichment in 1-4 specific tissue types.

Although several reports demonstrated that the DNA methylation of vtRNA promoters is well correlated to vtRNA transcript expression in various tissues (Ahn et al., 2018; Fort et al., 2018; Helbo et al., 2015; Helbo et al., 2017; Lee et al., 2014a; Lee et al., 2014b; Park et al., 2017b; Sallustio et al., 2016; Treppendahl et al., 2012), the association between the DNA methylation and chromatin accessibility at their promoters was not previously investigated.

To assess the expression of vtRNAs in cancer tissues, we analyzed a 500 bp region, containing the full vtRNA genes, the proximal promoter region and the RNA polymerase lII A and B box elements located downstream of the transcription start site (TSS) (Figure 1). Indeed, Vilalta et al., 1994 demonstrated that the deletion of the 300bp 5´-flanking region of the rat vtRNA transcript largely reduced its transcription greater than 30-fold (Vilalta et al., 1994). As depicted in Figure 1, this region bears epigenetic marks of transcriptionally active chromatin (Li et al., 2007), including histone H3K27ac, H3K4me1, H3K4me3 at the promoters and DNaseI hypersensitivity clusters around the TSS in the cell types compiled by UCSC browser (Goldman et al., 2013). Interestingly, vtRNA2-1 is the only vtRNA immersed in a CpG island and vtRNA1-1 displays the strongest epigenetic marks of transcriptional activity (Figure 1). The analysis of the ATAC-seq values of this 500 bp region of all the TCGA primary tumors reveals that the four human vtRNAs have different levels of chromatin accessibility (Figure 2A). Moreover, the average DNA methylation of the promoters (Figure 2B) in primary tumor samples supports this finding (32 tissues, Extended Data: Figure S1 and Tables S2 and S4). VtRNA1-1 promoter has the highest chromatin accessibility and the lowest DNA methylation, showing the smallest dispersion of values among the samples (Ave. 4.1, 0.6 SD) relative to the other three vtRNAs. Meanwhile, vtRNA1-2 (Ave. 2.9, 1.6 SD), vtRNA1-3 (Ave. 1.7, 1.2 SD) and vtRNA2-1 (Ave. 2.5, 1.4 SD) have broader ranges of ATAC-seq and DNA methylation values (Figure 2, Table 1, Extended Data: Tables S1 and S2 and Figure S1). A comparison among the vtRNAs shows that the average chromatin accessibility of their promoters is inversely proportional to their average DNA methylation (Table 1 and Figure 2) (relative ATAC-seq values are vtRNA1-1 (4.1) > vtRNA1-2 (2.9) > vtRNA2-1 (2.5) and accordingly, relative DNA methylation are vtRNA1-1 (0.09) < vtRNA1-2 (0.4) < vtRNA2-1 (0.5)). Yet vtRNA1-3 average low chromatin accessibility is not accompanied by a comparatively denser promoter methylation (values 1.7 and 0.22 respectively). Nevertheless, the intra-sample correlation between the ATAC-seq and DNA methylation for the vtRNAs is between -0.24 and -0.78 (see rs-values in Table 1 and Extended Data: Figure S2), suggesting that regardless of the relatively smaller average effect of DNA methylation on vtRNA1-3 promoter accessibility, all the vtRNAs are regulated by DNA methylation to some extent (Table 1). In addition, the CpG DNA methylation negatively correlates with ATAC-seq values in the 500bp of the promoters of the four vtRNAs (Table 1 and Extended Data: Figure S2). The strength of the correlation for each vtRNA is higher for the two vtRNAs that have a broader spectrum of DNA methylation, i.e. vtRNA1-2 and vtRNA2-1 (rs = -0.74 and rs = -0.78 respectively) (Table 1 and Extended Data: Figure S2B and S2D). Meanwhile, vtRNA1-1 and vtRNA1-3 promoters, which have narrower range of DNA methylation (vtRNA1-1 has also a small range of chromatin accessibility) present a smaller association between both values (rs = -0.24 and rs = -0.54 respectively) (Table 1 and Extended Data: Figure S2A and S2C). Technical differences in ATAC-seq and DNA methylation array approaches, as well as naturally occurring non-linear correlation between average DNA methylation and chromatin structure may account for these differences (Corces et al., 2018; Liu et al., 2018). The low ATAC-seq value of the pseudogenes vtRNA3-1P (Ave. -1.2, 0.6 SD) and vtRNA2-2P (Ave. 0.1, 0.6 SD), represents a proof of concept of the analyses (Extended Data: Figure S3 and Tables S1-S5). The same approach was performed for the Y-RNAs genes RNY1 (Ave. 3.7, 0.4 SD), RNY3 (Ave. 2.1, 0.6 SD, RNY4 (Ave. 1.9, 0.6 SD) and RNY5 (Ave. 3.1, 0.4 SD) (Canella et al., 2010), which code for ubiquitously expressed RNA polymerase III transcripts (Extended Data: Figure S3 and Tables S1).

Figure 1. Genomic view of human vtRNA genes epigenetic features.

Genomic view of 1.5 kb region of human vtRNA genes in UCSC Genome browser (GRCh37/hg19) centered at the 500bp bin highlighted in yellow, which was used for ATAC-seq and CpG methylation analyses. VtRNAs of locus 1 cluster (vtRNA1-1, vtRNA1-2, vtRNA1-3) and of locus 2 (vtRNA2-1) are in the sense and antisense orientation respectively. Several Gene annotation and ENCODE Project tracks for seven cell lines (GM12878, H1-hESC, HSMM, HUVEC, K562, NHEK, NHL) are displayed: DNA accessibility (DNaseI hypersensitivity clusters (color intensity is proportional to the maximum signal strength)), DNA methylation (for CpG islands length greater than 200 bp), histone modification (H3K27Ac, H3K4me1, H3K4me3 marks), conservation of the region in 100 Vertebrates (log-odds Phylop scores). The vertical viewing range of the tracks displays the default settings of the browser for each variable.

Figure 2. Chromatin accessibility and DNA methylation of the vtRNA promoters.

Analysis of vtRNA 500bp promoter region. A. Average chromatin accessibility ATAC-seq values for the 385 primary tumors available at the Pan-Cancer TCGA dataset (23 tissues types). B. Average DNA methylation beta-values for the 8403 primary tumors (32 tissues types). Dashed horizontal lines denote unmethylated (bottom gray, average beta-value ≤ 0.2), 50% methylated (middle red, average beta-value = 0.5) and highly methylated (top gray, average beta-value ≥ 0.8) promoters. The box plots show the median and the lower and upper quartile, and the whiskers the 2.5 and 97.5 percentile of the distribution.

Chromatin accessibility, DNA methylation and transcription factor binding at the four vtRNA promoters are correlated

Aiming to investigate a possible transcriptional co-regulation of vtRNA genes, we determined the pairwise correlations in chromatin accessibility and DNA methylation among the four vtRNA promoters. We found that all pairs of vtRNAs, except vtRNA1-2 and vtRNA2-1, show positive correlations in both datasets (Figure 3A and 3B) (rs-values 0.23-0.41 for ATAC-seq and rs-values 0.20-0.50 for DNA methylation). The unsupervised clustering of these epigenetic marks indicates that the three vtRNAs clustered at vtRNA1 locus (vtRNA1-1, vtRNA1-2 and vtRNA1-3), which are the vault particle associated vtRNAs, are more similar than vtRNA2-1 (Figure 3A and 3B).

Figure 3. Comparative chromatin accessibility and transcription factor occupancy at the vtRNA loci.

A–B Matrix of pairwise Spearman correlations (two-way hierarchical clustering distance measured by Euclidean and Ward clustering algorithms) of vtRNA 500 bp promoter region for ATAC-seq data (385 primary tumors samples across 23 tissues) (A) and DNA methylation average beta-values data (8403 primary tumors samples across 32 tissues) (B). C. Venn diagram of transcription factors identified as ChIP-seq Peaks by ENCODE 3 project in the cell line K562 (Venny 2.1; https://bioinfogp.cnb.csic.es/tools/venny/index.html). The region for TFs assignment was defined as ±3000 bp from the vtRNA transcript sequence boundaries. D. Interaction cluster of the core 23 TFs common to all vtRNAs performed with STRING (Szklarczyk et al., 2017). The colored labels of the top 3 enriched KEGG pathway terms (FDR < 0.05) are indicated.

Given that the chromatin accessibility of a DNA region is partially controlled by transcription and chromatin remodeler’s factors (Corces et al., 2018; Klemm et al., 2019), we analyzed transcription factors (TFs) occupancy at a ± 3000 bp region centered at each vtRNA gene. Since there is no such study of the TCGA samples, we investigated the CHIP-seq data of the K562 chronic myelogenous leukemia cell line, which has the richest TF data at ENCODE dataset (Euskirchen et al., 2007; Goldman et al., 2013). We observed that the four vtRNAs share a common core of 23 TFs (26%). The three members of the vtRNA1 cluster share 16 additional TFs (18%), reaching 39 common TFs (44%), while only 2-4 TFs with vtRNA2-1 (Figure 3C and 3D and Extended Data: Table S6). We finally asked whether this common core of TFs is linked to a specific biological process. The STRING enrichment analysis of the core of 23 TFs common to all vtRNAs identified HTLV-I infection (ATF1, ATF3, E2F1, EP300, NFATC3, TBP), Hepatitis B (E2F1, EP300, NFATC3) and Transcriptional misregulation in cancer (AFF1, ATF1, MAX) as the top three enriched KEGG pathway terms (FDR < 0.05) (Figure 3D and Extended Data: Table S6).

Taken together these findings indicate that the transcriptional activity of the vtRNAs is gene specific, being vtRNA1-1 the most accessible and possibly the most expressed. In addition, the data suggest that the regulatory status of the vtRNA promoters could be coordinately controlled, and the three vtRNA1s are more co-regulated among themselves. Yet, the core of TFs shared by the four vtRNAs is associated with viral infection and cancer related terms in the myeloid cell line K562.

The vtRNA promoters present gene and tissue specific patterns of chromatin accessibility in primary tumors

In order to investigate the expression of vtRNAs in different tissues we analyzed the Pan-Cancer TCGA ATAC-seq data discriminating the tissue of origin. As expected, vtRNA1-1 has the highest promoter chromatin accessibility among tissues and the smallest variation (Figure 4A, Extended Data: Tables S1-S5 and Figure S1). Remarkably, low-grade glioma tumor samples (LGG) show a global reduction in chromatin accessibility at the four vtRNA gene promoters (Figure 4A). An opposite pattern is seen in adrenocortical carcinoma tissue (ACC), where the vtRNAs have the highest concerted chromatin accessibility (Figure 4A). Individually, vtRNAs promoter accessibility is maximum and minimum in THCA (4.7) and LGG (3.3) for vtRNA1-1, ACC (4.2) and LGG (0.29) for vtRNA1-2, ACC (3.7) and LGG (0.24) for vtRNA1-3 and KIRC (4.2) and LGG (1.0) for vtRNA2-1 (Figure 4A). VtRNA3-1P reaches a maximum of promoter chromatin accessibility in Testicular Germ Cell Tumors (TGCT) (0.65) with an extremely low average ATAC-seq value and its minimum in PRAD (-1.5) (Extended Data: Figure S3C). VtRNA2-2P reaches a maximum and minimum promoter chromatin accessibility in Lung Adenocarcinoma (LUAD) (1.0) and ACC (-0.34) respectively (Extended Data: Figure S3C). Remarkably, although vtRNA1-1 and vtRNA1-3 have the highest and lowest promoter accessibility in the majority of tissues respectively, vtRNA1-2 and vtRNA2-1 are the most variable (vtRNA1-1 SD = 0.34, vtRNA1-2 SD = 1.4, vtRNA1-3 SD = 0.75 and vtRNA2-1 SD = 0.78) (Figure 4A). In addition, in 12 tissues the chromatin accessibility of vtRNA1-2 is higher than vtRNA2-1 (57%) whereas in nine tissues the chromatin accessibility of vtRNA1-2 is lower than vtRNA2-1 (43%) (Figure 4A).

Figure 4. Chromatin accessibility and DNA methylation of the vtRNA promoters in tumor and normal tissues.

A. ATAC-seq values for vtRNA promoters in 21 primary tumor tissues. B. Average beta-values of promoter DNA methylation for vtRNAs in 32 primary tumors tissues. C. Average beta-values of promoter DNA methylation for vtRNAs in 16 normal adjacent tissues. The charts show the average value and standard deviation of each vtRNAs for the tissues with at least five samples available at Pan-Cancer TCGA dataset.

Using the vtRNAs promoter DNA methylation data, we extended the analysis to 32 primary tumor tissues, incorporating nine more tissues than ATAC-seq samples (DLBC, KICH, OV, PAAD, READ, SARC, THYM, UCS, UVM) (Figure 4B, Extended Data: Figure S1) and 16 normal adjacent tissues (Figure 4C, Extended Data: Figure S1). Again, the DNA methylation profile mirrors the chromatin accessibility in the individual tissue types (Figure 4). Remarkably, the association between the average vtRNA promoter chromatin accessibility and DNA methylation of the vtRNA1 cluster in each tissue type is higher than that observed for the average of all tissues (vtRNA1-1 rs = -0.28, vtRNA1-2 rs = -0.90 vtRNA1-3 rs = -0.74 and vtRNA2-1 rs = -0.58) (Table 1 and Extended Data: Figure S4). The opposite finding for vtRNA2-1 may be explained by the impact of the previously recognized epigenetic variation in the locus (Silver et al., 2015), which may prevail over the tissue specific variation for this gene.

In primary tumor tissues, promoter DNA methylation of vtRNA2-1 is higher than vtRNA1-2 in 17 tissues (53%) and the opposite is observed in the remaining 15 tissues (47%) (Figure 4B). Meanwhile in normal tissues, promoter DNA methylation of vtRNA2-1 is higher than vtRNA1-2 in 4 tissues (25%) and lower in 12 tissues (75%) (Figure 4C), indicating that although vtRNA1-2 is more methylated in normal samples, vtRNA2-1 gains methylation in neoplastic tissue, becoming more methylated than vtRNA1-2 (in 7 tissue types and vtRNA1-2 loose methylation: BLCA, BRCA, CHOL, HNCSC, LUAD, LUSC, UCEC). Indeed, the Fisher's exact test identifies a mirrored profile of chromatin accessibility between normal (4 tissues with high vtRNA1-2 and 12 tissues with high vtRNA2-1 from 16 total tissues) and tumor tissues (17 tissues with high vtRNA1-2 and 15 tissues with high vtRNA2-1 from 32 total tissues) for vtRNA1-2 and vtRNA2-1 (p-value = 0.07). Despite this deregulation in transformed tissues, we wondered if the tissue specific differences in vtRNA promoters were explained by their pre-existing status in the normal tissues. The correlation among average vtRNA promoter DNA methylation in normal and tumor tissues suggest that the variation in chromatin accessibility among tissue types is already established in the normal tissue counterparts (vtRNA1-1 rs = 0.53, vtRNA1-2 rs = 0.84, vtRNA1-3 rs = 0.47 and vtRNA2-1 rs = 0.75, Extended Data: Figure S5).

VtRNA promoter´s DNA methylation is associated with tumor stage and tissue of origin

The assessment of differential vtRNA expression from normal to tumor condition at the TCGA cohort can only be performed using DNA methylation dataset since no ATAC-seq analyses were performed in the normal tissues. The average beta-value of CpG sites in the 500bp vtRNAs promoter of normal and tumor tissue evidenced gene specific patterns of deregulation in neoplastic tissues (Figure 5). As depicted in Figure 4, vtRNA1-1 and vtRNA1-3 promoter DNA methylation and chromatin accessibility is low and high tissue-wide respectively, in both normal and tumor tissues (yet, there are highly methylated tumor outliers as LGG, that lack normal counterpart) (Figure 4 and Figures 5A and 5C). Conversely, vtRNA1-2 and vtRNA2-1 revealed significant differences in average promoter DNA methylation among the tissue categories (Figures 5B and 5D). In agreement with previous results, the methylation of vtRNA1-2 promoter decreases from normal to tumor samples, while vtRNA2-1´s increases (Figures 5B and 5D and Figure 4).

Figure 5. Promoter DNA methylation of vtRNAs in Normal and Tumor tissues of Pan-Cancer TCGA dataset.

Average beta-values of promoter DNA methylation of vtRNA1-1 (A), vtRNA1-2 (B), vtRNA1-3 (C) and vtRNA2-1 (D), assessed in 746 normal and 8403 tumor tissues respectively. The box plots show the median line and lower and upper quartile and the whiskers the 2.5 and 97.5 percentile. Horizontal grey striped and red dotted lines denote unmethylated (average beta-value ≤ 0.2), 50% methylated (average beta-value = 0.5) and highly methylated (average beta-value ≥ 0.8) promoters. One-way ANOVA multiple test analysis with Sidak as posthoc test was performed. **** p-value < 0.0001.

We next asked whether the promoter status of the vtRNAs is deregulated during the normal and tumor transformation of different tissue types. This comparison was restricted to the 16 tissues with normal samples (with data of at least five samples, Figure 6 and Extended Data: Figure S1, Table S5). The results revealed that vtRNA1-1 is globally unmethylated in normal and tumor from different tissue origins, but presents small but statistically significant changes in 4 tissue types (OG profile in BLCA, THCA UCEC and TSG in LUSC), whose biological impact remains to be determined (Figure 6A). Similarly, vtRNA1-3 is globally unmethylated, showing larger variability and significant upregulation in BRCA and PRAD, compatible with a TSG function (TSG trends for LIHC p-value = 0.08 and LUSC p-value = 0.09, Figure 6C). In agreement with the literature, vtRNA2-1 presents a TSG pattern of deregulation in PRAD, LUSC and BRCA (Cao et al., 2013; Fort et al., 2018; Romanelli et al., 2014), and an OG pattern in KIRP that has not been previously described (Figure 6D). Finally, a statistically significant dysregulation of vtRNA1-2 promoter DNA methylation across almost all cancer types (13 out of 16 tissues analyzed) poses it as a candidate OG in cancer (Figure 6B).

Figure 6. VtRNA promoters DNA methylation in Normal vs. Tumor samples.

Average beta-values of promoter DNA methylation for vtRNA1-1 (A) vtRNA1-2 (B), vtRNA1-3 (C) and vtRNA2-1 (D). Acronyms indicate the tissue condition (normal (N) and tumor (T)). Blue and red colors indicate an increase or reduction in promoter methylation in tumor vs their normal tissues counterparts. The box plots show the median and the lower and upper quartile, and the whiskers the 2.5 and 97.5 percentile of the distribution. Horizontal, lines denote the methylation level of the promoters: grey striped bottom and top for unmethylated (average beta-value ≤ 0.2) or highly methylated (average beta-value ≥ 0.8) respectively, and red dotted for 50% methylated (average beta-value = 0.5)) promoters. One-way ANOVA multiple test analysis with Sidak as posthoc test was performed. * p-value < 0.05; ** p-value < 0.01; **** p-value < 0.0001.

Overall, the data shows different DNA methylation changes in the vtRNAs promoters upon malignant transformation in several tissues. In addition, the deregulation of each vtRNA occurs mostly in a gene specific direction, where vtRNA1-2 has an OG like epigenetic de-repression while vtRNA2-1 (and to a lesser extent vtRNA1-1, vtRNA1-3) has a TSG like repressive methylation in tumor tissues.

Low DNA methylation at the promoter of vtRNA1-1, vtRNA1-2 and vtRNA2-1 is associated with low patient overall survival

Seeking for a possible clinical significance of the chromatin accessibility change of the vtRNA promoters, we studied its association with the overall survival of patients (Figure 7 and Extended Data: Tables S7 and S8). The analysis of patients stratified by vtRNA promoter accessibility quartiles did not show differences in overall survival (data available at Extended Data: Table S8). However, a significantly lower patient survival probability when patient tumors have lower promoter DNA methylation is observed for vtRNA1-1 (p-value = 0.003), vtRNA1-2 (p-value = 0.004) and vtRNA2-1 (p-value = 0.002) (patient stratified in quartiles of the promoter DNA methylation cohort values, Figure 7). Therefore, a lower promoter DNA methylation of vtRNA1-1, vtRNA1-2 and vtRNA2-1, a surrogate of their high expression, might be associated with poor patient overall survival cancer-wide.

Figure 7. Association between patient Overall survival and vtRNA promoters DNA methylation status in Pan-Cancer TCGA dataset.

Kaplan-Meier curves of overall patient survival probability over the time (4000 days) discriminating the primary tumors into two cohorts with relatively higher (75^th) and lower (25^th) DNA methylation at each vtRNA promoter. Promoter DNA methylation data of the 8060 primary tumors was used to stratify the patients in two quartiles, based on the expression of vtRNA1-1 (A) vtRNA1-2 (B), vtRNA1-3 (C) and vtRNA2-1 (D). Patient survival probability of the two groups was analyzed using Log-rank (Mantel-Cox) test. The dotted lines represent the 95% confidence interval for each curve.

Genome-wide chromatin accessibility correlations with vtRNA promoters reveal their link to specific cancer related functions

Seeking to get insight into the function of the vtRNAs from their transcription regulatory marks, we searched for the genes co-regulated at chromatin level. We calculated the correlation of the ATAC-seq promoter values of each vtRNA and all the individual genes in the 385 primary tumor cohort and selected those showing a Spearman correlation rs ≥ 0.4. Using STRING software analysis, we then investigated their connection with Biological Process categories and KEGG pathways (Szklarczyk et al., 2017). A few enriched terms for the genes correlated with vtRNA1-2, vtRNA1-3 and vtRNA2-1 is identified (FDR < 1×10^-3) (Extended Data: Tables S9 and S10). The 196 protein coding genes correlated with vtRNA1-2 are enriched in the Biological Process “Epithelial cell differentiation” (FDR < 1×10^-4) (Extended Data: Table S10). Meanwhile, the 358 protein coding genes correlated with vtRNA2-1 are enriched in the Biological Processes “Immune system process”, “Inflammatory response” and “Immune response” (FDR < 1×10^-5) and the KEGG pathways term “Cytokine-cytokine receptor interaction” (FDR < 1×10^-4) (Extended Data: Table S10). The only six protein coding genes correlated with vtRNA1-3 are enriched in the “Thyroid cancer” KEGG pathway term (FDR < 1×10^-4). Remarkably, although STRING analysis of the vtRNA1-1 associated 37 protein coding genes did not find enriched process or pathways, 36 of the 37 correlated genes are located at chromosome 5q region. Using the Cluster Locator algorithm (Pazos Obregón et al., 2018), we found a statistically significant non-random clustering behavior for the chromosomal location of the genes co-regulated with vtRNA1-1 (p-value < 1×10^-10). Indeed, when the analysis was extended to 173 protein coding genes with Spearman correlation rs ≥ 0.3, 139 genes of 173 were situated at chromosome 5q region (p-value < 1×10^-10). These genes are positioned in particular regions at chr5q31, chr5q35, chr5q23, chr5q14, chr5q32, chr5q34 (FDR < 1×10^-4) (Kuleshov et al., 2016) (Extended Data: Table S10). We did not find a similar phenomenon for the other three vtRNAs. In summary, vtRNA1-3, vtRNA1-2 and vtRNA2-1 transcription may be co-regulated with genes belonging to specific biological processes located elsewhere, while vtRNA1-1 transcription may be controlled by a locally specialized chromatin domain at chromosome 5q region lacking apparent functional relatedness.

Immune Subtypes associated with the vtRNAs profile

Taking into account that vtRNAs have been previously related to the immune response (Golec et al., 2019; Li et al., 2015; Nandy et al., 2009), and that we found an association of vtRNA2-1 with immune related terms, we sought to investigate a putative relation between the vtRNAs expression and the six Immune Subtypes defined by Thorsson et al. (Thorsson et al., 2018), compiled at Pan-Cancer TCGA (Extended Data: Tables S11-S14). These six Immune Subtypes are named because of the foremost immune characteristic, comprising wound healing, IFN-g dominant, lymphocyte depleted, inflammatory, immunologically quiet and TGF-b dominant types (Thorsson et al., 2018). As was expected, mirrored patterns were observed for promoter ATAC-seq and DNA methylation data of the six groups (Figure 8). The immunologically quiet subtype shows a concerted shutdown of all vtRNA promoters (Figure 8), which is largely composed by LGG samples showing low vtRNA ATAC-seq and high promoter methylation average values in LGG tumors (Figures 4A and 4B). The immunological quiet subtype exhibits the lowest leukocyte fraction (very low Th17 and low Th2) and the highest macrophage:lymphocyte ratio (high M2 macrophages) (Thorsson et al., 2018). The remaining subtypes present a similar vtRNA profile except for the inversion of vtRNA1-2/vtRNA2-1 accessibility ratio, which is higher than 1 for the wound healing and IFN-g dominant subtypes while is lower than 1 for the lymphocyte depleted and inflammatory subtypes (Figure 8A). As expected, a mirrored pattern was observed for vtRNA1-2 and vtRNA2-1 at DNA methylation data (Figure 8B). Since the wound healing and IFN-g dominant subtype have high proliferation rate opposite to the lymphocyte depleted and inflammatory subtypes (Thorsson et al., 2018), we evaluated the correlation between vtRNA1-2 and vtRNA2-1 chromatin status and proliferation rate or wound healing status defined by Thorsson et al. (Extended Data: Tables S12 and S14) (Thorsson et al., 2018). A negative correlation between vtRNA1-2 promoter DNA methylation and the proliferation rate (rs = -0.44) and wound-healing features of the tumors is revealed (rs = -0.48) (Extended Data: Tables S12 and S14). These findings suggest that vtRNA transcriptional regulation may be associated with tumor immunity.

Figure 8. Chromatin accessibility and DNA methylation of vtRNA promoters in the Six Immune Subtypes analyzed in the PANCAN TCGA dataset.

A. ATAC-seq values expressed as log2 normalized values for vtRNA promoters grouped in six Immune Subtype: wound healing (105 samples), IFN-g dominant (93 samples), inflammatory (110 samples), lymphocyte depleted (43 samples), immunologically quiet (9 samples), TGF-b dominant (4 samples). B. DNA methylation average beta-values for vtRNA promoters grouped in six Immune Subtype: wound healing (1963 samples), IFN-g dominant (2137 samples), inflammatory (2126 samples), lymphocyte depleted (933 samples), immunologically quiet (383 samples), TGF-b dominant (151 samples) (Thorsson et al., 2018). The charts show the average and standard deviation for each vtRNA promoter in each Immune Subtype group.

Underlying data

All raw data used in this study can be downloaded from Xena Browser, a tool to explore functional genomic data sets (Goldman et al., 2020) (https://xenabrowser.net/), and UCSC genome browser (Kent et al., 2002) (https://genome.ucsc.edu/).

Extended data

Zenodo: Pan-Cancer chromatin analysis of the human vtRNA genes - Supplementary Figures, http://doi.org/10.5281/zenodo.4784353 (Fort & Duhagon, 2021a)

This project contains the following extended data:

Zenodo: Pan-Cancer chromatin analysis of the human vtRNA genes - Supplementary Tables, http://doi.org/10.5281/zenodo.4784473 (Fort & Duhagon, 2021b)

This project contains the following extended data:

Data are available under the terms of the Creative Commons Zero "No rights reserved" data waiver (CC0 1.0 Public domain dedication).