We performed Illumina sequencing starting with total RNA extracted from VeroE6 cells and from SARS-CoV-2 viral preparations cultured on VeroE6 cells (Figure 1A; Supplementary Table S1). To improve efficiency and quantitative assessment of small RNA-seq, defined here as RNA of <200 nucleotides in length, we built two libraries for each sample. The first library was treated with a demethylase mixture (DM) (Zheng et al., 2015) which removed many Watson-Crick face methylations in tRNA that impede reverse transcription in library construction while the second library was left untreated. As described previously, the DM-treated libraries are useful for quantitative assessments of transcript abundance, whereas the untreated samples are useful for modification analysis (Clark et al., 2016). As expected, sequencing reads of the VeroE6 cells mostly mapped to tRNAs, followed by those from 5S to 5.8S rRNA, a small amount of SRP RNA, and others such as spliceosomal RNA (snRNA) and Y RNA (Figure 1B). In the viral preparations, a substantial proportion of reads mapped to viral genomic RNA as expected, and tRNA and SRP RNA are present at almost equally high proportions, followed by a small amount of rRNA (Figure 1B). Although one cannot exclude RNA in exosomes or extracellular RNA not associated with vesicles, our data clearly show RNAs that are differently present in our cell free preparations as compared to those in cells. Our viral preparations used for the sequencing experiments contain high levels of infectious particles outside the cell, strongly suggesting that the sequenced RNA are derived, to a large extent, from cell free virions. For example, our results show a ∼150-fold enrichment of the SRP RNA over tRNA in the viral preparation samples vs. the cell samples, which suggests that we eliminated most if not all of the cellular debris. These results indicate that SARS-CoV-2 virions also package tRNA and SRP RNA in significant proportions.

We next examined the selectivity of packaged tRNA at two levels. Mammalian genomes contain many tRNA isodecoder genes that share the anticodon, but possess different body sequences; all tRNA isodecoders with the same anticodon belong to a single tRNA isoacceptor family (Goodenbour and Pan, 2006; Schimmel, 2018). Abundance of isodecoders was summed for each isoacceptor family and used to calculate the fraction of tRNA reads for each anticodon. The isoacceptor abundance fraction for each sample (three biological replicates of uninfected VeroE6 cells and six distinct viral culture supernatant preparations), was compared to the mean of Vero E6 cells. All three VeroE6 cell replicates were nearly identical, as the heat map shows close to zero values in all cases (Figure 1C). To avoid exaggerated representation of low abundant tRNAs by ratioed comparison, we subtracted the tRNA fraction in each viral preparation to its counterpart in the VeroE6 cells, so that the differences were readily identified for more abundant tRNAs (Figure 1C; Supplementary Table S2). We found several isoacceptor families that are significantly enriched across all six isolates. They include tRNA^Glu(TTC), tRNA^Lys(TTT), tRNA^Leu(AAG), tRNA^Ser(AGA), tRNA^Ser(GCT), and tRNA^Ser(TGA). These results indicate that SARS-CoV2 virions selectively incorporate tRNA isoacceptors.

Our downstream analysis, thus, focused on those six tRNAs enriched in the viral preparations. First, we analyzed the tRNA at the isodecoder level for all six tRNAs. Among the seven tRNA^Glu(TTC) isodecoders, four could be detected in the viral preparations. However, only two isodecoders represent almost all tRNA^Glu(TTC) in the viral preparations, even though neither is the most abundant isodecoder in VeroE6 cells (Figure 2A). In contrast, the single dominant tRNA^Leu(AAG), tRNA^Lys(TTT), tRNA^Ser(AGA) isodecoders in cells are also the ones in the viral preparations (Figures 2B–D). For tRNA^Ser(GCT) and tRNA^Ser(TGA), two isodecoders each are present at appreciable levels, in each case, the isodecoder at the highest level is also the one in the viral preparations (Figures 2E,F).

We examined the read pileup of the most abundant isodecoder in the viral preparations. By experimental design, our tRNA-seq results always start from the 3′ end of the tRNA and show a decline toward the 5′ end with sharp drops at certain tRNA modifications, an expected behavior for full-length tRNAs (Zheng et al., 2015). Three types of results are observed: first, the pileup decreases faster in the viral preparation tRNA compared to the cellular tRNA, this group includes tRNA^Glu(TTC) (Figure 3A). The pronounced drop of the viral preparation tRNA in the anticodon loop region is consistent with GluTTT-75 in the virion being a 3′ half tRNA fragment with the 5′ end in the anticodon loop, because there is no known RT stopping modifications in this tRNA. In the second type, the read pileup decreases at about the same rate, this group includes tRNA^Leu(AAG) and tRNA^Lys(TTT) (Figures 3B,C). The similar drop off is consistent with the tRNA in the viral preparations as the full-length tRNA like those in cells, and the sharp drop offs corresponds to the N2,2-dimethyl-G at position 26 (m² ₂G26) in tRNA^Leu(AAG) which is difficult to remove by the demethylase because it is buried in the tRNA structure (Dai et al., 2017) and 2-methylthio-6-carbamoylthreonine at position 37 (ms²t⁶A37) in tRNA^Lys(TTT) (Machnicka et al., 2013) which does not react with the demethylase. In the third type, the read pileup decreases slower in the viral preparation tRNA, this group includes all three tRNA^Ser, and the sharp drop offs correspond to the m² ₂G26 modification (Figures 3D–F). This result is consistent with the tRNA^Ser in the viral preparations having lower modification levels in the anticodon stem-loop region which can include N6-methyl-N6-threonylcarbamoyladenosine at position 37 (m⁶t⁶A37) in tRNA^Ser (Machnicka et al., 2013).

We next looked for RNA modification differences among the tRNAs from cells and viral preparations. In tRNA-seq, certain modifications can be identified by their “mutation” signatures in the sequencing data. Vero cells are derived from Chlorocebus sabaeus kidney, their tRNA modifications have not been reported in the literature. We analyzed the mutation signatures of the cellular tRNAs in the libraries without demethylase treatment and validated the methylations using the results from demethylase treatment (Clark et al., 2016) to provide a comprehensive analysis of Vero cell tRNA modifications (Table 2). Overall, the C. sabaeus tRNA modification patterns detected by sequencing are very similar to those from the human HEK293T cells (Table 1 in Clark et al. (2016)). A minor difference is the m² ₂G26 modification which is present in tRNA^Val in C. sabaeus Vero cells but not in HEK293T cells. We detected inosine modification at the wobble anticodon position (I34) in all tRNAs that are A34 in the genome which include tRNA^Leu(AAG) and tRNA^Ser(AGA) (Figure 4A; Supplementary Figure S1A; Table 2).

We then compared the mutation levels between the tRNAs enriched in viral preparations and their counterparts in cells. For high confidence analysis we applied a filter of ≥50 read coverage at nucleotide positions of interest. The modification sites passing this filter among the tRNAs in viral preparations only include m¹A58 (tRNA nomenclature) in tRNA^Leu(AAG), tRNA^Lys(TTT), and tRNA^Glu(TTC), and I34 in tRNA^Leu(AAG) and tRNA^Ser(AGA). For tRNA^Leu(AAG) and tRNA^Lys(TTT), the mutation fraction at m¹A58 is higher in the tRNA from the viral preparations than the VeroE6 cell tRNA (Figure 4B; Supplementary Figure S1B), suggesting preferential packaging of m¹A modified tRNA. Among the tRNA^Glu(TTC) isodecoders, m¹A58 level is variable in cells. Only two of the five abundant tRNA^Glu(TTC) isodecoders have high modification levels in VeroE6 cells, but only isodecoders with low modification fractions are present in the viral preparations (Figure 4C). Since tRNA^Glu(TTC) in the virions are likely tRNA fragments, this result is consistent with low m¹A modified tRNA^Glu(TTC) being the preferred source of tRNA^Glu(TTC) fragments in cells. tRNAs in cells and in the viral preparations are >90% modified with I34 in both tRNA^Leu(AGA) and tRNA^Ser(AGA) (Figure 4A; Supplementary Figure S1A).

We next carried out large RNA-seq (>200 nt) of the viral preparations to characterize the viral RNA and its candidate modifications. We first removed small RNAs by size-selection, followed by chemical fragmentation and library construction. As expected, most of the reads mapped to the SARS-CoV-2 genome (Wuhan reference). We used a mutation threshold of >90% to identify single nucleotide polymorphisms (SNPs) in these samples (Figure 5A). These SNPs did not change upon our enzymatic or chemical treatment described below (data not shown). Our samples are derived from distinct viral isolates obtained from patients infected during the first pandemic wave in the spring of 2020, at the time when most of the circulating SARS-CoV-2 viruses still had only a few sequence changes (Gonzalez-Reiche et al., 2020). We measured the read counts of the SARS-CoV-2 genomic RNA and the 18S + 28S rRNA in our viral preparations (Figure 5B). We found an average ratio of SARS-CoV-2/rRNA of ∼9.5, or a molar ratio of SARS-CoV-2/rRNA of ∼2. Given that the cellular ribosomes and viral particles produced upon infection (Sender et al., 2021) is higher than 100:1, our results show a >200-fold enrichment of SARS-CoV-2 viral RNA over rRNA in our culture supernatant preparations.

We also measured the ratio of reads from the SRP RNA, which was the most abundant RNA in the viral preparations in the small RNA-seq data. We counted the read coverage for the six SARS-CoV-2 genomes and those mapping to the C. sabaeus SRP from the Rfam database (Kalvari et al., 2021). We then used a normalization factor that accounts for the length difference of SARS-CoV-2 genome (29,903 nt) and SRP RNA (300 nt) to get a ratio of SRP to SARS-CoV-2 RNA (Figure 5C). We obtained a ratio of .3–.55 for the six viral preparations. Given the intermediate size of the SRP RNA that may incur losses during the size-selection step, we estimate that a virion contains up to one SRP RNA transcript.

We also asked whether our cell free viral preparations contain SARS-CoV-2 subgenomic RNA. Subgenomic RNAs are generated in the infected cells by joining a 5′ leader sequence to each of the SARS-CoV-2 protein coding genes (Supplementary Figure S2). We found 2,000–3,000 reads that cover the junction region of the 5′ leader with Orf1a, which corresponds to the full-length SARS-CoV-2 RNA. Normalizing the reads covering other genes to this junction, we found up to 2% of subgenomic RNAs containing S, Orf3a, Orf7a, and E genes in all six isolates (Figure 5D). Altogether, up to 10% of the viral RNAs in the cell free viral preparations are subgenomic RNAs. This low level is consistent with our viral preparations comprising primarily of SARS-CoV-2 virions, and most virions in the culture containing the full-length genomic RNA. The subgenomic RNAs may represent in part virions devoid of the full-length genomic RNA (i.e., defective viral particles).

Finally, we carried out experiments to identify novel modifications of SARS-CoV-2 vRNA. To facilitate this identification, we added two consecutive steps in the large RNA library construction, one with demethylase treatment (DM), and the other with 1-cyclohexyl-(2-morpholinoethyl)carbodiimide (CMC) treatment (Figure 6A). Demethylase treatment generally detects Watson-Crick methylations such as m¹A in our tRNA studies, whereas CMC treatment is a standard procedure for detecting pseudouridine (Ψ) modification in mRNA (Gilbert et al., 2016). We first carried out a threshold analysis in four pairwise comparisons of ±DM, ±CMC, ±DM (plus CMC treatment), ±CMC (plus DM treatment) using the filters for mutation difference of 5%, stop difference of 15%, deletion difference of 2%, and insertion difference of 2%. We filtered next the candidate sites from the threshold analysis to only include those that show the same signature change in at least four of the six viral preparations. Only the mutation signature of five sites passed these two filters (Figure 6B). These sites fall into four groups: 1) U8323 and U20331 are DM-sensitive without CMC, CMC-sensitive without DM, not detected in DM with CMC, nor in CMC with DM. These two sites may represent N3-methyl-U (m³U) derivatives. 2) A29517 is DM-sensitive without CMC, CMC-sensitive without DM, the mutation signals are reversed in DM with CMC, or in CMC with DM. This site may represent N1-methyl-A (m¹A) derivatives. 3) U3877 is DM-insensitive without CMC, CMC-sensitive without DM, DM-sensitive with CMC, not detected in CMC with DM. This site may represent N1-methyl or N3-methyl-pseudouridine (m¹Ψ, m³Ψ) derivatives. 4) A29780 is DM-insensitive without CMC, CMC-sensitive without DM, DM-sensitive with CMC, not detected in CMC with DM. We do not recognize an A modification at this time that would generate such signatures. No stop signature was observed for any of these five sites (Supplementary Figure S3).

It is intriguing to note that these five candidate modification sites in SARS-CoV-2 do not follow the expected pattern of changes for well-characterized RNA modifications. In all five cases, the mutation levels are at most 20%, which may reflect either low levels of modification at these sites or under-counting the modification levels because of the unknown nature of these modifications. Future studies are needed to pinpoint the exact nature of these modifications.

Residual nasopharyngeal swab specimens were collected after completion of the diagnostic process as part of the Mount Sinai Pathogen Surveillance Program. To culture SARS-CoV-2 isolates, .1 ml of viral transport media was inoculated into one well of a six-well plate seeded with a confluent monolayer of VeroE6 cells. Culture supernatants were harvested when CPE (cytopathic effect) became visible, aliquotted, and stored at −80°C. All work related to SARS-CoV-2 culture was performed in a BSL3 biocontainment facility by trained personnel and in accordance with the research registration approved by the Institutional Biosafety Committee (IBC).

We cultured six distinct SARS-CoV-2 isolates representing the early lineages of the pandemic (Table 1). After isolation of the clinical isolates on VeroE6 cells, we determined the infectious viral titers for each of the viral culture supernatants by plaque assay (Table 1). All six isolates displayed medium plaque phenotypes. We shared aliquots of the viral stocks analyzed in this study with the NIH BEI repository in the early summer of 2020 (see Table 1 for specifics).

The viral culture supernatants were spun at 3,000 rpm for 10 min to remove particulates. vRNA from each viral preparation and total RNA from VeroE6 cells was performed using QIAamp Viral RNA Kits (Qiagen) following the manufacturer’s instructions. We did not filter the supernatants as that could result in some viruses being absorbed to the membrane and/or causing a loss of viral infectivity.

Reference RNA sequences from Chlorocebus sabaeus that included non-coding RNA and tRNA were downloaded from Rfam database (https://rfam.xfam.org/, Kalvari et al., 2021). C. sabeaus tRNA sequences from Rfam were processed through tRNAScan-SE (http://lowelab.ucsc.edu/tRNAscan-SE/, Lowe and Chan, 2016), only sequences with high confidence (i.e., tRNAScan score ≥50) were used as reference. Following this, tRNA sequences were appended by adding CCA at their 3′ end as well as removing intron sequences (Supplementary Table S1). These processed C. sabeaus cytosolic tRNA, mitochondrial tRNA, and non-coding RNA sequences such as SRP from Rfam were combined with the Wuhan SARS-CoV-2 genome sequence (MN908947.3) to generate a custom reference database.

Raw reads following sequencing were designated reads 1 and reads 2 and were merged together using bbmerge.sh present within the bbmap package (https://github.com/BioInfoTools/BBMap), which results in merged fastq files. These merged fastq were converted to fasta file format using reformat.sh present within the bbmap package. These fasta files were aligned to our custom reference genome using bowtie2 (http://bowtie-bio.sourceforge.net/bowtie2/index.shtml, Langmead et al., 2009) with the following parameters: -f -p 10 --local —no-unal. The aligned reads were then used to determine RNA sequence abundance using custom python script. RNA modifications were detected based on aligned reads using samtools sort (http://www.htslib.org/, Li et al., 2009) feature sort the reads in a bam file format. Then IGVtools count (https://software.broadinstitute.org/software/igv/igvtools) count feature was utilized to output a wig files using the following parameters: -z 5 -w 1 -e 250 —bases. The resulting wig files were processed using a custom python script to identify nucleotide mutations as well as coverage of aligned reads.