The m6A landscape of polyadenylated nuclear (PAN) RNA and its related methylome in the context of KSHV replication.

Polyadenylated nuclear (PAN) RNA is a long noncoding transcript involved in Kaposi's sarcoma-associated herpesvirus (KSHV) lytic reactivation and regulation of cellular and viral gene expression. We have previously shown that PAN RNA has dynamic secondary structure and protein binding profiles that can be influenced by epitranscriptomic modifications. N6-methyladenosine (m6A) is one of the most abundant chemical signatures found in viral RNA genomes and virus-encoded RNAs. Here, we combined antibody-independent next-generation mapping with direct RNA sequencing to address the epitranscriptomic status of PAN RNA in KSHV infected cells. We showed that PAN m6A status is dynamic, reaching the highest number of modifications at the late lytic stages of KSHV infection. Using a newly developed method, termed selenium-modified deoxythymidine triphosphate (SedTTP)-reverse transcription (RT) and ligation assisted PCR analysis of m6A (SLAP), we gained insight into the fraction of modification at identified sites. By applying comprehensive proteomic approaches, we identified writers and erasers that regulate the m6A status of PAN, and readers that can convey PAN m6A phenotypic effects. We verified the temporal and spatial subcellular availability of the methylome components for PAN modification by performing confocal microscopy analysis. Additionally, the RNA biochemical probing (SHAPE-MaP) outlined local and global structural alterations invoked by m6A in the context of full-length PAN RNA. This work represents the first comprehensive overview of the dynamic interplay that takes place between the cellular epitranscriptomic machinery and a specific viral RNA in the context of KSHV infected cells.

INTRODUCTION

Kaposi's sarcoma-associated herpesvirus (KSHV) is etiologically linked to several malignancies, including Kaposi's sarcoma (KS), primary effusion lymphoma (PEL), and multicentric Castleman's disease (Carbone et al. 2015; Lurain et al. 2018; Ramaswami et al. 2021). Like all herpesviruses, the KSHV infectivity cycle is bipartite, consisting of latent and lytic stages of infection. During latency, KSHV persists as episomal DNA and expresses only a few of its genes. The full repertoire of viral genes expression occurs in a highly ordered manner during the KSHV lytic replication (Arias et al. 2014). Polyadenylated nuclear (PAN) RNA is a multifunctional long noncoding (lnc) RNA produced at high levels during KSHV lytic reactivation (Rossetto and Pari 2012), although it is also expressed during KSHV latency in cultured cells (Rossetto et al. 2013). PAN RNA regulates cellular and viral gene expression by recruiting chromatin-modifying complexes to cellular and viral DNA (Rossetto and Pari 2014; Purushothaman et al. 2016; Campbell and Izumiya 2020). Another study showed that PAN facilitates late viral mRNA export to the cytoplasm, increasing the levels of viral mRNAs available for the expression of late lytic products (Schifano et al. 2017). PAN also acts as a scaffold for viral proteins ORF50 and ORF57, and RNA polymerase II at chromatin, leading to the execution of the viral lytic program (Campbell and Izumiya 2020). We have previously shown that PAN RNA has dynamic secondary structure and protein binding profiles, which change depending upon biological context (Sztuba-Solinska et al. 2017). This work constitutes the most extensive structural characterization of viral lncRNA and its interactome inside the living cells and virions, providing a broad framework for understanding PAN's roles in KSHV infection.

Recent studies have shown that the structure (Roost et al. 2015; Jacob et al. 2017; Liu et al. 2017; Coker et al. 2019; Waduge et al. 2019; Charette and Gray 2000), metabolism (Helm and Alfonzo 2014; Dai et al. 2018; Hao et al. 2019), and function (Charette and Gray 2000; Schwartz et al. 2014a; Trixl and Lusser 2019; Porman et al. 2020) of coding and noncoding RNAs can be modulated by over 100 different types of chemical modifications, which are collectively referred to as the epitranscriptome. N6-methyladenosine (m6A) is the most abundant of those, as it has been found in over 25% of transcripts expressed in mammalian cells (Ye 2018), as well as in viral RNA genomes (Gokhale et al. 2016; Lichinchi et al. 2016a,b; Tirumuru et al. 2016; McIntyre et al. 2018) and virus-encoded RNAs (Kennedy et al. 2016; Baquero-Perez et al. 2021). m6A has been shown to impact nearly every aspect of RNA biology, from splicing to translation, and all the way through stability and decay (Zhou et al. 2018, 2019; Lee et al. 2020). The breadth of m6A impact has been attributed to creation of new sites for protein recognition, in part via local structural changes to the modified RNA (Liu et al. 2015, 2017; Ashraf et al. 2019). However, the biological significance of m6A, especially in the context of a specific viral transcript, has not been defined.

Previous studies have outlined the mechanism regulating the m6A status of KSHV coding RNAs (Ye et al. 2017; Hesser et al. 2018; Ye 2018; Baquero-Perez et al. 2019; Barker et al. 2020). It has been shown that the modification is mediated by a complex of methyltransferases (METTL3/14) and Wilms tumor 1 associated protein (WTAP), which are recruited to target mRNAs via RNA-binding motif 15 protein (RBM15) (Knuckles et al. 2018). The fat mass and obesity-associated protein (FTO) reverses this process and acts as an m6A eraser (Jia et al. 2011; Bartosovic et al. 2017). Additionally, the biological functions of m6A can be mediated by reader proteins (Hesser et al. 2018). For example, heterogeneous nuclear ribonucleoprotein C (HNRNPC), responsible for RNA splicing and stability, selectively binds m6A, where the modification acts as a “switch,” destabilizing an RNA helix and exposing a single-stranded HNRNPC binding motif (Liu et al. 2015). The structural influence of this “m6A-switch” also dictates the binding of YTH domain-containing proteins. In the context of KSHV infected cells, the YTHDC1 binding potentiates viral mRNA splicing (Ye et al. 2017), while YTHDF2 overexpression facilitates the degradation of KSHV coding transcripts (Tan et al. 2017). In addition, Staphylococcal nuclease domain-containing protein 1 (SND1) has been reported to bind and stabilize the KSHV ORF50 mRNA, supporting viral replication (Baquero-Perez et al. 2019). Nonetheless, the complexity of the KSHV transcriptome and the likely related epitranscriptome expands beyond coding RNAs and includes many noncoding transcripts, with PAN leading in abundance and significance for KSHV replication (Arias et al. 2014; Chavez-Calvillo et al. 2018; Tagawa et al. 2020). As of today, the m6A landscape and its dynamics for any of these viral noncoding transcripts have not been elucidated.

In this study, we applied an antibody-independent next-generation sequencing-based approach and direct RNA sequencing to tackle the PAN m6A landscape during the latent and lytic stages of KSHV replication. We developed a novel approach, SedTTP-RT and ligation assisted PCR analysis (SLAP), to quantify the modified fraction, which is a critical parameter for investigating m6A biological significance. Using target-specific RNA cross-linking, affinity capture and mass-spectrometry, we delineated the cellular methylome that guides m6A installation on PAN RNA and potentially conveys its phenotypic effect. Furthermore, confocal microscopy analyses identified the temporal and spatial subcellular availability of these components for PAN m6A modification. We also addressed the influence of m6A on the local and global secondary structure of full-length PAN RNA by performing biochemical RNA structure probing (SHAPE-MaP). This study provides the first comprehensive insight into the m6A status of a specific viral lncRNA, revealing the sophisticated interplay that exists between this viral lncRNA and cellular epitranscriptomic machinery, and creates a paradigm for future studies in the field of host-pathogen interactions.

RESULTS

PAN RNA is most extensively modified during the late lytic stages of KSHV infection

Most m6A mapping studies have been focusing on analyzing the global epitranscriptome of a specific organism or a virus with little insight into the position and frequency of modification on an individual transcript. The critical involvement of PAN RNA in the modulation of viral replication and cellular processes prompted us to examine whether this viral noncoding transcript undergoes m6A modifications, and if its epitranscriptomic status is dynamic. We designed a set of 10 antisense oligonucleotide probes to affinity capture PAN RNA during uninduced stage, that corresponds to KSHV latency (0 h post-induction, h pi), immediate early (8 h pi), early (24 h pi), and late lytic stages of replication (48 and 72 h pi) (Supplemental Fig. 1a), which are defined by distinct gene expression programs in the KSHV infectivity cycle (Arias et al. 2014). The region of the KSHV genome from which PAN RNA is transcribed overlaps the K7/survivin locus (Borah et al. 2011). To differentiate between the epitranscriptomic imprints of PAN RNA and K7 mRNA, we designed an antisense oligonucleotide probe that targets the K7 transcript in a region upstream of the overlapping sequence (Supplemental Fig. 1b; Yakushko et al. 2011). Additionally, the analysis included the antisense oligonucleotide probes specific to U1 small nuclear (sn) RNA to control for the affinity capture efficiency and specificity (data not shown) (Engreitz et al. 2014).

Previously it has been reported that the inclusion of 4-selenothymidine-5′-triphosphate during reverse transcription (4SedTTP RT) can identify the position of m6A on target RNA at single-nucleotide resolution (Hong et al. 2018). The presence of this deoxythymidine triphosphate (dTTP) derivative that carries a selenium atom replacing an oxygen at position 4 weakens its ability to base-pair with m6A, creating a stop at the position immediately downstream from the modified site (+1), while maintaining normal A-T base-pairing. We verified the accuracy of this approach by performing reverse transcription (RT) in the presence and absence of 4SedTTP on a 40 nt long transcript that carried a single m6A. The analysis confirmed the presence of a modification-induced truncation (Supplemental Fig. 2a). Next, we sought to determine the compatibility of the 4SedTTP RT method with next-generation sequencing analysis. We used two control transcripts, one carrying m6A marks at positions 13 and 82 and another unmodified transcript. The transcripts were directed to the 4SedTTP RT and control RT reactions, which instead of 4SedTTP, included thymine triphosphate (dTTP). To account for the m6A positional effect, we included random hexamer primers in all RT reactions. Out of 8394 total reads obtained for the modified transcript in 4SedTTP RT reaction, 1487 carried RT stops at the expected positions, resulting in the estimated m6A ratio of 5.6% (m6ARATIO, defined as the proportion of reads supporting RT stops at a position out of all the reads overlapping it). No reads indicating truncations were detected for control RT reactions performed on modified (dTTP-m6A), unmodified transcript (dTTP-A), and the 4SedTTP RT reaction performed on unmodified transcript (4SedTTP-A). The m6A fold change (m6AFC, defined as the log2 transformed m6A ratio in the treated samples [+4SedTTP] divided by the m6A ratio in the nontreated samples [−4SedTTP]) for two sites on modified transcript (4SedTTP-m6A), was estimated at 8.5 and 8.1, respectively (Supplemental Fig. 2b). The m6AFC calculated for 4SedTTP-m6A/dTTP-m6A libraries showed no significant correlation between the results, while that same ratio calculated for 4SedTTP-A/dTTP-A libraries indicated significant correlation (ρ = 0.99%, pairwise P-value <10−13). This indicated that RT stops, other than the ones induced by modification, were random and nonspecific.

After verifying the accuracy of 4SedTTP RT, we used the method to map m6A on PAN RNA captured at uninduced (0 h pi) and lytic stages of KSHV replication (8–72 h pi). We identified a total of five peaks that mapped to m6A at nt positions 18 and 203, 672, 1041, and 1048 (Fig. 1). These results were reproducible across two distinct biological replicates. The peak representing the modification at nt 18 showed up first on PAN RNA expressed during uninduced stage and became distinct at late lytic stage of KSHV infection (72 h pi). The peak corresponding to m6A at nt 203 was weak at immediate early lytic stage (8 h pi) and became pronounced at both late lytic stages of KSHV infection (48–72 h pi). The peak corresponding to m6A at nt 672 was detectable during uninduced stage and became well-defined at immediate early (8 h pi) and late lytic stages of KSHV infection (48–72 h pi), while at early lytic stage (24 h pi) it was not detectable. The peaks at nt 1041 and 1048 were both the most distinct at late lytic stages (48–72 h pi) (Fig. 1). These data demonstrated that the PAN m6A landscape is dynamic and changes during the KSHV infectivity cycle. They also indicated that PAN RNA expressed at the late lytic stages of infection is the most extensively modified, as it can carry up to five m6A sites.

FIGURE 1.

The m6A landscape of PAN RNA expressed at latent and lytic stages of KSHV replication. (A) The graphs represent m6ARATIO (y-axis) against PAN RNA sequence (x-axis). The peaks correspond to the RT stop 1 nt downstream from m6A. The threshold line (red) for calling m6A peaks is set at 0.2 and it was calculated by dividing the average m6ARATIO values calculated for all PAN nucleotides divided by the total number of nucleotides. (B) The charts represent the signal-to-noise ratio for sequences overlapping m6A peaks within a specified nucleotide window (x-axis). (C) The average m6ARATIO values for each modification at specified time points (n = 2). Time points of KSHV infection are color-coded as follows: uninduced stage/latency—0 h pi in blue; lytic stages: immediate early—8 h pi in orange; early—24 h pi in gray; late—48 h pi in yellow; late—72 h pi in green.

Direct RNA sequencing (DRS) confirms m6A modifications on PAN RNA

The DRS has been adopted to examine the m6A status of full-length RNAs (McIntyre et al. 2019; Lorenz et al. 2020; Parker et al. 2020). It has been shown that the base-calling of reads corresponding to m6A results in a higher mismatch frequency (MF) as compared to the unmodified reads (Liu et al. 2019). Also, base quality (BQ), insertion and deletion (INDEL) frequency, and current intensity can be affected by the presence of m6A (Liu et al. 2019). We sought to determine these values for PAN RNA expressed during uninduced stage (0 h pi) and at the late lytic stage of replication (48 h pi), during which the transcript carried the largest number of modified residues (Fig. 1) and compared them to two control transcripts. We found that the mismatch frequency in the 100% modified transcript was significantly higher than that of the unmodified transcript (P-value for MF = 1.868 × 10−135, t-test, P-value < 0.05) (Fig. 2A). A similar correlation was noted for PAN expressed during KSHV late lytic replication versus uninduced stage (P-value for MF = 1.586 × 10−47, t-test, P-value <0.05). The analysis of mismatch directionality for m6A in 0% and 100% modified control transcripts, and PAN expressed during uninduced (0 h pi) and late lytic stages (48 h pi) of KSHV replication indicated that the modification did not invoke biases towards any specific type of mismatch (Supplemental Fig. 3a). Also, we did not observe that the INDEL frequency, base quality scores, or current intensity were significantly different among analyzed samples (Supplemental Fig. 3b–d).

FIGURE 2.

The results of DRS data analysis for PAN RNA. (A) Mismatch frequency (MF). (B) The presence of m6A on PAN RNA expressed during uninduced stage of KSHV replication (0 h pi, blue), and late lytic stage (48 h pi, orange) induces the distortion in electric current. Similar current distortions are noted for all m6A on 100% modified PAN transcript (yellow) versus unmodified transcript (gray).

Previous studies have shown that DRS raw current intensity signals can be used to identify the current changes invoked by the presence of m6A with single-nucleotide resolution (Parker et al. 2019; Lorenz et al. 2020). In particular, the EpiNano, nanopore base-calling software, trained with m6A modified and unmodified synthetic sequences, has been shown to recognize m6A with ∼90% accuracy in vitro, while in vivo modifications have been identified with ∼87% accuracy and an m6A recovery percentage of around 32% (defined as the number of positive methylation signals divided by the total number of adenosines) (Liu et al. 2019). We verified these findings by analyzing the current intensity signals obtained for 100% modified and unmodified PAN transcripts. Out of 259 m6A sites present in 100% modified control transcripts, we detected 183 current distortions, estimating the recovery percentage at 64% (Supplemental Fig. 3e). No current alterations were noted for the unmodified transcript. We applied this approach to validate whether the presence of m6A on PAN RNA expressed during uninduced (0 h pi) and late lytic stage (48 h pi) of KSHV replication triggers current distortions. We mapped a total of ∼8000 reads corresponding to PAN RNA at every stage of viral replication and assessed the recovery percentage for each m6A site. For PAN RNA expressed during uninduced stage, we observed the highest recovery percentages at sites 18 and 672 (18.4% and 11.2%, respectively), while the other sites had lower recovery rates (203–1.0%, 1041–3.3%, 1048–4.2%, respectively). For PAN RNA expressed during the late lytic KSHV replication, all PAN m6A sites had high recovery percentages (sites 18–35.3%, 203–29.5%, 672–15.3%, 1041–40.9%, and 1048–45.3%).

To validate the precision of the m6A-calling by EpiNano and visualize the results, we used the MotifSeq (SquiggleKit pipeline). The MotifSeq takes a primary query sequence as input (.fasta file), converts it to a normalized signal trace, and performs a signal level local alignment using a dynamic programming algorithm (Ferguson and Smith 2019). The normalization of signal traces occurs by subtracting the expected current values from the observed currents gathered by MotifSeq. This approach is particularly useful for detecting m6A, as signal traces for modified residues have been shown to hold a greater disparity between expected and observed currents (McIntyre et al. 2019). We compared the current values obtained for the unmodified PAN in vitro transcript with current values obtained for 100% modified PAN in vitro transcript, PAN RNA expressed during uninduced (0 h pi), and late lytic (48 h pi) stages of KSHV replication. Two m6A at nt positions 18 and 672 on PAN RNA expressed during uninduced stage triggered current alternations (Fig. 2B). For PAN RNA expressed during late lytic replication (48 h pi), all m6A sites were manifested with current changes. A similar pattern of distortion was differentiating the modified and unmodified transcripts.

The modification frequency at identified m6A on PAN RNA reaches the highest level during the late lytic stages of KSHV infection

The modification abundance at a specific site is a variable function that can be instrumental for its functional importance in the context of specific transcripts. Unfortunately, high throughput epitranscriptomic mapping approaches, which involve multiple chemical handling steps, often result in a higher than desired background that obscures that relevant measurement. Also, the quantification of m6A stoichiometry based on DRS readouts has limitations as detection of modifications is partly dependent on sequence context and stoichiometry prediction is heavily affected by choice of re-squiggling algorithms (Begik et al. 2020). To gain insight into KSHV replication state-dependent stoichiometry of m6A in PAN RNA, we developed a new method, termed selenium-modified deoxythymidine triphosphate reverse transcription and ligation assisted PCR analysis (SLAP), which allows quantitative RT-PCR-based determination of the modification fraction (Supplemental Fig. 4). This method relies on the use of selenium-modified deoxythymidine triphosphate (4SedTTP) that invokes reverse transcriptase (RT) stops 1 nt downstream from m6A, and at the same time allows the formation of canonical pairs with unmodified adenine. The m6A-specific truncation effect is achieved due to the weaker electronegativity of Se, which reduces the probability of hydrogen bond formation (Vercoutere et al. 2003; Kumar et al. 2007). Also, the base-stacking interactions are disturbed by the methyl group at the N6-position of m6A because of the steric hindrance. The resulting truncated m6A-specific products and full-length products corresponding to unmodified sites are site-selectively ligated to adapter oligonucleotides (Supplemental Table 1) to introduce a common forward primer binding site. This forward primer, together with the reverse primer that binds to the upstream sequence shared by both modified and unmodified products, is used to generate two distinct products in the same sample during single PCR reaction. This strategy yields the quantitative information for individual m6A residues.

We first evaluated the accuracy and sensitivity of the SLAP method by setting up a calibration curve using 30 fmole of in vitro synthesized 100% m6A modified and unmodified PAN transcripts (standards), which were pooled at specific percentages (from 2.5% to 100% m6A modified) and spiked with 1 µg of total RNA extracts from BCBL-1 cells (Fig. 3A). For these combined standards we achieved a linear regression of R2 = 0.988 (Fig. 3B). We estimated from the calibration curve that SLAP allows for quantification of the modification fraction at as low as a 2.5% level. We noticed the presence of a very weak m6A-specific product in the positive reaction with only the unmodified transcript (+0%, Fig. 3A), which was regarded as the background, and subsequently a twofold threshold above that background was used to estimate the modification frequency. Minor background products were also noted for negative (−) control reactions but only for some m6A sites, that is, at nt 18, suggesting the occurrence of nonspecific-to-m6A truncations in some RT reactions. Interestingly, the sample including a 100% modified transcript resulted in an estimation of m6A frequency at a 92% level, which suggests a slight underestimation of the actual modification fraction (Fig. 3B).

FIGURE 3.

Selenium-modified deoxythymidine triphosphates reverse transcription and ligation assisted PCR analysis (SLAP) quantifies the modification frequency on PAN RNA. (A) Native PAGE shows PCR products derived from the modified (100% m6A) and unmodified (0% m6A) RNA standards, combined at the indicated percentages and subjected to SLAP analysis. The negative (−) and positive (+) reactions are indicated on top, together with the percentage of modified RNA. (B) The graph represents the SLAP standard calibration curve created based on quantification of results from panel A. The x-axis represents m6A fractions tested, while the y-axis represents corresponding densitometric measurements (n = 3). (C) Native PAGE for representative replicates (n = 3) of the SLAP analysis performed for nt 18, 203, 672, 1041, and 1048 on PAN RNA. Negative controls for the SLAP analysis include unmodified adenines on PAN RNA at nt 366 and 410 (n = 3). The time points of infection are indicated on top (0–72 h pi). The position of products specific for modified and unmodified sites is labeled. L stands for 50 bp DNA ladder. (D) The column graphs represent the average modification frequency on PAN at each site and the time point of infection (n = 3). Standard deviations for frequency measurements are indicated.

Next, we applied the SLAP analysis to quantify the frequency of m6A at each identified site on PAN RNA and two nonmodified adenine residues on PAN RNA at sites 366 and 410 (negative controls) (Fig. 3C). Due to the proximity of m6A at nt 18 to the 5′ terminus, analysis of that site required the use of an extended adapter containing the forward primer binding site, which would allow size-based differentiation of modified and unmodified products. As a result, the m6A-specific product corresponding to that site is longer (127 nt) than the product corresponding to unmodified residue (100 nt). Overall frequency of modification at all sites was relatively low for PAN RNA expressed during uninduced and early lytic stages of KSHV replication reaching the highest levels during late lytic stages (48–72 h pi). In particular, the m6A at position 18 was estimated at 31.4 ± 3.5% during uninduced stage, decreasing at immediate early stage to 0 ± 1.7%, followed by an increase through lytic and late lytic phases to 48.0 ± 2.9%, 56.6 ± 2.7%, and 57.9 ± 2.6%, respectively. The modification level at nt 203 was estimated at 3.4 ± 4.8% during uninduced, 59.7 ± 2.5% during immediate early, and 53.5 ± 4.0% at early stages of viral replication, followed by an increase to 52.2 ± 5.2% and 51.4 ± 2.8% during the late lytic stages (48–72 h pi). Position 672 was modified at 7.4 ± 5.0% during uninduced stage, 19.6 ± 3.7% at immediate early, 0 ± 5.1% at early lytic stages, followed by an increase to 16.3 ± 4.2% and 29.7 ± 3.2% during late lytic stages of KSHV replication. Nucleotide 1041 was modified at 10.8 ± 5.2% during uninduced stage, 1.0 ± 1.6% during immediate early, 0 ± 1.3% during early lytic stages, and showed the highest frequency during the late lytic stages, reaching a frequency of 57.4 ± 5.3% and 55.9 ± 4.9%, respectively. The modification frequency at nt 1048 was detectable at 10.4 ± 7.5% level during uninduced stage, reaching a frequency of 0 ± 0.7% and 9.5 ± 1.6% for immediate early and early lytic stages. However, the frequency of m6A at that site again increased during late lytic stages of KSHV replication, 28.2 ± 9.8% and 33.4 ± 3.7%, respectively (48–72 h pi).

The m6A modifications in the human transcriptome have been primarily localized within two canonical motifs, Gm6AC (∼70%) and Am6AC (∼30%) (Wei and Moss 1977; Wei et al. 1976), which have been ascribed to the METTL3/14 complex (Dominissini et al. 2012; Meyer et al. 2012), although another study identifies the GGm6AGG motif as present at the 41% frequency (Wongsurawat et al. 2018). A recent study performed on the KSHV transcriptome has reported that approximately 60% of all m6A is contained within a GGm6AC[G/U] motif (Baquero-Perez et al. 2019). We performed a motif search to investigate the m6A consensus sequence specific to PAN RNA. The overall PAN m6A consensus sequence was designed as [A/U/C][A/G/C]m6AC[A/U/C], encompassing previously identified canonical motifs, which is suggestive of the involvement of the METTL3/14 complex (Supplemental Fig. 5). In particular, the m6A at position 18 was located within a Cm6AC motif, three modifications at nt positions 203, 1041, 1048 were located at an Am6AC motif, while the m6A at nt 672 was associated with a Gm6AC motif. We also analyzed the m6A consensus for cellular lncRNAs that have been previously reported as m6A modified. Here, we found that the majority of m6A modifications were located within an A/Gm6AC consensus sequence (Supplemental Fig. 5).

PAN RNA associates with m6A writers, readers, and erasers

The RNA m6A status is dynamically coordinated by the catalytic activity of writer and eraser proteins, while the phenotypic effect of m6A is facilitated by various groups of reader proteins (Bokar et al. 1997; Liu et al. 2014; Schwartz et al. 2014b; Balacco and Soller 2019). To identify the methylome that is involved in the regulation of PAN m6A status, we performed cross-linking of PAN RNA-protein complexes during uninduced (0 h pi) and late lytic (48 h pi) stages of KSHV replication, followed by affinity capture and mass spectrometry analysis (RAP MS) (McHugh and Guttman 2018). RAP MS relies on 254 nm UV light cross-linking to create covalent bonds between RNA directly interacting with proteins (Brimacombe et al. 1988; Hockensmith et al. 1993). The biggest caveat of UV cross-linking is its low efficiency, which has been estimated at ∼5% for any given ribonucleoprotein complex. Thus, we also performed formaldehyde (FA) cross-linking, which occurs through the formation of methylol derivatives between amino groups of cytosines, guanine, and adenines or the imino groups of thymines with amino acids such as lysine, arginine, tryptophan, and histidine (Solomon and Varshavsky 1985). Both analyses were performed in three independent experimental replicates to account for methodological bias (Fig. 4A). We assessed the relative abundance of PAN RNA by RT qPCR and showed that the target was efficiently captured at both time points (over 100-fold enrichment for both uninduced and late lytic stages, Supplemental Fig. 6). The proteins that were captured in cross-linked samples and not present in the non-cross-linked samples were compared statistically between methods and the replicate samples. This categorical analysis of PAN-associating proteins and their mass spectrometry MASCOT scores showed that the number and type of proteins bore a strong correlation between the two cross-linking methods (FA and UV, R2 = 0.9894) and among individual biological replicates (R2= 0.9827 for replicates I and II [FA], R2= 0.9895 for replicates I and II [UV], R2= 0.9954 for II and III [FA], R2= 0.9803 for II and III [UV], R2= 0.9901 for I and III [FA], R2= 0.9918 for I and III [UV], Fig. 4A). We also performed principal component analysis (PCA) on the data obtained from RAP MS with UV and FA cross-linking and found that the number and type of PAN-associating proteins shared 99% overlap between both cross-linking methods, and among the three distinct biological replicates (Fig. 4B). For PAN RNA expressed during the KSHV late lytic replication (48 h pi), we identified RBM15, a protein crucial for the recruitment of METTL3/14 writer complex (Ma et al. 2007) and demethylase FTO. We also identified three m6A reader proteins, that is, HNRNPC, YTHDF2, and SND1. Importantly, none of these proteins were found to associate with PAN RNA expressed during uninduced stage of KSHV replication. Based on these data, we established the lncRNA-protein interactions network (LPI), including the identified PAN methylome components, and predicted interacting proteins derived from STRING, the functional protein association network database (Supplemental Fig. 6b; Szklarczyk et al. 2017). The network identified five nodes (proteins) and eight edges (interactions) associated with PAN modification. We classified RBM15 as a central hub for the PAN m6A methylome, which was predicted to interact with the METTL3/14 writer complex, accessory proteins, that is, RNA binding protein 15 (RBM15), Wilms’ tumor 1-associating protein (WTAP), and identified reader proteins, that is, SND1, YTHDF2, HNRNPC. The FTO, PAN-associating eraser was predicted to bind METTL3, RBM15B, and two identified reader proteins, that is, SND1 and YTHDF2. Overall, the LPI network was interconnected within each functional subset of PAN m6A methylome and between writers, erasers, and readers.

FIGURE 4.

The RAP MS analysis for PAN RNA m6A methylome and effect on PAN m6A methylation frequencies. (A) Venn diagrams representing the overlap between results generated from formaldehyde (FA) and UV cross-linked PAN-protein complexes. The three distinct biological replicates are indicated with colors (I-III). (B) Principal component analysis (PCA) of proteins identified in RAP MS (UV, x-axis) and FA experiments (y-axis). Proteins involved in regulation of the PAN m6A status are marked with purple dots, and they include the writer complex recruiter RBM15, readers HNRNPC, YTHDF2, SND1, and the FTO eraser. (C) Normalized expression of PAN RNA assessed by RT qPCR in non-cross-linked total RNA samples extracted during KSHV latent (0 h pi) and lytic (48 h pi) stages of replication, and after the PAN RNA–protein cross-linking and affinity capture. (D) Results of western blot analyses performed on whole cell BCBL-1 lysates (left) and captured PAN RNA-protein complexes (right) during uninduced (0 h pi) and lytic (8–72 h pi) stages of KSHV replication. β-actin was used as a control.

We further verified the associations of PAN RNA with the methylome components by performing a western blot analysis. We first assessed the expression of individual methylome proteins in BCBL-1 whole-cell lysates and showed that their expression levels were not significantly affected by the progression of KSHV replication. After performing FA cross-linking and affinity capture of PAN RNA-protein complexes, we assessed the relative abundance of target RNA and showed that it was efficiently captured at all time points (Fig. 4C). We confirmed that all proteins identified by RAP MS analysis interact with PAN RNA (Fig. 4D). RBM15 and SND1 were shown to interact with PAN at all lytic stages of KSHV infection (8–72 h pi), but not during uninduced stage (0 h pi). FTO was captured as interacting with PAN during all time points (0–72 h pi), while two reader proteins, HNRNPC and YTHDF2, were detected at the late lytic stage of KSHV infection (48–72 h pi and 72 h pi, respectively). METTL3 was not captured as directly interacting with PAN RNA, perhaps because the METTL3/14 writer complex is known to interact with the target RNA via RBM15 association (Patil et al. 2016). We also assessed the normalized expression levels of mRNAs encoding the cross-linked proteins and showed that their expression levels varied throughout the progression of the KSHV infectivity cycle. Only YTHDF2 mRNA showed significantly increased expression during the lytic stages (Supplemental Fig. 7a); however, that increase was not reflected in the expression of YTHDF2 protein (Fig. 4D). RBM15 mRNA expression level decreased during immediate early (8 h pi) and late lytic (48 h pi) stages of infection, while displaying a slight increase during early (24 h pi) and last late lytic stages (72 h pi). METTL3, FTO, and SND1 mRNA expression levels showed a slight increase during most lytic stages of KSHV replication (Ye et al. 2017; Baquero-Perez et al. 2019), while HNRNPC mRNA stayed at a relatively unchanged level, with some decrease observed during immediate early stage (8 h pi) of KSHV replication.

The knockdowns of methylome components affect PAN RNA expression levels and its m6A landscape

To assess the involvement of METTL3, RBM15, and FTO in the regulation of PAN m6A status, we performed the siRNA-directed knockdowns (KD) of individual proteins in BCBL-1 cells, followed by induction of KSHV lytic replication. The RT qPCR analysis indicated nearly complete expressional ablations of RBM15 (∼83%), METTL3 (∼89%), and FTO (∼92%,) at mRNA levels for each time point of infection (Fig. 5A). Western blots further verified the efficiency of knockdowns for all three proteins (Fig. 5B).

FIGURE 5.

The knockdowns of m6A methylome components affect PAN RNA expression levels. The normalized mRNA expression (y-axis) of (A) RBM15, METTL3, and FTO in BCBL-1 knockdowns (KD) cell lines. The x-axis corresponds to uninduced (0 h pi) and lytic (8–72 h pi) stages of KSHV infection. Control samples represent latent BCBL-1 cells treated with scrambled siRNAs. (B) Western blots representing the efficiency of knockdowns on the level of protein expression. Tubulin was used as a control. (C) PAN RNA expression in knockdown cell lines versus control at specified time points of KSHV replication. (D) ORF50 (left) and ORF57 (right) mRNA expression levels in knockdown versus wild-type BCBL-1 cells (control) at indicated time points post-induction. MALAT1 and β-actin were used for data normalization as their expression levels were not significantly affected by the knockdowns and progression of KSHV infectivity cycle (n = 4). (E) The representative native PAGE indicating frequency of modification at m6A sites in BCBL-1 RMB15, METTL3, and FTO knockdown (KD) cells. L represents 50 bp DNA ladder. The analyzed time points of infection (0 h pi, 48 h pi), negative (−4Sed) and positive (+4Sed) reactions are indicated on top. The position of products specific for modified and unmodified sites is indicated to the left. The column graphs represent the average modification frequency for each site at 0 h pi and 48 h pi time point of infection in specified knockdown cell lines. Standard deviations for frequency measurements are indicated.

A previous study has shown that the functional inhibition of FTO by BCBL-1 treatment with meclofenamic acid, an FTO inhibitor, and FTO ablation with shRNAs, increases the expression of KSHV lytic genes on both mRNA and protein level at early lytic stage of infection. In contrast, the impediment of m6A by treatment with DAA, a METTL3 inhibitor, and METTL3 ablation with shRNAs, abolished the lytic genes’ expression (Ye et al. 2017). We wanted to assess how the METTL3, RBM15, and FTO ablations affect PAN RNA expression levels, as this has not been previously addressed, and can provide insight into the functional significance of m6A for this lncRNA. For the FTO knockdown, the PAN RNA expression level initially increased at uninduced (0 h pi) and immediate early stages of KSHV infection (8 h pi), followed by a sharp decline at early and both late lytic stages (24–72 h pi) (Fig. 5C). For METTL3 knockdown, there was a strong decrease in PAN expression at most time points of infection (0–48 h pi) except for the last late lytic stage (72 h pi). In the case of the RBM15 knockdown, PAN RNA expression was decreased at all time points of viral replication. We also assessed the expression levels of replication and transcription activator (RTA/ORF50) mRNA and mRNA transcript accumulation (Mta/ORF57) in knockdown cell lines, as these key transcripts code for early viral products that influence PAN RNA expression and stability (Fig. 5D; Song et al. 2001; Majerciak and Zheng 2015). Also, the metabolism of both transcripts has been shown to be affected by the knockdown of methylome complexes (Ye et al. 2017; Hesser et al. 2018; Baquero-Perez et al. 2019). We anticipated that the expressional pattern of ORF50 and ORF57 mRNAs would follow PAN RNA expression levels. Interestingly, FTO knockdown slightly increased the expression of ORF50 and ORF57 mRNAs compared to wild-type BCBL-1 cells. Since FTO knockdown caused the decline in PAN RNA expression at the early (8 h pi) and both late lytic stages of KSHV replication (48–72 h pi) (Fig. 5C), these effects cannot be explained by merely the influence of ORF50 and ORF57 expression levels. Also, both RBM15 and METTL3 knockdowns reduced ORF50 and ORF57 mRNA expression from two- to fivefold, depending on the analyzed lytic stage of infection. However, METTL3 knockdown caused an increase of PAN RNA expression at the last late lytic stage (72 h pi).

Next, we applied the SedTTP RT with next-generation sequencing and SLAP methods to examine PAN m6A landscape in knockdown cell lines and to validate the involvement of the methylome components in PAN modification. As anticipated, the knockdown of METTL3 and RBM15 completely abolished the installation of m6A on PAN RNA at all identified sites, and at both uninduced (0 h pi) and late lytic (48 h pi) stages of infection (Supplemental Fig. 8). Meanwhile, the FTO knockdown resulted in an m6A pattern like the one observed in wild-type BCBL-1 cells (Fig. 1). Peaks of slightly lower m6ARATIO were observed at nt 18 and 672 during latency (0 h pi), and at all lytic stages of KSHV replication (8–72 h pi). We verified these effects by performing the SLAP analysis on PAN RNA expressed in each of the knockdown cell lines at uninduced (0 h pi) and late lytic (48 h pi) stages (Fig. 5E). Again, the m6A-specific products were not detected for METTL3 and RBM15 knockdowns at both analyzed stages of KSHV infection, and for FTO knockdowns, each site previously identified as modified resulted in generation of m6A-specific products. In general, for FTO knockdowns, we observed a slight increase in m6A frequency during uninduced stage, and a decrease in modification frequency at the late lytic stage (Supplemental Fig. 8). These results correlate with the observed increased levels of PAN RNA expression during uninduced and immediate early stages in FTO knockdown cell line, followed by a sharp decline observed at late lytic stages, and likely can account for the higher levels of m6A peaks. Most importantly, however, they emphasize the importance of FTO in the dynamic regulation of PAN m6A landscape. For site 18, the modification frequency reached 35.3 ± 1.8% during uninduced stage (WT—31.4 ± 3.5%), and 15.2 ± 1.3% during late lytic stage (WT—56.6 ± 2.7%). For site 203, the modification levels were 13.0 ± 6.3% (WT—3.4 ± 4.8%) and 10.5 ± 1.2% (WT—52.2 ± 5.2%) for latent and late lytic stages, respectively. For the residue at position 672, the modification levels were at 17.5 ± 12.7% (WT—7.4 ± 5.0%) and 13.8 ± 3.1% (WT—16.3 ± 4.2%), respectively. However, the modification abundance at sites 1041 and 1048 increased for both time points of KSHV infection and reached 24.2 ± 2.7% (WT—10.8 ± 5.2%) and 33.2 ± 1.3% (WT—57.4 ± 5.3%) for site 1041, and 35.4 ± 4.2% (WT—10.4 ± 7.5%) and 35.5 ± 1.0% (WT—28.2 ± 9.8%) for site 1048 (0 h pi and 48 h pi, respectively).

The m6A modification affects PAN RNA secondary structure at local and global levels

We performed selective 2′-hydroxyl acylation analyzed by primer extension and mutational profiling analysis (SHAPE-MaP) of PAN RNA expressed in wild-type, METTL3, and FTO knockdown BCBL-1 cell lines to gain insight into the structural impact of m6A. We chose to probe the structure of PAN RNA expressed during the late lytic stage of infection (48 h pi), as at that time point the transcript carried the largest number of modified residues (Fig. 1) at the highest modification frequency (Fig. 3). It has been shown that m6A methylation occurs at sites that are mostly unpaired (Kierzek and Kierzek 2003; Schwartz et al. 2013). Hence, we anticipated that the expressional ablation of METTL3 would trigger more double-stranded and constrained structures, while the knockdown of FTO would invoke more structurally flexible, that is, unpaired regions in PAN RNA. Since probing of PAN RNA secondary structure in wild-type BCBL-1 cells (Supplemental Fig. 9) would be expected to yield the ensemble of reactivity values at m6A reflecting varying modification frequency and distribution, we compared the reactivity profiles (ΔSHAPE analysis) for PAN RNA expressed in METTL3 and FTO knockdown cells (Fig. 6A; Smola et al. 2015; Sztuba-Solinska et al. 2017). As expected, we observed that the ablation of writer invoked lower reactivity values at m6A, while the knockdown of eraser induced strong structural signatures marked by increased exposure of m6A residues to the modifying reagent. The noted structural changes were of two type—local, defined as within 20 nt to the identified m6A (nt 25–32, 195–210, 703–709, and 1039–1050), and global, defined as outside the 20 nt range (nt 120–130, 151–158, 820–827, 900–950) (Fig. 6A). The most pronounced global changes invoked by the ablation of METTL3 were noted within the expression and nuclear retention element (ENE), which is adjacent to two m6A at nt 1041, and 1048 (Conrad et al. 2007). The PAN ENE forms a hairpin containing a U-rich internal loop flanked by short helices, which is involved in the sequestration of poly(A) tail, leading to the formation of a stabilizing triple helix (Conrad et al. 2007; Mitton-Fry et al. 2010). Interestingly, in METTL3 knockdown cells, the ENE U-rich residues were split between two hairpins (Fig. 6B). Such extensive secondary structure rearrangement was not noted for PAN RNA expressed in the FTO knockdown BCBL-1 cells, which was predicted to invoke mostly local structural rearagnements.

FIGURE 6.

The structural analysis of PAN RNA in METTL3 and FTO knockdown BCBL-1. (A) PAN RNA SHAPE-MaP reactivity profiles with blue trace corresponding to reactivity values obtained in METTL3 KD cells, red trace representing reactivity values from FTO KD cells. X-axis indicates nucleotide position on PAN RNA, y-axis represents SHAPE-MaP reactivity values. Nucleotides 1–14 and 1056–1077 were not covered in the analysis due to primer annealing. The green boxes indicate the local and global structural changes invoked by the methylome ablations. The position of the expression and nuclear retention (ENE) motif is indicated. (B) The secondary structure of PAN ENE motif (green) probed in wild-type, FTO and METTL3 knockdown (KD) cells. Nucleotides are colored by SHAPE reactivity as follows: gray residues display low reactivity (<0.2), pink residues representing medium reactivity (0.2–0.8), and red values representing high reactivity (>0.8). (C) Fraction of structural changes estimated for PAN RNA SHAPE profile obtained in wild-type BCBL-1 cells vs wild-type (background), FTO, and METTL3 knockdown cells. No change is marked in blue, local change is indicated in orange, while global change is colored in gray. (D) ROC curves generated for PAN RNA SHAPE profiles representing the m6A influence over PAN secondary structure. The x-axis corresponds to false positive rate (FPR); the y-axis represents true positive rate (TPR). The AUC column chart (right) indicates the level of agreement between nucleotide reactivity and pairing status. A value of 0.5 indicates that nucleotide reactivity reports no base-pairing information.

To further classify the PAN RNA secondary structure changes triggered by the METTL3 and FTO knockdowns, we performed classSNitch analysis (Fig. 6C,D; Sztuba-Solinska et al. 2017; Woods and Laederach 2017). The classSNitch categorizes RNA structural changes for each nucleotide based on the comparison of wild-type and experimental conditions by ranking eight features, for example, Pearson correlation coefficient, contiguousness, change variance, etc., to determine local, global, or no changes to RNA secondary structure (Woods and Laederach 2017). We compared PAN SHAPE reactivity profiles obtained in METTL3 and FTO knockdown conditions to the wild-type reactivity profile. The wild-type reactivity profile, which served as a base input, was initially assessed for the background levels of global and local changes by comparing two independently obtained wild-type profiles. It showed that 99% of residues did not exhibit structural changes, while 6.8% and 3.2% of residues displayed local and global changes, respectively. The PAN RNA reactivity profile obtained in FTO knockdown cells included 67% residues that did not exhibit structural changes, while 32.5% and 0.5% residues displayed local and global changes, respectively. The PAN RNA reactivity profile in METTL3 knockdown cells was characterized by 39.2% residues that displayed no structural changes, 57% residues with local changes, and 3.8% residues that showed global changes (Fig. 6C). Overall, PAN RNA profile obtained in METTL3 knockdown cells displayed over sevenfold greater level of global changes than FTO KD. These results indicate that m6A can impact not only the local PAN RNA secondary structure but also its global folding, and as such, can hold substantial control over the transcript's biology. Additionally, we generated the receiver operator curve (ROC, comparison between RandomForest RNA classifier, wild-type, and knockdowns) for each PAN RNA SHAPE profile to estimate the accuracy of secondary structure prediction. The PAN RNA secondary structure profile obtained in FTO knockdown cells had ROC curve with an estimated area under the curve value of 0.86 (AUC, measurement of accuracy) (Fig. 6D). In contrast, the PAN RNA profile obtained in METTL3 knockdown had ROC curve with an estimated AUC value of 0.77 (Fig. 6D). AUC values ranging from 0.6–0.7 are regarded as having low accuracy.

PAN colocalizes with the m6A methylome components to varying degrees during the KSHV lytic replication

PAN RNA has been considered as a mainly nuclear lncRNA (Borah et al. 2012), with some reports suggesting that a fraction of the PAN transcript can leak to the cytoplasm to support specific intermolecular interactions (Arias et al. 2014; Sztuba-Solinska et al. 2017). Similar observations have been made for other predominantly nuclear lncRNAs, for example, metastasis-associated lung adenocarcinoma transcript 1 (MALAT1) (Yang et al. 2013). We performed immunofluorescence analysis to assess the temporal and spatial subcellular distribution of our previously identified methylome enzymes, their potential for colocalization with PAN, and dependance on KSHV replication stage. We manually scored the immunofluorescence pattern of on average 10 cells from a number of microscopic fields and at each time point of infection, and found that PAN RNA was detectible as early as 8 h post lytic induction at the immediate early stage of KSHV infection (8 h pi). PAN RNA accumulated to the higest level at the early lytic stage (24 h pi), and was located mainly within the nucleus (Fig. 7; Supplemental Fig. 10).

FIGURE 7.

Confocal microscopy analysis of PAN and the indicated methylome components performed in BCBL-1 cells at uninduced (0 h pi) and lytic stages (8–72 h pi) of KSHV replication. The confocal microscopy micrographs are accompanied by thresholded Pearson Correlation Coefficients (PCC) with colocalization that were assessed for each bioreplicate, and fluorescence intensity profiles for indicated magnified cells and identified interactions: (A) PAN-RBM15, (B) PAN-METTL3, (C) PAN-FTO, (D) PAN-HNRNPC, (E) PAN-SND1, and (F) PAN-YTHDF2. DAPI staining of the nuclei is in blue, the methylome components are in pink and PAN RNA in green (n = 3). The fluorescence intensity profiles include lines color-coded according to staining; the red dashed line outlines the nucleus.

Recently, m6A-dependent colocalization of MALAT1 with RBM15 has been described, where RBM15 appeared to form small foci scattered over the nucleus and nuclear speckles (Wang et al. 2021). Our results are consistent with this observation, as we noted that RBM15 formed minute puncta that were contained within the nucleus. We calculated thresholded PCC using the JACoP plugin after selecting the “Costes automatic thresholding” option. The plugin sets thresholds for each image using the algorithm of Costes et al. (2004) and calculates PCC for pixels brighter than the threshold, and PCC for pixels fainter than one of both of the thresholds (Costes et al. 2004; Barlow et al. 2010). The PAN and RBM15 colocalization was most noticeable at early (24 h pi, thresholded PCC = 0.83, 0.89) and late lytic (72 h pi, thresholded PCC = 0.91, 0.89) stages of viral lytic replication (Fig. 7A). Similarily, METTL3 formed a small foci within the subnuclear space (Liu et al. 2014; Ping et al. 2014), and colocalization with PAN RNA was most evident during late lytic stage of KSHV infection (72 h pi, thresholded PCC = 0.93, 0.91) (Fig. 7B). We noted that the FTO localized abundantly to the nucleus (Bartosovic et al. 2017; Yang et al. 2018), and colocalized with PAN RNA to the highest degree during early (24 h pi, thresholded PCC = 0.94, 0.90) and late lytic (48 h pi, thresholded PCC = 0.89, 0.79) stages of viral replication (Fig. 7C). In the case of the m6A readers, HNRNPC resided abundantly within the nucleus, as previously reported (Anantha et al. 2013), and colocalized with PAN to the highest degree at early (24 h pi, thresholded PCC = 0.81, 0.83) and late lytic (48 h pi, thresholded PCC = 0.89) stages of infection (Fig. 7D). Both SND1 and YTHDF2 were found to have profuse cytoplasmic localization, with small fractions of both enzymes being subnuclear. PAN and SND1 colocalization was detected during early lytic stage (24 h pi, thresholded PCC = 0.9, 0.83) in the nucleus (Fig. 7E), while PAN and YTHDF2 colocalization was noted at late lytic stage of KSHV infection (72 h pi, thresholded PCC = 0.93, 0.89), and occurred within the nucleus (Fig. 7F). These data highlight that the progression of KSHV replication does not influence the localization of the m6A methylome components identified to influence PAN RNA modification, and that all of them are temporarily and spatially available to regulate or mediate the phenotypic effect of PAN m6A.

DISCUSSION

The intricate nature of the KSHV transcriptome facilitates the expression of numerous long noncoding RNAs, which regulate viral and cellular gene expression by means that are not fully understood (Campbell et al. 2014b; Schifano et al. 2017; Chavez-Calvillo et al. 2018; Tagawa et al. 2020). These noncoding transcripts themselves endure regulation by various mechanisms, including epitranscriptomic modifications. Recent studies have revealed the critical impact of m6A modification on the progression of KSHV latent and lytic stages of replication (Yang et al. 2018; Zaccara et al. 2019). No studies, however, addressed the dynamics of the m6A landscape and its abundance in the context of specific noncoding viral transcript, and during the progression of viral replication.

In our study, we provide the first comprehensive insight into the m6A status of KSHV-encoded PAN RNA, the key regulator of viral lytic reactivation (Rossetto et al. 2013; Withers et al. 2018; Hiura et al. 2020). Using two independent mapping approaches, that is, the 4SedTTP RT with next-generation sequencing and direct RNA sequencing, we revealed that PAN RNA can carry up to 5 m6A sites. These modifications localize within PAN RNA secondary structure motifs that we have previously shown to display low SHAPE and low Shannon entropy (Sztuba-Solinska et al. 2017), which are characteristic of biologically significant RNA motifs (Siegfried et al. 2014; Smola et al. 2016; Dethoff et al. 2018; Huber et al. 2019). We noted that PAN m6A status is dynamic and changes during the progression of KSHV replication, achieving the highest level during the late lytic stages. This process might operate as a fine-tuning mechanism, coordinating the function and biology of PAN RNA, as it has been demonstrated for other transcripts (Liu et al. 2015; Park et al. 2019; Zhou et al. 2019; Lee et al. 2020). PAN RNA is a highly multifunctional transcript involved in the interactions with viral (Rossetto et al. 2013; Campbell et al. 2014a) and cellular proteins the (Rossetto and Pari 2012), modulating the chromatin and episomal DNA state during viral lytic reactivation (Strahan et al. 2017; Hiura et al. 2020), and facilitating the nuclear export of late viral mRNAs (Withers et al. 2018). The dynamic m6A modifications to PAN RNA might provide the means to orchestrate its functions with the immediate need of the pathogen. We established the PAN m6A consensus sequence at the [A/U/C][A/G/C]m6AC[A/U/C], and noted that the motif shows higher sequence variability compared to the one we identified for other lncRNAs, but in general encompasses all canonical motifs that have been previously identified for m6A and METTL3/METTL14 complex (Baquero-Perez et al. 2019).

Using a newly developed method that we termed selenium-modified deoxythymidine triphosphates (SedTTP) RT and ligation assisted PCR analysis (SLAP), we gained insight into the modification frequency at identified m6A sites. We reasoned that the quantification of the extent of modification could advance the understanding of how changes in PAN modification pattern, the site and fraction, are modulated during the KSHV replication. This is also a crucial parameter for investigating the biological significance of m6A. For example, it has been shown that only one out of fourteen m6A detected on HOX antisense intergenic lncRNA (HOTAIR), which displayed the most consistent modification, plays a critical role in HOTAIR-mediated cellular proliferation, colony formation, and HOTAIR nuclear retention (Porman et al. 2020). Also, one out of four m6A detected in MALAT1, which shows over 80% modification frequency in two cell lines (Liu et al. 2013), is attributed to the m6A-induced structural change associated with m6A reader binding (Zhou et al. 2016). In our studies, the measurement of the m6A frequency at each site identified on PAN RNA was the highest during the late lytic stages of KSHV infection, with three m6A at nt positions 18, 203, and 1048 reaching over 50% modification rate.

We also tackled the identification of methylome components that regulate PAN m6A status. Using a combination of proteomic approaches and confocal microscopy analysis, we showed that the writer METTL3, writer recruiter RBM15, and eraser FTO, are not only temporarily and spatially available for PAN modification but also directly interact with the target RNA. Knockdown of any of these enzymes affected PAN RNA expression to a various degree. FTO ablation caused an increase in PAN expression during uninduced and immediate early stages of KSHV replication followed by a sharp decline at early and both late lytic stages. These effects were not reflected by the comparable expressional patterns of ORF50 and ORF57 mRNAs. For METTL3 knockdown, there was a decrease in PAN expression at most time points of infection, except for the last late lytic stage where we observed a surge in PAN RNA levels, which again, was not reflected in the increased expression of ORF50 and ORF57 mRNAs at that time point of KSHV infection. Both ORF50 and ORF57 are known regulators of PAN RNA expression and stability (Song et al. 2001; Majerciak and Zheng 2015). Also, their expression levels and metabolism were shown to be affected by the knockdown of METTL3 and FTO (Ye et al. 2017; Baquero-Perez et al. 2019). This raised the possibility that the observed perturbed expression of PAN RNA could result from the indirect effect of ORF50 and ORF57 mRNA dysregulation in knockdown cell lines. However, since ORF50 and ORF57 mRNA expression patterns do not follow the same kinetics pattern that we established for PAN RNA in tested knockdown conditions, other factors related to m6A-directed regulation of PAN RNA biology should be considered. The significance of PAN RNA m6A modification likely reaches beyond just the potential modulation of its interactions with proteins but may also include the influence over PAN RNA stability, structure, and turnover.

Next, we focused on identifying the m6A readers of PAN RNA, that could facilitate the phenotypic effect of the modification. Using RNA affinity capture with mass spectrometry and western blot analyses, we identified three PAN m6A readers, that is, SND1, HNRNPC, and YTHDF2. SND1 has been reported to bind and stabilize m6A modified viral ORF50 mRNA (Baquero-Perez et al. 2019). That study also identified the SND1-binding motif as a U-tract immediately followed by the Gm6AC motif. We found similar motif near m6A at position 672 on PAN RNA. Notably, U-rich motifs found adjacent to m6A are also bound by RBM15 and HNRNPC (Liu et al. 2015; Patil et al. 2016), both of which we found to associate with PAN RNA. Another PAN m6A reader, the YTHDF2, has been shown to recruit the CCR4-NOT deadenylase complex and promote the degradation of modified transcripts (Du et al. 2016). In KSHV replication, the knockdown of YTHDF2 has been shown to trigger elevated expression of viral transcripts and release of virions, which is suggestive of YTHDF2 functioning as a restriction factor (Tan and Gao 2018). Another study, however, reported a contradictory result, indicating that YTHDF2 mediates viral gene expression and virion production in iSLK.219 cells (Hesser et al. 2018). We showed that YTHDF2 directly associates with PAN RNA at the late lytic stages of KSHV infection (48 and 72 h pi), during which our target transcript is the most extensively modified. YTHDF2 has been shown to promote mRNA degradation by binding to m6A-modified RNAs (Hou et al. 2021) and recruitment of the CCR4-NOT deadenylase complex (Du et al. 2016; Park et al. 2019). Thus, it is feasible that the YTHDF2 might be involved in the regulation of PAN RNA turnover (Wang et al. 2014; Park et al. 2019; Huang et al. 2020).

We also delineated the structural effect of m6A in the context of full-length PAN RNA. The m6A modification is known to diminish the ability of RNA to form duplexes, favoring the linear, unpaired form of RNA. This reduction in structure provides greater access to RNA by RNA-binding proteins that associate with single-stranded motifs (Liu et al. 2015, 2017; Wu et al. 2018). Another study has shown that the m6A-directed structural effect strongly depends on the modified RNA's structural context, and when neighboring by a 5′ bulge, m6A stabilizes rather than disrupts base-pairing (Liu et al. 2018). In our study we showed that the knockdown of FTO and METTL3 results in the alternation of the local and global secondary structure of PAN RNA. The most striking structural change was invoked by METTL3 ablation, that disrupted the formation of a U-rich loop, which is a part of the ENE, motif involved in the sequestration of PAN poly(A) tail and regulation of transcripts stability (Conrad et al. 2007; Mitton-Fry et al. 2010). In our studies, the region downstream from the ENE was found to carry two m6A (nts 1041, 1048), which displayed moderate and high modification frequency, respectively (up to 33 and 57%) indicative of their functional importance. It is likely that the disrupted modification of these sites can not only affect the secondary structure of the ENE, but also influence the trajectory and interaction of poly(A) tail with the U-rich motif involved in the triple helix formation. These structural effects can perhaps account for the generally observed decreased level of PAN expression in METTL3 and RBM15 knockdown cells at most stages of KSHV replication. Also, these observations stay in agreement with previous studies showing that m6A can influence RNA stability (Wang et al. 2014), structure, and protein binding (Liu et al. 2015). The m6A holds a powerful grip over RNA biology, that can either result in the disruption of functionally critical RNA motifs or unmasking new ones, which in turn, can influence the cellular fate of the target transcript, including its subcellular localization, stability, and degradation. Protein binding can also induce changes to RNA structure, and since m6A can regulate RNA structural access to their respective protein partners, the observed structural changes of PAN RNA can be the result of either an altered m6A modification pattern or RNA–protein interactions influenced by m6A modifications. Our future studies will aim at addressing the influence of PAN RNA m6A modification over its stability and intermolecular interaction and KSHV replication to identify the functionality of this dynamic chemical signature.

Our work represents the most comprehensive overview of the dynamic interplay that takes place among the cellular epitranscriptomic machinery and a specific viral long noncoding RNA. Our findings provide thought-provoking clues regarding the influence of m6A over PAN RNA structure and stability. PAN is an unusual lncRNA, which has evolved multiple posttranscriptional strategies that lead to its increased stability and abundant accumulation (Conrad et al. 2007; Mitton-Fry et al. 2010); however, as all RNAs, it must undergo decay (Conrad 2016). Since we do not know the protein binders that would associate to and resolve the PAN RNA triple helix leading to the release of its poly(A) tail and subsequent deadenylation, perhaps m6A can afford a conceivable factor regulating that process. Only recently has it been shown that the PAN RNA expression is upregulated in cells with the knockdown of nonsense-mediated decay (NMD) factors (Zhao et al. 2020). Considering that our data show PAN association with YTHDF2, which is known to interact with the NMD complex and trigger m6A-mediated RNA decay (Park et al. 2019), and that METTL3 and FTO knockdowns affect PAN expression levels and the secondary structure of the ENE motif, the m6A modification can indeed hold clues to the conundrum of PAN RNA turnover.

MATERIAL AND METHODS

Cell lines and culture conditions

The KSHV positive body cavity-based lymphoma cell line (BCBL-1, a generous gift from Dr. Denise Whitby, NCI Frederick) was seeded at 2 × 105 cells/mL and grown in RPMI-1640 medium (Thermo Fisher 11875085) supplemented with 10% fetal bovine serum (FBS), 1% Penicillin/Streptomycin and 1% L-Glutamine at 37°C in 5% CO2. The induction of KSHV lytic infection was performed by treating 2 × 107 BCBL-1 cells with sodium butyrate (NaB, EMD Millipore 654833) to a final concentration of 0.3 mM, and collecting cells at 8-, 24-, 48-, and 72-h post-induction (h pi).

In vitro RNA transcription

The first 100 nt of the PAN sequence with T7 promoter was synthesized using oligonucleotides listed in Supplemental Table 1 and served as a template for in vitro transcription of the unmodified transcript (negative control). The RNA with the equivalent sequence but carrying two m6A sites was purchased from Integrated DNA Technologies (IDT, positive control). In vitro transcription of the unmodified PAN transcript was performed using the MEGAscript T7 Transcription Kit (Thermo Fisher AM13345). The 100% m6A modified full-length PAN transcript was synthesized using the MEGAscript T7 Transcription Kit following the manufacturer's protocol but substituting adenosine-triphosphate with N6-methyladenosine-5′-triphosphate (Jena Bioscience NU-1101S). RNAs were purified using the MEGAclear Transcription Clean-Up Kit (Thermo Fisher AM1908).

Probe generation

Ten biotinylated antisense probes, each 120 nt long, were designed to anneal to the PAN RNA sequence and synthesized per previously published method (Engreitz et al. 2013). One 45 nt long probe annealing outside the K7/PAN loci overlap was designed to perform the depletion of the K7 mRNA sequence (Arias et al. 2014). Three 50 nt long probes were purchased (IDT) and used to capture human U1 RNA (positive control). DNA templates for PAN-specific probes were obtained through reverse transcription (RT) of 100 ng total RNA from the lytically-induced BCBL-1 using oligonucleotides listed in Supplemental Table 1. Reverse transcription (RT) reactions (total volume of 20 µL) included 250 ng random hexamer primers (IDT 51-01-18-26), 10 mM dNTP, 1× First-strand buffer, 0.1 M DTT, 2 µL RNase OUT (40U final, Thermo Fisher 10777-019), and 1 µL of SuperScript II RT (200U final, Thermo Fisher 18064014). The left and right tag sequences were added during PCR amplification (Supplemental Table 1). The products were purified using the PureLink PCR purification kit (Thermo Fisher K310002) and used as templates for a second round of amplification to add the T7 promoter sequence in the antisense orientation. DNA products were purified and used for in vitro transcription. 10 µg of resulting RNA was used for RT using the biotinylated tag primer (Supplemental Table 1). The biotinylated probes were purified using the RNA Clean & Concentrator-5 purification kit (Zymo R1013).

PAN RNA affinity capture

Total RNA was extracted using TRIzol (Thermo Fisher 15596026), followed by DNase I treatment (Thermo Fisher AM1907), and purification using RNA Clean & Concentrator-5. K7 RNA was depleted by affinity capture with a biotinylated antisense K7 oligonucleotide probe (Supplemental Table 1). The remaining RNA fraction was purified using SILANE beads (Thermo Fisher 37002D) and fragmented using the COVARIS 8 sonicator with the following conditions: 250 sec at 75 peak power with 26.66% duty factor, 1000 cycles, with 20 average power at 4°C. PAN affinity capture was performed using 5 µg of total RNA, 15 pmol pool of the antisense probes (Supplemental Table 1), and streptavidin-coated magnetic beads (Thermo Fisher 65002).

Quantitative reverse transcription PCR (RT qPCR)

RT qPCR was performed using SuperScript II RT with random hexamers following the manufacturer's protocol. A 1:10 dilution of the RT reaction was used to perform quantitative real-time PCR (qPCR) assays using Perfecta SYBR Green FastMix Low ROX (Quanta BioSciences 95071). The BestKeeper ver. 1 was used to perform the repeated pairwise correlation and regression analysis for the level of RNA expression of housekeeping genes, that is, β-actin, glyceraldehyde 3-phosphate dehydrogenase (GADPH), U1 small nuclear (sn) RNA, and Metastasis Associated Lung Adenocarcinoma Transcript 1 (MALAT1) lncRNA during the latent and lytic stages of KSHV replication. As previously shown (Borah et al. 2011), the expression levels of β-actin and MALAT1 RNA were not affected by the progression of KSHV replication. The analysis of MALAT1 and β-actin transcripts expression resulted in the coefficient of correlation values of R = 0.98 and R = 0.94, respectively. Subsequently, these RNAs were used for data normalization (Supplemental Fig. 7b,c).

4-Selenothymidine-5′-triphosphate reverse transcription (4SedTTP RT)

A total of 2.5 µg of fragmented and affinity captured PAN RNA was used for the 20 µL ligation reaction that included 20 pmol of RiL-19 oligonucleotide (Supplemental Table 1), 2 µL 10× T4 ligation buffer, 1.8 µL 100% DMSO, 100 mM ATP, 8 µL PEG-8000, 0.3 µL of RNase Inhibitor (12U final, NEB M0307L) and 1.8 µL of T4 RNA ligase (20U final, NEB M0204S), and incubated at room temperature for 90 min. Reverse transcription reaction was initiated by annealing 20 pmol of AR-17 oligonucleotide (Supplemental Table 1) to the sample and incubating the reaction for 5 min at 75°C and slowly cooling to 25°C. Next, 2 µL of 4SedTTP reaction buffer (500 mM Tris-HCl pH 8.0, 500 mM KCl, 50 mM MgCl2, 100 mM DTT), 2 µL of 800 µM dATP, dCTP, dGTP, 4SedTTP (final concentration of 80 µM for each), 2 µL RNase Inhibitor (80U final, NEB M0307L) and 2 µL SuperScript III (400U final, Thermo Fisher 18080044) were added to each reaction to a final volume of 20 µL. The RT reactions were incubated at 42°C for 40 min, followed by incubation at 85°C for 15 min. The 100 nt long unmodified and m6A modified in vitro transcripts were treated in parallel with the experimental samples (Supplemental Table 1). Negative controls included fragmented, affinity captured PAN RNA (at 0–72 h pi) that was ligated to AR-17 oligonucleotide and directed to control RT reactions.

Deep-sequencing data deconvolution

The libraries were prepared by the RNA ligation method (Engreitz et al. 2015) and paired-end sequenced on Illumina HiSeq 2500 platform (Psomagen USA), resulting in ∼1 million reads per sample. The trimmed high-quality reads were mapped against the PAN sequence using Bowtie2 (Langmead and Salzberg 2012). Mapped reads for each library pair (+4SedTTP, −4SedTTP) were scaled to the data set with the lowest number of mapped reads. Two metrics were used to analyze 4SedTTP RT sequencing data: 1) m6ARATIO, which describes the proportion of reads supporting RT stops at a position out of all the reads overlapping it, as per Equation 1. where +4SedTTP corresponds to RT reactions including 4SedTTP, and −4SedTTP corresponds to RT reactions including dTTP; 2) m6A fold change (m6AFC), which is the log2 transformed m6A ratio in the treated samples (+4SedTTP) divided by the m6A ratio in the nontreated samples (−4SedTTP) (Schwartz et al. 2014a), as per Equation 2. The reported m6AFC value corresponds to the RT stop 1 nt downstream from the m6A position.

Direct RNA sequencing analysis

The 500 ng of total RNA extracted from BCBL-1 cells at the KSHV latent (0 HPI), and late lytic infection (48 HPI) were subjected to ribosomal RNA depletion using the RiboCop rRNA Depletion Kit (Lexogen 037). The direct RNA sequencing library preparation was performed according to the manufacturer's directions (ONT SQK-RNA002). Base-calling was performed using Albacore 1.2.1 [-f FLO-MIN106 -k SQK-RNA001 -r -n 0 -o fastq, fast5]. Reads present in the “pass” folder were used in the subsequent analysis. EpiNano software was utilized to determine the position of m6A (Liu et al. 2019). The SquiggleKit was used to view current intensities (Ferguson and Smith 2019); SAMtools and BCFtools (Li et al. 2009) were used for evaluating the number of insertions/deletions (INDELs), mismatches, and base quality scores.

PAN RNA cross-linking (UV and FA) followed by antisense purification with mass spectrometry (RAP-MS)

A total of 2 × 107 BCBL-1 cells per time point were cross-linked with 2% formaldehyde (VWR M134) in prewarmed 1× phosphate buffered saline (PBS) and incubated for 10 min at 37°C (Engreitz et al. 2013). About 2 × 107 BCBL-1 cells per time point were cross-linked with UV light at 254 nm wavelength with 8000 × 100 µJ/cm2 in 1× PBS (McHugh and Guttman 2018). These reactions were quenched with glycine to a final concentration of 500 mM and washed thrice with 1× PBS followed by PAN affinity capture. 60 pmol of antisense probes were added to 300 µL 1× hybridization buffer (20 mM Tris-HCl pH 7.5, 7 mM EDTA, 3 mM EGTA, 150 mM LiCl, 1% NP-40, 0.1% sodium deoxycholate, 3 M guanidine thiocyanate, 2.5 mM tris[2-carboxyethyl] phosphine [TCEP]), 5 µL protease inhibitor cocktail set (EDM Millipore 539134), and 10 µL murine RNase Inhibitor (200U, NEB M0314L) and incubated for 2 h at 37°C, after which they were washed six times with 1× PBS at 45°C. The elution of PAN-protein complexes was performed using 0.5 µL Benzonase nuclease (125U final, EDM Millipore 71206) in 1 mL of 0.35 M NaCl for 2 h at 37°C (Engreitz et al. 2013, 2014; McHugh and Guttman 2018).

MS data processing and normalization

Captured proteins were analyzed using the Ultra-High-Performance Liquid Chromatography directed to Q Exactive Hybrid Quadrupole-Orbitrap Mass Spectrometry (YPED Proteomics at Yale School of Medicine). The Orbitrap raw files were analyzed using MASCOT (ver. 2.6) and the UniProt canonical human protein database supplemented with common contaminants. Samples were searched, allowing for a fragment ion mass tolerance of 10 ppm and cysteine carbamidomethylation (static) and methionine oxidation (variable). A 1% false discovery rate for both peptides and proteins was applied. Up to two missed cleavages per peptide were allowed, and at least four peptides were required for protein identification and quantitation with a fragment mass tolerance of 0.02 Da and a monoisotopic weight. Proteins that were specifically captured as interacting with PAN by both formaldehyde (FA) and UV cross-linking and that were not identified in non-cross-linked samples were directed for further analysis. Data normalization was performed using the R package NOMAD (normalization of mass spectrometry data [Murie et al. 2018]) followed by PCA (principal component analysis) to provide a dimensionality-reduction solution to identify the proteins that have a high confidence (P > 0.05) for FA and UV cross-linking methods. The PCA analysis was performed utilizing two different R packages, prcomp (v.3.6.2) and EigenR (v1.0.0). We first generated the eigenvalues and eigenvectors of the covariance matrix related to FA (Principal Component 1) and UV (Principal Component 2) samples. PC1 (FA) was determined to be the most discriminating factor, and at least 60% of variance in the peak data and 57% of variance in the full spectrum data when using normalization methods. PC2 (UV) was the second most discriminating factor, where at least 45% of variance peak data and 44% of variance in full spectrum data.

Sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) and western blotting

A total of 2 × 107 BCBL-1 cells were washed with 1× PBS and resuspended in nuclear lysis buffer I (10 mM HEPES pH 7.4, 20 mM KCl, 1.5 mM MgCl2, 0.5 mM EDTA, 1 mM TCEP, 0.5 mM phenylmethylsulphonyl fluoride [PSMF] and 0.1% NP-40) supplemented with 10 µL of murine RNase Inhibitor (200U final), and 5 µL Protease Inhibitor Cocktail Set (EMD Millipore 539134) in a total volume of 1 mL. Cells were sonicated for 180 sec at 4°C, 15 peak power, 10 duty factor, 200 cycles at 1.5 average power with the wavelength adaptor no. 500534. DNase I (NEB M0303) treatment was performed according to the manufacturer's protocol, followed by 30 min 16,000g centrifugation at 4°C. The whole-cell protein lysates and the proteins cross-linked to PAN were analyzed by 10% SDS-PAGE. Proteins were transferred to a polyvinylidene fluoride (PDVF) membrane using a semidry blotting system (Thermo Fisher IB24002). The membranes were washed with 1× PBS with 0.1% Tween 20 (PBST) and blocked with 4% (w/v) nonfat milk in 1× PBST for 1 h at room temperature. The membranes were incubated with a specific primary antibody (Supplemental Table 2) for 1 h at room temperature, washed thrice with 1× PBST, and incubated with the horseradish peroxidase (HRP)-conjugated secondary antibody (Supplemental Table 2) for 45 min at room temperature. Blots were developed using the enhanced chemiluminescence substrate (Bio-Rad 1705061).

siRNA knockdowns

A total of 1 × 107 BCBL-1 cells were seeded in 30 mL fresh RPMI-1640 medium supplemented with 1% L-glutamine without antibiotics and without FBS to synchronize cells at G0 before transfection. For transfection, 80 pmol target-specific siRNAs (Santa Cruz Biotechnology, siMETTL3 sc-92172, siFTO sc-75002, siRBM15 sc-76365, and control siRNA-A sc-37007) were combined with 1 ml OPTI-MEM (Thermo Fisher 31985062) and incubated for 5 min at room temperature. Lipofectamine RNAiMAX (Thermo Fisher 13778) was combined with 1 ml OPTI-MEM and incubated for 5 min at room temperature. Both siRNAs/OPTI-MEM and Lipofectamine/OPTI-MEM solutions were combined, incubated for 5 min at room temperature, and added to the cells followed by 5 h incubation at 37°C and 5% CO2. Following the incubation, the second siRNA-directed knockdown procedure was performed as described above. Twenty-four hours after transfection, 10% FBS was added to induce the cell cycle reentry, followed by NaB treatment to the final concentration of 0.3 mM. The efficiency of siRNA knockdowns was verified by RT qPCR and western blotting.

Selenium-modified deoxythymidine triphosphates (SedTTP)-RT and ligation assisted PCR (SLAP)

A total of 2 × 107 BCBL-1 cells per time point were used for total RNA extraction. The 5′ end phosphorylation was performed by treating 10 µg total RNA with 0.5 µL of RNase Inhibitor (20U final, NEB M0307L), 1 µL of 10× T4 PNK reaction buffer, 0.1 mM ATP, 1 µL of 72 HPI polynucleotide kinase (10U final) in 10 µL total volume for 30 min at 37°C. To ligate the RNA-5 oligonucleotide (Supplemental Table 1), 1 µL of 10× T4 RNA Ligase Reaction Buffer, 100 pmol RNA-5 oligo, 0.1 mM ATP, 1 µL of RNase Inhibitor (10U final), 3 µL of 100% DMSO, and 1 µL of T4 RNA ligase I (10U final, NEB M0437M) were combined in 20 µL total volume and incubated for 16 h at 16°C. A total of 13.5 µL of the ligation reaction was combined with 10× annealing buffer (250 mM Tris-HCl, pH 7.4, 480 mM KCl) and 0.5 pmol target-specific RT primer in 15.5 µL, and the reaction was incubated for 2 min at room temperature. A total of 2 µL of 4SedTTP reaction buffer (500 mM Tris-HCl pH 8.0, 500 mM KCl, 50 mM MgCl2, 100 mM DTT), 2 µL of 800 µM dATP, dCTP, dGTP, 4SedTTP (final concentration of 80 µM for each), 2 µL RNase Inhibitor (80U final, NEB M0307L), and 1 µL AMV RT (12U final, Thermo Fisher 18080044) were added to each ligation reaction to the final volume of 20 µL and incubated for 1 h at 42°C, followed by the addition of 1 µL RNase H (5U final, NEB M0297L) and incubation for 20 min at 37°C. To anneal the adapters, 1.5 pmol adapter/splint oligos mixture (1.5 µM each) was added to the RT mixture, followed by incubation for 3 min at 75°C. To ligate the adapter, 4 µL of T4 DNA ligase (40U final, NEB M0202L), 1× T4 DNA ligase reaction buffer, and 2 µL 100% DMSO were added to the above mixture in a total volume of 25 µL and incubated for 16 h at 16°C. The reactions were purified using SILANE beads. PCR was performed with 0.2 µL Platinum Taq High Fidelity DNA Polymerase (2U final, Thermo Fisher 11304011) in 25 µL total, and the products were analyzed on 10% polyacrylamide gel electrophoresis (PAGE) and stained with 1× ethidium bromide solution overnight at 4°C (Zhang et al. 2019). Densitometric analysis of the PAGE gels was performed using ImageJ v1.52a (Abràmoff et al. 2006) on inverted TIFF files.

Selective 2-hydroxyl acylation analyzed by primer extension and mutational profiling (SHAPE-MaP) of PAN RNA

Each reaction included 1 × 107 BCBL-1 cells washed with 1× PBS, collected by centrifugation at 300g for 5 min, and resuspended in 900 µL of RPMI-1640 medium. The SHAPE (+) reaction included the addition of 100 µL of 100 mM 1-methyl-7-nitroisatoic anhydride (1M7, DC Chemicals DC8649) resuspended in 100% DMSO to a final concentration of 10 mM to 900 µL of cells, followed by the incubation for 5 min at 37°C. For the SHAPE (−) reaction, 100 µL of 100% DMSO was added to 900 µL of cells and incubated for 5 min at 37°C. The denatured control (dc) reaction included 1 µg/µL of total RNA extracted from BCBL-1 cells, which was combined with 5 µL 100% formamide and 1 μL 10× denaturing buffer (500 mM HEPES pH 8.0, 40 mM EDTA) in a total volume of 9 μL and incubated for 1 min at 95°C. An amount of 1 µL of 100 mM 1M7 was subsequently added to 9 µL of the sample and incubated for 1 min at 95°C. SHAPE-MaP libraries were constructed by dividing the PAN RNA sequence into three overlapping zones, including nt 15–343 for zone 1, nt 311–707 for zone 2, and nt 694–1060 for zone 2 (Supplemental Table 1). Libraries were sequenced using MiSeq Nano V2 (Illumina) paired-end sequencing 2 × 250 bp with on average 100,000 reads per sample. SHAPE-MaP data were analyzed using the ShapeMapper v1.2 pipeline (21).

Fluorescence in situ hybridization, costaining, and image quantification

A total of 1 × 107 BCBL1 cells were collected at the indicated time points post-induction, fixed with 4% paraformaldehyde (Electron Microscopy Sciences 15852-15s), permeabilized with 0.5% Triton X, and washed thrice with 1× PBS. PAN specific oligonucleotide probes were labeled with a DIG-Oligonucleotide Tailing Kit (Sigma Roche, 2nd generation 03353583910) and mixed at the 1:1:1 ratio. Probes were denatured for 5 min at 95°C before their addition to the hybridization solution (10% w/v dextran sulfate, Sigma D-6001), 10% v/v deionized formamide (VWR 97062-008), 2× nuclease-free saline sodium citrate, followed by overnight incubation at 37°C. After 1 h blocking, primary mouse anti-DIG (Sigma Roche 11333062910) was added at a dilution of 1:200 in 1% BSA and incubated for 1.5 h at room temperature. After three complete washes, the slides were incubated with a secondary Alexa Fluor 488-conjugated goat anti-mouse antibody (Thermo Fisher A-11029) at a dilution of 1:500 in 1% BSA. Methylome detection including METTL3, RBM15, FTO, YTHDF2, HNRNPC, SND1 was performed by using specific primary anti-rabbit antibodies in 1% BSA (Supplemental Table 2). After three complete washes, the secondary Alexa Fluor 647-conjugated goat anti-rabbit antibody (Thermo Fisher A-21244) was added at a dilution of 1:500 in 1% BSA. An amount of 1 µg/mL DAPI (Thermo Fisher 62248) was used to stain nuclei for 30 min at 37°C. The cells were mounted using VECTASHIELD Antifade Mounting Medium H-1200 (Vector Laboratories) and observed under Nikon ECLIPSE Ti confocal microscope. Image J 1.52a was used to analyze FISH and IF images. To determine the raw fluorescence intensity in the stacked images, the embedded line tool and Plot Profile function were selected in Image J. The line traces quantitatively representing the raw fluorescence signals across all channels and data for the different stains were exported to Microsoft Excel, normalized, and plotted. PCC was calculated as a thresholded coefficient in Volocity 5.3 using the JACoP plugin (http://imagejdocu.tudor.lu/doku.php?id_plugin:analysis:jacop_2.0:just_another_colocalization_plugin:start) by ignoring pixel pairs with a gray level beneath the threshold of either image (Barlow et al. 2010). PCC was then calculated only from those pixels that are above the threshold in both channels (Barlow et al. 2010). These coefficients are referred to as thresholded PCC.

Motif analysis

The identification of m6A consensus motifs within the PAN RNA sequence and cellular lncRNA sequences was performed with Multiple EM for Motif Elicitation (MEME, v5.1.1) using 8 nt windows surrounding the m6A site.

RiboSNitch analysis

The analysis of PAN RNA structural changes invoked by the ablation of METTL3 and FTO was performed by using SNPfold (Halvorsen et al. 2010) with set default parameters, followed by the variant analysis of global and local structural change by classSNitch (Woods and Laederach 2017).

Data visualization

Graphics were created using Microsoft Excel, Biorender (biorender.com), BEG Venn Diagram (bioinformatics.psb.ugent.be/webtools/Venn), Adobe Photoshop (ver 22.2), ImageJ (ver. 1.51), Geneious Prime 11.1.5 (https://www.geneious.com), Cytoscape ver. 3.8.2, and R studio ver. 3.6.2 (R Core Team 2018).

DATA DEPOSITION

Data supporting reported results can be found at NCBI Sequence Read Archive (SRA) submission numbers SUB9396779, SUB9396663, and SUB9421156. All data are reported in the text and Supplemental Material. At home-developed R scripts are available at https://github.com/jzs0165. The PAN RNA sequence can be found at accession number U50139.1.

SUPPLEMENTAL MATERIAL

Supplemental material is available for this article.