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 m6A modification status, not RNA binding of individual proteins. Specifically, there are 7105 m6A-modified genes among ~16,000 expressed genes in HeLa cells. Among 6814 mRNA with translation efficiency data acquired with Ribo-seq, 4424 m6A-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 m6A-modified mRNA could be indirect. Our previous m6A-QTL studies have already shown that various RBPs can promote or suppress translation through m6A [21], which could be responsible for functional outcomes of ~3000 m6A-modified mRNA not bound by YTHDFs analyzed by Zaccara et al. (Fig. 1a, yellow bar). Analyzing all m6A-modified mRNAs clearly does not represent the effects of only YTHDF1 or YTHDF2.
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 m6A 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 m6A sites (Additional file 1: Fig. S1a). Analyzing only the translated genes, we did not find significant correlations between numbers of m6A 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 m6A-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 m6A-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 m6A 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.
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 m6A and non-m6A 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 log2(enrichment) from these two datasets to define the P-body transcriptome. Unlike P-body depleted (“Depl”) transcripts, m6A 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 m6A-modified transcripts in the P-body-enriched fraction are more stabilized after YTHDF1-3 triple knockdown (Fig. 2c). We also found that m6A-modified transcripts are significantly enriched in the immunoprecipitated fractions in both datasets (Fig. 2d); we categorized transcripts based on numbers of m6A peaks and confirmed that groups with more m6A 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 (log2(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 m6A-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 m6A-dependent process, and that the relationship between m6A 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, m6A 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 m6A dependent.