DHX36 binding induces RNA structurome remodeling and regulates RNA abundance via m6A reader YTHDF1

Abstract

RNA structure constitutes a new layer of gene regulatory mechanisms. RNA binding proteins can modulate RNA secondary structures, thus participating in post-transcriptional regulation. The DEAH-box helicase 36 (DHX36) is known to bind and unwind RNA G-quadruplex (rG4) structure but the transcriptome-wide RNA structure remodeling induced by DHX36 binding and the impact on RNA fate remain poorly understood. Here, we investigate the RNA structurome alteration induced by DHX36 depletion. Our findings reveal that DHX36 binding induces structural remodeling not only at the localized binding sites but also on the entire mRNA transcript most pronounced in 3'UTR regions. DHX36 binding increases structural accessibility at 3'UTRs which is correlated with decreased post-transcriptional mRNA abundance. Further analyses and experiments uncover that DHX36 binding sites are enriched for N6-methyladenosine (m⁶A) modification and YTHDF1 binding; and DHX36 induced structural changes may facilitate YTHDF1 binding to m⁶A sites leading to RNA degradation. Altogether, our findings uncover the structural remodeling effect of DHX36 binding and its impact on RNA abundance through regulating m⁶A dependent YTHDF1 binding.

Fig. 1. In vivo RNA structurome profiling unveils DHX36 depletion-induced global increase of RNA structures.

A Left: schematic illustration of the CRISPR-Cas9 mediated generation of DHX36 knockout (KO) in HEK293T cells. 169 bp of exon 1 of the DHX36 gene was deleted by two single-guide RNAs and CRISPR-Cas9. Right: western blot confirmed the inactivation of DHX36 in KO cells with GAPDH as the loading control. Source data are provided as a Source Data file. B Schematic illustration of in vivo structure mapping in WT and DHX36-KO cells. Orange circle, SHAPE modification; NGS next-generation sequencing, RT reverse transcription. C Sequencing gel for 5.8S rRNA showing the RT stops induced by NAI modification. Source data are provided as a Source Data file. D The average SHAPE reactivity score and E Gini index of reactivity scores of all mRNAs (n = 7792) in WT and DHX36-KO. The most abundant mRNA per gene was analyzed. Wilcoxon signed-rank test was used to calculate the statistical significance. The average fold change (FC) between the KO vs WT is shown. The boxes indicate median (center), Q25, and Q75 (bounds of box), the smallest value within 1.5 times interquatile range below Q25 and the largest value within 1.5 times interquatile range above Q75 (whiskers). F Top: the binned average Gini index across the length of 5’UTRs (5 bins), CDSs (10 bins), and 3’UTRs (5 bins) of all mRNAs. Bottom: The binned ΔGini (KO-WT) across the above regions. G The binned average reactivity and ΔReactivity (KO-WT) across the 5’UTR (25 bins), CDS (50 bins) and 3’UTR (25 bins). The shaded area represents 95% confidence intervals (CIs) of the average ΔGini or ΔReactivity of each bin calculated by paired two-sided Student’s t-test.

Fig. 2. DHX36 binding induces mRNA structural loss and gained accessibility.

A The workflow for analyzing PAR-CLIP and Structure-seq data to identify DHX36 binding-induced structure remodeling. Created in BioRender. Zhang, Y. (2022) BioRender.com/p72o682. B The distribution of DHX36 binding sites in 5’UTRs, CDSs, and 3’UTRs. C The transcriptomic distribution of DHX36 binding sites. D rG4 formation was predicted in DHX36 binding sites, and the distribution of each subtype is shown. E The average SHAPE reactivity and Gini index of the reactivity scores of DHX36-bound mRNAs (n = 1787). The most abundant mRNA per gene was analyzed. The significance was calculated by a two-sided Wilcoxon signed-rank test. The fold change (FC) between the KO vs WT is shown. F Comparison of ΔReactivity and ΔGini (KO-WT) between DHX36-bound and unbound mRNAs. The significance was calculated by a two-sided Wilcoxon rank-sum test. G The binned average reactivity (top) and ΔReactivity (bottom) across the length of 5’UTRs (25 bins), CDSs (50 bins), and 3’UTRs (25 bins) of DHX36-bound mRNAs. The shaded area was plotted in the same way as Fig. 1G. H The average reactivity of 5’UTR, CDS, and 3’UTR of DHX36-bound mRNAs (n = 1787). The significance was calculated by a two-sided Wilcoxon signed-rank test. The boxes indicate median (center), Q25 and Q75 (bounds of box), the smallest value within 1.5 times interquatile range below Q25 and largest value within 1.5 times interquatile range above Q75 (whiskers). I Scatterplot showing regional ΔReactivity vs. GC content of the designated regions. r and P represent the correlation coefficient and statistical significance calculated by the Pearson correlation test. J The proportion of DHX36 depletion-induced DRRs located within (+DHX36) and outside (−DHX36) DHX36 binding sites. K The distribution of DHX36-induced DRRs across different regions of DHX36-bound mRNAs. L The distribution of DRRs across the mRNAs with DHX36 binding to 5’UTR (top), CDS (middle) or 3’UTR (bottom) only. The line plot shows the location of DHX36 binding sites, and the heatmap shows the DRR enrichment. M ΔReactivity (KO-WT) of PYCR1 and CCNA2 mRNA. DRRs and DHX36 binding sites are highlighted in the shaded areas and green arrows, respectively. ΔR and P represent the average difference in reactivity and the significance calculated by the two-sided Wilcoxon signed-rank test.

Fig. 3. DHX36 induces localized structural loss of 3’UTR binding sites.

A Average SHAPE reactivity (top), BPP (middle), and nucleotide composition (bottom) across DHX36 binding sites. Position 0 represents the exact crosslink site identified by PAR-CLIP data. B The percentage of DHX36 binding sites with or without DRRs. C rG4 formation was predicted in DRR-containing DHX36 binding sites, and the distribution of each subtype is shown. D Comparison of average SHAPE reactivity across DHX36 binding sites with and without rG4s in WT and DHX36-KO. E Comparison of average SHAPE reactivity across DHX36 binding sites located in 5’UTR (top), CDS (middle), and 3’UTR (bottom) in WT and DHX36-KO. F The percentage of DHX36-bound DRRs in 5’UTR, CDS, and 3’UTR. G Barplot showing the SHAPE reactivity of the DHX36 binding site within mRNA JUN in WT and DHX36-KO. H Illustration of the folded structure of the DHX36 binding site within mRNA JUN in WT and DHX36-KO. Nucleotides were color-coded based on the SHAPE reactivity scores. The black arrow indicates the exact crosslink site. I Heatmap depicting the coefficients of the explanatory variables in the linear regression models, including localized GC content of DHX36 binding sites, regional GC content, and ln(length) of DHX36-bound mRNA regions. Rows represent the models designed for four kinds of response variables: the localized ΔReactivity (KO-WT) of all, 5’UTR, CDS, and 3’UTR binding sites.

Fig. 4. DHX36 binding-induced 3’UTR structural change is associated with post-transcriptional regulation of mRNA abundance.

A Schematic illustration showing the method for defining DHX36-induced post-transcriptional regulatory effect. The post-transcriptional fold change (pFC) of mRNA abundance was calculated by dividing the fold change in chromatin fraction by that in the whole cell. Created in BioRender. Zhang, Y. (2021) BioRender.com/a92l088. B Scatterplot showing the definition of pUGs and pDGs based on the pFC of mRNA abundance. C Venn diagram showing the overlapping between DHX36-bound mRNAs and pUGs or pDGs. Statistical significance was calculated by the Chi-square test. D Comparison of pFCs of DHX36-bound mRNAs and all protein-coding genes. E Comparison of ΔReactivity of pDG (n = 149) and pUG targets (n = 517). The boxes indicate median (center), Q25 and Q75 (bounds of box), the smallest value within 1.5 times interquatile range below Q25 and largest value within 1.5 times interquatile range above Q75 (whiskers). A two-sided Wilcoxon rank-sum test was used to calculate the statistical significance in (D, E). F–I Scatterplots showing the correlations of mRNA abundance pFCs with reactivity changes of (F) entire mRNAs, G 5’UTRs, H CDSs, and I 3’UTRs of DHX36-bound mRNAs. r and P represent the correlation coefficient and statistical significance calculated by the Pearson correlation test. J Heatmaps showing the hierarchical clustering on the DHX36-bound mRNAs based on their abundance changes and the average reactivity changes of 5’UTR (left), CDS (mid), and 3’UTR (right). For each cluster, the Pearson correlation between the structure changes and the abundance changes was displayed in a scatterplot. The error bands in (F–J) represent the 95% CIs. K Average SHAPE reactivity across DHX36 binding sites within the designated regions of pUGs and pDGs. The colored horizontal lines represent the means of average reactivity across the binding sites. Position 0 represents the exact crosslink sites identified by PAR-CLIP data. L ΔReactivity of DHX36 binding sites within the designated regions of pUGs (n = 517) and pDGs (n = 149). Error bars represent 95% CI of the mean of ΔReactivity. P values were calculated by unpaired two-sided Student’s t-test.

Fig. 5. DHX36-induced RNA structural accessibility facilitates m⁶A-dependent YTDHF1 binding and RNA degradation.

A Schematic illustration showing the investigation of DHX36/m⁶A/YTHDF1 co-regulation of mRNA abundance via RNA structure remodeling. Created in BioRender. Zhang, Y. (2022) BioRender.com/r15e384. B Significant enrichment of m⁶A and YTHDF1 motifs in the 3’UTR binding sites of DHX36. C, D Enrichment of m⁶A and YTHDF1 binding within the true set (n = 1) of DHX36 binding sites compared to 100 sets (n = 100) of random sites. The randomly selected sites in each set are in equal numbers, equal length, and from the same 3’UTR with the DHX36 binding sites. Data are presented as mean values ± SD. Statistical significance was calculated by the Monte Carlo method. E Comparison of pFCs between DHX36-bound mRNAs with m⁶A or m⁶A/YTHDF1 within 3’UTRs and those without m⁶A/YTHDF1. Statistical significance was calculated by the Wilcoxon rank-sum test. F Identification of DHX36/m⁶A/YTHDF1 degradation targets by comparing DHX36 pUG targets with YTHDF1 degradation targets. G Venn diagram showing significant overlapping between 3’UTR YTHDF1 binding sites and DRR-containing DHX36 binding sites within 3’UTRs. Statistical significance in (F, G) was calculated by hypergeometric test. H Comparison of average reactivity and BPP of the regions surrounding m⁶A sites within DHX36-bound mRNAs in WT and DHX36-KO cells. Position 0 and green area denote m⁶A residues and m⁶A motifs, respectively. I Histogram showing the top five m⁶A motifs of YTHDF1 binding sites in WT and DHX36-KO HE293T cells. J Comparison of YTHDF1 binding peak numbers in WT and DHX36-KO cells. K The difference in the normalized ΔSignal (KO-WT) of YTHDF1 peaks within DHX36 binding sites.

Fig. 6. DHX36 facilitates YTHDF1 binding to decrease mRNA stability on selected targets.

A Genome tracks showing RNA abundance in chromatin fraction and whole cell, DHX36 binding sites, m⁶A sites, YTHDF1 binding sites, and SHAPE reactivity of the DHX36 binding sites in mRNAs EZR, MEPCE, PHF23, and ZNF768. B The folded structures of the DHX36 binding sites within the above mRNAs in WT and DHX36-KO cells. m⁶A and YTHDF1 binding sites are highlighted using yellow lines. Nucleotides are color-coded based on the SHAPE reactivity scores. C The abundances of the above mRNAs in WT and DHX36-KO cells were quantified by RT-qPCR (n = 3 biological replicates). D mRNA stability of the above mRNAs was determined by quantifying the mRNA abundance 0, 4, and 8 h after actinomycin D treatment (n = 3 biological replicates). E Schematic illustration showing the construction of the EGFP reporters. F The stability of the EGFP mRNAs fused with DHX36 3’UTR binding site of each selected mRNA was determined by quantifying the EGFP mRNA abundance 0, 2, and 6 h after actinomycin D treatment (n = 3 biological replicates). G The half-life of the above mRNAs in WT and YTHDF1-KO. H RIP assay was performed in WT and DHX36-KO cells with a YTHDF1 antibody. The enrichment of the immunoprecipitated YTHDF1 and GAPDH proteins was assessed by Western blotting. IgG was used as a negative control. I The enrichment of the representative mRNAs in the above RIP was detected by RT-qPCR. The relative enrichment was normalized by input (n = 3 biological replicates). J Western blot confirmed the unaltered protein levels of YTHDF1 in WT and DHX36-KO cells. Data are presented as mean values ± SD in (C, D, F, G, I). The statistical significances in (C, D, F, G, I) were calculated by a two-sided Student’s t-test. Source data are provided as a Source Data file.

Fig. 7. Schematic illustration of the functional mechanism of DHX36-induced structure change in regulating the activity of YTHDF1 binding and target mRNA degradation.

DHX36 binds to and resolves the structured sites. The increased structural accessibility of DHX36-bound m⁶A sites promotes m⁶A- dependent YTHDF1 binding to promote target mRNA degradation. Created in BioRender. Zhang, Y. (2023) BioRender.com/u61b678.