Changes of Small Non-coding RNAs by Severe Acute Respiratory Syndrome Coronavirus 2 Infection

Introduction

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a beta coronavirus belonging to the sarbecovirus subgenus of Coronaviridae family (Zhu et al., 2020). It is a positive-sense single-stranded RNA virus with a genome length of ∼30 kb. By the middle of January 2022, the ongoing coronavirus disease 2019 (COVID-19) pandemic, caused by SARS-CoV-2, has respectively caused more than 320 million infectious cases and over five million deaths globally (World Health, 2021).

Small non-coding RNAs (sncRNAs) have diverse functions through various regulatory mechanisms. They virtually participate in all biological pathways, including cell proliferation, differentiation, apoptosis, autophagy, and tissue remodeling. sncRNAs are also essential to regulate host responses to viral infections (Choudhuri, 2010; Beermann et al., 2016; Romano et al., 2017; Rajput et al., 2020; Wu et al., 2020). Among sncRNAs, the most widely studied sncRNAs are microRNAs (miRNAs), which are 18–24 nt in length, carry 5′ monophosphate and 3′ hydroxyl (3′-OH) ends, and generally regulate genes via the argonaute (AGO) platform (Schwarz et al., 2004; Fabian and Sonenberg, 2012).

Other than miRNAs, piwi-interacting RNAs (piRNAs), small nucleolar RNAs (snoRNAs), and tRNA-derived RNA fragments (tRFs) are also important members of sncRNAs (Dozmorov et al., 2013). Currently, there is very limited information on whether or how SARS-CoV-2 regulates the sncRNA expression, except the reports on SARS-CoV-2-impacted miRNAs (Mallick et al., 2009; Hasan et al., 2014; Khan et al., 2020).

Using T4 polynucleotide kinase (T4 PNK)‐RNA‐seq, a modified next-generation sequencing (NGS), we found that tRFs and piRNAs were the two most abundant sncRNAs in nasopharyngeal swabs (NPS) samples of the SARS-CoV-2-positive group. However, only tRFs were significantly enhanced in SARS-CoV positive samples. Generally, tRFs are generated by specific cleavages within pre-tRNAs or mature tRNAs (Lee et al., 2009). Compared with other sncRNAs, tRFs are relatively new members. However, their importance in diseases, such as cancer, infectious diseases, neurodegenerative diseases, and metabolic diseases, was quickly acknowledged after the discovery (Wang et al., 2013; Selitsky et al., 2015; Shen et al., 2018; Sun et al., 2018; Zhu et al., 2018; Choi et al., 2020; Qin et al., 2020; Wu et al., 2021). tRFs are classified mainly into tRF-1 series, tRF-3 series, and tRF-5 series (Fu et al., 2009). tRF-1 series are usually those from the 3′-trailer sequences of pre-tRNA, while tRF-3 and tRF-5 series are aligned to the 3′- and 5′- end of the mature tRNAs respectively. Among SARS-CoV-2-impacted tRFs, the most impacted tRFs belonged to tRF5s. In addition, the impacted tRF profile seemed to be SARS-CoV-2 specific, which is consistent with what we and others found previously on the changes in tRF signatures being virus-dependent (Wang et al., 2013; Selitsky et al., 2015), implicating tRFs as potential prognosis and diagnosis biomarkers. The impacted tRFs were also observed in SARS-CoV-2 infected human alveolar type II-like epithelial cells expressing human angiotensin-converting enzyme 2 (A549-ACE2) and human small airway epithelial cells (SAECs) in the air-liquid interface (ALI) culture.

In addition to host-derived ncRNAs, viral genomes can also encode ncRNAs. These viral ncRNAs vary in length and have diverse biological functions, including the regulation of viral replication, viral persistence, host immune evasion, host inflammatory response, and cell transformation (Tycowski et al., 2015). For example, SARS-CoV-encoded small RNAs contribute to SARS-CoV-induced lung injury (Morales et al., 2017), and SARS-CoV-2-encoded miRNAs enhance inflammation (Cheng et al., 2021). In this study, we revealed several new small viral RNA (svRNA) fragments, with the length of 25 nt, 33 nt, and 36 nt, by T4 PNK‐RNA‐seq. Among svRNAs derived from CoV-2 (sv-CoV2), a svRNA spanning from site 346 to site 382 of nsp1(sv-CoV2-346) had the highest expression.

In summary, this is the first report demonstrating the altered tRFs by SARS-CoV-2. T4 PNK pretreatment also enabled small RNA seq to reveal additional new sv-CoV2. In the future, we will characterize the biogenesis and function mechanisms of these new sncRNAs associated with SARS-CoV-2 infection.

Materials and Methods

Nasopharyngeal Swab Specimens

NPS were collected from patients who visited outpatient clinics of the University of Texas Medical Branch (UTMB) for SARS-CoV-2 screening in April 2020. NPS samples in universal viral transport media were transported to the Molecular Pathology laboratory, directed by Dr. Jianli Dong, and subjected to SARS-CoV-2 test using Abbott m2000 SARS-CoV-2 RT-PCR assay. The limit of detection (LOD) of detection assays is 100 viral genome copies/ml.

Thirteen anonymous NPS samples were used in this study, including seven SARS-CoV-2 negative (51.7 ± 13.7 years old) and six SARS-CoV-2 positive (49.2 ± 10.5 years old) samples. The protocol was approved by the Institutional Review Boards (IRB) of UTMB at Galveston, under the IRB protocol # 02-089 and 03-385.

RNA Isolation

After the SARS-CoV-2 validation, 1 ml of NPS sample from each individual was subjected to RNA extraction using the mirVana PARIS kit (Invitrogen, MA, United States) according to the manufacturer’s protocol. At the elution step, samples were incubated on the column for 5 min at 65°C, and the RNA was eluted with 45 µL nuclease-free water. To extract RNAs from cells, TRIzol reagents (Thermo Fisher Scientific, MA, United States) were used for total RNA preparation, as described (Choi et al., 2020), followed by qRT-PCR.

T4 PNK-RNA-seq and Data Analyses

To study whether other sncRNAs than miRNAs are impacted by SARS-CoV-2, we used T4 PNK-RNA-seq, a modified NGS, to get sncRNA profiles for samples derived from NSP or cultured cells, similarly as described in (Honda et al., 2015; Giraldez et al., 2019). A flowchart of the T4 PNK-RNA-seq is shown in Supplementary Figure S1A and data have been deposited in GEO (GSE193555). Basically, we treated sample RNAs with T4 PNK before the library construction and small RNA-seq to make the 3′-end of RNAs homogenous with -OH, as the ligation of the 3′-end of sncRNAs with sequencing barcodes requires the presence of 3′-OH and not all sncRNAs have 3-OH ends. The seq was done in the NGS Core of UTMB. In brief, the RNA samples were pretreated with 10 units of T4 PNK using 14 µL extracted RNAs in a final reaction volume of 50 µL and incubated at 37°C for 30 min, and then were heat-inactivated at 65°C for 20 min. The RNA was purified and concentrated within 6 µL nuclease-free water using Zymo RNA Clean and Concentrator kit (Irvine, CA, United States). Ligation-based small RNA libraries were prepared with an RNA input of 6 µL using NEB Next Multiplex Small RNA Library Prep Set for Illumina (Ipswich, MA, United States). Libraries were sequenced using the Illumina NextSeq 550 Mid-Output sequencing run. About 7,680 Mb of sequence data was generated.

To analyze the seq data, adaptor sequences were first removed using Cutadpat and reads with a length of more than 15 bp were extracted. We further filtered out RNAs with counts of less than 10 and all rRNA sequences, using the remainders as cleaned input reads. In terms of the mapping databases, we prepared tRF5 and tRF3 databases using the same sequences derived from different tRNAs [sequences downloaded from tRNA genes using the Table Browser of the UCSC genome browser (Karolchik et al., 2004)]. We also prepared tRF1 sequences using genome locations of tRNAs. Our in-house small RNA database includes 1) these tRFs, 2) miR/snoR sequences downloaded from the UCSC genome browser, and 3) piRNA sequences downloaded from piRBase (http://www.regulatoryrna.org/database/piRNA/). The cleaned input reads were mapped to our inhouse small RNA database using bowtie2 (v2.4.1) allowing two mismatches (option N -1). After we mapped the cleaned input reads to the small RNA database, the unmapped sequences were then mapped to the hg38 genome using the bowtie2 pre-built index (GRCh38_noalt_as) to detect all human sequences. The unmapped sequences to the human genome were then mapped to the SARS-CoV-2 reference genome (NC_045512) using the same parameters.

Raw read counts were normalized with the DEseq2 median of ratios method. Differentially expressed genes were determined by p-value < 0.05, fold change >2, and mean of normalized counts >10 in either Control (CN) or SARS-CoV-2 group. Unsupervised hierarchical clustering was performed using the Pearson correlation coefficient. A flowchart of the sequencing data analyses is summarized in Supplementary Figure S1B.

Cell Culture and Viruses

African green monkey kidney epithelial cells (Vero E6) were obtained from ATCC (Manassas, VA, United States) and maintained in a high-glucose DMEM (Gibco, MA, United States) supplemented with 10% fetal bovine serum (FBS), 10 units/ml penicillin, and 10 μg/ml streptomycin. The human alveolar type II-like epithelial cells expressing human angiotensin-converting enzyme 2 (A549-ACE2) cells were a kind gift from Dr. Shinji Makino and were cultured in DMEM (Gibco, MA, United States) containing 10% FBS, 10 units/ml penicillin, and 10 μg/ml streptomycin.

Small airway epithelial cells (SAECs), isolated from the normal distal portion of the lung in the 1 mm bronchiole area, were purchased from Lonza (Basel, Switzerland) to generate cells in the air-liquid interface (ALI) culture. The cells were cultured and differentiated using Complete PneumaCult™-Ex plus medium and PneumaCult™-ALI-S Maintenance medium (Stemcell Technologies, Vancouver, Canada), respectively, according to the manufacturer’s instructions.

Briefly, the cells at passage two (P2) were expanded in the T-25 flask using the complete PneumaCultTM-Ex plus medium, with a medium change every other day. For ALI cultures, the cells (P3) were seeded into Corning Costar 12 mm transwell inserts (Corning, NY, United States) at a concentration of 11 × 10⁴ cells/insert in 0.5 ml medium/insert, and another 1 ml/well medium was added to the basal chamber. Cells were submerged cultured in Complete PneumaCultTM-Ex plus medium, with a medium change every other day. After reaching ∼100% confluency, ALI was initiated by removing the apical medium and replacing the PneumaCultTM-Ex plus medium in the basal compartment with PneumaCult™-ALI-S Maintenance medium. The basal compartment medium was changed every other day. It took about 21 days to complete cell differentiation.

SARS-CoV-2 (USA-WA1/2020 strain) was obtained from the World Reference Center for Emerging Viruses and Arboviruses (WRCEVA) at the UTMB. Viral stocks were prepared by propagation in Vero E6 cells. Viral titers were determined by plaque assay as described (Blanco-Melo et al., 2020). All experiments using live SARS-CoV-2 were performed in a biosafety level 3 (BSL-3) laboratory with redundant fans in the biosafety cabinets. All personnel wore powered air-purifying respirators (Breathe Easy, 3M) with Tyvek suits, aprons, booties, and double gloves. All cell cultures, cell lines or primary cultured cells, and viruses have been approved for use by the Institutional Biosafety Committee of UTMB (NOU# 2018056 and NOU# 2020043).

Viral Infection

To infect A549-ACE2 cells in monolayer culture, the cells were seeded into the 24-well plate 24 h prior to the infection to allow the cells to reach 80–90% confluence in the following day. For infection, the cells were incubated with viruses in DMEM media with 10% FBS at a multiplicity of infection of 0.1 (MOI = 0.1). After 1 h incubation, cells were washed with PBS three times to remove the remaining viruses and cultured in fresh media containing 10% FBS. The cells were collected on day 4 post-infection.

Regarding the infection of SAECs in ALI culture, the infection was performed when hallmarks of excellent differentiation were evident, such as extensive apical coverage with cilia. Prior to infection, the apical side of the cells was washed five times with PBS, and the basal surface was washed once with PBS. Viruses were diluted to the specified MOI in 200 µL MEM medium and inoculated onto the apical surface of the ALI cultures. After a 2-h incubation at 37°C with 5% CO₂, unbound viruses were removed by washing the surface with PBS three times. The cells were collected on day 1 or 3 post-infection. The SARS-CoV-2 S gene was detected using qRT-PCR with primers as follows: S forward primer, 5′ CCT​ACT​AAA​TTA​AAT​GAT​CTC​TGC​TTT​ACT; reverse primer, 5′’ CAA​GCT​ATA​ACG​CAG​CCT​GTA.

qRT-PCR and RT-PCR

To evaluate sncRNAs expression, qRT-PCR was performed, as described previously (Choi et al., 2020; Wu et al., 2021). A schematic summary of tRF quantification by qRT-PCR is shown in the left panel of Figure 1A. Briefly, the total RNA was treated with T4 PNK, and then ligated to a 3′-RNA linker using T4 RNA ligase from Thermo Fisher Scientific (Waltham, MA, United States). The product was used as a template for reverse transcription (RT) with a linker-specific reverse primer using TaqMan Reverse Transcription Reagents from Thermo Fisher Scientific. The cDNA was subjected to SYBR Green qPCR using iTaq™ Universal SYBR Green Supermix kit from Bio-Rad (Hercules, CA, United States) and primers specific to the 5′-end of tRFs and RNA linker. U6 was used for normalization. The addition of a 3′-RNA linker enables the detection of tRF5s without the signal interference from its corresponding parent tRNAs, possibly because 1) the 3-end of tRNA is usually attached with an amino acid, preventing RNA linker attachment, and 2) reverse transcription annealing temperature sets tRFs, not tRNAs, to be annealed by the primer, as tRNA cloverleaf structure requires a specific denaturing temperature before annealing (right panel of Figure 1A) (Choi et al., 2020; Wu et al., 2021). The primers and 3′-RNA linker sequences are listed in Table 1.

FIGURE 1. The schematic summaries on tRF quantification by qRT-PCR (A) and the detection of sv-CoV2-346 by RT-PCR (B). (C) The sequence information on sv-CoV2-346 and associated primers and the 3′ RNA linker for its detection.

TABLE 1. Sequence information for tRF5s, RNA linker, RT primer, qPCR primers and PCR primers.

To validate the seq data of CoV2-encoded small RNAs (CoV2-346), RT-PCR was performed, using RT and PCR primers listed in Table 1. The overall experimental design to detect CoV2-346 by RT-PCR is illustrated in Figure 1B, with detailed seq information of primers and the RNA linker shown in Figure 1C. In brief, RNAs, pretreated with or without T4 PNK, were ligated to a 3′-RNA linker. The RT was done using primers complementary to the RNA linker, followed by PCR using forward primers annealing to 5′-end of CoV2-346 and reverse primers annealing to the last 4 nt of CoV2-346 and the first 15 nt of RNA linker.

Northern Blot

Northern hybridization for tRF5-GluCTC was performed as described (Wang et al., 2013). Briefly, 3 µg RNA was loaded on 15% denaturing polyacrylamide gel with 7 M ureas and then transferred to a positively charged nylon membrane (Amersham Biosciences, NJ, United States). The membrane was hybridized with a ³²P-labeled DNA probe in ULTRAhyb-Oligo solution (Life Technologies, Grand Island, NY, United States), followed by membrane washing and image development. The ³²P-labeled DNA probe for tRF5-GluCTC was 5′-CGC​CGA​ATC​CTA​ACC​ACT​AGA​CCA​CCA-3′.

Statistical Analysis

The experimental results were analyzed using Graphpad Prism 5 software. To compare the sncRNAs expression of NPS between SARS-CoV-2 negative and positive groups, an unpaired two-tailed Mann-Whitney U test was used. To compare the sncRNAs expression in SARS-CoV-2 infected cells and mock-infected cells, an unpaired two-tailed t-test was employed. A p-value < 0.05 was considered to indicate a statistically significant difference. Single and two asterisks represent a p-value of <0.05 and <0.01, respectively. Means ± standard errors (SE) are shown.

Results

T4 PNK-RNA-seq Revealed SARS-CoV-2-Impacted sncRNAs in NPS Samples.

To identify SARS-CoV-2-impacted sncRNAs, we initialized T4 PNK-RNA-seq for the NPS samples from four SARS-CoV-2 positive patients with their ages at 54.3 ± 4.0 years old and four SARS-CoV-2 negative subjects, with matched age at 50.5 ± 10.2 years old. The seq data were analyzed similarly as described in (Wang et al., 2013; Liu et al., 2018). In brief, the sequences with length >15 bp and reads >10 were mapped to the in-house small RNA database containing tRFs, miR/snoRs, and piRs to address redundant tRNA sequences across the genome after removing rRNAs. Unmapped sequences were then mapped to the hg38 human genome to identify other human-derived sequences and their composition. We found that piRNAs and tRFs were the two most abundant sncRNAs in SARS-CoV-2 positive samples. The top-10 ranked piRNAs and tRFs in the SARS-CoV-2 positive group are listed in Supplementary Tables S1, S2, respectively. As shown in Supplementary Table S2, all tRFs were derived from the 5′-ends of tRNAs, therefore tRF5s. Compared with the tRFs and piRNAs, the overall reads of miRNAs were much less (Figure 2A). We also compared the sncRNA profiles between SARS-CoV-2 positive and negative samples. As shown in Figure 2A, while the tRFs consisted of about 14% of all sncRNA counts in the control group, tRFs counts became 42% in the COVID-19 group, demonstrating a significant increase by COVID-19. In contrast, the overall expression of miRNAs and piRNAs was not impacted by COVID-19 (Figure 2A).

FIGURE 2. Impacted sncRNAs by SARS-CoV-2 infection in patients. (A) The relative sequencing frequency of miRNAs, tRFs, piRNAs, and snoRNAs was calculated by normalizing their raw reads with the DEseq2 median ratio method. (B) The volcano plot showed that sncRNAs were differentially expressed between and control group (CN) and SARS-CoV-2 patient group (SARS-COV-2). (C) Heatmap for unsupervised clustering of the differently expressed tRFs with >20 mean of normalized counts in any groups according to Pearson correlation. Data are shown as means ± standard error (SE). The single asterisk represents p values of <0.05.

Differential expression analysis for individual sncRNAs was also performed for SARS-CoV-2 negative and positive groups. As shown in Figure 2B, there were more up-regulated tRFs than down-regulated tRFs, while SARS-CoV-2 down-regulated snoRNAs were more than up-regulated ones. We also listed sncRNAs, which were significantly altered by SARS-CoV-2 in Tables 2–4. The cutoff was set as a fold change >2, with the significance of p < 0.05 in changes by SARS-CoV-2, and the mean of normalized counts >10 in the negative or positive group. The differentially expressed tRFs, miRNAs, and snoRNAs were listed in Tables 2–4, respectively.

TABLE 2. Changes in tRFs by SARS-CoV-2.

TABLE 3. Changes in miRNAs by SARS-CoV-2.

TABLE 4. Changes in snoRNAs by SARS-CoV-2.

As shown in Table 2, tRFs were significantly impacted by SARS-CoV-2 in NPS. Among those, 2 tRFs belong to the tRF1 series, 28 tRFs were tRF5s, and 53 tRFs were tRF3s. However, top-ranked SARS-CoV-2-impacted tRFs all belong to the tRF5. In Figure 2C, the mean of normalized counts >20 in the control (CN) or SARS-CoV-2 positive (SARS-CoV-2) group were selected to plot the heatmap and their sequences were listed in Table 5.

TABLE 5. tRF5s with mean of normalized counts > 20 in Control (CN) or COVID-19 groups.

Experimental Validation of SARS-CoV-2-Impacted tRFs

To validate the seq data, we used modified qRT-PCR to detect the expression of tRF5-GluCTC, tRF5-LysCTT, and tRF5-ValCAC, three top-ranked tRF5s in SARS-CoV-2-positive patients according to the seq data, as we previously described (Choi et al., 2020; Wu et al., 2021). Compared with Northern blot validation, the modified qRT-PCR with T4 PNK pretreatment and 3′-end RNA linker ligation provides the possibility to validate as many tRFs as possible for NPS samples, which usually have a limited yield of RNAs. Our results demonstrated that tRF5-GluCTC, tRF5-LysCTT, and tRF5-ValCA were significantly increased in the SARS-CoV-2 group (Figures 3A–C). Unlike SARS-CoV-2, which could induce tRF5-ValCAC, respiratory syncytial virus (RSV), a negative-sense RNA virus, does not induce tRF5-ValCAC infected cells (Wang et al., 2013), suggesting the change in tRF profile in response to viral infections is virus-specific.

FIGURE 3. Changes in the expression of tRF5s in NPS samples by SRAS-CoV-2. qRT-PCR was performed to detect tRF5-GluCTC (A), tRF5-LysCTT (B), tRF5-ValCAC (C), tRF5-CysGCA (D), tRF5-GlnCTG (E), and tRF5-GlyGCC (F) in the NPS from SARS-COV-2 and control (CN) patients. U6 was used as an internal control. Unpaired two-tailed Mann-Whitney U tests were performed for statistical comparisons. Single and double asterisks represent p values of <0.05 and <0.01, respectively. Data are shown as means ± SE.

Other than the three tRF5s mentioned above, tRF5-CysGCA, tRF5-GlnCTG and tRF5-GlyGCC were also chosen for the validation, as these three tRFs are highly inducible by RSV with the function tRF5-GlyGCC and tRF5-GlnCTG being vital in promoting RSV replication (Choi et al., 2020; Zhou et al., 2017). Although the function of tRF5-CysGCA in RSV is not clear in viral infection, it is important in regulating stress responses and neuroprotection (Ivanov et al., 2014). We validated that tRF5-CysGCA and tRF5-GlnCTG were also significantly enhanced in the SARS-CoV-2 group, compared with the control (CN) group (Figures 3D,E). However, SARS-CoV-2 did not affect tRF5-GlyGCC expression (Figure 3F). Given the fact that RSV induces significantly tRF5-GlyGCC (Wang et al., 2013), the result of Figure 3F supported virus-specific induction of tRFs and tRFs as potential biomarkers of viral infection.

Impacted tRFs in SARS-CoV-2-Infected Cells

A549-ACE2 is a common cell model to study coronaviruses, such as SARS-CoV and SARS-CoV-2 (Mossel et al., 2005; Blanco-Melo et al., 2020; Weston et al., 2020; Buchrieser et al., 2021). Herein, we studied whether the induction of tRFs can be recapitulated in SARS-CoV-2-infected A549-ACE2. As shown in Figures 4A–E, A549-ACE2 cells, at day 4 post-infection (p.i.) of SARS-CoV-2, had dramatically induction of tRF5-GluCTC, tRF5-LysCTT, tRF5-ValCAC, tRF5-CysGCA, and tRF5-GlnCTG. The northern blot also confirmed the induction of tRF5-GluCTC (Figure 4F), which was the most abundant tRF5 among the tested four tRF5s, confirming the liability of our newly developed qRT-PCR.

FIGURE 4. SARS-CoV-2-impacted tRF5s in A549-ACE2 cells. (A–E) A549-ACE2 cells were infected with SARS-CoV-2 at an MOI of 0.1 for 4 days, qRT-PCR was performed to detect tRF5-GluCTC (A), tRF5-LysCTT (B), tRF5-ValCAC (C), tRF5-CysGCA (D), and tRF5-GlnCTG (E). A northern blot was performed to detect tRF5-GluCTC (F). Mock-infected cells (those without SARS-CoV-2 infection) were used as control (CN). Unpaired two-tailed t-tests were performed for statistical comparisons. Single and double asterisks represent p values of <0.05 and <0.01, respectively. Data, shown as means ± SE, are representative of three independent experiments.

We also used primary SAECs in ALI culture, a commonly acknowledged physiology airway infection model (Bhowmick and Gappa-Fahlenkamp, 2016; Chandorkar et al., 2017), to confirm SARS-CoV-2-affected tRFs. SAECs, after a few weeks of ALI culture, have been shown to establish a pseudostratified cell layer that resembles the small airway epithelium as found in vivo (Hiemstra et al., 2018). Moreover, SAECs in ALI cultures have been found to express the receptor for SARS-CoV-2, therefore, a physiologically relevant cell model to study SARS-CoV-2 (Zhang et al., 2020; Zhu et al., 2020; Schweitzer et al., 2021). Therefore, we also studied tRF5s expression in SAECs in ALI culture. As shown in Figures 5A,B, the cilia were oriented towards the upper transwell compartment, after the cells were in ALI culture for 21 days. The differentiated cultures were infected with SARS-CoV-2 at an MOI of 0.1 for 1 or 3 days, followed by viral S gene quantification using qRT-PCR (Figure 5C). Our qRT-PCR confirmed the expression change in tRF5-GlnCTG and tRF5-ValCAC, two tRFs with relatively low abundance in SARS-CoV-2 positive NPS samples and infected A549-ACE2 cells, can also be detected in SAECs in ALI culture (Figures 5D, E). Since the cleaved tRNAs account for a very tiny portion of parent tRNAs, the difference in the induction folds of tRF5-GlnCTG and tRF5-ValCAC should not be determined by the abundance of their parent tRNAs, but possibly resulted from the distinct sensitivities of their parent tRNAs to the cleavage during the infection. Overall, in this study, we established two cell models, A549-AEC2 cells in monolayer culture and SAECs in ALI culture, to characterize SARS-CoV-2-induced tRFs in the future.

FIGURE 5. SARS-CoV-2-affected tRF5s in ALI-cultured SAECs. (A) Histological examination SAECs in ALI culture. After SAECs were in ALI culture for 21 days, the insert was fixed with 4% polyoxymethylene, followed by tissue processing, sectioning, and H&E staining. (B) Transmission electron microscopy (TEM) of SAECs in ALI culture. (C) ALI-cultured SEACs were apically infected with SARS-CoV-2 at an MOI of 0.1 for 1 or 3 days, qRT-PCR was performed to detect SARS-CoV-2 S gene expression. GAPDH was used as an internal control. (D,E) On day 3 postinfection, tRF5-GlnCTG (D) and tRF5-ValCAC (E) expressions were measured using qRT-PCR. Mock-infected cells were used as control (CN). Unpaired two-tailed t-tests were performed for statistical comparisons. Single asterisks represent p values of <0.05. Data are shown as means ± SE and are representative of three independent experiments.

Identification of SARS-CoV-2-Encoded svRNAs

Viral-derived sncRNAs are also an important family of sncRNAs (Tycowski et al., 2015). To investigate whether SARS-CoV-2-encoded svRNAs are produced in the context of SARS-CoV-2 infection, the seq data were also aligned to the complete genome sequence of SARS-CoV-2 isolate Wuhan-Hu-1 (NC_045512.2). Several SARS-CoV-2-derived svRNAs were identified in SARS-CoV-2 positive samples. The eight most abundant SARS-CoV-2-encoded svRNAs are listed in Table 6.

TABLE 6. SARS-CoV-2-encoded svRNAs.

We further analyzed the sequences of svRNAs. Since only RNA <200 bp were selected for the cDNA library, our results should not give svRNAs larger than 200 bp. We found that the length of mapped svRNAs ranged from 17 to 75 nt. Interestingly, svRNAs with the length of 25 nt, 33 nt, and 36 nt were enriched (Figure 6A, two representatives are shown). In Figure 6B, the locations of the top 8 svRNAs along the SARS-CoV-2 genome are shown.

FIGURE 6. Sequence information of virus-derived sncRNAsA. (A) Length distribution of viral small RNAs from two representative patient samples. (B) Two representative visual inspections of the small RNA sequences aligning with the SARS-CoV-2 genome. The names of viral genes and the genome positions (nt) are indicated.

Experimental Validation of SARS-CoV-2-Encoded svRNAs

Among CoV-2-derived svRNAs (sv-CoV2), a 36 nt sv-CoV2, derived from genomic site 346 to site 382 of nsp1 (sv-CoV2-346) had the highest expression. To further validate the presence of sv-CoV2-346, NPS RNAs from two control samples and two COVID-19 samples were treated with T4 PNK and linked to a 3′ RNA linker, and then the RT-PCR was performed. The RT-PCR was also performed without the T4 PNK treatment and RNA linker addition so that the importance of such treatments can be demonstrated. The overall workflow is illustrated in Figures 1B,C. The specific 55 nt RT-PCR products of sv-CoV2-346 were observed in SARS-CoV-2 samples, but not in the control samples, when samples were pretreated with T4 PNK and ligated with an RNA linker (Figure 7A). The length reflected the 36 nt sv-CoV2-346 along with the 3′ RNA linker. In addition, we found that the RT-PCR product of one SARS-CoV-2 sample was more than another one, which was consistent with their read frequency in Seq-data. The presence of sv-CoV2-346 was confirmed in SARS-CoV-2-infected A549-ACE2 cells using RT-PCR (Figure 7B).

FIGURE 7. Experimental validation of sv-CoV2-346. (A) RT-PCR was performed to detect sv-CoV2-346 in NPS samples of two control (CN, SARS-CoV-2 negative) and two SARS-CoV-2 patients. (B) The presence of svCoV2-346 in SARS-CoV-2-infected A549-ACE2 cells. (C) Detection of sv-CoV2-346 needs the pretreatment of T4 PNK. The 3′-end of sv-CoV2-346 does not have –OH, as samples without pretreatment of T4 PNK did not result in sv-CoV2-346 bands. All experiments were independently repeated twice.

Coronavirus-encoded svRNAs have been previously reported to be 18–22 nt long and therefore, share similar lengths with miRNAs (Morales et al., 2017; Cheng et al., 20212021). Coronavirus-encoded svRNAs with lengths longer than 30 nt have not been identified. Herein, we think that the identification of additional SARS-CoV-2-derived svRNAs was benefited from the treatment of T4 PNK and RNA linker ligation at their 3′-end. As shown in Figure 7C, both patient or infected cell samples, without such treatments, did not result in the band presence, supporting the lack of 3′-OH end of sv-CoV2-346 and the necessity of specific T4 PNK treatment for sv-CoV2-346 detection.

Herein, we also initialed to characterize sv-CoV2-346 by predicting the secondary RNA structure of svRNAs. Besides sv-CoV2-346, there were two other svRNAs, sv-CoV2-299 and svCoV2-404, near the region where sv-CoV-346 was derived (Table 7). sv-CoV2-299, sv-CoV2-346, and sv-CoV2-404 were derived from nucleotide 299 to 328, 346 to 382, and 400 to 443, respectively. Therefore, we took the viral genome spanning from 289 to 485, which covers all these three regions with some nt extension on both ends and predicted its RNA secondary structure using RNAfold web server based on minimum free energy to have a clue of biogenesis mechanisms (Hofacker, 2003) (Figure 8A). We found that nucleotides 299, 328, 400, 443, and 382 are all located on loops, implying the cleavage at these five sites along with the single-stranded RNA by ribonuclease (Figure 8A). Only nucleotide 346 was in the middle of the stem (Figure 8A). Interestingly, we found that 68 nt long svRNAs (sv-CoV2-314) overlapped with sv-CoV2-346 (Table 7). We, therefore, took the genome section spanning from 314 to 382 and run the secondary and tertiary structures of sv-CoV2-314 using RNAfold web server and RNAComposer, respectively (Hofacker, 2003; Popenda et al., 2012) (Figures 8B,C). This 68 nt fragment contained three hairpin loops and was folded into an L-shaped-like tertiary structure, and nucleotide 346 was located within the bottom loop (Figures 8B,C). The secondary and tertiary structures of sv-CoV2-314 were similar to tRNA, and nucleotide 346 location was similar to the cleavage site of tRFs (Figures 8D,E). Thus, we speculated that 68 nt sv-CoV2-314 may be the precursor of 36 nt sv-CoV2-346 and the virus may use endonuclease involved in tRF biogenesis to generate viral small RNAs fragments.

TABLE 7. RNA sequence information for those with high raw reads and mapping to viral nucleotides 289–485.

FIGURE 8. Structure prediction of viral small RNA. (A) The secondary structure of viral nucleotides 289–485. The secondary (B) and tertiary (C) structure of sv-CoV2-314. The secondary (D) and tertiary (E) structure of a tRNA example GluCTC. Arrows indicated cleavage sites.

Discussion

In this study, we identified the change in sncRNA expression by SARS-CoV-2. Among sncRNAs, miRNAs have been getting lots of attention in the studies (Vaz et al., 2010; Plieskatt et al., 2014; Max et al., 2018; Hou and Yao, 2021). Standard barcode-ligation-based small RNA-seq are usually designed to capture miRNAs, which usually have 3′-OH ends (Hafner et al., 2008). It is increasingly acknowledged that the 3′-ends of other types of sncRNAs are heterogeneous (Honda et al., 2015), resulting in unsuccessful sequencing barcode ligation in the standard small RNA-seq. In our study, T4 PNK-RNA-seq was employed to profile sncRNAs with heterogeneous ends. Sequencing data revealed that piRNAs and tRFs had higher global expression than miRNAs in NPS (Figure 2A). sncRNAs may carry various unidentified modifications, which are insensitive to T4 PNK treatment. Therefore, T4 PNK-RNA-seq may leave some SARS-CoV-2-impacted sncRNAs unidentified. Herein, the consistency among the seq data, qRT-PCR result, and NB data of tRF5-GluCTC suggested the reliability of T4 PNK-RNA-seq for tRF5 detection.

Notable, among differently expressed tRF5s with mean of normalized counts >20 in control (CN) or SARS-CoV-2 groups (Table 5), we found four tRF5s: tRF5-GluTTC-1-2, tRF5-GlyCCC-6-1, tRF5-LysCTT61 and tRF5-SecTCA-2-1, were not classic tRF5s. While their 3′-ends commonly stop around the anticodon region like classical tRF5s, they lack the first 10-15 nt of the tRNA 5′end. Since they span the complete region of the D loop and the first half of the anticodon loop, we subgrouped and named them as tRF5DC (Supplementary Figure S2). Interestingly, tRF5DC and classic tRF5s were derived from the different tRNA isoacceptors tRNA, suggesting different biogenesis mechanisms of these two type tRFs.

This study further supported that tRFs induction is virus-specific. Previously, we and others have shown that RSV, hepatitis B virus (HBV), and hepatitis C virus (HCV) infections lead to different tRF profile changes (Wang et al., 2013; Selitsky et al., 2015; Zhou et al., 2017; Choi et al., 2020). Compared with tRF induction by RSV, we found that SARS-CoV-2 could induce tRF5-ValCAC, while RSV cannot. On the other hand, we found that tRF5-GlyGCC, which is significantly inducible by RSV, was not induced by SARS-CoV-2. The virus-specific changes in tRFs suggest them to be promising biomarkers for viral infections.

Other than host-derived sncRNAs, sncRNAs can also be derived from viruses. svRNAs have been reported to be involved in the regulation of viral replication, viral persistence, host immune evasion, and cellular transformation (Tycowski et al., 2015). SARS-CoV-2-encoded sncRNAs have been demonstrated by two independent groups (Cheng et al., 20212021; Fu et al., 2021). Using T4-PNK-RNA-seq, several novel svRNAs in SARS-CoV-2 NPS samples were identified. One of them, sv-CoV2-346, was verified to be present in SARS-CoV-2-infected A549-ACE2 cells as well. Due to the limited NPS RNA samples, the leftover RNAs after sequencing were not enough for stem-loop qRT-PCR to validate svRNAs. Thus, we detected sv-CoV2-346 by RT-PCR using the same cDNA generated by the RT step for the qRT-PCR assays for tRF5s. Our RT-PCR revealed a sv-CoV2-346 specific band around 55 nt and a non-specific band. In the future, we will study the relationship between SARS-CoV-2 svRNAs and viral loads/disease severity.

The most widely studied viral sncRNAs are miRNAs-like svRNAs. Both DNA and RNA viruses could encode miRNAs-like svRNAs via Dicer-dependent miRNAs biogenesis pathways (Tycowski et al., 2015; Grundhoff and Sullivan, 2011). Among RNA viruses, cytoplasmic restricted RNA viruses were thought incapable of producing miRNA-like svRNAs. However, accumulating evidence indicates cytoplasmic RNA viruses, such as H5N1 influenza virus, enterovirus 71(EV71), West Nile virus (WNV), SARS-CoV, and SARS-CoV-2, also encode viral miRNAs (Morales et al., 2017; Cheng et al., 20212021; Fu et al., 2021; Perez et al., 2010; Li et al., 2018; Weng et al., 2014; Hussain et al., 2012). These cytoplasmic RNA viruses generate viral miRNAs via multiple non-canonical miRNAs biogenesis mechanisms. Dicer, not Drosha, is involved in WNV and EV71 viral miRNAs generation (Weng et al., 2014; Hussain et al., 2012). H5N1 influenza virus and SARS-CoV encode viral miRNAs in a Dicer-and Drosha-independent way (Morales et al., 2017; Li et al., 2018). Besides viral miRNAs, the induction of functional svRNAs, which do not look like miRNAs, was reported for cytoplasmic RNA viruses (Perez et al., 2010). However, the knowledge on how cytoplasmic RNA viruses produce svRNAs is limited. One of the interesting observations of newly discovered sv-CoV2 is that sv-CoV2-314 may have a similar tertiary structure as tRNAs and function as the potential precursor of 36 nt svCoV2-346 (Figure 8). Whether a tRNA-like shape (three-leafed clover) of svRNAs makes them as prone as tRNAs to the cleavage during SARS-CoV-2 infection will be investigated in the future. Recently, the cleavage of tRNAs has been reported to be regulated by nt modifications and tRNA anticodon loop is enriched with modification (Blanco et al., 2014; Ranjan and Rodnina, 2016). Therefore, it is possible that the anticodon and/or D loop experience nt modification changes in SARS-CoV-2 infection, resulting in the cleavage. It has been also previously reported that the cleavage of tRNAs is mediated by specific ribonuclease(s) (Lee et al., 2009; Wang et al., 2013; Zhou et al., 2017; Choi et al., 2020). Therefore, it is also possible SARS-CoV-2 favors the activation of certain enzymes to enrich the corresponding sncRNA population. How viruses use the host proteins to control the biogenesis of sncRNAs is still unclear and awaits investigation.

In summary, we investigated COVID-19-impacted sncRNAs comprehensively using the NPS samples by T4 PNK-RNA-seq and modified qRT-PCR method. We are aware that our study has some limitations. For example, T4 PNK-RNA-seq may not catch all types of sncRNAs. In addition, our NPS samples were all from outpatient clinics, which set the limitation to study the correlation of tRF changes with the disease severity. We are closely working with our Molecular Diagnosis Laboratory to obtain the samples from outpatient, inpatient, and ICU services so that whether tRFs serve as disease biomarkers can be determined. Our recent publication on the correlation of tRF changes with Alzheimer’s disease severity supports tRFs to be promising disease biomarkers (Wu et al., 2021). Other than biomarkers, the studies of tRFs in viral infections may also provide insight into therapeutic/prophylactic strategy development. In RSV infection, we found some induced tRFs promote viral replication (Wang et al., 2013; Zhou et al., 2017; Choi et al., 2020). Therefore, any mechanisms associated with their biogenesis and function would not only reveal potential targets to control viral replication but also benefit the ncRNA research community.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: Gene Expression Omnibus accession number: GSE193555.

Ethics Statement

The studies involving human participants were reviewed and approved by University of Texas Medical. IRB committee. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

WW and XB wrote the manuscript; WW, E-JC, BW, KZ, contributed to the experiment performance and data analyses. LT, PH, RG, II, LH, and IL contribute to the deep seq data analyses and related methodology establishment. GH and JD developed IBM protocols for sample collection and NSP sample handling. JD, TW, and XB contributed to overall project management and experimental design.

Funding

This work was supported by grants from the NIH R01 AI116812, R21AG069226, and Research Pilot Grant from the Institute for Human Infections and Immunity (IHII) at UTMB to XB and NIH grants R01AI127744 (TW), R21 AI140569 (TW), and R21 EY029112 (TW).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank Darby L. Buck for the manuscript editing.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2022.821137/full#supplementary-material

Supplementary Figure 1 | The pipeline of T4 PNK-RNA-seq and data analyses. (A) The flowchart of T4 PNK-RNA-seq. (B) The flowchart of the sequencing data analyses.

Supplementary Figure 2 | Biogenesis of tRF5s and tRF5DC. tRF5s are generated by the cleavage on the anticodon loop (black arrows); tRF5DCs are derived from the cleavage on the D arm (red or blue arrow) and anticodon loop (black arrows).

References