Human cancer cells compensate the genes unfavorable for translation by N⁶-methyladenosine modification and enhance their translation efficiency

Introduction

RNA modification is omnipresent in all living organisms and cells. More than one hundred types of RNA modifications are known (1), of which the two most prevalent in the animal kingdom are A-to-I RNA editing and N⁶-methyladenosine (m⁶A). A-to-I RNA editing is most commonly observed in coding mRNAs from more primitive organisms (2-5) or repetitive sequences in mammals (6-8), whereas m⁶A is widespread in the transcriptomes of many animal species (9-17). In humans, tens of thousands of m⁶A sites have been identified in cell lines, such as HeLa cells (13,14); typically, these sites are enriched around stop codons and have the consensus sequence context GRAC (R = A or G; A = methylated A) (14). Proteins interacting with m⁶A include the m⁶A writers (METTL3 and METTL14), readers (YTH domain proteins) and erasers (18). Functional studies have shown that m⁶A can affect mRNA stability (10,13) or promote translation of host genes (10,14) according to the different readers that bind to m⁶A sites. Interestingly, in human HeLa cells, m⁶A sites bound by the reader protein YTHDF1 increase the translation efficiency of host genes (14), whereas m⁶A sites bound by another reader (IGF2BP) can either facilitate translation or stabilize host mRNAs (10). A handful of studies have discussed the potential relationship between the m⁶A modification and cancer (10,19-21), but no study has systematically investigated how m⁶A can affect the translation of global oncogenes (rather than particular oncogenes). Furthermore, the overall relationship (or overlap) between m⁶A genes and oncogenes is unreported. We believe that this basic information is important for determining whether m⁶A genes and oncogenes are mutually favored or avoided. Importantly, even researchers who mention the enhanced translation efficiency of oncogenes caused by m⁶A have not given a “biological” reason why translation of the target genes should be elevated. Specifically, why should m⁶A facilitate the translation of oncogenes if the oncogenes are already optimized for a high translation efficiency? With the development of next-generation sequencing (NGS) techniques, NGS big data and many bioinformatics tools have been commonly applied to cancer studies (22,23). In this study, by utilizing the YTHDF1 target genes and the normal or si-METTL3 NGS data from HeLa cells generated by a previous work (14), we demonstrate that the m⁶A genes in HeLa cells originally are unfavorable for translation due to their significantly longer mRNA lengths and stronger codon usage bias. However, the translation of these target genes is compensated by the m⁶A modification, because the translation efficiency of the m⁶A genes is significantly higher in normal cells (with m⁶A) than in si-METTL3 cells (no m⁶A). Furthermore, we found enrichment of m⁶A genes in human oncogenes, and these oncogenes were even less suitable for translation. Our results suggest that the compensation of translation by m⁶A may originally have been designed for those oncogenes to help cancer cell growth. Alternatively, the m⁶A modification at least facilitates the translation of target genes that originally are unfavorable for translation. This strategy may be used by cancer cells to increase the protein quantity of a particular set of genes and eventually achieve rapid cell growth. This study deepened our understanding of the role played by m⁶A modification in cancer cells and revealed why m⁶A genes and some oncogenes need methylation to enhance their translation. We also provided novel insights into a potential method to suppress oncogenes in cancer cells.

Methods

Data collection

We collected the m⁶A peaks bound by the reader protein YTHDF1 reported in a previous study (14). The list of oncogenes was downloaded from the latest version of the Cancer Gene Census website (CGC, https://cancer.sanger.ac.uk/census/). NGS mRNA-Seq and Ribo-Seq data from normal or si-METTL3 HeLa cells were obtained from the same study (14). Adenosines within the GRAC (R = A or G; A = methylated A) motif in m⁶A peaks were defined as m⁶A sites in HeLa cells.

Annotation of m⁶A sites

We annotated the m⁶A sites using the hg19 human genome downloaded from the UCSC Genome Browser (genome.ucsc.edu). If an m⁶A site hit multiple isoforms of the same gene, then the transcript with the longest CDS (canonical transcript) was retained. The canonical transcript of each gene was defined by the SnpEff software (24). An m⁶A modification that did not hit any genes was annotated as intergenic.

Processing of the next-generation sequencing data

We aligned the NGS reads (mRNA-Seq and Ribo-Seq from the normal and si-METTL3 cells) to the hg19 reference genome using STAR (25). The uniquely mapped reads were kept for downstream analysis. The read counts of each gene in each sample were calculated by htseq-count (26). In the gene expression analysis, the canonical transcript of each gene was chosen, and all reads that overlapped with exon regions were counted.

Calculating differences in translation efficiency

We counted the reads within CDS regions of the canonical transcripts of each gene. The translation efficiency (TE) was defined as the ratio of normalized Ribo-Seq and mRNA-Seq read counts. Here, using the mRNA-Seq and Ribo-Seq counts from the normal (Control) and si-METTL3 (Treated) conditions, we employed xtail (27) to detect differences in the translation efficiency between the Control and Treated samples. The TE of the Control and Treated conditions and the log₂TE fold change (FC) were given by the software. Genes with a log₂TE FC <0 represent genes with a downregulated TE following si-METTL3 treatment.

Gene ontology (GO) enrichment

The GO analysis was performed using DAVID (28). Highly expressed genes (raw read count >100) in normal HeLa cells were used as background genes.

Conservation analysis

The conservation level of genomic positions is measured by the phyloP score (downloaded from the UCSC Genome Browser, genome.ucsc.edu). Briefly, sites with higher conservation levels have higher phyloP scores. For the comparison of the conservation levels of m⁶A+ and m⁶A− sites in coding regions, sites in different codon positions were compared separately.

Codon usage bias

The protocol used to calculate codon usage bias was described in an earlier study (29). The codon bias of a gene is calculated by the deviation (chi-square) of the A/T content of synonymous codons from that of the intronic regions. Higher deviation indicates a stronger codon bias for a gene. The codon bias of a particular codon is the correlation coefficient between the codon frequency within a synonymous codon family and the deviation (Chi-square value) of each gene (29). A higher correlation coefficient suggests stronger bias for a codon.

Statistical analysis

All statistical analyses were conducted in the R environment (http://www.R-project.org/).

Results

m⁶A methylome in human HeLa cells

We retrieved the m⁶A peaks (YTHDF1 target) identified in HeLa cells from a previous study (14) and extracted all adenosine sites within the GRAC (R = A or G; A = methylated A) motif according to the instructions in the literature. Adenosines located in the GRAC context in m⁶A peaks were regarded as m⁶A sites (Figure 1A). In total, we obtained 18,276 unique m⁶A sites. We annotated these m⁶A sites according to the reference genome hg19, and the canonical transcript of each gene was chosen if a site was located in multiple isoforms (Methods). The majority of m⁶A sites were located in CDS regions, and the 3'UTRs (untranslated regions) also contained a large fraction of m⁶A sites (Figure 1B). Apart from a few m⁶A sites in intergenic regions, most of the sites were assigned to 6,025 unique human genes. The m⁶A genes were significantly enriched in transcription factors based on the GO enrichment analysis (Figure 1C, only GO terms with FDR values <0.05 were listed), which agreed well with known concepts.

Figure 1 Landscape of m⁶A sites (YTHDF1 targets) reported in human HeLa cells. (A) The 4-mer motif around m⁶A sites in HeLa cells; (B) the genomic annotations of m⁶A sites in HeLa cells; (C) gene ontology (GO) of highly expressed m⁶A genes in HeLa cells; (D) conservation level (phyloP score) of m⁶A+ and m⁶A− sites in the coding regions of m⁶A genes. “pos” represents the position of the m⁶A sites on a codon. The Wilcoxon rank sum test was used to calculate the P values; (E) enrichment of codons and amino acids at m⁶A+ sites. “Random shuffle” represents the enrichment observed by shuffling the m⁶A+ sites among all GRAC motifs.

The m⁶A modification was previously reported to increase the translation efficiency of host genes (14,17). Since m⁶A events (most of which are located around stop codons) have been reported to help recruit translation initiation factors and facilitate translation, the exact position of an m⁶A site on a mRNA may not be important as long as the methylation event takes place on this mRNA (in the CDS, UTRs or around stop codons). To investigate whether the particular m⁶A position was important, we sought potential differences between m⁶A sites (m⁶A+) and comparable non-m⁶A sites (m⁶A−). The m⁶A− sites were defined as unmethylated adenosines within GRAC motifs in m⁶A genes (in HeLa cells). We found that the m⁶A+ sites in coding regions largely subjected to natural selection were not more conserved than the m⁶A− sites at the genome level (Figure 1D and Methods). This result agreed with an earlier study, which reported that generally human m⁶A sites were non-conserved (30). The particular m⁶A positions may not be important; otherwise, they would be preserved by natural selection and exhibit high conservation levels. Furthermore, we compared codons containing m⁶A+ and m⁶A− sites in CDSs. Due to the constraint of the GRAC sequence context, the m⁶A sites were only found in a small set of codons, and the m⁶A+ sites did not show any striking enrichment compared to the m⁶A− sites (Figure 1E). Again, we could not deduce any putative function of the m⁶A positions from this result, suggesting that the m⁶A modifications might not exert their function at the “intra-gene” level. Instead, as many previous studies have revealed (14,17), the major function of m⁶A is to increase the translation efficiency of host genes.

Crosstalk between m⁶A genes and oncogenes

m⁶A methylation events are known to promote host gene translation (14,17), and this mechanism does not rely much on the exact position of the methylation sites. Next, we investigated whether human cancer cells could enhance the translation of oncogenes via m⁶A modification.

To address this question, first we searched for known human oncogenes from the Cancer Gene Census (CGC, https://cancer.sanger.ac.uk/census/). We downloaded 719 human oncogenes from the latest version of the CGC website. A total of 335 of these 719 oncogenes were methylated in HeLa cells (Figure 2A). Furthermore, considering the mRNA expression levels, we extracted genes with a raw read count >100 in HeLa cells (see Methods for detail). A total of 4,968 of the 6,025 m⁶A genes and 478 of the 719 oncogenes were highly expressed in HeLa cells. The 4,968 m⁶A genes and 478 oncogenes included 283 overlapping genes (Figure 2B). Thus, more than half of the highly expressed oncogenes were methylated. Then, we examined the gene ontology of the highly expressed oncogenes. We found that these oncogenes were enriched in the transcriptional regulation and metabolism categories (Figure 2C), which agreed with the known features of m⁶A genes. Next, we calculated the number and density of m⁶A sites in oncogenes and other genes. To exclude the potential bias caused by mRNA length (as longer genes tend to bear more m⁶A sites by chance), we ranked all genes into five groups with decreasing mRNA length. Within each bin, we compared the number of m⁶A sites (Figure 2D, top) and the density of m⁶A sites (Figure 2D, bottom) in oncogenes versus other genes. The m⁶A density is defined as m⁶A sites per adenosine, which canceled the bias introduced by gene length. Strikingly, our results show that both the m⁶A number and density is generally higher in oncogenes than other genes (Figure 2D). This pattern indicates a potential functional role of the m⁶A modification in oncogenes.

Figure 2 Crosstalk between m⁶A genes and human oncogenes. (A) Overlap between all m⁶A genes and human oncogenes; (B) overlap between highly expressed m⁶A genes and oncogenes in HeLa cells; (C) gene ontology (GO) of highly expressed oncogenes in HeLa cells; (D) the numbers (top) and densities (bottom) of m⁶A sites in oncogenes and other genes. To exclude the potential bias caused by mRNA length, we ranked all genes into five groups with decreasing mRNA length. Within each bin, we compared the number of m⁶A sites (top) and the density of m⁶A sites (bottom) in oncogenes versus other genes. Wilcoxon rank sum test was used to calculate the P values; (E) fractions of m⁶A sites located in the CDSs and 3'UTRs of the oncogenes or non-oncogenes. Fisher’s exact test was used to calculate the P values. In this Figure, m⁶A+ represents m⁶A genes, and m⁶A− represents non-m⁶A genes.

We began to search for differences between the oncogenes and non-oncogenes (termed other genes) or between the m⁶A genes and non-m⁶A genes. First, we compared the m⁶A site distribution of the oncogenes and the remaining genes among the methylated gene set. We calculated the proportion of m⁶A sites that were located in CDSs and 3'UTRs. Intriguingly, the oncogenes showed a remarkably higher fraction of m⁶A sites in the CDSs and 3'UTRs than the other genes (Figure 2E). Since YTHDF1 binding of m⁶A sites facilitate translation via recruitment of initiation factors and circularization of host mRNAs (14), the m⁶A modifications on CDSs and 3'UTRs are likely to assist with the recruitment and circularization processes (than those modifications in regions such as the 5'UTR).

m⁶A genes, including oncogenes, are unfavorable for translation

Circularization of mRNA is important for translation. The possibility that the mRNA length may be an important factor that influences circularization is intuitive. We speculated that the longer genes had more difficulty with circularization and therefore were less favorable for translation (this assumption is tested in the next section). Among the highly expressed m⁶A genes in HeLa cells, we globally profiled the relationship between the mRNA length (the genes were divided into bins) and the fraction of m⁶A sites in the CDSs and 3'UTRs. Interestingly, the fraction of m⁶A sites in the CDSs and 3'UTRs increased with the mRNA length (Figure 3A, P<2.2e-16). This result suggests that longer m⁶A genes have a greater need to promote their translation by methylation. Next, we compared the mRNA lengths of the m⁶A genes versus non-m⁶A genes (Figure 3B) or oncogenes versus other genes (Figure 3C). The results are as follows: (I) m⁶A genes are significantly longer than non-m⁶A genes (Figure 3B), and (II) among these two gene sets, the oncogenes are significantly longer than the other genes (Figure 3C). If longer genes are indeed unfavorable for translation, then the m⁶A genes (and especially the oncogenes among them) will suffer from a disadvantage in translation.

Figure 3 m⁶A genes, including oncogenes, are unfavorable for translation. (A) Correlation between the mRNA length and the fraction of m⁶A sites in the CDS and 3'UTR. The highly expressed genes in HeLa cells were divided into twenty bins with increasing mRNA lengths (X-axis); (B) the mRNA lengths of the m⁶A genes and non-m⁶A genes. The Wilcoxon rank sum test was used to calculate the P values; (C) the mRNA lengths of oncogenes and non-oncogenes and of m⁶A and non-m⁶A genes. The Wilcoxon rank sum test was used to calculate the P values; (D) codon bias of codons enriched in m⁶A genes compared to non-m⁶A genes. The codons are listed on the right in the order of codon bias; (E) codon bias of m⁶A genes and non-m⁶A genes. The Wilcoxon rank sum test was used to calculate the P values; (F) codon bias of oncogenes and non-oncogenes in the m⁶A or non-m⁶A genes. The Wilcoxon rank sum test was used to calculate the P values. m⁶A+ represents m⁶A genes, and m⁶A− represents non-m⁶A genes.

Another factor that can influence the translation rate is (synonymous) codon usage bias. During translation elongation, the rate-limiting step of the decoding process is waiting for the corresponding tRNA of each codon. A relationship should exist between codon usage and the translation efficiency. We followed the method used an early study (29) to calculate the codon bias of each human gene and each codon (see Methods for details). For m⁶A genes versus non-m⁶A genes in HeLa cells, we found that the codons enriched in the m⁶A genes had a stronger bias (Figure 3D). At the gene level, the m⁶A genes had a stronger bias than the non-m⁶A genes (Figure 3E), whereas the oncogenes showed almost no difference in codon bias compared to that of the other genes (Figure 3F). If a stronger codon bias is unfavorable for translation (tested in the following section), then the m⁶A genes should find a solution to neutralize this disadvantage.

m⁶A methylation facilitates the translation of host genes, including oncogenes

We have proposed that a longer mRNA length and stronger codon bias may be unfavorable for translation and that m⁶A genes may suffer from these disadvantages. Here, using mRNA-Seq and Ribo-Seq NGS data from HeLa cells generated by a previous study (14), we examined the correlation between the translation efficiency and CDS length (Figure 4A) or codon bias (Figure 4B). Both variables show a negative correlation with the translation efficiency (P<2.2e-16). These correlations verified our assumption that both long mRNA (CDS) lengths and strong codon bias were unfavorable for translation. Moreover, the length and codon bias did not show any correlations (Figure 4C), proving that these two factors might contribute independently to the lower translation efficiency.

Figure 4 m⁶A methylation facilitates the translation of host genes, including oncogenes. (A) Correlation between the translation efficiency and CDS length of highly expressed genes in normal HeLa cells; (B) correlation between the translation efficiency and codon bias of highly expressed genes in normal HeLa cells; (C) correlation between codon bias and the CDS length of highly expressed genes in normal HeLa cells; (D) the log2 fold-change of the translation efficiency (TE) in si-METTL3 versus normal HeLa cells. m⁶A genes and non-m⁶A genes were compared separately. The Wilcoxon rank sum test was used to calculate the P values; (E) dot plot displaying the translation efficiency (TE) of m⁶A genes and non-m⁶A genes in si-METTL3 versus normal HeLa cells; (F) genes belonging to both m⁶A genes and oncogenes; genes with the most decreased translation efficiencies in the si-METTL3 condition are displayed. In this Figure, m⁶A+ represents m⁶A genes, and m⁶A− represents non-m⁶A genes.

A question arises that since m⁶A genes (especially the oncogenes) are less suitable for translation, can they compensate for these disadvantages by m⁶A modification? We fully utilized data from normal (Control) and si-METTL3 (Treated) HeLa cells. Knock down (si-) of the m⁶A writer gene METTL3 largely reduced the transcriptome-wide m⁶A level (14). We compared the translation efficiency of m⁶A and non-m⁶A genes in normal or si-METTL3 HeLa cells. As expected, the m⁶A genes but not the non-m⁶A genes showed a reduced translation efficiency in the si-METTL3 condition (Figure 4D,E). Notably, the translation efficiency of the m⁶A genes was remarkably lower than that of the non-m⁶A genes when the m⁶A writer was removed, but the translation efficiencies of the m⁶A genes were still lower, even with help from m⁶A (in the normal condition) (Figure 4D). Take together with our previous results, we propose that the m⁶A genes are unfavorable for translation due to their longer lengths and stronger codon bias and that they indeed have low translation efficiencies. With the help of the m⁶A modification, translation of the target genes is elevated.

Since we showed that m⁶A genes were enriched in oncogenes (Figure 2D), we searched for oncogenes that benefited from m⁶A modification. We listed the oncogenes with the most decreased translation efficiencies in the si-METTL3 versus normal condition (Figure 4F). This set of oncogenes increased their translation efficiencies through m⁶A methylation and might play important roles in cancer cell oncogenesis.

Discussion

m⁶A methylation participates in many biological processes, of which one of the most well-studied functions is facilitating the translation of host genes via the reader protein YTHDF1 (14). Although several studies have mentioned the role of m⁶A in oncogenesis (10,19,20), they are either case studies of particular genes or do not systematically investigate the translation of target genes. In this work, we utilized the YTHDF1 target genes and normal or si-METTL3 NGS data (mRNA-Seq and Ribo-Seq) in HeLa cells generated by a previous work (14) and elucidated the potential function of the m⁶A modification in cancer cells.

We found that the m⁶A genes were enriched in oncogenes when compared to non-m⁶A genes. We observed remarkably longer mRNA lengths for m⁶A genes, especially the oncogenes among them. We also observed stronger codon usage bias for m⁶A genes than for non-m⁶A genes. Lines of evidence revealed that a longer mRNA length and stronger codon bias were unfavorable for translation, because these two factors were significantly negatively correlated with the translation efficiencies of the genes. This finding provides a simple explanation for why m⁶A genes (or the oncogenes among them) need the m⁶A modification to enhance their translation. Indeed, the translation efficiency of m⁶A genes is significantly elevated in normal HeLa cells (Control) compared to that in si-METTL3 HeLa cells (Treated) where methylation is removed, whereas the translation of non-m⁶A genes is almost unchanged in the control versus treated cells. In other words, the unfavorable features for translation of m⁶A genes (or the oncogenes among them) are compensated by the m⁶A modification (Figure 5). Recruitment of initiation factors does not directly resolve the codon bias problem, but the increased translation initiation rate will definitely enhance the global translation efficiency of host genes. We should note that the impact of codon usage bias on the mRNA translation efficiency is still debatable (31). However, we do not wish to contradict any previous reports. Conservatively, we declare that the translation efficiency of a gene is negatively correlated with its codon bias, at least for the HeLa cell data used in this study.

Figure 5 Model illustrating that the features unfavorable for mRNA translation are compensated by m⁶A methylation. Short mRNAs are favorable for loop formation and translation initiation, whereas long mRNAs are not. Genes with biased codon usage may have a low elongation efficiency due to tRNA availability. With help from the m⁶A modification and recognition by the m⁶A reader proteins, translation initiation factors are recruited, and global translation of host genes is enhanced. The eIF (shown in dark green) in the diagram represents a general but not a particular eukaryotic initiation factor.

Notably, compensation by m⁶A modification may originally have been designed for oncogenes, since we have observed enrichment of m⁶A genes in oncogenes. If methylated oncogenes obtain higher translation efficiencies and eventually facilitate cancer cell proliferation, this strategy or mechanism may explain how cancer cells/tissues achieve rapid cell growth. However, direct evidence for this assumption is still lacking. Hopefully, detailed experimental validation will be carried out in the future.

The main contribution of this study is to unveil the enrichment of m⁶A genes in oncogenes and to clarify why the m⁶A target is needed to enhance their translation efficiencies by m⁶A methylation (Figure 5). Understanding this relationship is important and should be interesting for the fields of RNA modification, translational regulation and cancer studies.

Conclusions

Our results demonstrate that HeLa cells compensate genes unfavorable for translation by m⁶A modification and enhance their translation efficiencies (Figure 5). These m⁶A target genes are enriched in oncogenes, and the enhancement of translation of these genes may be related to cancer cell oncogenesis.

Acknowledgments

Funding: This research was financially supported by the National Natural Science Foundation of China (Grant no. 31770213).

Footnote

Conflicts of Interest: Both authors have completed the ICMJE uniform disclosure form (available at http://dx.doi.org/10.21037/tcr.2019.03.04). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Cantara WA, Crain PF, Rozenski J, et al. The RNA Modification Database, RNAMDB: 2011 update. Nucleic Acids Res 2011;39:D195-201. [Crossref] [PubMed]
2. Buchumenski I, Bartok O, Ashwal-Fluss R, et al. Dynamic hyper-editing underlies temperature adaptation in Drosophila. PLoS Genet 2017;13:e1006931. [Crossref] [PubMed]
3. Liscovitch-Brauer N, Alon S, Porath HT, et al. Trade-off between Transcriptome Plasticity and Genome Evolution in Cephalopods. Cell 2017;169:191-202.e11. [Crossref] [PubMed]
4. Yu Y, Zhou H, Kong Y, et al. The Landscape of A-to-I RNA Editome Is Shaped by Both Positive and Purifying Selection. PLoS Genet 2016;12:e1006191. [Crossref] [PubMed]
5. Zhang R, Deng P, Jacobson D, et al. Evolutionary analysis reveals regulatory and functional landscape of coding and non-coding RNA editing. PLoS Genet 2017;13:e1006563. [Crossref] [PubMed]
6. Porath HT, Knisbacher BA, Eisenberg E, et al. Massive A-to-I RNA editing is common across the Metazoa and correlates with dsRNA abundance. Genome Biol 2017;18:185. [Crossref] [PubMed]
7. Ramaswami G, Lin W, Piskol R, et al. Accurate identification of human Alu and non-Alu RNA editing sites. Nat Methods 2012;9:579-81. [Crossref] [PubMed]
8. Ramaswami G, Zhang R, Piskol R, et al. Identifying RNA editing sites using RNA sequencing data alone. Nat Methods 2013;10:128-32. [Crossref] [PubMed]
9. Haussmann IU, Bodi Z, Sanchez-Moran E, et al. m(6)A potentiates Sxl alternative pre-mRNA splicing for robust Drosophila sex determination. Nature 2016;540:301-4. [Crossref] [PubMed]
10. Huang H, Weng H, Sun W, et al. Recognition of RNA N(6)-methyladenosine by IGF2BP proteins enhances mRNA stability and translation. Nat Cell Biol 2018;20:285-95. [Crossref] [PubMed]
11. Kan L, Grozhik AV, Vedanayagam J, et al. The m(6)A pathway facilitates sex determination in Drosophila. Nat Commun 2017;8:15737. [Crossref] [PubMed]
12. Lence T, Akhtar J, Bayer M, et al. m(6)A modulates neuronal functions and sex determination in Drosophila. Nature 2016;540:242-7. [Crossref] [PubMed]
13. Wang X, Lu Z, Gomez A, et al. N6-methyladenosine-dependent regulation of messenger RNA stability. Nature 2014;505:117-20. [Crossref] [PubMed]
14. Wang X, Zhao BS, Roundtree IA, et al. N(6)-methyladenosine Modulates Messenger RNA Translation Efficiency. Cell 2015;161:1388-99. [Crossref] [PubMed]
15. Xu K, Yang Y, Feng GH, et al. Mettl3-mediated m(6)A regulates spermatogonial differentiation and meiosis initiation. Cell Res 2017;27:1100-14. [Crossref] [PubMed]
16. Zheng Y, Nie P, Peng D, et al. m6AVar: a database of functional variants involved in m6A modification. Nucleic Acids Res 2018;46:D139-45. [Crossref] [PubMed]
17. Zhou J, Wan J, Gao X, et al. Dynamic m(6)A mRNA methylation directs translational control of heat shock response. Nature 2015;526:591-4. [Crossref] [PubMed]
18. Zhao BS, Nachtergaele S, Roundtree IA, et al. Our views of dynamic N(6)-methyladenosine RNA methylation. RNA 2018;24:268-72. [Crossref] [PubMed]
19. Liu J, Eckert MA, Harada BT, et al. m(6)A mRNA methylation regulates AKT activity to promote the proliferation and tumorigenicity of endometrial cancer. Nat Cell Biol 2018;20:1074-83. [Crossref] [PubMed]
20. Li Z, Weng H, Su R, et al. FTO Plays an Oncogenic Role in Acute Myeloid Leukemia as a N(6)- Methyladenosine RNA Demethylase. Cancer Cell 2017;31:127-41. [Crossref] [PubMed]
21. Lin S, Choe J, Du P, et al. The m(6)A Methyltransferase METTL3 Promotes Translation in Human Cancer Cells. Mol Cell 2016;62:335-45. [Crossref] [PubMed]
22. Cao Y, Cong H, Wang XD, et al. Circular RNAs in cancer: old tree with new flowers. Translational Cancer Research 2017;6:843-51. [Crossref]
23. Lee CY, Chiu YC, Wang LB, et al. Common applications of next-generation sequencing technologies in genomic research. Translational Cancer Research 2013;2:33-45.
24. Cingolani P, Platts A. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly (Austin) 2012;6:80-92. [Crossref] [PubMed]
25. Dobin A, Davis CA, Schlesinger F, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 2013;29:15-21. [Crossref] [PubMed]
26. Anders S, Pyl PT, Huber W. HTSeq--a Python framework to work with high-throughput sequencing data. Bioinformatics 2015;31:166-9. [Crossref] [PubMed]
27. Xiao Z, Zou Q, Liu Y, et al. Genome-wide assessment of differential translations with ribosome profiling data. Nat Commun 2016;7:11194. [Crossref] [PubMed]
28. Jiao X, Sherman BT. DAVID-WS: a stateful web service to facilitate gene/protein list analysis. Bioinformatics 2012;28:1805-6. [Crossref] [PubMed]
29. Akashi H. Inferring weak selection from patterns of polymorphism and divergence at "silent" sites in Drosophila DNA. Genetics 1995;139:1067-76. [PubMed]
30. Liu Z, Zhang J. Most m6A RNA modifications in protein-coding regions are evolutionarily unconserved and likely nonfunctional. Mol Biol Evol 2017; [Epub ahead of print]. [PubMed]
31. Yang JR, Chen X, Zhang J. Codon-by-codon modulation of translational speed and accuracy via mRNA folding. PLoS Biol 2014;12:e1001910. [Crossref] [PubMed]