The translatome of glioblastoma

Abstract

Glioblastoma (GB), the most common and aggressive brain tumor, demonstrates intrinsic resistance to current therapies, resulting in poor clinical outcomes. Cancer progression can be partially attributed to the deregulation of protein translation mechanisms that drive cancer cell growth. In this study, we present the translatome landscape of GB as a valuable data resource. Eight patient-derived GB sphere cultures (GSCs) were analyzed using ribosome profiling and messenger RNA (mRNA) sequencing. We investigated inter-cell-line differences through differential expression analysis at both the translatome and transcriptome levels. Translational changes post-radiotherapy were assessed at 30 and 60 min. The translation of non-coding RNAs (ncRNAs) was validated using in-house and public mass spectrometry (MS) data, whereas RNA expression was confirmed by quantitative PCR (qPCR). Our findings demonstrate that ribosome sequencing provides more detailed information than MS or transcriptional analyses. Transcriptional similarities among GSCs correlate with translational similarities, aligning with previously defined subtypes such as proneural and mesenchymal. Additionally, we identified a broad spectrum of open reading frame types in both coding and non-coding mRNA regions, including long non-coding RNAs (lncRNAs) and pseudogenes undergoing active translation. Translation of ncRNAs into peptides was independently confirmed by in-house data and external MS data. We also observed that translational regulation of histones (downregulated) and splicing factors (upregulated) occurs in response to radiotherapy. These data offer new insights into genome-wide protein synthesis, identifying translationally regulated genes and alternative translation initiation sites in GB under normal and radiotherapeutic conditions, providing a rich resource for GB research. Further functional validation of differentially expressed genes after radiotherapy is needed. Understanding translational control in GB can reveal mechanistic insights and identify currently unknown biomarkers, ultimately enhancing the diagnosis and treatment of this aggressive brain cancer.

Keywords: glioblastoma; non‐coding RNA; radioresistance; radiotherapy; translatome.

Fig. 1

Outlining procedure for ribosome profiling in eight glioblastoma patient‐derived GSCs. (A) Schematic overview experimental approach. GSCs were analyzed with ribosome profiling and mRNA sequencing. With ribosome profiling we were able to capture actively translated genes at a given time point and mRNA sequencing was used to map the identified open reading frames (ORFs). The lower part shows the flowchart of computational analysis of Ribo‐seq data. RibORF [37] was used to filter mRNA sequences being translated by ribosomes. RiboReadCountByFrame is our laboratory‐developed bioinformatic script to identify newly translated ncRNAs and to generate gene count numbers per translated gene and consists of the following steps; first, unnecessary adaptor sequences are removed, followed by sorting based on specific barcode sequences for tracking experimental branches. Unique Molecular Identifiers are then extracted to cluster similar genetic sequences, which are mapped to a standard human genome to identify their location. Redundant sequences are filtered out, and ribosomal RNA is removed. After ensuring sample quality, ORFs are identified, and potential new protein‐building regions are predicted. In preparing the mRNA‐seq reference file, genetic mapping to a standard human genome is conducted to create a comprehensive list of genetic instructions. This includes identifying transcripts, merging data with known genetic sequences, and discovering new types of genetic sequences [129]. Quantifying genetic data involves counting the number of reads associated with each gene using the newly assembled transcriptome for mRNA‐seq and specialized code for ribosome profiling. Data analysis includes examining patterns of genetic reads around start and stop codons, comparing ribosome profiling and RNA sequencing data with protein data, identifying different types of protein‐building regions and start codons, assessing coding potential of newly identified regions, and exploring non‐coding genes that are transcribed. DE, differential expression; GSC, glioma sphere cells; mRNA, messenger RNA; n, represents the number of patient derived GSC samples; nt, nucleotide; ORF, open reading frame; RNAseq, RNA sequencing; rRNA, ribosomal RNA. The mRNA and ribosome profiling data are based on one technical replicate per GSC sample. (B) Read densities around start and stop codons of canonical ORFs. Good quality ribosome footprints show a clear three‐nucleotide periodicity after nuclease digestion and follow a typical high‐low‐low number of reads pattern (depicted in blue‐orange‐green), arising from the translocation of ribosomes along the mRNA one codon at a time. CDS, coding sequence; ORF, open reading frame; P, P‐site (i.e., binds to the tRNA holding the growing polypeptide chain of amino acids); t‐RNA, transfer‐RNA; UTR, untranslated region. (C) Percentage of reads in each nucleotide (nt) of codons, respectively first (P‐site), second, and third nucleotide of all ribosome profiling samples combined, illustrating the clear three‐nucleotide periodicity expected from ribosome footprints after appropriate nuclease digestion. P‐values were determined by Student's t‐test (***P value < 0.001; ****P value < 0.0001). (D) Identified transcriptome (mRNA sequencing) and translatome (ribosome profiling) genes in four GSC cultures under normal circumstances (n = 4, i.e., GSC34, GSC2, GSC20, and GSC28), compared with previously published proteome (mass spectrometry) data of these four identical cell lines [38]. The majority of identified genes show overlap in all three levels of gene expression. Ribosome profiling revealed 10 times as many translated genes compared to mass spectrometry (12 439 genes vs. 1385 genes). Proteome, translatome, and transcriptome data are based on one technical replicate per GSC.

Fig. 2

GSCs contain numerous open reading frame types and start codons. (A) Percentage of identified unique ORF types in all GSC samples (n = 24, i.e., three different conditions per GSC sample of eight GSC samples in total, based on one technical replicate per sample). Conditions per GSC sample are control (t = 0), 30 min after 2 Gy radiotherapy (t = 1) and 60 min after 2 Gy radiotherapy (t = 2), respectively. Of the identified ORFs 17% are previously annotated (main ORF), 44% are variants of CDSs (alt ORFs), and 39% are candidate ORFs. ORF, open reading frame; uORF, upstream ORF. (B) Correlation of ORF type abundance per GSC sample (n = 8, control t = 0) and between GSC samples. The relative distribution of ORF types is consistent between GSCs, with canonical, noncoding, and truncation ORFs being most prevalent. There is no relevant association in ORF type abundance per GB subtype. Subtypes are defined and color‐coded in both TCGA transcriptome analysis [64] and single‐cell RNA‐sequencing analysis [42]. AC, astrocyte‐like; CL, classical; MES, mesenchymal; NPC, neural‐progenitor‐like; OPC, oligodendrocyte‐progenitor‐like; PN, proneural. (C) Start codon types per ORF type, wherein canonical ORFs mainly use AUG whereas other ORF types use a variety of start codons. (D) ORF length (amino acids, aa) per ORF type. Overall, candidate ORFs contain shorter ORF lengths (non‐coding median length 312 aa; uORFs median length 318 aa) compared with canonical (median length 912 aa).

Fig. 3

Ribosome profiling reveals non‐protein coding biotypes being translated into GB. (A) Donut chart with the number of non‐coding genes transcribed (n = 3158) and translated in the control GSC samples (n = 8). About 15% of the transcribed non‐coding genes are being translated (n = 483). ncRNA, non‐coding RNA. (B) Translated non‐protein coding biotypes in GSCs, respectively pseudogenes (n = 198), short ncRNAs (n = 62), lncRNAs (n = 202) and newly identified ncRNAs (n = 21). (C) Non‐protein coding biotypes presented per patient. As illustrated, the non‐protein biotypes contain the same inter‐patient distribution patterns. (D) Multiple sequence alignment of all identified novel ncRNAs on amino acid level shows that ORFs are aligned based on size (i.e. complexity, depicted in shades of blue) rather than on amino acid similarity (depicted in shades of red). Asterisk (*) in the figure represents where the novel ncRNAs are located in the plot. (E) Histogram showing the validation of the translation of ncRNAs by MS analysis based on our own data as well as data from the Human PeptideAtlas database. Around 30% of the riboseq identified ncRNAs could be validated by external data.

Fig. 4

GSCs can be classified into four subgroups based on gene expression on translational and transcriptional levels, matching previously published GB subtypes. (A) Schematic overview of the analysis, eight patients at t = 0 (control) were correlated by unsupervised clustering to identify subgroups on transcriptional and translational level. (B) Dimension reduction analysis using tSNE of the transcriptome (mRNA‐seq, upper panel) and translatome (Ribo‐seq, lower panel) of the GSCs resulted in four subgroups (subgroup 1 GSC2, GSC34; subgroup 2 VU593, VU598; subgroup 3 VU609, VU591; subgroup 4 GSC20, GSC28, respectively). Correlation of the patients on both transcriptional and translational levels resulted in similar subgroups. Subgroups 1 and 4 were classified respectively as the Proneural (PN) and Mesenchymal subtypes, according to previously published TCGA GB subtype data [64] (Fig. S6A). Classification of subgroups based on public data using the TCGA transcriptome analysis [64] or single‐cell RNA‐sequencing analysis [42] is shown for each GSC cell line. AC, astrocyte‐like; CL, classical; MES, mesenchymal; MES, mesenchymal‐like; OPC, oligodendrocyte‐progenitor‐like; PN, proneural. (C) Corresponding heatmap of 61 differentially expressed (DE) genes at t = 0 (control) at translatome level between subgroups. (D) Top 10 GO terms (DAVID functional analysis) matching the DE genes of C, matching subgroup 1 (PN). The orange dots show the number of genes that are common between the GO term's gene set and the respective DE gene set. The gray bars are the −log of the P‐value determined by DAVID functional analysis (P‐value < 0.0001).

Fig. 5

Direct effect of radiotherapy results in alterations of histone and splicing factor dynamics in GB. (A) Schematic overview of the experiment, we used eight GSC cell lines in total and were radiated once with 2 Gy. Ribosome profiling and mRNA sequencing was performed on respectively control (t = 0), 30 min after irradiation (t = 1) and 60 min after irradiation (t = 2). Gy, gray; RT, radiotherapy. (B) MA plot illustrating 43 genes were significantly differentially expressed over time after radiation, (adjusted P‐value < 0.05), respectively, six genes up‐ and 37 downregulated (NS, non‐significant; P, adjusted P‐value). P‐values were determined by Wald's test. (C) Corresponding heatmap of 43 DE genes over time, wherein mainly histone genes are lower expressed after irradiation and splicing factors are higher expressed over time. Gy, gray; NS, non‐significant; RT, radiotherapy; t ₀, control; t ₁, 30 min after radiotherapy; t ₂, 60 min after radiotherapy.