The mechanism underlying redundant functions of the YTHDF proteins

Abstract

The YTH N⁶-methyladenosine RNA binding proteins (YTHDFs) mediate the functional effects of N⁶-methyladenosine (m⁶A) on RNA. Recently, a report proposed that all YTHDFs work redundantly to facilitate RNA decay, raising questions about the exact functions of individual YTHDFs, especially YTHDF1 and YTHDF2. We show that YTHDF1 and YTHDF2 differ in their low-complexity domains (LCDs) and exhibit different behaviors in condensate formation and subsequent physiological functions. Biologically, we also find that the global stabilization of RNA after depletion of all YTHDFs is driven by increased P-body formation and is not strictly m⁶A dependent.

Background

RNA modification represents an important layer of post-transcriptional regulation of RNA metabolism. Among more than 170 distinct chemical modifications identified on cellular RNA [1], N⁶-methyladenosine (m⁶A) is the most prevalent internal mRNA modification in mammals. One of the major pathways through which m⁶A exerts its function is the preferential binding of “reader” proteins to methylated transcripts. Proteins containing the YT521-B homology (YTH) domain, including YTHDF1, 2, and 3 and YTHDC1 and 2 in mammals, are direct m⁶A readers possessing a dedicated m⁶A-binding domain [2,3,4,5]. The binding of YTHDF1 to m⁶A-modified mRNAs was shown to induce their translation, which has been linked to various physiologically relevant processes [6,7,8,9]. YTHDF2 accelerates the decay of its mRNA targets via recruiting the CCR4-NOT deadenylation complex [4, 10]. YTHDF3 facilitates translation and accelerates decay of its mRNA targets [11, 12]. While Ythdf2 knockout is embryonically lethal in mice [13], Ythdf1 knockout mice develop normally within the first 3 months but exhibit defects in long-term learning and memory [8]. Knockdown of either YTHDF2 or YTHDF3 delays somatic cell reprogramming while YTHDF1 knockdown does not affect this process [14]. These reports suggest that the physiological roles of YTHDF proteins are different, and their functions could be specific to certain cellular contexts.

We have reported that the triple knockdown of YTHDF proteins leads to the highest mRNA stabilization compared to the single and double knockdown of YTHDFs [11]. The molecular mechanism behind this observation remains elusive. More recently, one report suggested that YTHDF proteins redundantly function in the decay of methylated RNA and can compensate for each other [15] and that this redundancy accounts for the significant RNA stabilization after triple knockdown of YTHDFs. The synergistic effect on mRNA decay has also been observed by others with triple knockdown of YTHDF proteins [16, 17]; however, various previous reports have also indicated a role for YTHDF1 in translation promotion [18,19,20], which made us speculate that a different mechanism may explain the transcriptome stabilization effect observed with YTHDF1-3 triple knockdown.

Results and discussion

We first examined the notion that YTHDF proteins are not involved in translation regulation but rather act redundantly to destabilize RNA. This model was presented alongside the concordant idea that YTHDF proteins share highly similar RNA targets, protein partners, and biological functions [15].

YTHDF1 promotes translation of its target transcripts

We identified two key methodological choices by Zaccara et al. [15] in their analyses which differ from our methodology and likely contribute to divergent conclusions. First, they analyzed the effects of individual YTHDF proteins by grouping RNA by m⁶A modification status, not RNA binding of individual proteins. Specifically, there are 7105 m⁶A-modified genes among ~16,000 expressed genes in HeLa cells. Among 6814 mRNA with translation efficiency data acquired with Ribo-seq, 4424 m⁶A-modified mRNA were used by Zaccara et al. for their YTHDF1 knockdown analysis (Fig. 1a), which far exceeds the ~753 high confidence transcripts directly bound by YTHDF1 identified with photoactivatable ribonucleoside-enhanced crosslinking and immunoprecipitation (PAR-CLIP) [3, 4, 11, 15] (Fig. 1a). While Zaccara et al. examined the effects of YTHDF1 knockdown on translation or stability of ~60% of the transcriptome, YTHDF1 mainly binds only ~10% of the transcriptome. Thus, the effects of YTHDF1 knockdown on other m⁶A-modified mRNA could be indirect. Our previous m⁶A-QTL studies have already shown that various RBPs can promote or suppress translation through m⁶A [21], which could be responsible for functional outcomes of ~3000 m⁶A-modified mRNA not bound by YTHDFs analyzed by Zaccara et al. (Fig. 1a, yellow bar). Analyzing all m⁶A-modified mRNAs clearly does not represent the effects of only YTHDF1 or YTHDF2.

Fig. 1

YTHDF1 and YTHDF2 have different molecular functions. a UpSet plots comparing transcripts analyzed with m⁶A-centric and RBP target-centric approaches. b Venn diagram showing RNA analyzed by Zaccara et al. and RNA analyzed in this work. c Cumulative plots showing changes in RNA translation efficiency (TE, left) and abundance (right) after YTHDF1 knockdown using data from Zaccara et al. P values were determined by a Mann-Whitney-Wilcoxon test. d Cumulative plots showing changes in RNA translation efficiency (TE, left) and abundance (right) after YTHDF2 knockdown using data from Zaccara et al. P values were determined by a Mann-Whitney-Wilcoxon test. e Sequence alignment of YTHDF proteins showing differences in the low complexity domain (LCD) region. Conservation scores were calculated by Jalview [22]. f Heatmaps showing the similarity scores calculated by the BLOSUM62 algorithm of YTHDF LCDs (left) and full-length proteins (right). g Higher-order structures formed by LCDs of YTHDFs under an electron microscopy (EM). Scale bar: 200 nm. h Cytoplasmic RNP granule protein localization map of the cell generated by t-distributed stochastic neighbor embedding (t-SNE) from the CELL MAP project [23]. YTHDF1 and eIFs are delineated in red, YTHDF2 and CNOTs are delineated in blue, and YTHDF3 is delineated in yellow. For boxplots, the center line represents the median, the box limits show the upper and lower quartiles, and whiskers represent 1–99%. P values were determined by a Mann-Whitney-Wilcoxon test

Second, Zaccara et al. analyzed the effects of YTHDFs on mRNA abundance by including only actively translated genes. For their analysis of transcript abundance, the transcripts were pre-filtered to have non-zero read counts in the ribosome-protected fragment (RPF) samples (Fig. 1b). This restrained their analysis to only include actively translated RNA. By studying all detected transcripts (with a sum of > 10 read counts across all samples), we found that the abundance of transcripts with more m⁶A sites tends to decrease more upon YTHDF1 knockdown (Additional file 1: Fig. S1a). In contrast, only YTHDF2 knockdown and triple knockdown cause stabilization of transcripts with more m⁶A sites (Additional file 1: Fig. S1a). Analyzing only the translated genes, we did not find significant correlations between numbers of m⁶A sites and changes in RNA abundance after YTHDF1 or YTHDF3 knockdown (Additional file 1: Fig. S1a). Interrogating functions of YTHDF proteins only based on actively translated RNA may lead to incorporation of indirect effects in the analysis, especially when aiming to elucidate their roles in translation and decay. Moreover, drawing conclusions about the relative roles of YTHDF proteins in translation and decay based on data from cells treated with a translation inhibitor could compromise the analysis.

To clarify the functional effects of YTHDF proteins on their target mRNA, we grouped transcripts according to their binding by YTHDF1 and YTHDF2 in HeLa cells from published PAR-CLIP datasets [3, 4]. Analyzing T-C mutations which are caused by direct protein binding, we showed that YTHDF1 and YTHDF2 have different mRNA targets (Additional file 1: Fig. S1b,c). Applying these groupings to RNA-seq and ribosome profiling data from knockdown experiments of YTHDF1 and YTHDF2 revealed significant differences. Knockdown of YTHDF1 decreases translation efficiency only of YTHDF1 unique targets and YTHDF1/2 shared targets, while mRNA abundance is not significantly altered for any group of genes (Fig. 1c). In contrast, YTHDF2 knockdown leads to more significant stabilization of its RNA targets (DF1/2 shared and DF2 unique) (Fig. 1d). Refining YTHDF targets to m⁶A-modified mRNA leads to the same conclusion (Additional file 1: Fig. S1d). We conclude through these careful analyses that the major effect of YTHDF1 is to promote the translation of its RNA targets, while YTHDF2 plays a greater role in mRNA stability. In contrast, analyzing all m⁶A-modified actively translating genes led to the conclusion that they do not affect mRNA translation [15]. Moreover, a recent report using individual YTHDF proteins fused to RNA editors (TRIBE/STAMP), respectively, to map mRNA binding by YTHDF proteins. This study showed that individual mRNAs could be bound by more than one YTHDF protein over their lifetime and that there are YTHDF1 unique target mRNAs as also shown in our analysis [24]. Evolutionarily, Drosophila only has one YTHDF ortholog, and it promotes translation of its target mRNA transcripts while having little effect on their abundance [25]. These observations reinforce the involvement of YTHDF proteins in mRNA translation regulation, not just decay.

YTHDF1 and YTHDF2 bind different protein partners and form distinct higher-order structures

Zaccara et al. also reported that YTHDF proteins bind similar sets of proteins [15]. Amino acid sequences of proteins determine their higher-order structures and molecular functions. If YTHDF proteins share the same protein partners, RNA targets, and biological functions, as was proposed [15], they should have highly conserved amino acid sequences. However, sequence alignment of human YTHDF proteins shows that they differ in their low-complexity domains, which might lead to distinct features in condensate formation (Fig. 1e). Homology analysis involving the calculation of distance scores by BLOSUM62 suggests that YTHDF2 is the most different, while YTHDF1 and YTHDF3 are more similar to each other (Fig. 1e, bottom right panel). The similarity calculations between YTHDF proteins in their low-complexity domains (LCDs) are within the 49–52% range while those for full-length proteins are 58–69%, indicating that they mainly differ in their LCDs (Fig. 1f). Indeed, we found that the LCDs of YTHDF2 form fibril-like structures distinct from structures formed by YTHDF1 and YTHDF3 under electron microscopy (Fig. 1g and Additional file 1: Fig. S1e). Proteins with low-complexity domains may condense with different protein partners based on their intrinsic aggregate forming properties. RBPs in distinct RNP granules are known to have varied functional outcomes and bind different RNA substrates [26].

The different amino acid sequences and higher-order structures of YTHDF proteins indicate that YTHDF proteins should not share highly similar protein partners. To clarify protein interactions between YTHDFs and either translation or decay machineries, we analyzed a recent protein localization map generated by Bio-ID from HEK293 cells [23]. YTHDF1 tends to display similar subcellular localization in cytosolic RNP granules as eukaryotic initiation factors (eIFs) while YTHDF2 colocalizes better with CNOTs, with YTHDF3 in the margin between CNOTs and eIFs (Fig. 1h and Additional file 1: Fig. S1f). These features are in accordance with their reported functions to promote mRNA translation or facilitate decay. We next examined interacting proteins of YTHDF1, 2, and 3. Using the same proteomics dataset [27], Zaccara et al. reported protein partners by comparing proteins enriched with C-terminal BirA* fusion of YTHDF1 to those with the N-terminal BirA* fusion of YTHDF2 and YTHDF3, respectively [15], while we compared results generated using the C-terminal BirA* fusions of YTHDF1,2,3. CNOT proteins are not found as shared high-confidence protein interactors of C-terminal BirA* fusions of YTHDF1,2,3 (Additional file 1: Fig. S1g). The report constructed the Bio-ID database comparing the overlap between N-terminal fusion and C-terminal fusion of the same protein. The results showed overlap ratio varying between 45 and 92% [27]. Thus, using data obtained from YTHDF proteins with the BirA* fusion at the same terminus provides more appropriate comparison of their protein partners. Therefore, our and others’ results indicate YTHDF1 and YTHDF2 proteins have distinct protein partners and form different higher-order structures.

YTHDF2 knockdown or knockout suffices to cause stabilization of its target transcripts, as reported by Zaccara et al. and others [28,29,30], while YTHDF1 or YTHDF3 single knockdown does not cause significant alteration of RNA abundance. This was attributed to a higher level of YTHDF2 in HeLa cells when compared to YTHDF1 or YTHDF3 using Ribo-seq translation efficiency results [15]. However, we analyzed relative protein levels in HeLa and HEK293T cells directly and found that the YTHDF1 and YTHDF3 levels are higher in HeLa while YTHDF2 is more abundant in HEK293T (Additional file 1: Fig. S1h). If YTHDF1 and YTHDF3 are more abundant than YTHDF2 in HeLa cells, YTHDF2 should not suffice to compensate for YTHDF1 or YTHDF3 knockdown.

YTHDF1-3 triple knockdown leads to increased P-body formation and global mRNA stabilization

If each YTHDF protein has a distinct structure and function, how can we explain the synergistic mRNA stabilization observed upon knockdown of all three YTHDF proteins by us, Zaccara et al. and others [11, 15, 16]? To answer this, we performed mRNA-seq with spike-in calibration. A slight decrease in mRNA abundance was observed after individual knockdown of each YTHDF protein (Fig. 2a, siDF1, siDF2, siDF3); however, triple knockdown of YTHDF1-3 caused global stabilization of the whole transcriptome (Fig. 2a, siDF1-3). The majority of cytosolic mRNA is stabilized by triple knockdown independent of m⁶A methylation or YTHDF binding (Fig. 2a). Differential gene analyses of datasets obtained from Zaccara et al. also show much more significant alteration of the transcriptome (~5-19x more transcripts with adjusted P values < 0.05 for differential expression) with YTHDF1-3 triple knockdown than with knockdown of any single YTHDF protein (Additional file 1: Fig. S2a). These results suggest the perturbation of a more fundamental process regulating mRNA stability in the cytosol after YTHDF triple knockdown.

Fig. 2

Increased P-body formation accounts for the global RNA stabilization after YTHDF1-3 triple knockdown. a Boxplots showing changes in RNA abundance after single knockdown or triple knockdown of YTHDFs in HeLa cells. Gray boxes, non-targets (NT) in each analysis (for siDF1, non-DF1 bound transcripts; for siDF2, non-DF2 bound transcripts; for siDF3, non-DF3 bound transcripts; for siDF1-3, non-methylated transcripts). Colored boxes, YTHDF1 targets (brown), YTHDF2 targets (red), YTHDF3 targets (blue), and m⁶A modified transcripts (purple). b Fluorescence microscopy analysis of P-body numbers after single knockdown or triple knockdown of YTHDFs in HeLa cells. Numbers of DCP1A foci per cell were quantified with CellProfiler 3.0. Number of cells analyzed: siCtrl (32), siMETTL3 (49), siDF1 (57), siDF2 (56), siDF3 (46), siDF1/2/3 (41). c Boxplots showing changes in RNA abundance after triple knockdown of YTHDFs, grouped by both m⁶A modification status and P-body enrichment (depleted: “Depl”, enriched: “Enr”).Red boxes denote m⁶A methylated transcripts and gray boxes denote unmethylated transcripts. d Scatter plots showing enrichment of m⁶A-modified transcripts in the P-body enriched pool. Top: EDC3 dataset, bottom: EDC4 dataset. Red dots denote m⁶A methylated transcripts and gray dots denote unmethylated transcripts. e Boxplots showing changes in RNA abundance after DDX6 knockdown in control (left) or YTHDF1-3 (right) knockdown HeLa cells grouped by both m⁶A modification status and P-body enrichment in WT cells (depleted: “Depl”, enriched: “Enr”). Red boxes denote m⁶A methylated transcripts and gray boxes denote unmethylated transcripts. f Boxplots showing changes in RNA abundance after DDX6 overexpression in control (left) or YTHDF1-3 (right) knockdown HeLa cells grouped by both m⁶A modification status and P-body enrichment in WT cells (depleted: “Depl”, enriched: “Enr”). Red boxes denote m⁶A methylated transcripts and gray boxes denote unmethylated transcripts. For boxplots, the center line represents the median, the box limits show the upper and lower quartiles, and whiskers represent 1.5 × interquartile range. P values were determined by a Mann-Whitney-Wilcoxon test

We found that depletion of all YTHDF proteins caused increased numbers of processing bodies (P-bodies) in cultured HeLa cells by DCP1A staining (Fig. 2b and additional file 1: Fig. S2b). Cytosolic P-bodies are hubs for RNA processing with reports suggesting roles in facilitating decay or RNA stabilization [31,32,33]. We decided to characterize the P-body-associated transcriptome in HeLa cells in order to further assess how P-body perturbation may affect m⁶A and non-m⁶A methylated transcripts. We performed RIP-seq using two individual antibodies against EDC3 and EDC4, respectively. The enrichments of transcripts in P-bodies in the two datasets correlate well (Additional file 1: Fig. S2c). Thus, we used the averaged log₂(enrichment) from these two datasets to define the P-body transcriptome. Unlike P-body depleted (“Depl”) transcripts, m⁶A methylation does not cause a more significant stabilization effect for P-body enriched (“Enr”) transcripts after YTHDF1-3 triple knockdown (Fig. 2c). This indicates that YTHDF1-3 triple knockdown preferentially stabilizes P-body enriched transcripts. Moreover, the m⁶A-modified transcripts in the P-body-enriched fraction are more stabilized after YTHDF1-3 triple knockdown (Fig. 2c). We also found that m⁶A-modified transcripts are significantly enriched in the immunoprecipitated fractions in both datasets (Fig. 2d); we categorized transcripts based on numbers of m⁶A peaks and confirmed that groups with more m⁶A peaks are more enriched in P-bodies (Additional file 1: Fig. S2d).

To study whether P-body dynamics account for the transcriptome stabilization observed after YTHDF1-3 triple knockdown, we knocked down or over-expressed DDX6, a protein shown to be essential for P-body assembly [34]. This led to significant alterations of P-body numbers stained with an anti-DCP1A antibody (Additional file 1: Fig. S2e). RNA sequencing with spike-in normalization showed that abolishment of P-bodies leads to global destabilization of RNA (log₂(FoldChange) < 0) as would be expected (Fig. 2e). Upon YTHDF1-3 triple knockdown, P-body enriched transcripts are more stabilized compared to the P-body depleted transcripts (Fig. 2c). Although P-bodies enrich m⁶A-modified transcripts, our results suggest that the stabilization of transcripts following depletion of all three YTHDF proteins could be a result of increased P-body formation in cells rather than an m⁶A-dependent process, and that the relationship between m⁶A and stabilization after triple knockdown could be confounded by the role and dysregulation of P bodies.

We performed YTHDF protein triple knockdown in DDX6-depleted cells with small interfering RNA (siRNA). The global stabilization effect of YTHDF triple knockdown was almost abolished in DDX6-depleted cells (Fig. 2e). Over-expression of DDX6 caused global stabilization of cellular mRNA, and YTHDF triple knockdown exaggerated this effect (Fig. 2f). Of note, m⁶A modification does not cause more significant stabilization of mRNA in P-body enriched fractions in both DDX6 knockdown and overexpression (Fig. 2e and f). Conversely, P-body enriched mRNA are always more stabilized, suggesting that P-bodes are playing a critical role in global stabilization following YTHDF triple knockdown. Collectively, our results show that the global stabilization of mRNA after depletion of all YTHDF proteins is a result of increased P-body formation and is not strictly m⁶A dependent.

Conclusions

The effects of m⁶A on RNA fate are heterogeneous and depend heavily on biological context. The context-dependent roles of m⁶A reader proteins enable the m⁶A-mediated multifaceted regulation of multiple biological processes. YTHDF1 and YTHDF2 proteins have notable sequence differences in their low-complexity regions and form different LLPS granules. Using YTHDF1 target transcripts instead of all methylated mRNA for analysis, we show that YTHDF1 indeed promotes translation in HeLa cells.

In this study, we confirm that depletion of all three YTHDF proteins exhibits a synergistic effect to stabilize mRNA [11, 15,16,17]. However, we show that this effect is not strictly m⁶A dependent. The triple knockdown of all three YTHDF proteins leads to increased cellular P-body formation and global stabilization of most mRNAs, regardless of their methylation status. Therefore, in addition to the individual functions of these YTHDF proteins that are affected by their relative levels, these abundant proteins also participate in a scaffolding role to maintain cellular RNP granules. Depriving these proteins may expose cellular mRNAs to induce P-body formation for global stabilization of the whole transcriptome.

Materials and methods

siRNA and plasmid transfections

AllStars negative control siRNA (QIAGEN, 1027281) was used as control siRNA in knockdown experiments. Cells were transfected by using Lipofectamine RNAiMAX (Invitrogen 13778075) for siRNAs (human YTHDF1: QIAGEN SI00764715, human YTHDF2: QIAGEN SI00764757, human YTHDF3: QIAGEN SI04133339, human METTL3: QIAGEN SI05020414 and human DDX6: Dharmacon J-006371-05-0002) according to the manufacturer’s protocols.

Amino acid sequence alignment

Canonical amino acid sequences of YTHDF1-3 from Homo sapiens were retrieved from UniProt Knowledgebase [35]. Sequence alignments were performed with Jalview (version 2.11.2.0) [22]. Conservation scores were calculated with default settings and phylogenetic tree was calculated with the built-in function “Tree” with “Neighbor Joining” and the BLOSUM62 method. Numbers on the tree denote distances in the virtual space depicting similarities between proteins.

Negative staining transmission electron microscopy

A protein solution (5 μL) of the prion-like domain of each YTHDF was loaded on an EM grid for 10 seconds and excess solution was removed via blotting with filter paper. The grid was then washed with water and stained with 5 μL uranyl acetate (2%) for 15 s. All negative staining samples were imaged on a JOEL 1400 microscope.

Western blot

Protein samples were prepared from respective zebrafish embryos by lysis in RIPA buffer (ThermoFisher Scientific 89900) containing 1 × Halt™ Protease and Phosphatase Inhibitor Cocktail (ThermoFisher Scientific 78441). Protein concentration was measured by NanoDrop 8000 Spectrophotometer (ThermoFisher Scientific). Lysates of equal total protein concentration were heated at 90°C in 1 × loading buffer (Bio-Rad 1610747) for ten minutes. Denatured protein was loaded into 4–12% NuPAGE Bis-Tris gels (Invitrogen NP0335BOX) and transferred to PVDF membranes (ThermoFisher Scientific 88585). Membranes were blocked in Tris-Buffered Saline, 0.1% Tween® 20 (TBST) with 3% BSA (MilliporeSigma A7030) for 30 min at room temperature, incubated in a diluted primary antibody solution at 4 °C overnight, and then washed and incubated in a dilution of secondary antibody conjugated to HRP for 1 h at room temperature. Protein bands were detected using SuperSignal West Dura Extended Duration Substrate kit (ThermoFisher Scientific 34075) on a FluroChem R (Proteinsimple).

Fluorescence microscopy

For imaging of P-bodies, HeLa cells were fixed with 4% paraformaldehyde and permeabilized with 0.3% Triton X-100. The blocked coverglass (ThermoFisher Scientific 155409PK) was incubated with an anti-DCP1A-AlexaFluor 488 conjugate (Abcam ab208275) at 4 °C overnight. After three washes with DPBS, the nucleus was counterstained with Hoechst 33342 (Abcam ab228551), and the coverglass was kept in DPBS at 4 °C before imaging. Samples were imaged on a Leica SP8 laser scanning confocal microscope at the University of Chicago. P-body numbers in each cell were quantified with Cellprofiler 3.0 [36] with a custom workflow.

RNA-seq library construction and bioinformatic analysis

Library preparation was performed using a SMARTer Stranded Total RNA-Seq Kit v2 (TaKaRa, 634417) following the manufacturer’s protocols. Sequencing was carried out at the University of Chicago Genomics Facility on an Illumina NextSeq machine in single-end mode with 75 base pairs (bp) per read. Raw reads were trimmed with cutadapt (version 1.10) [37] and then aligned to the human genome and transcriptome (hg38) using HISAT (version 2.1.0) [38] with the parameter ‘--rna-strandness R’. Annotation files (RefSeq, 2020-04-01, in gtf format) were downloaded from NCBI.

P-body isolation and analysis

The protocol was adapted from a previous publication [39]. For each immunoprecipitation, HeLa cells from ten 15-cm dishes were collected on ice and flash-frozen with liquid nitrogen before being kept at – 80 °C until use. The pellet was thawed on ice for 5 min, re-suspended in 1 mL SG lysis buffer (50 mM Tris-HCl pH 7.4, 100 mM KCl, 0.5% NP40, cOmplete mini EDTA-free protease inhibitor (MilliporeSigma 11836170001), 1 U/μl of RNasin Plus RNase Inhibitor (Promega N2611)), and passed through a 25-gauge 5/8 needle attached to a 1-ml syringe 7 times. After lysis, the lysates were spun at 1000 × g for 5 min at 4 °C to pellet cell debris. Fifty microliters and 950 μl of the supernatants were transferred to new microcentrifuge tubes for isolating total and P-body RNAs, respectively. For isolating total RNA, TRIzol^TM Reagent (Invitrogen 15596026) was added to the system and RNA was extracted following the manufacturer’s protocol. Following isopropanol precipitation, the RNA pellet was re-suspended in 50 μl RNase-free H₂O.

The following steps were performed to isolate mammalian P-body cores and extract their RNA: (1) the 950 μl supernatant was spun at 18,000 × g for 20 min at 4 °C to pellet P-body cores. (2) The resulting supernatant was discarded, and the pellet was re-suspended in 1 ml SG lysis buffer. (3) Steps 1 and 2 were repeated to enrich for SG cores. (4) The resulting pellet was then re-suspended in 300 μl of SG lysis buffer and spun at 850 × g for 2 min at 4 °C. (5) The supernatant which represents the mammalian SG core enriched fraction was transferred to a new tube. (6) The enriched fraction was pre-cleared twice by adding 60 μL equilibrated DEPC-treated Protein A Dynabeads (ThermoFisher Scientific 10001D) and nutating at 4 °C for 30 min. Dynabeads were removed using a magnet. (7) Twenty micrograms of EDC3 or EDC4 antibody was added to the enriched fraction and nutated at 4 °C overnight to affinity purify P-body cores. (8) The solution was spun at 18,000 × g for 20 min at 4 °C and the supernatant was discarded to remove any unbound antibody. (9) The pellet was then re-suspended in 500 μl SG lysis buffer and 100 μl of equilibrated DEPC-treated Protein A Dynabeads was added. (10) The sample was nutated for 3 h at 4 °C. (11) The Dynabeads were washed three times with wash buffer 1 (20 mM Tris-HCl pH 8.0, 200 mM NaCl, 1 U/μl of RNasin Plus RNase Inhibitor) for 5 min, once with wash buffer 2 (20 mM Tris-HCl pH 8.0, 500 mM NaCl, 1 U/μl of RNasin Plus RNase Inhibitor) for 5 min, and once with wash buffer 3 (SG lysis buffer + 2 M Urea, 1 U/μl of RNasin Plus RNase Inhibitor) for 2 min at 4 °C. (12) The beads were resuspended in 200 μl of 100 μg/ml Protease K solution (1X TE buffer, 2M Urea, 1 U/μl of RNasin Plus RNase Inhibitor) and incubated for 15 min at 37 °C. (13) TRIzol^TM Reagent (Invitrogen) was added to the samples and RNA was extracted following the manufacturer’s protocol. Following isopropanol precipitation, the RNA pellet was re-suspended in 20 μl RNase-free H2O. After processing data similarly to the RNA-seq workflow described above through alignment, reads on each NCBI annotated gene were called using the DESeq2 package in R [40]. The fold changes from the DESeq2 output were used as fold enrichment fold in P-bodies.

RNA seq analysis with spike-in strategy

The same number of cells were counted, and total RNA was isolated with TRIzol^TM Reagent (Invitrogen 15596026), according to the manufacturer’s protocol. An amount of ERCC RNA spike-in control (Invitrogen 4456740) proportional to the total cell number was added to each purified total RNA sample before library preparation. After RNA-seq data alignment, reads on each NCBI annotated gene were converted to attomoles by dividing by the sum of reads aligned to the ERCC spike-in. Average log₂(Fold changes) between attomole amounts were calculated and analyzed.

Analysis of publicly available RNA-seq data

Raw fastq files were downloaded from Gene Expression Omnibus (accession number GSE134380) [15, 41] and quality checked using FastQC v0.11.5 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Adapters were trimmed using Cutadapt [37]. Reads were then mapped to the hg38 human genome [42] using HISAT2 v. 2.1.0 [43] with the option --rna-strandedness R. Output sam files were converted to bam files, sorted, and indexed using samtools v. 1.7 [44]. Read counts mapping to each gene were obtained using htseq-count [45] with reference to the hg38 annotation gtf file with options -s reverse, -t exon, -f bam, -i gene_id, -m interserction-nonempty, and -r pos. Before further analysis, the sum of counts across all samples was computed for each gene in R v. 4.0.3 using the rowSums function, and genes with 10 or fewer mapped reads were removed. These were taken to represent “expressed” genes. Differential expression analysis of expressed genes was performed with DESeq2 [40] in R v. 4.0.3. Data on N⁶-methyladenosine sites and differential expression of “translated” genes in HeLa was obtained directly from processed data files published on Gene Expression Omnibus from the study of interest. Further analyses were performed and plots were created in R v. 4.0.3 using the following packages: biomaRt v. 2.44.4 [46], ggplot2 v. 3.3.5 [47], dplyr v. 1.0.7 [48], ggpubr v. 0.4.0 [49], and forcats v. 0.5.1 [50]. When performing analyses involving methylation status, unnamed genes were removed since the published file describing m⁶A sites contained only gene names and m⁶A site numbers. For boxplots, outliers were not shown but were included in statistical analyses.