Global loss of cellular m6A RNA methylation following infection with different SARS-CoV-2 variants

  1. Tanmoy Mondal1,8
  1. 1Department of Laboratory Medicine, Institute of Biomedicine, University of Gothenburg, Gothenburg 41345, Sweden;
  2. 2Virus and Immunity Unit, Institut Pasteur, Université Paris Cité, CNRS UMR 3569, Paris, France;
  3. 3GA08/NRG-Laboratory, Department of Molecular Reproduction, Development and Genetics, Indian Institute of Sciences, Bengaluru 560012, India;
  4. 4Department of Infectious Diseases, Sahlgrenska Academy, University of Gothenburg, Gothenburg 41345, Sweden;
  5. 5Department of Family Medicine, Gangnam Severance Hospital, Yonsei University College of Medicine, Gangnam-gu, Seoul 120-752, Korea;
  6. 6Department of Internal Medicine, Gangnam Severance Hospital, Yonsei University College of Medicine, Gangnam-gu, Seoul 120-752, Korea;
  7. 7Department of Medical Biochemistry and Biophysics, Karolinska Institute, Stockholm 17177, Sweden;
  8. 8Department of Clinical Chemistry, Sahlgrenska University Hospital, University of Gothenburg, Gothenburg 41345, Sweden

  1. ↵9 These authors contributed equally to this work.

  • Corresponding author: tanmoy.mondal{at}gu.se
    •  
    Next Section

    Abstract

    Insights into host–virus interactions during SARS-CoV-2 infection are needed to understand COVID-19 pathogenesis and may help to guide the design of novel antiviral therapeutics. N6-Methyladenosine modification (m6A), one of the most abundant cellular RNA modifications, regulates key processes in RNA metabolism during stress response. Gene expression profiles observed postinfection with different SARS-CoV-2 variants show changes in the expression of genes related to RNA catabolism, including m6A readers and erasers. We found that infection with SARS-CoV-2 variants causes a loss of m6A in cellular RNAs, whereas m6A is detected abundantly in viral RNA. METTL3, the m6A methyltransferase, shows an unusual cytoplasmic localization postinfection. The B.1.351 variant has a less-pronounced effect on METTL3 localization and loss of m6A than did the B.1 and B.1.1.7 variants. We also observed a loss of m6A upon SARS-CoV-2 infection in air/liquid interface cultures of human airway epithelia, confirming that m6A loss is characteristic of SARS-CoV-2-infected cells. Further, transcripts with m6A modification are preferentially down-regulated postinfection. Inhibition of the export protein XPO1 results in the restoration of METTL3 localization, recovery of m6A on cellular RNA, and increased mRNA expression. Stress granule formation, which is compromised by SARS-CoV-2 infection, is restored by XPO1 inhibition and accompanied by a reduced viral infection in vitro. Together, our study elucidates how SARS-CoV-2 inhibits the stress response and perturbs cellular gene expression in an m6A-dependent manner.

    The ongoing pandemic of COVID-19 disease is caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). Since 2019, the virus has spread over the world, and novel variants have emerged in succession, associated with changes in transmissibility and disease severity capacity (Tao et al. 2021). The SARS-CoV-2 genome consists of a single-stranded positive genomic RNA of approximately 30,000 nucleotides (Wu et al. 2020). Although different vaccines have been developed to mitigate SARS-CoV-2 spreading, it is still necessary to gain knowledge on the basic mechanisms underlying SARS-CoV-2 infection and host response to assist in the development of improved therapeutic options.

    N6-Methyladenosine modification (m6A) is a prevalent internal RNA modification (Baquero-Perez et al. 2021). It is involved in the regulation of a broad range of biological processes including cell differentiation (Geula et al. 2015), mRNA stability, translation, liquid-phase separation, and stress granule formation, among others (Zaccara et al. 2019). The writing of m6A on mRNA is mediated by a methyltransferase complex comprising the core catalytic subunits METTL3 and METTL14 and the adapter subunits WTAP, VIRMA, HAKAI, and RBM15/B (Ping et al. 2014; Patil et al. 2016; Xiao et al. 2016; Yue et al. 2018; Bawankar et al. 2021). The m6A-modified RNAs are recognized by cytoplasmic YTH domain-containing proteins, namely, YTHDF1–YTHDF3, which regulate mRNA stability and stress granule formation. In addition, the nuclear m6A reader proteins YTHDC1 and HNRNPA2B1 play critical roles in splicing and nuclear export (Alarcón et al. 2015; Xiao et al. 2016). Finally, m6A modifications can also be erased by the dioxygenases FTO and ALKBH5, which specifically demethylate m6A RNA (Zhang et al. 2017; Wei et al. 2022).

    The m6A RNA modification has been identified in several viral genomic RNA, being first described for DNA viruses, such as the simian virus 40 (SV40), herpes simplex virus, and adenovirus 2, and also in retroviruses, such as Rous sarcoma virus and human immunodeficiency virus 1 (HIV-1) (Lu et al. 2018). Moreover, a role for m6A modification has been described for RNA viruses, such as dengue, West Nile fever, yellow fever, Zika, and hepatitis C viruses, in which it was involved in the suppression of viral gene expression and replication (Lichinchi et al. 2016; Kim and Siddiqui 2021). The SARS-CoV-2 genome is m6A-modified by host proteins, and the m6A modification is important for promoting viral replication and for limiting the host immune response (Li et al. 2021; Liu et al. 2021; Zhang et al. 2021b). Depletion of the cytoplasmic m6A readers, the YTHDF proteins, and the writer METTL3 suppress SARS-CoV-2 (and HCoV-OC43) replication. METTL3 inhibition mediated by a small-molecule inhibitor induces viral RNA synthesis suppression, suggesting that m6A modification is needed for optimal viral RNA expression (Burgess et al. 2021).

    Though recent studies reported the presence of m6A in the SARS-CoV-2 RNA genome (Kim et al. 2020; Burgess et al. 2021; Li et al. 2021; Liu et al. 2021; Zhang et al. 2021b), it has not been explored in detail how the host's m6A mRNA profile is affected during infection. To address this question, we studied the effect of SARS-CoV-2 infection in Vero cells and air/liquid interface (ALI) cultures of human airway epithelia, which are highly permissive to SARS-CoV-2 infection (Robinot et al. 2021; Wei et al. 2021b).

    Previous SectionNext Section

    Results

    Infection by different SARS-CoV-2 variants leads to deregulation of RNA catabolism–related genes

    To identify variant-specific differences in gene expression perturbations, we infected Vero cells with three SARS-CoV-2 variants: B.1, B.1.1.7, and B.1.351. At 48 h postinfection, spike positivity was detected in 94%–95% of infected cells for all the variants (Supplemental Fig. S1A). We isolated total RNA 48 h postinfection and evaluated gene expression changes using RNA sequencing (RNA-seq). We first checked viral RNA levels postinfection and observed that RNA-seq reads mapping to the viral genome were comparable for the three SARS-CoV-2 variants (Supplemental Fig. S1B). The expression levels analyzed per viral gene were also similar for the three variants (Supplemental Fig. S1C). Viral reads contributed 1.2%–2.8% of the total reads recovered from infected cells (adding reads from the viral and host genomes) (Supplemental Fig. S1D), which is consistent with an earlier report (Blanco-Melo et al. 2020). Next, we performed a differential gene expression analysis to identify genes with altered expression postinfection in Vero cells (Fig. 1A; Supplemental Data S1). We identified a considerable number of up-regulated (998) and down-regulated (950) genes that were modulated across the three variants during infection (Fig. 1B). We then performed a principal component analysis (PCA) and a regression analysis to evaluate similarities in gene expression patterns across the three variants. We observed variant-specific patterns in gene expression, with gene expression patterns being more similar between the B.1.1.7 and B.1.351 variants (Supplemental Fig. S1E,F). Gene expression changes postinfection with three SARS-CoV-2 variants overlapped with several publicly available data sets deposited since the outbreak of COVID-19 (Supplemental Fig. S1G). We also compared the gene expression data set post-SARS-CoV-2 infection in three different variants with one of the publicly available data sets (Riva et al. 2020), in which differential gene expression was analyzed 24 h post-SARS-COV-2 infection in the Vero cells, and we observed a moderate but significant positive correlation between our data and publicly available data (Supplemental Fig. S1H). Pathway analysis revealed that altered genes were enriched in pathways associated with protein localization to the membrane, RNA catabolic processes, and cilium assembly (Fig. 1C; Supplemental Data S2). This was evident in the pathway analysis performed on up- and down-regulated genes separately after viral infection. Up-regulated genes were enriched for pathways related to several RNA catabolic processes, especially in B.1 and B.1.1.7 variants, whereas down-regulated genes were enriched with cilium assembly and organization (Supplemental Fig. S1I; Supplemental Data S2). We then chose to systematically look at the genes involved in RNA catabolic processes that were frequently deregulated across variants during SARS-CoV-2 infection. To do so, we selected RNA catabolism–related pathways, which were enriched in deregulated genes in at least two variants and visualized the interactions of those genes using network analysis (Fig. 1D). We observed that the mean connectivity of nodes in the RNA catabolism pathway–related network was much higher compared with random networks (Supplemental Fig. S1J). In particular, the network analysis suggested widespread interactions (as detected by connectivity degree) of m6A-related genes with different RNA catabolic processes (Supplemental Fig. S1K). We observed that m6A-related genes were frequently deregulated during SARS-CoV-2 infection in Vero cells (Fig. 1E), although we did not observe any significant change in the expression of the main m6A writers METTL3 and METTL14. We validated the differential expression of m6A eraser FTO, which was down-regulated, and the up-regulation of the m6A-related genes SPEN and YTHDF1 (Dossin et al. 2020), which was observed after infection with the three variants. Similar changes were observed after infection with the B.1.617.2 (Delta) variant, which was available during the final compilation of the data (Fig. 1F). Although we observed variant-specific gene expression changes, our data suggest widespread deregulation of genes associated with RNA catabolic processes and m6A modification–related pathways for all the SARS-CoV-2 variants tested.

    Figure 1.
    View larger version:
    Figure 1.

    Gene expression changes postinfection with different variants of SARS-CoV-2. (A) Volcano plots showing gene expression changes and significance (−log10 scale) in SARS-CoV-2-infected versus noninfected cells. Dashed lines indicate the significance threshold (adjusted P-value < 0.01) and log2 fold-change threshold (abs log2 fold-change > 1). Differentially expressed genes are highlighted in red, and labels indicate some of the common differentially expressed genes across all conditions. (B) Venn diagram comparison of up- and down-regulated genes across variants. (C) Top enriched Gene Ontology (GO) term biological processes after infection. The size of the dots represents the enrichment of genes with a GO term, colored according to their significance level. (D) ClueGO clustering and visualization of common terms associated with m6A-related genes across any pair of SARS-CoV-2 variants. m6A genes that are deregulated across variants are highlighted in red. (E) Heat map of m6A-related genes, with gray boxes indicating nondifferentially expressed genes and colors indicating the log2 fold-change values of significantly deregulated genes in at least two variants. (F) Relative mRNA expression of m6A-related genes in noninfected Vero cells and Vero cells infected with B.1, B.1.1.7, B.1.351, and B.1.617.2 (Delta) variants of SARS-CoV-2. TBP and POL2RG were used to normalize the qPCR data. Data are shown as mean ± SD of three replicates (n = 3). Statistics: two-tailed paired t-test; (*) P < 0.05, (**) P < 0.01, (***) P < 0.001.

    Variant-specific changes in cellular RNA m6A levels after viral infection

    A change in m6A-related genes during viral infection prompted us to check the effect of infection on cellular m6A levels. To this end, isolated RNA from noninfected Vero cells and cells infected with three SARS-CoV-2 variants (B.1, B.1.1.7, and B.1.351) was supplemented with spike-in bacterial RNA as a control (spike-in) for m6A RNA immunoprecipitation (RIP) followed by sequencing (m6A RIP-seq) (Fig. 2A). We observed a loss of m6A peaks in the cellular RNA postinfection with different SARS-CoV-2 variants, with the loss of m6A peaks being more marked in the B.1 and B.1.1.7 infections compared with a B.1.351 infection (Supplemental Fig. S2A; Supplemental Data S3). The global distribution of m6A peak density across transcripts in noninfected cells was similar to that reported previously, with strong enrichments near the start and stop codons (Meyer et al. 2012). Spike-in normalized relative m6A peak density showed m6A loss across the whole transcript length following SARS-CoV-2 infection in Vero cells, with more pronounced effects in B.1 and B.1.1.7 infection compared with B.1.351 infection (Fig. 2B). The input RNA signal across the transcript length was largely unaltered between noninfected and different SARS-CoV-2-infected conditions, suggesting input RNA does not contribute to the observed m6A loss (Fig. 2B; Supplemental Fig. S2B). This global loss of m6A was also evident by the reduced number of genes with m6A peaks following SARS-CoV-2 infection (Supplemental Fig. S2C). A decrease in the spike-in normalized m6A signal was visible over the commonly lost m6A peaks (Supplemental Fig. S2D) and also at the level of individual transcripts (Fig. 2C) after infection with three SARS-CoV-2 variants. Identified m6A peaks in noninfected cells were enriched with the known “GGAC” motif, whereas the dominant motif was more degenerate in infected cells (Fig. 2D). Although the number of m6A peaks was drastically reduced, we also detected new m6A peaks that were gained over many transcripts after infection with all three variants. The number of gained m6A peaks was highest for the B.1.351 variant (Fig. 2E; Supplemental Fig. S2C). When we investigated the genomic location of the lost and gained peaks during infection, we observed that gained peaks were found in different genomic locations, with major gained peaks located at intergenic and TSS regions (Supplemental Fig. S2F). The most enriched motifs (despite less-significant P-values) in the gained m6A peaks were different from the canonical “GGAC” motif (Supplemental Fig. S2E). We also checked by scatter plots the relation between expression (log2 fold-changes: infected vs. noninfected) and m6A enrichment (infected vs. noninfected) across the transcripts and found that change in expression cannot simply explain alteration in m6A enrichment (Supplemental Fig. S2G). This was further evident in the input RNA signal, which was mostly unaffected over the common lost peaks in infected cells even though the m6A signal was reduced over those peaks (Supplemental Fig. S2D), suggesting the involvement of additional mechanisms that drive m6A loss during SARS-CoV-2 infection.

    Figure 2.
    View larger version:
    Figure 2.

    Variant-specific changes in cellular RNA m6A level after viral infection. (A) Flow chart describing the m6A RIP-seq protocol. (B) Metagene analysis showing spike-in normalized relative m6A peak density at host genes in noninfected and SARS-CoV-2-infected Vero cells (variants as specified by colors). (C) Genome browser screenshots showing the ratio of m6A RIP/input for three different genes in noninfected and SARS-CoV-2-infected (variants as specified) Vero cells. Ratio tracks were calculated using spike-in normalized counts per million (CPM). m6A peak regions identified using the MACS peak caller software in noninfected cells are highlighted using gray boxes. (D) Identified motifs from de novo motif analysis of peaks in noninfected and infected cells. (E) The total number of peaks classified as retained, gained, or lost postinfection with the SARS-CoV-2 variants are indicated and compared with the noninfected condition. (F–H) METTL3 localization. (F) METTL3 and α-tubulin coimmunostaining in Vero cells that were either noninfected or infected with SARS-CoV-2 variants (as denoted). Scale bar is 20 μm. (G) METTL3 intensities (in relative fluorescence units [RFU]) derived from the immunostainings performed in F. METTL3 RFU in the nucleus (red) and cytoplasm (blue) were estimated using ImageJ, using DAPI marking the nucleus and α-tubulin staining the cytoplasm as references. Data are shown as mean ± SD. Data presented from multiple experiments with the total number of cells counted, n = <70. One-way ANOVA test was performed; (**) P < 0.01, (***) P < 0.001, (****) P < 0.0001. (H) Distribution of METTL3 in the nucleus and cytoplasm calculated from the percentage of RFU intensities measured as described in G. Data are shown as mean ± SD. (I,J) METTL3–METTL14 interaction. (I) Proximity ligation assay (PLA) in noninfected and SARS-CoV-2-infected Vero cells (variants as indicated) depicting METTL3 and METTL14 PLA foci in the nucleus (marked by DAPI). The background control shows PLA with only the METTL3 antibody. Scale bar is 20 μm. (J) Quantification of METTL3–METTL14 PLA foci as detected in I. The number of PLA foci/nuclei are shown as mean ± SD. Data presented from multiple experiments with the total number of cells counted, n = <100. One-way ANOVA test was performed; (****) P < 0.0001.

    The global loss of m6A peaks during infection suggests an altered function of the key cellular enzyme complex (METTL3/METTL14), which deposits m6A modifications. The METTL3/METTL14 complex is normally localized in the nucleus and deposits m6A modification cotranscriptionally (Huang et al. 2019). We determined whether METTL3/METTL14 function could be compromised during infection owing to altered cellular localization. We observed that during SARS-CoV-2 infection, METTL3 was partially relocalized from the nucleus to the cytoplasm (Fig. 2F). Importantly, cytoplasmic localization of METTL3 was more evident with the B.1 and B.1.1.7 variants compared with B.1.351 (Fig. 2F–H), which is consistent with a greater reduction in m6A level in B.1 and B.1.1.7 infections. To verify that the METTL3 antibody used was specific and did not detect a viral antigen in infected cells, we infected both control (Control-sh) and METTL3 knockdown (KD) Vero cells with the B.1 variant (Supplemental Fig. S2H,I). We observed a drastic reduction in METTL3 signal upon METTL3 KD in noninfected cells. After infection with B.1, the Control-sh cells showed METTL3 localization in both nuclear and cytoplasmic compartments; however, in the METTL3 KD condition, this staining was absent, suggesting that the METTL3 antibody was indeed specific (Supplemental Fig. S2J). The nuclear localization of METTL14 was mostly unaffected during SARS-CoV-2 infection (Supplemental Fig. S2K). As m6A loss was marked in infected cells, we wondered if partial METTL3 cytoplasmic localization was sufficient to compromise METTL3/METTL14 functional complex formation. To test this, we performed a proximity ligation assay (PLA) to detect the METTL3/METTL14 complex in the noninfected and SARS-CoV-2-infected cells. We observed that the METTL3/METTL14 PLA signal was decreased in infected cells, with an effect that was more marked in B.1 and B.1.1.7 compared with B.1.351, in which METTL3 localization was also less affected (Fig. 2I,J).

    SARS-CoV-2 genomic RNA contains m6A modification

    Recently several studies reported the presence of m6A modification on SARS-CoV-2 genomic RNA (Burgess et al. 2021; Li et al. 2021; Liu et al. 2021; Zhang et al. 2021b). As we observed that the key m6A depositing enzyme METTL3 showed partial relocalization to the cytoplasm where SARS-CoV-2 genome replication occurs, we investigated the m6A profile in viral RNAs. Our strand-specific m6A RIP sequencing (m6A RIP-seq) data allowed us to profile m6A in both the positive (genomic) and the negative (replicative intermediates) strands of viral RNA for the three SARS-CoV-2 variants. Consistent with previous reports, we detected m6A-enriched regions at several positions of the positive-strand SARS-CoV-2 genome for all three variants, with broader m6A peaks detected in the N gene region (Fig. 3A; Supplemental Data S3; Burgess et al. 2021; Li et al. 2021; Liu et al. 2021; Zhang et al. 2021b). Using liquid chromatography–tandem mass spectrometry (LC-MS/MS), we then determined the presence of a relative number of m6A residues in the SARS-CoV-2 RNA genome. In these experiments, viral genomic RNA was isolated from infected Vero cell supernatant, and the ratio of A/m6A was measured using a synthetic m6A-containing RNA oligo as a standard. Using the parallel reaction-monitoring (PRM) mode, we monitored LC-MS/MS profiles for m/z 268.0–136.0 and m/z 282.0–150.1, which correspond to A and m6A, respectively (Sun et al. 2021). We observed that each SARS-CoV-2 genome contains, on average, 10 m6A modifications (Supplemental Fig. S3A), consistent with a recent report (Li et al. 2021). We observed m6A-containing peaks were also present in the negative strand of the viral RNA, mainly in the N/ORF10 region (Fig. 3B). We then validated the presence of m6A in both the positive and negative strands by strand-specific m6A RIP-qPCR for the B.1.1.7 variant (Fig. 3C). Further, using SARS-CoV-2-infected patient RNA samples from throat/nose swabs, we could detect m6A enrichment in the N gene region in both the positive and negative strands of SARS-CoV-2 RNAs (Fig. 3D). We verified up-regulation of selected interferon-stimulated genes in SARS-CoV-2-infected patient throat/nose swab RNA as previously reported (Supplemental Fig. S3B; Gao et al. 2021; Lorè et al. 2021b). The presence of m6A peaks in both strands of SARS-CoV-2 RNA suggests an important functional role for such modifications (Burgess et al. 2021; Li et al. 2021; Liu et al. 2021; Zhang et al. 2021b). Inspection of publicly available viral-RNA protein interaction data shows that m6A reader protein partners are enriched in viral RNA interacting proteins (Schmidt et al. 2021). The m6A readers previously shown to interact with SARS-CoV-2 RNA were searched for known interacting protein partners in the STRING database. We found that these interacting proteins are enriched in RNA metabolism and viral process pathways (Supplemental Fig. S3C, top) and that these proteins show overlap with the SARS-CoV-2 RNA–protein interactome (Supplemental Fig. S3C, bottom). These observations suggest that SARS-CoV-2 may use m6A residues in its genome to recruit other RNA-binding proteins to the viral RNA using m6A readers as intermediates. This is consistent with the recently reported role of m6A reader YTHDF proteins in SARS-CoV-2 infection (Burgess et al. 2021). Given m6A has been implicated in other RNA viruses, the detailed characterization of interactions between m6A readers and SARS-CoV-2 RNA will further elucidate the role of m6A reader proteins in host–virus interactions (Lichinchi et al. 2016; Kim and Siddiqui 2021).

    Figure 3.
    View larger version:
    Figure 3.

    SARS-CoV-2 genomic RNA contains m6A modification. (A) Identified m6A peaks on SARS-CoV-2 positive strands from B.1, B.1.1.7, and B.1.351 variants are shown. The m6A peaks in the positive strand of SARS-CoV-2 RNA from publicly available data. (B) m6A peaks on the SARS-CoV-2 negative strand. The presence of m6A peaks in the N gene region is highlighted, and normalized coverage is shown on the panel below. (C,D) m6A-RIP qPCR data showing enrichment of N gene region in both the positive strand and the negative strand of viral RNA. (C) Twenty-four hours postinfection with B.1.1.7 in Vero cells; data are represented as a percentage of input, and IgG was used as the negative control. Statistics: two-tailed paired t-test, n = 3; (*) P < 0.05. (D) m6A-RIP qPCR data showing enrichment of the N gene region in both the positive strand and the negative strand of viral RNA from the infected patient samples (n = 5). Data are represented as a percentage of input. IgG was used as a negative control. Statistics: two-tailed paired t-test; (**) P < 0.01.

    m6A loss modulates cellular gene expression in SARS-CoV-2-infected cells

    To better understand the effects of m6A modulation on cellular gene expression during viral infection, we first identified the set of transcripts showing a change in m6A level during infection with the three SARS-CoV-2 variants compared with noninfected cells. We then compared the genes with lost m6A (shows loss of at least one m6A peak compared with noninfected cells) or with gained m6A (no m6A peak present in noninfected cells, but one or more m6A peaks gained after infection) with our set of differentially expressed genes (up- and down-regulated during infection). We found that the list of up- and down-regulated genes overlapped with that of genes with lost and gained m6A. In general, down-regulated genes showed higher overlap with lost m6A genes across the three variants, whereas the up-regulated genes were overrepresented in the gained m6A gene category (Fig. 4A). To better understand how loss and gain of m6A contributed to gene expression, we performed a cumulative distribution function (CDF) analysis that measured RNA abundance using RNA-seq data in genes with either lost or gained m6A peaks. We found that lost m6A genes showed a decrease in abundance, whereas gained m6A genes were increased in abundance after SARS-CoV-2 infection for the three variants studied (Supplemental Fig. S4A). Similar findings were obtained when genes with m6A gain or loss were compared with genes with no change in m6A level (retained genes) (Supplemental Fig. S4B). We further categorized genes in noninfected cells based on m6A content to either non-m6A, low (one peak), medium (two peaks), or high (three or more peaks) m6A. We observed higher-m6A-containing genes were more susceptible to viral infection, with a decreased expression of these genes upon infection with SARS-CoV-2 variants (Fig. 4B). This suggests that m6A modification on these genes contributes to their expression and that their expression is compromised owing to the global loss of m6A during viral infection. We chose a few genes with varying levels of m6A in noninfected cells and validated their down-regulation during viral infection (Fig. 4C). Expression of these genes was also decreased in METTL3 KD cells, indicating an m6A-dependent expression (Supplemental Fig. S4C). The exception was EID1, which did not change expression in absence of METTL3, consistent with this gene having no m6A modification (Supplemental Fig. S4C). Infection with the B.1.617.2 (Delta) variant also down-regulated the expression of the five tested genes (Fig. 4C). We checked the localization of METTL3 after B1.617.2 infection and observed that METTL3 was partially relocalized from the nucleus to the cytoplasm in B.1.617.2 infection as well (Supplemental Fig. S4D). We then performed a pathway analysis for the up- and down-regulated genes with lost m6A. We observed that down-regulated genes with lost m6A were mostly involved in cilium organization, whereas up-regulated genes with lost m6A were primarily part of pathways related to covalent chromatin/histone modifications (Supplemental Fig. S4E). We validated the down-regulation of some of the known cilia-related genes such as HDAC6 and PKHD1 (Fig. 4C; Zhang et al. 2004; Ran et al. 2015) consistent with a recent report suggesting widespread loss of motile cilia after SARS-CoV-2 infection (Robinot et al. 2021). The genes with lost m6A involved in covalent chromatin modification and up-regulated during infection, such as KDM6A (Supplemental Data S2), were recently identified as proviral genes in a CRISPR screen (Wei et al. 2021b). Thus, our data suggest that loss of m6A during SARS-CoV-2 infection establishes a pattern of host cell gene expression that favors viral replication.

    Figure 4.
    View larger version:
    Figure 4.

    Gene expression after viral infection and global m6A loss. (A, top) Overlap diagrams of differentially expressed and m6A-modified genes after infection. Each set shows the number of up-regulated and down-regulated genes associated with lost or gained m6A modifications. The remaining genes are shown in gray. Hypergeometric test P-values depict the calculated probability of overlap between differentially expressed and m6A lost/gained genes. (A, bottom) Bar-graph summarizing the percentage of up-/down-regulated genes overlapping with either lost or gained m6A genes. (B) Log2 fold-change distributions of differentially expressed genes postinfection with three different variants, categorized according to their m6A level: non, low, medium/high. Statistical significance was calculated using the Wilcoxon test. (ns) Nonsignificant, (*) P < 0.05, (**) P < 0.01, (***) P < 0.001, (****) P < 0.0001. (C) Relative mRNA expression of genes with varying levels of m6A (FAM111A and PKHD1, high m6A; HDAC6, medium m6A; RNASEL, low m6A; EID1, no m6A) after viral infection with the variant indicated compared with noninfected cells. qPCR data were normalized to TBP and POL2RG. Statistics: two-tailed t-test; (*) P < 0.05, (**) P < 0.01, (***) P < 0.001, n = 2.

    m6A modification is known to be implicated in alternative splicing (Wei et al. 2021a). Considering the global loss of m6A during SARS-CoV-2 infection, we looked for differential exon use (DEU) in the RNA-seq data. We identified about 790 DEU events common to cells infected with the three SARS-CoV-2 variants (Supplemental Data S4). For instance, we observed that exon 1 of the COL6A2 collagen gene is differentially included during infection, which is associated with the loss of m6A peaks surrounding the newly included exon 1 (Supplemental Fig. S4F). We also found that DEU-containing genes had on average more m6A peaks compared with random m6A-positive genes in noninfected cells (Supplemental Fig. S4G). Further, we observed that m6A-containing DEU genes showed a loss of m6A peaks during infection (Supplemental Fig. S4H).

    Treatment with selinexor restores METTL3 cellular localization during SARS-CoV-2 infection

    During the SARS-CoV-2 infection, we observed partial relocalization of METTL3 to the cytoplasm from its original nuclear location. Considering that m6A modification of cellular mRNA is deposited cotranscriptionally in the nucleus (Huang et al. 2019), we aimed to restore METTL3 nuclear localization to see if we could antagonize the effects of SARS-CoV-2 infection. We observed that the expression of XPO1 (or exportin 1), one of the major nuclear export proteins, was up-regulated during infection (Fig. 5A; Supplemental Fig. S5A). Therefore, we set to determine whether XPO1 was involved in the relocalization of METTL3 during viral infection. Using PLA, we observed that XPO1 and METTL3 interacted with each other in Vero cells (Fig. 5B). We then treated these cells during infection with selinexor, a well-characterized XPO1 inhibitor (Kashyap et al. 2021). We observed that selinexor treatment resulted in the restoration of METTL3 localization to the nucleus of B.1-infected cells (Fig. 5C,D). We then checked if the restoration of METTL3 localization could rescue METTL3/METTL14 complex formation, which was compromised in infected cells. We observed that the METTL3/METTL14 PLA signal was robustly detected in B.1- and B.1.1.7-infected cells treated with selinexor but not in infected cells treated with DMSO (Fig. 5E; Supplemental Fig. S5B). These results suggested that change in METTL3 localization perturbed the formation of the METTL3/METTL14 complex and that promoting METTL3 nuclear localization by selinexor was sufficient to drive the formation of the METTL3/METTL14 complex in infected cells. Importantly, selinexor treatment was also effective in decreasing SARS-CoV-2 infection, as measured by spike protein positivity for all the four variants tested (B.1, B.1.1.7, B.1.351, B.1.617.2), without decreasing cell viability (Supplemental Fig. S5C,D). Of note, the efficacy of selinexor in preventing SARS-CoV-2 infection has also been reported recently (Kashyap et al. 2021).

    Figure 5.
    View larger version:
    Figure 5.

    Treatment with selinexor restores METTL3 cellular localization during viral infection. (A) Relative expression of the XPO1 gene evaluated post-SARS-CoV-2 infection (variants as specified) using RT-qPCR. Data were normalized to TBP and POL2RG. Two-tailed t-test; (*) P < 0.05, (**) P < 0.01, n = 3. (B, top) PLA in noninfected Vero cells depicting interacting XPO1 and METTL3 foci. (Bottom) Background control for PLA with only XPO1 antibody. Scale bar is 20 μm. (C,D) Rescue of METTL3 localization. (C) METTL3 and α-tubulin coimmunostaining in Vero cells postinfection (24 h) with B.1 SARS-CoV-2 that were either treated with DMSO or selinexor (150 nM). Scale bar is 20 μm. (D, top) METTL3 intensities (RFU) from the immunostaining performed in C. METTL3 (RFU) in the nucleus and cytoplasm were estimated using ImageJ with help of DAPI marking the nucleus and α-tubulin staining the cytoplasm. Data are shown as Mean ± SD. Data presented from multiple experiments with number of cells counted, n = <70. One-way ANOVA test was performed; (**) P < 0.01, (***) P < 0.001. (D, bottom) Percentage of METTL3 distribution in nucleus and cytoplasm calculated from intensities from the upper panel. (E, top) PLA depicts interacting METTL3 and METTL14 foci in the nucleus (marked by DAPI) of Vero cells infected with B.1 SARS-CoV-2 variants and treated with either DMSO or selinexor (150 nM). Scale bar is 20 μm. (E, bottom) Quantification of METTL3–METTL14 PLA foci as detected in the above panel. The number of PLA foci/nuclei are shown as mean ± SD. Data presented from multiple experiments with number of cells counted, n = >100. One-way ANOVA test was performed; (****) P < 0.0001. (F) Immunostaining showing G3BP1 localization after DMSO or selinexor (150 nM) treatment in B.1-infected cells (24 h postinfection). White arrow highlights some of the G3BP1 foci. Scale bar is 50 μm. (G) Immunostaining showing G3BP1 and m6A/YTHDF2 localization after selinexor treatment in B.1-infected cells (24 h postinfection). White arrow highlights some of the G3BP1 foci overlapping with either m6A or YTHDF2. Scale bar is 10 μm. (H) G3BP1 and YTHDF2 RIP showing enrichment of m6A-positive genes in Vero cells. IgG was used as a negative control. Statistics: two-tailed t-test; (***) P < 0.001, (****) P < 0.0001. n = 4. (I) G3BP1 and YTHDF2 RIP showing enrichment of m6A-positive genes in Vero cells treated with either METTL3 inhibitor STM2457 (5 μM) or DMSO. IgG was used as a negative control. Two-tailed t-test; (*) P < 0.05, (**) P < 0.01; n = 3. (J) Relative mRNA expression of genes with varying levels of m6A after selinexor treatment in B.1.1.7-infected cells. Treatment with selinexor partially restores the expression of m6A-modified genes (FAM111A, HDAC6, RNASEL, and PKHD1), but not the non-m6A gene (EID1). qPCR data were normalized to TBP and POL2RG. Statistics: two-tailed t-test; (*) P < 0.05, (**) P < 0.01, (***) P < 0.001, (****) P < 0.0001; n = 3. (K) m6A-RIP qPCR showing recovery of m6A levels at selected genes in selinexor-treated B.1.1.7-infected cells. DMSO was used as treatment control and IgG was used as a negative control for RIP. Two-tailed t-test; (*) P < 0.05, (**) P < 0.01, (***) P < 0.001, (****) P < 0.0001; n = 3.

    m6A-modified RNA and m6A reader YTHDF proteins have been implicated in stress granule formation by promoting the recruitment of stress granule proteins such as G3BP1 (Fu and Zhuang 2020). Stress granule formation is known to be compromised in SARS-CoV-2-infected cells, and KD of the stress granule protein G3BP1 enhances SARS-CoV-2 infection (Zheng et al. 2021). We reasoned that the inactivation of stress granule formation by SARS-CoV-2 could be mediated by the global loss of m6A in cellular mRNAs that we uncovered in infected cells. Therefore, we tested if selinexor-mediated reversal of METTL3 nuclear localization could promote stress granule formation in infected cells by reinstalling m6A modification in cellular mRNAs. We observed the appearance of stress granules in selinexor-treated, but not in DMSO-treated (control), SARS-CoV-2-infected cells (Fig. 5F; Supplemental Fig. S5E). Immunostaining showed that stress granules in selinexor-treated infected cells often overlapped with m6A-modified RNAs and with the m6A reader protein YTHDF2 (Fig. 5G; Supplemental Fig. S5F). We then used an RIP assay to check the interaction of four validated mRNAs that show loss of m6A and down-regulation postinfection with the m6A reader YTHDF2 and the stress granule protein G3BP1. Indeed, we found that these mRNAs coprecipitated with both G3BP1 and YTHDF2 (Fig. 5H). Inhibiting the METTL3 enzyme by the small-molecule inhibitor STM2457 (Yankova et al. 2021) decreased the interaction of the four mRNAs with G3BP1 and YTHDF2, suggesting that m6A is required for such interaction (Fig. 5I). METTL3 inhibition did not alter the expression of G3BP1 and YTHDF2 (Supplemental Fig. S5G). Treatment with selinexor, which restored METTL3 nuclear localization during infection, also rescued the expression of the four candidate down-regulated genes, which included the HDAC6 and PKHD1 genes related to cilia formation (Fig. 5J). We also observed that selinexor treatment could reinstate the m6A modifications that were lost during infection in these four mRNAs (Fig. 5K). Therefore, rescuing nuclear localization of METTL3 by selinexor promoted stress granule formation in infected cells, restored cellular gene expression, and inhibited SARS-CoV-2 infection in vitro, highlighting the possibility of targeting the m6A pathway as an antiviral strategy.

    Global loss of m6A modification in primary human airway epithelial cells upon SARS-CoV-2 infection

    Next, we aimed to validate the global loss of m6A methylation during SARS-CoV-2 infection in human cell infection models. To this goal, we tested METTL3 localization post-SARS-CoV-2 infection in the bronchial epithelial cell line, BEAS-2B. BEAS-2B cells were positive for spike protein, and METTL3 showed partial relocalization from the nucleus to the cytoplasm in these infected cells (Fig. 6A; Supplemental Fig. S6A). Further, we also tested METTL3 localization postinfection in primary human bronchial epithelial (HBE) cells grown in monolayer culture, and we observed, similar to BEAS-2B, partial cytoplasmic relocalization of METTL3 (Supplemental Fig. S6B,C).

    Figure 6.
    View larger version:
    Figure 6.

    Loss of m6A in human airway epithelia after SARS-CoV-2 infection. (A, left) METTL3 localization in BEAS-2B cells postinfection with SARS-CoV-2. The scale bar is 20 μm. (A, right) The percentage distribution of METTL3 in the nucleus and cytoplasm was calculated using ImageJ with help of DAPI marking the nucleus and α-tubulin staining the cytoplasm. Data are shown as mean ± SD. Data presented from multiple experiments with the number of cells counted, n = <300. Unpaired t-test; (****) P < 0.0001. (B) Flow chart describing the experimental design for SARS-CoV-2 infection of reconstructed human bronchial/nasal epithelium at ALI. (C) Total number of m6A peaks classified as retained, gained, or lost, after 4- or 7-d post infection (dpi) with SARS-CoV-2, compared with noninfected HBE. (D) Bar plots representing the number of genes that lost, gained, or retained m6A in HBE after 4 and 7 dpi. (E) Metagene plots showing spike-in normalized relative m6A peak density distribution at host genes in noninfected and SARS-CoV-2-infected HBE cells at 4 and 7 dpi. (F) Boxplot showing log2 ratio of m6A signal over input obtained from spike-in normalized CPM data at lost m6A peak regions (±125 bp from m6A peak submit) taken from panel C at 4 and 7 dpi in HBE. Statistical significance was calculated using the Wilcoxon test; (****) P < 0.0001. (G) Genome browser screenshots showing the ratio of m6A RIP/input for three different genes in noninfected and SARS-CoV-2-infected HBE at 4 dpi. Ratio tracks were calculated using spike-in normalized CPM. m6A peak region identified using MACS peak caller in noninfected cells are highlighted using the gray box. (H) log2 fold-change distributions of differentially expressed genes in HBE at 4 dpi, categorized according to their m6A level: non, low, medium/high. Statistical significance was calculated using the Wilcoxon test; (**) P < 0.01, (***) P < 0.001, (****) P < 0.0001. (I) Total number of m6A peaks classified as retained, gained, or lost, after 4 dpi with SARS-CoV-2 compared with noninfected HNE. (J) The number of genes that lost, gained, or retained m6A in HNE post-SARS-CoV-2 infection. (K) Metagene plots showing spike-in normalized relative m6A peak density distribution in noninfected and SARS-CoV-2-infected HNE. (L) Boxplot showing log2 ratio of m6A signal over input obtained from spike-in normalized CPM data at lost m6A peak regions (±125 bp from m6A peak submit) taken from panel I at 4 dpi in HNE. Statistical significance was calculated using the Wilcoxon test; (****) P < 0.0001.

    We then used a reconstructed HBE as the SARS-CoV-2 infection model, in which cells are differentiated by the culture at the ALI, forming a pseudostratified ciliated epithelium before viral infection as described previously (Robinot et al. 2021). RNA was extracted from HBE at 4 and 7 days postinfection (dpi), and we observed high levels of SARS-CoV-2 replication as measured by N gene expression (Supplemental Fig. S6D). Extracted RNA from HBE at 4 and 7 dpi was supplemented with spike-in bacterial RNA and used in m6A RIP-seq experiments (Fig. 6B). SARS-CoV-2 reads contributed 0.48% and 0.47% of the total reads at 4 and 7 dpi, respectively (Supplemental Fig. S6E). Consistent with our observations in infected Vero cells, we observed a marked loss in the number of m6A peaks postinfection (Fig. 6C; Supplemental Fig. S6F; Supplemental Data S3). The top enriched motif in the m6A peaks from uninfected HBE contained the “GGAC” sequence, whereas the top motif was more degenerate in infected HBE (Supplemental Fig. S6G). A few transcripts showed a gain of m6A during infection in the HBE model, but the fraction of transcripts with gained m6A peaks was much lower compared with that in Vero cells (Fig. 6C). We observed a drastic decrease in the number of m6A-positive genes in HBE on both 4 and 7 dpi (Fig. 6D,E). m6A enrichment (IP/input) and spike-in normalized m6A signal over the lost peaks were significantly reduced in the SARS-CoV-2-infected cells, whereas input RNA was unchanged over these lost peaks (Fig. 6F; Supplemental Fig. S6H). Genome browser visualization of three chosen genes showed the loss of m6A peaks at 4 dpi (Fig. 6G). We could detect m6A peaks in both the positive and negative strands of viral RNA at 4 dpi. We also detected m6A peaks in the positive strand of viral RNA at 7 dpi, but m6A peaks were absent in the negative strand at day 7, possibly because the infection reaches a plateau or has already started to decrease at this stage, as described previously (Supplemental Fig. S6I; Supplemental Data S3; Robinot et al. 2021). Genes related to the defense response to viruses and type I interferon signaling pathways were up-regulated upon SARS-CoV-2 infection in HBE, whereas genes related to cilium organization were the most down-regulated (Supplemental Fig. S6J; Supplemental Data S1,S2). In the HBE infection model, genes containing one or more m6A peaks were prone to down-regulation postinfection, similar to the findings obtained in Vero cells (Fig. 6H). Down-regulated genes with lost m6A were again enriched in the pathways of cilium organization and cilium assembly (Supplemental Fig. S6K), suggesting that m6A loss could be a major contributor to the loss of cilia as previously reported in the HBE model (Robinot et al. 2021). These data collectively confirm that SARS-CoV-2 infection causes a general perturbation of m6A-dependent gene expression in human primary epithelial cells.

    We further evaluated changes in cellular m6A level following SARS-CoV-2 infection in reconstructed human nasal epithelia (HNE) cultures in ALI conditions, which also support robust SARS-CoV-2 infection (Fig. 6B; Samelson et al. 2022). We verified SARS-CoV-2 infection in HNE 4 dpi using the expression of the N gene (Supplemental Fig. S7A). SARS-CoV-2 reads contributed 0.79% of total reads in the HNE (Supplemental Fig. S7B). Infected HNE showed a drastic loss of m6A peaks 4 dpi and a decrease in the number of m6A-positive genes as well (Fig. 6I–K). As observed in HBE, m6A enrichment (IP/input) and spike-in normalized m6A signal over the lost peaks were significantly reduced in HNE, whereas input RNA was unchanged (Fig 6F; Supplemental Fig. S6H; Fig 6L; Supplemental Fig. S7C). Importantly, infected HNE showed a reduced level of nuclear METTL3 with partial relocalization to the cytoplasm (Supplemental Fig. S7D,E). We observed that infected HNEs were positive for spike protein and showed disturbed staining of cilia marker TUBB4A as reported earlier (Supplemental Fig. S7D; Robinot et al. 2021). Taken together, our observations in multiple transformed and primary cell culture models suggest that SARS-CoV-2 infection causes a global loss of m6A in cellular RNA.

    Previous SectionNext Section

    Discussion

    We report that SARS-CoV-2 infection causes a global loss of m6A in cellular RNAs while enabling m6A addition in viral RNAs, highlighting how this virus usurps mRNA modification pathways to promote its replication. Mechanistically, we found that SARS-CoV-2 infection induced partial relocalization of the m6A methyltransferase METTL3 from the nucleus to the cytoplasm, which compromised the formation of the METTL3/METTL14 complex in the nucleus and could account for the overall decrease in m6A modification in infected cells. The latter findings are consistent with those of Zhang et al. (2021b), who reported a cytoplasmic localization of both METTL3 and METTL14 along with other m6A reader proteins in SARS-CoV-2-infected Vero cells. The extent of METTL3 cytoplasmic relocalization and its effect on METTL3/METTL14 complex function and m6A levels depended on the viral variants. The B.1 and B.1.1.7 variants induced a more pronounced METTL3 relocalization and a more profound loss of m6A than the B.1.351 variant, suggesting that differences in variant fitness may depend not only on their capacity to escape from the innate and adaptive immune responses (Alefishat et al. 2022) but also on their capacity to perturb gene expression in infected cells. We also observed global loss of m6A peaks in primary human airway epithelial cells during SARS-CoV-2 infection, confirming the generality of our findings. In previous reports, however, a global loss in m6A was not detected in the Huh7 and A459 cell lines following SARS-CoV-2 infection, even though a decrease in m6A peaks was observed near the stop codon region in Huh7-infected cells (Burgess et al. 2021; Liu et al. 2021). The differential effect of SARS-CoV-2 infection on m6A regulation in Huh7 and A549 cells compared with Vero, HBE, and HNE infection models may result from different infection rates and/or different tissue origins of the cell types used (Saccon et al. 2021). Although we observed a global decrease in m6A peaks in infected cells, a fraction of transcripts nevertheless showed a gain in m6A peaks upon infection. These gained m6A modifications may be deposited by residual METTL3/METTL14 in the nucleus or by METTL16, which was recently described to also have m6A writing capacity in mRNAs (Su et al. 2022). Future investigations will be needed to understand how the variant-specific effect on m6A contributes to viral adaptation and how a subset of m6A peaks can be gained during infection.

    Gene expression analysis also highlighted common pathways dysregulated by all the SARS-CoV-2 variants tested. We observed in particular that genes associated with RNA catabolism tended to be up-regulated in infected Vero cells. Several genes associated with the m6A pathway, including m6A readers and m6A erasers, were deregulated as well during infection. Change in the RNA splicing pattern of m6A-related genes was recently reported after the acute depletion of METTL3, which was interpreted as a feedback loop to compensate for the acute depletion of m6A (Wei et al. 2021a). In SARS-CoV-2-infected cells, depletion of m6A following METTL3 relocalization and loss of functional METTL3/METTL14 complex in the nucleus probably activates a similar feedback loop, which could explain why genes involved in RNA catabolism and m6A modification were preferentially deregulated. Our combined analysis of differentially expressed genes and change in m6A peak abundance suggests that, in general, the loss of m6A peaks correlates with decreased RNA expression. However, it is important to consider that loss of m6A was more widespread than mRNA decrease in the infected cells and, thus, that loss of m6A did not always lead to decreased gene expression. Further investigation is needed to decipher the mechanisms that dictate how the loss of m6A influences gene expression during viral infection. Although we have pointed out DEU events in the subset of genes with m6A loss, further study is required to understand the effect of m6A loss on the other cellular genes that show no change in expression. We observed that m6A-containing genes were more prone to down-regulation during infection in both Vero and HBE infection models. The SARS-CoV-2 Nsp1 protein has been shown to inhibit the nuclear export of cellular RNAs during infection and to promote host mRNA decay (Burke et al. 2021; Zhang et al. 2021a). The nuclear export of cellular RNA is known to be regulated by m6A in a YTHDC1-dependent and nuclear transcription factor, X-box binding 1 (NFX1)–dependent manner (Roundtree et al. 2017). Viral protein Nsp1 also inhibits NFX1 function, resulting in the retention of cellular mRNAs in the nucleus during SARS-CoV-2 infection (Zhang et al. 2021a). It will be interesting to determine if the loss of m6A and Nsp1-mediated inhibition of NFX1 act synergistically to cause nuclear retention of cellular mRNAs during infection and if the retained mRNAs might be prone to degradation.

    Using a strand-specific m6A RIP-seq approach, we found that both the viral genomic RNA (positive-strand) and replicative negative-strand RNAs carried m6A modification. In the HBE model, we detected a robust m6A peak signal in genomic RNA at 4 and 7 dpi. The m6A peaks were also present in the replicative negative strand at 4 dpi, but not at day 7, presumably because viral replication had already abated at this time point. In a previous report by Liu et al. (2021), m6A peaks could be detected on negative-strand of SARS-CoV-2 using m6A RIP-seq, but not with the m6A cross-linking and immunoprecipitation (CLIP) technique. The investigators suggested lack of m6A detection could be owing to the limited coverage of the negative strand in CLIP data (Liu et al. 2021). Although our data from Vero and HBE infection models suggest the presence of m6A peaks in the negative strand, future experiments using m6A-CLIP in the human cell infection model will be required to precisely identify the specific residues modified by m6A. Viral RNA m6A modification is known to have a proviral effect by allowing escape from RIG-I binding and limiting the induction of inflammatory gene expression (Li et al. 2021). As RIG-I is known to be activated by viral double-stranded replicative intermediates (Yamada et al. 2021), it will be interesting to check if the m6A modifications detected in the replicative negative strand of SARS-CoV-2 contribute to an escape mechanism from RIG-I. Although the majority of m6A peaks were similar across the three SARS-CoV-2 variants studied, there were m6A peaks that were specific to particular variants. Further detailed studies will be required to understand if the differences in m6A modification across variants contribute to changes in pathogenicity by differential recruitment of m6A reader proteins or by modulating the RIG-I-dependent escape mechanisms described above.

    To explore new therapeutic options against SARS-CoV-2 infection, we have exploited the altered localization of the m6A writer METTL3. We succeeded in restoring METTL3 localization during infection by inhibiting XPO1 with selinexor. Selinexor treatment also increased the expression of specific genes that showed down-regulation with concomitant loss of m6A during infection. XPO1 inhibition by selinexor was previously proposed to be effective in limiting SARS-CoV-2 infection by promoting the nuclear relocalization of the ACE2 receptor (Kashyap et al. 2021). Our study provides a further mechanistic explanation for the effectiveness of XPO1 inhibition during SARS-CoV-2 infection. Although we provide evidence that altered METTL3 localization could be mediated by XPO1, we do not rule out other possible mechanisms that might affect METTL3 localization during SARS-CoV-2 infection. In particular, cytokines that are secreted by infected cells might influence METTL3 localization in both infected and bystander cells, resulting in broad perturbations in cellular m6A, a notion that requires further investigation. Cellular m6A-modified RNAs and m6A reader YTHDF proteins promote stress granule formation, which contributes to the antiviral defense (Fu and Zhuang 2020). The SARS-CoV-2 N-protein has been shown to phase-separate with the stress granule protein G3BP1 and thereby prevent stress granule formation in the infected cells (Wang et al. 2021). We propose an additional mechanism for the inhibition of stress granule formation through the depletion of m6A on cellular RNA. The recovery of stress granules formation in selinexor-treated infected cells may indeed result from restored m6A levels in cellular RNA, which could be visualized by the colocalization of G3BP1, m6A, and YTHDF2. We observed that selinexor treatment also had functional consequences, by restoring the gene expression and m6A enrichment during infection. The recovery of stress granules may have also contributed to stabilizing these RNAs, as association with stress granule proteins such as G3BP1 has been shown to regulate the stability of cellular RNAs (Laver et al. 2020; Somasekharan et al. 2020). Some of the validated m6A-dependent RNAs showed an association with G3BP1 and YTHDF2, suggesting that these RNAs could indeed be guided to stress granules in an m6A-dependent manner.

    Collectively, our findings highlight how SARS-CoV-2 perturbs the m6A RNA modification pathway to deregulate cellular RNAs and limit stress granule formation. Change in the cellular localization of METTL3 in infected cells resulted in a global loss of m6A in cellular mRNAs, whereas viral RNAs remained m6A modified. We propose that rescuing METTL3 localization during infection could be explored as a novel antiviral strategy against SARS-CoV-2.

    Previous SectionNext Section

    Methods

    Cell culture, SARS-CoV-2 infection, and XPO1 inhibitor treatment

    SARS-CoV-2 was isolated from samples with viral genomes sequenced at the diagnostic laboratory (Ringlander et al. 2021). Isolates from the B.1, B.1.1.7, B.1.351, and B.1.617.2 pango lineages were used to infect Vero CCL-81 cells (ATCC) grown in DMEM, 1% penicillin–streptomycin, and 2% fetal calf serum at 37°C and in 5% CO2. Details of the infection method in Vero cells, primary HBE cells (Lonza CC-2540S), and bronchial epithelial cell line BEAS-2B (ATCC) and XPO1 treatment are described in the Supplemental Material. SARS-CoV-2 infection of reconstructed human bronchial and nasal epithelia is also described in the Supplemental Material.

    Immunofluorescence staining

    Immunofluorescence staining was performed in Vero cells, primary HBE cells, and bronchial epithelial cell line BEAS-2B cells using the method as described in the Supplemental Material.

    HNE MucilAir cultures were fixed and cut using a scalpel blade, and immunofluorescence staining was performed as described in the Supplemental Material.

    RNA-seq and m6A RIP-seq

    RNA was isolated from both SARS-CoV-2-infected (48 h postinfection) and noninfected (mock) Vero cells using TRIzol reagent (Thermo Fisher Scientific 15596026) and direct-zol RNA miniprep (ZYMO Research R2050). For RNA extraction of primary epithelial cultures, cells were washed in cold PBS and then lyzed in 150 µL of TRIzol reagent (Thermo Fisher Scientific) added to the apical side of the insert. RNA was purified using the direct-zol miniprep kit (Zymo Research ZR2080), according to the manufacturer's specifications. HBE RNA was isolated at 4 and 7 dpi, whereas HNE RNA was isolated at 4 dpi along with the corresponding noninfected controls.

    RNA isolated either from Vero cells (15 μg) or from epithelial culture (5 μg) was supplemented with 10 ng or 3 ng of bacterial RNA, respectively, as a spike-in control, before fragmentation using RNA fragmentation reagents (Thermo Fisher Scientific AM8740). The fragmented RNA was used in RNA-seq and m6A RIP-seq experiments.

    Sequencing libraries for total RNA-seq/input for m6A RIP-seq were prepared from 10 ng of the fragmented RNA using SMARTer stranded total RNA-seq kit V2, pico input mammalian (Takara Bio). m6A RIP was performed with fragmented RNA as previously described (Zeng et al. 2018) with m6A antibody (Synaptic Systems 202003) incubated with 15 μg of Vero or 5 μg of epithelial (HBE and HNE) RNA. m6A RIP RNA was used to make sequencing libraries using the SMARTer stranded total RNA-seq kit V2, pico input mammalian (Takara Bio). All the libraries were single-end sequenced (1× 88 bp) on an Illumina NextSeq 2000 platform at the BEA core facility, Stockholm, Sweden.

    RT-qPCR and m6A RIP-qPCR

    The method followed for RT-qPCR and m6A RIP-qPCR are described in the Supplemental Material, and the primers used in this study are listed in Supplemental Table S1.

    Processing of RNA-seq data

    Single-end sequencing reads from SMARTer stranded total RNA-seq kit v2 were analyzed with FastQC for quality control (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/), and adapters were removed using Trim Galore! v0.6.6 (https://github.com/FelixKrueger/TrimGalore) with a minimal length threshold of 20 bp. Trimmed reads from Vero infected and noninfected cells were aligned to the following reference genomes: concatenated chlSab2 (Chlorosebus sabeus), wuhCor1 (SARS-CoV-2), plus Escherichia coli str. K-12 substr. MG1655 (spike-in), or the SARS-CoV-2 + E. coli genomes, obtained from the UCSC Genome Browser (http://hgdownload.soe.ucsc.edu/goldenPath/). Alignments were made using HISAT2 v2.2.1 (Kim et al. 2019) with the parameters (-U ‐‐rna-strandness R). Sequencing data from human bronchial and nasal epithelial noninfected and infected cells were mapped to the concatenated human (hg38) + SARS-CoV-2 (wuhCor1) + E. coli (spike-in) reference genomes. Duplicate alignments were labeled using MarkDuplicates from Picard v2.23.4 (http://broadinstitute.github.io/picard), and marked alignment files were further processed using Sambamba v0.7.1 (Tarasov et al. 2015; Kim et al. 2019), keeping only mapping reads separated by strand after duplicate removal.

    Analysis of RNA-seq data for differential gene expression

    Differential gene expression following SARS-CoV-2 infection and Gene Ontology (GO) analysis is described in detail in the Supplemental Material.

    m6A RIP-seq data analysis

    m6A peak calling on SARS-CoV-2 viral genome and infected cell samples was performed with callpeak from MACS2 v2.2.6 (Zhang et al. 2008) on IP and input processed alignments mapped to SARS-CoV-2 wuhCor1 or the concatenated genomes (chlSab2 + wuhCor1) or (hg38 + wuCor1), using the parameters “‐‐no-model ‐‐keep-dup auto ‐‐call-summits” and effective genome sizes according to each genome with a P-value cutoff of 0.05 for SARS-CoV-2 genome and a P-value cutoff of 0.01 for the Vero and human data sets. Samples processed using MACS2 were scaled by default to the library size to account for differences during peak calling.

    The called peaks were annotated according to the nearest genomic feature using annotatePeaks.pl from HOMER v4.11 (Heinz et al. 2010). Retained, gained, and lost peaks between noninfected and infected cells were identified using BEDTools intersect from BEDTools v2.29.2 (Quinlan and Hall 2010). Briefly, peaks were classified as “retained” depending on whether a peak was intersected between noninfected and infected samples; “lost” if a peak was present in noninfected cells but not intersecting with peaks of infected samples; and “gained” if a peak was present in infected cells but not in noninfected cells. Motif analysis on m6A called peaks was performed with findMotifsGenome.pl from HOMER and visualized with the universalmotif R package (https://bioconductor.org/packages/universalmotif/) using R computing language (R core team 2022). Changes in overall m6A enrichment across transcripts after infection were calculated by extracting the spike-in normalized m6A/input signals from bigWig files within the annotated gene coordinates in the chlsab2 genomes using the rtracklayer package (Lawrence et al. 2009). Further correlation between m6A enrichment and gene expression changes (log2 fold-changes infected vs. noninfected) were then generated. Potential N6,2′-O-dimethyladenosine (m6Am) signals were removed from m6A RIP-seq data using the method as described in the Supplemental Material.

    Normalization of m6A RIP-seq data using spike-ins

    To control for systematic variations across m6A RIP experiments, the amount of spike-in bacterial RNA was estimated by counting the total number of reads uniquely mapped to the E. coli K-12 reference genome using Sambamba v0.7.1 (Tarasov et al. 2015). E. coli spike-in counts were further used to calculate scaling factors for each batch of m6A RIP-seq samples (Supplemental Table S2). Computed scaling factors were then used in metagene density distribution analyses (described below) to normalize the density of m6A peaks to the spike-in content in each sample. Genome-wide coverage tracks were calculated using normalized counts per million (CPM) and further adjusted according to the computed spike-in scaling factors for m6A and input RNA signals using the bamCoverage –scaleFactor parameter from deepTools v3.3.2 (Ramírez et al. 2016). m6A RIP/input ratio tracks were then calculated using the bigwigCompare function from deepTools.

    Metagene analysis

    To analyze the genome-wide distribution of m6A, a metagene analysis of m6A peak density distribution was performed by overlapping the peak coordinates with the following genomic features: 5′ UTR, CDS, and 3′ UTR obtained from the GTF genome annotation files from UCSC (for chlsab2) or GENCODE v36 (for hg38), considering the longest isoform for each gene. Each transcript was scaled to fixed-size metagene bins according to each reference genome. m6A peak density distribution profiles were generated after mapping the m6A peaks to the metagene coordinates using the plyranges R package (Lee et al. 2019). To compare multiple conditions, the relative m6A density distributions were calculated using the relative density function from the ggmulti package (https://cran.r-project.org/web/packages/ggmulti/index.html). The relative density function calculates the sum of the density estimate area of all conditions, where the total sum is scaled to a maximum of one and the area of each condition is proportional to its own count. The relative m6A-RIP peak distributions were further normalized using the calculated spike-in factors.

    Enrichment of m6A peaks

    To calculate the enrichment of RIP signals across m6A peak regions, we extracted the number of reads in m6A-RIP (m6A signal) and input RNA alignments 125 bp around peak summit coordinates using the ScoreMatrix function from the genomation package using the parameter rpm = TRUE to obtain normalized CPM counts to account for differences of library sizes (Akalin et al. 2015). The m6A-RIP and input RNA signals around 125-bp peak summits were further normalized using spike-in factors. The log2 ratio of m6A signal/input RNA was then calculated using the enrichmentMatrix function from the same package.

    Previous SectionNext Section

    Data access

    All raw and processed sequencing data generated in this study have been submitted to the NCBI Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) under accession number GSE188477. The source code used for the analysis of data in this study is available at GitHub (https://github.com/AkramMendez/m6a_sarscov2) and as Supplemental Code.

    Previous SectionNext Section

    Competing interest statement

    The authors declare no competing interests.

    Previous SectionNext Section

    Acknowledgments

    We thank the core facility at Novum, BEA, Bioinformatics and Expression Analysis, which is supported by the board of research at the Karolinska Institute and the research committee at the Karolinska hospital for help with sequencing, and the Proteomic core facility at GU for their support with LC-MS/MS. The computations and data handling were enabled by resources in project SNIC-2022-22-85 provided by the Swedish National Infrastructure for Computing (SNIC) at UPPMAX, partially funded by the Swedish Council through grant agreement no. 2018-05973. We acknowledge Anders Sjölander for assistance concerning technical and implementational aspects of the UPPMAX resources. This work was funded by grants from the Swedish Research Council (Vetenskapsrådet) to T.M. (2018-02224), Sweden-South Korea COVID-19 grant from the Swedish Research Council (Vetenskapsrådet) to T.M. (2020-06311 and 2022-05965), and project grant to T.M. from Svenska Läkaresällskapet; Kungl. Vetenskaps- och Vitterhets-Samhället (KVVS) grant to T.M.; grants COROCHIP and PFR7 from the Urgence COVID-19 Fundraising Campaign of Institut Pasteur to L.A.C.; Bollan scholarship to R.V.; and grants from the National Research Foundation of Korea to J.-J.D. (NRF-2021K1A4A7A02098793) and J.Y.C. (NRF-2022K1A3A1A47088296).

    Previous SectionNext Section

    Footnotes

    • [Supplemental material is available for this article.]

    • Article published online before print. Article, supplemental material, and publication date are at https://www.genome.org/cgi/doi/10.1101/gr.276407.121.

    • Freely available online through the Genome Research Open Access option.

    • Received November 20, 2021.
    • Accepted February 17, 2023.

    This article, published in Genome Research, is available under a Creative Commons License (Attribution-NonCommercial 4.0 International), as described at http://creativecommons.org/licenses/by-nc/4.0/.

    Previous Section
     

    References

    1. ↵
      Akalin A, Franke V, Vlahoviček K, Mason CE, Schübeler D. 2015. genomation: a toolkit to summarize, annotate and visualize genomic intervals. Bioinformatics 31: 1127–1129. doi:10.1093/bioinformatics/btu775
      CrossRefMedlineGoogle Scholar
    2. ↵
      Alarcón CR, Goodarzi H, Lee H, Liu X, Tavazoie S, Tavazoie SF. 2015. HNRNPA2B1 is a mediator of m6A-dependent nuclear RNA processing events. Cell 162: 1299–1308. doi:10.1016/j.cell.2015.08.011
      CrossRefMedlineGoogle Scholar
    3. ↵
      Alefishat E, Jelinek HF, Mousa M, Tay GK, Alsafar HS. 2022. Immune response to SARS-CoV-2 variants: a focus on severity, susceptibility, and preexisting immunity. J Infect Public Health 15: 277–288. doi:10.1016/j.jiph.2022.01.007
      CrossRefGoogle Scholar
    4. ↵
      Baquero-Perez B, Geers D, Díez J. 2021. From A to m6A: the emerging viral epitranscriptome. Viruses 13: 1049. doi:10.3390/v13061049
      CrossRefMedlineGoogle Scholar
    5. ↵
      Bawankar P, Lence T, Paolantoni C, Haussmann IU, Kazlauskiene M, Jacob D, Heidelberger JB, Richter FM, Nallasivan MP, Morin V, et al. 2021. Hakai is required for stabilization of core components of the m6A mRNA methylation machinery. Nat Commun 12: 3778. doi:10.1038/s41467-021-23892-5
      CrossRefGoogle Scholar
    6. ↵
      Blanco-Melo D, Nilsson-Payant BE, Liu W-C, Uhl S, Hoagland D, Møller R, Jordan TX, Oishi K, Panis M, Sachs D, et al. 2020. Imbalanced host response to SARS-CoV-2 drives development of COVID-19. Cell 181: 1036–1045.e9. doi:10.1016/j.cell.2020.04.026
      CrossRefMedlineGoogle Scholar
    7. ↵
      Burgess HM, Depledge DP, Thompson L, Srinivas KP, Grande RC, Vink EI, Abebe JS, Blackaby Wesley P, Hendrick A, Albertella MR, et al. 2021. Targeting the m6A RNA modification pathway blocks SARS-CoV-2 and HCoV-OC43 replication. Genes Dev 35: 1005–1019. doi:10.1101/gad.348320.121
      Abstract/FREE Full Text
    8. ↵
      Burke JM, St Clair LA, Perera R, Parker R. 2021. SARS-CoV-2 infection triggers widespread host mRNA decay leading to an mRNA export block. RNA 27: 1318–1329. doi:10.1261/rna.078923.121
      Abstract/FREE Full Text
    9. ↵
      Dossin F, Pinheiro I, Żylicz JJ, Roensch J, Collombet S, Le Saux A, Chelmicki T, Attia M, Kapoor V, Zhan Y, et al. 2020. SPEN integrates transcriptional and epigenetic control of X-inactivation. Nature 578: 455–460. doi:10.1038/s41586-020-1974-9
      CrossRefMedlineGoogle Scholar
    10. ↵
      Fu Y, Zhuang X. 2020. m6A-binding YTHDF proteins promote stress granule formation. Nat Chem Biol 16: 955–963. doi:10.1038/s41589-020-0524-y
      CrossRefGoogle Scholar
    11. ↵
      Gao X, Liu Y, Zou S, Liu P, Zhao J, Yang C, Liang M, Yang J. 2021. Genome-wide screening of SARS-CoV-2 infection-related genes based on the blood leukocytes sequencing data set of patients with COVID-19. J Med Virol 93: 5544–5554. doi:10.1002/jmv.27093
      CrossRefGoogle Scholar
    12. ↵
      Geula S, Moshitch-Moshkovitz S, Dominissini D, Mansour AA, Kol N, Salmon-Divon M, Hershkovitz V, Peer E, Mor N, Manor YS, et al. 2015. Stem cells. m6A mRNA methylation facilitates resolution of naïve pluripotency toward differentiation. Science 347: 1002–1006. doi:10.1126/science.1261417
      Abstract/FREE Full Text
    13. ↵
      Heinz S, Benner C, Spann N, Bertolino E, Lin YC, Laslo P, Cheng JX, Murre C, Singh H, Glass CK. 2010. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol Cell 38: 576–589. doi:10.1016/j.molcel.2010.05.004
      CrossRefMedlineGoogle Scholar
    14. ↵
      Huang H, Weng H, Zhou K, Wu T, Zhao BS, Sun M, Chen Z, Deng X, Xiao G, Auer F, et al. 2019. Histone H3 trimethylation at lysine 36 guides m6A RNA modification co-transcriptionally. Nature 567: 414–419. doi:10.1038/s41586-019-1016-7
      CrossRefGoogle Scholar
    15. ↵
      Kashyap T, Murray J, Walker Christopher J, Chang H, Tamir S, Hou B, Shacham S, Kauffman MG, Tripp RA, Landesman Y. 2021. Selinexor, a novel selective inhibitor of nuclear export, reduces SARS-CoV-2 infection and protects the respiratory system in vivo. Antiviral Res 192: 105115. doi:10.1016/j.antiviral.2021.105115
      CrossRefGoogle Scholar
    16. ↵
      Kim G-W, Siddiqui A. 2021. Hepatitis B virus X protein recruits methyltransferases to affect cotranscriptional N6-methyladenosine modification of viral/host RNAs. Proc Natl Acad Sci 118: e2019455118. doi:10.1073/pnas.2019455118
      Abstract/FREE Full Text
    17. ↵
      Kim D, Paggi JM, Park C, Bennett C, Salzberg SL. 2019. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat Biotechnol 37: 907–915. doi:10.1038/s41587-019-0201-4
      CrossRefMedlineGoogle Scholar
    18. ↵
      Kim D, Lee J-Y, Yang J-S, Kim JW, Kim VN, Chang H. 2020. The architecture of SARS-CoV-2 transcriptome. Cell 181: 914–921.e10. doi:10.1016/j.cell.2020.04.011
      CrossRefMedlineGoogle Scholar
    19. ↵
      Laver JD, Ly J, Winn AK, Karaiskakis A, Lin S, Nie K, Benic G, Jaberi-Lashkari N, Cao WX, Khademi A, et al. 2020. The RNA-binding protein Rasputin/G3BP enhances the stability and translation of its target mRNAs. Cell Rep 30: 3353–3367.e7. doi:10.1016/j.celrep.2020.02.066
      CrossRefGoogle Scholar
    20. ↵
      Lawrence M, Gentleman R, Carey V. 2009. rtracklayer: an R package for interfacing with genome browsers. Bioinformatics 25: 1841–1842. doi:10.1093/bioinformatics/btp328
      CrossRefMedlineGoogle Scholar
    21. ↵
      Lee S, Cook D, Lawrence M. 2019. plyranges: a grammar of genomic data transformation. Genome Biol 20: 4. doi:10.1186/s13059-018-1597-8
      CrossRefGoogle Scholar
    22. ↵
      Li N, Hui H, Bray B, Gonzalez GM, Zeller M, Anderson KG, Knight R, Smith D, Wang Y, Carlin AF, et al. 2021. METTL3 regulates viral m6A RNA modification and host cell innate immune responses during SARS-CoV-2 infection. Cell Rep 35: 109091. doi:10.1016/j.celrep.2021.109091
      CrossRefGoogle Scholar
    23. ↵
      Lichinchi G, Zhao BS, Wu Y, Lu Z, Qin Y, He C, Rana TM. 2016. Dynamics of human and viral RNA methylation during Zika virus infection. Cell Host Microbe 20: 666–673. doi:10.1016/j.chom.2016.10.002
      CrossRefMedlineGoogle Scholar
    24. ↵
      Liu J, Xu Y-P, Li K, Ye Q, Zhou H-Y, Sun H, Li X, Yu L, Deng Y-Q, Li R-T, et al. 2021. The m6A methylome of SARS-CoV-2 in host cells. Cell Res 31: 404–414. doi:10.1038/s41422-020-00465-7
      CrossRefGoogle Scholar
    25. ↵
      Lorè NI, De Lorenzo R, Rancoita PMV, Cugnata F, Agresti A, Benedetti F, Bianchi ME, Bonini C, Capobianco A, Conte C, et al. 2021. CXCL10 levels at hospital admission predict COVID-19 outcome: hierarchical assessment of 53 putative inflammatory biomarkers in an observational study. Mol Med 27: 129. doi:10.1186/s10020-021-00390-4
      CrossRefMedlineGoogle Scholar
    26. ↵
      Lu W, Tirumuru N, St. Gelais C, Koneru PC, Liu C, Kvaratskhelia M, He C, Wu L. 2018. N6-Methyladenosine–binding proteins suppress HIV-1 infectivity and viral production. J Biol Chem 293: 12992–13005. doi:10.1074/jbc.RA118.004215
      Abstract/FREE Full Text
    27. ↵
      Meyer KD, Saletore Y, Zumbo P, Elemento O, Mason CE, Jaffrey SR. 2012. Comprehensive analysis of mRNA methylation reveals enrichment in 3′ UTRs and near stop codons. Cell 149: 1635–1646. doi:10.1016/j.cell.2012.05.003
      CrossRefMedlineGoogle Scholar
    28. ↵
      Patil DP, Chen C-K, Pickering BF, Chow A, Jackson C, Guttman M, Jaffrey SR. 2016. m6A RNA methylation promotes XIST-mediated transcriptional repression. Nature 537: 369–373. doi:10.1038/nature19342
      CrossRefMedlineGoogle Scholar
    29. ↵
      Ping X-L, Sun B-F, Wang L, Xiao W, Yang X, Wang W-J, Adhikari S, Shi Y, Lv Y, Chen Y-S, et al. 2014. Mammalian WTAP is a regulatory subunit of the RNA N6-methyladenosine methyltransferase. Cell Res 24: 177–189. doi:10.1038/cr.2014.3
      CrossRefMedlineGoogle Scholar
    30. ↵
      Quinlan AR, Hall IM. 2010. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26: 841–842. doi:10.1093/bioinformatics/btq033
      CrossRefMedlineGoogle Scholar
    31. ↵
      Ramírez F, Ryan DP, Grüning B, Bhardwaj V, Kilpert F, Richter Andreas S, Heyne S, Dündar F, Manke T. 2016. deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Res 44: W160–W165. doi:10.1093/nar/gkw257
      CrossRefMedlineGoogle Scholar
    32. ↵
      Ran J, Yang Y, Li D, Liu M, Zhou J. 2015. Deacetylation of α-tubulin and cortactin is required for HDAC6 to trigger ciliary disassembly. Sci Rep 5: 12917. doi:10.1038/srep12917
      CrossRefMedlineGoogle Scholar
    33. ↵
      R Core Team. 2022. R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna. https://www.R-project.org/.
      Google Scholar
    34. ↵
      Ringlander J, Olausson J, Nyström K, Härnqvist T, Jakobsson HE, Lindh M. 2021. Recurrent and persistent infection with SARS-CoV-2: epidemiological data and case reports from Western Sweden, 2020. Infect Dis 53: 900–907. doi:10.1080/23744235.2021.1957143
      CrossRefGoogle Scholar
    35. ↵
      Riva L, Yuan S, Yin X, Martin-Sancho L, Matsunaga N, Pache L, Burgstaller-Muehlbacher S, De Jesus PD, Teriete P, Hull MV, et al. 2020. Discovery of SARS-CoV-2 antiviral drugs through large-scale compound repurposing. Nature 586: 113–119. doi:10.1038/s41586-020-2577-1
      CrossRefMedlineGoogle Scholar
    36. ↵
      Robinot R, Hubert M, de Melo GD, Lazarini F, Bruel T, Smith N, Levallois S, Larrous F, Fernandes J, Gellenoncourt S, et al. 2021. SARS-CoV-2 infection induces the dedifferentiation of multiciliated cells and impairs mucociliary clearance. Nat Commun 12: 4354. doi:10.1038/s41467-021-24521-x
      CrossRefGoogle Scholar
    37. ↵
      Roundtree IA, Luo G-Z, Zhang Z, Wang X, Zhou T, Cui Y, Sha J, Huang X, Guerrero L, Xie P, et al. 2017. YTHDC1 mediates nuclear export of N6-methyladenosine methylated mRNAs. eLife 6: e31311. doi:10.7554/eLife.31311
      CrossRefMedlineGoogle Scholar
    38. ↵
      Saccon E, Chen X, Mikaeloff F, Rodriguez JE, Szekely L, Vinhas BS, Krishnan S, Byrareddy SN, Frisan T, Végvári Á, et al. 2021. Cell-type-resolved quantitative proteomics map of interferon response against SARS-CoV-2. iScience 24: 102420. doi:10.1016/j.isci.2021.102420
      CrossRefGoogle Scholar
    39. ↵
      Samelson AJ, Tran QD, Robinot R, Carrau L, Rezelj VV, Kain AM, Chen M, Ramadoss GN, Guo X, Lim SA, et al. 2022. BRD2 inhibition blocks SARS-CoV-2 infection by reducing transcription of the host cell receptor ACE2. Nat Cell Biol 24: 24–34. doi:10.1038/s41556-021-00821-8
      CrossRefGoogle Scholar
    40. ↵
      Schmidt N, Lareau CA, Keshishian H, Ganskih S, Schneider C, Hennig T, Melanson R, Werner S, Wei Y, Zimmer M, et al. 2021. The SARS-CoV-2 RNA–protein interactome in infected human cells. Nat Microbiol 6: 339–353. doi:10.1038/s41564-020-00846-z
      CrossRefMedlineGoogle Scholar
    41. ↵
      Somasekharan SP, Zhang F, Saxena N, Huang JN, Kuo I-C, Low C, Robert B, Adomat H, Stoynov N, Foster L, et al. 2020. G3BP1-linked mRNA partitioning supports selective protein synthesis in response to oxidative stress. Nucleic Acids Res 48: 6855–6873. doi:10.1093/nar/gkaa376
      CrossRefGoogle Scholar
    42. ↵
      Su R, Dong L, Li Y, Gao M, He PC, Liu W, Wei J, Zhao Z, Lei G, Han L, et al. 2022. METTL16 exerts an m6A-independent function to facilitate translation and tumorigenesis. Nat Cell Biol 24: 205–216. doi:10.1038/s41556-021-00835-2
      CrossRefGoogle Scholar
    43. ↵
      Sun H, Li K, Zhang X, Liu J, Zhang M, Meng H, Yi C. 2021. m6Am-seq reveals the dynamic m6Am methylation in the human transcriptome. Nat Commun 12: 4778. doi:10.1038/s41467-021-25105-5
      CrossRefGoogle Scholar
    44. ↵
      Tao K, Tzou PL, Nouhin J, Gupta RK, de Oliveira T, Kosakovsky Pond SL, Fera D, Shafer RW. 2021. The biological and clinical significance of emerging SARS-CoV-2 variants. Nat Rev Genet 22: 757–773. doi:10.1038/s41576-021-00408-x
      CrossRefMedlineGoogle Scholar
    45. ↵
      Tarasov A, Vilella AJ, Cuppen E, Nijman IJ, Prins P. 2015. Sambamba: fast processing of NGS alignment formats. Bioinformatics 31: 2032–2034. doi:10.1093/bioinformatics/btv098
      CrossRefMedlineGoogle Scholar
    46. ↵
      Wang J, Shi C, Xu Q, Yin H. 2021. SARS-CoV-2 nucleocapsid protein undergoes liquid-liquid phase separation into stress granules through its N-terminal intrinsically disordered region. Cell Discov 7: 5. doi:10.1038/s41421-020-00240-3
      CrossRefGoogle Scholar
    47. ↵
      Wei G, Almeida M, Pintacuda G, Coker H, Bowness JS, Ule J, Brockdorff N. 2021a. Acute depletion of METTL3 implicates N6-methyladenosine in alternative intron/exon inclusion in the nascent transcriptome. Genome Res 31: 1395–1408. doi:10.1101/gr.271635.120
      Abstract/FREE Full Text
    48. ↵
      Wei J, Alfajaro MM, DeWeirdt PC, Hanna RE, Lu-Culligan WJ, Cai WL, Strine MS, Zhang S-M, Graziano VR, Schmitz CO, et al. 2021b. Genome-wide CRISPR screens reveal host factors critical for SARS-CoV-2 infection. Cell 184: 76–91.e13. doi:10.1016/j.cell.2020.10.028
      CrossRefMedlineGoogle Scholar
    49. ↵
      Wei J, Yu X, Yang L, Liu X, Gao B, Huang B, Dou X, Liu J, Zou Z, Cui X-L, et al. 2022. FTO mediates LINE1 m6A demethylation and chromatin regulation in mESCs and mouse development. Science 376: 968–973. doi:10.1126/science.abe9582
      CrossRefGoogle Scholar
    50. ↵
      Wu F, Zhao S, Yu B, Chen Y-M, Wang W, Song Z-G, Hu Y, Tao Z-W, Tian J-H, Pei Y-Y, et al. 2020. A new coronavirus associated with human respiratory disease in China. Nature 579: 265–269. doi:10.1038/s41586-020-2008-3
      CrossRefMedlineGoogle Scholar
    51. ↵
      Xiao W, Adhikari S, Dahal U, Chen Y-S, Hao Y-J, Sun B-F, Sun H-Y, Li A, Ping X-L, Lai W-Y, et al. 2016. Nuclear m6A reader YTHDC1 regulates mRNA splicing. Mol Cell 61: 507–519. doi:10.1016/j.molcel.2016.01.012
      CrossRefMedlineGoogle Scholar
    52. ↵
      Yamada T, Sato S, Sotoyama Y, Orba Y, Sawa H, Yamauchi H, Sasaki M, Takaoka A. 2021. RIG-I triggers a signaling-abortive anti-SARS-CoV-2 defense in human lung cells. Nat Immunol 22: 820–828. doi:10.1038/s41590-021-00942-0
      CrossRefGoogle Scholar
    53. ↵
      Yankova E, Blackaby W, Albertella M, Rak J, de Braekeleer E, Tsagkogeorga G, Pilka ES, Aspris D, Leggate D, Hendrick AG, et al. 2021. Small-molecule inhibition of METTL3 as a strategy against myeloid leukaemia. Nature 593: 597–601. doi:10.1038/s41586-021-03536-w
      CrossRefGoogle Scholar
    54. ↵
      Yue Y, Liu J, Cui X, Cao J, Luo G, Zhang Z, Cheng T, Gao M, Shu X, Ma H, et al. 2018. VIRMA mediates preferential m6A mRNA methylation in 3′UTR and near stop codon and associates with alternative polyadenylation. Cell Discov 4: 10. doi:10.1038/s41421-018-0019-0
      CrossRefGoogle Scholar
    55. ↵
      Zaccara S, Ries RJ, Jaffrey SR. 2019. Reading, writing and erasing mRNA methylation. Nat Rev Mol Cell Biol 20: 608–624. doi:10.1038/s41580-019-0168-5
      CrossRefMedlineGoogle Scholar
    56. ↵
      Zeng Y, Wang S, Gao S, Soares F, Ahmed M, Guo H, Wang M, Hua JT, Guan J, Moran MF, et al. 2018. Refined RIP-seq protocol for epitranscriptome analysis with low input materials. PLoS Biol 16: e2006092. doi:10.1371/journal.pbio.2006092
      CrossRefGoogle Scholar
    57. ↵
      Zhang M-Z, Mai W, Li C, Cho S-Y, Hao C, Moeckel G, Zhao R, Kim I, Wang J, Xiong H, et al. 2004. PKHD1 protein encoded by the gene for autosomal recessive polycystic kidney disease associates with basal bodies and primary cilia in renal epithelial cells. Proc Natl Acad Sci 101: 2311–2316. doi:10.1073/pnas.0400073101
      Abstract/FREE Full Text
    58. ↵
      Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, Nusbaum C, Myers RM, Brown M, Li W, et al. 2008. Model-based Analysis of ChIP-seq (MACS). Genome Biol 9: R137. doi:10.1186/gb-2008-9-9-r137
      CrossRefMedlineGoogle Scholar
    59. ↵
      Zhang S, Zhao BS, Zhou A, Lin K, Zheng S, Lu Z, Chen Y, Sulman EP, Xie K, Bögler O, et al. 2017. m6A demethylase ALKBH5 maintains tumorigenicity of glioblastoma stem-like cells by sustaining FOXM1 expression and cell proliferation program. Cancer Cell 31: 591–606.e6. doi:10.1016/j.ccell.2017.02.013
      CrossRefMedlineGoogle Scholar
    60. ↵
      Zhang K, Miorin L, Makio T, Dehghan I, Gao S, Xie Y, Zhong H, Esparza M, Kehrer T, Kumar A, et al. 2021a. Nsp1 protein of SARS-CoV-2 disrupts the mRNA export machinery to inhibit host gene expression. Sci Adv 7: eabe7386. doi:10.1126/sciadv.abe7386
      FREE Full Text
    61. ↵
      Zhang X, Hao H, Ma L, Zhang Y, Hu X, Chen Z, Liu D, Yuan J, Hu Z, Guan W. 2021b. Methyltransferase-like 3 modulates severe acute respiratory syndrome coronavirus-2 RNA N6-methyladenosine modification and replication. mBio 12: e0106721. doi:10.1128/mBio.01067-21
      CrossRefGoogle Scholar
    62. ↵
      Zheng Z-Q, Wang S-Y, Xu Z-S, Yu-Zhi F, Wang Y-Y. 2021. SARS-CoV-2 nucleocapsid protein impairs stress granule formation to promote viral replication. Cell Discov 7: 38. doi:10.1038/s41421-021-00275-0
      CrossRefGoogle Scholar
    | Table of Contents
    OPEN ACCESS ARTICLE

    This Article

    1. Published in Advance March 1, 2023, doi: 10.1101/gr.276407.121 Genome Res. 33: 299-313 © 2023 Vaid et al.; Published by Cold Spring Harbor Laboratory Press

    Article Category

    ORCID

    1. Profile for Vaid, R.http://orcid.org/0000-0002-2074-5080
    2. Profile for Mendez, A.http://orcid.org/0000-0001-9195-3808
    3. Profile for Thombare, K.http://orcid.org/0000-0002-0183-3376
    4. Profile for Burgos-Panadero, R.http://orcid.org/0000-0003-0811-2675
    5. Profile for Robinot, R.http://orcid.org/0000-0002-3651-0171
    6. Profile for Reinius, B.http://orcid.org/0000-0002-7021-5248
    7. Profile for Chakrabarti, L. A.http://orcid.org/0000-0002-1895-3630
    8. Profile for Mondal, T.http://orcid.org/0000-0001-6228-8001

    Share

    Preprint Server



    Navigate This Article

    Current Issue

    1. March 2026, 36 (3)
    1. Current Issue
    1. Alert me to new issues of Genome Research